Issue 
A&A
Volume 623, March 2019



Article Number  A4  
Number of page(s)  21  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834344  
Published online  26 February 2019 
Secular spinaxis dynamics of exoplanets
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, LAL, Université de Lille,
75014
Paris,
France
email: melaine.saillenfest@obspm.fr
Received:
28
September
2018
Accepted:
7
January
2019
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 landbased life, but orbital and rotational parameters of exoplanets are still poorly constrained. Numerical explorations usually realised in this situation are therefore in heavy contrast with the uncertain nature of the available data.
Aims. We aim to provide an analytical formulation of the longterm spinaxis 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 enable a quick and straightforward exploration of the spinaxis dynamics.
Methods. The longterm orbital solution is decomposed into quasiperiodic series and the spinaxis 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 the angular momentum deficit).
Results. This method gives accurate results when the orbital evolution is well known. The detailed structure of the chaotic zones for the solar system planets can be retrieved from simple analytical formulas. For lessconstrained planetary systems, the maximal extent of the chaotic regions can be computed, requiring only the mass, the semimajor 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 spinorbit secular resonances (unless the precession rate is affected by the presence of massive satellites).
Key words: planets and satellites: dynamical evolution and stability / celestial mechanics / planets and satellites: general
© M. Saillenfest et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
From the works by Laskar & Robutel (1993) and Laskar et al. (1993a), we know that the longterm 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^{o} to 60^{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 spinaxis evolution of the Earth (Néron de Surgy & Laskar 1997; Laskar et al. 2004a; Li & Batygin 2014a), Venus (Correia et al. 2003; Correia & Laskar 2003), and Mars (Laskar et al. 2004b). On the other hand, paleorecords for 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 highobliquity planets undergo severe seasonal variations (e.g. Spiegel et al. 2009) even for a stable obliquity. Nevertheless, 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 therefore of more than simple academic interest. This motivated further studies of the structure of the chaotic region (Li & Batygin 2014b) and some additional numerical experiments (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 attainable 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 spinaxis dynamics). Several authors have already tackled this problem (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 timedependent perturbation of the spin axis. Applied to extrasolar planets, this method requires a nominal value to be chosen for both the rotational parameters and the unknown orbital elements. The parameter space to be explored is thus very wide, meaning that even elaborate numerical explorations require some degree of arbitrariness. Consequently, the use of numerical integrations at this stage could appear to be 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 suggestions for a possible analytical treatment were actually given by Laskar (1996) and partially exploited by Atobe et al. (2004), Li & Batygin (2014b), 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 sensitivity of the spinaxis 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 hypotheticalterrestrial planets in the habitable zone of known exoplanetary systems. Their analytical developments, though, were limited to the lowestorder 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 Kepler62f and Kepler186f. This time, their calculations were mostly designed to support and explain numerical results, and they did not try to substantially improve the theory by Atobe et al. (2004).
In this context, the goal of this article is twofold: (i) we aim to provide a general analytical formalism for studying the longterm spinaxis dynamics of (exo)planets and (ii) to 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 tothe characteristics of their spin dynamics. This would allow for the determination of the exoplanets that merit further, detailed study (in particular if a complex chaotic spin dynamics is expected) and those that have a necessarily very simple spin dynamics, thus negating the need for 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 spinaxis dynamics and shows how it can be expanded in terms of the orbital motion parameters. The secular resonances at all orderscan then be isolated and used to delimit the chaotic regions. Subsequently, 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, this method allows for a preliminary classification of the “nonresonant” exoplanets, for which no chaos can appear and for which the obliquity variations are constrained by an analytical bound.
2 Analytical model of the longterm spin dynamics
2.1 Development of the Hamiltonian
We 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 spinorbit resonance. Considering only the lowestorder 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). This Hamiltonian can be written (1)
where the conjugate coordinates are X (cosine of obliquity) and − ψ (minus the precession angle). The quantity α is referred to as the precession constant (contrary to previous studies, we prefer to exclude here the eccentricity e appearing in the denominator from the definition of α). Following the derivation proposed by Néron de Surgy & Laskar (1997), we obtain (2)
where is the gravitational constant, m_{0} is the mass of the star, a is the semimajor axis of the planet in orbit around the star, ω is its spin angular velocity, and A ≤ B ≤ C are its momenta of inertia. The Hamiltonian (Eq. (1)) depends explicitly on time t through the eccentricity e and the functions (3)
in which q = sin(I∕2)cos Ω and p = sin(I∕2)sin Ω, where I and Ω are the orbital inclination and the longitude of ascending node of the planet, respectively. In the following, we write η ≡ sin(I∕2). One should note that if there is only one planet in the system (twobody problem), the obliquity is constant and the precession angle circulates with constant angular velocity .
We consider a scenario where the eccentricity and the inclination of the planet are small, such that we are allowed to develop the Hamiltonian in series of e and η. In the following, we present the terms up to order 3, but the method presented here can be generalised to any order (as we see below, the third order is the first one at which the eccentricity begins to play a substantial role). Using the fact that , we obtain (4)
We nowconsider the secular orbital dynamics of the planet, resulting from the perturbations by the other planets, to be quasiperiodic. 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 see below, this holds very well for the solar system (this methodology was first proposed by Laskar 1996; it is used e.g. by Li & Batygin 2014b). In this case, we can write (5)
where ϖ is the longitude of pericentre of the planet in orbit around the star. The angles θ_{j} and ϕ_{j} evolve linearly with frequencies μ_{j} and ν_{j}, that is, (6)
whereas the amplitudes E_{j} and S_{j} are real constants of order e and η 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, μ_{j} and ν_{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 thus obtain (7) (8)
and from Eq. (4), (9)
In order to obtain an autonomous Hamiltonian, we introduce the momenta Θ_{j} and Φ_{j} conjugate to θ_{j} and ϕ_{j}. The system has now N + M + 1 degrees of freedom. The new Hamiltonian(that we still denote ) can be written (10)
where we suppose that . The different parts are, respectively, (11) (12) (13)
We note that there can be no resonance among the angles ϕ_{j} and θ_{j} because they come from the quasiperiodic solution of the orbital dynamics. By definition, they are thus already “integrated”.
2.2 One perturbing term: Colombo’s top
From Eqs. (11)–(12), we conclude that at lowestorder to the perturbation, resonant angles can only be of the form σ = ϕ_{j} + ψ. We now consider a single resonance with the term j = p. We introduce the resonant canonical coordinates by the linear transformation (15)
Assuming that the system is far from any other resonance, the longterm dynamics at first order to the perturbation is given by averaging the Hamiltonian over all angles but σ. Dropping the constant terms, we get (16)
As shown in Appendix A, the Hamiltonian has the same form in the case of a 1:1 spinorbit resonance, with a slightly different expression of the coefficients however. This would only slightly shift the position of the secular resonances considered here. In this case, the tidal damping (responsible for this capture in spinorbit 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 dτ = α dt to Eq. (16), we obtain the following Hamiltonian (that we still denote ): (18)
The dynamical system with Hamiltonian function (Eq. (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 γ and β.
First of all, it is enough to study the case γ ≥ 0 and β ≥ 0 since we get the negative cases by the transformations Σ →−Σ and σ → σ + π, respectively. When studying the equilibrium points of the system (Henrard & Murigande 1987, or Appendices B.1 and B.2), we find that the phase space can have two different geometries according to the value of γ and β (see Fig. 2 for the phase portraits). The boundary between these two regions of the parameter space is the curve (20)
Below (regions A, B, C of Fig. 1), the dynamical system (Eq. (18)) has four equilibrium points, which we denote (a, b, c, d). The equilibrium points c and d merge along the curve , and disappear above it (region D of Fig. 1). Their respective positions are (21)
The values Σ_{a,b,c,d} have explicit expressions in terms of γ and β, as given in Appendix B.1 (they correspond to the different roots of a quartic equation). Since the resonant angle ψ + ϕ_{p} is equal to 0 or π for each of these equilibrium points, they all correspond to configurations where the spin axis, the normal to the orbit (reduced to its pth harmonic), and the normal to the reference plane are all 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 (22)
We are now interested in the width of the resonance, that is, the interval of Σ 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 β (as used by Atobe et al. 2004 or Li & Batygin 2014b), but this approximation is no longer valid when β grows. Since an analytical expression can be derived even in the general case, we use it here. The computations (Appendix B.5) lead to the following extreme values of Σ spanned by the resonance: (24)
They are defined whenever Σ_{c} itself is defined (regions A, B, C of Fig. 1). At this point, it is important to note that the coordinates (Σ, σ) are singular at Σ = ±1 since the problem actually takes place on the sphere (Henrard & Murigande 1987). We must therefore carefully study the meaning of the limits (Eq. (24)) when one of them crosses ± 1. This leads to two other limits in the parameter space, as the curves (25)
(see Appendices B.3 and B.4), delimiting the regions A–B and B–C of Fig. 1, respectively. We note that and intersect at (γ;β) = (1 ;0), and intersect at , and and intersect at (0 ; 1∕2). Contrary to , the boundaries and 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 [Σ_{−}; Σ_{+}] in region A; in [Σ_{−};+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 Σ = 1.
When β → 0, the width of the island tends to 0 and all the level curves in the (σ, Σ) plane tend to be horizontal. In this limit, the resonance width is small and almost independent of γ (pendulum approximation). When γ → 0, the phase portrait in the (σ, Σ) plane tends to be symmetric with respect to the line Σ = 0.
From Eq. (19), we note that (26)
Remembering that S_{p} is the amplitude of a given term in the inclination series (Eq. (5)), this means that the parameter region above the line β =2γ in Fig. 1 (that is, S_{p} = 1) cannot be reached in this problem.
As a simple rule of thumb, one can consider that β controls the resonance width, and that γ controls the location of the resonance centre (see the Hamiltonian in Eq. (18)). From Eq. (26), we know that β 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 byBoué & Laskar (2010) in their scenario for tilting the spinaxis 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 β depends on the precession constant α (Eq. (2)) as well, meaning that “small” mutual inclinations do not guarantee that the resonances are thin. For the terrestrial planets of the solar system, some values of α produce firstorder resonances larger than 70^{o}, even if the mutual orbital inclinations are modest (see Sect. 2.3).
Contrary to β, the magnitude of γ cannot be easily traced back from the planetary architecture: the firstorder secular spinorbit resonances of any planet can be located anywhere between 0^{o} and 180^{o} of obliquity. The large majority of them actually lie in [0^{o};90^{o}] because most of the frequencies ν_{p} are negative(the explanation for this property is given in Sect. 3.1).
Fig. 1 Parameter space of Colombo’s top Hamiltonian (Eq. (18)). There is no separatrix (and thus no resonance) in region D, delimited by the curve (Eq. (20), thick black line). The curves and (Eq. (25), thin black lines) delimit regions A–B and B–C, respectively. In region A, the resonance island lies between Σ_{−} and Σ_{+} (Eq. (24)); in region B, it lies between Σ_{−} and + 1; in region C, it lies between − 1 and + 1 (hence the 180^{o} full width given by the colour scale). The problem becomes unphysical above β = 2γ (blue line) since it corresponds to an amplitude S_{p} > 1 in the inclination series. 
Fig. 2 Examples of phase portraits for the four regions of Fig. 1. The equilibrium points are labelled as in Eq. (21). The level curves of the Hamiltonian are drawn with black lines out of the resonance island and with red lines inside the resonance. The parameters chosen are, from A to D: (γ, β) = (0.4, 0.1); (0.55, 0.15); (0.05, 0.7); and (0.8, 0.3). 
2.3 Overlap of firstorder resonances
Going back to the full Hamiltonian (Eq. (10)), the main chaotic regions of the system can be estimated as the overlap of the firstorder 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 quasiperiodic representation (Eq. (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 quasiperiodic 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 β coefficients, resulting in virtually no secular variation of the obliquity (see the shades of grey in Fig. 1). A nochaos criterion can be obtained even if the amplitudes S_{i} are not known: we simply have to verify that γ_{i} ≈ ν_{i}∕α > 1 ∀ i = 1, …, M. On the contrary, if γ_{i} < 1 for at least two i, the existence of a chaotic region is possible, but not guaranteed. These properties are 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 quasiperiodic 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 firstorder resonances for the terrestrial planets (light and darkred 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 firstorder resonances overlap (there are almost no lightred 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 firstorder resonances. The following section is therefore dedicated to second and thirdorder resonances.
Fig. 3 Top row: estimate of the chaotic regions of the spin dynamics as the superposition of secular spinorbit resonances. The orbitalevolution of each planet is approximated by the synthetic representation of Laskar (1990), as detailed in Appendix F. Lightred and darkred regions represent the firstorder resonances and their overlaps, respectively(Colombo’s top Hamiltonian, Sect. 2.2); lightblue and darkblue regions represent the secondorderresonances and their overlap (Sect. 2.4); and green regions represent the overlap of thirdorder resonances. The nonoverlapping thirdorder resonances are not indicated because they are very thin and thus unimportant for a global picture of the dynamics. Second row: as a comparison, the system given by Eq. (1) is integrated numerically with the same orbital model (quasiperiodic decomposition of Laskar 1990), and a frequency map analysis is performed to locate the chaotic zones. The colour scale goes from black (no chaos), to red (strong chaos). Third row: same maps obtained from a more detailed model in which the orbital evolution is directly taken from a numerical integration (adapted from Laskar & Robutel 1993). In the weakly chaotic zones, the dots are shifted vertically according to the level of chaos; in the strongly chaotic zones, they are plotted in boldface. Bottom row: same as top row, but the longterm orbital evolution of each planet is approximated by the LagrangeLaplace system (Sect. 3.1). The eight planets of the solar system are included, with the initial conditions of Bretagnon (1982). 
2.4 Higher order resonances
Outside of firstorder resonances, we can use a nearidentity canonical change of coordinates in order to suppress the angular dependency at first order (as already used in a similar context by Li & Batygin 2014b). Let us consider an intermediary Hamiltonian , such that the current coordinates are obtained from the new ones through its flow at time 1. The Hamiltonianin the new coordinates is then (27)
In these expressions, Poisson’s brackets are defined as (29)
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 must fulfil the homological equation (30)
in which is the zerothorder term of the multidimensional Fourier decomposition of (average of over all angles), which is here equal to zero. By matching the terms of the Fourier decomposition of and one by one, the solution to the homological equation is (31)
Injecting this function into the expression of the new Hamiltonian (Eq. (28)), we get as required, and the secondorder term is given in Appendix C. The only possible resonant angles at second order have the form σ = ϕ_{j} + ϕ_{k} + 2ψ. 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 secondorder resonances are computed in Appendix C and given in the first line of Table 1.
Outside of both firstorder and secondorder resonances, the same method can be used to compute the location and width of thirdorder resonances. An intermediary Hamiltonian of the form is used, in which must satisfy a second homological equation (see Appendix D). The possible resonant angles, as well as the centre and half width of all the thirdorder resonances are given in Table 1. As before, the same method can be used in the case of synchronous rotation, by adding the term to (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 thirdorder resonance. Resonances of order 2 or 3 with a centre located inside a resonance of lower order are not considered (in other words, loworder resonances are assumed to be 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 highorder resonances in this context. Indeed, numerous frequencies of the quasiperiodic representation are quite close to each other, which implies that the corresponding resonances overlap massively. Hence, even if the second and thirdorder resonances are thin, most of them are located one after another, resulting in large chaotic zones. The chaotic diffusion is however slower for higher order resonances. Moreover, 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 be facilitated by small modulations of the resonance locations.
We notethat thirdorder 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. 1993b, 2004a), has two different origins: one is a firstorder resonance with ν_{23}, and the other is a thirdorder resonance between μ_{1}, μ_{10} and ν_{12} (the numbering refers to Laskar et al. 2004a). The amplitude of the first resonance, which is the one emphasised in the literature, is larger by a factor of thirty. This resonance is present in the synthetic representation used in this paper (see the term with frequency − 50.30212″ yr^{−1} in Table F.1). In the top row of Fig. 3, it appears as a thin isolated resonance (upperleft 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 mixedtyperesonance 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 the whole planetary system can be damped even if only one planet dissipates energy with the star. Therefore, 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 loweccentricity and lowinclination regime, the longterm rotational dynamics ofplanets 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 quasiperiodic representation) and (ii) the rotation parameters (α coefficient). In this section, we outline how these quantities can be estimated from physical and dynamical arguments, even with scarce data.
Critical angle, location, and half width of every second and thirdorder secular spinorbit resonance.
3.1 The LagrangeLaplace system
Regarding the longterm orbital dynamics, one can use nominal orbital elements (either bestfit or assumed ones) and numerically integrate the equations of motion. The resulting 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 quasiperiodic 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. In fact we see below that at this level of precision, the LagrangeLaplace system is already a goodenough approximation of the orbital dynamics, up to moderate eccentricities and inclinations (and without strong effects coming from meanmotion resonances).It was used for the same purpose by Atobe et al. (2004) in the case of a massless hypothetical terrestrial planet.
The LagrangeLaplace system is the lowestorder model of the longterm orbital dynamics: it uses a development of the Hamiltonianat 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 (32)
where the index k = 1, 2, …, N represents a given planet of the system. Writing z and ζ the vectors of all z_{k} and ζ_{k}, the equations of motion in the LagrangeLaplace approximation are (33)
where the real matrices A and B are only functions of the masses and semimajor axes. They can be retrieved from the lowestorder terms in eccentricity and inclination of the orbital Hamiltonian. Organising the planets by increasing semimajor axes, we get from Laskar & Robutel (1995): (34)
in which , and the functions C_{2}(α) and C_{3} (α) are expressed in terms of the Laplace coefficients : (36)
(see Laskar & Robutel 1995, Laskar et al. 2012 or Murray & Dermott 1999). The equations of motion for z and ζ 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) and has the form of quasiperiodic series (Eq. (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 (37)
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), ζ_{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 μ_{j} and ν_{j} of the quasiperiodic 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 thus have 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 the bottom row of Fig. 3. 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 surrounding the elements of an exoplanetary system^{1}. Using this approach with the nominal orbital elements given by Brasser et al. (2014) and Deitrick et al. (2018), we retrieve from analytical formulas their maps showing the possible obliquity variations of HD 40307 g and Kepler62 f, in terms of the locations and widths of the secular spinorbit 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 LagrangeLaplace system, associated with the development of the Hamiltonian (Sect. 2), is enough to obtain the level of detail required for studying the longterm 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 the ascending node. From now on, we suppose that only the masses, the semimajor axes, and the eccentricities are known for all planets of the system. In this case, the LagrangeLaplace matrices A and B can still be computed since they only depend on the masses and the semimajor axes. We thus obtain the two sets of frequencies μ_{j} and ν_{j}. Because the initial conditions z_{i}(0) and ζ_{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 tothe available data.
For the eccentricity, this amounts to maximising the modulus of a sum of complex numbers with unknown phase. The result is therefore simply the sum of the moduli: (38)
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 ζ_{i} (0) are unknown. It is therefore necessary to introduce additional arguments, either from a physical or dynamical point of view. 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 multiplanet system seems simplistic. We opt here for the hypothesis by Laskar & Petit (2017) of equipartition of the angular momentum deficit (AMD) among the secular degrees of freedom. As these latter authors 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, equipartition of the AMD allows one 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: (39)
Contrary to Laskar & Petit (2017), we use here the Hamiltonian decomposition of Laskar & Robutel (1995), where the integrable part is the Sunplanet twobody problem. This is also the decomposition chosen when expressing the matrices A and B of the LagrangeLaplace system (Eqs. (34)–(35)). The AMD equipartition hypothesis amounts to considering that the total AMD of the system, (41)
Hence, even if the individual orbital inclinations are not known, the soobtained value of C gives a bound for them. For instance, the maximum possible value of the inclination of the kth planet is given by (43)
In our case, we are trying to maximise the quantity (44)
obtained from Eq. (37) with unknown Ω_{i}, using the constraint from Eq. (42). This constraint can be rewritten (45)
where η_{i} = sin(I_{i}∕2), and Z = C_{p}, whereas the quantity to be maximised can be written (46)
where . The coefficients c_{i} and b_{i} are all positive, and 0 ≤ η_{i} ≤ 1. Equation (45) forms a hyperellipsoid, whereas the quantity to be maximised (Eq. (46)) forms a hyperplane. Except from particular cases that we 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: (47)
where λ > 0 by definition of η_{i}. We get the value of λ from the imposed value of Z: (48)
The maximum of Y with the constraint Z is thus (49)
Going back to the original notations, this finally gives (50)
However, Eq. (47) does not take into account the condition that all the η_{i} are smaller than 1. In practice, if the value obtained for λ implies that one or several η_{i} are larger than 1, we simply have to fix them to 1 and use the same resolution method iteratively^{2} with the remaining η_{i} (changing the definition of Z and Y accordingly).
Table 2 shows the comparison of the amplitudes obtained for the Earth with the full LagrangeLaplace system, and their maximisation supposing that the mutual inclinations and longitudes of the 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. Therefore, taking all the maximum amplitudes at once, as if they formed one single representation, gives a large upper bound for the chaotic zones. Moreover, due to the large amplitudes of the series obtained, the smallwidth approximation for second and thirdorder resonances does not necessarily hold.
For simple systems, such 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 therefore enough to give a clear view of the dynamics. Using the maximised amplitudes gives the largest possible widths of the resonances, which, in turn shows the maximum obliquity variations and their locations. We fully retrieve their results. In contrast to these simple and ordereddynamics, Fig. 5 shows as an illustration of the maximised chaotic regions for the GJ 3293 system. The available orbital elements are taken from AstudilloDefru et al. (2017), who pointed out that GJ 3293 d is in the habitable zone. In the spindown 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, however, that the chaotic zones are maximised here according to the available orbital data. We also predict a rich obliquity dynamics for the Trappist1 planets, but due to the confirmed strong effects of meanmotion resonances (see e.g. Quarles et al. 2017), the use of the LagrangeLaplace model is probably inadequate in this case. Building an orbital theory specific to this system would be beyond the scope of this paper.
Secular representation of the Earth orbital dynamics given by the LagrangeLaplace theory.
Fig. 4 Estimate of the chaotic regions of the spinaxis dynamics for the Earth, where the longterm orbital dynamics is approximated by the LagrangeLaplace system. The same colour code as Fig. 3 is used. Left panel: the complete set of initial conditions is used (same as Fig. 3, bottom row). Right panel: the orbital elements (I, ϖ, Ω) are assumedto be unknown for all planets while the remaining ones are taken from Bretagnon (1982). Accordingly, the coefficients (E_{k}, S_{k}) of the quasiperiodic series are maximised according to the estimated AMD value (see Sect. 3.2). 
3.3 Maximisation of α
In Sects. 3.1 and 3.2, we saw how to obtain a quasiperiodic approximation of the longterm orbital motion of a planet, and how to estimate bounds for its coefficients if some of the orbital parameters are unknown. However, in order to study its longterm spin dynamics, we still lack an estimate of its precession constant α. Even in the solar system, the precession constants of the planets are not very well known. The estimate of α for an extrasolar planet would require observations that are very hard to obtain (its rotation period and a model of interior), and would be specific to one exoplanet. In order to keep this study as general as possible, we do not try to obtain a single value for the precession constant; instead, we look for its upper bound from general physical considerations.
After the sphere, the simplest shape model for a rotating planet is given by the Maclaurin ellipsoid (Chandrasekhar 1969). This latter describes the equilibrium shape of a selfgravitating homogeneous body in rotation with constant angular velocity. The rotational symmetry is imposed (circular equator), leading to the formula (51)
where ω and ρ are the rotation velocity and the density of the body and ɛ is the eccentricity of its ellipsoidal figure (in any plane containing the rotation axis). Studying as a function of the ellipsoid eccentricity, f is zero for ɛ = 0 and ɛ = 1, and it has one maximum at ɛ_{0} ≈ 0.929956 with the value f_{0} ≈ 0.224666. This implies that there is no such equilibrium figure possible for rotation velocities larger than (52)
Converting the ellipsoid eccentricity in terms of momenta of inertia, we obtain the relation (53)
Injecting this into the expression of the precession constant (Eq. (2)), we obtain the maximum value (54)
For rotation velocities close to ω_{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 α_{max} and the homogeneous approximation is quite crude anyway, allowing us to stick to Eq. (54). Planets are expected to spin much more slowly than ω_{max} (Eq. (52)), including giant gaseous planets (Batygin 2018). In the remainder of the article, ω_{max} is referred to as the “rotational breakup” velocity, meaning that above this limit the planet spins so fast that it experiences major structural damage.
Using the average density of the Earth, we obtain a minimum rotation period of about 2.4 h, leading to a maximum precession constant of about 230″ yr^{−1}. The true value for the Earth is 20 yr^{−1}, or 50″ yr^{−1}, 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 α_{max} ≈ 0.75″ yr^{−1}, while the true value is 0.20″ yr^{−1}, but it increases to 0.83″ yr^{−1} 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 have been observed so far. When using our model to study the spin dynamics of exoplanets, we must therefore 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 ρ 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 massradius distribution (see Fig. 6). The density is, in any case, not the major source of uncertainty in our method.
Fig. 5 Chaotic regions of the spinaxis dynamics for exoplanets of the GJ 3293 system, maximised with the method presented in this paper. The colour code is the same as in previous figures. The bounds for their precession constants are, from left to right panels: α_{max} = 66, 20, 6.2 and 1.2 deg yr^{−1} (corresponding to maximum rotation periods of a few hours). For each exoplanet, the horizontal line shows the precession constant corresponding to a rotation period equal to the orbital period (obtained from the method of Sect. 3.3). They are, from left to right: 13, 31, 48 and 123 days. 
3.4 Classification of nonresonant exoplanets
In Sect. 2.3, we saw that the overlap of firstorder secular resonances, leading to the largest chaotic zones of the spin dynamics, can only be produced if the ratio ν_{i} ∕α is smaller than 1 for at least two frequencies ν_{i} of the inclination quasiperiodic representation. Assuming that α is bounded by α_{max}, we can deduce that there can be no substantial chaotic zone if ν_{i}∕α_{max} > 1 whatever the frequency ν_{i}. Moreover, if the mutual inclinations are small (as we assume they are), this implies that the obliquity is almost constant. Our bound for α (Eq. (54)) can thus be used for a preliminary classification of the exoplanets, while the bounds for the amplitudes (Eqs. (38) and (50)) are required for a more specific application to one exoplanet (they allow to constrain both γ and β, see Fig. 1).
Table 3 shows the 94 planets classified strictly nonresonant with this criterion, using the LagrangeLaplace matrix to estimate the frequencies (Sect. 3.1) and Eq. (54) as a bound of α. All the exoplanets^{3} with knownmass, semimajor axis and eccentricity were analysed (taking m sinI instead of the mass if a realmass estimate was unavailable). On 20180307, 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 in this category. Such large semimajor 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 firstorder resonances (Sect. 2.3) because there is only one forced frequency in inclination (Sect. 3.1). However, this single term allows us 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 γ and β (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 β is small as well (Eq. (26)), resulting in quite flat level curves for Colombo’s top Hamiltonian (Eq. (18)). The maximum obliquity variations are achieved around Σ = 0, and at leading order, they are equal to (55)
This approximation holds very well for small values of S_{1}. For the exoplanet WASP81 c, which has the lowest bound for S_{1} in Table 3, we obtain a maximum obliquity excursion of 0.5^{o}. This limit is very close to what is obtained by plotting the level curves of the Hamiltonian (Eq. (18)), and it holds as long as α and S_{1} are below their estimated bounds. Such a good constraint cannot be achieved for every planet in twoplanet systems, however, 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 therefore also be used as an orderofmagnitude estimate. This results in a maximum of about 20^{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).
Fig. 6 Empirical massradius relationship obtained from the exoplanets with known mass and radius (http://exoplanet.eu). The exoplanets are supposed to be rocky up to 1 Earth mass, gaseous beyond 200 Earth masses, and of intermediate composition in between. Three power laws are used: 1∕3, 1∕2, and − 0.06 from left to right, in general accordance with e.g. Seager et al. (2007) and Weiss et al. (2013). 
4 Conclusion
The spinaxis 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 present an analytical formulation of the longterm spinaxis 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 quasiperiodic series. This method is therefore valid as long as the orbital chaos, if any, takes place on a much larger timescale than the spinaxis evolution. The spinaxis Hamiltonian is then expanded in powers of the eccentricity and inclination amplitudes of the orbital series. Here, we provide 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 with 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 α (gathering the physical characteristics of the planet under study) and to the quasiperiodic 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 for the previous numerical results to be analytically retrieved with good precision. Numerical integrations prove therefore 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. 2004b) 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 one to stick to the simple analytical formulas presented in this article; no further information would be brought by numerical integrations.
At this level of uncertainty, the LagrangeLaplace system provides a goodenough representation of the orbital motions (except for exoplanetary systems featuring highly excited orbits or strong effects of meanmotion resonances). The formulas obtained allow for an upper bound to be set on the amplitude of the eccentricity and inclination terms if the mutual orientations of the orbitsare 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 spinaxis 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 highly affected by meanmotion resonances (like Trappist1) 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 quasiperiodic series.
However, this method does not allow for tidal dissipations (playing an important role for exoplanets close to their star) to be considered, 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). Such a modelling amounts to making the precession constant α and/or the amplitudes of the orbital series gradually vary. This method does not include either the effects of libration around spinorbit resonances, even if a method exists that allows one 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 α. This bound is obtained from the flattening of the planet corresponding to its rotational breakup velocity. Since α governs the width and location of the resonances, this allows one to classify the exoplanets that cannot be subject to firstorder secular spinorbit 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 α is however invalidated by the possible presence of massive satellites (such as our Moon), but some exoplanets are so far from resonance that their classification is quite safe. This is the case for 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.
Exoplanets from http://exoplanet.eu on 20180307 classified nonresonant with the method detailed in this study.
Acknowledgements
We thank the anonymous referee for her or his detailed review. Funding and support from Paris Sciences et Lettres (PSL) university through the project origins and conditions for the emergence of life (OCAV) is acknowledged.
Appendix A Case of a 1:1 spinorbit 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 spinaxis dynamics, we can retrieve Colombo’s top Hamiltonian (Sect. 2.2), allowing us to use the same approach as in the nonresonant case. As before, though, we do not consider the effect of the tidal dissipation on the obliquity. This is therefore only valid for systems for which the tidal damping of the obliquity acts on a larger timescale than the spinaxis dynamics.
We usethe same method as Correia et al. (2003). We denote λ the mean longitude of the planet in orbit around the star, and ℓ its rotation angle. The mean longitude λ is measured from the equinox at a reference epoch (for instance J2000), whereas the rotation angle ℓ 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 ℓ − λ during the average over the mean longitude and the fast rotation angles (see Néron de Surgy & Laskar 1997), the corresponding “semiaveraged” Hamiltonian is (A.1)
where we neglected terms of order e(B − A)∕C. The momenta L = Cω and Y = LX are conjugate to ℓ and − ψ, respectively. The momentum Λ, conjugate to λ, has been added such that (mean motion). The resonant precession constant is defined as (A.2)
using the same notation as Eq. (2). We note that the angle λ + ψ appearing in the Hamiltonian corresponds to the mean longitude measured from the equinox of the date. Now we use the canonical change of coordinates (A.3)
The momentum Γ is an arbitrary constant of motion and the Hamiltonian becomes (A.4)
We now suppose that the dynamics of θ, corresponding to the “semisecular” timescale (either circulation or oscillation), is much faster than the evolution of the other degrees of freedom, corresponding to the secular timescale. We therefore consider for now that apart from (I, θ), all the variables are fixed (adiabatic approximation). The equations of motion are (A.5)
Using the definitions of I, Y, and α_{r}, the first equation gives (A.6)
resulting, for any value of X, in two equilibrium points: θ = ψ and ψ + π∕2 mod π. We note that θ = ψ is an elliptic equilibrium while θ = ψ + π∕2 is hyperbolic. Injecting this into the second equation, we obtain (A.7)
in which the small terms correspond to the precession of the spin axis (α and α_{r}) and the precession of the orbit ( and ). The equilibrium condition, corresponding to the exact resonance, is thus ω ≈ n. Considering that the planet is locked in synchronous rotation, we have therefore θ = ψ and ω ≈ 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: (A.8)
where this time, we use X as conjugate momentum of − ψ (the Hamiltonian is thus divided by the constant L). In the expression of α and α_{r}, we must replace ω by n. We get here oneextra term with respect to Eq. (1), due to the spinorbit resonance. Using the same method as in Sect. 2.2, the Hamiltonian in case of a firstorder secular spinorbit resonance is (A.9)
which must be compared to Eq. (16). This Hamiltonian has the same general form and can be reduced to Colombo’s top. We can therefore apply the same method of resolution (redefining the constants accordingly).
Appendix B Characteristic quantities of Colombo’s top
B.1 Equilibrium points
From Eq. (18), the equations of motion are (B.1)
Apart from the coordinate singularity at Σ = ±1, the first equation implies that when σ = 0 or π. Injecting this into the second equation, we get (B.2)
where β ≥ 0 by hypothesis. The resolution of this equation requires to square left and righthand terms, loosing the information^{4} about the sign of cosσ. We obtain a quartic equation in Σ: (B.3)
This discriminant is equal to zero for the particular cases γ = 0 or β =0, for which the polynomial can be factored into, respectively, (B.5)
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 γ > 0 and β > 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 (B.6)
where numerous intermediary variables are required in order to obtain compact expressions: (B.7)
We notethat Σ_{c,d} are real solutions only when Δ_{4} ≥ 0 (see below for the limit in terms of γ and β). The corresponding values of σ are (B.8)
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 Eq. (B.4) corresponds to a bifurcation. Its position can be computed by solving the equation Δ_{4} = 0, which corresponds to solving a cubic equation either in γ^{2} or β^{2}. Choosing to solve it in terms of β, the discriminant is (B.9)
meaning that there is only one real solution. This solution is (B.10)
which is the boundary (Eq. (20)).
B.3 Second boundary (A/B)
The other two boundaries can be obtained by studying the level curves of the Hamiltonian passing through Σ = ±1 (which is singular using the coordinates Σ and σ, but it does not matter here).
Let us begin with the + 1 case, for which the Hamiltonianhas a value of − 1∕2 + γ. We now look for this specific level curve along the axes σ = 0 and σ = π. This leads to the equation (B.11)
for which Σ = +1 is a solution. By reorganising the terms, taking the square (thus loosing the information about the sign of cos σ), and dividing by (Σ − 1), we get (B.12)
which is a cubic equation in Σ. Its determinant is (B.13)
Once again, it is zero for β = 0. Moreover the solutions for γ = 0 can be easily computed. In these two particular cases, the polynomial can be factored into, respectively, (B.14)
showing the solutions and their multiplicities. For γ > 0 and β > 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 Δ_{3} = 0, which amounts to solving a quadratic equation in β^{2} or a quartic equation in γ. Choosing to solve it in terms of β, the only positive solution is (B.15)
which is the boundary (Eq. (25)).
B.4 Third boundary (B/C)
Let us nowstudy the level curve of the Hamiltonian passing in Σ =−1, which has value − 1∕2 − γ. The procedure is the same as for the second boundary, and the new formulas are obtained simply by replacing γ with − γ. There is however an ambiguity because there are two positive solutions β^{2} (as a function of γ) which cancel the determinant. The one corresponding to the bifurcation is the largest, that is, (B.16)
which is the boundary (Eq. (25)).
B.5 Separatrices
The position at σ = 0 or π of the separatrix emerging from the hyperbolic point (Σ, σ) = (Σ_{c}, π) defines the boundaries of the resonant region (see Fig. 1). Writing , the equations to solve are (B.17)
The resolution of this equation requires to square left and righthand terms, loosing the information^{5} about the sign of cosσ. We obtain a quartic equation in Σ, (B.18)
in which Σ_{c} is a double root. It can thus be divided by , leading to the quadratic equation (B.19)
This equation always has two real solutions, provided that Σ_{c} exists (i.e. in zones A, B, or C). These solutions are (B.20)
where we replaced f with Eq. (18) in terms of Σ_{c}.
Appendix C Secondorder resonances
Using the intermediary Hamiltonian (Eq. (31)), the Hamiltonian in the new coordinates is obtained term by term from Eq. (28). The two first terms are simple: we have (given at Eq. (11)) and by definition of . The secondorder term is more complex since it requires computation of Poisson’s brackets. Using the fact that and reorganising the terms adequately, we obtain (C.1)
Since by hypothesis there is no firstorder resonance in the system, the only possible resonant angles at second order are of the form σ = ϕ_{j} + ϕ_{k} + 2ψ. Let us perform the canonical change of coordinates (C.2)
Assuming that σ is the only resonant angle, the dynamics at second order is given by averaging over all other angles (this is another change of coordinates close to identity). The momenta Γ_{1} and Γ_{2} become arbitrary constants of motion that we conveniently choose equal to zero. Dropping the unnecessary constants, the resonant Hamiltonian is thus (C.4)
in which X must be replaced by − 2Σ. Highorder resonances are relatively thin, so it is enough to consider the dynamics in the neighbourhood of the resonance centre at first order: (C.5)
or equivalently X_{0} = − 2Σ_{0} = − (ν_{j} + ν_{k})∕(2α). Considering that , we therefore obtain (C.6)
in which we dropped the unnecessary constants, and where (C.7)
By injecting the momentum Σ instead of X and by using the modifiedtime dτ = −4αdt, we obtain (C.8)
This is the Hamiltonian of a pendulum of centre Σ_{0} and half width . In terms of the obliquity cosine X, the resonance has position X_{0} and half width .
Appendix D Thirdorder 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 , such that the new coordinates are given by its flow at time 1. The Hamiltonian in the new coordinates is then (D.1)
The firstorder part of required to suppress the angular dependency at order 1 can be directly taken from Eq. (31). We write (D.3)
(average plus oscillating part) for which the expression is given by Eq. (C.1). The homological equation for order 2 is then (D.4)
which defines the Hamiltonian . This leads to (D.5)
We must now compute the remainders at order 3. First of all, we can simplify their expressions by taking into account that, by definition: , and . We have then (D.6)
The expression of the coefficients D_{j} is very complex. We do not give them here as they are of no interest at this stage (the angles ϕ_{j} + ψ are nonresonant by hypothesis).
As shown in Appendix C, in the pendulum approximation, the halfwidth of any possible resonance is two times the square root of its coefficient divided by α, 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 LagrangeLaplace 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 × N matrix with elements (b_{ij}). The ith “Geršgorin disc” (i = 1, 2, …, N) is the closed disc of the complex plane centred at b_{ii} and with radius (E.1)
Fig. E.1 Geršgorin discs in the complex plane corresponding to the LagrangeLaplace matrix B for three planets. The centre of the circles are the diagonal entries of B (red spots). Every eigenvalue of B lies on the real line, inside the union of all the discs. 
Theorem (Geršgorin 1931). Any eigenvalue of B lies inside at least one of the discs, i = 1, 2, …, N.
Corollary. All the eigenvalues of B are located inside the union of the 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. E.1.
Appendix F Orbital solution used for the inner solar system
In order to apply our method to a given planet, we first need a quasiperiodic approximation of its longterm 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 and (Tables VI and VII of Laskar 1990) by the matrix corresponding to the linear part of the solution (Table V of Laskar 1990). In the series obtained, the terms with the same combinationof frequencies are then merged together, finally resulting in 56 terms in eccentricity and 60 terms in inclination.
These series are given in Table F.1 for the inner planets, under the form: (F.1)
with N = 56 and M = 60. They are used in Fig. 3 of the present work.
Quasiperiodic representation of the orbital dynamics of Mercury.
References
 Armstrong, J. C., Barnes, R., DomagalGoldman, S., et al. 2014, Astrobiology, 14, 277 [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Atobe, K., Ida, S., & Ito, T. 2004, Icarus, 168, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K. 2018, AJ, 155, 178 [NASA ADS] [CrossRef] [Google Scholar]
 Boué, G., & Laskar, J. 2006, Icarus, 185, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Boué, G., & Laskar, J. 2010, ApJ, 712, L44 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., Ida, S., & Kokubo, E. 2014, MNRAS, 440, 3685 [NASA ADS] [CrossRef] [Google Scholar]
 Bretagnon, P. 1982, A&A, 114, 278 [NASA ADS] [Google Scholar]
 Canup, R. M., & Asphaug, E. 2001, Nature, 412, 708 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, J. A., & Winn, J. N. 2010, ApJ, 716, 850 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (New Haven: Yale University Press) [Google Scholar]
 Colombo, G. 1966, AJ, 71, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M. 2014, A&A, 570, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003, Icarus, 163, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Laskar, J., & de Surgy O. N. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Deitrick, R., Barnes, R., Quinn, T. R., et al. 2018, AJ, 155, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Geršgorin, S. 1931, Bulletin de l’Académie des Sciences de l’URSS, Classe des sciences mathématiques et naturelles, 6, 749 [Google Scholar]
 Hartmann, W. K., & Davis, D. R. 1975, Icarus, 24, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Hays, J. D., Imbrie, J., & Shackleton, N. J. 1976, Science, 194, 1121 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Henrard, J., & Murigande, C. 1987, Celest. Mech., 40, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1988, A&A, 198, 341 [NASA ADS] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1994, A&A, 287, L9 [NASA ADS] [Google Scholar]
 Laskar, J. 1996, Celest. Mech. Dyn. Astron., 64, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 2008, Icarus, 196, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Petit, A. C. 2017, A&A, 605, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., & Robutel, P. 1993, Nature, 361, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Laskar, J., Joutel, F., & Robutel, P. 1993a, Nature, 361, 615 [Google Scholar]
 Laskar, J., Joutel, F., & Boudin, F. 1993b, A&A, 270, 522 [NASA ADS] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004a, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., Correia, A. C. M., Gastineau, M., et al. 2004b, Icarus, 170, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, G., & Batygin, K. 2014a, ApJ, 795, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Li, G., & Batygin, K. 2014b, ApJ, 790, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Barnes, J. W., & Chambers, J. E. 2012, Icarus, 217, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Lock, S. J., Stewart, S. T., Petaev, M. I., et al. 2018, J. Geophys. Res. Planets, 123, 910 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Néron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [Google Scholar]
 Peale, S. J. 1969, AJ, 74, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Quarles, B., Quintana, E. V., Lopez, E., Schlieder, J. E., & Barclay, T. 2017, ApJ, 842, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., Kuchner, M., HierMajumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Shan, Y., & Li, G. 2018, AJ, 155, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596 [Google Scholar]
 Varga, R. S. 2004, Geršgorin and His Circles (Berlin: SpringerVerlag) [CrossRef] [Google Scholar]
 Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501 [NASA ADS] [CrossRef] [Google Scholar]
 Weertman, J. 1976, Nature, 261, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, J.W., Dong, S., Zhu, Z., et al. 2016, Proc. Natl. Acad. Sci., 113, 11431 [NASA ADS] [CrossRef] [Google Scholar]
Some subtle dynamical effects are not reproduced by the LagrangeLaplace system, like the s_{6} + g_{5} − g_{6} firstorder resonance mentioned in Sect. 2.4.
All Tables
Critical angle, location, and half width of every second and thirdorder secular spinorbit resonance.
Secular representation of the Earth orbital dynamics given by the LagrangeLaplace theory.
Exoplanets from http://exoplanet.eu on 20180307 classified nonresonant with the method detailed in this study.
All Figures
Fig. 1 Parameter space of Colombo’s top Hamiltonian (Eq. (18)). There is no separatrix (and thus no resonance) in region D, delimited by the curve (Eq. (20), thick black line). The curves and (Eq. (25), thin black lines) delimit regions A–B and B–C, respectively. In region A, the resonance island lies between Σ_{−} and Σ_{+} (Eq. (24)); in region B, it lies between Σ_{−} and + 1; in region C, it lies between − 1 and + 1 (hence the 180^{o} full width given by the colour scale). The problem becomes unphysical above β = 2γ (blue line) since it corresponds to an amplitude S_{p} > 1 in the inclination series. 

In the text 
Fig. 2 Examples of phase portraits for the four regions of Fig. 1. The equilibrium points are labelled as in Eq. (21). The level curves of the Hamiltonian are drawn with black lines out of the resonance island and with red lines inside the resonance. The parameters chosen are, from A to D: (γ, β) = (0.4, 0.1); (0.55, 0.15); (0.05, 0.7); and (0.8, 0.3). 

In the text 
Fig. 3 Top row: estimate of the chaotic regions of the spin dynamics as the superposition of secular spinorbit resonances. The orbitalevolution of each planet is approximated by the synthetic representation of Laskar (1990), as detailed in Appendix F. Lightred and darkred regions represent the firstorder resonances and their overlaps, respectively(Colombo’s top Hamiltonian, Sect. 2.2); lightblue and darkblue regions represent the secondorderresonances and their overlap (Sect. 2.4); and green regions represent the overlap of thirdorder resonances. The nonoverlapping thirdorder resonances are not indicated because they are very thin and thus unimportant for a global picture of the dynamics. Second row: as a comparison, the system given by Eq. (1) is integrated numerically with the same orbital model (quasiperiodic decomposition of Laskar 1990), and a frequency map analysis is performed to locate the chaotic zones. The colour scale goes from black (no chaos), to red (strong chaos). Third row: same maps obtained from a more detailed model in which the orbital evolution is directly taken from a numerical integration (adapted from Laskar & Robutel 1993). In the weakly chaotic zones, the dots are shifted vertically according to the level of chaos; in the strongly chaotic zones, they are plotted in boldface. Bottom row: same as top row, but the longterm orbital evolution of each planet is approximated by the LagrangeLaplace system (Sect. 3.1). The eight planets of the solar system are included, with the initial conditions of Bretagnon (1982). 

In the text 
Fig. 4 Estimate of the chaotic regions of the spinaxis dynamics for the Earth, where the longterm orbital dynamics is approximated by the LagrangeLaplace system. The same colour code as Fig. 3 is used. Left panel: the complete set of initial conditions is used (same as Fig. 3, bottom row). Right panel: the orbital elements (I, ϖ, Ω) are assumedto be unknown for all planets while the remaining ones are taken from Bretagnon (1982). Accordingly, the coefficients (E_{k}, S_{k}) of the quasiperiodic series are maximised according to the estimated AMD value (see Sect. 3.2). 

In the text 
Fig. 5 Chaotic regions of the spinaxis dynamics for exoplanets of the GJ 3293 system, maximised with the method presented in this paper. The colour code is the same as in previous figures. The bounds for their precession constants are, from left to right panels: α_{max} = 66, 20, 6.2 and 1.2 deg yr^{−1} (corresponding to maximum rotation periods of a few hours). For each exoplanet, the horizontal line shows the precession constant corresponding to a rotation period equal to the orbital period (obtained from the method of Sect. 3.3). They are, from left to right: 13, 31, 48 and 123 days. 

In the text 
Fig. 6 Empirical massradius relationship obtained from the exoplanets with known mass and radius (http://exoplanet.eu). The exoplanets are supposed to be rocky up to 1 Earth mass, gaseous beyond 200 Earth masses, and of intermediate composition in between. Three power laws are used: 1∕3, 1∕2, and − 0.06 from left to right, in general accordance with e.g. Seager et al. (2007) and Weiss et al. (2013). 

In the text 
Fig. E.1 Geršgorin discs in the complex plane corresponding to the LagrangeLaplace matrix B for three planets. The centre of the circles are the diagonal entries of B (red spots). Every eigenvalue of B lies on the real line, inside the union of all the discs. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.