Secular spin-axis dynamics of exoplanets

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).


Introduction
From the works by Laskar & Robutel (1993) and Laskar et al. (1993a), 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 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 spin-axis evolution of the Earth (Néron de Surgy & Laskar 1997;Laskar et al. 2004a;Li & Batygin 2014a), Venus , 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 high-obliquity 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 spin-axis 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 time-dependent 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 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 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 spin-axis 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 to the 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 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. 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 "non-resonant" exoplanets, for which no chaos can appear and for which the obliquity variations are constrained by an analytical bound.

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 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). This Hamiltonian can be written 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 where 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, ω 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 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 (two-body problem), the obliquity is constant and the precession angle circulates with constant angular velocity αX/(1 − e 2 ) 3/2 . 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 C = η 2Ω = O(η 2 ), we obtain We now consider 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 A4, page 2 of 21 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, 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(Laskar , 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 As shown in Appendix A, the Hamiltonian has the same form in the case of a 1:1 spin-orbit 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τ = a α dt to Eq. (16), we obtain the following Hamiltonian (that we still denote F ):  (18)). There is no separatrix (and thus no resonance) in region D, delimited by the curve C 1 (Eq. (20), thick black line). The curves C 2 and C 3 (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.
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 Below C 1 (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 C 1 , and disappear above it (region D of Fig. 1). Their respective positions are 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 for β 1 : and for γ 1 : 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. 2004or 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: 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 (see Appendices B.3 and B.4), delimiting the regions A-B and B-C of Fig. 1, respectively. We note that C 1 and C 2 intersect at (γ; β) = (1 ; 0), C 1 and C 3 intersect at (1/8 ; 3 √ 3/8), and C 2 and C 3 intersect at (0 ; 1/2). Contrary to C 1 , the boundaries C 2 and 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 [Σ − ; Σ + ] 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 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 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 β 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 first-order 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 first-order 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).

Overlap of first-order 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 quasi-periodic 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 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 β coefficients, resulting in virtually no secular variation of the obliquity (see the shades of grey in Fig. 1). A no-chaos 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 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 (1) is integrated numerically with the same orbital model (quasi-periodic 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 long-term orbital evolution of each planet is approximated by the Lagrange-Laplace system (Sect. 3.1). The eight planets of the solar system are included, with the initial conditions of Bretagnon (1982).
A4, page 6 of 21 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 therefore dedicated to second-and third-order resonances.

Higher order resonances
Outside of first-order resonances, we can use a near-identity 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 X = ε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 theñ wherẽ In these expressions, Poisson's brackets are defined as 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 X must fulfil the homological equatioñ in which H 1 is the zeroth-order term of the multidimensional Fourier decomposition of H 1 (average of H 1 over all angles), which is here equal to zero. By matching the terms of the Fourier decomposition of H 1 and X 1 one by one, the solution to the homological equation is Injecting this function into the expression of the new Hamiltonian (Eq. (28)), we getH 1 = 0 as required, and the secondorder termH 2 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 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 X = εX 1 + ε 2 X 2 is used, in which 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 −α r (1 + X) 2 /2 to 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 resonance. 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 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 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 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 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. 1993b(Laskar et al. , 2004a, has two different origins: one is a first-order resonance with ν 23 , and the other is a third-order 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 (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 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.

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 (α coefficient). In this section, we outline how these quantities can be estimated from physical and dynamical arguments, even with scarce data. Table 1. Critical angle, location, and half width of every second-and third-order secular spin-orbit resonance.
Notes. In the fourth line, the symbol P stands for:

The Lagrange-Laplace system
Regarding the long-term orbital dynamics, one can use nominal orbital elements (either best-fit 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, andDeitrick 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. In fact we see below 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 z k = e k exp(i k ) and ζ k = sin 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 Lagrange-Laplace approximation arė 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): and in which n 2 j a 3 j = G(m 0 + m j ), and the functions C 2 (α) and C 3 (α) are expressed in terms of the Laplace coefficients b (k) s : (see Laskar & Robutel 1995, Laskar et al. 2012or 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 quasi-periodic 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 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 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 thus have M = N − 1 terms, whereas the eccentricity A4, page 8 of 21 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 Kepler-62 f, in terms of the locations and widths of the secular spin-orbit resonances (see Fig. 8 by Brasser et al. 2014 andFigs. 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.

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 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 µ 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 to the 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: 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 multi-planet 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(Laskar , 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: where 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 decomposition 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, is equal to 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 kth planet is given by In our case, we are trying to maximise the quantity obtained from Eq. (37) with unknown Ω i , using the constraint from Eq. (42). This constraint can be rewritten where η i = sin(I i /2), c i = 2Λ i 1 − e 2 i and Z = C p , whereas the quantity to be maximised can be written where b i = Q −1 ki . The coefficients c i and b i are all positive, and 0 η i 1. Equation (45) forms a hyper-ellipsoid, 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: where λ > 0 by definition of η i . We get the value of λ from the imposed value of Z: The maximum of Y with the constraint Z is thus Going back to the original notations, this finally gives 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 Lagrange-Laplace 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 small-width approximation for second-and third-order 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 ordered dynamics, Fig. 5 shows as an illustration of 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 2 From the form of the constraint (Eq. (45)), decreasing one η i to 1 implies that at least one of the remaining η i should increase; from the solution (Eq. (47)), this actually means that all the remaining η i increase. Notes. The eight planets of the Solar System are included. In the third column, each amplitude is maximised in the case where both the mutual inclinations and the longitudes are unknown, assuming the equipartition of AMD between secular degrees of freedom (since the Solar System is hierarchically AMD stable, the AMD of the inner and outer parts were taken separately, see Laskar & Petit 2017). The initial conditions and physical parameters are taken from Bretagnon (1982).   Fig. 3, bottom row). Right panel: the orbital elements (I, , Ω) are assumed to 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). 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 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 beyond the scope of this paper.

Maximisation of α
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 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 self-gravitating homogeneous body in rotation with constant angular velocity. The rotational symmetry is imposed (circular equator), leading to the formula 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 f = ω 2 /(2πGρ) 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 Converting the ellipsoid eccentricity in terms of momenta of inertia, we obtain the relation Injecting this into the expression of the precession constant (Eq. (2)), we obtain the maximum value 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 mass-radius distribution (see Fig. 6). The density is, in any case, not the major source of uncertainty in our method.

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 only be produced if the ratio ν i /α is smaller than 1 for at least two frequencies ν i of the inclination quasi-periodic representation. Assuming that α is bounded by α max , we can deduce that there can be no substantial chaotic zone if ν i /α max > 1 whatever A4 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 non-resonant with this criterion, using the Lagrange-Laplace matrix to estimate the frequencies (Sect. 3.1) and Eq. (54) as a bound of α. All the exoplanets 3 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). On 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 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 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 3 http://exoplanet.eu leading order, they are equal to 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 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 two-planet 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 order-of-magnitude 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).

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 present 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 therefore 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. 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 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 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;Laskar et al. 2004b) but not yet for exoplanets because their initial Notes. The first column gives the name of the exoplanet; the second column gives the rank of the exoplanet (sorted by increasing semi-major axis) and the total number of planets in the system; the third column gives the semi-major axis value; the fourth column gives the maximum value of the precession constant estimated from Eq. (54); the fifth column gives the minimum ratio of the eigenfrequencies of the Lagrange-Laplace system and α max in absolute value (γ parameter of Colombo's top); the sixth column gives the maximum amplitude of the series decomposition obtained from the maximisation (Eq. (50)), allowing to obtain a maximum bound for the β parameter of Colombo's top.
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 Lagrange-Laplace system provides a good-enough representation of the orbital motions (except for exoplanetary systems featuring highly excited orbits or strong effects of mean-motion 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 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 highly 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 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 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 α 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.

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 us to use the same approach as in the non-resonant 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 spin-axis dynamics.
We use the same method as . 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 "semi-averaged" Hamiltonian is 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λ = n (mean motion). The resonant precession constant is defined as 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 The momentum Γ is an arbitrary constant of motion and the Hamiltonian becomes We now suppose that the dynamics of θ, 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 therefore consider for now that apart from (I, θ), all the variables are fixed (adiabatic approximation). The equations of motion arė (A.5) Using the definitions of I, Y, and α r , the first equation giveṡ 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 obtaiṅ in which the small terms correspond to the precession of the spin axis (α and α r ) and the precession of the orbit (A and B). 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: 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 one extra term with respect to Eq. (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 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).
cos σ . Apart from the coordinate singularity at Σ = ±1, the first equation implies thatΣ = 0 when σ = 0 or π. Injecting this into the second equation, we get where β 0 by hypothesis. The resolution of this equation requires to square left-and right-hand terms, loosing the information 4 about the sign of cos σ. We obtain a quartic equation in Σ: This discriminant is equal to zero for the particular cases γ = 0 or β = 0, for which the polynomial can be factored into, respectively, 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 where numerous intermediary variables are required in order to obtain compact expressions: We note that Σ c,d are real solutions only when ∆ 4 0 (see below for the limit in terms of γ and β). The corresponding values of σ are 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 meaning that there is only one real solution. This solution is which is the boundary C 1 (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 Hamiltonian has a value of −1/2 + γ. We now look for this specific level curve along the axes σ = 0 and σ = π. This leads to the equation 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 which is a cubic equation in Σ. Its determinant is ∆ 3 = β 2 −β 4 + 1 4 − 5γ − 2γ 2 β 2 + γ(1 − γ) 3 . (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, 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 β 2 = 1 8 1 − 20γ − 8γ 2 + (1 + 8γ) 3/2 , (B.15) which is the boundary C 2 (Eq. (25)).

B.4. Third boundary (B/C)
Let us now study 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, which is the boundary C 3 (Eq. (25)).

(B.19)
This equation always has two real solutions, provided that Σ c exists (i.e. in zones A, B, or C). These solutions are where we replaced f with Eq. (18) in terms of Σ c .
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 α, 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 1931or 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 i j ). The ith "Geršgorin disc" G i (i = 1, 2, . . . , N) is the closed disc of the complex plane centred at b ii and with radius Theorem (Geršgorin 1931). Any eigenvalue of B lies inside at least one of the G i discs, i = 1, 2, . . . , N.
Corollary. All the eigenvalues of B are located inside the union of the 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. 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 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 and ζ • i (Tables VI and VII of Laskar 1990) by the matrix 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 Table F.1 for the inner planets, under the form: