Issue 
A&A
Volume 662, June 2022



Article Number  L3  
Number of page(s)  16  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202243327  
Published online  03 June 2022 
Letter to the Editor
The origin of chaos in the Solar System through computer algebra
IMCCE, CNRS, Observatoire de Paris, Université PSL, Sorbonne Université, 77 avenue DenfertRochereau, 75014 Paris, France
email: federico.mogavero@obspm.fr
Received:
14
February
2022
Accepted:
8
April
2022
The discovery of the chaotic motion of the planets in the Solar System dates back more than 30 years. Still, no analytical theory has satisfactorily addressed the origin of chaos so far. Implementing canonical perturbation theory in the computer algebra system TRIP, we systematically retrieve the secular resonances at work along the orbital solution of a forced longterm dynamics of the inner planets. We compare the time statistic of their halfwidths to the ensemble distribution of the maximum Lyapunov exponent and establish dynamical sources of chaos in an unbiased way. New resonances are predicted by the theory and checked against direct integrations of the Solar System. The image of an entangled dynamics of the inner planets emerges.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability / chaos
© F. Mogavero and J. Laskar 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
The chaotic longterm behaviour of the planetary orbits in the inner Solar System (ISS) emerged when the numerical integration of analytically averaged equations of motion revealed a maximum Lyapunov exponent (MLE) of about (5 million yr)^{−1} (Laskar 1985, 1989). Previously, the existence of secular resonances among the precession frequencies of planet perihelia and nodes had been shown to generate small divisors, which prevent the representation of the orbits as quasiperiodic series (Laskar 1984, 1988). Investigating the origin of chaos, Laskar measured the libration period of the Fourier harmonic θ_{2 : 1} = 2(g_{3} − g_{4})−(s_{3} − s_{4}), which involves the fundamental frequencies of Earth and Mars, to be about 4.6 Myr (Laskar 1990a). He then proposed the libration–circulation transitions of the corresponding argument as a source of the observed MLE. Evidence of a large chaotic zone bridging the resonances θ_{2 : 1} and θ_{1 : 1} = (g_{3} − g_{4})−(s_{3} − s_{4}) was later given (Laskar 1992). Nevertheless, when the integration of the full equations of motion confirmed the MLE of the averaged ones, the chaotic nature of the θ_{2 : 1} dynamics was questioned (Sussman & Wisdom 1992). As a consequence, a claim remains in literature that an undisputed dynamical mechanism for the observed chaos is missing (Lecar et al. 2001; Murray & Holman 2001; Hayes 2007). In the meantime, the alternating librations of θ_{2 : 1} and θ_{1 : 1} have been confirmed by direct integrations (Laskar et al. 2004) and supported by geological records (Ma et al. 2017; Olsen et al. 2019; Zeebe & Lourens 2019).
The highdimensional dynamics of the inner planets probably discouraged systematic analytical studies of its resonant structure. A couple of analyses have focused on the longterm motion of Mercury, by freezing all the other planets on quasiperiodic orbits (Lithwick & Wu 2011; Batygin et al. 2015), but this simplification leads to predictions that conflict with the findings of realistic models (Mogavero & Laskar 2021). Nevertheless, a Trojan horse against the curse of dimensionality affecting the ISS dynamics is offered by computer algebra, which allows the formal manipulation of the analytical series of celestial mechanics and the implementation of canonical perturbation theory in particular. Computer algebra has produced some remarkable results, such as the reproduction of Delaunay’s monumental lunar theory (Deprit et al. 1970) and the application of the Kolmogorov–Arnold–Moser theory to the threebody problem (Robutel 1995; Locatelli & Giorgilli 2000), in addition to the demonstration of the chaotic behaviour of the Solar System itself (Laskar 1985, 1989). Still, its use in celestial mechanics may seem limited given its potential.
We have recently proposed a forced model of the longterm dynamics in the ISS (Mogavero & Laskar 2021). It allows the secular phase space of the Solar System to be restricted in a consistent way to the eight degrees of freedom (DOFs) dominated by the inner planets. In this study we employ the computer algebra system TRIP (Gastineau & Laskar 2011, 2021) to carry out an unbiased analysis of the Fourier harmonics that constitute its Hamiltonian (Appendix A).
2. Forced secular inner Solar System
The longterm dynamics of the Solar System planets essentially consists of the slow precession of their perihelia and nodes, driven by secular, that is, orbitaveraged, gravitational interactions (Laskar 1990a; Laskar et al. 2004). The precession frequencies of the outer planet orbits are practically constant over billions of years when compared to those of the ISS (Laskar 1990a; Laskar et al. 2004; Hoang et al. 2021). Built on these facts, the model of forced secular ISS consists in predetermining a quasiperiodic secular solution for the giant planets, with the inner ones moving in the resulting timedependent gravitational potential (Mogavero & Laskar 2021). The quasiperiodic form of the giant planet orbits is established through frequency analysis (Laskar 2005) of the orbital solution of a comprehensive model of the Solar System (Laskar et al. 2004). The low planetary masses and the absence of strong meanmotion resonances in the ISS allow us to simply consider firstorder secular averaging of the Nbody Hamiltonian. This corresponds to Gauss’s dynamics of Keplerian rings (Gauss 1818), which we correct for the leading secular contribution of general relativity. The pertinence of our model has been thoroughly demonstrated (Mogavero & Laskar 2021). It matches reference orbital solutions of the Solar System over timescales shorter than or comparable to the Lyapunov time, with an average discrepancy in the fundamental frequencies of only a few hundredths of an arcsecond per year over the next 20 Myr. Moreover, it correctly reproduces the MLE and the statistics of the high eccentricities of Mercury over the next 5 billion years (Gyr).
Dynamical model. The secular Hamiltonian of the Solar System planets at first order in planetary masses and corrected for the leading contribution of general relativity reads (Mogavero & Laskar 2021)
The planets are indexed in order of increasing semimajor axis, are the planetary masses, a_{k} and e_{k} are the (secular) semimajor axes and eccentricities, respectively, G is the gravitational constant, and c is the speed of light. The vectors r_{k} are the planet heliocentric positions, and the bracket operator represents the averaging over the mean longitudes of the planets λ_{k}, which results from the suppression of the nonresonant Fourier harmonics of the Nbody Hamiltonian at first order in planetary masses (Mogavero & Laskar 2021). The semimajor axes a_{k} are constants of motion in the secular dynamics. A suitable set of canonically conjugate momentumcoordinate pairs of variables for the secular dynamics are the Poincaré rectangular coordinates in complex form, , with
where Λ_{k} = μ_{k}[G(m_{0} + m_{k})a_{k}]^{1/2}, m_{0} and μ_{k} = m_{0}m_{k}/(m_{0} + m_{k}) are the Sun mass and the reduced masses of the planets, respectively, ℐ_{k} are the planet inclinations, ϖ_{k} are the longitudes of the perihelia, and Ω_{k} are the longitudes of the nodes (Poincaré 1896; Laskar 1991; Laskar & Robutel 1995).
The model of forced ISS consists of the choice of an explicit quasiperiodic time dependence for the orbits of the outer planets (Mogavero & Laskar 2021),
for k ∈ {5, 6, 7, 8}, where t denotes the time, and are complex amplitudes, m_{kℓ} and n_{kℓ} are integer vectors, and ϕ(t) = ω_{o}t, with ω_{o} = (g_{5}, g_{6}, g_{7}, g_{8}, s_{6}, s_{7}, s_{8}) representing the septuple of the constant fundamental frequencies of the outer orbits (Laskar 1990a). Gauss’s dynamics of the forced ISS is obtained by substituting this predetermined time dependence into Eq. (1),
The resulting Hamiltonian system consists of two DOFs for each inner planet, corresponding to the x_{k} and y_{k} variables, respectively. The forced secular ISS is thus characterised by eight DOFs. As a result of the forcing from the outer planets, its orbital solutions live in a 16dimensional phase space since no trivial integrals of motion, such as the total energy or angular momentum, exist.
The development of the twobody perturbing function (Laskar & Robutel 1995), when implemented in TRIP, allows the Hamiltonian ℋ to be systematically expanded by exploiting the low eccentricities and inclinations of the planets (Mogavero & Laskar 2021). This development provides truncated Hamiltonians ℋ_{2n} that are multivariate polynomials of total degree 2n in the Poincaré variables of the inner planets (Appendix B). At the lowest degree, ℋ_{2} produces a forced LaplaceLagrange (LL) dynamics, which can be analytically integrated by introducing complex proper mode variables, . By introducing actionangle variables through and , the truncated Hamiltonians expressed in the proper modes can be expanded as finite Fourier series,
where I = (X, Ψ) and θ = (χ, ψ) are the eightdimensional vectors of the action and angle variables, respectively, and (k, ℓ) is the wave vector of a given harmonic. There are 2748 harmonics with a nonnull amplitude at degree four and more than ten million at degree ten (Mogavero & Laskar 2021).
3. Maximum Lyapunov exponent
Computing the MLE is fundamental to the determination of the origin of chaos, as its value can be compared to the halfwidth of the leading resonant harmonics of the Hamiltonian, which constitute the dynamical sources of chaoticity (Chirikov 1979). The nonnull probability of unstable orbital evolutions in the ISS (Laskar & Gastineau 2009) makes the definition of the MLE as an infinitetime limit (Oseledec 1968; Benettin et al. 1980) not pertinent. We numerically compute a finitetime MLE (FTMLE) employing the standard algorithm of Benettin et al. (1980) (faster chaoticity detectors have been developed starting from the fast Lyapunov indicator of Froeschlé et al. 1997). The FTMLE is timeasymptotically a stochastic function of the initial conditions of the system, and its computation acquires full physical significance for an ensemble of orbital solutions (Mogavero & Laskar 2021). Manipulation of the truncated Hamiltonians ℋ_{2n} in TRIP allows us to systematically derive the equations of motion and the corresponding variational equations, which we integrate through an Adams PECE method of order 12 and a timestep of 250 yr. Figure 1 shows the FTMLE expressed as an angular frequency over the next 5 Gyr for different degrees of truncation of the Hamiltonian. In each case, the FTMLE is computed for 128 stable (i.e. noncollisional) solutions, with initial conditions very close to the nominal values of Gauss’s dynamics and for different initial tangent vectors (Appendix C). The figure shows the [5th, 95th] percentile range of the probability distribution function (PDF) of the FTMLE estimated from each ensemble of solutions. We also report the PDF from the full Hamiltonian ℋ, computed in Mogavero & Laskar (2021), along with the FTMLE of its nominal solution 𝒮 for nine different initial tangent vectors, to manifest its asymptotic behaviour. In a few hundred million years, each FTMLE becomes independent of the initial tangent vector, the renormalisation time, and the norm chosen for the phasespace vectors. Its distribution only reveals the intricate dependence on the initial position of the system in the phase space. At 5 Gyr, the FTMLE roughly ranges from 0.15 to 0.5″ yr^{−1} with a 90% probability. Figure 1 shows that the truncated Hamiltonians ℋ_{2n} reproduce the asymptotic distribution of the FTMLE of the full dynamics ℋ, even at the lowest degree. It suggests, in particular, that the dynamical interaction of the Fourier harmonics at degree four constitutes the primary source of the observed FTMLE.
Fig. 1. Finitetime MLE of the forced secular ISS over 5 Gyr for different degrees of truncation of the Hamiltonian and the corresponding Lyapunov time FTMLE^{−1}. The bands represent the [5th, 95th] percentile range of the PDF estimated from ensembles of 128 (ℋ_{2n}) and 1148 (ℋ) stable orbital solutions with very close initial conditions. The black lines denote the nominal solution 𝒮 for nine different initial tangent vectors. 
4. Leading harmonics
Significant insight into the ISS dynamics is provided by the knowledge of the leading harmonics of the Hamiltonian, that is, those that drive the system trajectory in the action space the most, without necessarily being resonant. The action variables in turn control the essential longterm variation in the frequencies of motion (and therefore the activation of a specific web of resonances) through the mean frequencies that derive from the integrable part of the Hamiltonian in Eq. (5). The contribution of the harmonic k, ℓ to the action vector I(t) is
We ranked the 69 339 harmonics of ℍ_{6} according to the time median of the relative Euclidean norm δI_{k, ℓ}(t) = ‖ΔI_{k, ℓ}(t)‖/‖I(t)‖ along the nominal solution 𝒮 of Gauss’s dynamics, spanning 5 Gyr (Appendix D). We report the first 30 harmonics in Table 1. As usual, each harmonic is identified by the corresponding combination of frequency labels (Laskar 1990a). Given the absence of harmonics of order six, Table 1 confirms the leading role of those of order two and four, which enter the truncated Hamiltonian ℍ_{2n} at degree four. The harmonics θ_{1 : 1} and σ_{1 : 1} = (g_{1} − g_{5})−(s_{1} − s_{2}) appear among the very leading terms, confirming their dynamical relevance, as suggested in previous studies (Laskar 1990a, 1992; Lithwick & Wu 2011; Boué et al. 2012; Batygin et al. 2015). Among the top terms, there are harmonics that couple more extensively the DOFs of the inner planets, with the remarkable examples of (g_{1} − g_{4})+(s_{1} − s_{4}), associating the proper modes of Mercury and Mars and (s_{1} − s_{2})+(s_{3} − s_{4}), which concatenates all the inclination DOFs. A nonnegligible role of the Saturndominated eccentricity mode g_{6} also emerges. Even though these harmonics are not all necessarily resonant, they suggest that the dynamical entanglement of all the DOFs is significant along the nominal solution 𝒮. This consideration is supported by the study of partial Hamiltonians constructed from a limited number of Fourier harmonics (Appendix E). As Fig. E.1 shows, while the harmonics in Table 1 allow the asymptotic FTMLE of Fig. 1 to be robustly reproduced, simplified Hamiltonians based on the selection of specific DOFs may provide inconsistent predictions, with too low values down to a nonchaotic dynamics. This notably indicates that the resonant nature of a Fourier harmonic should only be established along the orbital solution of a realistic Hamiltonian.
Leading Fourier harmonics of ℍ_{6} along the nominal solution 𝒮 spanning 5 Gyr.
5. Resonant harmonics
Figure 1 and Table 1 suggest that the nonlinear interaction of the Fourier harmonics at degree four constitutes the primary source of chaos. In the framework of canonical perturbation theory, we constructed a Lie transform to define a change of variables that eliminates the terms of degree four from the truncated Hamiltonian ℍ_{2n} (Appendix F). In the Lietransformed Hamiltonian, , which is in Birkhoff normal form to degree four, these terms are replaced with the chain of highorder harmonics that arise from their dynamical interaction. Very importantly, the aim of this procedure is not to set up successive analytical approximations of the dynamics, which would be a vain goal. We simply define new canonical variables that let the interactions of the terms of degree four appear explicitly in the amplitudes of the Fourier harmonics at higher degrees.
To reveal the resonant harmonics of the Hamiltonian , we first retrieved those that present episodes of libration along the Lietransformed nominal solution, . Following Mogavero & Laskar (2021), we defined timedependent fundamental frequencies for the inner orbits (Appendix F). For each harmonic of , we then evaluated the corresponding combination of frequencies along the nominal solution. When this combination becomes null at least once over the time span of the solution, we consider the harmonic to be librating. This procedure filters out a majority of the Fourier harmonics as they never librate.
The appearance of libration episodes, and thus the existence of libration islands in the phase space, does not guarantee the resonant nature of a harmonic, as this is connected to the presence of both stable and unstable manifolds. Therefore, we employed the classic divideetimpera approach of Chirikov (1979) and considered a reduced Hamiltonian for each librating harmonic:
The resonant variables p, φ are related to the Lietransformed actionangle variables I′,θ′ by and φ = k ⋅ θ′ + Ωt, where we denote Ω = ℓ ⋅ ω_{o} (Appendix G). We therefore inspected the fictitious oneDOF dynamics that would be generated if k, ℓ were the only harmonic appearing in . The topology of the reduced phase space (p, φ) depends on seven integrals of motion, whose values relate to the position of the system in the fulldimensional action space. For a given librating harmonic, our study considers the oneparameter family of the reduced phase spaces that arise at each point along the nominal solution , that is, (Appendix G). This family of phase spaces is therefore spanned by the time t_{0}. As an example, Fig. 2 shows, for different t_{0} times along , the phase spaces of the harmonics θ_{1 : 1}, θ_{2 : 1}, and θ_{3 : 2}, with θ_{m : n} = m(g_{3} − g_{4})−n(s_{3} − s_{4}), as deduced from the Hamiltonian . The times were chosen across the first transition of θ_{2 : 1} from libration to circulation, which occurs at about 340 Myr along the nominal solution , and corresponds to the point where its FTMLE starts to increase in Fig. 1. The figure shows the level curves of the reduced Hamiltonian reconstructed by numerical integration as well as its fixed points, computed semianalytically with TRIP. For a resonant harmonic, separatrices emerge from the hyperbolic fixed points and enclose libration islands. The level curve corresponding to the value of the resonant variables at time t_{0} along is also shown and used to define the temporary libration or rotation state of the harmonic.
Fig. 2. Reduced phase spaces (p, φ) of the harmonics θ_{1 : 1}, θ_{2 : 1}, and θ_{3 : 2} along the Lietransformed nominal solution as deduced from the Lietransformed Hamiltonian . The phase spaces correspond to the t_{0} times indicated in the upperright boxes. The black dots reproduce the level curves of the reduced Hamiltonian, and the orange dots represent its fixed points. The separatrices are shown by red dots, and the level curve corresponding to the resonant variables at time t_{0} is shown in blue. The righthandside axes report the projection of the mean frequencies along the normal k to the resonance plane. 
Figure 2 shows how a resonant phase space can significantly differ from that of a simple pendulum, which is the universal model of nonlinear resonance and is often invoked to perform analytical computations. While the pendulum approximation turns out to be well suited for the majority of the resonant harmonics, it does not apply in general to the leading resonances, which are our main interest. We therefore discarded this approximation and performed systematic algebraic manipulations with TRIP to retrieve the fixed points of the reduced phase spaces through the roots of a univariate polynomial in the action p (Appendix G). Relying on a polynomial solver (Appendix G, Bini 1996; Bini & Fiorentino 2000; Bini & Robol 2014), we established the resonant nature of each harmonic along the nominal solution in a numerically robust and efficient way. We retrieved, in particular, the characteristic frequencies of the motion ω_{ell} and ω_{hyp} around the elliptic and hyperbolic fixed points, respectively. These frequencies are related to the square root of the eigenvalues of the variational matrix that governs the linear stability of the fixed points. Via a similar procedure, we retrieved the extrema of the action p along a level curve of the reduced Hamiltonian, which are employed to determine the librationrotation state of the harmonic and the resonance halfwidth, Δω. Following Chirikov (1979), we considered the projection of the frequency vector ω′=dθ′/dt along the normal to the resonance plane k ⋅ ω′+Ω = 0, that is, ω = ∥k∥^{−1}k ⋅ ω′. The reduced dynamics in Eq. (7) induces a variation in this projection equal to . Building on ideas behind frequency analysis (Laskar 1993), we then defined the resonance halfwidth Δω from the variation in the rotational frequency of the angle φ across the separatrix (Appendix G). In the pendulum approximation we have, up to a constant factor close to one,
which is the Chirikov (1979) halfwidth (recall that ω_{ell} = ω_{hyp} for the pendulum).
6. Origin of chaos
To establish a meaningful connection between the resonant phase spaces and the observed FTMLE, we assumed that the time distribution of physical observables along the nominal solution spanning 5 Gyr is representative of their ensemble distribution (that is, over a set of stable orbital solutions with very close initial conditions) at some large time of the order of billions of years. Based on this sort of finitetime ergodic assumption, we systematically retrieved for each librating harmonic of the samples of the fixed point frequencies ω_{hyp}, ω_{ell} and of the resonance halfwidth Δω along the nominal solution . Table 2 shows the first 30 resonant harmonics, along with their time statistic, as deduced from the truncation of the Lie transform at degree ten. We report for each frequency its time median as well as the 5th and 95th percentiles as subscript and superscript, respectively. The resonances are ranked according to their median halfwidth. To measure the overall dynamical impact of the terms, we denote as τ_{res} the fraction of time a harmonic is resonant and as τ_{libr} the fraction of libration states in the resonant case. Table 2 shows the harmonics that are resonant for more than 1% of the time.
First 30 resonant harmonics of along the nominal solution spanning 5 Gyr.
It is perhaps hopeless to obtain, for highdimensional dynamical systems, a precise formula connecting the FTMLE to the halfwidth of the leading resonances or to their fixed point frequencies. Nevertheless, studies of the periodically forced pendulum (Chirikov 1979; Holman & Murray 1996; Murray & Holman 1997; Li & Batygin 2014) suggest that these quantities should differ only by a factor of order unity from one another, that is,
To fix ideas, we also assumed that Eq. (8) holds beyond pendulum approximation. For the top resonances of Table 2, the statistical intervals of both Δω and ω_{hyp} extensively overlap the range spanned by the asymptotic FTMLE of the nominal solution in Fig. 1. Moreover, they are significantly compatible with the longterm ensemble distribution of the FTMLE. A statistical connection between the leading resonances of the Hamiltonian and the MLE thus emerges, even in the absence of relations more precise than Eqs. (8) and (9). It indicates the dynamical sources of the observed chaos. One may appreciate the low number of resonances with Δω ≳ 0.1″ yr^{−1}, which constitutes a sounding explanation for the large width of the FTMLE distribution, while the dense set of terms with Δω < 0.1″ yr^{−1} explains why the exponent never goes below this value. We point out the absence of resonances of order two. The harmonic g_{1} − g_{5}, in particular, is known to be involved in the very high eccentricities of Mercury (Boué et al. 2012) and is not expected to become resonant along a stable solution.
The resonances marked by a star in Table 2 only involve the proper modes 3 and 4, while those distinguished by a dagger exclusively concern the proper modes 1, 2, and 5. The first group of harmonics is composed of θ_{1 : 1}, θ_{2 : 1}, θ_{3 : 1}, and θ_{3 : 2} and is associated with the chaotic zone proposed in Laskar (1992). The second group belongs to the frequency subspace spanned by the combinations g_{1} − g_{5}, g_{2} − g_{5}, and s_{1} − s_{2} and includes the wellknown resonance σ_{1 : 1}. We first point out that the libration frequencies of 0.28 and 0.12″ yr^{−1} associated with θ_{2 : 1} and σ_{1 : 1} in Laskar (1990a) are consistent with the corresponding statistics of the elliptic fixed point frequency ω_{ell}. We then show that each of these two sets of harmonics constitutes a source of chaoticity. The emergence of chaos from resonant oneDOF phase spaces, as in Fig. 2, can be stated in two ways. As time changes and the system moves along its nominal trajectory, the phasespace region swept by the separatrices forms a chaotic zone. In an equivalent way, a chaotic zone results, at fixed time t_{0}, from separatrixsplitting when one restores the harmonics of that are suppressed in the reduced Hamiltonian. In any case, chaos derives from the interaction of resonant harmonics. Figure 3 represents such an interaction in terms of resonance overlap (Chirikov 1979) for the groups of starred and daggered harmonics. It shows the resonance planes in the frequency space along with asymmetric resonance layers, defined by the time median of the signed halfwidths Δω^{+}, Δω^{−} (Appendix G). The reported halfwidths are meaningful in the region visited by the nominal solution, which is also shown. The existence of the chaotic zone suggested in Laskar (1992) is indeed supported by a statistically robust overlap. The resonance θ_{2 : 1}, in particular, appears right in the centre of the resonance network, in agreement with its large τ_{res} value. Its amplitude results to a large extent from the interaction of θ_{1 : 1} and θ_{1 : 0} = g_{3} − g_{4} at degree four (see Table 1). The resonance overlap in the g_{1} − g_{5}, g_{2} − g_{5}, s_{1} − s_{2} subspace is shown in Fig. 3 in terms of integer combinations of σ_{1 : 1} and (g_{1} − g_{2})+(s_{1} − s_{2}). It appears to be more relevant, at least for our nominal solution, than the restriction to the subspace spanned by g_{1} − g_{5} and s_{1} − s_{2} investigated in previous studies (Lithwick & Wu 2011; Batygin et al. 2015). A relevant dynamical role of the proper mode g_{2}, hinted at in Lithwick & Wu (2011) and Batygin et al. (2015), is confirmed. Figure 3 also suggests a longterm diffusion mainly perpendicular to the resonance planes, while Arnold’s diffusion along them seems negligible (Laskar 1993), at least for the present regular (i.e. noncollisional) solution.
Fig. 3. Overlaps of secular resonances from Table 2. Left panel: starred resonances θ_{m : n} = m(g_{3} − g_{4}) − n(s_{3} − s_{4}) in the frequency subspace spanned by g_{4} − g_{3} and s_{4} − s_{3}. Right panel: daggered resonances as integer combinations of (g_{1} − g_{5})−(s_{1} − s_{2}) and (g_{1} − g_{2})+(s_{1} − s_{2}) (the harmonic ℱ_{66} = 3g_{1} − g_{2} − 2g_{5} − (s_{1} − s_{2}) is not shown in Table 2). The dashed lines represent the resonance centres, and the coloured regions correspond to the time median of the signed halfwidths. The nominal solution , which spans 5 Gyr, is shown after lowpass filtering with a cutoff frequency of (5 Myr)^{−1}, the cross indicating the current position of the ISS. Frequencies are in arcsec yr^{−1}. 
Right after the top resonances, harmonics with Δω < 0.1″ yr^{−1} significantly couple all the proper modes, confirming the implications of Table 1. Coupling resonances such as (g_{2} − g_{4})+(s_{2} − s_{4}), (g_{1} − g_{3})+(s_{2} − s_{3}) and (g_{1} − g_{4})+(s_{1} − s_{4}), along with highorder ones resulting from their interaction with leading harmonics, strongly suggest that resonance overlap takes place extensively in the fulldimensional frequency space. Even though we have isolated lowdimensional sources of chaoticity, the largescale chaos in the ISS is probably best understood as a highdimensional phenomenon.
The statistical nature of these findings suggests that the harmonics of Table 2 should be found statistically in resonance over an ensemble of sufficiently long orbital solutions of the Solar System, independently of the precise proper modes and dynamical model employed. A vivid example is provided by the term (g_{1} − g_{2})+(s_{1} − s_{2}) that we found among the top resonances: it enters libration at around 1 Gyr in Lithwick & Wu (2011) but was not mentioned for the much shorter solutions of previous studies (Laskar 1990a, 1992; Sussman & Wisdom 1992; Laskar et al. 2004). More generally, we show in Fig. 4 examples of libration episodes for a number of harmonics in Table 2 along the direct integrations of the Solar System presented in Laskar & Gastineau (2009) and spanning 5 Gyr, when one uses the proper modes z^{•}, ζ^{•} defined in Laskar (1990a). Despite the different dynamical model and proper modes, one systematically confirms the librations predicted by the theory by inspecting just a few solutions. We point out the remarkable case of (g_{1} − g_{3})+(s_{2} − s_{3}), which will surely be librating over the next few tens of millions of years. Clearly, the chaotic zone of the ISS may extend to resonances that are barely, or not at all, active along our nominal solution and thus do not appear in Table 2. In any case, Fig. 4 suggests that, despite being deduced from the forced secular dynamics, the web of resonances pictured here should underlie the chaotic dynamics of the inner planets in every realistic model of the Solar System.
Fig. 4. Combination of angle variables for different harmonics of Table 2 as a function of time along direct integrations of the Solar System presented in Laskar & Gastineau (2009) (n = 0 is the nominal solution). The angle variables correspond to the proper modes z^{•}, ζ^{•} defined in Laskar (1990a). The time series are lowpassfiltered with a cutoff frequency of (500 kyr)^{−1}. 
The polynomial F_{1} generally has complex coefficients as the phases of the quasiperiodic decomposition of the giant planet orbits in Eq. (3) are explicitly contained in the complex amplitudes.
In the general case, the integral must be split over subsets of [0, 2π] for which does not change sign, and the contributions summed in absolute value to obtain the period of the motion. We simply do not compute the rotational frequency when φ is not monotonic, as this constitutes a minority of cases.
The harmonic θ_{1 : 1} is resonant for only 0.3% of the time at degree six, so the corresponding PDFs present spikes. We thus had to adapt the bandwidth choice to get a smoother estimation for Δω. We also point out that θ_{1 : 1} possesses no negative signed halfwidth at degree ten (see its phase space at 450 Myr in Fig. 2).
Acknowledgments
We are truly indebted to M. Gastineau for his support with TRIP. We acknowledge the precious help of L. Robol for the integration of MPSolve in TRIP, and the fruitful discussions with N. H. Hoàng and M. Saillenfest. F.M. has been supported by a PSL postdoctoral fellowship and by a grant of the French Agence Nationale de la Recherche (AstroMeso ANR19CE31000201). This project has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo885250). This work was granted access to the HPC resources of MesoPSL financed by the Region ÎledeFrance and the project Equip@Meso (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
References
 Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120 [CrossRef] [Google Scholar]
 Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J. M. 1980, Meccanica, 15, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D. 1996, Numer. Algorithms, 13, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D., & Fiorentino, G. 2000, Numer. Algorithms, 23, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Bini, D. A., & Robol, L. 2014, J. Comput. Appl. Math., 272, 276 [CrossRef] [Google Scholar]
 Boué, G., Laskar, J., & Farago, F. 2012, A&A, 548, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chirikov, B. V. 1960, J. Nucl. Energy, 1, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Chirikov, B. V. 1969, Research Concerning the Theory of Nonlinear Resonance and Stochasticity, (Novosibirsk: Institute of Nuclear Physics) [in Russian], CERNTrans7140 (1971) [Google Scholar]
 Chirikov, B. V. 1979, Phys. Rep., 52, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Deprit, A. 1969, Celest. Mech., 1, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Deprit, A., Henrard, J., & Rom, A. 1970, Science, 168, 1569 [NASA ADS] [CrossRef] [Google Scholar]
 Froeschlé, C., Gonczi, R., & Lega, E. 1997, Planet. Space Sci., 45, 881 [CrossRef] [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194 [Google Scholar]
 Gastineau, M., & Laskar, J. 2021, TRIP 1.4.120, TRIP Reference Manual, IMCCE, Paris Observatory, https://www.imcce.fr/trip [Google Scholar]
 Gauss, C. 1818, Werke, 3, 331 [Google Scholar]
 Hairer, E., Nørsett, S., & Wanner, G. 1993, Solving Ordinary Differential Equations I – Nonstiff Problems, 2nd edn. (Berlin, Heidelberg: SpringerVerlag) [Google Scholar]
 Hayes, W. B. 2007, Nat. Phys., 3, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Henon, M., & Heiles, C. 1964, AJ, 69, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Hoang, N. H., Mogavero, F., & Laskar, J. 2021, A&A, 654, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Holman, M. J., & Murray, N. W. 1996, AJ, 112, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, G. 1966, PASJ, 18, 287 [NASA ADS] [Google Scholar]
 Laskar, J. 1984, PhD Thesis, Observatoire de Paris, France [Google Scholar]
 Laskar, J. 1985, A&A, 144, 133 [NASA ADS] [Google Scholar]
 Laskar, J. 1988, A&A, 198, 341 [NASA ADS] [Google Scholar]
 Laskar, J. 1989, Nature, 338, 237 [Google Scholar]
 Laskar, J. 1990a, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1990b, in Modern Methods in Celestial Mechanics, eds. D. Benest, & C. Froeschle, 89 [Google Scholar]
 Laskar, J. 1991, NATO Adv. Study Inst. (ASI) Ser. B, 272, 93 [NASA ADS] [Google Scholar]
 Laskar, J. 1992, in Chaos, Resonance, and Collective Dynamical Phenomena in the Solar System, ed. S. FerrazMello, IAU Symp., 152, 1 [NASA ADS] [Google Scholar]
 Laskar, J. 1993, Phys. D Nonlinear Phenom., 67, 257 [Google Scholar]
 Laskar, J. 2005, in Hamiltonian Systems and Fourier Analysis: New Prospects for Gravitational Dynamics, eds. D. Benest, C. Froeschlé, & E. Lega (Cambridge Scientific Publishers Ltd), 93 [Google Scholar]
 Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lecar, M., Franklin, F. A., Holman, M. J., & Murray, N. J. 2001, ARA&A, 39, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Li, G., & Batygin, K. 2014, ApJ, 790, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Wu, Y. 2011, ApJ, 739, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Locatelli, U., & Giorgilli, A. 2000, Celest. Mech. Dyn. Astron., 78, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, C., Meyers, S. R., & Sageman, B. B. 2017, Nature, 542, 468 [NASA ADS] [CrossRef] [Google Scholar]
 Mogavero, F., & Laskar, J. 2021, A&A, 655, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (Taylor & Francis) [Google Scholar]
 Murray, N., & Holman, M. 1997, AJ, 114, 1246 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, N., & Holman, M. 2001, Nature, 410, 773 [NASA ADS] [CrossRef] [Google Scholar]
 Olsen, P. E., Laskar, J., Kent, D. V., et al. 2019, Proc. Natl. Acad. Sci., 116, 10664 [NASA ADS] [CrossRef] [Google Scholar]
 Oseledec, V. I. 1968, Trans. Moscow Math. Soc., 19, 197 [Google Scholar]
 Parzen, E. 1962, Ann. Math. Stat., 33, 1065 [Google Scholar]
 Poincaré, H. 1896, Comptes rendus hebdomadaires de l’Académie des sciences de Paris, 123, 1031 [Google Scholar]
 Prince, P., & Dormand, J. 1981, J. Comput. Appl. Math., 7, 67 [CrossRef] [Google Scholar]
 Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenblatt, M. 1956, Ann. Math. Stat., 27, 832 [Google Scholar]
 Silverman, B. W. 1986, Density Estimation for Statistics and Data Analysis (London: Chapman & Hall/CRC) [Google Scholar]
 Sussman, G. J., & Wisdom, J. 1992, Science, 257, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, G. H., & Ford, J. 1969, Phys. Rev., 188, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J. 1980, AJ, 85, 1122 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Zeebe, R. E., & Lourens, L. J. 2019, Science, 365, 926 [NASA ADS] [CrossRef] [Google Scholar]
 Zurbenko, I. G., & Smith, D. 2018, WIREs Comput. Stat., 10, 1419 [Google Scholar]
Appendix A: The need for computer algebra. TRIP
When it comes to elucidating the causes of a chaotic dynamics in celestial mechanics, one typically resorts to the analysis of the Fourier harmonics of a given Hamiltonian (Chirikov 1960, 1969; Walker & Ford 1969; Wisdom 1980). The problem is first reduced to a set of partial oneDOF Hamiltonians that describe the dynamics arising from each harmonic independently of one another. The interaction of pairs of harmonics is then addressed through the resonance overlap criterion (Chirikov 1960; Walker & Ford 1969; Chirikov 1979) and Poincaré’s sections (Henon & Heiles 1964) to reveal the dynamical sources of chaos. In practice, the choice of the harmonics to study may rely on a somewhat intuitive perception of their relevance, supported by observational data, numerical investigations, and the findings of previous research.
When the Hamiltonian system under investigation has a few DOFs, the above approach is quite straightforward to apply and very reliable (e.g. Wisdom et al. 1984). Conversely, when dealing with several DOFs, as in case of the Solar System planets, the curse of dimensionality makes its implementation especially difficult. First of all, the number of loworder harmonics that are potentially resonant rapidly grows with the dimension of the phase space: an a priori choice of the pertinent harmonics may thus lead to a biased analysis. More importantly, the resonant nature of each harmonic, even in a oneDOF model, still depends on the system position in the fulldimensional phase space. Some authors have invoked the quasiperiodicity of certain variables to reduce the analysis to a lowerdimension phase space (Lithwick & Wu 2011; Boué et al. 2012; Batygin et al. 2015). However, the validity of such an assumption strongly depends on the problem under investigation, and it cannot provide the basis of a general approach. In the ISS, for instance, the fundamental frequencies of the planet orbits all vary in an intricate way in a few tens of millions of years, with the exception of the eccentricity mode g_{2}, which has slower frequency variations (Laskar 1990a; Laskar et al. 2004; Hoang et al. 2021). A quasiperiodic approximation for the corresponding DOFs is thus not consistent over timescales much longer than 10 Myr (Mogavero & Laskar 2021).
To overcome the difficulties associated with the highdimensional phase space of the ISS, we set up an unbiased systematic analysis of the Fourier harmonics of the secular planetary Hamiltonian. To this end, we employed TRIP, a computer algebra system dedicated to the perturbation series of celestial mechanics, which has been developed over the past 30 years at the Institut de mécanique céleste et de calcul des éphémérides (Laskar 1990b; Gastineau & Laskar 2011, 2021). The main objects of its symbolic kernel are the Poisson series, that is, multivariate Fourier series whose coefficients are multivariate Laurent series^{1},
where and are complex and real variables, respectively, and . The angle variables enter the series through their complex exponential, which is encoded in TRIP as a dedicated variable. The complex coefficients C_{k, ℓ} can be rational functions or numerical coefficients of different types: double, quadruple, and multipleprecision floatingpoint numbers; fixed or multipleprecision integers or rational numbers; or double and quadrupleprecision floatingpoint intervals. The symbolic kernel implements several operations on the Poisson series: sum, product, division, differentiation, integration, exponentiation, substitution of variables, and selection of specific terms. At the heart of these functionalities, TRIP embeds a monomial truncated product, which allows the product of two series to be truncated according to a given condition on the partial or total degree of one or more variables (truncation based on the size of the coefficients C_{k, ℓ} is also available). Truncations are formal objects in TRIP and permit perturbative computations to be performed at a minimal cost. TRIP also provides special functionalities for celestial mechanics, such as Poisson brackets and the manipulation of formal Lie series, which allow the algorithms of canonical perturbation theory (Hori 1966; Deprit 1969) to be implemented. A numerical kernel, based on vectors, matrices, and multidimensional tables, provides a stateoftheart evaluation of the Poisson series. It also includes a number of routines typical of more general computer algebra systems, such as algorithms for the zeros of univariate polynomials, the solutions of systems of algebraic equations, and the analysis of time series (e.g. fast Fourier transform, frequency analysis, integration, and interpolation). Moreover, TRIP natively integrates polynomial dynamical systems through Adams PECE and DOPRI8 methods (Hairer et al. 1993; Prince & Dormand 1981). Finally, TRIP provides a dedicated procedural language that supports loops, conditional statements, and function and structure definition, which permits the symbolic and numerical kernels to be interfaced in autonomous programs. When unsupported functionalities are needed, it allows external C and Fortran routines to be called and allows communication with other computer algebra systems, such as Maple and Mathematica.
Appendix B: Development of the secular planetary Hamiltonian
The secular Hamiltonian of the entire Solar System in Eq. (1) can be expanded in series of the Poincaré complex variables and truncated at a given total degree 2n, where n is a positive integer (Laskar 1991; Laskar & Robutel 1995). This results in a polynomial Hamiltonian , where groups all the monomials of same total degree 2p in the variables . The expansion straightforwardly provides the truncated Hamiltonian for the forced ISS,
At the lowest degree, the ℋ_{2} Hamiltonian describes an integrable, forced LL dynamics for the inner planets. The corresponding equations of motion are analytically solved by introducing complex variables , with conjugate momentumcoordinate pairs, that correspond to the proper (i.e. normal) modes of the free oscillations (Mogavero & Laskar 2021). By switching to the new set of canonical variables (u_{k}, v_{k}), the truncated Hamiltonian in Eq. (B.1) transforms to
where ℍ_{(2p)} groups all the terms of same total degree 2p (Mogavero & Laskar 2021). Actionangle variables that trivially integrate the LL problem are thus introduced as
for k ∈ {1, 2, 3, 4} and allow the truncated Hamiltonian in Eq. (B.2) to be written as a finite Fourier series,
where are complex amplitudes, with , and where we employ a compact notation for the actionangle variables, that is, I = (X, Ψ) and θ = (χ, ψ). At the lowest degree, one has the LL Hamiltonian , with ω_{LL} = −(g_{LL}, s_{LL}), and the Hamiltonian in Eq. (B.4) is therefore in quasiintegrable form. In the framework of canonical perturbation theory, the explicit time dependence of the Hamiltonian (coming from the forcing of the outer planets) can be absorbed in a phasespace extension, resulting in
The dynamics of the additional angles is consistently given by , while that of the conjugated actions Φ is irrelevant.
Appendix C: Nominal solution and choice of initial conditions
We estimated the statistical properties of the chaotic dynamics in the ISS from the orbital solution 𝒮 of Gauss’s dynamics (Hamiltonian ℋ in Eq. (4)) that is numerically integrated over 5 Gyr starting from the nominal initial conditions given in Mogavero & Laskar (2021, Appendix D). This solution is denoted as the nominal solution throughout the Letter and consists of the values of the canonical variables sampled with a timestep Δt of 1000 yr, that is, with 𝒮(t_{m}) = (I(t_{m}),θ(t_{m})) and t_{m} = mΔt. We considered the time distribution of a physical observable along the nominal solution 𝒮 as representative of its ensemble distribution (that is, over a set of stable orbital solutions with very close initial conditions) at some large time of the order of billions of years.
Following Mogavero & Laskar (2021), we chose the initial conditions for the ensemble computation of the FTMLE by taking the relative variation in each coordinate of the nominal phasespace vector of Gauss’s dynamics as a normal random variable, with zero mean and a standard deviation of 10^{−9}. In other words, the initial conditions are distributed according to
for k ∈ {1, 2, 3, 4}, where represents the initial conditions of the nominal solution 𝒮, are standard normal deviates, and σ = 10^{−9}. An analogous expression holds for the variables y_{i}. The initial conditions thus follow a multivariate Gaussian distribution centred at the nominal initial conditions of ℋ, with a diagonal covariance matrix. The initial tangent vectors that we used in the application of the Benettin et al. (1980) algorithm are also sampled from a multivariate Gaussian distribution, with null mean and an identity covariance matrix. We chose a renormalisation time of 5 Myr and the Euclidean norm for the phasespace vectors.
Appendix D: Leading harmonics
There can be many different ways of ranking the Fourier harmonics of Eq. (B.4), depending on the dynamical quantities that one aims to track. As we intend to reconstruct the resonant structure of the ISS, we are interested in the longterm variation in the frequencies of motion . From Eq. (B.4), Hamilton’s equations give
where (k, ℓ)∈ℤ^{8} × ℤ^{7} and the star means that the null wave vector (0, 0) is excluded from the summation. As the Hamiltonian ℍ_{2n} is close to integrable, the actions I(t) vary over a timescale much longer than that of the angles θ(t). Therefore, Eq. (D.1) suggests that the frequencies ω_{2n}(t) rapidly fluctuate around the slowvarying mean frequencies , only depending on the action variables. If the motion were quasiperiodic, the harmonics appearing in the summation of Eq. (D.1) could be systematically eliminated by a series of canonical transformations in the form of Lie transforms. The mean frequencies would then give the (constant) frequencies of motion as a function of the transformed action variables. For a chaotic dynamics, the presence of resonant harmonics shifts the shorttime average of ω_{2n} with respect to . Nevertheless, the mean frequencies still provide the overall longterm behaviour of the frequencies of motion. We thus chose to rank the Fourier harmonics of the Hamiltonian according to their contribution to the system trajectory in the action space, which in turn drives the variation in the mean frequencies with time.
Hamilton’s equation gives the contribution of the harmonic k, ℓ to the action vector I(t) as
By construction, the sum of the contributions from all the harmonics exactly reconstructs the system displacement in the action space at a given time, that is, ΔI(t) = I(t)−I(0). The integral was computed numerically in TRIP for each harmonic along the sampled nominal solution 𝒮. The harmonic ranking was established by considering the relative Euclidean norm of the contribution,
The set of values is treated as a sample from an underlying PDF and characterised via statistical estimators. The harmonics in Table 1 are listed according to their median contribution, while the 5th and 95th percentiles show the dispersion of δI_{k, ℓ}.
Appendix E: Partial Hamiltonians
The identification of the leading harmonics of the Hamiltonian allows the FTMLE to be reproduced from a small set 𝒰 of Fourier harmonics, by considering
Partial Hamiltonians can be systematically constructed in TRIP through the selection of specific terms of ℍ_{6}. We integrated the corresponding variational equations to obtain the FTMLE distribution of an ensemble of 128 stable solutions with initial conditions very close to the nominal ones of Gauss’s dynamics. Figure E.1 shows the [5th, 95th] percentile range of the FTMLE distribution from the 30 leading harmonics of Table 1. This simplified dynamics reproduces the asymptotic FTMLE distribution of the Hamiltonians ℋ and ℋ_{2n} very well, confirming the dynamical relevance of the leading harmonics. Very importantly, the FTMLE prediction turns out to be robust with respect to the addition of further harmonics from our ranking.
Fig. E.1. Finitetime MLE and corresponding Lyapunov time FTMLE^{−1} from different partial Hamiltonians (Eq. (E.1)). The bands represent the [5th, 95th] percentile range of the PDF estimated from an ensemble of 128 stable orbital solutions with very close initial conditions. The blue colour corresponds to the Fourier harmonics of Table 1, the green colour to the harmonics of the same table that only include the angular DOFs 1, 2, and 5, and the red colour to those harmonics that only involve the angular DOFs 3 and 4. The hatched regions correspond to the case where the choice of the DOFs is made across the entire ℍ_{6} Hamiltonian. The grey region stands for the full ℋ Hamiltonian, as shown in Fig. 1. 
The construction of partial Hamiltonians allows the extensive entanglement of the DOFs already shown in Table 1 to be highlighted. Figure E.1 reports the FTMLE distribution from the leading harmonics of Table 1 that only involve the angles of the proper modes 1, 2, and 5, that is, the harmonics 𝒰 = {ℱ_{6}, ℱ_{8}, ℱ_{20}, ℱ_{27}, ℱ_{29}, ℱ_{30}} only involving combinations of g_{1}, g_{2}, g_{5}, s_{1}, and s_{2}. This partial Hamiltonian includes, in particular, the Fourier harmonics considered in the simplified longterm dynamics of Mercury in Batygin et al. (2015). The predicted Lyapunov time of about 100 Myr is clearly incompatible with the full system dynamics. The prediction does not improve when the selection of the harmonics only involving g_{1}, g_{2}, g_{5}, s_{1}, and s_{2} is made across the entire ℍ_{6} Hamiltonian (as shown by the hatched green region). The harmonics of Table 1 that only involve the angles of the proper modes 3 and 4, that is, the harmonics 𝒰 = {ℱ_{2}, ℱ_{3}, ℱ_{5}, ℱ_{25}} only involving g_{3}, g_{4}, s_{3} and s_{4}, produce an asymptotic FTMLE that is compatible with the full system dynamics (as shown by the red region in Fig. E.1). Nevertheless, when the selection of such harmonics is made across the entire ℍ_{6} Hamiltonian, the resulting FTMLE turns out to be a monotonically decreasing function of time, suggesting an integrable dynamics (hatched red region in Fig. E.1). These numerical experiments show that freezing a priori some DOFs may dramatically change the behaviour of the system, resulting in predictions that may largely depend on the specific choice of the Fourier harmonics. Clearly, the surprising behaviours of Fig. E.1 may depend on the fact that the initial conditions of the simplified dynamics are not adjusted according to the choice of the harmonics. Proper initial conditions should be computed from the Lie transform that allows the harmonics of the Hamiltonian ℍ_{6} that do not appear in to be eliminated. However, this adjustment is often omitted when simplified models are constructed in literature. Due to this fact, we decided to keep the nominal initial conditions of Gauss’s dynamics in these numerical experiments. In any case, taking the ensemble of the DOFs of the inner system into account seems essential. The resonant nature of a Fourier harmonic, in particular, should only be established along the orbital solution of a realistic Hamiltonian.
Appendix F: Lie transform
We employed canonical perturbation theory to construct a change of variables that eliminates the Fourier harmonics of degree 4 from the truncated Hamiltonian, that is, the harmonics with a nonnull wave vector appearing in ℍ_{(4)} in Eq. (B.2). Such a transformation of variables can be canonically defined as the time1 flow of a generating Hamiltonian S that satisfies the homologic equation
where the braces represent the Poisson bracket and where we employ the extended phasespace formalism of Eq. (B.5) to deal with a timeindependent Hamiltonian. The resulting change of variables is given by the formal Lie transform
where L_{S} ⋅ ={S, ⋅} is the Lie derivative associated with the generating function S. The exponential of an operator A acting on the phasespace functions is formally expressed as
with A^{0} defined as the identity operator and A^{q} = AA^{q − 1} for q ≥ 1. The transformed Hamiltonian reads
where the Lie transform is formally evaluated at the new canonical variables I′,θ′;Φ′,ϕ′. The solution to the homologic Eq. (F.1) is expressed as
As ℍ_{(4)} possesses a finite number of nonnull harmonics, the Fourier series in Eq. (F.5) is finite. Moreover, all the denominators being constant and different from zero, the generating function S is a welldefined analytical function.
Because the generating function S in Eq. (F.5) does not depend on the actions Φ, the only term appearing in ℍ′ that depends on the transformed variables Φ′ is . The values of Φ′ are therefore irrelevant from a dynamical point of view (as are those of Φ), and the transformation Φ′=E^{−LS}Φ can be discarded. Moreover, from {S, ϕ}=0 it follows that ϕ′=ϕ, in agreement with the previous point.
To implement the exponential operator in TRIP, the formally infinite summation in Eq. (F.3) must be truncated. As all the terms in the generating function S are of degree 4, the Lie derivative L_{S} raises the degree of the Hamiltonian terms by 2. Therefore, to truncate the transformed Hamiltonian at degree 2n, it is sufficient to compute the summation up to the order n − 1 and then discard the terms of degree higher than 2n, that is,
Once the Hamiltonian is computed, the corresponding nominal solution is obtained from the truncation of the operator E^{−LS} at the same order n − 1, that is, . Computing the change in the nominal solution is essential for consistently evaluating the amplitudes of the Hamiltonian terms resulting from the Lie transform.
F.1. Convergence
Since the generating Hamiltonian S is an analytical function, the convergence of the series generated by the Lie transforms E^{±LS} depends on the size of its Fourier coefficients (Morbidelli 2002). To address the convergence of our Lie transform, we first plotted the relative increments of the transformed nominal solution as a function of time and truncation degree (see Fig. F.1). We used the Euclidean norm for the phasespace vectors, as usual. To highlight the longterm trend of the time series, which may be hidden by large shorttime oscillations, we applied the lowpass KolmogorovZurbenko (KZ) filter with three iterations of the moving average and a cutoff frequency of (5 Myr)^{−1} (Zurbenko & Smith 2018; Mogavero & Laskar 2021). As also shown by the time mean of the increments , the contributions to decay exponentially with the degree of truncation, up to degree ten at least. That said, the decrease is in practice quite slow. Figure F.2 then addresses the convergence of the Lietransformed Hamiltonian along the transformed nominal solution . The relative increments are plotted as a function of time and truncation degree 2n. We also report the convergence of the original ℍ_{2n} Hamiltonian along the nominal solution 𝒮 for comparison. Again, the contributions to decrease with the degree of truncation, up to degree ten at least, even though the decay is slow when compared to that of the original ℍ_{2n} Hamiltonian. Even if the asymptotic convergence cannot be guaranteed, the absence of any divergence in Figs. F.1 and F.2 suggests the stability of the estimations based on the truncation of the Lie transform at a finite degree 2n ≤ 10. Some oscillations of these estimations can be expected as a result of the slow (finitedegree) convergence of the series.
Fig. F.1. Convergence of the Lietransformed nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are lowpassfiltered with a cutoff frequency of (5 Myr)^{−1}. 
Fig. F.2. Convergence of the Lietransformed Hamiltonian along the nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are lowpassfiltered with a cutoff frequency of (5 Myr)^{−1}. The convergence of the original ℍ_{2n} Hamiltonian along the nominal solution 𝒮 is also shown. 
F.2. Librating harmonics
To retrieve the set of librating harmonics, we defined timedependent fundamental frequencies along a sampled orbital solution, as we did in Mogavero & Laskar (2021) in the case of the frequency g_{1}. Following Eq. (B.3), the instantaneous (angular) frequencies of the proper modes (u_{i}(t),v_{i}(t)) are given by . We sampled the time derivatives via the Lagrange fivepoint formula, as very high numerical precision is unnecessary. Since the instantaneous frequencies present large shorttime fluctuations that hide their longterm evolution, we applied the lowpass KZ filter (Zurbenko & Smith 2018; Mogavero & Laskar 2021) to the time series and defined . We used three iterations of the moving average and a cutoff frequency of (5 Myr)^{−1}. This specific value was chosen to filter out the nonchaotic part of the dynamical spectrum from the time series. When the fundamental frequencies were computed for a given orbital solution, the librating harmonics were defined from the combinations of frequencies that become null at least once over its time span. We find around 1000, 15 000, and 150 000 librating harmonics among those appearing in along the Lietransformed nominal solution at degree six, eight, and ten, respectively.
Appendix G: Reduced Hamiltonians
To retrieve the resonant harmonics of the Lietransformed Hamiltonian , that is,
we considered the dynamics generated by each Fourier component separately. If we assume that at time t = t_{0} along the nominal solution the system is at position in the phase space, the dynamics that would be generated by the harmonic k, ℓ, starting at , if all the other harmonics are ignored, is given by the reduced Hamiltonian
where τ ≥ t_{0} is the time of the fictitious dynamics, Ω = ℓ ⋅ ω_{o}, and we consider k ∈ ℤ^{8} ∖ {0}. The resonant nature of the harmonic depends on the closeness of the actions to the resonance condition
which defines a hyperplane of normal vector k in the space of the mean frequencies and an algebraic variety in the space of the actions I′ (at degree four this is a hyperplane as depends linearly on I′). According to the reduced Hamiltonian, the actions I′ oscillate along the direction k, that is,
Following the Chirikov (1979) derivations, we performed the explicit reduction of the Hamiltonian in Eq. (G.2) to one DOF. We considered the canonical transformation (I′,θ′) → (p, φ) defined by the timedependent generating function
where 𝕄 is a real 8 × 8 matrix, , and T denotes the transposition operator. One has
We chose the first row of the matrix 𝕄 to be k^{T}, and the remaining seven rows to be orthonormal vectors perpendicular to k, that is, and for 2 ≤ i, ȷ ≤ 8. It follows that
with p = (p_{1}, …, p_{8})^{T}. We finally chose ν = (Ω, 0, …, 0)^{T} so that φ_{1} = k ⋅ θ′+Ωτ. The reduced Hamiltonian governing the new canonical variables (p, φ) is given by
The reduced dynamics only affects the momentum p_{1}, and are all integrals of motion. These constants can be set to zero by choosing and employing the initial condition so that p_{1}(τ = t_{0}) = 0 and
As the action variables I′ are nonnegative by definition, the dynamics of the momentum p_{1} is generally restricted to a subset of the real line that is bounded from above, below, or both, depending on the wave vector k.
G.1. Reduced phase spaces
A different choice considered by Chirikov (1979) consists in taking on the resonance surface given by Eq. (G.3), with the additional constraint that is parallel to the direction k, that is,
With this choice, the integrals of motion have null values and the action variables are given by
To simplify the notation, we omit the subscripts of the resonant variables (p_{1}, φ_{1}) from now on. With Chirikov’s choice and under the assumption of sufficiently small oscillations of the momentum p, the reduced Hamiltonian can be expanded around and provides, at first order, the universal description of a nonlinear resonance in terms of pendulum dynamics (Chirikov 1979). In this case the vector represents the centre of oscillation in the action space.
When truncating the Lie transform in Eq. (F.4) at degree 2n, the two constraints in Eq. (G.10) can be rewritten as a polynomial equation of degree n − 1 in the variable p(t_{0}), which thus possesses up to n − 1 real solutions. As a generalisation of the pendulum approximation, when no real solutions for p(t_{0}) exist, or when the corresponding vectors have some negative components, generally no hyperbolic fixed points appear in the reduced phase space (p, φ). The dynamics is thus globally characterised by rotation states of the angle φ. Depending on the Fourier amplitude , libration islands may still exist, for example close to the boundaries of the p variable. However, in the absence of hyperbolic points, they are not enclosed by a separatrix. These reduced phase spaces are therefore considered as nonresonant. When a real solution for p(t_{0}) exists and the corresponding point lies in the physical action space (that is, its components are all nonnegative), hyperbolic fixed points generally appear in the phase space. The separatrices that emerge from these points separate resonant libration states of φ from rotation states. The position of the fixed points along the momentum axis is offset with respect to p = 0 (corresponding to the action point ) by an amount that depends on the Fourier coefficient .
G.2. Fixed points
The topology of a reduced phase space depends on the position of the system in the eightdimensional action space, that is, it depends on through Eq. (G.9) or Eq. (G.10). An explicit study of the oneDOF Hamiltonian in Eq. (G.8) can be carried out when a given point is considered. The fundamental idea in this work is to retrieve the topology of the reduced phase spaces along the nominal solution , that is, . For each Fourier harmonic k, ℓ, one obtains a oneparameter family of phase spaces spanned by the time t_{0}. We studied their topology through the retrieval of the fixed points. Their position can be analytically obtained in the pendulum approximation; however, while the pendulum model works fine for the majority of the harmonics, it turns out not to be pertinent for the leading resonant ones, which are our main interest. Fortunately, as the following derivations show, computer algebra allows the fixed points of a general reduced phase space to be systematically retrieved, making the pendulum approximation unnecessary.
In light of Eqs. (G.9) and (G.11), the reduced Hamiltonian can be written as
where we rename , , and as ℏ, f_{0}, and f_{1}, respectively, to keep a simpler notation. From Eq. (B.3) it follows that one can write , with
where we denote and . The function μ is thus a real multivariate monomial of degree in the square roots of the action variables, while F_{1} can be expressed as a multivariate polynomial of degree (2n − k−ℓ)/2 in the action variables^{2}, the order of the harmonic k, ℓ (i.e. k+ℓ) being an even integer. The fixed points of the Hamiltonian ℏ satisfy the system of equations
where d stands for derivation with respect to p. The function df_{0} is a multivariate polynomial of degree n − 1 in the action variables and thus a univariate polynomial of the same degree when expressed in the variable p. Because of the square roots appearing in f_{1} and in its derivative df_{1}, the system of equations (G.14) is nonalgebraic. Nevertheless, systematic manipulations of these equations allow its solutions to be found through the roots of a univariate polynomial in the variable p. Searching for solutions that correspond to nonnull values of the action variables I′, we assumed μ ≠ 0 and write the first equation as . From Eq. (G.13) it follows that the square root dependence in df_{1} can be eliminated through multiplication by the function
that is, the function can be expressed as a polynomial of degree in the action variables. Through multiplication by η and taking the square, the second equation of system (G.14) implies . The first equation then allows the dependence on the variable φ to be eliminated from the second equation, providing the system
The second equation is now a univariatepolynomial equation in the variable p of degree . Its complex solutions can be systematically retrieved through a univariate polynomial solver. Among all the solutions, one has to keep the real ones that correspond to positive values of the action variables I′ in Eqs. (G.9) or (G.11). For each solution p^{⋆}, the value of the angle φ^{⋆} is then given by , where the sign is chosen to restore the sign loss that occurs in squaring the second equation of the system (G.14).
The linear stability of a fixed point is assessed from the sign of the eigenvalues of the variational matrix
evaluated at the fixed point. The matrix 𝕍 is the product of the symplectic matrix 𝕁 = ((0, −1),(1, 0)) with the Hessian of the reduced Hamiltonian, and its eigenvalues λ are given by the equation λ^{2} = −det(𝕍). In the case of an elliptic fixed point (i.e. when λ^{2} < 0), the angular frequency of the small oscillations is given by , while in the case of a hyperbolic fixed point (i.e. when λ^{2} > 0) we denote .
G.3. Level curves
We systematically characterised the level curves of the reduced Hamiltonian by retrieving their extrema. The tangent direction t to the level curve ℏ(p, φ) = ℏ_{0} can be expressed as
The extrema of the variable p along the level curve satisfy the equation ∂ℏ/∂φ = 0, while the extrema of φ are given by ∂ℏ/∂p = 0. To distinguish between the minima and maxima of the variables, one studies the acceleration vector along the level curve, which is given by
It follows that −(∂ℏ/∂p)(∂^{2}ℏ/∂φ^{2}) is positive for a minimum of p and negative in the case of a maximum. Similarly, an extremum of φ is a minimum when −(∂ℏ/∂φ)(∂^{2}ℏ/∂p^{2}) is a positive quantity and a maximum if it is negative.
As an example, the extrema of the momentum p along the level curve ℏ(p, φ) = ℏ_{0} are given by the system of equations
Similarly to the computation of the fixed points, the solutions of Eq. (G.20) can be searched for among those of the system of equations
which involves the solution of a univariatepolynomial equation of degree 2n in the variable p.
G.4. Resonance halfwidth
Following Chirikov (1979), we denote as ω the projection of the frequency vector ω′=dθ′/dt along the normal to the resonance plane k ⋅ ω′+Ω = 0, that is, ω = ∥k∥^{−1}k ⋅ ω′. The reduced dynamics in Eq. (G.12) induces a variation in this projection equal to . Based on ideas from frequency analysis (Laskar 1993), we then associated with the level curve of the reduced Hamiltonian passing through the point (p_{0}, φ_{0}) its rotational frequency ν(p_{0}, φ_{0}), that is, the time mean of the instantaneous frequency ,
where the time integration is meant along the trajectory starting at (p_{0}, φ_{0}) at time τ = 0. The rotational frequency of a separatrix is null, as for the librations of the angle φ inside a resonant island, while ν is different from zero for rotations. Therefore, we established the halfwidth of a resonant harmonic from the variation in the rotational frequency across its separatrix. To associate an intrinsic halfwidth with each harmonic as in the Chirikov (1979) definition, we defined once and for all the twofold halfwidth as
where (p_{±}, φ_{±}) are the locations of the maximum and minimum, respectively, of the momentum p along the separatrix of the resonance (it is implicitly assumed that the line curve through (p_{±}, φ_{±} + π/2) corresponds to a rotation of the angle φ). When applied to the simple pendulum, definition (G.23) gives
where is a constant factor, being the complete elliptic integral of the first kind. As β ≈ 0.95, Eq. (G.24) is practically the Chirikov (1979) halfwidth (recall that ω_{ell} = ω_{hyp} for the pendulum).
To compute the frequency ν of a rotating line curve in a numerically efficient way, we write ν(p_{0}, φ_{0}) = 2π/𝒯(p_{0}, φ_{0}), where the signed period 𝒯 is given by
with . Equation (G.25) assumes that φ is a monotonic function of time along the level curve through (p_{0}, φ_{0}) (i.e. or ), and p = p(φ; p_{0}, φ_{0}) then represents the momentum as a (singlevalued) function of the angle^{3}. To compute the signed period without numerical integration of the corresponding trajectory, we reconstructed the curve p = p(φ; p_{0}, φ_{0}) geometrically. We started by retrieving the minimum, p_{min}, and maximum, p_{max}, of the momentum along the level line (p_{0}, φ_{0}) as solutions of system (G.21). We then linearly sampled values of p in the interval [p_{min}, p_{max}] and retrieved corresponding values for φ through Eq. (G.12), which consists of a linear trigonometric equation in the angle variable. We finally used spline interpolation to uniformly sample the curve along the angle axis and compute the integral in Eq. (G.25).
The twofold definition in Eq. (G.23) attributes a different halfwidth to each side of the resonance along the momentum axis. While the pendulum dynamics is symmetric, the two halfwidths Δω_{±} are not equal in absolute value in the general case. Moreover, their signs depend on the topology of the reduced phase space and are not necessarily opposite (see, as an example, the resonance θ_{1 : 1} at 450 Myr in Fig. 2). The time statistic of the resonance halfwidth along the sampled nominal solution (given in Table 2) is constructed from the sample
where m spans the subset of times the harmonic is resonant. The two sides of the resonance are thus equally weighted in absolute value. More informative statistics are the signed halfwidths, Δω^{+} and Δω^{−}, derived from the samples
whose medians define the asymmetric resonance layers in Fig. 3.
G.5. Time sampling
To retrieve the samples of the fixed point frequencies ω_{hyp}, ω_{ell} and of the resonance halfwidth Δω along the sampled nominal solution , we first applied the lowpass KZ filter to the time series of the actionangle variables,
with three iterations of the moving average and a cutoff frequency of (5 Myr)^{−1}. Similarly to our definition of a librating harmonic in Appendix F, we are not interested in harmonics that are resonant over timescales much shorter than the Lyapunov time. Moreover, lowpass filtering allows the time series to be resampled with a timestep Δt′ = 250 kyr, which is much bigger than the original one, Δt = 1 kyr Appendix C. This resampling is fundamental for actually being able to perform the systematic retrieval of the reduced phasespace topology over the entire nominal solution, which spans 5 Gyr, and for all the librating harmonics. We point out that, because of the linearity of the KZ filter, filtering the angle variables θ′(t) is equivalent to filtering the corresponding instantaneous frequencies dθ′(t)/dt, as done in Appendix F to define timedependent fundamental frequencies.
G.6. Polynomial solver
To find numerical approximations to the roots of a univariate polynomial with complex coefficients, TRIP implements the fixedprecision algorithm of Bini (1996) and extends it to quadruple and multiprecision floating point coefficients. Based on the EhrlichAberth iteration (e.g. Bini 1996), which allows simultaneous approximations of all the roots, the algorithm guarantees a posteriori error bounds. This allows, in particular, the real roots we are interested in to be isolated in a robust way. The algorithm deals with polynomials of high degree and with very large or small coefficients. Such polynomials can be generated in the symbolic preprocessing of systems of equations, as in our case. When a lack of convergence to some of the roots is detected after a fixed maximum number of iterations, or in the rare case of an overflow, we temporarily switched to the more recent implementation of MPSolve 3 (Bini & Fiorentino 2000; Bini & Robol 2014), which we integrated in TRIP and which provides a guaranteed approximation of the roots with any desired number of digits.
Appendix H: Transformed proper modes and time distributions
The proper modes and introduced in Appendix B change under the Lie transform according to
While the transformed proper modes are formally determined by Eq. (H.1), the need to truncate the Lie transform in actual computations produces different definitions, that is,
Even though the Lie transform defines a canonical change of variables close to identity, its slow (finitedegree) convergence shown in Appendix F suggests that the contributions from highdegree terms are actually nonnegligible. Indeed, harmonics of order eight and ten appear among the leading resonances of Table 2, and their amplitudes are largely the result of the Lie transform. The truncated proper modes in Eq. (H.2) may thus be expected to differ from one another in an appreciable way. As the resonances in Table 2 are given in the proper modes truncated at degree ten, we show in Fig. G.1 some examples of the dependence of the time statistic on the degree of truncation of the Lie transform. For the resonant harmonics σ_{1 : 1}, θ_{2 : 1}, θ_{1 : 1}, and (g_{1} − g_{3})+(s_{2} − s_{3}), we report in the corresponding column the time PDF of the hyperbolic fixed point frequency ω_{hyp} (first line), the resonance halfwidth Δω (second line), and the signed halfwidths Δω^{−}, Δω^{+} (third line), for degrees of truncation from six to ten. Each time distribution thus results from the truncated Hamiltonian along the corresponding nominal solution , which spans 5 Gyr. We show the kernel density estimation (Rosenblatt 1956; Parzen 1962) of each PDF (we used the Gaussian kernel and the Silverman (1986) rule of thumb to select the bandwidth^{4}) as well as the observed median and the [5th, 95th] percentile range. Figure G.1 shows that the time distributions of different degrees largely overlap one another. Their medians, in particular, differ by a factor of two at most. These examples indicate that, as the nominal solution varies within the chaotic zone, the resonant nature of a Fourier harmonic can be stated statistically, in a way that is, within the framework of this study, largely independent of the precise truncation of the proper modes in Eq. (H.2). The dependence on the truncation degree does not modify the general picture of the resonance web given in Table 2. Truncation at degree ten, used for that table, allows the appearance of higherorder harmonics to be accounted for. In this context, we point out that one may expect libration episodes of a resonant harmonic to emerge in a statistical framework independently of the precise proper modes and even the dynamical model employed, as shown in Fig. 4. This consideration is clearly restricted to sets of proper modes that differ by a quasiidentity change of variables from those that integrate the LL dynamics (e.g. u, v in this work or z^{•}, ζ^{•} in Laskar 1990a), and to dynamical models that faithfully reproduce the predictions of Nbody numerical integrations.
Fig. G.1. Time PDF of the hyperbolic fixed point frequency ω_{hyp} (first line), resonance halfwidth Δω (second line), and signed halfwidths Δω^{−}, Δω^{+} (third line), for four resonant harmonics listed in Table 2 and degree of truncation 2n from six to ten. Each time distribution results from the truncated Hamiltonian along the nominal solution spanning 5 Gyr. Kernel density estimation of the PDF is shown as a solid line, and the observed median and the [5th, 95th] percentile range are represented by a vertical dashed line and a shaded region, respectively. Frequencies are in arcsec yr^{−1}. 
We remark that we do not provide confidence intervals for the PDF estimation in Fig. G.1. Estimation of confidence intervals must take into account that the time samples are correlated over a length comparable to the Lyapunov time of the dynamics. Confidence intervals for correlated data can be obtained via movingblock bootstrap (Hoang et al. 2021), for example. While one may expect the variance of the PDF estimation to be rather large, we aim here to show the compatibility of the distributions across different degrees of truncation, and this is already indicated by the PDF estimation itself. More precisely, we do not intend to state exact time distributions for ω_{hyp} or Δω (which would be limited by the assumptions behind our approach, by the way) but simply to identify the ranges over which these dynamical quantities may statistically vary. As an example, the overlap between the asymptotic Lyapunov exponent in Fig. 1 and the ranges spanned by ω_{hyp} and Δω in the case of the leading resonances θ_{1 : 1} and θ_{2 : 1} is clearly robust across different degrees of truncation. In this sense, the medians and percentile ranges in Table 2 are not meant to represent strict values for the corresponding observables, but more simply to synthesise distributions such as those shown in Fig. G.1.
All Tables
All Figures
Fig. 1. Finitetime MLE of the forced secular ISS over 5 Gyr for different degrees of truncation of the Hamiltonian and the corresponding Lyapunov time FTMLE^{−1}. The bands represent the [5th, 95th] percentile range of the PDF estimated from ensembles of 128 (ℋ_{2n}) and 1148 (ℋ) stable orbital solutions with very close initial conditions. The black lines denote the nominal solution 𝒮 for nine different initial tangent vectors. 

In the text 
Fig. 2. Reduced phase spaces (p, φ) of the harmonics θ_{1 : 1}, θ_{2 : 1}, and θ_{3 : 2} along the Lietransformed nominal solution as deduced from the Lietransformed Hamiltonian . The phase spaces correspond to the t_{0} times indicated in the upperright boxes. The black dots reproduce the level curves of the reduced Hamiltonian, and the orange dots represent its fixed points. The separatrices are shown by red dots, and the level curve corresponding to the resonant variables at time t_{0} is shown in blue. The righthandside axes report the projection of the mean frequencies along the normal k to the resonance plane. 

In the text 
Fig. 3. Overlaps of secular resonances from Table 2. Left panel: starred resonances θ_{m : n} = m(g_{3} − g_{4}) − n(s_{3} − s_{4}) in the frequency subspace spanned by g_{4} − g_{3} and s_{4} − s_{3}. Right panel: daggered resonances as integer combinations of (g_{1} − g_{5})−(s_{1} − s_{2}) and (g_{1} − g_{2})+(s_{1} − s_{2}) (the harmonic ℱ_{66} = 3g_{1} − g_{2} − 2g_{5} − (s_{1} − s_{2}) is not shown in Table 2). The dashed lines represent the resonance centres, and the coloured regions correspond to the time median of the signed halfwidths. The nominal solution , which spans 5 Gyr, is shown after lowpass filtering with a cutoff frequency of (5 Myr)^{−1}, the cross indicating the current position of the ISS. Frequencies are in arcsec yr^{−1}. 

In the text 
Fig. 4. Combination of angle variables for different harmonics of Table 2 as a function of time along direct integrations of the Solar System presented in Laskar & Gastineau (2009) (n = 0 is the nominal solution). The angle variables correspond to the proper modes z^{•}, ζ^{•} defined in Laskar (1990a). The time series are lowpassfiltered with a cutoff frequency of (500 kyr)^{−1}. 

In the text 
Fig. E.1. Finitetime MLE and corresponding Lyapunov time FTMLE^{−1} from different partial Hamiltonians (Eq. (E.1)). The bands represent the [5th, 95th] percentile range of the PDF estimated from an ensemble of 128 stable orbital solutions with very close initial conditions. The blue colour corresponds to the Fourier harmonics of Table 1, the green colour to the harmonics of the same table that only include the angular DOFs 1, 2, and 5, and the red colour to those harmonics that only involve the angular DOFs 3 and 4. The hatched regions correspond to the case where the choice of the DOFs is made across the entire ℍ_{6} Hamiltonian. The grey region stands for the full ℋ Hamiltonian, as shown in Fig. 1. 

In the text 
Fig. F.1. Convergence of the Lietransformed nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are lowpassfiltered with a cutoff frequency of (5 Myr)^{−1}. 

In the text 
Fig. F.2. Convergence of the Lietransformed Hamiltonian along the nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are lowpassfiltered with a cutoff frequency of (5 Myr)^{−1}. The convergence of the original ℍ_{2n} Hamiltonian along the nominal solution 𝒮 is also shown. 

In the text 
Fig. G.1. Time PDF of the hyperbolic fixed point frequency ω_{hyp} (first line), resonance halfwidth Δω (second line), and signed halfwidths Δω^{−}, Δω^{+} (third line), for four resonant harmonics listed in Table 2 and degree of truncation 2n from six to ten. Each time distribution results from the truncated Hamiltonian along the nominal solution spanning 5 Gyr. Kernel density estimation of the PDF is shown as a solid line, and the observed median and the [5th, 95th] percentile range are represented by a vertical dashed line and a shaded region, respectively. Frequencies are in arcsec yr^{−1}. 

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.