Secular spin-axis dynamics of exoplanets

TL;DR

This paper develops an analytical model for the long-term spin-axis dynamics of exoplanets, enabling quick exploration of their climate stability and obliquity evolution based on limited orbital data.

Contribution

It introduces a formula linking exoplanet spin-axis behavior to physical and dynamical parameters, with bounds for uncertain data, improving analysis efficiency.

Findings

Accurate retrieval of Solar System planets' chaotic zones.

Maximal chaotic region extent from minimal orbital data.

Classification of exoplanets outside major spin-orbit resonances.

Abstract

Context: Seasonal variations and climate stability of a planet are very sensitive to the planet obliquity and its evolution. This is of particular interest for the emergence and sustainability of land-based life, but orbital and rotational parameters of exoplanets are still poorly constrained. Numerical explorations usually realised in this situation are thus in heavy contrast with the uncertain nature of the available data. Aims: We aim to provide an analytical formulation of the long-term spin-axis dynamics of exoplanets, linking it directly to physical and dynamical parameters, but still giving precise quantitative results if the parameters are well known. Together with bounds for the poorly constrained parameters of exoplanets, this analysis is designed to allow a quick and straightforward exploration of the spin-axis dynamics. Methods: The long-term orbital solution is decomposed in quasi-periodic series and the spin-axis Hamiltonian is expanded in powers of eccentricity and inclination. Chaotic zones are measured by the resonance overlap criterion. Bounds for the poorly known parameters of exoplanets are obtained from physical grounds (rotational breakup) and dynamical considerations (equipartition of AMD). Results: This method gives accurate results when the orbital evolution is well known. The chaotic zones for planets of the Solar System can be retrieved in details from simple analytical formulas. For less constrained planetary systems, the maximal extent of the chaotic regions can be computed, requiring only the mass, the semi-major axis and the eccentricity of the planets present in the system. Additionally, some estimated bounds of the precession constant allow to classify which observed exoplanets are necessarily out of major spin-orbit secular resonances (unless the precession rate is affected by the presence of massive satellites).

Tables

Table 1. Table 1: Critical angle, location and half width of every second and third-order secular spin-orbit resonance.

Table 2. Table 2: Secular representation of the Earth orbital dynamics given by the Lagrange-Laplace theory.

Table 3. Table 3: Exoplanets from http://exoplanet.eu at date 2018-03-07 classified non-resonant with the method detailed in this study.

Table 4. Table 4: Quasi-periodic representation of the orbital dynamics of Mercury.

Table 5. Table 5: Quasi-periodic representation of the orbital dynamics of Venus.

Table 6. Table 6: Quasi-periodic representation of the orbital dynamics of the Earth.

Table 7. Table 7: Quasi-periodic representation of the orbital dynamics of Mars.

Equations

$$\mathcal{H}(X,-\psi ,t)$$

$$\displaystyle-\sqrt{1-X^{2}}\big{(}\mathcal{A}(t)\sin\psi+\mathcal{B}(t)\cos\psi\big{)}$$

$$+2X\mathcal{C}(t),$$

$$\alpha =\frac{3\mathcal{G}{m}_0}{2\omega a^3}\frac{2C-A-B}{2C}.$$

$$\left\{\begin{aligned} \mathcal{A}(t)&=\frac{2\big{(}\dot{q}+p\,\mathcal{C}(t)\big{)}}{\sqrt{1-p^{2}-q^{2}}}\\ \mathcal{B}(t)&=\frac{2\big{(}\dot{p}-q\,\mathcal{C}(t)\big{)}}{\sqrt{1-p^{2}-q^{2}}}\\ \end{aligned}\right.\hskip 14.22636pt,\hskip 14.22636pt\mathcal{C}(t)=q\dot{p}-p\dot{q},$$

$$\mathcal{A}$$

$$\mathcal{B}$$

$$e\exp (i\varpi )$$

$$\eta \exp (i\Omega )$$

$${\theta }_j(t)={\mu }_{j}t+{\theta }_j^{(0)}\text{and}{\varphi }_j(t)={\nu }_{j}t+{\varphi }_j^{(0)},$$

$$e^2=\sum _{j=1}^N{E}_j^2+2\sum _{j<k}^N{E}_j{E}_k\cos ({\theta }_j-{\theta }_k),$$

$$\mathcal{C}=\sum _{j=1}^M{\nu }_j{S}_j^2+\sum _{j<k}^M({\nu }_j+{\nu }_k)S_j{S}_k\cos ({\varphi }_j-{\varphi }_k),$$

$$\mathcal{A}\sin \psi +\mathcal{B}\cos \psi =2\sum _{j=1}^M{\nu }_j{S}_j\cos ({\varphi }_j+\psi )$$

$$+\sum _{j=1}^M\left[{\nu }_j{S}_j^3-2S_j\left(\sum _{k=1}^M{\nu }_k{S}_k^2\right)\right]\cos ({\varphi }_j+\psi )$$

$$-\sum _{j=1}^M\sum _{\begin{array}{c}k=1\\ k\ne j\end{array}}^M{\nu }_k{S}_j^2{S}_k\cos (2{\varphi }_j-{\varphi }_k+\psi )$$

$$-2\sum _{i=1}^M\sum _{\begin{array}{c}j<k\\ j,k\ne i\end{array}}^M{\nu }_i{S}_i{S}_j{S}_k\cos (-{\varphi }_i+{\varphi }_j+{\varphi }_k+\psi )$$

$$+\mathcal{O}({\eta }^4).$$

$$\mathcal{H}={\mathcal{H}}_0+\epsilon {\mathcal{H}}_1+{\epsilon }^2{\mathcal{H}}_2+{\epsilon }^3{\mathcal{H}}_3+\mathcal{O}({\epsilon }^4),$$

$${\mathcal{H}}_0(X,\Theta ,\Phi )=-\frac{\alpha }{2}{X}^2+\sum _{j=1}^N{\mu }_j{\Theta }_j+\sum _{j=1}^M{\nu }_j{\Phi }_j,$$

$$\epsilon {\mathcal{H}}_1(X,-\psi ,\varphi )=-2\sqrt{1-X^2}\sum _{j=1}^M{\nu }_j{S}_j\cos ({\varphi }_j+\psi ),$$

$${\epsilon }^2{\mathcal{H}}_2(X,\theta ,\varphi )$$

$$-\frac{3\alpha }{2}{X}^2\sum _{j<k}^N{E}_j{E}_k\cos ({\theta }_j-{\theta }_k)$$

$$+2X\sum _{j<k}^M({\nu }_j+{\nu }_k)S_j{S}_k\cos ({\varphi }_j-{\varphi }_k),$$

$${\epsilon }^3{\mathcal{H}}_3(X,-\psi ,\varphi )=$$

$$-\sqrt{1-X^2}\sum _{j=1}^M\left[{\nu }_j{S}_j^3-2S_j\left(\sum _{k=1}^M{\nu }_k{S}_k^2\right)\right]\cos ({\varphi }_j+\psi )$$

$$+\sqrt{1-X^2}\sum _{j=1}^M\sum _{\begin{array}{c}k=1\\ k\ne j\end{array}}^M{\nu }_k{S}_j^2{S}_k\cos (2{\varphi }_j-{\varphi }_k+\psi )$$

$$+2\sqrt{1-X^2}\sum _{i=1}^M\sum _{\begin{array}{c}j<k\\ j,k\ne i\end{array}}^M{\nu }_i{S}_i{S}_j{S}_k\cos (-{\varphi }_i+{\varphi }_j+{\varphi }_k+\psi ).$$

$$\{\begin{array}{rl}\sigma & =\psi +{\varphi }_p\\ \xi & =-{\varphi }_p\end{array}\text{and}\{\begin{array}{rl}\Sigma & =-X\\ \Xi & =-X-{\Phi }_p.\end{array}$$

$$\mathcal{F}(\Sigma ,\sigma )=-\frac{1}{2}\mathfrak{a}\alpha {\Sigma }^2+\mathfrak{b}\Sigma +\mathfrak{c}\sqrt{1-{\Sigma }^2}\cos \sigma ,$$

$$\mathfrak{a}$$

$$\mathfrak{b}$$

$$\mathfrak{c}$$

$$\mathcal{F}(\Sigma ,\sigma )=-\frac{1}{2}{\Sigma }^2+\gamma \Sigma +\beta \sqrt{1-{\Sigma }^2}\cos \sigma ,$$

$$\gamma =\frac{\mathfrak{b}}{\mathfrak{a}\alpha }\text{and}\beta =\frac{\mathfrak{c}}{\mathfrak{a}\alpha }.$$

$${\mathcal{C}}_1=\left\{\gamma ,\beta ⩾0:{\gamma }^{2/3}+{\beta }^{2/3}=1\right\}.$$

$$\{\begin{array}{rl}{\Sigma }_a& \in [0,+1],{\sigma }_a=0\to \text{ elliptic (A,B,C,D)}\\ {\Sigma }_b& \in [-1,0],{\sigma }_b=\pi \to \text{ elliptic (A,B,C,D)}\\ {\Sigma }_c& \in [0,+1],{\sigma }_c=\pi \to \text{ hyperbolic (A,B,C)}\\ {\Sigma }_d& \in [0,+1],{\sigma }_d=\pi \to \text{ elliptic (A,B,C)}.\end{array}$$

$$\text{for }\beta ⩽1:\{\begin{array}{rl}& \underset{\gamma \to 0}{\lim }{\Sigma }_a=\underset{\gamma \to 0}{\lim }{\Sigma }_c=0\\ & \underset{\gamma \to 0}{\lim }{\Sigma }_d=-\underset{\gamma \to 0}{\lim }{\Sigma }_b=\sqrt{1-{\beta }^2}\end{array}$$

$$\text{for }\beta ⩾1:\underset{\gamma \to 0}{\lim }{\Sigma }_a=\underset{\gamma \to 0}{\lim }{\Sigma }_b=0,$$

Full text

*Aims. *We aim to provide an analytical formulation of the long-term spin-axis dynamics of exoplanets, linking it directly to physical and dynamical parameters, but still giving precise quantitative results if the parameters are well known. Together with bounds for the poorly constrained parameters of exoplanets, this analysis is designed to allow a quick and straightforward exploration of the spin-axis dynamics.

*Methods. *The long-term orbital solution is decomposed in quasi-periodic series and the spin-axis Hamiltonian is expanded in powers of eccentricity and inclination. Chaotic zones are measured by the resonance overlap criterion. Bounds for the poorly known parameters of exoplanets are obtained from physical grounds (rotational breakup) and dynamical considerations (equipartition of AMD).

*Results. *This method gives accurate results when the orbital evolution is well known. The chaotic zones for planets of the Solar System can be retrieved in details from simple analytical formulas. For less constrained planetary systems, the maximal extent of the chaotic regions can be computed, requiring only the mass, the semi-major axis and the eccentricity of the planets present in the system. Additionally, some estimated bounds of the precession constant allow to classify which observed exoplanets are necessarily out of major spin-orbit secular resonances (unless the precession rate is affected by the presence of massive satellites).

1 Introduction

From the works by Laskar & Robutel (1993) and Laskar et al. (1993b), we know that the long-term dynamics of the terrestrial planets of the Solar System feature wide chaotic regions allowing large variations of their obliquity. In particular, Mars is currently in a chaotic region extending from $0^{\text{o}}$ to $60^{\text{o}}$ obliquity, whereas the Earth is located in a stable region thanks to the presence of the Moon, resulting in obliquity variations of only a few degrees. Subsequent studies detailed both the past and future spin-axis evolution the Earth (Néron de Surgy & Laskar, 1997; Laskar et al., 2004b; Li & Batygin, 2014b), Venus (Correia et al., 2003; Correia & Laskar, 2003) and Mars (Laskar et al., 2004a). On the other hand, paleorecords on Earth show that even the very slight variations of its orbit and obliquity led to major climate changes (e.g. Weertman, 1976; Hays et al., 1976). This implies that life on Earth would be very different to what it is now if the Earth had evolved in a large chaotic zone, as it would have without the stabilising effect of the Moon. This conclusion is reinforced by the fact that high-obliquity planets undergo severe seasonal variations (e.g. Spiegel et al., 2009) even for a stable obliquity (but conditions suitable for the very emergence of life could still be achieved: in some extreme cases, the amount of liquid water on the surface may even be favoured by large obliquity variations, as reported by Armstrong et al., 2014).

As the formation of the Moon is thought to have resulted from an accidental collision event (e.g. Hartmann & Davis, 1975; Canup & Asphaug, 2001; Lock et al., 2018), a “moonless” Earth was a possible (or even likely) outcome of the planetary formation process. In the broader context of exoplanets, the dynamics of such a moonless Earth is thus more than of academic interest. This motivated further works about the structure of the chaotic region (Li & Batygin, 2014a) and some additional numerical studies (Lissauer et al., 2012).

When it comes to exoplanets, we must face the problem of the incomplete and imprecise nature of both dynamical and physical data. Except from very favourable cases, like a fast precession and an important flattening inducing detectable transit depth modulations (Carter & Winn, 2010; Correia, 2014), the spin orientation and the flattening of exoplanets are far from being reachable by observations. They must therefore be taken as completely free parameters, in the spirit of the work by Laskar & Robutel, 1993 for the Solar System planets. Moreover, the orbital properties of exoplanets are not well known either, especially the respective orientations of the orbits (including the mutual inclinations, which play a crucial role in the spin-axis dynamics). Several authors tackled this problem already (e.g. Brasser et al., 2014; Deitrick et al., 2018; Shan & Li, 2018): they used numerical integrations of the planetary system in order to build the time-dependent perturbation of the spin axis. Applied to extrasolar planets, this method requires to choose a nominal value for both the rotational parameters and the unknown orbital elements. The parameter space to be explored is thus very wide, so that even elaborate numerical explorations require some degree of arbitrariness. Consequently, the use of numerical integrations at this stage could appear a bit in contradiction with the very incomplete nature of the data. However, the secular problem (averaged over rotational and orbital motions) is not as complex as it could appear. Some hints about a possible analytical treatment were actually given by Laskar (1996) and partially exploited by Atobe et al. (2004), Li & Batygin (2014a) and Shan & Li (2018). In the exoplanetary case, a refined analytical theory would be very convenient, since it would give in a direct way the sensibility of the spin-axis dynamics to the various known and unknown parameters, instead of giving a list of possible outcomes. Associated with some bounds for the unknown parameters, such a theory would give a clear range for these possible outcomes.

This was the approach used by Atobe et al. (2004), with the aim of finding the probability of small obliquity variations for hypothetical terrestrial planets in the habitable zone of known exoplanetary systems. Their analytical developments, though, were limited to the lowest-order approximations (their parameter space was indeed enormous since, in this case, the planet itself was hypothetical). The recent work by Shan & Li (2018) also contains an analytical part applied to the particular case of exoplanets Kepler-62f and Kepler-186f. This time, their calculations were mostly designed to precise and explain numerical results, so they did not try to bring any substantial improvement to the theory by Atobe et al. (2004).

In this context, the goal of this article is twofold: i) provide a general analytical formalism for studying the long-term spin-axis dynamics of (exo)planets and ii) clarify what kind of information about the spin can be obtained from typical observed exoplanetary systems, that is, with numerous unknown physical and dynamical parameters. In particular, it is crucial for future studies to have a simple way to classify the observed exoplanets according to the characteristics of their spin dynamics. This would allow to determine which exoplanets are worth to be studied in more details (in particular if a complex chaotic spin dynamics is expected) and which ones have necessarily a very simple spin dynamics, making unnecessary any further numerical or analytical study. Such a general analysis will give both a qualitative view of the system if it is poorly known (in the continuation of Atobe et al., 2004), and a quantitative description of the dynamics if it is well known (as an analytical counterpart of Laskar & Robutel, 1993).

This article is organised as follows: Sect. 2 recalls the Hamiltonian of the secular spin-axis dynamics and shows how it can be expanded in terms of the orbital motion parameters. The secular resonances at all orders can then be isolated and used to delimit the chaotic regions. Then, Sect. 3 shows how an incomplete set of orbital elements can still be used to constrain the orbital solution of an exoplanet. Combined with the rotational breakup limit, it allows to make a preliminary classification of the “non-resonant” exoplanets, for which no chaos can appear and the obliquity variations are constrained by an analytical bound.

2 Analytical model of the long-term spin dynamics

2.1 Development of the Hamiltonian

Let us consider a system composed of a star and several planets. We study the rotational dynamics of one planet among them. For now, we consider that this planet is far from any spin-orbit resonance. Considering only the lowest-order term of the torque from the star expanded in Legendre polynomials, the Hamiltonian of rotation averaged over orbital and rotational motions is given for instance by Laskar & Robutel (1993) and detailed by Néron de Surgy & Laskar (1997). It can be written

[TABLE]

where the conjugate coordinates are $X$ (cosine of obliquity) and $-\psi$ (minus the precession angle). The quantity $\alpha$ is called the “precession constant” (contrary to previous studies, we prefer to exclude here the eccentricity $e$ appearing in denominator from the definition of $\alpha$). Following the derivation proposed by Néron de Surgy & Laskar (1997), we obtain

[TABLE]

In this expression, $\mathcal{G}$ is the gravitational constant; $m_{0}$ is the mass of the star; $a$ is the semi-major axis of the planet in orbit around the star; $\omega$ is its spin angular velocity, and $A\leqslant B\leqslant C$ are its momenta of inertia. The Hamiltonian (1) depends explicitly on time $t$ through the eccentricity $e$ and the functions

[TABLE]

in which $q=\sin(I/2)\cos\Omega$ and $p=\sin(I/2)\sin\Omega$, where $I$ and $\Omega$ are respectively the orbital inclination and the longitude of ascending node of the planet. In the following, we will write $\eta\equiv\sin(I/2)$. One can note that if there is only one planet in the system (two-body problem), the obliquity is constant and the precession angle circulates with constant angular velocity $\alpha X/(1-e^{2})^{3/2}$.

Let us suppose that the eccentricity and the inclination of the planet are small, such that we can develop the Hamiltonian in series of $e$ and $\eta$. In the following, we present the terms up to order $3$, but the method presented here can be generalised to any order (as we will see, the third order is the first one at which the eccentricity begins to play a substantial role). Using the fact that $\mathcal{C}=\eta^{2}\dot{\Omega}=\mathcal{O}(\eta^{2})$, we obtain

[TABLE]

Let us now suppose that the secular orbital dynamics of the planet, resulting from the perturbations by the other planets, is quasi-periodic. This amounts to considering that the chaos present in the orbital secular system acts on a much larger timescale than the spin dynamics under study. As we will see below, this holds very well for the Solar System (this methodology was first proposed by Laskar, 1996; it is used for instance by Li & Batygin, 2014a). In this case, we can write

[TABLE]

where $\varpi$ is the longitude of pericentre of the planet in orbit around the star. The angles $\theta_{j}$ and $\phi_{j}$ evolve linearly with frequencies $\mu_{j}$ and $\nu_{j}$, that is,

[TABLE]

whereas the amplitudes $E_{j}$ and $S_{j}$ are real constants of order $e$ and $\eta$ or smaller. Such series can be obtained either from analytical theories or from frequency analysis of numerical solutions (Laskar, 1988, 1990). In a general integrable case, $\mu_{j}$ and $\nu_{j}$ are integer combinations of the fundamental frequencies of the orbital dynamics (usually noted $g_{k}$ and $s_{k}$), and the series contain an infinite number of terms. Arranging the terms by decreasing amplitude, we consider here a truncation with $N$ terms for the eccentricity and $M$ terms for the inclination. We get then

[TABLE]

[TABLE]

and from (4),

[TABLE]

In order to obtain an autonomous Hamiltonian, let us introduce the momenta $\Theta_{j}$ and $\Phi_{j}$ conjugate to $\theta_{j}$ and $\phi_{j}$. The system has now $N+M+1$ degrees of freedom. The new Hamiltonian (that we still denote $\mathcal{H}$) can be written

[TABLE]

where we suppose that $\mathcal{O}(e)=\mathcal{O}(\eta)=\mathcal{O}(\varepsilon)$. The different parts are respectively

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

We note that there can be no resonance among the angles $\phi_{j}$ and $\theta_{j}$ because they come from the quasi-periodic solution of the orbital dynamics. By definition, they are thus already “integrated”.

2.2 One perturbing term: Colombo’s top

From (11-12), we conclude that at lowest-order to the perturbation, resonant angles can only be of the form $\sigma=\phi_{j}+\psi$. Let us consider a single resonance with the term $j=p$. We introduce the resonant canonical coordinates by the linear transformation

[TABLE]

Assuming that the system is far from any other resonance, the long-term dynamics at first order to the perturbation is given by averaging the Hamiltonian over all angles but $\sigma$. Dropping the constant terms, we get

[TABLE]

with

[TABLE]

As shown in Appendix A, the Hamiltonian has the same form in the case of a $1$:$1$ spin-orbit resonance, with though a slightly different expression of the coefficients. This would thus only shift a bit the position of the secular resonances considered here. In this case, the tidal damping (responsible for this capture in spin-orbit resonance) is supposed to act on a much larger timescale than the spin dynamics studied here, such that our approach still holds. Finally, applying the modified time $\mathrm{d}\tau=\mathfrak{a}\,\alpha\,\mathrm{d}t$ to (16), we obtain the following Hamiltonian (that we still denote $\mathcal{F}$):

[TABLE]

where

[TABLE]

The dynamical system with Hamiltonian function (18) is well known. It was thoroughly studied by Henrard & Murigande (1987), who called it “Colombo’s top” in memory of Colombo (1966). In the following, we detail its characteristics of interest here, in terms of the two constant parameters $\gamma$ and $\beta$.

First of all, it is enough to study the case $\gamma\geqslant 0$ and $\beta\geqslant 0$ since we get the negative cases by the transformations $\Sigma\rightarrow-\Sigma$ and $\sigma\rightarrow\sigma+\pi$, respectively. When studying the equilibrium points of the system (Henrard & Murigande, 1987, or Appendix B.1 and B.2), we find that the phase space can have two different geometries according to the value of $\gamma$ and $\beta$ (see Fig. 2 for the phase portraits). The boundary between these two regions of the parameter space is the curve

[TABLE]

Below $\mathscr{C}_{1}$ (regions A,B,C of Fig. 1), the dynamical system (18) has four equilibrium points, that we denote $(a,b,c,d)$. The equilibrium points $c$ and $d$ merge along the curve $\mathscr{C}_{1}$, and disappear above it (region D of Fig. 1). Their respective positions are

[TABLE]

The values $\Sigma_{a,b,c,d}$ have explicit expressions in terms of $\gamma$ and $\beta$, as given in Appendix B.1 (they correspond to the different roots of a quartic equation). Since the resonant angle $\psi+\phi_{p}$ is equal to [math] or $\pi$ for each of these equilibrium points, they all correspond to configurations where the spin axis, the normal to the orbit (reduced to its $p$th harmonic) and the normal to the reference plane are in the same plane. As such, they are commonly called “Cassini’s states” and labelled $(1,2,3,4)$ after Peale (1969), corresponding to the equilibrium points $(c,a,b,d)$.

We note the limiting cases

[TABLE]

and

[TABLE]

We are now interested in the width of the resonance, that is the interval of $\Sigma$ enclosed in the separatrix emerging from the hyperbolic fixed point $c$ and containing the fixed point $a$. A pendulum approximation can be obtained for small values of $\beta$ (as used by Atobe et al., 2004 or Li & Batygin, 2014a), but this approximation is no longer valid when $\beta$ grows. Since an analytical expression can be derived even in the general case, we will use it here. The computations (Appendix B.5) lead to the following extreme values of $\Sigma$ spanned by the resonance:

[TABLE]

They are defined whenever $\Sigma_{c}$ itself is defined (regions A,B,C of Fig. 1). At this point, it is important to note that the coordinates $(\Sigma,\sigma)$ are singular at $\Sigma=\pm 1$ since the problem actually takes place on the sphere (Henrard & Murigande, 1987). We must thus study carefully the meaning of the limits (24) when one of them crosses $\pm 1$. This leads to two other limits in the parameter space, as the curves

[TABLE]

(see Appendix B.3-B.4), delimiting the regions A-B and B-C of Fig. 1, respectively. We note that $\mathscr{C}_{1}$ and $\mathscr{C}_{2}$ intersect at $(\gamma;\beta)=(1\,;0)$, $\mathscr{C}_{1}$ and $\mathscr{C}_{3}$ intersect at $(1/8\,;3\sqrt{3}/8)$, and $\mathscr{C}_{2}$ and $\mathscr{C}_{3}$ intersect at $(0\,;1/2)$. Contrary to $\mathscr{C}_{1}$, the boundaries $\mathscr{C}_{2}$ and $\mathscr{C}_{3}$ do not correspond to actual bifurcations of the dynamical system, but only to the limits where the resonant island contains the north pole of the sphere (region B) and both poles of the sphere (region C). Hence, the resonance lies in $[\Sigma_{-};\Sigma_{+}]$ in region A; in $[\Sigma_{-};+1]$ in region B; and in $[-1;+1]$ in region C (see Fig. 2). In region D, there is no more resonance (the separatrix disappears), but the obliquity can still vary substantially. In Fig. 1, the colour shades in region D show the oscillation amplitude of the obliquity as the trajectory passes through $\Sigma=1$.

When $\beta\rightarrow 0$, the width of the island tends to [math] and all the level curves in the $(\sigma,\Sigma)$ plane tend to be horizontal. In this limit, the resonance width is small and almost independent of $\gamma$ (pendulum approximation). When $\gamma\rightarrow 0$, the phase portrait in the $(\sigma,\Sigma)$ plane tends to be symmetric with respect to the line $\Sigma=0$.

From (19), we note that

[TABLE]

Remembering that $S_{p}$ is the amplitude of a given term in the inclination series (5), this means that the parameter region above the line $\beta=2\gamma$ in Fig. 1 (that is, $|S_{p}|=1$) cannot be reached in this problem.

As a simple rule of thumb, one can consider that $\beta$ controls the resonance width, and that $\gamma$ controls the location of the resonance centre (see the Hamiltonian at Eq. 18). From Eq. (26), we know that $\beta$ is proportional to the amplitude $S_{p}$. Therefore, the resonance widths are larger in hot planetary systems, for which the mutual inclinations are large. This was exploited by Boué & Laskar (2010) in their scenario for tilting the spin-axis of Uranus. On the contrary, the secular system cannot produce any obliquity variation if the mutual inclinations are exactly zero. One must keep in mind that $\beta$ depends on the precession constant $\alpha$ (Eq. 2) as well, so that “small” mutual inclinations do not guarantee that the resonances are thin. For the terrestrial planets of the Solar System, some values of $\alpha$ produce first-order resonances larger than $70^{\text{o}}$, even if the mutual orbital inclinations are modest (see Sect. 2.3).

Contrary to $\beta$, the magnitude of $\gamma$ cannot be easily traced back from the planetary architecture: the first-order secular spin-orbit resonances of any planet can be located anywhere between $0^{\text{o}}$ and $180^{\text{o}}$ of obliquity. The large majority of them actually lie in $[0^{\text{o}};90^{\text{o}}]$ because most of the frequencies $\nu_{p}$ are negative (the explanation for this property is given in Sect. 3.1).

2.3 Overlap of first-order resonances

Going back to the full Hamiltonian (10), the main chaotic regions of the system can be estimated as the overlap of the first-order resonances taken separately (Chirikov’s criterion). With the analytical expression of their respective widths given in Sect. 2.2, the computation of the overlapping regions is straightforward for any quasi-periodic representation (5) of the orbital motion.

When comparing two secular terms, we note that no chaotic zone can form if at least one of them is in region D, because there is no separatrix in D. This is well verified by Poincaré sections. A direct consequence of this property is that if all the terms of the quasi-periodic series are in region D, the secular dynamics of the obliquity cannot be chaotic at first order. Moreover, if all the $S_{i}$ coefficients are small (low inclination regime), so are the corresponding $\beta$ coefficients, resulting in virtually no secular variation of the obliquity (see the shades of grey in Fig. 1). Actually, a no-chaos criterion can be obtained even if the amplitudes $S_{i}$ are not known: we just have to check that $\gamma_{i}\approx\nu_{i}/\alpha>1\ \forall\ i=1...M$. On the contrary, if $\gamma_{i}<1$ for at least two $i$, the existence of a chaotic region is possible, but not guaranteed. These properties will be discussed further in Sect. 3.4.

This method can be easily checked in the case of the Solar System, since the secular spin dynamics of all planets have been studied in the literature. Moreover, a very accurate quasi-periodic approximation of the orbital dynamics can be obtained, since the properties and initial conditions of the planets are very well known. As an example, the upper row of Fig. 3 shows the widths and overlaps of first-order resonances for the terrestrial planets (light and dark-red regions). This figure was produced by applying the previous analytical formulas to the orbital series of Laskar (1990), containing more than 50 terms in both eccentricity and inclination (see Appendix F). We note that most of the first-order resonances overlap (there are almost no light-red regions). This results in wide chaotic zones even if the individual amplitudes $S_{i}$ are small. However, the full extent of the chaotic regions given by the frequency analysis (Fig. 3, second and third rows) cannot be retrieved by only considering first-order resonances. The following section is thus dedicated to second and third-order resonances.

2.4 Higher-order resonances

Outside of first-order resonances, we can use a canonical change of coordinates close to identity in order to suppress the angular dependency at first order (as already used in a similar context by Li & Batygin, 2014a). Let us consider an intermediary Hamiltonian $\mathcal{X}=\varepsilon\mathcal{X}_{1}$, such that the current coordinates are obtained from the new ones through its flow at time $1$. The Hamiltonian in the new coordinates is then

[TABLE]

where

[TABLE]

In these expressions, Poisson’s brackets are defined as

[TABLE]

where the pairs $(p_{j},q_{j})$ are conjugate variables, $p_{j}$ being the momentum and $q_{j}$ the coordinate. In order to suppress the angular dependency at order 1, the Hamiltonian $\mathcal{X}$ must fulfil the homological equation

[TABLE]

in which $\overline{\mathcal{H}}_{1}$ is the 0th-order term of the multidimensional Fourier decomposition of $\mathcal{H}_{1}$ (average of $\mathcal{H}_{1}$ over all angles), which is here equal to zero. By matching the terms of the Fourier decomposition of $\mathcal{H}_{1}$ and $\mathcal{X}_{1}$ one by one, the solution to the homological equation is

[TABLE]

Injecting this function into the expression of the new Hamiltonian (28), we get $\tilde{\mathcal{H}}_{1}=0$ as required, and the second order term $\tilde{\mathcal{H}}_{2}$ is given in Appendix C. The only possible resonant angles at second order have the form $\sigma=\phi_{j}+\phi_{k}+2\psi$. The width of second and higher order resonances is quite small, so that their separatrices can be computed assuming that $X$ is near the exact resonance (pendulum approximation). Accordingly, the centre and half width of the second-order resonances are computed in Appendix C and given in the first line of Table 1.

Outside of both first-order and second-order resonances, the same method can be used to compute the location and width of third-order resonances. An intermediary Hamiltonian of the form $\mathcal{X}=\varepsilon\mathcal{X}_{1}+\varepsilon^{2}\mathcal{X}_{2}$ is used, in which $\mathcal{X}_{2}$ must satisfy a second homological equation (see Appendix D). The possible resonant angles, as well as the centre and half width of all the third-order resonances are given in Table 1. As before, the same method can be used in the case of synchronous rotation, by adding the term $-\alpha_{r}(1+X)^{2}/2$ to $\mathcal{H}_{0}$ (Appendix A). This would only slightly shift the resonances.

With these values, it is straightforward to compute the overlap regions of every possible second and third-order resonances. Resonances of order 2 or 3 with a centre located inside a resonance of lower order are not considered (in other words, low-order resonances are supposed unaffected by higher-order ones). The result is shown in the top row of Fig. 3 for the inner Solar System (blue and green zones). We obtain a much better match with the frequency map analysis, showing the importance of high-order resonances in this context. Indeed, numerous frequencies of the quasi-periodic representation are quite close to each other, which implies that the corresponding resonances overlap massively. Hence, even if the second and third-order resonances are thin, most of them are all located one after another, resulting in large chaotic zones. The chaotic diffusion is though slower for higher-order resonances. In the real Solar System, in which the secular orbital frequencies are actually not fixed but vary slowly (Laskar, 1990), the diffusion of the obliquity will besides be facilitated by small modulations of the resonances locations.

We note that third-order resonances include terms mixing both eccentricity and inclination (last line of Table 1). The existence of these terms shows that the $s_{6}+g_{5}-g_{6}$ resonance, which is known to play an important role in the future obliquity evolution of the Earth (Laskar et al., 1993a, 2004b), has two different origins: one is a first-order resonance with $\nu_{23}$, and the other is a third-order resonance between $\mu_{1}$, $\mu_{10}$ and $\nu_{12}$ (the numbering refers to Laskar et al., 2004b). The amplitude of the first resonance, which is the one emphasised in the literature, is larger by a factor thirty. This resonance is present in the synthetic representation used in this paper (see the term with frequency $-50.30212\ ^{\prime\prime}$/yr in Tables 4 to 7). In the top row of Fig. 3, it appears as a thin isolated resonance (upper-left corner of the graphs).

Finally, since tidal dissipations are much more efficient in decreasing the eccentricity than the inclination, one can imagine a planet with an initially chaotic obliquity wandering in a mixed-type-resonance overlap region, becoming frozen out of resonance when the amplitudes $E_{i}$ decrease due to tidal dissipation. As shown by Laskar et al., 2012, the eccentricity amplitudes of all the planetary system can be damped even if only one planet dissipates energy with the star. Hence the eccentricity modes of an external planet can be damped even if it is not itself subject to tidal dissipation.

3 Application to exoplanetary systems

In the previous sections, we saw that in the low-eccentricity and low-inclination regime, the long-term rotational dynamics of planets can be studied very efficiently by a simple analytical model. However, even if numerous exoplanetary systems are known nowadays, most of the information required to characterise the rotation of their planets remains poorly constrained. This information can be split into two groups: i) the orbital dynamics (amplitudes and frequencies of the quasi-periodic representation) and ii) the rotation parameters ($\alpha$ coefficient). In this section, we will see how these quantities can be estimated from physical and dynamical arguments, even with scarce data.

3.1 The Lagrange-Laplace system

Regarding the long-term orbital dynamics, one can use nominal orbital elements (either best-fit or assumed ones) and integrate numerically the equations of motion. The so-obtained solution can then be used directly (as did for instance Brasser et al., 2014, and Deitrick et al., 2018), or put in the form of quasi-periodic series and used as shown above. However, this method puts a heavy contrast between the very uncertain nature of the orbital elements used and the refined numerical solution applied. Actually, we will see that at this level of precision, the Lagrange-Laplace system is already a good-enough approximation of the orbital dynamics, up to moderate eccentricities and inclinations (and without strong effects coming from mean-motion resonances). It was used for the same purpose by Atobe et al. (2004) in the case of a massless hypothetical terrestrial planet.

The Lagrange-Laplace system is the lowest-order model of the long-term orbital dynamics: it uses a development of the Hamiltonian at second order of the eccentricities and inclinations, which is itself averaged over the fast angles (secular model at first order to the mutual perturbations). Let us write

[TABLE]

where the index $k=1,2...N$ represents a given planet of the system. Writing $\mathbf{z}$ and $\boldsymbol{\zeta}$ the vectors of all $z_{k}$ and $\zeta_{k}$, the equations of motion in the Lagrange-Laplace approximation are

[TABLE]

where the real matrices $A$ and $B$ are only functions of the masses and semi-major axes. They can be retrieved from the lowest-order terms in eccentricity and inclination of the orbital Hamiltonian. Organising the planets by increasing semi-major axes, we get from Laskar & Robutel (1995):

[TABLE]

and

[TABLE]

in which $n_{j}^{2}a_{j}^{3}=\mathcal{G}(m_{0}+m_{j})$, and the functions $C_{2}(\alpha)$ and $C_{3}(\alpha)$ are expressed in terms of the Laplace coefficients $b_{s}^{(k)}$:

[TABLE]

(see Laskar & Robutel, 1995, Laskar et al., 2012 or Murray & Dermott, 1999). The equations of motion for $\mathbf{z}$ and $\boldsymbol{\zeta}$ are decoupled and linear, such that the solution can be obtained by diagonalising the matrices $A$ and $B$. Its expression can be taken directly from Laskar et al. (2012): it has the form of quasi-periodic series (5) as required by our model. The frequencies $g_{k}$ and $s_{k}$ are the eigenvalues of $A$ and $B$, and the amplitude of the term $k$ for the planet $j$ is

[TABLE]

where $(P,Q)$ are the matrices composed of the eigenvectors of $(A,B)$, the matrices $(P^{-1},Q^{-1})$ are their inverses, and $z_{i}(0)$, $\zeta_{i}(0)$ are the initial conditions of planet $i$. In this case, we note that there is a single term for each proper frequency of the system; the frequencies $\mu_{j}$ and $\nu_{j}$ of the quasi-periodic representation are thus directly equal to one of the $g_{k}$ and $s_{k}$, respectively.

From the conservation of total orbital angular momentum, one of the inclination proper frequencies $s_{k}$ is identically equal to zero (matrix $B$ has rank deficiency of order 1). The inclination series have thus $M=N-1$ terms, whereas the eccentricity series have $N$ terms. Moreover, all the inclination proper frequencies are negative (this can be shown from the Geršgorin circles theorem, see Appendix E).

The result, in terms of chaotic zones for the spin dynamics, is shown in Fig. 3, bottom row. Although the match with the numerical maps (Fig. 3, second and third rows) is not as good as when we used the synthetic representation of the orbital dynamics (Fig. 3, top row), the estimate obtained is still remarkably good considering the uncertainties of the elements of an exoplanetary system222Some subtle dynamical effects are not reproduced by the Lagrange-Laplace system, like the $s_{6}+g_{5}-g_{6}$ first-order resonance mentioned in Sect. 2.4.. Using this approach with the nominal orbital elements given by Brasser et al. (2014) and Deitrick et al. (2018), we retrieve analytically their maps showing the possible obliquity variations of HD 40307 g and Kepler-62 f, in terms of the locations and widths of the secular spin-orbit resonances (see Fig. 8 by Brasser et al., 2014 and Figs. 5, 6, 10, 11 by Deitrick et al., 2018). The differences of oscillation amplitude that they observe are a natural consequence of the initial position of the planet with respect to the resonance centre (see Fig. 2). This shows that the Lagrange-Laplace system, associated with the development of the Hamiltonian (Sect. 2), is enough to obtain the level of detail required for studying the long-term rotation of exoplanets up to moderate eccentricities and inclinations. The use of a more elaborate model would add no substantial information, owing to the large uncertainties of the exoplanetary system under study.

3.2 Maximisation of $E_{k}$ and $S_{k}$

Unfortunately, several orbital elements remain unknown for most of the observed exoplanetary systems. The unknown elements usually include the mutual inclinations and the relative longitudes of ascending node. From now on, we suppose that only the masses, the semi-major axes and the eccentricities are known for all planets of the system. In this case, the Lagrange-Laplace matrices $A$ and $B$ can still be computed since they only depend on the masses and the semi-major axes. We thus obtain the two sets of frequencies $\mu_{j}$ and $\nu_{j}$. Because the initial conditions $z_{i}(0)$ and $\zeta_{i}(0)$ appearing in Eq. (37) are not fully known, the goal here is to obtain the maximum possible value of the amplitudes $E_{k}$ and $S_{k}$ according to the available data.

For the eccentricity, it amounts to maximise the modulus of a sum of complex numbers with unknown phase. The result is thus simply the sum of the moduli:

[TABLE]

using the fact that by definition, $|z_{i}(0)|=e_{i}$. The problem is more complex for the inclinations, since both the amplitudes and the phases of the initial conditions $\zeta_{i}(0)$ are unknown. It is thus necessary to introduce additional arguments, either from physical or from dynamical grounds. Guided by statistics on the orbital excitation due to close encounters, Atobe et al. (2004), while dealing mostly with systems with a single observed planet, imposed $I=e/2$ for each of them. This law is also in agreement with statistical distributions of observed exoplanetary systems (Xie et al., 2016). In our case, though, the application of this statistical result as a strict rule for each planet of a multi-planet system seems a bit simplistic. We will opt here for the hypothesis by Laskar & Petit (2017) of equipartition of the Angular Momentum Deficit (AMD) among the secular degrees of freedom. As they point out, this hypothesis is motivated both by theoretical arguments on chaotic diffusion in the secular dynamics (Laskar, 1994, 2008) and by the aforementioned correlations in observed distributions. As shown below, this allows to smooth the statistical law over all the planets contained in the system. Let us introduce the “coplanar AMD” of a planetary system, that is, the AMD it would have if it was strictly coplanar:

[TABLE]

where

[TABLE]

Contrary to Laskar & Petit (2017), we use here the Hamiltonian decomposition of Laskar & Robutel (1995), where the integrable part is the Sun-planet two-body problem. This is also the one chosen when expressing the matrices $A$ and $B$ of the Lagrange-Laplace system (Eqs. 34-35). The AMD equipartition hypothesis amounts to considering that the total AMD of the system,

[TABLE]

is equal to

[TABLE]

Hence, even if the individual orbital inclinations are not known, the so-obtained value of $C$ gives a bound for them. For instance, the maximum possible value of the inclination of the $k$th planet is given by

[TABLE]

In our case, we are trying to maximise the quantity

[TABLE]

obtained from (37) with unknown $\Omega_{i}$, using the constraint (42). This constraint can be rewritten

[TABLE]

where $\eta_{i}=\sin(I_{i}/2)$, $c_{i}=2\Lambda_{i}\sqrt{1-e_{i}^{2}}$ and $Z=C_{p}$, whereas the quantity to be maximised can be written

[TABLE]

where $b_{i}=\left|Q^{-1}_{ki}\right|$. The coefficients $c_{i}$ and $b_{i}$ are all positive, and $0\leqslant\eta_{i}\leqslant 1$. The constraint (45) forms an hyper-ellipsoid, whereas the quantity to be maximised (46) forms an hyperplane. Except from particular cases that we will dismiss here, there is thus only one solution for the maximisation of $Y$, which corresponds to the tangency of the plane and the ellipsoid. This implies that the two gradients are collinear:

[TABLE]

where $\lambda>0$ by definition of $\eta_{i}$. We get the value of $\lambda$ from the imposed value of $Z$:

[TABLE]

The maximum of $Y$ with the constraint $Z$ is thus

[TABLE]

Going back to the original notations, this finally gives

[TABLE]

However, Eq. (47) does not take into account the condition that all the $\eta_{i}$ are smaller than $1$. In practice, if the value obtained for $\lambda$ implies that one or several $\eta_{i}$ are larger than $1$, we just have to fix them to $1$ and use the same resolution method iteratively333From the form of the constraint (45), decreasing one $\eta_{i}$ to $1$ implies that at least one of the remaining $\eta_{i}$ should increase; from the solution (47), this actually means that all the remaining $\eta_{i}$ increase. with the remaining $\eta_{i}$ (changing the definition of $Z$ and $Y$ accordingly).

Table 2 shows the comparison of the amplitudes obtained for the Earth with the full Lagrange-Laplace system, and their maximisation supposing that the mutual inclinations and longitudes of node are unknown. As shown in Fig. 4, a chaotic map can be obtained using these maximum values. However, we note that each maximisation is specific to one single amplitude since it implies a distinct set of inclination values. Taking all the maximum amplitudes at once as if they formed one single representation gives thus a large upper bound for the chaotic zones. Moreover, due to the large amplitudes of the series obtained, the small-width approximation for second and third-order resonances does not necessarily hold.

For simple systems as those studied by Brasser et al. (2014), Deitrick et al. (2018) or Shan & Li (2018), the resonances are thin and well separated. A picture of the resonant regions (which do not overlap, in these cases) is thus enough to give clear view of the dynamics. Using the maximised amplitudes gives the largest possible widths of the resonances, which, in turns, show the maximum obliquity variations and their locations. We fully retrieve their results. In contrast to these simple and ordered dynamics, Fig. 5 shows as an illustration the maximised chaotic regions for the GJ 3293 system. The available orbital elements are taken from Astudillo-Defru et al. (2017), who pointed out that GJ 3293 d is in the habitable zone. In the spin-down process toward synchronous rotation due to tidal dissipative effects from the star (thus decreasing the precession constant), planet d is the most likely to suffer from large obliquity changes. One must remember, though, that the chaotic zones are here maximised according to the available orbital data. We also predict a rich obliquity dynamics for the Trappist-1 planets, but due to the confirmed strong effects of mean-motion resonances (see e.g. Quarles et al., 2017), the use of the Lagrange-Laplace model is probably inadequate in this case. Building an orbital theory specific to this system would be out the scope of this paper.

3.3 Maximisation of $\alpha$

In Sects. 3.1 and 3.2, we saw how to obtain a quasi-periodic approximation of the long-term orbital motion of a planet, and how to estimate bounds for its coefficients if all the orbital parameters are not known. However, in order to study its long-term spin dynamics, we still lack an estimate of its precession constant $\alpha$. Even in the Solar System, the precession constants of the planets are not very well known. The estimate of $\alpha$ for an extrasolar planet would require observations that are very hard to obtain (its rotation period and a model of interior), and which would be specific to one exoplanet. In order to keep the study as general as possible, we will not try here to obtain a single value for the precession constant: instead, we will look for an upper bound of it from general physical considerations.

After the sphere, the simplest shape model for a rotating planet is given by the Maclaurin ellipsoid (Chandrasekhar, 1969). It describes the equilibrium shape of a self-gravitating homogeneous body in rotation with constant angular velocity. The rotational symmetry is imposed (circular equator), leading to the formula

[TABLE]

where $\omega$ and $\rho$ are the rotation velocity and the density of the body and $\epsilon$ is the eccentricity of its ellipsoidal figure (in any plane containing the rotation axis). Studying $f=\omega^{2}/(2\pi\mathcal{G}\rho)$ as a function of the ellipsoid eccentricity, $f$ is zero for $\epsilon=0$ and $\epsilon=1$, and it has one maximum at $\epsilon_{0}\approx 0.929956$ with value $f_{0}\approx 0.224666$. This implies that there is no such equilibrium figure possible for rotation velocities larger than

[TABLE]

Converting the ellipsoid eccentricity in terms of momenta of inertia, we obtain the relation

[TABLE]

Injecting it into the expression of the precession constant (2), we obtain the maximum value

[TABLE]

For rotation velocities close to $\omega_{\text{max}}$, it is known that there exist equilibrium ellipsoidal figures with three unequal axes that have a lower total energy, called Jacobi ellipsoids (Chandrasekhar, 1969). However, we only need an order of magnitude for $\alpha_{\text{max}}$ and the homogeneous approximation is anyway quite crude, allowing us to stick to the expression (54). Planets are expected to spin much more slowly than $\omega_{\text{max}}$ (Eq. 52), including giant gaseous planets (Batygin, 2018). In the remaining part of the article, we allow us to abusively refer to the “rotational breakup” velocity.

Using the average density of the Earth, we obtain a minimum rotation period of about $2.4$ hours, leading to a maximum precession constant of about $230\leavevmode\nobreak\ ^{\prime\prime}$/yr. The true value for the Earth is $20\leavevmode\nobreak\ ^{\prime\prime}$/yr, or $50\leavevmode\nobreak\ ^{\prime\prime}$/yr if we include the additional effects of the Moon (Laskar & Robutel, 1993). This remains well below our bound, but the difference with this “effective” precession constant due to the presence of satellites actually constitutes the largest source of uncertainty. Close satellites increase the effective flattening of the planet, whereas far satellites increase the effective torque from the star (Boué & Laskar, 2006). In both cases, this increases the precession constant to be used in our model. This effect is particularly problematic for Saturn, because our upper bound gives $\alpha_{\text{max}}\approx 0.75\leavevmode\nobreak\ ^{\prime\prime}$/yr, while the true value is $0.20\leavevmode\nobreak\ ^{\prime\prime}$/yr, but it increases to $0.83\leavevmode\nobreak\ ^{\prime\prime}$/yr if we take the satellites into account (Ward & Hamilton, 2004). Hence, the introduction of Saturn’s satellites makes the precession constant exceed our upper bound. This problem is unavoidable for exoplanets, because the observation of satellite systems is hard and none has been observed so far. When using our model to study the spin dynamics of exoplanets, we must thus always keep in mind that the presence of numerous or massive satellites could modify our conclusions for borderline cases like Saturn.

Moreover, we must assume a density $\rho$ for the exoplanets if their radius has not been measured, which adds even more uncertainty. If the radius is unknown, an order of magnitude of the density can be estimated through an empirical law adjusted to the observed mass-radius distribution (see Fig. 6). The density is anyway not the major source of uncertainty of our method.

3.4 Classification of non-resonant exoplanets

In Sect. 2.3, we saw that the overlap of first-order secular resonances, leading to the largest chaotic zones of the spin dynamics, can be produced only if the ratio $\nu_{i}/\alpha$ is smaller than $1$ at least for two frequencies $\nu_{i}$ of the inclination quasi-periodic representation. Assuming that $\alpha$ is bounded by $\alpha_{\text{max}}$, we can deduce that there can be no substantial chaotic zone if $\nu_{i}/\alpha_{\text{max}}>1$ whatever the frequency $\nu_{i}$. Moreover, if the mutual inclinations are small (as we assume they are), this implies that the obliquity is almost constant. Our bound for $\alpha$ (Eq. 54) can thus be used for a preliminary classification of the exoplanets, while the bounds for the amplitudes (Eqs. 38,50) are required for a more specific application to one exoplanet (they allow to constrain both $\gamma$ and $\beta$, see Fig. 1).

Table 3 shows the 94 planets classified strictly non-resonant with this criterion, using the Lagrange-Laplace matrix to estimate the frequencies (Sect. 3.1) and Eq. (54) as a bound of $\alpha$. All the exoplanets from http://exoplanet.eu with known mass, semi-major axis and eccentricity were analysed (taking $m\sin I$ instead of the mass if a real-mass estimate was unavailable). At date 2018-03-07, this represents 143 systems with more than one planet, which contain 353 planets in total (plus the Solar System). For some of them, the frequency ratios are so far from $1$ that their classification is quite safe, even when considering the numerous sources of error inherent to our method, and in particular, the possible presence of satellites. This mostly concerns planets that are far from their star, like Uranus and Neptune, and no terrestrial exoplanet has been observed yet is in this category. Such large semi-major axes (third column of Table 3) imply that most of the planets listed in Table 3 are also unaffected by orbital and rotational tidal dissipation resulting from the interaction with their central star.

Most of the exoplanetary systems known so far contain only two planets. In this case, there can be no chaotic region anyway coming from the overlap of first-order resonances (Sect. 2.3) because there is only one forced frequency in inclination (Sect. 3.1). However, this single term allows to go one step further and compute the variation range of the obliquity at first order. This is obtained by looking at the interval of parameters $\gamma$ and $\beta$ (Fig. 1) allowed for the exoplanet. A very simple formula can be derived if the inclination amplitude $S_{1}$ (which is the only amplitude of the decomposition) is small. Indeed, this implies that $\beta$ is small as well (Eq. 26), resulting in quite flat level curves for Colombo’s top Hamiltonian (18). The maximum obliquity variations are achieved around $\Sigma=0$, and at leading order, they are equal to

[TABLE]

This approximation holds very well for small values of $S_{1}$. For the exoplanet WASP-81 c, which has the lowest bound for $S_{1}$ in Table 3, we obtain a maximum obliquity excursion of $0.5^{\text{o}}$. This limit is very close to what is obtained by plotting the level curves of the Hamiltonian (18), and it holds as long as $\alpha$ and $S_{1}$ are below their estimated bounds. Such a good constraint cannot be achieved for every planet in two-planet systems, though, since our bound on $S_{1}$ from the AMD is sometimes not very informative. It is even dramatic for planets perturbed by a very massive companion: for example the maximum inclination amplitude of HD 92788 c is higher than $1$, indicating that all inclinations are possible.

The maximum excursion of the obliquity is harder to obtain if there are more than two planets in the system, since the dynamics is ruled by the superposition of several forcing terms. However, if the maximum maximised amplitude is small (last column of Table 3), the superposition of all the terms is unlikely to bring the obliquity over the limit given at Eq. (55). It can thus be used as well as an order-of-magnitude estimate. This results in a maximum of about $20^{\text{o}}$ for Uranus and Neptune, showing the crude nature of our maximisation (as shown by previous works, it is very hard to tilt Uranus and Neptune by the mean of planetary perturbations, see e.g. Boué & Laskar, 2010).

4 Conclusion

The spin-axis dynamics of a planet plays a major role in its climate setting, and, by extension, in its suitability for life. However, the rotation properties of exoplanets are still very poorly constrained. In this paper, we presented an analytical formulation of the long-term spin-axis dynamics of a planet, allowing to link known and unknown parameters to its obliquity evolution and to provide a global picture of the dynamics in a straightforward way.

At first, the orbital solution is modelled by quasi-periodic series. This method is thus valid as long as the orbital chaos, if any, takes place on a much larger timescale than the spin-axis evolution. The spin-axis Hamiltonian is then expanded in powers of the eccentricity and inclination amplitudes of the orbital series. We provided all terms up to order 3 but the development can be conducted to higher orders.

A clear picture of the phase space structure is given by the obliquity ranges associated to the various resonant regions. The resonant dynamics at order 1 can be characterised analytically in terms of two parameters, which are linked to the precession constant $\alpha$ (gathering the physical characteristics of the planet under study) and to the quasi-periodic representation of the orbit. The pendulum approximation is only used at order 2 and beyond, for which the resonances are thin enough. The regions of resonance overlap at all orders are identified as chaotic. In some cases (as for the terrestrial planets of the Solar System), these chaotic regions allow wide excursions of the obliquity. The method presented here allows to retrieve analytically the previous numerical results with a good precision. Numerical integrations prove thus to be necessary only if detailed statistics on the obliquity evolution are required. This is very informative for Solar System planets (as shown by Néron de Surgy & Laskar, 1997; Correia & Laskar, 2003; Laskar et al., 2004a) but not yet for exoplanets because their initial conditions and physical parameters are still poorly known. Hence, the uncertainty of our results remains largely dominated by our lack of knowledge of the exoplanetary systems rather than by the approximations inherent to our method. This allows to stick to the simple analytical formulas presented here.

At this level of uncertainty, the Lagrange-Laplace system provides a good-enough representation of the orbital motions (excepted for exoplanetary systems featuring highly excited orbits or strong effects of mean-motion resonances). The formulas obtained allow to set an upper bound for the amplitude of the eccentricity and inclination terms if the mutual orientations of the orbits are unknown. On the other hand, the AMD equipartition hypothesis can be used, if required, to place a bound on the inclination from the eccentricity values. Through our analytical model of the spin-axis dynamics, these maximum amplitudes provide the maximum extent of the chaotic zones. For example, a large chaotic region is expected for exoplanet GJ 3293 d for rotational velocities above the synchronous rotation. Systems very affected by mean-motion resonances (like Trappist-1) can still be studied using the method described here, but with the prior construction of a synthetic representation for the orbital motion, written in the form of a quasi-periodic series.

However, this method does not allow to consider tidal dissipations (playing an important role for exoplanets close to their star), which could be modelled as an adiabatic process acting on a much longer timescale than the obliquity variations (see Néron de Surgy & Laskar, 1997). This amounts to make the precession constant $\alpha$ and/or the amplitudes of the orbital series gradually vary. This method does not include either the effects of libration around spin-orbit resonances, even if a trick allows to take into account a possible locking in synchronous rotation (Appendix A).

Finally, under the hypothesis of hydrostatic equilibrium, we can set a bound to the precession constant $\alpha$. This bound is obtained from the flattening of the planet corresponding to its rotational breakup velocity. Since $\alpha$ governs the width and location of the resonances, this allows to classify the exoplanets that cannot be subject to first-order secular spin-orbit resonances. Among the sufficiently known systems with more than one planet, we found $94$ planets in this category ($26\%$ of our sample). If they belong to exoplanetary systems with low mutual inclinations (as it is expected in most cases for orbital stability), this implies that their obliquity is almost constant. This bound for $\alpha$ is though invalidated by the possible presence of massive satellites (as our Moon), but some exoplanets are so far from resonance that their classification is quite safe. This is the case of Uranus and Neptune.

Considering the high efficiency of the analytical method proposed here, an obliquity stability map could be designed easily in the future for each new exoplanet discovered, and in particular for those classified as “habitable”. However, such a stability map should always be computed again if any additional planet is found in the system. Indeed, it would shift the existing frequencies (especially if the new planet is massive), and add one frequency in both the inclination and eccentricity series, multiplying the possibilities of resonance. On the other hand, the total AMD of the system would increase, resulting in wider maximised chaotic zones.

Acknowledgements.

We thank the anonymous referee for her or his detailed review.

Appendix A Case of a $1\!:\!1$ spin-orbit resonance

Numerous exoplanets are observed very close to their star, in a place where the tidal frictions are strong enough to efficiently lock them in synchronous rotation. In this section, we show that if the librations around the synchronous rotation are much faster than the secular spin-axis dynamics, we can retrieve Colombo’s top Hamiltonian (Sect. 2.2), allowing to use the same approach as in the non-resonant case. As before, though, we will not consider the effect of the tidal dissipation on the obliquity. This is thus only valid for systems for which the tidal damping of the obliquity acts on a larger timescale than the spin-axis dynamics.

We will use the same method as Correia et al. (2003). Let us write $\lambda$ the mean longitude of the planet in orbit around the star, and $\ell$ its rotation angle. The mean longitude $\lambda$ is measured from the equinox at a reference epoch (for instance J2000), whereas the rotation angle $\ell$ is measured from the equinox of the date up to a fixed point of the equator (principal axis A). If we keep the angles of the form $\ell-\lambda$ during the average over the mean longitude and the fast rotation angles (see Néron de Surgy & Laskar 1997), the corresponding “semi-averaged” Hamiltonian is

[TABLE]

where we neglected terms of order $e(B-A)/C$. The momenta $L=C\omega$ and $Y=LX$ are conjugate to $\ell$ and $-\psi$, respectively. The momentum $\Lambda$, conjugate to $\lambda$, has been added such that $\dot{\lambda}=n$ (mean motion). The resonant precession constant is defined as

[TABLE]

using the same notation as Eq. (2). We note that the angle $\lambda+\psi$ appearing in the Hamiltonian corresponds to the mean longitude measured from the equinox of the date. Let us use the canonical change of coordinates

[TABLE]

The momentum $\Gamma$ is an arbitrary constant of motion and the Hamiltonian becomes

[TABLE]

We will now suppose that the dynamics of $\theta$, corresponding to the “semi-secular” timescale (either circulation or oscillation), is much faster than the evolution of the other degrees of freedom, corresponding to the secular timescale. We thus consider for now that except $(I,\theta)$, all the variables are fixed (adiabatic approximation). The equations of motion are

[TABLE]

Using the definition of $I$, $Y$ and $\alpha_{r}$, the first equation gives

[TABLE]

resulting, for any value of $X$, to two equilibrium points: $\theta=\psi$ and $\psi+\pi/2\mod\pi$. We note that $\theta=\psi$ is an elliptic equilibrium while $\theta=\psi+\pi/2$ is hyperbolic. Injecting this into the second equation, we obtain

[TABLE]

in which the small terms correspond to the precession of the spin axis ($\alpha$ and $\alpha_{r}$) and the precession of the orbit ($\mathcal{A}$ and $\mathcal{B}$). The equilibrium condition, corresponding to the exact resonance, is thus $\omega\approx n$. Considering that the planet is locked in synchronous rotation, we have thus $\theta=\psi$ and $\omega\approx n$. According to the adiabatic approximation, this will be verified whatever the value of the slow variables, such that we can inject them into the full Hamiltonian:

[TABLE]

where this time, we use $X$ as conjugate momentum of $-\psi$ (the Hamiltonian is thus divided by the constant $L$). In the expression of $\alpha$ and $\alpha_{r}$, we must replace $\omega$ by $n$. We get here one extra term with respect to (1), due to the spin-orbit resonance. Using the same method as in Sect. 2.2, the Hamiltonian in case of a first-order secular spin-orbit resonance is

[TABLE]

which must be compared to (16). This Hamiltonian has the same general form and it can be reduced to Colombo’s top. We can thus apply the same method of resolution (redefining the constants accordingly).

Appendix B Characteristic quantities of Colombo’s top

B.1 Equilibrium points

From (18), the equations of motion are

[TABLE]

Apart from the coordinate singularity at $\Sigma=\pm 1$, the first equation implies that $\dot{\Sigma}=0$ when $\sigma=0$ or $\pi$. Injecting this into the second equation, we get

[TABLE]

where $\beta\geqslant 0$ by hypothesis. The resolution of this equation requires to square left and right-hand terms, loosing the information666After having computed one solution $\Sigma_{0}$, this information is retrieved by checking the sign of $\Sigma_{0}/(\gamma-\Sigma_{0})$. about the sign of $\cos\sigma$. We obtain a quartic equation in $\Sigma$:

[TABLE]

with discriminant

[TABLE]

It is zero for the particular cases $\gamma=0$ or $\beta=0$, for which the polynomial can be factored into, respectively,

[TABLE]

showing the corresponding solutions and their multiplicities. They constitute equilibrium points of the system whenever they are real and in the interval $[-1;1]$.

For $\gamma>0$ and $\beta>0$, the discriminant can be either negative (two equilibrium points), zero (three equilibrium points among which one double root), or positive (four equilibrium points). The corresponding solutions can be written analytically according to the general resolution of quartic equations. They are namely

[TABLE]

where numerous intermediary variables are required in order to get compact expressions:

[TABLE]

We note that $\Sigma_{c,d}$ are real solutions only when $\Delta_{4}\geqslant 0$ (see below for the limit in terms of $\gamma$ and $\beta$). The corresponding values of $\sigma$ are

[TABLE]

The points $a$, $b$ and $d$ are elliptic fixed points, whereas the point $c$ is hyperbolic.

B.2 First boundary (BC/D)

The zero value of (68) corresponds to a bifurcation. Its position can be computed by solving the equation $\Delta_{4}=0$, which corresponds to solving a cubic equation either in $\gamma^{2}$ or $\beta^{2}$. Choosing to solve it in terms of $\beta$, the discriminant is

[TABLE]

meaning that there is only one real solution. This solution is

[TABLE]

which is the boundary $\mathscr{C}_{1}$ (20).

B.3 Second boundary (A/B)

The other two boundaries can be obtained by studying the level curves of the Hamiltonian passing through $\Sigma=\pm 1$ (which is singular using the coordinates $\Sigma$ and $\sigma$, but it does not matter here).

Let us begin with the $+1$ case, for which the Hamiltonian has value $-1/2+\gamma$. We now look for this specific level curve along the axes $\sigma=0$ and $\sigma=\pi$. This leads to the equation

[TABLE]

for which $\Sigma=+1$ is a solution. By reorganising the terms, taking the square (thus loosing the information about the sign of $\cos\sigma$), and dividing by $(\Sigma-1)$, we get

[TABLE]

which is a cubic equation in $\Sigma$. Its determinant is

[TABLE]

Once again, it is zero for $\beta=0$. Moreover the solutions for $\gamma=0$ can be easily computed. In these two particular cases, the polynomial can be factored into, respectively,

[TABLE]

showing the solutions and their multiplicities. For $\gamma>0$ and $\beta>0$, the discriminant can be either negative (one solution), zero (three solutions among which one double root), or positive (three solutions). The zero value corresponds to the limit we are looking for. Its position can be computed by solving the equation $\Delta_{3}=0$, which amounts to solving a quadratic equation in $\beta^{2}$ or a quartic equation in $\gamma$. Choosing to solve it in terms of $\beta$, the only positive solution is

[TABLE]

which is the boundary $\mathscr{C}_{2}$ (25).

B.4 Third boundary (B/C)

Let us now study the level curve of the Hamiltonian passing in $\Sigma=-1$, which has value $-1/2-\gamma$. The procedure is the same as for the second boundary, and the new formulas are obtained simply by replacing $\gamma$ by $-\gamma$. There is though an ambiguity because there are two positive solutions $\beta^{2}$ (as a function of $\gamma$) which cancel the determinant. The one corresponding to the bifurcation is the largest, that is,

[TABLE]

which is the boundary $\mathscr{C}_{3}$ (25).

B.5 Separatrices

The position at $\sigma=0$ or $\pi$ of the separatrix emerging from the hyperbolic point $(\Sigma,\sigma)=(\Sigma_{c},\pi)$ defines the boundaries of the resonant region (see Fig. 1). Writing $f=\mathcal{F}(\Sigma_{c},\pi)$, the equations to solve are

[TABLE]

The resolution of this equation requires to square left and right-hand terms, loosing the information777After having computed one solution $\Sigma_{0}$, this information is retrieved by checking the sign of $-\Sigma_{0}^{2}/2+\gamma\Sigma_{0}-f$. about the sign of $\cos\sigma$. We obtain a quartic equation in $\Sigma$,

[TABLE]

in which $\Sigma_{c}$ is a double root. It can thus be divided by $(\Sigma-\Sigma_{c})^{2}$, leading to the quadratic equation

[TABLE]

This equation has always two real solutions, provided that $\Sigma_{c}$ exists (that is, in zones A, B or C). These solutions are

[TABLE]

where we replaced $f$ by its expression (18) in terms of $\Sigma_{c}$.

Appendix C Second-order resonances

Using the intermediary Hamiltonian $\mathcal{X}=\varepsilon\mathcal{X}_{1}$ (31), the Hamiltonian in the new coordinates is obtained term by term from Eq. (28). The two first terms are simple: we have $\tilde{\mathcal{H}}_{0}=\mathcal{H}_{0}$ (given at Eq. 11) and $\tilde{\mathcal{H}}_{1}=0$ by definition of $\mathcal{X}$. The second-order term is more complex since it requires to compute Poisson’s brackets. Using of the fact that $\{\mathcal{X}_{1},\mathcal{H}_{0}\}=-\mathcal{H}_{1}$ and reorganising the terms adequately, we obtain

[TABLE]

Since by hypothesis there is no first-order resonance in the system, the only possible resonant angles at second order are of the form $\sigma=\phi_{j}+\phi_{k}+2\psi$. Let us perform the canonical change of coordinates

[TABLE]

and

[TABLE]

Assuming that $\sigma$ is the only resonant angle, the dynamics at second order is given by averaging $\tilde{\mathcal{H}}$ over all other angles (this is another change of coordinates close to identity). The momenta $\Gamma_{1}$ and $\Gamma_{2}$ become arbitrary constants of motion that we will conveniently choose equal to zero. Dropping the unnecessary constants, the resonant Hamiltonian is thus

[TABLE]

in which $X$ must be replaced by $-2\Sigma$. High-order resonances are quite thin, so it is enough to consider the dynamics in the neighbourhood of the resonance centre at first order:

[TABLE]

or equivalently $X_{0}=-2\Sigma_{0}=-(\nu_{j}+\nu_{k})/(2\alpha)$. Considering that $|X-X_{0}|=\mathcal{O}(\varepsilon)$, we have then

[TABLE]

in which we dropped the unnecessary constants, and where

[TABLE]

By injecting the momentum $\Sigma$ instead of $X$ and by using the modified time $\mathrm{d}\tau=-4\alpha\mathrm{d}t$, we obtain

[TABLE]

This is the Hamiltonian of a pendulum of centre $\Sigma_{0}$ and half width $\sqrt{|K|}$. In terms of the obliquity cosine $X$, the resonance has position $X_{0}$ and half width $2\sqrt{|K|}$.

Appendix D Third-order resonances

If there is no resonance at first and second orders, we can use a canonical change of coordinates close to identity in order to suppress the angular dependency at first and second orders. Let us consider an intermediary Hamiltonian $\mathcal{X}=\varepsilon\mathcal{X}_{1}+\varepsilon^{2}\mathcal{X}_{2}$, such that the new coordinates are given by its flow at time $1$. The Hamiltonian in the new coordinates is then

[TABLE]

where

[TABLE]

The first-order part of $\mathcal{X}$ required to suppress the angular dependency at order 1 can be directly taken from Eq. (31). Let us write

[TABLE]

(average plus oscillating part) for which the expression is given by (85). The homological equation for order 2 is then

[TABLE]

which defines the Hamiltonian $\mathcal{X}_{2}$. This leads to

[TABLE]

We must now compute the remainders at order 3. First of all, we can simplify their expressions by taking into account that, by definition: $\{\mathcal{X}_{1},\mathcal{H}_{0}\}=-\mathcal{H}_{1}$, $\{\mathcal{X}_{2},\mathcal{H}_{0}\}=-\widetilde{\mathcal{A}_{2}}$ and $\{\mathcal{X}_{1},\mathcal{H}_{1}\}=2(\mathcal{A}_{2}-\mathcal{H}_{2})$. We have then

[TABLE]

which gives

[TABLE]

where

[TABLE]

[TABLE]

and:

[TABLE]

The expression of the coefficients $D_{j}$ is very complex. We will not give it here since they have no interest at this stage (the angles $\phi_{j}+\psi$ are non resonant by hypothesis).

As shown in Appendix C, in the pendulum approximation, the half-width of any possible resonance is two times the square root of its coefficient divided by $\alpha$, and its position is given by the combination of the unperturbed frequencies. Accordingly, the possible resonances at order 3 are gathered in Table 1.

Appendix E Geršgorin circles

In order to prove that the Lagrange-Laplace matrix for the orbital inclinations has only negative or zero eigenvalues, one can use the Geršgorin circle theorem (see Geršgorin 1931 or Varga 2004). This theorem is recalled below, and we show how it applies to our matrix.

Definition.

Let $B$ be a complex $N\times N$ matrix with elements $(b_{ij})$. The $i$th “Geršgorin disc” $\mathscr{G}_{i}$ $(i=1,2..N)$ is the closed disc of the complex plane centred at $b_{ii}$ and with radius

[TABLE]

Theorem (Geršgorin 1931).

Any eigenvalue of $B$ lies inside at least one of the $\mathscr{G}_{i}$ discs, $i=1,2..N$.

Corollary.

All the eigenvalues of $B$ are located inside the union of the $\mathscr{G}_{i}$ discs, $i=1,2..N$.

In our case, the matrix $B$ is real (see Eq. 35). It has only real eigenvalues and one of them is identically equal to zero. Moreover, given the very particular form of this matrix, the centre of each Geršgorin disc is located on the real line, with an abscissa equal to the opposite of its radius. Therefore, all the eigenvalues of $B$ are negative or zero, as illustrated in Fig. 7.

Appendix F Orbital solution used for the inner Solar System

In order to apply our method to a given planet, we first need a quasi-periodic approximation of its long-term orbital dynamics.

In the case of the Solar System, the search for such series has been a challenge for centuries, eventually leading to very complete solutions (up to the degree of chaos inherent to the system). In the present work, we use the solution of Laskar (1990), obtained by multiplying the normalised proper modes $z_{i}^{\bullet}$ and $\zeta_{i}^{\bullet}$ (Tables VI and VII of Laskar 1990) by the matrix $\tilde{S}$ corresponding to the linear part of the solution (Table V of Laskar 1990). In the series obtained, the terms with the same combination of frequencies are then merged together, finally resulting in 56 terms in eccentricity and 60 terms in inclination.

These series are given in Tables 4-7 for the inner planets, under the form:

[TABLE]

with $N=56$ and $M=60$. They are used in Fig. 3 of the present work.