Spin dynamics of close-in planets exhibiting large TTVs

TL;DR

This paper investigates how gravitational interactions in multi-planet systems can lead to non-synchronous and chaotic spin states in close-in planets, using TTVs as a diagnostic tool.

Contribution

It demonstrates that planet-planet perturbations can significantly alter the spin states of close-in planets, challenging the common assumption of synchronous rotation.

Findings

KOI-227b can have non-synchronous rotation.

Kepler-88b's rotation can become chaotic.

TTVs effectively probe spin dynamics.

Abstract

We study the spin evolution of close-in planets in compact multi-planetary systems. The rotation period of these planets is often assumed to be synchronous with the orbital period due to tidal dissipation. Here we show that planet-planet perturbations can drive the spin of these planets into non-synchronous or even chaotic states. In particular, we show that the transit timing variation (TTV) is a very good probe to study the spin dynamics, since both are dominated by the perturbations of the mean longitude of the planet. We apply our model to KOI-227b and Kepler-88b, which are both observed undergoing strong TTVs. We also perform numerical simulations of the spin evolution of these two planets. We show that for KOI-227b non-synchronous rotation is possible, while for Kepler-88b the rotation can be chaotic.

Tables

Table 1. Table 1: Qualitative evolution of the spin as a function of the amplitude α 𝛼 \alpha and frequency ν 𝜈 \nu of the perturbation.

	$\sigma /\nu \ll 1$	$\sigma /\nu \sim 1$	$\sigma /\nu \gg 1$
$\alpha \ll 1$	synchronous resonance only
$\alpha \sim 1$	3 separated res.	chaos	modulated pendulum

Table 2. Table 2: Parameters for KOI-227   b 𝑏 b , c 𝑐 c used in this study.

Parameter	[unit]	$b$	$c$
$m$	[$M_{\oplus }$]	$11.09$	$43.89$
$R$	[$R_{\oplus }$]	$2.23$	–
$k_f$		$0.95$	–
$\zeta$		$0.333$	–
$C_{22,r}$		$1.4\times {10}^{-7}$	–
$a$	[AU]	$0.104483724$	$0.167413850$
$e$		$0.0756724$	$0.0173317$
$\lambda$	[deg]	$0$	$250.4696$
$\varpi$	[deg]	$-179.2092$	$119.1518$

Table 3. Table 3: Parameters for Kepler-88   b 𝑏 b , c 𝑐 c used in this study.

Parameter	[unit]	$b$	$c$
$m$	[$M_{\oplus }$]	$8.7$	$198.8$
$R$	[$R_{\oplus }$]	$3.78$	–
$k_f$		$0.45$	–
$\zeta$		$0.25$	–
$C_{22,r}$		0	–
$a$	[AU]	$0.095093133$	$0.152955525$
$e$		$0.05593$	$0.05628$
$i$	[deg]	$0.945$	$3.8$
$\lambda$	[deg]	$6.405$	$252.9$
$\varpi$	[deg]	$90.59$	$270.76$
$\Omega$	[deg]	$270$	$264.1$

Table 4. Table 4: Capture probabilities in spin-orbit resonances (in percent) for KOI-227   b 𝑏 b , using different τ 𝜏 \tau values.

$\eta -n$	$\tau$ (yr)
	${10}^1$	${10}^2$	${10}^3$	${10}^4$
$\nu$	$-$	$-$	3.8	2.3
$\nu /2$	$-$	57.3	74.2	31.4
$0$	100.0	42.7	22.0	56.7
$-\nu /2$	$-$	$-$	$-$	9.6

Equations

$$\ddot{\theta }=-\frac{6C_{22}}{\zeta }\frac{\mu }{r^3}\sin 2(\theta -f),$$

$$\ddot{\theta }=-\frac{6C_{22}}{\zeta }\frac{\mu }{a^3}\sum _{k\in \mathbb{Z}}{X}_k^{-3,2}(e)\sin (2\theta -k\lambda +(k-2)\varpi ),$$

$$\ddot{\theta }\approx -\frac{6C_{22}}{\zeta }\frac{\mu }{a^3}\sin 2(\theta -\lambda ).$$

$$\lambda ={\lambda }_0+nt,$$

$$\gamma$$

$$\sigma$$

$$\ddot{\gamma }\approx -\frac{{\sigma }^2}{2}\sin 2\gamma ,$$

$$\mathcal{H}={\mathcal{H}}_0({\Lambda }_i)+{\mathcal{H}}_1({\Lambda }_i,{\lambda }_i,D_i,{\varpi }_i),$$

$${\mathcal{H}}_0=-\frac{\mathcal{G}{m}_0{m}_1}{2a_1}-\frac{\mathcal{G}{m}_0{m}_2}{2a_2}=-\frac{{\mu }_1^2{\beta }_1^3}{2{\Lambda }_1^2}-\frac{{\mu }_2^2{\beta }_2^3}{2{\Lambda }_2^2},$$

$${\dot{\lambda }}_i$$

$${\dot{\Lambda }}_i$$

$${\dot{\lambda }}_i\approx \frac{{\mu }_i^2{\beta }_i^3}{{\Lambda }_i^3}\approx n_i{\left(\frac{{\Lambda }_i}{{\Lambda }_{i,0}}\right)}^{-3}\approx n_i-3n_i\frac{{\Lambda }_i-{\Lambda }_{i,0}}{{\Lambda }_{i,0}}.$$

$$\Lambda ={\Lambda }_0\left(1+\sum _j{A}_j\cos ({\nu }_{j}t+{\varphi }_j)\right),$$

$$\dot{\lambda }\approx n-3n\sum _j{A}_j\cos ({\nu }_{j}t+{\varphi }_j).$$

$$\lambda$$

$${\alpha }_j=3A_{j}n/{\nu }_j$$

$$\sin 2(\theta -\lambda )\approx \sin 2\gamma +2\sum _j{\alpha }_j\sin ({\nu }_{j}t+{\varphi }_j)\cos 2\gamma ,$$

$$\gamma =\theta -\left({\lambda }_0+nt+\sum _j{\alpha }_j\sin {\varphi }_j\right).$$

$$\frac{1}{a^3}$$

$$\displaystyle\ddot{\gamma}\approx-\frac{\sigma^{2}}{2}\Bigg{[}\!\!\!\!$$

$${\mathbf{e}}^{\mathrm{i}f}=\sum _{k\in \mathbb{Z}}{X}_k^{0,1}(e){\mathbf{e}}^{\mathrm{i}(k\lambda +(1-k)\varpi )}$$

$$f\approx {\lambda }_0+nt-\sum _j{\alpha }_j\left(\sin ({\nu }_{j}t+{\varphi }_j)-\sin {\varphi }_j\right)$$

$${\mathrm{TTV}}_k=t_k-kP\approx \sum _j\frac{{\alpha }_j}{n}\left(\sin ({\nu }_j{t}_k+{\varphi }_j)-\sin {\varphi }_j\right)-\frac{{\lambda }_0}{n},$$

$${\mathrm{TTV}}_k\approx \frac{\alpha }{n}\left(\sin (\nu t_k+\varphi )-\sin (\varphi )\right)-\frac{{\lambda }_0}{n},$$

$$\displaystyle\ddot{\gamma}\approx-\frac{\sigma^{2}}{2}\bigg{[}\sin 2\gamma\!\!\!\!$$

$$C_{22,r}\approx 1.4\times 10^{-7}.$$

$$C_{22,\text{sync.}}\approx 1.9\times 10^{-6},$$

$$C_{22,\text{sub/super}}\approx 6.6\times 10^{-7}.$$

$$C_{22,\text{sync.}}$$

$$C_{22,\text{sub/super}}$$

$${\ddot{\text{@vec}r}}_1$$

$${\ddot{\text{@vec}r}}_2$$

Full text

1 Introduction

The rotation of close-in planets is usually modified by tidal interactions with the central star and reaches a stationary value on time-scales 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 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_{0}$, and two companion planets with masses $m_{1}$ and $m_{2}$, such that $m_{1},m_{2}\ll m_{0}$. 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) 111Tidal effects drive the obliquity of the planets to zero degrees (e.g. Correia et al., 2003; Boué et al., 2016), so we expect that Kepler planets whose spin has been driven close to the synchronous rotation also have nearly zero obliquity..

We introduce $\theta$, the rotation angle of the planet with respect to an inertial line, whose evolution is described by (e.g. Murray & Dermott, 1999)

[TABLE]

where $\mu=\mathcal{G}(m_{0}+m)$ ($\mathcal{G}$ 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_{22}$ is the Stokes gravity field coefficient that measures the asymmetry in the equatorial axes of the planet, and $\zeta$ 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

[TABLE]

where $a$, $e$, $\lambda$, $\varpi$ 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 $X_{k}^{-3,2}(e)$ is of order $|k-2|$ in eccentricity. Thus, the leading term in this expansion corresponds to $k=2$, and we have $X_{2}^{-3,2}(e)=1+O(e^{2})$. If we only keep this leading term, Eq. (2) simplifies as

[TABLE]

If the planet remains unperturbed and follows a Keplerian orbit, the semi-major axis is constant and the mean longitude is given by

[TABLE]

where $n$ is the constant mean motion of the planet. Introducing

[TABLE]

we can rewrite Eq. (3) as

[TABLE]

which is the equation of a simple pendulum, with a stable equilibrium point at $\gamma=0$. This equilibrium corresponds to the exact synchronization ($\dot{\theta}=n$). The parameter $\sigma$ (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 $\lambda_{i}$. The action canonically conjugated with $\lambda_{i}$ is the circular angular momentum of planet $i$, $\Lambda_{i}=\beta_{i}\sqrt{\mu_{i}a_{i}}$, where $\beta_{i}=m_{i}m_{0}/(m_{0}+m_{i})$. We also introduce the angular momentum deficit (AMD) of the planets, $D_{i}=\Lambda_{i}\left(1-\sqrt{1-e_{i}^{2}}\right)$ which are the actions canonically conjugated with $\varpi_{i}$. The Hamiltonian of the three-body problem reads

[TABLE]

where $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ are the Keplerian and perturbative part of the Hamiltonian. We have

[TABLE]

and $\mathcal{H}_{1}$ is of order one in planet masses (over star mass). The equations of motion read

[TABLE]

However, it can be shown that the perturbative part has a weak dependency on $\Lambda_{i}$ (e.g. Delisle et al., 2014). At leading order in eccentricity, we have

[TABLE]

We assume that the orbital motion is not chaotic and thus quasi-periodic. We introduce the quasi-periodic decomposition of $\Lambda$

[TABLE]

where $\nu_{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_{2}$:$k_{1}$) but outside of it (see Appendix A), the dominant frequency is the frequency of circulation $k_{2}n_{2}-k_{1}n_{1}$. 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

[TABLE]

We thus have

[TABLE]

where the coefficients

[TABLE]

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^{\circ}$. With this approximation, we have at first order in $\alpha_{j}$

[TABLE]

with

[TABLE]

From Eq. (13) we deduce

[TABLE]

and finally, using the expressions of Eqs. (17), (19) in Eq. (3), we obtain

[TABLE]

The term ($\sin 2\gamma$) 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 $\alpha_{j}$) in Eq. (2), corresponding to a sub-synchronous resonance ($\dot{\theta}=n-\nu_{j}/2$), and a super-synchronous resonance ($\dot{\theta}=n+\nu_{j}/2$). 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 ($\dot{\theta}=kn/2$) appear, as do new resonances of the form $\dot{\theta}=kn/2\pm\nu_{j}/2$. Moreover, if the amplitudes $\alpha_{j}$ are not small, the series should be developed at a higher degree in $\alpha_{j}$, and resonances of the type $\dot{\theta}=kn/2\pm l\nu_{j}/2$ would appear (see Leleu et al., 2015, for the co-orbital case).

Different dynamical regimes can be observed depending on the values of $\alpha$, $\nu$, and $\sigma$. 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 2 describes qualitatively the dynamics of the spin of the planet, as a function of the amplitude ($\alpha$) and frequency ($\nu$) 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

[TABLE]

where the Hansen coefficient $X_{k}^{0,1}$ is of order $|k-1|$ in eccentricity. For the sake of simplicity, we only keep the leading order term $f\approx\lambda$. We thus have (see Eq. (15))

[TABLE]

We introduce $t_{k}$, the time of the $k$-th transit. We have $f(t_{k})=k2\pi$, thus the transit timing variations are given by

[TABLE]

where $P=2\pi/n$ is the orbital period. Therefore, the quasi-periodic decomposition of the TTV signal directly provides the amplitudes $\alpha_{j}$ and frequencies $\nu_{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))

[TABLE]

and (see Eq. (2))

[TABLE]

where we assume ($\nu\ll n$). The most interesting systems for our study (see Table 2) are those for which $\alpha$ is not negligible, such that the width of the resonances at $\dot{\theta}=n\pm\nu/2$ 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 $\alpha$ and $\nu$ (and of the phase $\phi$).

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 case 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_{22}\neq 0$ (Eq. (6)), while Kepler-88 $b$ is a gaseous planet for which the $C_{22}$ value is likely very close to zero (Campbell & Synnott, 1985).

4.1 Conservative evolution

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

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_{\oplus}$, a period of about 18 d, and TTVs with an amplitude of at least 10 hr have been observed (see Nesvorný et al., 2014). In terms of angular amplitude (see Eq. (15)), this corresponds to $\alpha\gtrsim 8^{\circ}$. The main TTV period is about 4.5 yr, thus $\nu/n\approx 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 ($\chi^{2}=37.6$). The mass of KOI-227 $b$ is $37.5\ M_{\oplus}$ for this solution, and its density would thus be $18.6\ \mathrm{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_{\oplus}$. The obtained solution has a $\chi^{2}$ of 39.1, which is still better than the best-fitting solution in other resonances ($\chi^{2}=51.5$ for the outer 3:2 resonance, and $\chi^{2}=82.2$ for the inner 3:2 resonance, see Nesvorný et al., 2014).

For this illustration, we adopted the mass of $11.09\ M_{\oplus}$ 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 ($\chi^{2}=53.3$). The obtained solution is given in Table 3.

Using the mass ($11.09\ M_{\oplus}$) and radius ($2.23\ R_{\oplus}$) of the planet, we estimate its permanent deformation (see Appendix C)

[TABLE]

This corresponds to $\sigma/n\approx 2.2\times 10^{-3}$ and $\sigma/\nu\approx 0.20$, which means that the sub/super-synchronous resonances ($\dot{\theta}=n\pm\nu/2$) 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_{22}$ increases due to the tidal deformation. We estimate the maximum deformation of the planet (see Appendices C and D) in the synchronous resonance

[TABLE]

and in the sub/super-synchronous resonances

[TABLE]

These values correspond to $\sigma/n\approx 8.3\times 10^{-3}$ and $\sigma/\nu\approx 0.77$ (synchronous resonance), and $\sigma/n\approx 4.9\times 10^{-3}$ and $\sigma/\nu\approx 0.45$ (sub/super-synchronous resonances). As $\sigma/\nu$ 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 ($\theta$). 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_{22}=1.4\times 10^{-7}$ (permanent deformation), $6.6\times 10^{-7}$ (maximum deformation in the sub/super-synchronous resonances), and $1.9\times 10^{-6}$ (maximum deformation in the synchronous resonance).

For $C_{22}=1.4\times 10^{-7}$ (see Fig. 1 top), we observe that the synchronous resonance is stable, as are the main super/sub-synchronous resonances ($\eta=n\pm\nu/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_{22}=6.6\times 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_{22}=1.9\times 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 $\eta=n\pm\nu$). 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 ($\eta=n\pm\nu/2$) are all possible, and should have comparable probabilities. However, for the highest $C_{22}$ 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)

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 $\alpha\approx 16^{\circ}$ (amplitude of 12 hr, 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 $\nu/n\approx 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 4 the orbital parameters obtained by Nesvorný et al. (2013). We note that the mutual inclination is very small (about $3^{\circ}$), and we neglect it in the following.

The bulk density of Kepler-88 $b$ is about $0.87\ \mathrm{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}\approx 0$). However, the tidally induced deformation of the planet, if it is captured in a resonance, is not negligible (see Appendices C and D)

[TABLE]

This corresponds to $\sigma/n\approx 0.018$ and $\sigma/\nu\approx 1.0$ (maximum deformation in the synchronous resonance), and $\sigma/n\approx 0.010$ and $\sigma/\nu\approx 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_{22}=2.2\times 10^{-6}$ (see Fig. 2 top), we observe that a large chaotic area is surrounding the synchronous and super/sub-synchronous resonances ($\eta=n\pm\nu/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_{22}=6.7\times 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_{22}$ 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 $\tau$, which corresponds to the relaxation time of the planet555$\tau=\tau_{v}+\tau_{e}$, where $\tau_{v}$ and $\tau_{e}$ are the viscous (or fluid) and Maxwell (or elastic) relaxation times, respectively. For simplicity, in this paper we consider $\tau_{e}=0$, since this term does not contribute to the tidal dissipation (for more details, see Correia et al., 2014)..

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)

[TABLE]

[TABLE]

where $\@vec{r}_{i}$ is the position vector of the planet $i$ (astrocentric coordinates), $\@vec{f}$ is the acceleration arising from the potential created by the deformation of the inner planet (Correia et al., 2014)

[TABLE]

and $\@vec{k}$ is the unit vector normal to the orbital plane of the inner planet. The torque acting to modify the inner planet rotation is

[TABLE]

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)

[TABLE]

where $k_{f}$ is the fluid second Love number for potential. The relaxation times $\tau$ are totally unknown for exoplanets, but if the tidal quality dissipation $Q-$factor can be estimated, then an equivalent $\tau$ can be obtained (see Correia et al., 2014). To cover all possible scenarios, in our numerical simulations we adopt a wide spectrum of $\tau$ values: $\log_{10}\tau_{\mathrm{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 $\tau=1.6$ d, and for Mars $\tau=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 ($\tau=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^{4}$ years ago is still going on, suggesting that the Earth’s mantle relaxation time is something like $\tau=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 3. The initial rotation period is set at 15.5 d, which corresponds to $\dot{\theta}/n\approx 1.14$. Since the libration frequency of KOI-227 $b$ is $\nu/n\approx 0.011$, the initial rotation rate is completely outside the resonant area, which is bounded by $\dot{\theta}/n\leq 1+\nu/n\approx 1.01$ (Fig. 1). For any $\tau$ 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 $\tau$ value we ran 1,000 simulations with slightly different initial rotation rates. The step $\Delta\dot{\theta}$ between each initial condition was chosen such that the moment at which each simulation crosses the resonance spreads equally over one eccentricity cycle. In Figure 3 we show some examples of evolution for different $\tau$ values. The evolution timescale changes with $\tau$ because the $Q-$factor was also modified.

For $\tau\leq 10^{-2}$ yr, the spin is in the low frequency regime ($\tau n\ll 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)

[TABLE]

The average eccentricity of KOI-227 $b$’s orbit is 0.057, which gives a value for the pseudo equilibrium of $\langle\dot{\theta}/n\rangle\approx 1.019$666The average rotation rate for an oscillating eccentricity is actually given by a value slightly higher than that obtained with expression (36) using the average value of the eccentricity (see Correia, 2011).. This value is already very close to the synchronous resonance, but since the libration width is $\sigma/n\approx 0.002$ (Eq. (26)), we are still outside the resonant area. Therefore, in this regime initial prograde rotations never cross any spin-orbit resonances.

For $10^{-1}\leq\tau\leq 10^{1}$ yr, the spin is in a transition of frequency regime ($\tau n\sim 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 $\tau$ 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 $\tau\geq 10^{2}$ yr, the spin is in the high frequency regime ($\tau n\gg 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 $\tau$ value. We observe that as $\tau$ increases, the number of captures in higher order resonances also increases. Indeed, for high $\tau$ values, the $C_{22}$ 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_{22}$ 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 $\tau=10^{4}$ yr, the sub-synchronous resonance can be reached after some wandering in this chaotic zone. In Figure 4 we show four examples of capture in each resonance for this $\tau$ value.

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

For gaseous planets we have $Q\sim 10^{4}$ (Lainey et al., 2009, 2012), which gives $\tau$ 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 4. The initial rotation period is set at 10 d, which corresponds to $\dot{\theta}/n\approx 1.1$. Since the libration frequency of Kepler-88 $b$ is $\nu/n\approx 0.017$, the initial rotation rate is completely outside the resonant area, which is bounded by $\dot{\theta}/n\leq 1+\nu/n\approx 1.02$ (Fig. 2). For any $\tau$ 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 Figure 5 we show some examples of evolution for different $\tau$ values.

As in the case of KOI-227 $b$, for $\tau\leq 10^{-2}$ yr, the spin is in the low frequency regime ($\tau n\ll 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 $\langle\dot{\theta}/n\rangle\approx 1.025$. The libration width for a maximum value of $C_{22}$ is $\sigma/n\approx 0.008$ (Eq. (27)), so in this regime initial prograde rotations never cross any spin-orbit resonances (Fig. 5, top). However, for $10^{-1}\leq\tau\leq 10^{0}$ yr, the equilibrium value is already inside the libration width of the synchronous resonance. Thus, for these $\tau$ values the spin can be captured in the synchronous resonance (Fig. 5, middle).

For $\tau\geq 10^{1}$ yr the spin is already in the high frequency regime, which means that capture in asynchronous resonances could be possible. However, the equilibrium $C_{22}\approx 6.7\times 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_{22}=0$) can also be chaotic when its orbit is perturbed by a companion planet.

In Figure 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_{22}$ 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_{22}>6.7\times 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 synchronised.

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 transit timing variations (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\ \mathrm{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 ($\tau$, 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).

Acknowledgements.

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.

Appendix A Near-resonance case

In this appendix we show how the quasi-periodic decomposition of $\Lambda_{i}$, $\lambda_{i}$ (see Eq. (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

[TABLE]

The secular evolution of the system is obtained by averaging this Hamiltonian over the fast angles $\lambda_{i}$ ($i=1,2$). This averaging transformation is obtained by a change of coordinates that is close to identity. We denote by $\Lambda_{i}^{\prime}$, $\lambda_{i}^{\prime}$, $D_{i}^{\prime}$, and $\varpi_{i}^{\prime}$ ($i=1,2$) the new coordinates, by $\mathcal{H}^{\prime}$ the new Hamiltonian, and by $W$ the generating Hamiltonian of the transformation. To first order in the planet masses, the change of coordinates reads

[TABLE]

and the Hamiltonian $\mathcal{H}^{\prime}$ is given by

[TABLE]

with

[TABLE]

and where the braces denote the Poisson brackets. We thus have

[TABLE]

which is called the homological equation and whose solution is

[TABLE]

where $n_{i}$ ($i=1,2$) are the unperturbed Keplerian mean-motions of the planets.

The secular Hamiltonian $\mathcal{H}^{\prime}$ no longer depends on the mean longitudes $\lambda_{i}^{\prime}$, which implies that the coordinates $\Lambda_{i}^{\prime}$ are constants of motion. The long-term variations of $D_{i}^{\prime}$ (i.e. eccentricities), $\varpi_{i}^{\prime}$, and $\lambda_{i}^{\prime}$, can be solved by using the Hamiltonian $\mathcal{H}^{\prime}$. 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 $D_{i}^{\prime}$ and $\varpi_{i}^{\prime}$ are constants and that $\lambda_{i}^{\prime}=\lambda_{i,0}^{\prime}+n_{i}t$. In order to obtain the real evolution of the system, we need to revert to the original coordinate system using Eqs. (38)- (41). In particular, we have

[TABLE]

Let us consider a system that is close to a $k_{2}$:$k_{1}$ resonance but outside of it. The combination $\nu=k_{2}n_{2}-k_{1}n_{1}$ 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

[TABLE]

with $k=(-k_{1},k_{2})$, $C_{k}=(-1)^{i+1}|C_{k}|\mathbf{e}^{\mathrm{i}\phi_{i}}$, and $C_{-k}$ is its complex conjugate (by construction). The degree of the resonance $k_{2}$:$k_{1}$ is denoted $q=k_{2}-k_{1}$. The coefficient $C_{k}$ is of order $q$ in eccentricity (D’Alembert rule), and can be written (e.g. Laskar & Robutel 1995)

[TABLE]

where $e=\max(e_{1},e_{2})$ and $c_{k}$ is of the order of unity. Finally, we have (see Eq. (16))

[TABLE]

where $m_{p}$ is the mass of the perturbing planet. We observe that when the system is very close to the resonance separatrix, $\nu\ll n$, the amplitude of oscillations increases (see Eq. (51)).

Appendix B Resonant case

In this appendix we show how the quasi-periodic decomposition of $\Lambda_{i}$, $\lambda_{i}$ (see Eq. (13)) can be obtained in the resonant case. The resonant case arises when the small divisor of Eq. (47) 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 $\mathcal{H}^{\prime}$ is constructed such that resonant terms (of the form $p(k_{2}\lambda_{2}-k_{1}\lambda_{1})$) are kept,

[TABLE]

The solution to the homological equation (Eq. (45)) is in this case

[TABLE]

and the resonant Hamiltonian reads

[TABLE]

where the coefficients $C_{l}$ are functions of $\Lambda_{i}^{\prime}$, $D_{i}^{\prime}$, $\varpi_{i}^{\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 $\Lambda^{\prime}$ (e.g. Delisle et al. 2014), and in first approximation we assume they are constant. According to Eq. (49), we have

[TABLE]

where $q=k_{2}-k_{1}$ is the degree of the resonance, and $c_{(-pk_{1},pk_{2})}$ is of the order of unity. In a first approximation, we only keep the terms of order $q$ in eccentricity, and obtain the Hamiltonian

[TABLE]

where $C_{0}$ has been dropped since we assumed it to be constant. We introduce the canonical change of coordinates

[TABLE]

Since $\mathcal{H}^{\prime}$ only depends on the angle $\psi$ and not on $\xi$, $\Gamma$ is a conserved quantity. We now develop the Keplerian part $\mathcal{H}_{0}^{\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$,

[TABLE]

with

[TABLE]

With this approximation, the Hamiltonian (56) reads

[TABLE]

which is the Hamiltonian of a simple pendulum. The frequency of libration is given by

[TABLE]

where $K$ is the complete elliptical integral of the first kind, and the amplitude of libration $\psi_{max}$ is in the range $[0,\pi]$. For small amplitude oscillations, we have

[TABLE]

and from Eqs. (57), (58), we obtain

[TABLE]

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 $\mathcal{H}^{\prime}$. We thus have $\Lambda_{i}\approx\Lambda_{i}^{\prime}$. Finally, from Eqs. (16), (64), (65), we deduce

[TABLE]

In particular, for $m_{1}\ll m_{2}$, we have $\alpha=\psi_{max}/k_{1}$ for the inner planet (1), and $\alpha=0$ for the outer planet (2). On the contrary, for $m_{1}\gg m_{2}$, we have $\alpha=0$ for the inner planet, and $\alpha=\psi_{max}/k_{2}$ 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 ($\sigma$) from the known properties of the planets. From Eq. (6), we have

[TABLE]

The internal structure factor $\zeta$ can be estimated from $k_{f}$ through the Darwin-Radau equation (e.g. Jeffreys 1976)

[TABLE]

We have $\zeta=2/5$ ($k_{f}=3/2$) for a homogenous sphere, $\zeta\approx 1/3$ for rocky planets, and $\zeta\approx 1/4$ for gaseous planets.

The deformation of the planet ($C_{22}$) 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)

[TABLE]

For gaseous planets, the permanent asymmetry is very weak ($C_{22,r}\approx 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 ($\tau$, 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. (4.2) with $\tau=0$, see also Correia & Rodríguez 2013)

[TABLE]

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. (2):

[TABLE]

with

[TABLE]

We assume in the following that the relaxation time $\tau$ is much longer that the variations of $C_{22,t,\text{inst.}}$, such that the actual $C_{22,t}$ of the planet is the mean value of $C_{22,t,\text{inst.}}$ (see Eq. (4.2)):

[TABLE]

From Eq. (C), 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_{22}$ 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

[TABLE]

when the amplitude of libration is small ($\sin 2\gamma\approx 0$, $\cos 2\gamma\approx 1$). Similarly, if the spin is locked in a sub/super-synchronous resonance, we have

[TABLE]

at the centre of the resonance.

Let us apply this reasoning to KOI-227 $b$ (rocky planet). We obtain $C_{22,r}=1.4\times 10^{-7}$ for the permanent deformation (see Eq. (69)). The tidal deformation is $C_{22,t}=2.6\times 10^{-6}$ at the centre of the synchronous resonance (see Eq. (74)), and $C_{22,t}=7.9\times 10^{-7}$ at the centre of the sub/super-synchronous resonance (see Eq. (75)). The maximum deformation ($C_{22}=C_{22,r}+C_{22,t}$) is thus $2.8\times 10^{-6}$ in the synchronous case, and $9.3\times 10^{-7}$ in the sub/super-synchronous case.

In the case of Kepler-88 $b$ (gaseous planet) we obtain $C_{22,r}=0$, $C_{22}=C_{22,t}=2.0\times 10^{-5}$ for the synchronous resonance, and $C_{22}=C_{22,t}=5.6\times 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 $\alpha$, and frequency $\nu\ll n$. We neglect here the semi-major axis variations since their contribution to the spin evolution is of order $\nu/n$ compared to the contribution of the mean longitude variations (see Eqs. (15), (19), and (2)).

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

[TABLE]

where (see Eq. (67))

[TABLE]

and (see Eq. (70))

[TABLE]

We introduce $h=\gamma+\alpha\sin(\nu t+\phi)$, such that

[TABLE]

with

[TABLE]

Equations (79), (80) can be developed in power series of $h$

[TABLE]

At first order (linearized equation), the forced solution is simply

[TABLE]

which is of order $\alpha$. At higher orders, odd harmonics of the forced frequency $\nu$ appears, and $h$ takes the form

[TABLE]

where $h_{k}$ is of order $\alpha^{k}$. We replace Eq. (86) in Eqs. (83) and (84), and truncate the resulting expression at a given order $N$ in $\alpha$. By identifying terms of frequency $k\nu$ ($k\in[1,N]$), we obtain $N$ polynomial equations on the coefficients $h_{1}$,…,$h_{N}$. Solving this set of equations allows us to determine the coefficients $h_{i}$, as well as the corresponding $C_{22}$ value. The forced oscillations of $\gamma$ are then given by

[TABLE]

with $\gamma_{1}=h_{1}-\alpha$, $\gamma_{k}=h_{k}$ ($k>1$).

Applying this reasoning to KOI-227 $b$, we obtain $C_{22,t}=1.8\times 10^{-6}$ (instead of $2.6\times 10^{-6}$) and $C_{22}=C_{22,t}+C_{22,r}=1.9\times 10^{-6}$ (instead of $2.8\times 10^{-6}$). The forced amplitude at frequency $\nu$ is $\gamma_{1}=17^{\circ}$, while the amplitudes of harmonics ($3\nu$, etc.) decrease rapidly ($\gamma_{k+2}/\gamma_{k}\sim 10^{-2}$). In the case of Kepler-88 $b$, we obtain $C_{22}=C_{22,t}=6.7\times 10^{-6}$ (instead of $2.0\times 10^{-5}$) and a forced amplitude of $36^{\circ}$ (at frequency $\nu$). As was true for KOI-227 $b$, for Kepler-88 $b$ the amplitudes of harmonics decrease rapidly.

D.2 Sub/super-synchronous case

We now consider the same effect (forced oscillations at frequency $\nu$ and its harmonics) but for the sub/super-synchronous resonances. We assume that the spin is locked in one of these resonances, such that

[TABLE]

with $\phi_{+}=\phi+\pi$, $\phi_{-}=\phi$. As for the synchronous resonance, we introduce $h=\gamma_{\text{forced}}+\alpha\sin(\nu t+\phi)$ such that

[TABLE]

As we did for the synchronous case, these expressions can be developed in power series of h. Replacing $h$ by

[TABLE]

and truncating at a given order $N$ in $\alpha$ ($h_{k}$ being of order $\alpha^{k}$), we obtain a set of $N$ polynomial equations on $h_{1}$, …, $h_{N}$. We then solve for $h_{i}$ and determine the corresponding $C_{22}$ value.

In the case of KOI-227 $b$, we obtain $C_{22,t}=5.2\times 10^{-7}$ (instead of $7.9\times 10^{-7}$) and $C_{22}=6.6\times 10^{-7}$ (instead of $9.3\times 10^{-7}$), with a forced amplitude of $\gamma_{1}=-6^{\circ}$ (at frequency $\nu$). The negative sign means that the forced oscillations are dephased by an angle $\pi$ with respect to the TTV signal. In the case of Kepler-88 $b$, we find $C_{22}=C_{22,t}=2.2\times 10^{-6}$ (instead of $5.6\times 10^{-6}$), with a forced amplitude of $-10^{\circ}$. As in the case of the synchronous resonance, the amplitudes of harmonics ($3\nu$, etc.) decrease rapidly for both planets ($\gamma_{k+2}/\gamma_{k}\sim 10^{-3}$).