© ESO, 2017

1. Introduction

The rotation of close-in planets is usually modified by tidal interactions with the central star and reaches a stationary value on timescales typically much shorter than the tidal evolution of orbits (e.g. Hut 1981; Correia 2009). As long as the orbit has some eccentricity, the rotation can stay in non-synchronous configurations, but tidal dissipation also circularizes the orbit, which ultimately results in synchronous motion (the orbital and rotation periods become equal).

Recently, Leconte et al. (2015) used simulations including global climate model (GCM) of the atmosphere of Earth-mass planets in the habitable zone of M-type stars to show that these planets might be in a state of asynchronous rotation (see also Correia et al. 2008). This asynchronous rotation is due to thermal tides in the atmosphere. This same effect was also invoked to explain the retrograde spin of Venus (see Correia & Laskar 2001, 2003). However, for close-in planets, the gravitational tides dominate the thermal tides, so synchronous rotation is believed to be the most likely scenario (Correia et al. 2008; Cunha et al. 2015).

In this paper we investigate another effect that can drive the spin of close-in planets to asynchronous rotation, namely planetary perturbations. Correia & Robutel (2013) showed that in the case of co-orbital planets, planet-planet interactions induce orbital perturbations that can lead to asynchronous spin equilibria, and even chaotic evolution of the spin of the planets. The planets librate around the Lagrangian equilibrium and have oscillations of their mean longitude that prevent the spin synchronization.

We generalize this study to other mean-motion resonances (2:1, 3:2, etc.) for which a similar libration of the mean longitude can be observed. We study both the resonant and the near-resonance cases. While no co-orbital planet has yet been observed, many planets have been found around other mean-motion resonances. Moreover, in some cases, strong planet-planet interactions have been observed, for instance GJ 876 using the radial velocity (RV) technique (Correia et al. 2010), or Kepler-88 (also known as the King of TTVs) using the transit timing variation (TTV) technique (Nesvorný et al. 2013). The TTV technique is particularly promising for such a study since the amplitude and period of the observed TTVs provide a direct measure of the amplitude and period of the orbital libration which affects the spin evolution.

In Sect. 2, we derive the equations of motion for the spin of a planet that is perturbed by another planet. In Sect. 3 we show that both the spin dynamics and the TTVs are dominated by the perturbations of the mean longitude of the planet, and the coefficients of their Fourier series are related to each other. In Sect. 4, we apply our analytical modelling to some observed planets showing significant TTVs and perform numerical simulations. Finally, we discuss our results in Sect. 5.

2. Spin dynamics

We consider a system consisting of a central star with mass m₀, and two companion planets with masses m₁ and m₂, such that m₁,m₂ ≪ m₀. We study the spin evolution of one of the planets, either the inner one or the outer one. The subscript 1 always refers to the inner planet, 2 refers to the outer one, and no subscript refers to the planet whose spin evolution is studied. For simplicity, we assume coplanar orbits and low eccentricities for both planets. We assume that the spin axis of the considered planet is orthogonal to the orbital plane (which corresponds to zero obliquity)¹.

We introduce θ, the rotation angle of the planet with respect to an inertial line, whose evolution is described by (e.g. Murray & Dermott 1999) $\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{6}{\mathit{C}}_{\mathrm{22}}}{\mathit{\zeta }}\frac{\mathit{\mu }}{{\mathit{r}}^{\mathrm{3}}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}\mathit{f}\mathrm{)}\mathit{,}$(1)where ( being the gravitational constant), r is the distance between the planet and the star (heliocentric coordinates), f is the true longitude of the planet, C₂₂ is the Stokes gravity field coefficient that measures the asymmetry in the equatorial axes of the planet, and ζ is the inner structure coefficient that measures the distribution of mass in the planet’s interior (see Appendix C).

Expanding expression (1) in Fourier series of the mean longitude, we obtain $\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{6}{\mathit{C}}_{\mathrm{22}}}{\mathit{\zeta }}\frac{\mathit{\mu }}{{\mathit{a}}^{\mathrm{3}}}\sum _{\mathit{k}\mathrm{\in }Z}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\sin \mathrm{(}\mathrm{2}\mathit{\theta }\mathrm{-}\mathit{k\lambda }\mathrm{+}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{2}\mathrm{\left)}\mathit{\varpi }\mathrm{\right)}\mathit{,}$(2)where a, e, λ, ϖ are the semi-major axis, the eccentricity, the mean longitude, and the longitude of periastron of the planet, respectively. We note that the Hansen coefficient ${\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ is of order | k − 2 | in eccentricity. Thus, the leading term in this expansion corresponds to k = 2, and we have ${\mathit{X}}_{\mathrm{2}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}$. If we only keep this leading term, Eq. (2) simplifies as $\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{\approx }\mathrm{-}\frac{\mathrm{6}{\mathit{C}}_{\mathrm{22}}}{\mathit{\zeta }}\frac{\mathit{\mu }}{{\mathit{a}}^{\mathrm{3}}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}\mathit{\lambda }\mathrm{)}\mathit{.}$(3)If the planet remains unperturbed and follows a Keplerian orbit, the semi-major axis is constant and the mean longitude is given by $\mathit{\lambda }\mathrm{=}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{nt,}$(4)where n is the constant mean motion of the planet. Introducing $\begin{array}{ccc}\mathit{\gamma }& \mathrm{=}& \mathit{\theta }\mathrm{-}\mathrm{(}{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{n}{\mathit{t}}^{\mathrm{)}}\mathit{,}\\ \mathit{\sigma }& \mathrm{=}& \mathit{n}\sqrt{\mathrm{12}\hspace{0.17em}{\mathit{C}}_{\mathrm{22}}\mathit{/}\mathit{\zeta }}\mathit{,}\end{array}$we can rewrite Eq. (3) as $\begin{array}{c}\mathrm{¨}\\ \mathit{\gamma }\end{array}\mathrm{\approx }\mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}\sin \mathrm{2}\mathit{\gamma ,}$(7)which is the equation of a simple pendulum, with a stable equilibrium point at γ = 0. This equilibrium corresponds to the exact synchronization (). The parameter σ (Eq. (6)) measures the width of the synchronous resonance, and also corresponds to the libration frequency at exact resonance.

We now consider planet-planet interactions that disturb the orbits, in particular a_i and λ_i. The action canonically conjugated with λ_i is the circular angular momentum of planet i, ${\mathrm{\Lambda }}_{\mathit{i}}\mathrm{=}{\mathit{\beta }}_{\mathit{i}}\sqrt{{\mathit{\mu }}_{\mathit{i}}{\mathit{a}}_{\mathit{i}}}$, where β_i = m_im₀/ (m₀ + m_i). We also introduce the angular momentum deficit (AMD) of the planets, ${\mathit{D}}_{\mathit{i}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}\right)$ which are the actions canonically conjugated with ϖ_i. The Hamiltonian of the three-body problem reads $\mathrm{\mathscr{H}}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{\right)}\mathrm{+}{\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathit{i}}\mathit{,}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}{\mathit{D}}_{\mathit{i}}\mathit{,}{\mathit{\varpi }}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(8)where ℋ₀ and ℋ₁ are the Keplerian and perturbative part of the Hamiltonian. We have ${\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{1}}}{\mathrm{2}{\mathit{a}}_{\mathrm{1}}}\mathrm{-}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{2}}}{\mathrm{2}{\mathit{a}}_{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\beta }}_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{2}{\mathrm{\Lambda }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\beta }}_{\mathrm{2}}^{\mathrm{3}}}{\mathrm{2}{\mathrm{\Lambda }}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}$(9)and ℋ₁ is of order one in planet masses (over star mass). The equations of motion read $\begin{array}{ccc}\mathit{\lambda ̇}\mathit{i}& \mathrm{=}& \frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathit{i}}}\mathrm{=}\frac{{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}{\mathit{\beta }}_{\mathit{i}}^{\mathrm{3}}}{{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{3}}}\mathrm{+}\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{1}}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathit{i}}}\mathit{,}\\ \mathrm{\Lambda ̇}\mathit{i}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }\mathrm{\mathscr{H}}}{\mathit{\partial }{\mathit{\lambda }}_{\mathit{i}}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }{\mathrm{\mathscr{H}}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\lambda }}_{\mathit{i}}}\mathrm{·}\end{array}$However, it can be shown that the perturbative part has a weak dependency on Λ_i (e.g. Delisle et al. 2014). At leading order in eccentricity, we have $\mathit{\lambda ̇}\mathit{i}\mathrm{\approx }\frac{{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}{\mathit{\beta }}_{\mathit{i}}^{\mathrm{3}}}{{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{3}}}\mathrm{\approx }{\mathit{n}}_{\mathit{i}}{\left(\frac{{\mathrm{\Lambda }}_{\mathit{i}}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}\right)}^{-3}\mathrm{\approx }{\mathit{n}}_{\mathit{i}}\mathrm{-}\mathrm{3}{\mathit{n}}_{\mathit{i}}\frac{{\mathrm{\Lambda }}_{\mathit{i}}\mathrm{-}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}\mathrm{·}$(12)We assume that the orbital motion is not chaotic and thus quasi-periodic. We introduce the quasi-periodic decomposition of Λ$\mathrm{\Lambda }\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{0}}\left(\mathrm{1}\mathrm{+}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{j}}\cos \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\right)\mathit{,}$(13)where ν_j can be any combination of the frequencies of the system. For a resonant system (see Appendix B), the frequency of libration in the resonance is the dominant term of the decomposition. For a system that is close to resonance (k₂:k₁) but outside of it (see Appendix A), the dominant frequency is the frequency of circulation k₂n₂ − k₁n₁. However, for this computation we can keep things general and use the generic decomposition of Eq. (13). We replace it in Eq. (12) to obtain the evolution of the mean longitude $\mathit{\lambda ̇}\mathrm{\approx }\mathit{n}\mathrm{-}\mathrm{3}\mathit{n}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{j}}\cos \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathit{.}$(14)We thus have $\begin{array}{ccc}\mathit{\lambda }& \mathrm{\approx }& {\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{nt}\mathrm{-}\sum _{\mathit{j}}\mathrm{3}{\mathit{A}}_{\mathit{j}}\frac{\mathit{n}}{{\mathit{\nu }}_{\mathit{j}}}\mathrm{(}\sin \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathrm{-}\sin {\mathit{\phi }}_{\mathit{j}}\mathrm{)}\\ & \mathrm{=}& \end{array}$(15)where the coefficients ${\mathit{\alpha }}_{\mathit{j}}\mathrm{=}\mathrm{3}{\mathit{A}}_{\mathit{j}}\mathit{n}\mathit{/}\mathit{\nu }j$(16)have the dimension of angles and correspond to the amplitudes of each term in the decomposition. For the sake of simplicity of the computations, we assume that these amplitudes remain small, but in principle they could reach 180°. With this approximation, we have at first order in α_j$\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}\mathit{\lambda }\mathrm{)}\mathrm{\approx }\sin \mathrm{2}\mathit{\gamma }\mathrm{+}\mathrm{2}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\sin \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\cos \mathrm{2}\mathit{\gamma ,}$(17)with $\mathit{\gamma }\mathrm{=}\mathit{\theta }\mathrm{-}\left({\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{nt}\mathrm{+}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\sin {\mathit{\phi }}_{\mathit{j}}\right)\mathit{.}$(18)From Eq. (13) we deduce $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{a}}^{\mathrm{3}}}& \mathrm{\approx }& \frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{0}}^{\mathrm{3}}}\left(\mathrm{1}\mathrm{-}\mathrm{6}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{j}}\cos \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\right)\\ & \mathrm{\approx }& \frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{0}}^{\mathrm{3}}}\left(\mathrm{1}\mathrm{-}\mathrm{2}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\cos \mathrm{(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\right)\end{array}$(19)and finally, using the expressions of Eqs. (17), (19) in Eq. (3), we obtain $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\gamma }\end{array}& \mathrm{\approx }& \mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}[\hspace{0.17em}\hspace{0.17em}\sin \mathrm{2}\mathit{\gamma }\\ & & \mathrm{+}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\right)\sin \mathrm{2}\left(\mathit{\gamma }\mathrm{+}\frac{{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}}{\mathrm{2}}\right)\\ & & \mathrm{+}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\right)\sin \mathrm{2}\left(\mathit{\gamma }\mathrm{-}\frac{{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{+}\mathit{\pi }}{\mathrm{2}}\right)]\mathrm{·}\end{array}$(20)The term (sin2γ) corresponding to the synchronous resonance still exists. However, for each frequency appearing in the quasi-periodic decomposition of the perturbed orbital elements (see Eqs. (13), (15)), two terms appear (at first order in α_j) in Eq. (20), corresponding to a sub-synchronous resonance (), and a super-synchronous resonance (). This splitting of the synchronous resonance was found in the case of co-orbital planets in Correia & Robutel (2013). We note that if we do not neglect the non-leading terms in the Fourier expansion (Eq. (2), the classical spin-orbit resonances () appear, as do new resonances of the form . Moreover, if the amplitudes α_j are not small, the series should be developed at a higher degree in α_j, and resonances of the type would appear (see Leleu et al. 2015, for the co-orbital case).

Different dynamical regimes can be observed depending on the values of α, ν, and σ. In particular, a chaotic evolution of the spin is expected when the separation between two resonances is of the order of the width of these resonances (Chirikov 1979). Table 1 describes qualitatively the dynamics of the spin of the planet, as a function of the amplitude (α) and frequency (ν) of the perturbing term.

3. TTV as a probe for the spin dynamics

The TTV of a planet is a very good probe that can be used to estimate the main frequencies appearing in the quasi-periodic decomposition of Eqs. (13), (15), and the associated amplitudes. As we did for the computation of the spin evolution, we consider coplanar planets with low eccentricities. We take the origin of the longitudes as the observer direction such that the transit occurs when f = 0 (f being the true longitude of the planet). The true longitude can be expressed as a Fourier series of the mean longitude ${\mathrm{e}}^{\mathrm{i}\mathit{f}}\mathrm{=}\sum _{\mathit{k}\mathrm{\in }Z}{\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{k\lambda }\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{k}\mathrm{\left)}\mathit{\varpi }\mathrm{\right)}}$(21)where the Hansen coefficient ${\mathit{X}}_{\mathit{k}}^{\mathrm{0}\mathit{,}\mathrm{1}}$ is of order | k − 1 | in eccentricity. For the sake of simplicity, we only keep the leading order term f ≈ λ. We thus have (see Eq. (15)) $\mathit{f}\mathrm{\approx }{\mathit{\lambda }}_{\mathrm{0}}\mathrm{+}\mathit{nt}\mathrm{-}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\mathrm{\left(}\sin \mathrm{\right(}{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathrm{-}\sin {\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathit{.}$(22)We introduce t_k, the time of the kth transit. We have f(t_k) = k2π, thus the TTVs are given by ${\mathrm{TTV}}_{\mathit{k}}\mathrm{=}{\mathit{t}}_{\mathit{k}}\mathrm{-}\mathit{kP}\mathrm{\approx }\sum _{\mathit{j}}\frac{{\mathit{\alpha }}_{\mathit{j}}}{\mathit{n}}\mathrm{\left(}\sin \mathrm{\right(}{\mathit{\nu }}_{\mathit{j}}{\mathit{t}}_{\mathit{k}}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathrm{-}\sin {\mathit{\phi }}_{\mathit{j}}\mathrm{)}\mathrm{-}\frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathit{n}}\mathit{,}$(23)where P = 2π/n is the orbital period. Therefore, the quasi-periodic decomposition of the TTV signal directly provides the amplitudes α_j and frequencies ν_j that we need in order to analyze the spin evolution. We are mainly interested in resonant and near-resonant systems. In these cases, one term is leading the expansion with a period much longer than the orbital period (libration period for the resonant case, and circulation period for the near-resonant case, see Appendices A, B). For these long periods, the amplitudes and frequency are well determined as long as the number of transits is sufficient to cover the whole period. In particular, there are no sampling/aliasing issues that could arise for periods that are close to the orbital period.

For the sake of simplicity we assume that one term is leading the TTVs, such that (see Eq. (23)) ${\mathrm{TTV}}_{\mathit{k}}\mathrm{\approx }\frac{\mathit{\alpha }}{\mathit{n}}\left(\sin \mathrm{(}\mathit{\nu }{\mathit{t}}_{\mathit{k}}\mathrm{+}\mathit{\phi }\mathrm{)}\mathrm{-}\sin \mathrm{\left(}\mathit{\phi }\mathrm{\right)}\right)\mathrm{-}\frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathit{n}}\mathit{,}$(24)and (see Eq. (20)) $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\gamma }\end{array}\mathrm{\approx }\mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}[\sin \mathrm{2}\mathit{\gamma }& \mathrm{+}& \mathit{\alpha }\sin \mathrm{2}\left(\mathit{\gamma }\mathrm{+}\frac{\mathit{\nu t}\mathrm{+}\mathit{\phi }}{\mathrm{2}}\right)\\ & \mathrm{+}& \mathit{\alpha }\sin \mathrm{2}\left(\mathit{\gamma }\mathrm{-}\frac{\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{+}\mathit{\pi }}{\mathrm{2}}\right)]\mathit{,}\end{array}$(25)where we assume (ν ≪ n). The most interesting systems for our study (see Table 1) are those for which α is not negligible, such that the width of the resonances at is not negligible compared to the synchronous resonance. Such systems can be locked in sub/super-synchronization or even show chaotic evolution of the spin (see Correia & Robutel 2013, for the co-orbital case). The TTVs provide a determination of the parameters α and ν (and of the phase φ).

4. Application to planets with TTVs

In this section we apply the results obtained in Sects. 2 and 3 to real planetary systems that show large TTVs, and we perform numerical simulations in the conservative case (Sect. 4.1) and dissipative case (Sect. 4.2). We have chosen two examples: KOI-227 b, which is trapped in a mean-motion resonance, and Kepler-88 b, which is near (but not trapped in) a mean-motion resonance. In both cases the perturber is not observed to transit but is inferred from the TTV signal. In addition, KOI-227 b is considered a rocky planet with a permanent equatorial asymmetry C₂₂ ≠ 0 (Eq. (6)), while Kepler-88 b is a gaseous planet for which the C₂₂ value is likely very close to zero (Campbell & Synnott 1985).

4.1. Conservative evolution

4.1.1. KOI-227  b (rocky planet, in resonance)

Fig. 1 Spin dynamics of KOI-227 b in the conservative case, with C22 = 1.4 × 10-7 (top, permanent deformation), 6.6 × 10-7 (middle, maximum deformation in the sub/super-synchronous resonances), and 1.9 × 10-6 (bottom, maximum deformation in the synchronous resonance). The left column shows the main frequency η of θ(t) for different initial values of (and with γ(0) = 0). The colour gives the derivative . Blue dots correspond to resonant motion, green dots to non-resonant regular motion, and red dots to chaotic motion. The vertical black line corresponds to the synchronization (η = n). The two grey lines correspond to the main sub/super-synchronous resonances (η = n ± ν/ 2). The dashed and dotted grey lines correspond to higher order sub/super-synchronous resonances (η = n + kν/ 2, k = ± 2, ± 3). The right column shows the same colour index (), but both γ(0) and are varied (2d maps). The vertical white line highlights the initial conditions taken in the left column (γ(0) = 0).

KOI-227 hosts at least two planets (see Nesvorný et al. 2014), but only one (KOI-227 b) is known to transit. This planet has a radius of 2.23 R_⊕, a period of about 18 d, and TTVs with an amplitude of at least 10 h have been observed (see Nesvorný et al. 2014). In terms of angular amplitude (see Eq. (15)), this corresponds to α ≳ 8°. The main TTV period is about 4.5 yr, thus ν/n ≈ 0.011.

Since the pertubing planet is not detected directly, the orbital parameters of the system cannot be completely solved for, due to degeneracies (see Nesvorný et al. 2014). Three possible families of solutions have been proposed by Nesvorný et al. (2014), corresponding to an outer 2:1 or 3:2 resonance or an inner 3:2 resonance between the observed planet and the perturber. For these three solutions, the planets must stay inside the resonance. The outer 2:1 configuration is favoured by the data but the two other configurations cannot be ruled out (see Nesvorný et al. 2014). The best-fitting solution is the 2:1 configuration (χ² = 37.6). The mass of KOI-227 b is 37.5 M_⊕ for this solution, and its density would thus be 18.6 g cm^-3, which seems very high. However, the orbital parameters, and masses are not very well constrained. As an example, we refitted the orbital parameters of the system, using the same TTV data as Nesvorný et al. (2014), but imposing the density of KOI-227 b to be the same as that of the Earth. The mass of KOI-227 b is thus set to 11.09 M_⊕. The obtained solution has a χ² of 39.1, which is still better than the best-fitting solution in other resonances (χ² = 51.5 for the outer 3:2 resonance, and χ² = 82.2 for the inner 3:2 resonance, see Nesvorný et al. 2014).

For this illustration, we adopted the mass of 11.09 M_⊕ for KOI-227 b. We also imposed the system to be coplanar in order to simplify the problem. We thus refitted the model imposing zero inclination between the planets, which also provides a good fit to the data (χ² = 53.3). The obtained solution is given in Table 2.

Using the mass (11.09 M_⊕) and radius (2.23 R_⊕) of the planet, we estimate its permanent deformation (see Appendix C) ${\mathit{C}}_{\mathrm{22}\mathit{,r}}\mathrm{\approx }\mathrm{1.4}\mathrm{\times }{\mathrm{10}}^{-7}\mathit{.}$(26)This corresponds to σ/n ≈ 2.2 × 10^-3 and σ/ν ≈ 0.20, which means that the sub/super-synchronous resonances () are well separated from the synchronous resonance, and that the planet could be locked in any of these resonances.

In addition to the permanent deformation, the tidal deformation could also play an important role in the spin dynamics of the planet. In particular, if the planet is captured in one of the spin-orbit resonances, the C₂₂ increases due to the tidal deformation. We estimate the maximum deformation of the planet (see Appendices C and D) in the synchronous resonance ${\mathit{C}}_{\mathrm{22}\mathit{,}\mathrm{sync}\mathrm{.}}\mathrm{\approx }\mathrm{1.9}\mathrm{\times }{\mathrm{10}}^{-6}\mathit{,}$(27)and in the sub/super-synchronous resonances ${\mathit{C}}_{\mathrm{22}\mathit{,}\mathrm{sub}\mathrm{/}\mathrm{super}}\mathrm{\approx }\mathrm{6.6}\mathrm{\times }{\mathrm{10}}^{-7}\mathit{.}$(28)These values correspond to σ/n ≈ 8.3 × 10^-3 and σ/ν ≈ 0.77 (synchronous resonance), and σ/n ≈ 4.9 × 10^-3 and σ/ν ≈ 0.45 (sub/super-synchronous resonances). As σ/ν approaches unity, the resonances get closer and closer, which may induce a chaotic evolution of the spin.

To study in more details the spin dynamics in these different resonances, we perform numerical simulations of the spin in the conservative case. We substitute in Eq. (1) the orbital solution given by a classical N-body integrator, and integrate it to obtain the evolution of the rotation angle (θ). Figure 1 shows a frequency analysis (using the NAFF algorithm, see Laskar 1988, 1990, 1993) of the spin of KOI-227 b in the conservative case, and assuming C₂₂ = 1.4 × 10^-7 (permanent deformation), 6.6 × 10^-7 (maximum deformation in the sub/super-synchronous resonances), and 1.9 × 10^-6 (maximum deformation in the synchronous resonance).

For C₂₂ = 1.4 × 10^-7 (see Fig. 1 top), we observe that the synchronous resonance is stable, as are the main super/sub-synchronous resonances (η = n ± ν/ 2). The widths of these three resonances are comparable, which indicates that the capture probability in any of these resonances should be similar. For C₂₂ = 6.6 × 10^-7 (see Fig. 1 middle), we observe that the three resonances are surrounded by a large chaotic area. However, a stable region is still visible in each of the three resonances. This means that the sub/super-synchronous resonances remain stable even if the tidal deformation increases to its maximum value after the resonant capture. Finally, for C₂₂ = 1.9 × 10^-6 (see Fig. 1 bottom), we observe a very large chaotic area that encompasses all the resonances. We still observe three areas of stability, which correspond to the synchronous resonance and to the second-order sub/super-synchronous resonances η = n ± ν). Stable capture in the synchronous resonance is thus still possible.

We conclude that in the case of KOI-227 b, stable captures in the synchronous, and sub/super-synchronous resonances (η = n ± ν/ 2) are all possible, and should have comparable probabilities. However, for the highest C₂₂ value some chaotic evolution can be expected before the rotation enters a stable island.

4.1.2. Kepler-88 b (gaseous planet, close to resonance)

Fig. 2 Same as Fig. 1 but for Kepler-88 b. We neglect the permanent deformation of the planet (C22,r = 0), and show the spin dynamics for C22 = 2.2 × 10-6 (top, maximum deformation in the sub/super-synchronous resonances), and 6.7 × 10-6 (bottom, maximum deformation in the synchronous resonance).

Kepler-88 b, also referred to as KOI-142 b, or as the King of TTVs, is the planet exhibiting the largest TTV known today. The TTV amplitude is α ≈ 16° (amplitude of 12 h, compared with an orbital period of 10.95 d, see Nesvorný et al. 2013), and the TTV period is about 630 d, thus we have ν/n ≈ 0.017. As is true for KOI-227 b, the perturber is not directly observed. However, the TTV signal and the transit duration variation (TDV) are sufficient in this case to obtain a unique orbital solution (see Nesvorný et al. 2013). We reproduce in Table 3 the orbital parameters obtained by Nesvorný et al. (2013). We note that the mutual inclination is very small (about 3°), and we neglect it in the following.

The bulk density of Kepler-88 b is about 0.87 g cm^-3 (Nesvorný et al. 2013), which means that this planet is mainly gaseous. Its permanent deformation is thus probably very small and we neglect it (C_22,r ≈ 0). However, the tidally induced deformation of the planet, if it is captured in a resonance, is not negligible (see Appendices C and D) $\begin{array}{ccc}& & \\ & & {\mathit{C}}_{\mathrm{22}\mathit{,}\mathrm{sub}\mathrm{/}\mathrm{super}}\mathrm{\approx }\mathrm{2.2}\mathrm{\times }{\mathrm{10}}^{-6}\mathit{.}\end{array}$This corresponds to σ/n ≈ 0.018 and σ/ν ≈ 1.0 (maximum deformation in the synchronous resonance), and σ/n ≈ 0.010 and σ/ν ≈ 0.59 (maximum deformation in the sub/super-synchronous resonances).

Figure 2 shows the spin dynamics of Kepler-88 b in the conservative case and with both estimates of the deformation. For C₂₂ = 2.2 × 10^-6 (see Fig. 2 top), we observe that a large chaotic area is surrounding the synchronous and super/sub-synchronous resonances (η = n ± ν/ 2). A stable area is visible at the centre of the synchronous resonance, but not in the super/sub-synchronous resonances. Therefore, permanent capture in non-synchronous resonances is not possible. For C₂₂ = 6.7 × 10^-6 (see Fig. 2, bottom), a very small island of stability is still present in the synchronous resonance. Therefore, stable capture in the synchronous resonance should be possible but might be difficult to achieve.

4.2. Numerical simulations with tidal dissipation

Tidal interactions with the star have a double effect on the planet: deformation and dissipation. The deformation occurs because the mass distribution inside the planet adjusts to the tidal potential. The dissipation occurs because this adjustment is not instantaneous, so there is a lag between the perturbation and the maximum deformation. As seen in Sect. 4.1, the deformation is very important, since different values of the C₂₂ can have very different implications for the spin dynamics.

In this section we also take into account the dissipative part of the tidal effect, which slowly modifies the spin rotation rate of the planet, and might drive it into the different configurations described in Sect. 4.1. In order to get a comprehensive picture of the spin dynamics of the considered planets, and especially to estimate capture probabilities in the different spin-orbit resonances, we run numerical simulations that take into account both the tidal deformation and the tidal dissipation.

Viscoelastic rheologies have been shown to reproduce the main features of tidal effects (for a review of the main models, see Henning et al. 2009). One of the simplest models of this kind 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). In this case, the planet can respond as an elastic solid or as a viscous fluid, depending on the frequency of the perturbation. The response of the planet to the tidal excitation is modelled by the parameter τ, which corresponds to the relaxation time of the planet².

We adopt here a Maxwell viscoelastic rheology using a differential equation for the gravity field coefficients (Correia et al. 2014). This method tracks the instantaneous deformation of the planet, and therefore allows us to correctly take into account the gravitational perturbations from the companion body. The complete equations of motion governing the orbital evolution of the system in an astrocentric frame are (Rodríguez et al. 2016) $\begin{array}{ccc}{\begin{array}{c}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{1}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\mu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{3}}}{{r}}_{\mathrm{1}}\mathrm{+}\mathrm{𝒢}{\mathit{m}}_{\mathrm{2}}\left(\frac{{{r}}_{\mathrm{2}}\mathrm{-}{{r}}_{\mathrm{1}}}{\mathrm{|}{{r}}_{\mathrm{2}}\mathrm{-}{{r}}_{\mathrm{1}}{\mathrm{|}}^{\mathrm{3}}}\mathrm{-}\frac{{{r}}_{\mathrm{2}}}{{\mathit{r}}_{\mathrm{2}}^{\mathrm{3}}}\right)\mathrm{+}{f}\mathit{,}\\ {\begin{array}{c}\mathrm{¨}\\ {r}\end{array}}_{\mathrm{2}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\mu }}_{\mathrm{2}}}{{\mathit{r}}_{\mathrm{2}}^{\mathrm{3}}}{{r}}_{\mathrm{2}}\mathrm{+}\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}\left(\frac{{{r}}_{\mathrm{1}}\mathrm{-}{{r}}_{\mathrm{2}}}{\mathrm{|}{{r}}_{\mathrm{1}}\mathrm{-}{{r}}_{\mathrm{2}}{\mathrm{|}}^{\mathrm{3}}}\mathrm{-}\frac{{{r}}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{3}}}\right)\mathrm{+}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}}{{\mathit{\mu }}_{\mathrm{1}}}{f}\mathit{,}\end{array}$where r_i is the position vector of the planet i (astrocentric coordinates), f is the acceleration arising from the potential created by the deformation of the inner planet (Correia et al. 2014) $\begin{array}{ccc}{f}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}{\mathit{\mu }}_{\mathrm{1}}{\mathit{R}}^{\mathrm{2}}}{\mathrm{2}{\mathit{r}}_{\mathrm{1}}^{\mathrm{5}}}{\mathit{J}}_{\mathrm{2}}{{r}}_{\mathrm{1}}\\ & & \mathrm{-}\frac{\mathrm{9}{\mathit{\mu }}_{\mathrm{1}}{\mathit{R}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{5}}}\left[{\mathit{C}}_{\mathrm{22}}\cos \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathrm{-}{\mathit{S}}_{\mathrm{22}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\right]{{r}}_{\mathrm{1}}\\ & & \mathrm{+}\frac{\mathrm{6}{\mathit{\mu }}_{\mathrm{1}}{\mathit{R}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{5}}}\left[{\mathit{C}}_{\mathrm{22}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathrm{+}{\mathit{S}}_{\mathrm{22}}\cos \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\right]{k}\mathrm{\times }{{r}}_{\mathrm{1}}\mathit{,}\end{array}$(33)and k is the unit vector normal to the orbital plane of the inner planet. The torque acting to modify the inner planet rotation is $\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{6}\mathrm{𝒢}{\mathit{m}}_{\mathrm{0}}}{\mathit{\zeta }{\mathit{r}}_{\mathrm{1}}^{\mathrm{3}}}\left[{\mathit{C}}_{\mathrm{22}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathrm{+}{\mathit{S}}_{\mathrm{22}}\cos \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\right]\mathit{.}$(34)The inner planet is deformed under the action of self rotation and tides. Therefore, the gravity field coefficients can change with time as the shape of the planet is continuously adapting to the equilibrium figure. According to the Maxwell viscoelastic rheology, the deformation law for these coefficients is given by (Correia et al. 2014) $\begin{array}{ccc}& & {\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathit{\tau }\mathit{J̇}\mathrm{2}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\frac{\mathit{\theta ̇}\mathrm{2}{\mathit{R}}^{\mathrm{3}}}{\mathrm{3}\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}}\mathrm{+}{\mathit{k}}_{\mathrm{f}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathrm{2}{\mathit{m}}_{\mathrm{1}}}{\left(\frac{\mathit{R}}{{\mathit{r}}_{\mathrm{1}}}\right)}^{\mathrm{3}}\mathit{,}\\ & & {\mathit{C}}_{\mathrm{22}}\mathrm{+}\mathit{\tau }\mathit{Ċ}\mathrm{22}\mathrm{=}{\mathit{C}}_{\mathrm{22}\mathit{,r}}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{{\mathit{m}}_{\mathrm{1}}}{\left(\frac{\mathit{R}}{{\mathit{r}}_{\mathrm{1}}}\right)}^{\mathrm{3}}\cos \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathit{,}\\ & & {\mathit{S}}_{\mathrm{22}}\mathrm{+}\mathit{\tau }\mathit{Ṡ}\mathrm{22}\mathrm{=}\mathrm{-}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{{\mathit{m}}_{\mathrm{1}}}{\left(\frac{\mathit{R}}{{\mathit{r}}_{\mathrm{1}}}\right)}^{\mathrm{3}}\sin \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathit{,}\end{array}$(35)where k_f is the fluid second Love number for potential. The relaxation times τ are totally unknown for exoplanets, but if the tidal quality dissipation Q-factor can be estimated, then an equivalent τ can be obtained (see Correia et al. 2014). To cover all possible scenarios, in our numerical simulations we adopt a wide spectrum of τ values: log ₁₀τ_yr = − 5 to 4, with step 1.

4.2.1. KOI-227 b (rocky planet, in resonance)

For rocky planets, such as the Earth and Mars, we have Q = 10 (Dickey et al. 1994) and Q = 80 (Lainey et al. 2007), respectively. We then compute for the Earth τ = 1.6 d, and for Mars τ = 14.7 d. 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 d). Although these values provide a good estimation 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 due to 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). Therefore, we conclude that the deformation timescale of rocky planets can range from a few days up to thousands of years.

In our numerical simulations we use the initial conditions from Table 2. The initial rotation period is set at 15.5 d, which corresponds to . Since the libration frequency of KOI-227 b is ν/n ≈ 0.011, the initial rotation rate is completely outside the resonant area, which is bounded by (Fig. 1). For any τ value, the rotation rate of the planet decreases due to tidal effects until it approaches this area, where multiple spin-orbit resonances are present.

Capture in resonance is a stochastic process. Therefore, for each τ value we ran 1000 simulations with slightly different initial rotation rates. The step between each initial condition was chosen such that the moment at which each simulation crosses the resonance spreads equally over one eccentricity cycle. In Fig. 3 we show some examples of evolution for different τ values. The evolution timescale changes with τ because the Q-factor was also modified.

Fig. 3 Examples of the spin evolution of KOI-227 b and corresponding global C22 value for different τ values. We use the initial conditions from Table 2, and the initial rotation period is set at 15.5 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.

For τ ≤ 10^-2 yr, the spin is in the low frequency regime (τn ≪ 1), which is usually known as the viscous or linear (e.g. Singer 1968; Mignard 1979). As a consequence, the rotation rate evolves into pseudo-synchronous equilibrium (e.g. Correia et al. 2014) $\mathit{\theta ̇}\mathit{/}\mathit{n}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{6}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathrm{𝒪}\mathrm{\left(}{\mathit{e}}^{\mathrm{4}}\mathrm{\right)}\mathit{.}$(36)The average eccentricity of KOI-227 b’s orbit is 0.057, which gives a value for the pseudo equilibrium of  ³. This value is already very close to the synchronous resonance, but since the libration width is σ/n ≈ 0.002 (Eq. (26)), we are still outside the resonant area. Therefore, in this regime initial prograde rotations never cross any spin-orbit resonances.

Fig. 4 Different final spin evolution of KOI-227 b and corresponding global C22 value for τ = 104 yr. We use the initial conditions from Table 2, and the initial rotation period is set at 15.5 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.

For 10^-1 ≤ τ ≤ 10¹ yr, the spin is in a transition of frequency regime (τn ~ 1). The rotation rate still evolves into a pseudo-synchronous equilibrium, but its value is below that provided by expression (36), and lies inside the libration width of the synchronous resonance (see Fig. 4 in Correia et al. 2014). As a result, for these τ values the spin is always captured in the synchronous resonance. Capture in higher order spin-orbit resonances is possible but with a very low probability (<1%); we did not obtain any examples in our simulations.

For τ ≥ 10² yr, the spin is in the high frequency regime (τn ≫ 1). In this regime the tidal torque has multiple equilibria that coincide with the spin-orbit resonances (see Fig. 4 in Correia et al. 2014). Therefore, capture in the asynchronous higher order resonances becomes a real possibility. In Table 4 we list the final distribution in the different resonances for each τ value. We observe that as τ increases, the number of captures in higher order resonances also increases. Indeed, for high τ values, the C₂₂ is able to retain its tidal deformation (Eqs. (27) and (28)) for longer periods of time, increasing the libration width of the individual resonances. When the C₂₂ reaches its maximum tidal deformation, some individual resonances overlap, which results in chaotic motion around these resonances, including synchronous resonance (Fig. 1). An interesting consequence is that for τ = 10⁴ yr, the sub-synchronous resonance can be reached after some wandering in this chaotic zone. In Fig. 4 we show four examples of capture in each resonance for this τ value.

4.2.2. Kepler-88 b (gaseous planet, close to resonance)

For gaseous planets we have Q ~ 10⁴ (Lainey et al. 2009, 2012), which gives τ values of a few minutes assuming that most of the dissipation arises in the convective envelope. However, the cores of these planets also experience tidal effects, which in some cases can be equally strong (Remus et al. 2012; Guenel et al. 2014). In addition, other tidal mechanisms, such as the excitation of inertial waves, are expected to take place, which also enhance the tidal dissipation (e.g. Ogilvie & Lin 2004; Favier et al. 2014). Therefore, the full deformation of gaseous planets may be of the order of a few years or even decades (for a review see Socrates et al. 2012).

In our numerical simulations we use the initial conditions from Table 3. The initial rotation period is set at 10 d, which corresponds to . Since the libration frequency of Kepler-88 b is ν/n ≈ 0.017, the initial rotation rate is completely outside the resonant area, which is bounded by (Fig. 2). For any τ value, the rotation rate of the planet decreases due to tidal effects until it approaches this area, where multiple spin-orbit resonances are present. In Fig. 5 we show some examples of evolution for different τ values.

Fig. 5 Examples of the spin evolution of Kepler-88 b and corresponding global C22 value for different τ values. We use the initial conditions from Table 3, and the initial rotation period is set at 10 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.

As in the case of KOI-227 b, for τ ≤ 10^-2 yr, the spin is in the low frequency regime (τn ≪ 1), and the rotation rate evolves into the pseudo-synchronous equilibrium (Eq. (36)). The average eccentricity of Kepler-88 b orbit is 0.065, which gives for the pseudo equilibrium . The libration width for a maximum value of C₂₂ is σ/n ≈ 0.008 (Eq. (27)), so in this regime initial prograde rotations never cross any spin-orbit resonances (Fig. 5, top). However, for 10^-1 ≤ τ ≤ 10⁰ yr, the equilibrium value is already inside the libration width of the synchronous resonance. Thus, for these τ values the spin can be captured in the synchronous resonance (Fig. 5, middle).

For τ ≥ 10¹ yr the spin is already in the high frequency regime, which means that capture in asynchronous resonances could be possible. However, the equilibrium C₂₂ ≈ 6.7 × 10^-6 of Kepler-88 b is large enough so that the libration zones of individual resonances merge (Eq. (29)). As a consequence, as explained in Sect. 4.1.2, a large chaotic zone around the synchronous resonance is expected. Indeed, in all simulations we observe a chaotic behaviour for the spin (Fig. 5, bottom). This result is very interesting as it shows that the rotation of gaseous planets (with a residual C₂₂ = 0) can also be chaotic when its orbit is perturbed by a companion planet.

In Fig. 2 we observe that a small stable synchronous island subsists at the middle of the chaotic zone. Therefore, we cannot exclude that after some chaotic wobble the spin finds a path into this stable region. Nevertheless, in our numerical experiments we never observed a simulation where the spin is permanently stabilized in the synchronous resonance. From time to time the rotation appears to enter the synchronous island, but then the C₂₂ grows to a value slightly higher than the theoretical estimation given by expression (29). When the maximum deformation is achieved the spin suddenly returns into the chaotic zone. Indeed, for C₂₂> 6.7 × 10^-6 the small resonant island might totally disappear and the spin might always remain chaotic.

5. Discussion

We show that close-in planets inside or close to orbital resonances undergo perturbations of their spins. For small eccentricities and weak orbital perturbations, the only spin equilibrium is the synchronous spin-orbit resonance, for which the rotation period equals the orbital revolution period. Tidal dissipation naturally drives the spin of the planet into this unique equilibrium, which is why close-in planets are usually assumed to be tidally synchronized.

When the planet-planet perturbations are significant, we demonstrate that new sub-synchronous and super-synchronous spin-orbit resonances appear, even for quasi-circular orbits. Moreover, for planets observed to transit, the TTVs provide the location (TTV period) and size (TTV amplitude) of these new resonances. For planets undergoing strong TTVs, the spin could be tidally locked in these asynchronous states or could even be chaotic.

We apply our modelling to KOI-227 b and Kepler-88 b, and run numerical simulations of the spin of these planets. We find that the spin of KOI-227 b has a non-negligible probability of being locked in an asynchronous resonance, while the spin of Kepler-88 b could be chaotic. In the case of KOI-227 b, we assume the planet to be mainly rocky since the bulk density found by Nesvorný et al. (2014) in the best fitting solution is very high (18.6 g cm^-3). However, the mass of the planet is not well constrained and this planet could have a non-negligible gaseous envelope. For Kepler-88 b (gaseous planet), in the most realistic cases, the spin is always locked in the (pseudo-)synchronous state. To observe a chaotic evolution we need the relaxation timescale (τ, see Sect. 4.2) to be at least 10 yr. This value is very typical for rocky planets, but probably overestimated for a gaseous planet such as Kepler-88 b. Nevertheless, these two cases illustrate very well what the spin evolution of smaller rocky planets with strong TTVs could be. We note that the planets with the strongest TTVs found in the literature (Nesvorný et al. 2012, 2013, 2014) are mostly giant planets, probably due to observational biases. Indeed, individual transits of small planets are very noisy, and these planets are usually fitted by phase-folding the light-curve. This makes small planets undergoing strong TTVs particularly challenging to detect.

For the sake of simplicity of the model we assume coplanar orbits, no obliquity, and low eccentricities for both planets. However, this is not a limitation for the application to observed systems, as numerical simulations including these effects can be performed. In particular, adding some inclination/obliquity will increase the number of degrees of freedom and probably ease the chaotic behaviour of the spin (Wisdom et al. 1984; Correia et al. 2015).

Acknowledgments

We thank the anonymous referee for his/her useful comments. We acknowledge financial support from SNSF and CIDMA strategic project UID/MAT/04106/2013. This work has been carried out in part within the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation.

References

Appendix A: Near-resonance case

In this appendix we show how the quasi-periodic decomposition of Λ_i, λ_i (see Eqs. (13), (15)) can be obtained in the near-resonance case. The non-resonant case is the easiest to deal with since the classical secular approximation can be used. The perturbative part of the Hamiltonian (see Eq. (8)) can be expanded in Fourier series of the mean longitudes ${\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{=}\sum _{\mathit{l}\mathrm{\in }{Z}^{\mathrm{2}}}{\mathit{C}}_{\mathit{l}}\mathrm{\left(}{\mathrm{\Lambda }}_{\mathit{i}}\mathit{,}{\mathit{D}}_{\mathit{i}}\mathit{,}{\mathit{\varpi }}_{\mathit{i}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathrm{(}{\mathit{l}}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{1}}\mathrm{+}{\mathit{l}}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}\mathrm{)}}\mathit{.}$(A.1)The secular evolution of the system is obtained by averaging this Hamiltonian over the fast angles λ_i (i = 1,2). This averaging transformation is obtained by a change of coordinates that is close to identity. We denote by ${\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}$, ${\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}$, ${\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}$, and ${\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}$ (i = 1,2) the new coordinates, by ℋ′ the new Hamiltonian, and by W the generating Hamiltonian of the transformation. To first order in the planet masses, the change of coordinates reads $\begin{array}{ccc}{\mathit{\lambda }}_{\mathit{i}}& \mathrm{=}& {\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\mathrm{\left\{}\mathit{W,}{\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}{\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathit{\partial W}}{\mathit{\partial }{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,}\\ {\mathrm{\Lambda }}_{\mathit{i}}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\mathrm{\left\{}\mathit{W,}{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathit{\partial W}}{\mathit{\partial }{\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,}\\ {\mathit{\varpi }}_{\mathit{i}}& \mathrm{=}& {\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\mathrm{\left\{}\mathit{W,}{\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}{\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathit{\partial W}}{\mathit{\partial }{\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,}\\ {\mathit{D}}_{\mathit{i}}& \mathrm{=}& {\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\mathrm{\left\{}\mathit{W,}{\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\right\}}\mathrm{=}{\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathit{\partial W}}{\mathit{\partial }{\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,}\end{array}$and the Hamiltonian ℋ′ is given by $\begin{array}{ccc}{\mathrm{\mathscr{H}}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{=}& {\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathit{,}\\ {\mathrm{\mathscr{H}}}_{\mathrm{1}}^{\mathrm{\prime }}& \mathrm{=}& {\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{+}\mathrm{\left\{}\mathit{W,}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{\right\}}\mathit{,}\end{array}$with ${\mathrm{\mathscr{H}}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}\mathrm{⟨}{\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{\mathrm{0}}\mathit{,}$(A.8)and where the braces denote the Poisson brackets. We thus have $\mathrm{\left\{}\mathit{W,}{\mathrm{\mathscr{H}}}_{\mathrm{0}}\mathrm{\right\}}\mathrm{=}\mathrm{⟨}{\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{⟩}\mathrm{-}{\mathrm{\mathscr{H}}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\sum _{\mathit{l}\mathrm{\in }{Z}^{\mathrm{2}}\mathrm{\backslash }\mathrm{0}}{\mathit{C}}_{\mathit{l}}{\mathrm{e}}^{\mathrm{i}\mathit{l}\mathit{.}\mathit{\lambda }}\mathit{,}$(A.9)which is called the homological equation and whose solution is $\mathit{W}\mathrm{=}\sum _{\mathit{l}\mathrm{\in }{Z}^{\mathrm{2}}\mathrm{\backslash }\mathrm{0}}\frac{{\mathit{C}}_{\mathit{l}}}{\mathrm{i}\mathit{l}\mathit{.}\mathit{n}}{\mathrm{e}}^{\mathrm{i}\mathit{l}\mathit{.}\mathit{\lambda }}\mathit{,}$(A.10)where n_i (i = 1,2) are the unperturbed Keplerian mean-motions of the planets.

The secular Hamiltonian ℋ′ no longer depends on the mean longitudes ${\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}$, which implies that the coordinates ${\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}$ are constants of motion. The long-term variations of ${\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}$ (i.e. eccentricities), ${\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}$, and ${\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}$, can be solved by using the Hamiltonian ℋ′. However, in this study, we are interested in the evolution of the system at shorter timescales and will neglect this secular evolution. We thus assume that ${\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}$ and ${\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}$ are constants and that ${\mathit{\lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{=}{\mathit{\lambda }}_{\mathit{i,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{n}}_{\mathit{i}}\mathit{t}$. In order to obtain the real evolution of the system, we need to revert to the original coordinate system using Eqs. (A.2)–(A.5). In particular, we have ${\mathrm{\Lambda }}_{\mathit{i}}\mathrm{=}{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{-}\sum _{\mathit{l}\mathrm{\in }{Z}^{\mathrm{2}}\mathrm{\backslash }\mathrm{0}}\frac{{\mathit{l}}_{\mathit{i}}{\mathit{C}}_{\mathit{l}}}{\mathit{l}\mathit{.}\mathit{n}}{\mathrm{e}}^{\mathrm{i}\mathit{l}\mathit{.}{\mathit{\lambda }}^{\mathrm{\prime }}}\mathit{.}$(A.11)Let us consider a system that is close to a k₂:k₁ resonance but outside of it. The combination ν = k₂n₂ − k₁n₁ is thus small (but not zero) which enhances the corresponding terms in the change of coordinates (small divisor). Keeping only this enhanced term, we have ${\mathrm{\Lambda }}_{\mathit{i}}\mathrm{\approx }{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{2}{\mathit{k}}_{\mathit{i}}}{\mathit{\nu }}\mathrm{\left|}{\mathit{C}}_{\mathit{k}}\mathrm{\right|}\cos \mathrm{(}\mathit{\nu t}\mathrm{+}{\mathit{\phi }}_{\mathit{i}}\mathrm{)}\mathit{,}$(A.12)with k = ( − k₁,k₂), C_k = ( − 1)^{i + 1} | C_k | e^iφ_i, and C_{− k} is its complex conjugate (by construction). The degree of the resonance k₂:k₁ is denoted q = k₂ − k₁. The coefficient C_k is of order q in eccentricity (D’Alembert rule), and can be written (e.g. Laskar & Robutel 1995) ${\mathit{C}}_{\mathit{k}}\mathrm{=}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{{\mathit{a}}_{\mathrm{2}}}{\mathit{e}}^{\mathit{q}}{\mathit{c}}_{\mathit{k}}\mathit{,}$(A.13)where e = max(e₁,e₂) and c_k is of the order of unity. Finally, we have (see Eq. (16)) $\begin{array}{ccc}\mathit{\nu }& \mathrm{=}& {\mathit{k}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{1}}\mathit{,}\\ \mathit{\alpha }& \mathrm{=}& \mathrm{6}{\mathit{k}}_{\mathit{i}}\mathrm{\left|}{\mathit{C}}_{\mathit{k}}\mathrm{\right|}\frac{{\mathit{n}}_{\mathit{i}}}{{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}{\mathit{\nu }}^{\mathrm{2}}}\mathrm{~}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}{\mathit{e}}^{\mathit{q}}{\left(\frac{{\mathit{n}}_{\mathit{i}}}{\mathit{\nu }}\right)}^{\mathrm{2}}\mathit{,}\end{array}$where m_p is the mass of the perturbing planet. We observe that when the system is very close to the resonance separatrix, ν ≪ n, the amplitude of oscillations increases (see Eq. (A.15)).

Appendix B: Resonant case

In this appendix we show how the quasi-periodic decomposition of Λ_i, λ_i (see Eq. (13)) can be obtained in the resonant case. The resonant case arises when the small divisor of Eq. (A.11) is too small and the averaging technique of Appendix A is no longer valid. However, the non-resonant terms can still be averaged out using the same procedure (see Appendix A). The resonant secular Hamiltonian ℋ′ is constructed such that resonant terms (of the form p(k₂λ₂ − k₁λ₁)) are kept, ${\mathrm{\mathscr{H}}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}\sum _{\mathit{p}\mathrm{\in }Z}{\mathit{C}}_{\mathrm{(}\mathrm{-}\mathit{p}{\mathit{k}}_{\mathrm{1}}\mathit{,p}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}{\mathrm{e}}^{\mathrm{i}\mathit{p}\mathrm{(}{\mathit{k}}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{)}}\mathit{.}$(B.1)The solution to the homological equation (Eq. (A.9)) is in this case $\mathit{W}\mathrm{=}\sum _{\mathit{l}\mathrm{\in }{Z}^{\mathrm{2}}\mathit{,l}\ne \mathrm{(}\mathrm{-}\mathit{p}{\mathit{k}}_{\mathrm{1}}\mathit{,p}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}\frac{{\mathit{C}}_{\mathit{l}}}{\mathrm{i}\mathit{l}\mathit{.}\mathit{n}}{\mathrm{e}}^{\mathrm{i}\mathit{l}\mathit{.}\mathit{\lambda }}\mathit{,}$(B.2)and the resonant Hamiltonian reads ${\mathrm{\mathscr{H}}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}\sum _{\mathit{p}\mathrm{\in }Z}{\mathit{C}}_{\mathrm{(}\mathrm{-}\mathit{p}{\mathit{k}}_{\mathrm{1}}\mathit{,p}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}{\mathrm{e}}^{\mathrm{i}\mathit{p}\mathrm{(}{\mathit{k}}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{)}}\mathit{,}$(B.3)where the coefficients C_l are functions of ${\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}$, ${\mathit{D}}_{\mathit{i}}^{\mathrm{\prime }}$, ${\mathit{\varpi }}_{\mathit{i}}^{\mathrm{\prime }}$. In the following we neglect the secular evolution of the eccentricities and longitudes of periastron as in the non-resonant case (see Appendix A). Moreover, the coefficients C_l have a weak dependency over Λ′ (e.g. Delisle et al. 2014), and in first approximation we assume they are constant. According to Eq. (A.13), we have ${\mathit{C}}_{\mathrm{(}\mathrm{-}\mathit{p}{\mathit{k}}_{\mathrm{1}}\mathit{,p}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}\mathrm{=}\frac{\mathrm{𝒢}{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{{\mathit{a}}_{\mathrm{2}}}{\mathit{e}}^{\mathrm{\left|}\mathit{p}\mathrm{\right|}\mathit{q}}{\mathit{c}}_{\mathrm{(}\mathrm{-}\mathit{p}{\mathit{k}}_{\mathrm{1}}\mathit{,p}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}\mathit{,}$(B.4)where q = k₂ − k₁ is the degree of the resonance, and c_{(− pk1,pk2)} is of the order of unity. In a first approximation, we only keep the terms of order q in eccentricity, and obtain the Hamiltonian ${\mathrm{\mathscr{H}}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\mathscr{H}}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}\mathrm{\left|}{\mathit{C}}_{\mathrm{-}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}}\mathrm{\right|}\cos \mathrm{(}{\mathit{k}}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{+}\mathit{\phi }\mathrm{)}\mathit{,}$(B.5)where C₀ has been dropped since we assumed it to be constant. We introduce the canonical change of coordinates

$\begin{array}{ccc}& & \mathit{\psi }\mathrm{=}{\mathit{k}}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{+}\mathit{\phi ,} \mathit{J}\mathrm{=}\frac{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}{\mathrm{\Lambda }}_{\mathrm{1}}^{\mathrm{\prime }}}{{\mathit{k}}_{\mathrm{1}}}\mathit{,}\\ & & \mathit{\xi }\mathrm{=}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathit{,} \mathrm{\Gamma }\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{2}}}{{\mathit{k}}_{\mathrm{1}}}{\mathrm{\Lambda }}_{\mathrm{1}}^{\mathrm{\prime }}\mathit{.}\end{array}$Since ℋ′ only depends on the angle ψ and not on ξ, Γ is a conserved quantity. We now develop the Keplerian part ${\mathrm{\mathscr{H}}}_{\mathrm{0}}^{\mathrm{\prime }}$ in power series of J. Constant terms can be neglected since they do not contribute to the dynamics. Moreover, the first order terms cancel out for a resonant system. We thus obtain, up to second order in J, ${\mathrm{\mathscr{H}}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\approx }\mathrm{-}{\mathit{K}}_{\mathrm{2}}{\mathit{J}}^{\mathrm{2}}\mathit{,}$(B.8)with ${\mathit{K}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\left({\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}\frac{{\mathit{n}}_{\mathrm{2}}}{{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}\right)\mathrm{·}$(B.9)With this approximation, the Hamiltonian (B.5) reads ${\mathrm{\mathscr{H}}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{K}}_{\mathrm{2}}{\mathit{J}}^{\mathrm{2}}\mathrm{+}\mathrm{\left|}{\mathit{C}}_{\mathrm{(}\mathrm{-}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathrm{)}}\mathrm{\right|}\cos \mathit{\psi ,}$(B.10)which is the Hamiltonian of a simple pendulum. The frequency of libration is given by $\mathit{\nu }\mathrm{=}\sqrt{\mathrm{2}{\mathit{K}}_{\mathrm{2}}\mathrm{\left|}{\mathit{C}}_{\mathit{k}}\mathrm{\right|}}\frac{\mathit{\pi }}{\mathrm{2}}\mathit{K}\left(\sin \left(\frac{{\mathit{\psi }}_{\max }}{\mathrm{2}}\right)\right)\mathit{,}$(B.11)where K is the complete elliptical integral of the first kind, and the amplitude of libration ψ_max is in the range [ 0,π ]. For small amplitude oscillations, we have $\begin{array}{ccc}\mathit{\nu }& \mathrm{=}& \sqrt{\mathrm{2}{\mathit{K}}_{\mathrm{2}}\mathrm{\left|}{\mathit{C}}_{\mathit{k}}\mathrm{\right|}}\mathrm{~}\mathit{n}\sqrt{\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}{\mathit{e}}^{\mathit{q}}}\mathit{,}\\ \mathit{\psi }& \mathrm{=}& {\mathit{\psi }}_{\max }\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\zeta }\mathrm{)}\mathit{,}\\ \mathit{J}& \mathrm{=}& \mathrm{-}\frac{\mathit{\nu }{\mathit{\psi }}_{\max }}{\mathrm{2}{\mathit{K}}_{\mathrm{2}}}\cos \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\zeta }\mathrm{)}\mathit{,}\end{array}$(B.12)and from Eqs. (B.6), (B.7), we obtain $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{1}}^{\mathrm{\prime }}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{1}}\mathit{J}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{k}}_{\mathrm{1}}\frac{\mathit{\nu }{\mathit{\psi }}_{\max }}{\mathrm{2}{\mathit{K}}_{\mathrm{2}}}\cos \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\zeta }\mathrm{)}\mathit{,}\\ {\mathrm{\Lambda }}_{\mathrm{2}}^{\mathrm{\prime }}& \mathrm{=}& {\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathit{J}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}{\mathit{k}}_{\mathrm{2}}\frac{\mathit{\nu }{\mathit{\psi }}_{\max }}{\mathrm{2}{\mathit{K}}_{\mathrm{2}}}\cos \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\zeta }\mathrm{)}\mathit{.}\end{array}$In principle, we should revert to the original coordinate system (coordinates without primes), as in the non-resonant case (see Appendix A). However, in the resonant case, there are no enhanced terms (with small divisors) in the change of coordinates W, since we kept these terms in the secular Hamiltonian ℋ′. We thus have ${\mathrm{\Lambda }}_{\mathit{i}}\mathrm{\approx }{\mathrm{\Lambda }}_{\mathit{i}}^{\mathrm{\prime }}$. Finally, from Eqs. (16), (B.13), (B.14), we deduce $\begin{array}{ccc}\mathit{\alpha }& \mathrm{=}& \frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{k}}_{\mathit{i}}{\mathit{n}}_{\mathit{i}}}{{\mathit{K}}_{\mathrm{2}}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}{\mathit{\psi }}_{\max }\\ & \mathrm{=}& \frac{\mathrm{1}\mathit{/}{\mathrm{\Lambda }}_{\mathit{i,}\mathrm{0}}}{{\mathit{k}}_{\mathrm{1}}\mathit{/}{\mathrm{\Lambda }}_{\mathrm{1}\mathit{,}\mathrm{0}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathit{/}{\mathrm{\Lambda }}_{\mathrm{2}\mathit{,}\mathrm{0}}}{\mathit{\psi }}_{\max }\\ & \mathrm{\approx }& \frac{{\mathit{n}}_{\mathit{i}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{/}{\mathit{m}}_{\mathit{i}}}{{\mathit{k}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{1}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{/}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{2}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{/}{\mathit{m}}_{\mathrm{2}}}{\mathit{\psi }}_{\max }\mathit{.}\end{array}$(B.15)In particular, for m₁ ≪ m₂, we have α = ψ_max/k₁ for the inner planet (1), and α = 0 for the outer planet (2). On the contrary, for m₁ ≫ m₂, we have α = 0 for the inner planet, and α = ψ_max/k₂ for the outer one. This is not surprising as the less massive planet undergoes the strongest perturbations.

Appendix C: Permanent and tidally induced deformation

In this appendix we show how to estimate the width of the synchronous resonance (σ) from the known properties of the planets. From Eq. (6), we have $\frac{\mathit{\sigma }}{\mathit{n}}\mathrm{=}\sqrt{\mathrm{12}\hspace{0.17em}{\mathit{C}}_{\mathrm{22}}\mathit{/}\mathit{\zeta }}\mathit{.}$(C.1)The internal structure factor ζ can be estimated from k_f through the Darwin-Radau equation (e.g. Jeffreys 1976) $\mathit{\zeta }\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathrm{5}}\sqrt{\frac{\mathrm{4}\mathrm{-}{\mathit{k}}_{\mathrm{f}}}{\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathrm{f}}}}\right)\mathrm{·}$(C.2)We have ζ = 2 / 5 (k_f = 3 / 2) for a homogenous sphere, ζ ≈ 1 / 3 for rocky planets, and ζ ≈ 1 / 4 for gaseous planets.

The deformation of the planet (C₂₂) has two components, the permanent deformation (C_22,r) due to the intrinsic mass repartition in the planet and the tidally induced deformation (C_22,t). The permanent asymmetry of the mass repartition can be roughly estimated from the mass and radius of a rocky planet, and using observations of the solar system planets (see Yoder 1995) ${\mathit{C}}_{\mathrm{22}\mathit{,r}}\mathrm{~}{\mathrm{10}}^{-6}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{\oplus }}}\right)}^{\mathrm{5}}{\left(\frac{\mathit{m}}{{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{·}$(C.3)For gaseous planets, the permanent asymmetry is very weak (C_22,r ≈ 0).

The tidally induced deformation corresponds to the adjustment of the planet’s mass distribution to the external gravitational potential. This deformation is not instantaneous, and the relaxation time (τ, see Sect. 4.2) depends on the planet’s internal structure. If the deformation were instantaneous, the C_22,t coefficient would be given by (Eq. (35) with τ = 0, see also Correia & Rodríguez 2013) ${\mathit{C}}_{\mathrm{22}\mathit{,t,}\mathrm{inst}\mathrm{.}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\cos \mathrm{2}\mathrm{(}\mathit{\theta }\mathrm{-}\mathit{f}\mathrm{)}\mathit{,}$(C.4)where we assume the obliquity to be negligible. This expression is very similar to Eq. (1), and using the same approximations, we obtain an expression very similar to Eq. (20): $\begin{array}{ccc}{\mathit{C}}_{\mathrm{22}\mathit{,t,}\mathrm{inst}\mathrm{.}}& \mathrm{=}& {\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}[\cos \mathrm{2}\mathit{\gamma }\\ & & \mathrm{+}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\right)\cos \mathrm{2}\left(\mathit{\gamma }\mathrm{+}\frac{{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}}{\mathrm{2}}\right)\\ & & \mathrm{+}\sum _{\mathit{j}}{\mathit{\alpha }}_{\mathit{j}}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\right)\cos \mathrm{2}\left(\mathit{\gamma }\mathrm{-}\frac{{\mathit{\nu }}_{\mathit{j}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{j}}\mathrm{+}\mathit{\pi }}{\mathrm{2}}\right)]\end{array}$(C.5)with ${\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{{\mathit{a}}_{\mathrm{0}}}\right)}^{\mathrm{3}}\mathrm{·}$(C.6)We assume in the following that the relaxation time τ is much longer that the variations of C_22,t,inst., such that the actual C_22,t of the planet is the mean value of C_22,t,inst. (see Eq. (35)): ${\mathit{C}}_{\mathrm{22}\mathit{,t}}\mathrm{=}⟨{\mathit{C}}_{\mathrm{22}\mathit{,}\hspace{0.17em}\mathit{t,}\hspace{0.17em}\mathrm{inst}\mathrm{.}}⟩\mathit{.}$(C.7)From Eq. (C.5), we deduce that as long as the spin is outside of any resonance, the tidal deformation average out (C_22,t = 0). Therefore, before the resonant capture, the C₂₂ coefficient of the planet reduces to its permanent deformation (C_22,r).

If the spin is locked in the synchronous resonance, the averaged tidally induced deformation reaches ${\mathit{C}}_{\mathrm{22}\mathit{,t}}\mathrm{=}{\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$(C.8)when the amplitude of libration is small (sin2γ ≈ 0, cos2γ ≈ 1). Similarly, if the spin is locked in a sub/super-synchronous resonance, we have ${\mathit{C}}_{\mathrm{22}\mathit{,t}}\mathrm{=}{\mathit{\alpha }}_{\mathit{j}}\left(\mathrm{1}\mathrm{\pm }\frac{{\mathit{\nu }}_{\mathit{j}}}{\mathit{n}}\right){\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$(C.9)at the centre of the resonance.

Let us apply this reasoning to KOI-227 b (rocky planet). We obtain C_22,r = 1.4 × 10^-7 for the permanent deformation (see Eq. (C.3)). The tidal deformation is C_22,t = 2.6 × 10^-6 at the centre of the synchronous resonance (see Eq. (C.8)), and C_22,t = 7.9 × 10^-7 at the centre of the sub/super-synchronous resonance (see Eq. (C.9)). The maximum deformation (C₂₂ = C_22,r + C_22,t) is thus 2.8 × 10^-6 in the synchronous case, and 9.3 × 10^-7 in the sub/super-synchronous case.

In the case of Kepler-88 b (gaseous planet) we obtain C_22,r = 0, C₂₂ = C_22,t = 2.0 × 10^-5 for the synchronous resonance, and C₂₂ = C_22,t = 5.6 × 10^-6 for the sub/super-synchronous resonances.

We note that all these values assume that the amplitude of libration in the resonance is vanishing. We show in Appendix D that forced oscillations are non-negligible and significantly reduce the tidal deformation.

Appendix D: Forced oscillations and tidal deformation

In this appendix we derive the amplitude of forced oscillations in the synchronous and sub/super-synchronous resonances, as well as their implications on the tidal deformation of the planet. We assume that a single term is dominating the planet’s TTVs, with amplitude α, and frequency ν ≪ n. We neglect here the semi-major axis variations since their contribution to the spin evolution is of order ν/n compared to the contribution of the mean longitude variations (see Eqs. (15), (19), and (20)).

Appendix D.1: Synchronous case

We first assume that the spin of the planet is locked in the synchronous resonance. From Eqs. (3) and (15) we deduce $\begin{array}{c}\mathrm{¨}\\ \mathit{\gamma }\end{array}\mathrm{\approx }\mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}\sin \mathrm{2}(\mathit{\gamma }\mathrm{+}\mathit{\alpha }\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)})\mathit{,}$(D.1)where (see Eq. (C.1)) ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}\mathrm{12}{\mathit{C}}_{\mathrm{22}}{\mathit{n}}^{\mathrm{2}}\mathit{/}\mathit{\zeta }$(D.2)and (see Eq. (C.4)) ${\mathit{C}}_{\mathrm{22}}\mathrm{=}{\mathit{C}}_{\mathrm{22}\mathit{,r}}\mathrm{+}{\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}⟨\cos \mathrm{2}(\mathit{\gamma }\mathrm{+}\mathit{\alpha }\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)})⟩\mathit{.}$(D.3)We introduce h = γ + αsin(νt + φ), such that $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{h}\end{array}& \mathrm{=}& \mathrm{-}\mathit{\alpha }{\mathit{\nu }}^{\mathrm{2}}\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathrm{-}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}\sin \mathrm{2}\mathit{h,}\\ {\mathit{\sigma }}^{\mathrm{2}}& \mathrm{=}& {\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{t,}\mathrm{0}}^{\mathrm{2}}⟨\cos \mathrm{2}\mathit{h}⟩\mathit{,}\end{array}$with $\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{12}{\mathit{C}}_{\mathrm{22}\mathit{,r}}{\mathit{n}}^{\mathrm{2}}\mathit{/}\mathit{\zeta ,}\\ & & {\mathit{\sigma }}_{\mathit{t,}\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathrm{12}{\mathit{C}}_{\mathrm{22}\mathit{,t}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{\mathit{n}}^{\mathrm{2}}\mathit{/}\mathit{\zeta }\mathit{.}\end{array}$Equations (D.4), (D.5) can be developed in power series of h$\begin{array}{ccc}& & \begin{array}{c}\mathrm{¨}\\ \mathit{h}\end{array}\mathrm{=}\mathrm{-}\mathit{\alpha }{\mathit{\nu }}^{\mathrm{2}}\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathrm{-}{\mathit{\sigma }}^{\mathrm{2}}\left(\mathit{h}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{h}}^{\mathrm{3}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{15}}{\mathit{h}}^{\mathrm{5}}\mathrm{+}\mathit{...}\right)\mathit{,}\\ & & {\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{t,}\mathrm{0}}^{\mathrm{2}}⟨\mathrm{1}\mathrm{-}\mathrm{2}{\mathit{h}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{h}}^{\mathrm{4}}\mathrm{+}\mathit{...}⟩\mathit{.}\end{array}$At first order (linearized equation), the forced solution is simply${\mathit{h}}_{\mathrm{lin}\mathrm{.}}\mathrm{=}\frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\nu }}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathit{t,}\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{\alpha }\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathit{,}$(D.10)which is of order α. At higher orders, odd harmonics of the forced frequency ν appears, and h takes the form $\mathit{h}\mathrm{=}\sum _{\mathit{k}\mathrm{\in }\mathrm{2}N\mathrm{+}\mathrm{1}}{\mathit{h}}_{\mathit{k}}\sin \mathit{k}\mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathit{,}$(D.11)where h_k is of order α^k. We replace Eq. (D.11) in Eqs. (D.8) and (D.9), and truncate the resulting expression at a given order N in α. By identifying terms of frequency kν (k ∈ [ 1,N ]), we obtain N polynomial equations on the coefficients h₁,...,h_N. Solving this set of equations allows us to determine the coefficients h_i, as well as the corresponding C₂₂ value. The forced oscillations of γ are then given by ${\mathit{\gamma }}_{\mathrm{forced}}\mathrm{=}\sum _{\mathit{k}\mathrm{\in }\mathrm{2}N\mathrm{+}\mathrm{1}}{\mathit{\gamma }}_{\mathit{k}}\sin \mathit{k}\mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathit{,}$(D.12)with γ₁ = h₁ − α, γ_k = h_k (k> 1).

Applying this reasoning to KOI-227 b, we obtain C_22,t = 1.8 × 10^-6 (instead of 2.6 × 10^-6) and C₂₂ = C_22,t + C_22,r = 1.9 × 10^-6 (instead of 2.8 × 10^-6). The forced amplitude at frequency ν is γ₁ = 17°, while the amplitudes of harmonics (3ν, etc.) decrease rapidly (γ_{k + 2}/γ_k ~ 10^-2). In the case of Kepler-88 b, we obtain C₂₂ = C_22,t = 6.7 × 10^-6 (instead of 2.0 × 10^-5) and a forced amplitude of 36° (at frequency ν). As was true for KOI-227 b, for Kepler-88 b the amplitudes of harmonics decrease rapidly.

Appendix D.2: Sub/super-synchronous case

We now consider the same effect (forced oscillations at frequency ν and its harmonics) but for the sub/super-synchronous resonances. We assume that the spin is locked in one of these resonances, such that $\mathit{\gamma }\mathrm{=}\mathrm{\pm }\frac{\mathit{\nu t}\mathrm{+}{\mathit{\phi }}_{\mathrm{\pm }}}{\mathrm{2}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{forced}}\mathit{,}$(D.13)with φ₊ = φ + π, φ₋ = φ. As for the synchronous resonance, we introduce h = γ_forced + αsin(νt + φ) such that $\begin{array}{ccc}& & \begin{array}{c}\mathrm{¨}\\ \mathit{h}\end{array}\mathrm{=}\mathrm{-}\mathit{\alpha }{\mathit{\nu }}^{\mathrm{2}}\sin \mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathrm{\mp }\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{2}}\sin \mathrm{(}\mathrm{2}\mathit{h}\mathrm{\pm }\mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{\left)}\mathrm{\right)}\mathit{,}\\ & & {\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\pm }{\mathit{\sigma }}_{\mathit{t,}\mathrm{0}}^{\mathrm{2}}⟨\cos \mathrm{(}\mathrm{2}\mathit{h}\mathrm{\pm }\mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{\left)}\mathrm{\right)}⟩\mathit{.}\end{array}$As we did for the synchronous case, these expressions can be developed in power series of h. Replacing h by $\mathit{h}\mathrm{=}\sum _{\mathit{k}\mathrm{\in }\mathrm{2}N\mathrm{+}\mathrm{1}}{\mathit{h}}_{\mathit{k}}\sin \mathit{k}\mathrm{(}\mathit{\nu t}\mathrm{+}\mathit{\phi }\mathrm{)}\mathit{,}$(D.16)and truncating at a given order N in α (h_k being of order α^k), we obtain a set of N polynomial equations on h₁, ..., h_N. We then solve for h_i and determine the corresponding C₂₂ value.

In the case of KOI-227 b, we obtain C_22,t = 5.2 × 10^-7 (instead of 7.9 × 10^-7) and C₂₂ = 6.6 × 10^-7 (instead of 9.3 × 10^-7), with a forced amplitude of γ₁ = − 6° (at frequency ν). The negative sign means that the forced oscillations are dephased by an angle π with respect to the TTV signal. In the case of Kepler-88 b, we find C₂₂ = C_22,t = 2.2 × 10^-6 (instead of 5.6 × 10^-6), with a forced amplitude of − 10°. As in the case of the synchronous resonance, the amplitudes of harmonics (3ν, etc.) decrease rapidly for both planets (γ_{k + 2}/γ_k ~ 10^-3).

All Tables

All Figures

Fig. 1 Spin dynamics of KOI-227 b in the conservative case, with C22 = 1.4 × 10-7 (top, permanent deformation), 6.6 × 10-7 (middle, maximum deformation in the sub/super-synchronous resonances), and 1.9 × 10-6 (bottom, maximum deformation in the synchronous resonance). The left column shows the main frequency η of θ(t) for different initial values of (and with γ(0) = 0). The colour gives the derivative . Blue dots correspond to resonant motion, green dots to non-resonant regular motion, and red dots to chaotic motion. The vertical black line corresponds to the synchronization (η = n). The two grey lines correspond to the main sub/super-synchronous resonances (η = n ± ν/ 2). The dashed and dotted grey lines correspond to higher order sub/super-synchronous resonances (η = n + kν/ 2, k = ± 2, ± 3). The right column shows the same colour index (), but both γ(0) and are varied (2d maps). The vertical white line highlights the initial conditions taken in the left column (γ(0) = 0).

In the text

	Fig. 2 Same as Fig. 1 but for Kepler-88 b. We neglect the permanent deformation of the planet (C22,r = 0), and show the spin dynamics for C22 = 2.2 × 10-6 (top, maximum deformation in the sub/super-synchronous resonances), and 6.7 × 10-6 (bottom, maximum deformation in the synchronous resonance).
In the text

	Fig. 3 Examples of the spin evolution of KOI-227 b and corresponding global C22 value for different τ values. We use the initial conditions from Table 2, and the initial rotation period is set at 15.5 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.
In the text

	Fig. 4 Different final spin evolution of KOI-227 b and corresponding global C22 value for τ = 104 yr. We use the initial conditions from Table 2, and the initial rotation period is set at 15.5 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.
In the text

	Fig. 5 Examples of the spin evolution of Kepler-88 b and corresponding global C22 value for different τ values. We use the initial conditions from Table 3, and the initial rotation period is set at 10 d. The green lines give the position of the spin-orbit resonances n, n ± ν/ 2, and n ± ν.
In the text