© F. Mogavero and J. Laskar 2022

1. Introduction

The chaotic long-term 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)⁻¹ (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 quasi-periodic series (Laskar 1984, 1988). Investigating the origin of chaos, Laskar measured the libration period of the Fourier harmonic θ_2 : 1 = 2(g₃ − g₄)−(s₃ − s₄), 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₃ − g₄)−(s₃ − s₄) 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 high-dimensional dynamics of the inner planets probably discouraged systematic analytical studies of its resonant structure. A couple of analyses have focused on the long-term motion of Mercury, by freezing all the other planets on quasi-periodic 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 three-body 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 long-term 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 long-term dynamics of the Solar System planets essentially consists of the slow precession of their perihelia and nodes, driven by secular, that is, orbit-averaged, 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 quasi-periodic secular solution for the giant planets, with the inner ones moving in the resulting time-dependent gravitational potential (Mogavero & Laskar 2021). The quasi-periodic 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 mean-motion resonances in the ISS allow us to simply consider first-order secular averaging of the N-body 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)

(1)

The planets are indexed in order of increasing semi-major axis, are the planetary masses, a_k and e_k are the (secular) semi-major 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 non-resonant Fourier harmonics of the N-body Hamiltonian at first order in planetary masses (Mogavero & Laskar 2021). The semi-major axes a_k are constants of motion in the secular dynamics. A suitable set of canonically conjugate momentum-coordinate pairs of variables for the secular dynamics are the Poincaré rectangular coordinates in complex form, , with

(2)

where Λ_k = μ_k[G(m₀ + m_k)a_k]^1/2, m₀ and μ_k = m₀m_k/(m₀ + 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 quasi-periodic time dependence for the orbits of the outer planets (Mogavero & Laskar 2021),

(3)

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) = ω_ot, with ω_o = (g₅, g₆, g₇, g₈, s₆, s₇, s₈) 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),

(4)

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 16-dimensional phase space since no trivial integrals of motion, such as the total energy or angular momentum, exist.

The development of the two-body 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, ℋ₂ produces a forced Laplace-Lagrange (LL) dynamics, which can be analytically integrated by introducing complex proper mode variables, . By introducing action-angle variables through and , the truncated Hamiltonians expressed in the proper modes can be expanded as finite Fourier series,

(5)

where I = (X, Ψ) and θ = (χ, ψ) are the eight-dimensional vectors of the action and angle variables, respectively, and (k, ℓ) is the wave vector of a given harmonic. There are 2748 harmonics with a non-null 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 half-width of the leading resonant harmonics of the Hamiltonian, which constitute the dynamical sources of chaoticity (Chirikov 1979). The non-null probability of unstable orbital evolutions in the ISS (Laskar & Gastineau 2009) makes the definition of the MLE as an infinite-time limit (Oseledec 1968; Benettin et al. 1980) not pertinent. We numerically compute a finite-time MLE (FT-MLE) 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 FT-MLE is time-asymptotically 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 FT-MLE expressed as an angular frequency over the next 5 Gyr for different degrees of truncation of the Hamiltonian. In each case, the FT-MLE is computed for 128 stable (i.e. non-collisional) 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 FT-MLE estimated from each ensemble of solutions. We also report the PDF from the full Hamiltonian ℋ, computed in Mogavero & Laskar (2021), along with the FT-MLE of its nominal solution 𝒮 for nine different initial tangent vectors, to manifest its asymptotic behaviour. In a few hundred million years, each FT-MLE becomes independent of the initial tangent vector, the renormalisation time, and the norm chosen for the phase-space vectors. Its distribution only reveals the intricate dependence on the initial position of the system in the phase space. At 5 Gyr, the FT-MLE roughly ranges from 0.15 to 0.5″ yr⁻¹ with a 90% probability. Figure 1 shows that the truncated Hamiltonians ℋ_2n reproduce the asymptotic distribution of the FT-MLE 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 FT-MLE.

Fig. 1. Finite-time MLE of the forced secular ISS over 5 Gyr for different degrees of truncation of the Hamiltonian and the corresponding Lyapunov time FT-MLE−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 long-term 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

(6)

We ranked the 69 339 harmonics of ℍ₆ 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₁ − g₅)−(s₁ − s₂) 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₁ − g₄)+(s₁ − s₄), associating the proper modes of Mercury and Mars and (s₁ − s₂)+(s₃ − s₄), which concatenates all the inclination DOFs. A non-negligible role of the Saturn-dominated eccentricity mode g₆ 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 FT-MLE 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 non-chaotic dynamics. This notably indicates that the resonant nature of a Fourier harmonic should only be established along the orbital solution of a realistic Hamiltonian.

5. Resonant harmonics

Figure 1 and Table 1 suggest that the non-linear 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 Lie-transformed Hamiltonian, , which is in Birkhoff normal form to degree four, these terms are replaced with the chain of high-order 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 Lie-transformed nominal solution, . Following Mogavero & Laskar (2021), we defined time-dependent 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 divide-et-impera approach of Chirikov (1979) and considered a reduced Hamiltonian for each librating harmonic:

(7)

The resonant variables p, φ are related to the Lie-transformed action-angle variables I′,θ′ by and φ = k ⋅ θ′  +  Ωt, where we denote Ω = ℓ ⋅ ω_o (Appendix G). We therefore inspected the fictitious one-DOF 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 full-dimensional action space. For a given librating harmonic, our study considers the one-parameter 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₀. As an example, Fig. 2 shows, for different t₀ times along , the phase spaces of the harmonics θ_1 : 1, θ_2 : 1, and θ_3 : 2, with θ_m : n = m(g₃ − g₄)−n(s₃ − s₄), 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 FT-MLE 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 semi-analytically 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₀ 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 Lie-transformed nominal solution as deduced from the Lie-transformed Hamiltonian . The phase spaces correspond to the t0 times indicated in the upper-right 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 t0 is shown in blue. The right-hand-side 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 non-linear 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 libration-rotation state of the harmonic and the resonance half-width, Δω. 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∥⁻¹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 half-width Δω 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,

(8)

which is the Chirikov (1979) half-width (recall that ω_ell = ω_hyp for the pendulum).

6. Origin of chaos

To establish a meaningful connection between the resonant phase spaces and the observed FT-MLE, 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 finite-time ergodic assumption, we systematically retrieved for each librating harmonic of the samples of the fixed point frequencies ω_hyp, ω_ell and of the resonance half-width Δω 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 half-width. 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.

It is perhaps hopeless to obtain, for high-dimensional dynamical systems, a precise formula connecting the FT-MLE to the half-width 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,

(9)

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 FT-MLE of the nominal solution in Fig. 1. Moreover, they are significantly compatible with the long-term ensemble distribution of the FT-MLE. 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⁻¹, which constitutes a sounding explanation for the large width of the FT-MLE distribution, while the dense set of terms with Δω < 0.1″ yr⁻¹ explains why the exponent never goes below this value. We point out the absence of resonances of order two. The harmonic g₁ − g₅, 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₁ − g₅, g₂ − g₅, and s₁ − s₂ and includes the well-known resonance σ_1 : 1. We first point out that the libration frequencies of 0.28 and 0.12″ yr⁻¹ 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 one-DOF phase spaces, as in Fig. 2, can be stated in two ways. As time changes and the system moves along its nominal trajectory, the phase-space region swept by the separatrices forms a chaotic zone. In an equivalent way, a chaotic zone results, at fixed time t₀, from separatrix-splitting 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 half-widths Δω⁺, Δω⁻ (Appendix G). The reported half-widths 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₃ − g₄ at degree four (see Table 1). The resonance overlap in the g₁ − g₅, g₂ − g₅, s₁ − s₂ subspace is shown in Fig. 3 in terms of integer combinations of σ_1 : 1 and (g₁ − g₂)+(s₁ − s₂). It appears to be more relevant, at least for our nominal solution, than the restriction to the subspace spanned by g₁ − g₅ and s₁ − s₂ investigated in previous studies (Lithwick & Wu 2011; Batygin et al. 2015). A relevant dynamical role of the proper mode g₂, hinted at in Lithwick & Wu (2011) and Batygin et al. (2015), is confirmed. Figure 3 also suggests a long-term 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. non-collisional) solution.

Fig. 3. Overlaps of secular resonances from Table 2. Left panel: starred resonances θm : n = m(g3  −  g4)  −  n(s3  −  s4) in the frequency subspace spanned by g4 − g3 and s4 − s3. Right panel: daggered resonances as integer combinations of (g1 − g5)−(s1 − s2) and (g1 − g2)+(s1 − s2) (the harmonic ℱ66 = 3g1 − g2 − 2g5 − (s1 − s2) 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 half-widths. The nominal solution , which spans 5 Gyr, is shown after low-pass 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⁻¹ significantly couple all the proper modes, confirming the implications of Table 1. Coupling resonances such as (g₂ − g₄)+(s₂ − s₄), (g₁ − g₃)+(s₂ − s₃) and (g₁ − g₄)+(s₁ − s₄), along with high-order ones resulting from their interaction with leading harmonics, strongly suggest that resonance overlap takes place extensively in the full-dimensional frequency space. Even though we have isolated low-dimensional sources of chaoticity, the large-scale chaos in the ISS is probably best understood as a high-dimensional 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₁ − g₂)+(s₁ − s₂) 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₁ − g₃)+(s₂ − s₃), 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 low-pass-filtered with a cutoff frequency of (500 kyr)−1.

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 post-doctoral fellowship and by a grant of the French Agence Nationale de la Recherche (AstroMeso ANR-19-CE31-0002-01). This project has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo-885250). This work was granted access to the HPC resources of MesoPSL financed by the Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

References

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 one-DOF 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 low-order 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 one-DOF model, still depends on the system position in the full-dimensional phase space. Some authors have invoked the quasi-periodicity of certain variables to reduce the analysis to a lower-dimension 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₂, which has slower frequency variations (Laskar 1990a; Laskar et al. 2004; Hoang et al. 2021). A quasi-periodic 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 high-dimensional 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¹,

(A.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 multiple-precision floating-point numbers; fixed or multiple-precision integers or rational numbers; or double- and quadruple-precision floating-point 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 state-of-the-art 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,

(B.1)

At the lowest degree, the ℋ₂ 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 momentum-coordinate 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

(B.2)

where ℍ_(2p) groups all the terms of same total degree 2p (Mogavero & Laskar 2021). Action-angle variables that trivially integrate the LL problem are thus introduced as

(B.3)

for k ∈ {1, 2, 3, 4} and allow the truncated Hamiltonian in Eq. (B.2) to be written as a finite Fourier series,

(B.4)

where are complex amplitudes, with , and where we employ a compact notation for the action-angle 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 quasi-integrable 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 phase-space extension, resulting in

(B.5)

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 FT-MLE by taking the relative variation in each coordinate of the nominal phase-space vector of Gauss’s dynamics as a normal random variable, with zero mean and a standard deviation of 10⁻⁹. In other words, the initial conditions are distributed according to

(C.1)

for k ∈ {1, 2, 3, 4}, where represents the initial conditions of the nominal solution 𝒮, are standard normal deviates, and σ = 10⁻⁹. 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 phase-space 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 long-term variation in the frequencies of motion . From Eq. (B.4), Hamilton’s equations give

(D.1)

where (k, ℓ)∈ℤ⁸ × ℤ⁷ 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 slow-varying mean frequencies , only depending on the action variables. If the motion were quasi-periodic, 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 short-time average of ω_2n with respect to . Nevertheless, the mean frequencies still provide the overall long-term 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

(D.2)

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,

(D.3)

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 FT-MLE to be reproduced from a small set 𝒰 of Fourier harmonics, by considering

(E.1)

Partial Hamiltonians can be systematically constructed in TRIP through the selection of specific terms of ℍ₆. We integrated the corresponding variational equations to obtain the FT-MLE 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 FT-MLE distribution from the 30 leading harmonics of Table 1. This simplified dynamics reproduces the asymptotic FT-MLE distribution of the Hamiltonians ℋ and ℋ_2n very well, confirming the dynamical relevance of the leading harmonics. Very importantly, the FT-MLE prediction turns out to be robust with respect to the addition of further harmonics from our ranking.

Fig. E.1. Finite-time MLE and corresponding Lyapunov time FT-MLE−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 FT-MLE 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 𝒰 = {ℱ₆, ℱ₈, ℱ₂₀, ℱ₂₇, ℱ₂₉, ℱ₃₀} only involving combinations of g₁, g₂, g₅, s₁, and s₂. This partial Hamiltonian includes, in particular, the Fourier harmonics considered in the simplified long-term 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₁, g₂, g₅, s₁, and s₂ is made across the entire ℍ₆ 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 𝒰 = {ℱ₂, ℱ₃, ℱ₅, ℱ₂₅} only involving g₃, g₄, s₃ and s₄, produce an asymptotic FT-MLE 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 ℍ₆ Hamiltonian, the resulting FT-MLE 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 ℍ₆ 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 non-null wave vector appearing in ℍ₍₄₎ in Eq. (B.2). Such a transformation of variables can be canonically defined as the time-1 flow of a generating Hamiltonian S that satisfies the homologic equation

(F.1)

where the braces represent the Poisson bracket and where we employ the extended phase-space formalism of Eq. (B.5) to deal with a time-independent Hamiltonian. The resulting change of variables is given by the formal Lie transform

(F.2)

where L_S  ⋅ ={S, ⋅} is the Lie derivative associated with the generating function S. The exponential of an operator A acting on the phase-space functions is formally expressed as

(F.3)

with A⁰ defined as the identity operator and A^q = AA^q − 1 for q ≥ 1. The transformed Hamiltonian reads

(F.4)

where the Lie transform is formally evaluated at the new canonical variables I′,θ′;Φ′,ϕ′. The solution to the homologic Eq. (F.1) is expressed as

(F.5)

As ℍ₍₄₎ possesses a finite number of non-null 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 well-defined 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^−L_SΦ 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,

(F.6)

Once the Hamiltonian is computed, the corresponding nominal solution is obtained from the truncation of the operator E^−L_S 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^±L_S 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 phase-space vectors, as usual. To highlight the long-term trend of the time series, which may be hidden by large short-time oscillations, we applied the low-pass Kolmogorov-Zurbenko (KZ) filter with three iterations of the moving average and a cutoff frequency of (5 Myr)⁻¹ (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 Lie-transformed 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 (finite-degree) convergence of the series.

Fig. F.1. Convergence of the Lie-transformed nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are low-pass-filtered with a cutoff frequency of (5 Myr)−1.

Fig. F.2. Convergence of the Lie-transformed 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 low-pass-filtered 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 time-dependent fundamental frequencies along a sampled orbital solution, as we did in Mogavero & Laskar (2021) in the case of the frequency g₁. 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 five-point formula, as very high numerical precision is unnecessary. Since the instantaneous frequencies present large short-time fluctuations that hide their long-term evolution, we applied the low-pass 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)⁻¹. This specific value was chosen to filter out the non-chaotic 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 Lie-transformed nominal solution at degree six, eight, and ten, respectively.

Appendix G: Reduced Hamiltonians

To retrieve the resonant harmonics of the Lie-transformed Hamiltonian , that is,

(G.1)

we considered the dynamics generated by each Fourier component separately. If we assume that at time t = t₀ 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

(G.2)

where τ ≥ t₀ is the time of the fictitious dynamics, Ω = ℓ ⋅ ω_o, and we consider k ∈ ℤ⁸ ∖ {0}. The resonant nature of the harmonic depends on the closeness of the actions to the resonance condition

(G.3)

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,

(G.4)

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 time-dependent generating function

(G.5)

where 𝕄 is a real 8 × 8 matrix, , and T denotes the transposition operator. One has

(G.6)

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

(G.7)

with p = (p₁, …, p₈)^T. We finally chose ν = (Ω, 0, …, 0)^T so that φ₁ = k ⋅ θ′+Ωτ. The reduced Hamiltonian governing the new canonical variables (p, φ) is given by

(G.8)

The reduced dynamics only affects the momentum p₁, and are all integrals of motion. These constants can be set to zero by choosing and employing the initial condition so that p₁(τ = t₀) = 0 and

(G.9)

As the action variables I′ are non-negative by definition, the dynamics of the momentum p₁ 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,

(G.10)

With this choice, the integrals of motion have null values and the action variables are given by

(G.11)

To simplify the notation, we omit the subscripts of the resonant variables (p₁, φ₁) 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 non-linear 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₀), which thus possesses up to n − 1 real solutions. As a generalisation of the pendulum approximation, when no real solutions for p(t₀) 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 non-resonant. When a real solution for p(t₀) exists and the corresponding point lies in the physical action space (that is, its components are all non-negative), 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 eight-dimensional action space, that is, it depends on through Eq. (G.9) or Eq. (G.10). An explicit study of the one-DOF 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 one-parameter family of phase spaces spanned by the time t₀. 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

(G.12)

where we rename , , and as ℏ, f₀, and f₁, respectively, to keep a simpler notation. From Eq. (B.3) it follows that one can write , with

(G.13)

where we denote and . The function μ is thus a real multivariate monomial of degree in the square roots of the action variables, while F₁ can be expressed as a multivariate polynomial of degree (2n − |k|−|ℓ|)/2 in the action variables², the order of the harmonic k, ℓ (i.e. |k|+|ℓ|) being an even integer. The fixed points of the Hamiltonian ℏ satisfy the system of equations

(G.14)

where d stands for derivation with respect to p. The function df₀ 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₁ and in its derivative df₁, the system of equations (G.14) is non-algebraic. 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 non-null 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₁ can be eliminated through multiplication by the function

(G.15)

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

(G.16)

The second equation is now a univariate-polynomial 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

(G.17)

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 λ² = −det(𝕍). In the case of an elliptic fixed point (i.e. when λ² < 0), the angular frequency of the small oscillations is given by , while in the case of a hyperbolic fixed point (i.e. when λ² > 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, φ) = ℏ₀ can be expressed as

(G.18)

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

(G.19)

It follows that −(∂ℏ/∂p)(∂²ℏ/∂φ²) is positive for a minimum of p and negative in the case of a maximum. Similarly, an extremum of φ is a minimum when −(∂ℏ/∂φ)(∂²ℏ/∂p²) 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, φ) = ℏ₀ are given by the system of equations

(G.20)

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

(G.21)

which involves the solution of a univariate-polynomial equation of degree 2n in the variable p.

G.4. Resonance half-width

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∥⁻¹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₀, φ₀) its rotational frequency ν(p₀, φ₀), that is, the time mean of the instantaneous frequency ,

(G.22)

where the time integration is meant along the trajectory starting at (p₀, φ₀) 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 half-width of a resonant harmonic from the variation in the rotational frequency across its separatrix. To associate an intrinsic half-width with each harmonic as in the Chirikov (1979) definition, we defined once and for all the two-fold half-width as

(G.23)

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

(G.24)

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) half-width (recall that ω_ell = ω_hyp for the pendulum).

To compute the frequency ν of a rotating line curve in a numerically efficient way, we write ν(p₀, φ₀) = 2π/𝒯(p₀, φ₀), where the signed period 𝒯 is given by

(G.25)

with . Equation (G.25) assumes that φ is a monotonic function of time along the level curve through (p₀, φ₀) (i.e. or ), and p = p(φ; p₀, φ₀) then represents the momentum as a (single-valued) function of the angle³. To compute the signed period without numerical integration of the corresponding trajectory, we reconstructed the curve p = p(φ; p₀, φ₀) geometrically. We started by retrieving the minimum, p_min, and maximum, p_max, of the momentum along the level line (p₀, φ₀) 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 two-fold definition in Eq. (G.23) attributes a different half-width to each side of the resonance along the momentum axis. While the pendulum dynamics is symmetric, the two half-widths Δω_± 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 half-width along the sampled nominal solution (given in Table 2) is constructed from the sample

(G.26)

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 half-widths, Δω⁺ and Δω⁻, derived from the samples

(G.27)

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 half-width Δω along the sampled nominal solution , we first applied the low-pass KZ filter to the time series of the action-angle variables,

(G.28)

with three iterations of the moving average and a cutoff frequency of (5 Myr)⁻¹. 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, low-pass 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 phase-space 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 time-dependent fundamental frequencies.

G.6. Polynomial solver

To find numerical approximations to the roots of a univariate polynomial with complex coefficients, TRIP implements the fixed-precision algorithm of Bini (1996) and extends it to quadruple- and multi-precision floating point coefficients. Based on the Ehrlich-Aberth 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

(H.1)

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,

(H.2)

Even though the Lie transform defines a canonical change of variables close to identity, its slow (finite-degree) convergence shown in Appendix F suggests that the contributions from high-degree terms are actually non-negligible. 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₁ − g₃)+(s₂ − s₃), we report in the corresponding column the time PDF of the hyperbolic fixed point frequency ω_hyp (first line), the resonance half-width Δω (second line), and the signed half-widths Δω⁻, Δω⁺ (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⁴) 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 higher-order 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 quasi-identity 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 N-body numerical integrations.

Fig. G.1. Time PDF of the hyperbolic fixed point frequency ωhyp (first line), resonance half-width Δω (second line), and signed half-widths Δω−, Δω+ (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 moving-block 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. Finite-time MLE of the forced secular ISS over 5 Gyr for different degrees of truncation of the Hamiltonian and the corresponding Lyapunov time FT-MLE−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 Lie-transformed nominal solution as deduced from the Lie-transformed Hamiltonian . The phase spaces correspond to the t0 times indicated in the upper-right 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 t0 is shown in blue. The right-hand-side 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(g3  −  g4)  −  n(s3  −  s4) in the frequency subspace spanned by g4 − g3 and s4 − s3. Right panel: daggered resonances as integer combinations of (g1 − g5)−(s1 − s2) and (g1 − g2)+(s1 − s2) (the harmonic ℱ66 = 3g1 − g2 − 2g5 − (s1 − s2) 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 half-widths. The nominal solution , which spans 5 Gyr, is shown after low-pass 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 low-pass-filtered with a cutoff frequency of (500 kyr)−1.
In the text

Fig. E.1. Finite-time MLE and corresponding Lyapunov time FT-MLE−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 Lie-transformed nominal solution . Left panel: Relative increments as a function of time and truncation degree 2n. Right panel: Corresponding time means . The time series are low-pass-filtered with a cutoff frequency of (5 Myr)−1.
In the text

	Fig. F.2. Convergence of the Lie-transformed 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 low-pass-filtered 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 half-width Δω (second line), and signed half-widths Δω−, Δω+ (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