# Spin-orbit coupling for close-in planets

Alexandre Correia, Jean-Baptiste Delisle

► To cite this version:

Alexandre Correia, Jean-Baptiste Delisle. Spin-orbit coupling for close-in planets. Astronomy and Astrophysics - A&A, 2019, 630, pp.A102. 10.1051/0004-6361/201936336. hal-03144466

HAL Id: hal-03144466

https://hal.science/hal-03144466

Submitted on 13 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spin-orbit coupling for close-in planets

Alexandre C. M. Correia<fnref n="1" />,<fnref n="2" /> and Jean-Baptiste Delisle<fnref n="2" />,<fnref n="3" />

Received 17 July 2019 / Accepted 20 August 2019

## ABSTRACT

We study the spin evolution of close-in planets in multi-body systems and present a very general formulation of the spin-orbit problem. This includes a simple way to probe the spin dynamics from the orbital perturbations, a new method for computing forced libration and tidal deformation, and general expressions for the tidal torque and capture probabilities in resonance. We show that planet-planet perturbations can drive the spin of Earth-size planets into asynchronous or chaotic states, even for nearly circular orbits. We apply our results to Mercury and to the KOI-1599 system of two super-Earths in a 3/2 mean motion resonance.

**Key words.** celestial mechanics – planets and satellites: general – chaos

## 1. Introduction

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

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

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

resonant rotation with the orbital libration frequency (Correia & Robutel 2013; Leleu et al. 2016; Delisle et al. 2017).

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

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

## 2. Orbital forcing

We consider a system composed by a central star of mass *M*, a companion planet of mass *m*, and additional companions in coplanar orbits. The planet is considered an oblate ellipsoid with mean radius *R* and gravity field coefficients given by $C_{22}$ and

<footnote n="1">¹ CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra, Portugal
e-mail: acor@uc.pt</footnote>

<footnote n="2">² ASD, IMCCE, Observatoire de Paris, PSL Université, 77 Av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote n="3">³ Observatoire de l'Université de Genève, 51 chemin des Maillettes, 1290 Sauverny, Switzerland</footnote>

<footnote n="1">¹ The expected libration width in resonance does not seem sufficient to cover this scenario. The current explanation is that Venus rotation is an equilibrium between gravitational and thermal atmospheric tides (Gold & Soter 1969; Dobrovolskis 1980; Correia & Laskar 2001).</footnote>
$S_{22}$. We furthermore assume that the spin axis of the planet is orthogonal to the orbital plane (zero obliquity). The torque of the star on the planet changes the rotation angle, $\theta$, whose variations are given by (e.g., Correia et al. 2014):

$$ \ddot{\theta} = -\frac{6GMmR^2}{Cr^3} [C_{22} \sin 2(\theta - f) + S_{22} \cos 2(\theta - f)], \quad (1) $$

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

The above expression can be rewritten in a more compact and usable form with complex numbers as

$$ \ddot{\theta} = \mathcal{B} \mathfrak{F} (Z x e^{-2i\theta}), \quad (2) $$

with

$$ Z = C_{22} - iS_{22}, \quad (3) $$

$$ x = \left(\frac{a}{r}\right)^3 e^{2if}, \quad (4) $$

$$ \mathcal{B} = \frac{6GMmR^2}{Ca^3} \approx 18n^2, \quad (5) $$

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

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

$$ X = r e^{if} = (r \cos f, r \sin f) \quad (6) $$

be the position vector of the planet, whose spin we want to study. The evolution of X can be obtained by integrating a system of N point-mass bodies (the star, the planet and the remaining companions). The orbital forcing on the spin (Eq. (4)) can then be directly computed from X as

$$ x = \frac{a^3 X^2}{|X|^5}. \quad (7) $$

Assuming quasi-periodic orbital motion, we can expand the previous expression as a Fourier series

$$ x = \sum_k x_k e^{ik \cdot v_t}, \quad (8) $$

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

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

### 3. Rigid bodies

As in most previous studies on spin-orbit coupling, let us assume for now that $Z = Z_0$, either because $|Z_t| \ll |Z_0|$ (case of Mercury) or because $Z_t$ is frozen from a time the planet was more fluid (case of the Moon). Then, Eq. (2) becomes<fnref n="2" />

$$ \ddot{\theta} = \mathcal{B} \mathfrak{F} \left( Z_0 \sum_k x_k e^{i(k \cdot v_t - 2\theta)} \right). \quad (9) $$

Each term in this expression individually behaves like a pendulum, with islands of rotational libration centered at

$$ \langle \dot{\theta} \rangle = \nu_k = \frac{k \cdot v}{2}, \quad (10) $$

and libration width

$$ \beta_k = \sqrt{2\mathcal{B}|Z_0x_k|} \approx 6n \sqrt{|Z_0|} \sqrt{|x_k|}. \quad (11) $$

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

$$ \bar{\beta} = \beta_{k_1} + \beta_{k_2} \approx \nu_{k_2} - \nu_{k_1} = \Delta v. \quad (12) $$

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

#### 3.1. Single planet case

In the classical case of a single planet in an elliptical orbit, the orbital forcing (Eq. (4)) is given by

$$ x = \sum_k X_k^{-3,2} (e) e^{2i\varpi} e^{iknt}, \quad (13) $$

<footnote n="2">² Many authors adopt a frame where $Z_0$ is a real number ($Z_0 = C_{22,0}$), which simplifies expression (1). This is correct for the rigid case, but it is no longer valid when we consider the tidal deformation of the planet.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Fast Fourier transform of the orbital forcing $x$ (Eq. (7)) for the two-body problem with eccentricity $e = 0.165$. The horizontal lines give the amplitudes $\beta_k$ for $|Z_0| = 1.5 \times 10^{-5}$, which corresponds to Mercury's value (Genova et al. 2013), and $|Z_0| = 2.9 \times 10^{-2}$, which corresponds to Hyperion's value (Thomas et al. 2007).</figcaption></figure>

where $X_k^{l,m}(e)$ are the Hansen coefficients (Hughes 1981) and $\omega$ is the longitude of the pericenter. Thus, resonant islands occur for $v_k = kn/2$, with amplitude $|x_k| = X_k^{-3,2}(e)$. The two largest terms correspond to $k = 2$ and $k = 3$, for which $X_2^{-3,2}(e) \approx 1$ and $X_3^{-3,2}(e) \approx 7e/2$, respectively (we neglect terms in $e^2$). Then

$$ \bar{\beta}/\Delta\nu \approx 12 \sqrt{|Z_0|} \left(1 + \sqrt{7e/2}\right). \qquad (14) $$

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

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

## 3.2. Multi-planet systems

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

<table><caption>**Table 1.** Parameters for the KOI-1599 system used in this study.</caption><thead><tr><th>Parameter</th><th>(Unit)</th><th>Inner</th><th>Outer</th></tr></thead><tbody><tr><td><i>m</i></td><td>(M<sub>⊕</sub>)</td><td>8.3</td><td>4.3</td></tr><tr><td><i>R</i></td><td>(R<sub>⊕</sub>)</td><td>1.9</td><td>1.9</td></tr><tr><td><i>a</i></td><td>(au)</td><td>0.11230</td><td>0.14727</td></tr><tr><td><i>e</i></td><td></td><td>0.0067</td><td>0.0177</td></tr><tr><td><i>λ</i></td><td>(deg)</td><td>217.4</td><td>251.5</td></tr><tr><td><i>ω</i></td><td>(deg)</td><td>40.6</td><td>220.1</td></tr></tbody></table>

**Notes.** The mass of the star is $1.02 M_{\odot}$. This orbital solution corresponds to the stable fit GEA I solution obtained by Panichi et al. (2019).

mutual perturbations are strong and need to be taken into account in the spin-orbit coupling analysis.

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

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

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

We need to estimate the rotational libration width $\beta_k$ for each frequency in Table 2 to determine the behavior of the spin. Since we know the orbital forcing amplitudes, the only
<table><caption>**Table 2.** Quasi-periodic decomposition of the orbital forcing function $x$ for the outer planet in the KOI-1599 system (Eq. (8)).</caption><thead><tr><th colspan="3">k · v</th><th rowspan="2">v<sub>k</sub>/n - 1</th><th rowspan="2">|x<sub>k</sub>|</th><th rowspan="2">φ<sub>k</sub> (°)</th></tr><tr><th>n</th><th>ν</th><th>g</th></tr></thead><tbody><tr><td>2</td><td>0</td><td>0</td><td>0.000000</td><td>0.93142</td><td>173.21</td></tr><tr><td>2</td><td>1</td><td>0</td><td>0.003316</td><td>0.25878</td><td>79.54</td></tr><tr><td>2</td><td>-1</td><td>0</td><td>-0.003316</td><td>0.24821</td><td>86.87</td></tr><tr><td>2</td><td>-2</td><td>0</td><td>-0.006632</td><td>0.04461</td><td>0.54</td></tr><tr><td>3</td><td>-1</td><td>-1</td><td>0.498128</td><td>0.04138</td><td>150.57</td></tr><tr><td>3</td><td>0</td><td>-1</td><td>0.501445</td><td>0.03730</td><td>-123.10</td></tr><tr><td>2</td><td>2</td><td>0</td><td>0.006632</td><td>0.02197</td><td>-14.13</td></tr><tr><td>3</td><td>-2</td><td>-1</td><td>0.494812</td><td>0.01212</td><td>64.24</td></tr><tr><td>3</td><td>1</td><td>-1</td><td>0.504761</td><td>0.01070</td><td>143.24</td></tr><tr><td>1</td><td>0</td><td>1</td><td>-0.501445</td><td>0.00753</td><td>-70.49</td></tr><tr><td>2</td><td>-3</td><td>0</td><td>-0.009949</td><td>0.00676</td><td>-85.79</td></tr></tbody></table>

**Notes.** We provide the dominating terms of $x$ by decreasing amplitude ($|x_k| > 5 \times 10^{-3}$). We have $v_k = \frac{1}{2}k \cdot v$ (Eq. (10)), and $\phi_k = \arg x_k$. $v = (n, v, g)$ is the fundamental frequencies vector of the orbital motion, where $n = 112.305$ rad yr$^{-1}$ is the mean motion, $v = 0.744857$ rad yr$^{-1}$ the libration frequency, and $g = -0.324472$ rad yr$^{-1}$ the precession rate.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Fast Fourier transform of the orbital forcing $x$ (Eq. (7)) for the outer planet in the KOI-1599 system. The horizontal lines give the amplitudes $\beta_k$ for $|Z_0| = 1.5 \times 10^{-6}$ and $|Z_0| = 1.0 \times 10^{-7}$.</figcaption></figure>

uncertainty is related to the permanent deformation of the planet, $Z_0$ (Eq. (11)). We can adopt the Earth's value, $|Z_0| = 1.5 \times 10^{-6}$ because KOI-1599 planets are likely rocky, as their densities are similar to that of the Earth. However, their surface gravity is a bit stronger, so maybe the inhomogeneities in the mass distribution are smaller, such that $|Z_0| \sim 10^{-7}$. These two values can be seen as maximal and minimal estimations for $|Z_0|$, respectively. In Fig. 2 we horizontally draw the libration width of the islands centered at each $v_k$ for both $|Z_0|$ values over the Fourier spectra of the orbital forcing $x$. We observe that, for $|Z_0| \sim 10^{-7}$, each island is isolated, so the rotation can be trapped in any of them, which includes the synchronous resonance, along with the super- and sub-synchronous resonances. On the other hand, for $|Z_0| \sim 10^{-6}$, the libration widths overlap, so we can expect chaotic rotation.

In order to study in more detail the spin dynamics, we have to run direct numerical simulations of the rotation angle of the outer planet. We then integrated Eq. (1) together with the orbital solution of the KOI-1599 system given by a classical N-body

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Frequency map analysis of KOI-1599 outer planet rotation angle assuming a rigid body with $|Z_0| = 1.5 \times 10^{-6}$ (top) or $|Z_0| = 1.0 \times 10^{-7}$ (bottom). Each dot corresponds to a grid with 181 × 201 different initial conditions of $\gamma(0)$ and $\dot{\gamma}(0)$, respectively, with $\gamma(t) = \theta(t) - \lambda(0)$. The color index gives the derivative $\partial\eta/\partial\dot{\gamma}(0)$, where $\eta$ is the main frequency of $\gamma(t)$. The blue dots correspond to resonant motion, green dots to nonresonant regular motion, and red dots to chaotic motion.</figcaption></figure>

integrator. In Fig. 3 we show a frequency map analysis of the rotation angle (Laskar 1993; Leleu et al. 2016) for the two $|Z_0|$ values in the neighborhood of the synchronous resonance. We conclude that our predictions for the spin of the outer planet based on the Chirikov (1979) general criterion are correct: for $|Z_0| \sim 10^{-7}$ the planet can present stable asynchronous rotation, while for $|Z_0| \sim 10^{-6}$ it likely presents chaotic rotation.

## 4. Tidal deformation and forced oscillations

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

One of the simplest viscoelastic models is to consider that the planet behaves like a Maxwell material, which is represented by a purely viscous damper and a purely elastic spring connected in series (e.g., Turcotte & Schubert 2002). More complex models, such as the Andrade model (Andrade 1910), can be thought of as the Maxwell model equipped with extra terms describing hereditary reaction of strain to stress (Efroimsky 2012). Therefore, the main features of tidal response are similar to the Maxwell model.
For simplicity, we adopt the Maxwell viscoelastic rheology. The viscous response of the planet to the tidal excitation is mod- eled by the parameter $\tau$, which corresponds to the relaxation time of the planet, while the elastic response in modeled by the parameter $\tau_e$. The deformation Z follows (Correia et al. 2014)

$$Z + \tau \dot{Z} = Z_0 + Z_e + \tau_e \dot{Z}_e, \quad (15)$$

with

$$Z_e = \mathcal{A} \bar{x} e^{2i\theta}, \quad \text{where} \quad \mathcal{A} = \frac{k_f}{4} \frac{M}{m} \left(\frac{R}{a}\right)^3, \quad (16)$$

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

We denote the average rotation rate of the planet by $\omega = \langle \dot{\theta} \rangle$. For simplicity, we assume that $\theta$ has a small forced amplitude of libration due to orbital forcing and a zero free amplitude of libration, such that

$$\theta = \omega t + \sum_k \theta_k e^{ik \cdot vt}, \quad (17)$$

$$\mathbf{e}^{2i\theta} = \mathbf{e}^{2i\omega t} \sum_{j \ge 0} \frac{(i \sum_k 2\theta_k \mathbf{e}^{ik \cdot vt})^j}{j!} = \sum_l y_l \mathbf{e}^{i(2\omega+l\cdot v)t}, \quad (18)$$

where all $\theta_k$ are small, and the series is absolutely convergent ($\sum_k |\theta_k|$ converges). Each coefficient $y_l$ can be expressed as a power series of the coefficients $\theta_k$, by identifying terms on both sides of Eq. (18).

From expressions (8), (16), and (18), we obtain

$$Z_e = \mathcal{A} \sum_k \left( \sum_l \bar{x}_{k+l} y_l \right) e^{i(2\omega - k \cdot v)t}, \quad (19)$$

and replacing this in expression (15), we get

$$Z = Z_0 + Z_t + C e^{-t/\tau}, \quad (20)$$

with

$$Z_t = \mathcal{A} \sum_k \left( \sum_l \bar{x}_{k+l} y_l \right) \frac{1 + i(2\omega - k \cdot v)\tau_e}{1 + i(2\omega - k \cdot v)\tau} e^{i(2\omega - k \cdot v)t}. \quad (21)$$

In the following, we neglect the term $Ce^{-t/\tau}$ because it vanishes after some time, and we only consider the solution $Z = Z_0 + Z_t$.

### 4.1. No resonant rotation

We assume that $\omega \neq \nu_k = k \cdot v/2$ for all $k$ (no resonant rotation). Then, over one rotational period we have (Eq. (21))

$$\langle Z_t \rangle = 0, \quad (22)$$

that is, there is no permanent tidal deformation.

### 4.2. Resonant rotation without forced oscillations

The situation changes if the rotation is locked in some spin-orbit resonance, that is, we assume that $\omega = \nu_r = r \cdot v/2$ (Eq. (10)). We can then rewrite expression (20) as

$$Z = \sum_k z_k e^{i(r-k)\cdot vt}, \quad (23)$$

with

$$z_k = Z_0 \delta_{r,k} + \mathcal{A} b_{r-k} \sum_l \bar{x}_{k+l} y_l, \quad (24)$$

and

$$b_k = \frac{1 + i(\mathbf{k} \cdot \mathbf{v})\tau_e}{1 + i(\mathbf{k} \cdot \mathbf{v})\tau}. \quad (25)$$

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

In the case of very weak or no forced oscillations, expres- sions (17) and (18) simplify as

$$\theta = \theta_0 + \frac{\mathbf{r} \cdot \mathbf{v}}{2} t, \quad (26)$$

$$\mathbf{e}^{2i\theta} = \mathbf{e}^{2i\theta_0} \mathbf{e}^{ir\cdot vt}, \quad (27)$$

that is, $y_0 = \mathbf{e}^{2i\theta_0}$, and $y_{k \neq 0} = 0$. Replacing in expression (24), we get for the deformation of the planet

$$z_k = Z_0 \delta_{r,k} + \mathcal{A} b_{r-k} \mathbf{e}^{2i\theta_0} \bar{x}_k. \quad (28)$$

Finally, to get the permanent tidal deformation, we can average over one rotational period

$$\langle Z_t \rangle = \langle Z \rangle - Z_0 = \mathcal{A} \mathbf{e}^{2i\theta_0} \bar{x}_r, \quad (29)$$

or

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

For instance, in the classical case of a single planet in an elliptical orbit (Sect. 3.1), we have (Eq. (13))

$$\bar{x}_k = X_k^{-3,2}(e) \mathbf{e}^{-2i\varpi}, \quad (31)$$

thus, for a planet trapped in the $\nu_k/n = k/2$ resonance

$$|\langle Z_t \rangle| = \mathcal{A} X_k^{-3,2}(e). \quad (32)$$

For the outer planet in the KOI-1599 system (Sect. 3.2), we have for the synchronous resonance (Table 2)

$$|x_{[2,0]}| \approx 0.93 \Rightarrow |\langle Z_t \rangle| \approx 2.5 \times 10^{-6}, \quad (33)$$

while for the super- and sub-synchronous resonances

$$|x_{[2, \pm 1]}| \approx 0.25 \Rightarrow |\langle Z_t \rangle| \approx 7.0 \times 10^{-7}. \quad (34)$$

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

### 4.3. Resonant rotation with forced oscillations

In case of resonant motion ($\omega = \nu_r$), replacing the deformation of the planet, $Z$ (Eq. (23)), in expression (2), we get

$$\ddot{\theta} = \mathcal{B} \mathcal{J} \left( \sum_k \left( \sum_{l,m} x_{k+l+m} \bar{y}_l z_m \right) \mathbf{e}^{ik \cdot vt} \right), \quad (35)$$

and from expression (17)

$$\ddot{\theta} = - \sum_k (\mathbf{k} \cdot \mathbf{v})^2 \theta_k \mathbf{e}^{ik \cdot vt}. \quad (36)$$
Then, by identifying terms in both series, we find

$$ \theta_k = \frac{1}{2i} (\alpha_k - \bar{\alpha}_{-k}), \qquad (37) $$

where

$$ \alpha_k = \frac{\mathcal{B}}{(\mathbf{k} \cdot \mathbf{v})^2} \sum_{l,m} \bar{x}_{l+m-k} y_l \bar{z}_m. \qquad (38) $$

We also note that $\theta_{-k} = \bar{\theta}_k$. Therefore, with $\phi_k = \arg \alpha_k$, we can rewrite expression (17) as

$$ \theta = v_r t + \theta_0 + \sum_{k \neq 0} |\alpha_k| \sin(\mathbf{k} \cdot \mathbf{v} t + \phi_k). \qquad (39) $$

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

We can estimate the order of magnitude of the expected forced oscillations. Since $|x_k| \le 1$, $|y_k| \le 1$ and $|z_k| \le \mathcal{Z} = \max(|Z_0|, \mathcal{A})$, we deduce from expression (37)

$$ |\theta_k| \le \frac{\mathcal{B}\mathcal{Z}}{(\mathbf{k} \cdot \mathbf{v})^2} \approx 18 \mathcal{Z} \left(\frac{n}{2v_k}\right)^2. \qquad (40) $$

For an Earth-size planet and a perturbing frequency on the order of the orbital frequency ($v_k \sim n$), we get $|\theta_k| \lesssim 10^{-5}$, that is, the amplitude of the forced oscillations remains very small. We then linearize expression (18) and find

$$ y_0 = \mathbf{e}^{i2\theta_0}, \quad \text{and} \quad y_{l \neq 0} \approx 2i \theta_l \mathbf{e}^{i2\theta_0} = (\alpha_l - \bar{\alpha}_{-l}) \mathbf{e}^{i2\theta_0}. \qquad (41) $$

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

For the sake of simplicity and readability, we restrict our study to the first order case in the following, as it already provides a good approximation in most cases. Then, from expressions (38) and (41) we have

$$ \alpha_k \approx \frac{18 n^2 \mathbf{e}^{i2\theta_0}}{(\mathbf{k} \cdot \mathbf{v})^2} \sum_m \bar{x}_{m-k} \bar{z}_m, \qquad (42) $$

where $\bar{z}_m$ is given by expression (24). In order to be captured in a given resonance ($\omega = v_r$), we have that $x_r \neq 0$ (Sect. 3). Since we always have $\mathbf{k} \neq \mathbf{0}$ (Eq. (39)), it is needed that at least another $x_{m+r} \neq 0$ to observe forced librations. Indeed, if there is only one orbital forcing term, as in the case of a single planet in a circular orbit, all $\alpha_k = 0$, and hence there are no forced librations.

### 4.3.1. Single planet case

We first consider a single planet in an elliptical orbit with small eccentricity. The two dominating terms of the orbital forcing are $x_2 \approx 1$ and $x_3 \approx 7e/2$ (Sect. 3.1), so we neglect the remaining

terms. Applying to the case of Mercury, we have $r = 3$, $v = n$, and $|Z_0| \gg \mathcal{A}$. Thus, expression (42) becomes

$$ \alpha_k \approx \frac{18 \mathbf{e}^{i2\theta_0}}{k^2} \bar{Z}_0 \bar{x}_{3-k}. \qquad (43) $$

The only $\alpha_k \neq 0$ is for $x_{3-k} = x_2$, that is, for $k = 1$. Then, the rotation angle becomes (Eq. (39))

$$ \theta \approx \frac{3}{2}nt + \theta_0 + 18|Z_0 \bar{x}_2| \sin(nt + \phi_1). \qquad (44) $$

This simple approximation gives a forced libration period equal to the orbital period, and a forced libration amplitude $|\alpha_1| \approx 18|Z_0| = 56''$, in perfect agreement with a full numerical simulation for Mercury (see Fig. 2 in Peale et al. 2007).

### 4.3.2. Multi-planet systems

We consider the outer planet in the KOI-1599 system, for which the dominating terms of the orbital forcing are $x_{[2,0]} = 0.93$, $x_{[2,1]} = 0.26$, and $x_{[2,-1]} = 0.25$ (Table 2), and we neglect the remaining terms. In this case, we have $|Z_0| \ll \mathcal{A}$, so replacing expression (24) in (42) gives

$$ \alpha_k \approx \frac{18 \mathcal{A} n^2 \mathbf{e}^{i2\theta_0}}{(\mathbf{k} \cdot \mathbf{v})^2} \sum_m b_{m-r} \bar{x}_{m-k} \sum_l x_{m+l} \bar{y}_l, \qquad (45) $$

which depends on the tidal model through $b_{m-r}$ (Eq. (25)). For simplicity, we further assume that $\tau_e = 0$ and $\tau v_k \gg 1$ (high frequency regime, see Sect. 5), such that $b_{m-r} \approx \delta_{r,m}$. Then

$$ \alpha_k \approx \frac{18 \mathcal{A} n^2 \mathbf{e}^{i2\theta_0}}{(\mathbf{k} \cdot \mathbf{v})^2} \bar{x}_{r-k} \sum_l x_{r+l} \bar{y}_l. \qquad (46) $$

Let us consider, for example, that the planet is captured in the $r = [2, 1]$ spin-orbit resonance, with frequency $v_{[2,1]} = n + v/2$. The only two $\alpha_k \neq 0$ are for $\bar{x}_{[2,1]} - k = \bar{x}_{[2,0]}$ and $\bar{x}_{[2,-1]}$, that is, for $\mathbf{k} = [0, 1]$ and $[0, 2]$, respectively, with associated frequencies $v$ and $2v$, where $v$ is the orbital libration frequency

$$ \alpha_v \approx 18 \mathcal{A} \mathbf{e}^{i2\theta_0} \left(\frac{n}{v}\right)^2 \bar{x}_{[2,0]} \sum_l x_{[2,1]+l} \bar{y}_l, \qquad (47) $$

$$ \alpha_{2v} \approx 18 \mathcal{A} \mathbf{e}^{i2\theta_0} \left(\frac{n}{2v}\right)^2 \bar{x}_{[2,-1]} \sum_l x_{[2,1]+l} \bar{y}_l. \qquad (48) $$

Since $|\alpha_{2v}/\alpha_v| \approx 0.067$, we can also neglect the contribution of $\alpha_{2v}$. Thus, the rotation angle (Eq. (39)) simply becomes

$$ \theta \approx (n + \frac{1}{2}v)t + \theta_0 + |\alpha_v| \sin(vt + \phi_v), \qquad (49) $$

and the rotation rate

$$ \dot{\theta} \approx n + v \left[ \frac{1}{2} + |\alpha_v| \cos(vt + \phi_v) \right]. \qquad (50) $$

From expression (41) we additionally have

$$ \sum_l x_{[2,1]+l} \bar{y}_l = x_{[2,1]} \bar{y}_{[0,0]} + x_{[2,0]} \bar{y}_{[0,-1]} + x_{[2,-1]} \bar{y}_{[0,-2]} \\ \approx \mathbf{e}^{-2i\theta_0} (x_{[2,1]} - x_{[2,0]} \alpha_v). \qquad (51) $$

Inserting in expression (47) and solving for $\alpha_v$, we finally get

$$ \alpha_v \approx \frac{18 \mathcal{A} (n/v)^2 \bar{x}_{[2,0]} \bar{x}_{[2,1]}}{1 + 18 \mathcal{A} (n/v)^2 |x_{[2,0]}|^2}. \qquad (52) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Forced libration of the rotation rate for the super-synchronous resonance ($v_r = n + v/2$) of the outer planet in the KOI-1599 system. We show the evolution of $(\theta - v_r)/n$ adopting $\tau_e = 0$ and $\tau n = 10^4$.</figcaption></figure>

With $\mathcal{A} = 2.7 \times 10^{-6}$ (Eq. (16)), and $(n/v) \approx 150$ (Table 2), we estimate $|\alpha_v| \approx 0.14$ rad. From expression (50) we also estimate $\Delta\theta/n = (v/n)|\alpha_v| \approx 9 \times 10^{-4}$. In Fig. 4 we show the forced libra- tions of the rotation rate for the $\mathbf{r} = [2, 1]$ spin-orbit resonance, which is obtained performing a numerical simulation of the KOI- 1599 system with $\tau n = 10^4$ (that is, $\tau v \sim 10^2$). We observe that the libration frequency of the spin is indeed $v$, and that the forced libration amplitude also matches the estimated value.

We can now estimate the correction that the forced libration introduces in the tidal deformation (Eq. (24)) as follows:

$$ \begin{aligned} \langle Z_t \rangle &\approx \mathcal{A}(\bar{x}_{[2,1]}Y_{[0,0]} + \bar{x}_{[2,0]}Y_{[0,-1]}) \\ &= \mathcal{A}e^{2i\theta_0}(\bar{x}_{[2,1]} - \bar{x}_{[2,0]}\bar{\alpha}_v) \\ &= \mathcal{A}e^{2i\theta_0}\bar{x}_{[2,1]} \left(1 - \frac{18 \mathcal{A}(n/v)^2 |x_{[2,0]}|^2}{1 + 18 \mathcal{A}(n/v)^2 |x_{[2,0]}|^2}\right). \end{aligned} \quad (53) $$

We thus obtain for the outer planet in the KOI-1599 system, $|\langle Z_t \rangle| \approx 3.6 \times 10^{-7}$, which is about half of the value without con- sidering the forced librations (Eq. (34)); this value corresponds to what is observed in the numerical simulations (Fig. 4). We conclude that the forced libration contributes to decreasing the tidal deformation (see also Delisle et al. 2017, Appendix D).

## 5. Tidal dissipation and capture probabilities

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

### 5.1. Tidal dissipation

Replacing the tidal deformation of the planet $Z_t$ (Eq. (21)) in the general expression of the torque (Eq. (2)), and averaging over the fundamental frequencies of the system we get

$$ \begin{aligned} \mathcal{T} &= \langle \mathcal{B} \mathfrak{I} (Z_t \times \mathbf{e}^{-2i\theta}) \rangle \\ &= \mathcal{A} \mathcal{B} \mathfrak{I} \left( \sum_{\mathbf{k},l,m} x_{\mathbf{k}+l} \bar{x}_{\mathbf{k}+m} y_{\mathbf{m}} \bar{y}_{l} \frac{1 + i(2\omega - \mathbf{k} \cdot \mathbf{v})\tau_e}{1 + i(2\omega - \mathbf{k} \cdot \mathbf{v})\tau} \right). \end{aligned} \quad (54) $$

In addition, assuming very small or no forced oscillations, that is, $\theta \approx \langle\theta\rangle = \theta_0 + \omega t$ (Eq. (17)), we have $y_0 = \mathbf{e}^{2i\theta_0}$ and $y_{k \neq 0} = 0$, so the average net torque can be simplified as

$$ \mathcal{T}(\omega) = -\mathcal{K} \sum_{\mathbf{k}} |x_{\mathbf{k}}|^2 \frac{(2\omega - \mathbf{k} \cdot \mathbf{v})\tau}{1 + (2\omega - \mathbf{k} \cdot \mathbf{v})^2 \tau^2}, \quad (55) $$

where

$$ \mathcal{K} = \mathcal{A}\mathcal{B}\left(1 - \frac{\tau_e}{\tau}\right) = k_f \frac{3GM^2R^5}{2Ca^6}\left(1 - \frac{\tau_e}{\tau}\right). \quad (56) $$

The initial rotation rate of the planets in unknown. A small number of large impacts at the end of the formation process can significantly change the spin rate or direction (Dones & Tremaine 1993; Kokubo & Ida 2007). However, since plan- ets contract from an initial larger protoplanet, in general it is expected that they rotate fast. Indeed, based on the solar sys- tem planets, MacDonald (1964) proposed an empirical relation $\omega_0 \propto m^{4/5}/R^2$, which leads to about 12 h for Earth-size planets. Therefore, in general we expect $\omega_0 \gg \mathbf{k} \cdot \mathbf{v}$, and since $\tau_e < \tau$, the initial rotation rate can only decrease (Eq. (55)).

As the planet arrives in a regime where $\omega \sim \mathbf{k} \cdot \mathbf{v}$ the situa- tion changes. Depending on $\tau$ and the amplitudes $|x_k|$, the torque may become positive, and hence create stable equilibria for the rotation rate. In Fig. 5 we draw the net torque as a function of $\tau$ for two different orbital forcing functions. The left-hand column represents the case of Mercury, for which $|x_k| = X_k^{-3.2}(e)$ with $e = 0.206$, while the right-hand column is for the outer planet in the KOI-1599 system, for which $|x_k|$ are given in Table 2.

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

### 5.2. Capture probabilities

The total torque on the rotation (Eq. (2)) can be written as

$$ \ddot{\theta} = \mathcal{B} \mathfrak{I} (Z_0 x \mathbf{e}^{-2i\theta}) + \mathcal{B} \mathfrak{I} (Z_t x \mathbf{e}^{-2i\theta}), \quad (57) $$

since $Z = Z_0 + Z_t$. The average of the second term is the tidal torque $\mathcal{T}(\omega)$ (Eq. (55)). The average of the first term is different from zero when the rotation rate is close to some resonant island, $\omega \approx v_r$, with $v_r = \mathbf{r} \cdot \mathbf{v}/2$ (Eq. (9)):

$$ \langle \mathcal{B} \mathfrak{I} (Z_0 x \mathbf{e}^{-2i\theta}) \rangle = \mathcal{B} \mathfrak{I} (Z_0 x_r \mathbf{e}^{-2i(\gamma-\phi)}), \quad (58) $$

where $\gamma = \theta - v_r t + \phi$, and $\phi = \frac{1}{2} \arg(Z_0 x_r)$ is a constant phase. Thus, the average of the total torque (Eq. (57)) becomes

$$ \dot{\gamma} = -\frac{\beta_r^2}{2} \sin 2\gamma + \mathcal{T}(v_r + \dot{\gamma}), \quad (59) $$

where $\beta_r^2 = 2\mathcal{B}|Z_0 x_r|$ (Eq. (11)). If $Z_0 = 0$, capture could only be possible at the equilibrium of the tidal torque, that is, when $\mathcal{T} = 0$. However, in a more general situation, the tidal torque can be counterbalanced by the conservative torque resulting from $\beta_r$. Therefore, even when $\mathcal{T} \neq 0$, each time the planet crosses a resonant island, there is some chance of being captured.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Normalized net torque $\mathcal{T}/\mathcal{K}$ (Eq. (55)) for different $\tau$ values and two orbital forcing functions: Mercury, for which $|x_k| = X_k^{-3,2}(e)$ and $e = 0.206$ (left column); the outer planet in the KOI-1599 system, for which the $|x_k|$ are given in Table 2 (right column).</figcaption></figure>

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

$$ \frac{dh}{dt} = 2 \mathcal{T} \frac{d\gamma}{dt}, \qquad (60) $$

where

$$ h = \dot{\gamma}^2 - \frac{\beta_r^2}{2} \cos 2\gamma \qquad (61) $$

is an integral of motion in the conservative problem (Sect. 3), which can be related to the total energy. The separatrix equation between dynamical regimes is $h = \beta_r^2/2$, where $h > \beta_r^2/2$ gives the trajectories in the circulation zone (outside the resonance) and $h < \beta_r^2/2$ the trajectories in the libration zone (captured in

resonance). The “energy” variation of the planet after a cycle around the resonance is then given by the function

$$ \Delta h(\gamma_1) = 2 \int_{\gamma_1}^{\gamma_2} \mathcal{T} \, d\gamma + 2 \int_{\gamma_2}^{\gamma_3} \mathcal{T} \, d\gamma, $$

where $\gamma_1$ is the $\gamma$ value when the planet crosses the separatrix between the circulation and libration zones. The values $\gamma_2$ and $\gamma_3$ are the first two following $\gamma$ values corresponding to $\dot{\gamma} = 0$. In order to be captured, after a cycle inside the resonance, the planet must remain within the libration zone, that is, the “energy” variation after a cycle must be negative $\Delta h(\gamma_1) < 0$. Assuming that $\gamma_1$ is uniformly distributed between $-\pi/2$ and $\pi/2$, capture inside the resonance occurs for “energy” variations within

$$ \{\mathcal{E} : \Delta h(-\pi/2) < \Delta h(\gamma_1) < 0\}, \qquad (62) $$

from a total of possibilities

$$ \{\mathcal{E}_T : \Delta h(-\pi/2) < \Delta h(\gamma_1) < \Delta h(\pi/2)\}, \qquad (63) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Capture probabilities in the presently observed 3/2 spin-orbit resonance for Mercury ($e = 0.206$ and $|Z_0| = 1.25 \times 10^{-5}$). The solid line gives the theoretical approximation given by expression (66), while the dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with $\theta_0$ differing by one degree.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Capture probabilities for the super-synchronous resonance $v_r = n + v/2$ for the outer planet in the KOI-1599 system with $|Z_0| = 10^{-7}$. The solid line gives the theoretical approximation given by expression (66) with all terms from Table 2, while the dashed line is obtained with only the first three dominating terms. The dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with $\theta_0$ differing by one degree.</figcaption></figure>

that is,

$$P_{\text{cap}} = \frac{\int_{E} d\gamma_1}{\int_{E_T} d\gamma_1} = \frac{\Delta h(-\frac{\pi}{2})}{\Delta h(-\frac{\pi}{2}) - \Delta h(\frac{\pi}{2})} \quad (64)$$

Usually the tidal torque is very weak, that is, $|\mathcal{T}(\dot{y})| \ll \beta_r^2$, which gives $\gamma_3 \simeq -\pi/2$, $\gamma_2 \simeq \pi/2$, and $h \simeq \beta_r^2/2$. Then, from expression (61) we have

$$\dot{y} = \operatorname{sign}(\dot{y})\beta_r \cos \gamma, \quad (65)$$

and replacing it in expression (64) we finally get

$$P_{\text{cap}} = 1 - \frac{T_{-}}{T_{+}}, \quad \text{with} \quad T_{\pm} = \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \mathcal{T}(v_r \pm \beta_r \cos \gamma) \, d\gamma. \quad (66)$$

This expression is valid for any kind of tidal torque. For the Maxwell rheology that we adopt in this paper, we use expression (55) for $\mathcal{T}$, which leads to

$$T_{\pm} = -\mathcal{K} \sum_{k} |x_k|^2 \left[ \frac{1}{\Delta_k^+} \ln \left( \frac{1 + \delta_k^+}{1 - \delta_k^+} \right) - \frac{1}{\Delta_k^-} \ln \left( \frac{1 + \delta_k^-}{1 - \delta_k^-} \right) \right], \quad (67)$$

with $\Delta_k^\pm = \sqrt{\eta^2 - (\lambda_k \pm i)^2}$, $\delta_k^\pm = (\pm(\eta - \lambda_k) - i)/\Delta_k^\pm$, $\eta = \pm 2\tau\beta_r$, and $\lambda_k = \tau(\mathbf{r} - \mathbf{k}) \cdot \mathbf{v}$. The $\pm$ in $\eta$ corresponds to $T_\pm$. Since both $\eta$ and $\lambda_k$ depend on $\tau$ and $|Z_0|$ (through $\beta_r$), the chances of capture in each resonant island change with these two parameters.

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

remain in the 3/2 spin-orbit equilibria even with $Z_0 = 0$ (Correia et al. 2014).

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

In Fig. 7, we additionally show the results from numerical simulations (Eqs. (1) and (15)). Contrary to the case of Mercury (Fig. 6), in this case the agreement between the theoretical and numerical capture probabilities is not so good. For $\tau n \ll 10^2$, this is because the theoretical probabilities (Eq. (66)) were estimated using the approximation (59), which assumes that the effect of the nearby resonances can be averaged to zero. However, although isolated resonant islands exist for $|Z_0| = 10^{-7}$, there is already a small chaotic zone around these islands (Fig. 3, bottom). Therefore, in the numerical simulations the rotation is chaotic for some time before it is stabilized in one of these islands (see Sect. 6) and the probabilities are altered. For $\tau n \gg 10^2$, in addition to the previous effect, there is a significant increase in the permanent tidal deformation $\langle Z_t \rangle > 10^{-7}$ (Eqs. (33) and (34)) that extends the size of the chaotic zone (Fig. 3, top). As a result, the rotation may remain permanently chaotic or evolve into some totally unexpected spin-orbit resonance (see next section for more details).
<table><caption>**Table 3.** Number of captures in each spin-orbit resonances (in percentage) for the KOI-1599 outer planet using different $\tau$ values.</caption><thead><tr><th>τn</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th><th>10<sup>-1</sup></th><th>10<sup>0</sup></th><th>10<sup>1</sup></th><th>10<sup>2</sup></th><th>10<sup>3</sup></th><th>10<sup>4</sup></th><th>10<sup>5</sup></th><th>10<sup>6</sup></th></tr></thead><tbody><tr><td>n + v</td><td>5.1</td><td>9.4</td><td>3.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>n + v/2</td><td>64.0</td><td>47.5</td><td>58.3</td><td>74.0</td><td>10.6</td><td>—</td><td>—</td><td>—</td><td>59.2</td><td>—</td><td>—</td></tr><tr><td>n</td><td>30.3</td><td>39.2</td><td>36.7</td><td>26.0</td><td>89.4</td><td>100</td><td>100</td><td>100</td><td>—</td><td>—</td><td>—</td></tr><tr><td>n - v/2</td><td>0.6</td><td>3.9</td><td>1.7</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>40.8</td><td>—</td><td>—</td></tr><tr><td>n - v</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>86.7</td></tr><tr><td>Chaotic</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>100</td><td>13.3</td></tr></tbody></table>

## 6. Complete evolution

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

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

In this section, we study the complete behavior of the outer planet in the KOI-1599 system (Table 1). This planet has an estimated mass of 4.3 $M_⊕$, that is, it can be classified as a super-Earth. An upper bound of ten Earth masses is commonly accepted to warrant that the planet is mostly rocky (e.g., Valencia et al. 2007). Therefore, the internal structure and rheology of this planet is likely similar to the terrestrial planets in the solar system. We adopt for the residual deformation $Z_0 = 1.0 \times 10^{-7}$. For the tidal deformation we use a Maxwell model (Eq. (15)), with $k_f = 0.9$, $\tau_e = 0$, and $C/mR^2 = 1/3$.

For rocky planets, such as the Earth and Mars, we have a dissipation factor $Q = 10$ (Dickey et al. 1994) and $Q = 80$ (Lainey et al. 2007), respectively. We then compute for the Earth $\tau = 1.6$ day, and for Mars $\tau = 14.7$ day. However, in the case of the Earth, the present Q-factor is dominated by the oceans, the Earth's solid body $Q$ is estimated to be 280 (Ray et al. 2001), which increases the relaxation time by more than one order of magnitude ($\tau = 46$ day). Although these values provide a good estimate for the average present dissipation ratios, they appear to be inconsistent with the observed deformation of the planets. Indeed, in the case of the Earth, the surface post-glacial rebound from the last glaciation about $10^4$ yr ago is still going on, suggesting that the Earth's mantle relaxation time is something like $\tau = 4400$ yr (Turcotte & Schubert 2002). Since $\tau$ can range from a few minutes up to thousands of years, we adopt $-6 \le \log_{10} \tau_{yr} \le 4$, which is equivalent to $-4 \le \log_{10} \tau n \le 6$.

In our numerical simulations we use the initial conditions from Table 1. The initial rotation period is set at 19.5 day, which corresponds to $\dot{\theta}/n \approx 1.05$. Since the libration frequency of the outer planet is $v/n \approx 0.0066$ (Table 2), the initial rotation rate

is completely outside the resonant area, which is bounded by $\dot{\theta}/n \le 1.01$ (Fig. 3). For any $\tau$ value, the rotation rate of the planet decreases owing to tidal effects until it approaches the region where multiple spin-orbit resonances are present. The shape of the tidal torque depends on the $\tau$ value (Fig. 5), so at this stage different scenarios can occur. In addition, capture in resonance is a stochastic process. Therefore, for each $\tau$ value we have run 180 simulations with $\theta_0$ differing by one degree. In Table 3 we list the final distribution of the spin in the different resonances for each $\tau$ value.

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

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

For $\tau n \gg 10^2$, the spin is in the “high frequency regime” ($\tau v \gg 1$). In this regime, the relaxation time is large and the planet is thus able to retain some tidal deformation $\langle Z_i \rangle \ne 0$ (Sect. 4). For the outer planet in the KOI-1599 system we have $A = 2.7 \times 10^{-6} \gg Z_0$, so for the main spin-orbit resonances the tidal deformation dominates over the permanent deformation. As a result, the phase space of the spin is modified with respect to the phase space drawn for $Z_0 = 10^{-7}$, it resembles more to that for $Z_0 \sim 10^{-6}$ (Fig. 3, top). We therefore expect chaotic motion for the spin in this regime. However, in our simulations we verify that, although chaotic rotation is indeed a possibility for some $\tau$ values, stable rotation is also possible in some asynchronous spin-orbit resonances (Table 3). One reason for this is that in
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Some examples of the spin evolution of the outer planet in the KOI-1599 system for $\tau n = 10^{-4}$ (low frequency regime). The blue lines give the position of the main spin-orbit resonances $n, n \pm v/2$, and $n \pm v$.</figcaption></figure>

this regime the tidal torque facilitates the capture in all these resonances (Fig. 5). The other reason is because the tidal deformation varies from resonance to resonance (Eq. (30)). For the synchronous resonance $|(Z_t)| \sim 10^{-6}$, so even if the rotation is somehow captured there, it becomes chaotic after some time. For the super-synchronous resonance $|(Z_t)| \approx 3.6 \times 10^{-7}$ (Eq. (53)), so it is not impossible to remain trapped there, provided that the spin finds a way in.

In Fig. 9 we show an example for $\tau n = 10^6$, which clarifies the behavior of the spin in the high frequency regime. In this figure we show the evolution of the rotation rate (top) together with the evolution of the shape (bottom). Before the spin reaches the resonant area, we have $\langle Z \rangle = Z_0 = 10^{-7}$, so in principle the rotation can be trapped in a stable spin-orbit resonance (Fig. 3, bottom). As soon as the spin enters the resonant area, the deformation increases and can experience values $|Z| > 10^{-6}$. Then, the rotation becomes chaotic (Fig. 3, top). However, just outside the chaotic zone there exists the sub-synchronous resonance with $\omega = n - v$. Somehow, after some time in the chaotic region, the spin finds a path into this stable region. Once there, the deformation of the planet decreases again to a value $|Z| \sim 10^{-7}$. Thus, the sub-synchronous resonance is no longer connected with the chaotic region and the spin can remain there forever (Fig. 3, bottom). In this example, it is remarkable that although the spin

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Spin evolution of the outer planet in the KOI-1599 system (top) and the corresponding deformation $|Z|$ (bottom) for $\tau n = 10^6$ (high frequency regime). The blue lines give the position of the main spin-orbit resonances $n, n \pm v/2$, and $n \pm v$.</figcaption></figure>

starts with an initial rotation rate faster than the synchronous rotation, it finally ends up with a rotation slower than that.

## 7. Discussion

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

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

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

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

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

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

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

## References

Adams, F. C., & Bloch, A. M. 2015, MNRAS, 446, 3676

Andrade, E. N. d. C. 1910, Proc. R. Soc. London A, 84, 1

Carpenter, R. L. 1964, AJ, 69, 2

Chirikov, B. V. 1979, Phys. Rep., 52, 263

Colombo, G. 1965, Nature, 208, 575

Colombo, G. 1966, AJ, 71, 891

Correia, A. C. M. 2015, A&A, 582, A69

Correia, A. C. M. 2018, Icarus, 305, 250

Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767

Correia, A. C. M., & Laskar, J. 2004, Nature, 429, 848

Correia, A. C. M., & Laskar, J. 2010, Icarus, 205, 338

Correia, A. C. M., & Robutel, P. 2013, ApJ, 779, 20

Correia, A. C. M., & Rodríguez, A. 2013, ApJ, 767, 128

Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271

Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50

Correia, A. C. M., Leleu, A., Rambaux, N., & Robutel, P. 2015, A&A, 580, L14

Darwin, G. H. 1879, Philos. Trans. R. Soc. London, 170, 1

Darwin, G. H. 1880, Philos. Trans. R. Soc. London, 171, 713

Delisle, J.-B., Correia, A. C. M., Leleu, A., & Robutel, P. 2017, A&A, 605, A37

Dickey, J. O., Bender, P. L., Faller, J. E., et al. 1994, Science, 265, 482

Dobrovolskis, A. R. 1980, Icarus, 41, 18

Dones, L., & Tremaine, S. 1993, Icarus, 103, 67

Efroimsky, M. 2012, Celest. Mech. Dyn. Astron., 112, 283

Genova, A., Iess, L., & Marabucci, M. 2013, Planet. Space Sci., 81, 55

Gold, T., & Soter, S. 1969, Icarus, 11, 356

Goldreich, P., & Peale, S. 1966, AJ, 71, 425

Goldreich, P., & Peale, S. 1967, AJ, 72, 662

Goldreich, P., & Soter, S. 1966, Icarus, 5, 375

Goldstein, R. M. 1964, AJ, 69, 12

Harbison, R. A., Thomas, P. C., & Nicholson, P. C. 2011, Celest. Mech. Dyn. Astron., 110, 1

Hay, H. C. F. C., & Matsuyama, I. 2019, ApJ, 875, 22

Henning, W. G., O'Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000

Henrard, J. 1993, in Dynamics Reported (New York: Springer Verlag), 117

Hughes, S. 1981, Celest. Mech., 25, 101

Hut, P. 1980, A&A, 92, 167

Kokubo, E., & Ida, S. 2007, ApJ, 671, 2082

Lainey, V., Dehant, V., & Pätzold, M. 2007, A&A, 465, 1075

Laskar, J. 1993, Physica D, 67, 257

Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632

Leleu, A. 2017, ArXiv e-prints [arXiv:1701.05585]

Leleu, A., Robutel, P., & Correia, A. C. M. 2016, Celest. Mech. Dyn. Astron., 125, 223

MacDonald, G. J. F. 1964, Revs. Geophys., 2, 467

Mignard, F. 1979, Moon Planets, 20, 301

Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (London: Taylor & Francis)

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press)

Noyelles, B., Frouard, J., Makarov, V. V., & Efroimsky, M. 2014, Icarus, 241, 26

Panichi, F., Migaszewski, C., & Goździewski, K. 2019, MNRAS, 485, 4601

Peale, S. J. 1969, AJ, 74, 483

Peale, S. J., Yseboodt, M., & Margot, J.-L. 2007, Icarus, 187, 365

Pettengill, G. H., & Dyce, R. B. 1965, Nature, 206, 1240

Ray, R. D., Eanes, R. J., & Lemoine, F. G. 2001, Geophys. J. Int., 144, 471

Showalter, M. R., & Hamilton, D. P. 2015, Nature, 522, 45

Singer, S. F. 1968, Geophys. J. R. Astron. Soc., 15, 205

Smith, B. A., Soderblom, L., Batson, R. M., et al. 1982, Science, 215, 504

Steffen, J. H., & Hwang, J. A. 2015, MNRAS, 448, 1956

Thomas, P. C., Armstrong, J. W., Asmar, S. W., et al. 2007, Nature, 448, 50

Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge: Cambridge University Press)

Valencia, D., Sasselov, D. D., & O'Connell, R. J. 2007, ApJ, 656, 545

Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137

Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants (Washington, D.C.: American Geophysical Union), 1