Timescales of Chaos in the Inner Solar System: Lyapunov Spectrum and Quasi-integrals of Motion

Numerical integrations of the Solar System reveal a remarkable stability of the orbits of the inner planets over billions of years, in spite of their chaotic variations characterized by a Lyapunov time of only 5 million years and the lack of integrals of motion able to constrain their dynamics. To open a window on such long-term behavior, we compute the entire Lyapunov spectrum of a forced secular model of the inner planets. We uncover a hierarchy of characteristic exponents that spans two orders of magnitude, manifesting a slow-fast dynamics with a broad separation of timescales. A systematic analysis of the Fourier harmonics of the Hamiltonian, based on computer algebra, reveals three symmetries that characterize the strongest resonances responsible for the orbital chaos. These symmetries are broken only by weak resonances, leading to the existence of quasi-integrals of motion that are shown to relate to the smallest Lyapunov exponents. A principal component analysis of the orbital solutions independently confirms that the quasi-integrals are among the slowest degrees of freedom of the dynamics. Strong evidence emerges that they effectively constrain the chaotic diffusion of the orbits, playing a crucial role in the statistical stability over the Solar System lifetime.


I. INTRODUCTION
The planetary orbits in the inner Solar System (ISS) are chaotic, with a Lyapunov time distributed around 5 million years (Myr) [1][2][3][4].Still, they are statistically very stable over a timescale that is a thousand times longer.The probability that the eccentricity of Mercury exceeds 0.7, leading to catastrophic events (i.e., close encounters, collisions, or ejections of planets), is only about 1% over the next 5 billion years (Gyr) [5][6][7].The dynamical halflife of Mercury orbit has recently been estimated at 30-40 billion years [4,7].A disparity of nearly four orders of magnitude between the Lyapunov time and the timescale of dynamical instability is intriguing, since the chaotic variations of the orbits of the inner planets cannot be constrained a priori.While the total energy and angular momentum of the Solar System are conserved, the disproportion of masses between the outer and inner planets implies that unstable states of the ISS are in principle easily realizable through exchanges of these quantities.The surprising stability of the ISS deserves a global picture in which it can emerge more naturally.
To our knowledge, the only study addressing the timescale separation in the long-term dynamics of the ISS is based on the simplified secular dynamics of a massless Mercury [8]: All the other planets are frozen on regular quasi-periodic orbits; secular interactions are expanded to first order in masses and degree 4 in eccentricities and inclinations; an a priori choice of the relevant terms of the Hamiltonian is made.The typical instability time of about 1 Gyr [8,9] is, however, too short and in significant contrast with realistic numerical integrations of the Solar System, which show a general increase of the instability rate with the complexity of the dynamical model [7].We have shown that truncating the secular Hamiltonian of the ISS at degree 4 in eccentricities and inclinations results in an even more stable dynamics, with an instability rate at 5 Gyr that drops by orders of magnitude when compared to the full system [10].From the perspective of these latest findings, the small probability of 1% of an instability over the age of the Solar System may be naturally regarded as a perturbative effect of terms of degree 6 and higher.Clearly, the striking stability of the dynamics at degree 4 is even more impressive in the present context, and remains to be explained.
A strong separation in dynamical timescales is not uncommon among classical quasi-integrable systems [e.g., 11,12].This is notably evinced by the Fermi-Pasta-Ulam-Tsingou (FPUT) problem, which deals with a chain of coupled weakly-anharmonic oscillators [13].Far from Kolmogorov-Arnold-Moser (KAM) and Nekhoroshev regimes (as is likely to be pertinent to the ISS, see Sect.III), one can generally state that the exponential divergence of close trajectories occurring over a Lyapunov time is mostly tangent to the invariant tori defined by the action variables of the underlying integrable problem, and hence contributes little to the diffusion in the action space [14,15].In other words, the Lyapunov time and the diffusion/instability time scale differently with the size of the terms that break integrability, and this can result in very different timescales [12].However, this argument is as general as poorly satisfactory in addressing quantitatively the timescale separation in a complex problem as the present one.Moreover, even though order-of-magnitude estimates of the chaotic diffusion in the ISS suggest that it may take hundreds of million years to reach the destabilizing secular resonance g 1 − g 5 [16], the low probability of an instability over 5 Gyr still remains unexplained [4].Establishing more pre-cisely why the ISS is statistically stable over a timescale comparable to its age is a valuable step in understanding the secular evolution of planetary systems through metastable states [4, 17][18].With its 8 secular degrees of freedom (d.o.f.), this system also constitutes a peculiar bridge between the low-dimensional dynamics often addressed in celestial mechanics and the systems with a large number of bodies studied in statistical mechanics: It cannot benefit from the straightforward application of standard methods of the two fields [e.g., 19,Appendix A].
This work aims to open a window on the long-term statistical behavior of the inner planet orbits.Section II briefly recalls the dynamical model of forced secular ISS introduced in Ref. [4].Section III presents the numerical computation of its Lyapunov spectrum.Section IV introduces the quasi-symmetries of the resonant harmonics of the Hamiltonian and the corresponding quasiintegrals (QIs) of motion.Section V establishes a geometric connection between the quasi-integrals and the slowest d.o.f. of the dynamics via a principal component analysis (PCA) of the orbital solutions.Section VI states the implications of the new findings on the long-term stability of the ISS.We finally discuss the connections with other classical quasi-integrable systems and the methods used in this work.

II. DYNAMICAL MODEL
The long-term dynamics of the Solar System planets consists essentially of the slow precession of their perihelia and nodes, driven by secular, orbit-averaged gravitational interactions [2,20].At first order in planetary masses, the secular Hamiltonian corrected for the leading contribution of general relativity reads [e.g. 4, 21] The planets are indexed in order of increasing semi-major axes (a i ) 8  i=1 , m 0 and m i are the Sun and planet masses, respectively, e i the eccentricities, G the gravitational constant and c the speed of light.The vectors r i are the heliocentric positions of the planets, and the bracket operator represents the averaging over the mean longitudes resulting from the elimination of the non-resonant Fourier harmonics of the N -body Hamiltonian [4,21].Hamiltonian (1) generates Gauss's dynamics of Keplerian rings [4,22], whose semi-major axes a i are constants of motion of the secular dynamics.
By developing the 2-body perturbing function [23,24] in the computer algebra system TRIP [25,26], the secular Hamiltonian can be systematically expanded in series of the Poincaré rectangular coordinates in complex form, where being the reduced masses of the planets, I i the inclinations, i the longitudes of the perihelia and Ω i the longitudes of the nodes [27].Pairs (x i , −jx i ) and (y i , −jy i ) are canonically conjugate momentum-coordinate variables.When truncating at a given total degree 2n in eccentricities and inclinations, the expansion provides Hamiltonians H 2n = H 2n [(x i , xi , y i , ȳi ) 8 i=1 ] that are multivariate polynomials.
Valuable insight into the dynamics of the inner planets is provided by the model of a forced ISS recently proposed [4].It exploits the great regularity of the long-term motion of the outer planets [2,20,28] to predetermine their orbits to a quasi-periodic form: for i ∈ {5, 6, 7, 8}, where t denotes time, xil and ỹil are complex amplitudes, m il and n il integer vectors, and ω o = (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ) represents the septuple of the constant fundamental frequencies of the outer orbits.Frequencies and amplitudes of this Fourier decomposition are established numerically by frequency analysis [29,30] of a comprehensive orbital solution of the Solar System [4, Appendix D].Gauss's dynamics of the forced ISS is obtained by substituting the predetermined time dependence in Eq. ( 1), so that H = H[(x i , y i ) 4 i=1 , t].The resulting dynamics consists of two d.o.f. for each inner planet, corresponding to the x i and y i variables, respectively.Therefore, the forced secular ISS is described by 8 d.o.f. and an explicit time dependence.As a result of the forcing from the outer planets, no trivial integrals of motion exist and its orbital solutions live in a 16-dimensional phase space.
A truncated Hamiltonian H 2n for the forced ISS is readily obtained by substituting Eq. ( 3) in the truncated Hamiltonian H 2n of the entire Solar System.At the lowest degree, H 2 generates a linear, forced Laplace-Lagrange (LL) dynamics.This can be analytically integrated by introducing complex proper mode variables When expressed in the proper modes, the truncated Hamiltonian can be expanded as a finite Fourier series: where I = (X, Ψ) and θ = (χ, ψ) are the 8-dimensional vectors of the action and angle variables, respectively, and we introduce the external angles φ(t) = −ω o t.The The quasi-periodic form of the outer orbits in Eq. (3) contains harmonics of order higher than one, that is, m il 1 > 1 and n il 1 > 1 for some i and l, where • 1 denotes the 1-norm.Therefore, the dynamics of H 2n and H 2n are not exactly the same [4].Still, the difference is irrelevant for the results of this work, so we treat the two Hamiltonians as equivalent from now on.Despite the simplifications behind Eqs. ( 1) and ( 3), the forced secular ISS has been shown to constitute a realistic model that is consistent with the predictions of reference integrations of the Solar System [2,5,6,20].It correctly reproduces the finite-time maximum Lyapunov exponent (FT-MLE) and the statistics of the high eccentricities of Mercury over 5 Gyr [4].Table I presents a summary of the different Hamiltonians and corresponding dynamics we consider in this work.

III. LYAPUNOV SPECTRUM
Ergodic theory provides a way, through the Lyapunov characteristic exponents (LCEs), to introduce a fundamental set of timescales for any differentiable dynamical system ż = F (z, t) defined on a phase space P ⊆ R P [31][32][33][34].If Φ(z, t) denotes the associated flow and z(t) = Φ(z 0 , t) the orbit that emanates from the initial condition z 0 , the LCEs λ 1 ≥ λ 2 ≥ • • • ≥ λ P are the logarithms of the eigenvalues of the matrix Λ(z 0 ) defined as where M(z 0 , t) = ∂Φ/∂z 0 is the fundamental matrix and T stands for transposition [32,33].Introducing the Jacobian J = ∂F /∂z, the fundamental matrix allows us to write the solution of the variational equations ζ = J(z(t), t)ζ as ζ(t) = M(z 0 , t)ζ 0 , where ζ(t) ∈ T z(t) P belongs to the tangent space of P at point z(t) and ζ 0 = ζ(0).The multiplicative ergodic theorem of Oseledec [31] states that if ρ is an ergodic (i.e.invariant and indecomposable) measure for the time evolution and has compact support, then the limit in Eq. ( 7) exists for ρ-almost all z 0 , and the LCEs are ρ-almost everywhere constant and only depend on ρ [32].Moreover, one has lim for ρ-almost all z 0 , where λ (1) > λ (2) > . . .are the LCEs without repetition by multiplicity, and z0 is the subspace of R P corresponding to the eigenvalues of Λ(z 0 ) that are smaller than or equal to exp λ (i) , with 8) is irrelevant [32,34].Once the LCEs have been introduced, a characteristic timescale can be defined from each positive exponent as λ −1 i .In the case of the maximum Lyapunov exponent λ 1 , the corresponding timescale is commonly called the Lyapunov time.
For a Hamiltonian system with p d.o.f.(i.e., P = 2p), the fundamental matrix is symplectic and the set of LCEs is symmetric with respect to zero, that is, If the Hamiltonian is time independent, a pair of exponents vanishes.In general, the existence of an integral of motion C = C(z) implies a pair of null exponents, one of them being associated with the direction of the tangent space that is normal to the surface of constant C [e.g.33].
The ISS is a clear example of a dynamical system that is out of equilibrium.Its phase-space density diffuses seamlessly over any meaningful timescale [5,28].Therefore, the infinite time limit in Eq. ( 7) is not physically relevant.The non-null probability of a collisional evolution of the inner planets [5,6,35,36] implies that such limit does not even exist as a general rule.Most of the orbital solutions stemming from the current knowledge of the Solar System are indeed asymptotically unstable [4,7].Physically relevant quantities are the finite-time LCEs (FT-LCEs), λ i (z 0 , t), defined from the eigenvalues The time dependence of the phase-space density translates in the fact that no ergodic measure is realized by the dynamics, and the FT-LCEs depend on the initial condition z 0 in a non-trivial way [37].
The FT-MLE of the forced secular ISS has been numerically computed over 5 Gyr for an ensemble of stable orbital solutions of the Hamiltonian H with initial conditions very close to their nominal values [38].Its long-term distribution is quite large and does not shrink over time [4,Fig. 3].At 5 Gyr, the probability density function (PDF) of the Lyapunov time peaks at around 4 Myr, it decays very fast below 2 Myr, while its 99th percentile reaches 10 Myr [4, Fig. 4].The significant width of the distribution relates to the aforementioned out-of-equilibrium dynamics of the ISS, as the FT-MLE of each orbital solution continues to vary over time.The dependence of the exponent on the initial condition is associated with the non-ergodic exploration of the phase space by the dynamics.As a remark, the fact that the lower tail of the FT-MLE distribution, estimated from more than 1000 solutions, does not extend to zero implies that invariant (KAM) tori are rare in a neighborhood of the nominal initial conditions (if they exist at all).This fact excludes that the dynamics is in a Nekhoroshev regime [12,39], in agreement with the indications of a multidimensional resonance overlapping at the origin of chaos [19,40].In such a case, the long dynamical half-life of the ISS should not be interpreted in terms of an exponentially-slow Arnold diffusion.
Computations of the FT-MLE of the Solar System planets have been reported for more than thirty years [1,3].However, the retrieval of the entire spectrum of exponents still represents a challenging task.Integrating an N -body orbital solution for the Sun and the eight planets that spans 5 Gyr requires the order of a month of wall-clock time [e.g.41].The computation by a standard method of the entire Lyapunov spectrum for a system with p d.o.f. also requires the simultaneous time evolution of a set of 2p tangent vectors [42].On the top of that, a computation of the exponents for an ensemble of trajectories is advisable for a non-ergodic dynamics [4].These considerations show how demanding the computation of the Lyapunov spectrum of the Solar System planets is.By contrast, a 5-Gyr integration of the forced ISS takes only a couple of hours for Gauss's dynamics (H) and a few minutes at degree 4 (H 4 ).This dynamical model thus provides a unique opportunity to compute all the FT-LCEs that are mainly related to the secular evolution of the inner orbits.
We compute the Lyapunov spectrum of the truncated forced ISS using the standard method of Benettin et al. [43], based on Gram-Schmidt orthogonalization.Manipulation of the truncated Hamiltonian H 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 years.Parallelization of the time evolution of the 16 tangent vectors, between two consecutive reorthonormalization steps of the Benettin et al. [43] algorithm, significantly reduces the computation time.Figure 1a shows the positive FT-LCEs expressed as angular frequencies over the next 10 Gyr for the Hamiltonian truncated at degree 4. The FT-LCEs are computed for 150 stable solutions, with initial conditions very close to the nominal values of Gauss's dynamics and random sets of initial tangent vectors [19,Appendix C].The figure shows the [5th, 95th] percentile range of the marginal PDF of each exponent estimated from the ensemble of solutions.For large times, the exponents of each solution become independent of the initial tangent vectors, the renormalization time, and the norm chosen for the phase-space vectors (see Appendix A and Fig. 8a).In this asymptotic regime, the Benettin et al. [43] algorithm purely retrieves the FT-LCEs as defined in Eq. (10), and the width of their distributions only reflects the out-of-equilibrium dynamics of the system.The convergence of our numerical computation is also assessed by verifying the symmetry of the spectrum stated in Eq. ( 9) (see Appendix A and Fig. 8b).
The spectrum in Fig. 1a has distinctive features.A set of intermediate exponents follow the MLE, ranging from 0.1 to 0.01 yr −1 , while the smallest ones fall below 0.01 yr −1 .Figure 1a reveals the existence of a hierarchy of exponents and corresponding timescales that spans two orders of magnitude, down to a median value of λ −1 8 ≈ 500 Myr.The number of positive exponents confirms that no integral of motion exists, as one may expect from the forcing of the outer planets.We also compute the spectrum for the Hamiltonian truncated at degree 6.As shown in Appendix A (Fig. 9), the asymptotic distributions of the exponents are very similar to those at degree 4.This result suggests that long-term diffusion of the phase-space density is very close in the two cases.The different instability rates of the two truncated dynamics mainly relates to the geometry of the instability boundary, which is closer to the initial position of the system for H 6 than for H 4 [7].
The relevance of the Lyapunov spectrum in Fig. 1a emerges from the fact that the existence of an integral of motion implies a pair of vanishing exponents.This is a pivotal point: By a continuity argument, the presence of positive exponents much smaller than the leading one constitutes a compelling indication that there are dynamical quantities whose chaotic decoherence over initially very close trajectories takes place over timescales much longer than the Lyapunov time.In the long term, such quantities should diffuse much more slowly than any LL action variable.Therefore, Fig. 1a suggests that the secular orbits of the inner planets are characterized by a slowfast dynamics that is much more pronounced than the well-known timescale separation arising from the LL integrable approximation.The existence of slow quantities, which are a priori complicated functions of the phasespace variables, is crucial in the context of finite-time stability, as they can effectively constrain the long-term diffusion of the phase-space density toward the unstable states.The next section addresses the emergence of these slow quantities from the symmetries of the Fourier harmonics that compose the Hamiltonian.

IV. QUASI-INTEGRALS OF MOTION
The emergence of a chaotic behavior of the planetary orbits can be explained in terms of the pendulum-like dynamics generated by each Fourier harmonic that composes the Hamiltonian in Eq. ( 6) [e.g.44].One can write H 2n (I, θ, t) = H 0,0 2n (I) +

M2n
i=1 F i (I, θ, t), with where (k i , i ) = (0, 0), M 2n is the number of harmonics in H 2n with a non-null wave vector, and c.c. stands a O is the order of the harmonic.b τ res is the fraction of time the harmonic is resonant.Only harmonics with τ res > 1% are shown.c 5th and 95th percentiles of the time distribution of ∆ω as subscripts and superscripts, respectively.
for complex conjugate.Chaos arises from the interaction of resonant harmonics, that is, those harmonics F i whose frequency combination k i • θ + i • φ(t) vanishes at some point along the motion.Using the computer algebra system TRIP, the harmonics of H 10 that enter into resonance along the 5-Gyr nominal solution of Gauss's dynamics have been systematically retrieved, together with the corresponding time statistics of the resonance half-widths ∆ω [19].The resonances have then been ordered by decreasing time median of their halfwidths.The resulting ranking of resonances is denoted as R 1 from now on.Table II recalls the 30 strongest resonances that are active for more than 1% of the 5-Gyr time span of the orbital solution.The wave vector of each harmonic is identified by the corresponding combination of frequency labels the order of each harmonic, defined as the even integer O = (k, ) 1 .The support of the asymptotic ensemble distribution of the FT-MLE shown in Fig. 1a overlaps in a robust way with that of the time distribution of the half-width of the strongest resonances.In other words, where ∆ω R1 stands for the half-width of the uppermost resonances of ranking R 1 .Equation ( 12) shows the dynamical sources of chaos in the ISS by connecting the top of the Lyapunov spectrum with the head of the resonance spectrum.Computer algebra allows us to establish such a connection in an unbiased way despite the multidimensional nature of the dynamics.We stress that such analysis is built on the idea that the time statistics of the resonant harmonics along a 5-Gyr ordinary orbital solution should be representative of their ensemble statistics (defined by a set of stable solutions with very close initial conditions) at some large time of the order of billions of years.This assumption was inspired by the good level of stationarity that characterizes the ensemble distribution of the MLE beyond 1 Gyr [4,19], and that extends to the entire spectrum in Fig. 1a.We remark that, strictly speaking, ranking R 1 is established on the Fourier harmonics of the Lie-transformed Hamiltonian H 2n [19,Appendix G].New canonical variables are indeed defined to transform H 2n in a Birkhoff normal form to degree 4. The goal is to let the interactions of the terms of degree 4 in H 2n appear more explicitly in the amplitudes of the harmonics of higher degrees in H 2n , the physical motivation being that the non-linear interaction of the harmonics at degree 4 constitutes the primary source of chaos [19].Keeping in mind the quasiidentity nature of the Lie transform, here we drop for simplicity the difference between the two Hamiltonians.Moreover, all the new analyses of this work involve the original variables of Eq. (5).

A. Quasi-symmetries of the resonant harmonics
In addition to the dynamical interactions responsible for the chaotic behavior of the orbits, Table II provides information on the geometry of the dynamics in the action variable space.Ranking the Fourier harmonics allows us to consider partial Hamiltonians constructed from a limited number m of leading terms [7,19], that is, The dynamics of a Hamiltonian reduced to a small set of harmonics is generally characterized by several symmetries and corresponding integrals of motion.We are interested in how these symmetries are progressively destroyed when one increases the number of terms taken into account in Eq. (13).Consider a set of m harmonics of H 2n and a dynamical quantity that is a linear combination of the action variables, that is, γ ∈ R 8 being a parameter vector.From Eq. ( 11), the partial contribution of the m harmonics to the time derivative of C γ along the flow of and Ċγ = Ċγ,M2n , where M 2n is the total number of harmonics with a non-null wave vector that appear in H 2n .Any quantity C γ with γ • k i = 0 is conserved by the one-d.o.f.dynamics generated by the single harmonic F i .In other words, such a quantity would be an integral of motion if F i were the only harmonic to appear in the Hamiltonian.Considering now m different harmonics, these do not contribute to the change of the quantity , that is, if the vector γ belongs to the orthogonal complement of the linear subspace of R 8 spanned by the wave vectors (k i ) m i=1 .We also consider the quantity Because of the explicit time dependence in the Hamiltonian, the partial contribution of a set of m harmonics to the time derivative of C γ along the flow of

) and one has Ċ
Dynamical quantities C γ or C γ that are unaffected by a given set of leading harmonics, that is, with null partial contribution in Eq. ( 15) or (17), are denoted as quasiintegrals of motion from now on.More specifically, we build our analysis on ranking R 1 , since the resonant harmonics are those responsible for changes that accumulate stochastically over long timescales, driving chaotic diffusion.
In the framework of the aforementioned considerations, the resonances listed in Table II possess three different symmetries.
a. First symmetry.The rotational invariance of the entire Solar System implies the d'Alembert rule [21,24,45,46].Moreover, the Jupiterdominated eccentricity mode g 5 is the only fundamental Fourier mode of the outer planet forcing to appear in Table II.The quantity with b. Second symmetry.We write the eccentricity and inclination parts of the harmonic wave vectors explicitly, that is, k = (k ecc , k inc ) with k ecc , k inc ∈ R 4 .One can visually check that the harmonics in Table II verify the relation ).Therefore, denoting γ 1 = (0 4 , 1 4 ), the quantity is conserved by these resonances.C inc is the angular momentum deficit (AMD) [47] contained in the inclination d.o.f.This symmetry can then be interpreted as a remnant of the conservation of the AMD of the entire (secular) Solar System.We remark that the AMD contained in the eccentricity d.o.f., C ecc = 4 i=1 X i , is not invariant under the leading resonances because of the eccentricity forcing mainly exerted by Jupiter through the mode g 5 .The conservation of C inc depends on two facts: the inclination modes s 6 , s 7 , s 8 of the external forcing do not appear in Table II; low-order harmonics like 2g 1 −s 1 −s 2 , 2g 1 − 2s 1 , and 2g 1 − 2s 2 are never resonant (even if they can raise large quasi-periodic contributions), so that two AMD reservoirs C ecc and C inc are decoupled in Table II.We recall that the absence of an inclination mode s 5 in the external forcing relates to the fixed direction of the angular momentum of the entire Solar System [2,21,46].
c. Third symmetry.The first two symmetries could be expected to some extent on the basis of physical intuition of the interaction between outer and inner planets.However, it is not easy to even visually guess the third one from Table II.Consider the 30 × 8 matrix K 30 whose rows are the wave vectors (k i ) 30  i=1 of the listed resonances.A singular value decomposition shows that the rank of K 30 is equal to 6. Therefore, the linear subspace V 30 = span(k 1 , k 2 , . . ., k 30 ) spanned by the wave vectors has dimension 6.A Gram-Schmidt orthogonalization allows us to determine two linearly independent vectors that span its orthogonal complement Since the second symmetry clearly requires that γ 1 ∈ V ⊥ 30 , the three quantities C inc , C 2 , C ⊥ 2 are not independent and one has indeed The additional symmetry can thus be interpreted in terms of a certain decoupling between the d.o.f. 3, 4 and 1, 2, representing in the proper modes the Earth-Mars and Mercury-Venus subsystems, respectively.
The aforementioned symmetries, that exactly characterize the resonances listed in We remark that, differently from C inc and C 2 , the quantity E 2n is a non-linear function of the action-angle variables.However, as far as stable orbital evolutions are concerned, the convergence of the series expansion of the Hamiltonian is sufficiently fast that the linear LL approximation E 2 = H 2 +g 5 1 8 •I = C γ3 , with γ 3 = −ω LL +g 5 1 8 , reproduces reasonably well E 2n along the flow of H 2n for n > 1.The vector γ 3 is used in Sect.V, together with γ 1 and γ 2 , to deal with the geometry of the linear action subspace spanned by the QIs.The explicit expressions of these vectors are given in Appendix B. We mention that, differently from γ 1 and γ 2 , the components of γ 3 are not integer and they have the dimension of a frequency.

B. Slow variables
The QIs of motion E 2n , C inc , C 2 are clearly strong candidates for slow variables once evaluated along the orbital solutions.In what follows, to assess the slowness of a dynamical quantity when compared to the typical variations of the action variables, we consider the variance of its time series along a numerical solution.
We define the dimensionless QIs where C 0 stands for the current total AMD of the inner planets, that is, the value of C ecc + C inc at time zero.We stress that, by introducing the unit vectors At degree 2, one also has E 2 = C γ3 /C 0 .We then consider the ensembles of numerical integrations of H 4 and H 6 , with very close initial conditions and spanning 100 Gyr in the future, that have been presented in Ref. [7].
The top row of Fig. 2 shows the time evolution over 5 Gyr of the dimensionless QIs and of two components of the dimensionless action vector I = I/C 0 along the nominal orbital solutions of the two ensembles.We subtract from each time series its mean over the plotted time span.The time series are low-pass filtered by employing the Kolmogorov-Zurbenko (KZ) filter with three iterations of the moving average [4,48].A cutoff frequency of 1 Myr −1 is chosen to highlight the long-term diffusion that can be hidden by short-time quasi-periodic oscillations.This is in line with our definition of quasi-integrals based on contribution from resonant harmonics only.Figure 2 clearly shows that the QIs are slowly-diffusing variables when compared to an arbitrary function of the action variables.The behavior of the QIs along the nominal orbital solutions of Fig. 2 is confirmed by a statistical analysis in Appendix C. Figure 10 shows the time evolution of the distributions of the same quantities as Fig. 2 over the stable orbital solutions of the entire ensembles of 1080 numerical integrations of Ref. [7]. Figure 11 details the growth of the QI dispersion over time.We remark that C 2 and E 2n show very similar time evolutions along stable orbital solutions, as can be seen in the top row Fig. 2.This is explained by the interesting observation that the components of the unit vectors γ 2 and γ 3 differ from each other by only a few percent, as shown in Appendix B. However, we stress that the two vectors are in fact linearly independent: C 2 does not depend on the actions X 1 and X 2 , while E 2n does.The two QIs move away from each other when high eccentricities of Mercury are reached, that is, for large excursions of the Mercury-dominated action X 1 .

C. Weak resonances and Lyapunov spectrum
A fundamental result from Table II is that the symmetries introduced in Sect.IV A are still preserved by resonances that have half-widths an order of magnitude smaller than those of the strongest terms.It is natural to extract from ranking R 1 the weak resonances that break the three symmetries.A new ranking of reso-nances R 2 is defined in this way.Table III II, only harmonics that are resonant for more than 1% of the 5-Gyr time span of the nominal solution of Gauss's dynamics are shown.The leading symmetry-breaking resonances have half-widths of about 0.01 yr −1 .For each QI, the dominant contribution comes from harmonics involving Fourier modes of the outer planet forcing other than g 5 : the Saturn-dominated modes g 6 , s 6 and the modes g 7 , s 7 mainly associated to Uranus.In the case of C inc , there is also a contribution that starts at about 0.006 yr −1 with F 8 = 4g 1 − g 2 − g 3 − s 1 − 2s 2 + s 4 and comes from high-order internal resonances, that is, resonances that involve only the d.o.f. of the inner planets.We remark that the decrease of the resonance half-width with the index of the harmonic in Table III is steeper for C inc than for E 2n , C 2 , and is accompanied by a greater presence of high-order resonances.This may notably explain why the secular variations of C inc are somewhat smaller in the top row of Fig. 2. We finally point out the important symmetry-breaking role of the modes g 7 , s 7 , representing the forcing mainly exerted by Uranus.Differently from what one might suppose, these modes cannot be completely neglected when addressing the long-term diffusion of ISS.This recalls the role of the modes s 7 and s 8 in the spin dynamics of Venus [49], and is basically a manifestation of the long-range nature of the gravitational interaction.
As we state in Sect.III, a pair of Lyapunov exponents would vanish if there were an exact integral of motion.In the presence of a weakly broken symmetry, one may expect a small positive Lyapunov exponent whose value relates to the half-width of the strongest resonances driving the time variation of the corresponding QI.Such an argument is a natural extension of the correspondence between the FT-MLE and the top of the resonance spectrum given in Eq. ( 12).Comparison of Table III with the Lyapunov spectrum in Fig. 1a shows that the time statistics of the half-widths of the symmetry-breaking resonances of ranking R 2 overlaps with the ensemble distribution of the three smallest FT-LCEs, that is, λ 6 , λ 7 , λ 8 .One can indeed write where ∆ω R2 stands for the half-width of the uppermost resonances of ranking R 2 .a Only harmonics that are resonant for more than one percent of time are shown, i.e., τ res > 1%. exponents: Equation ( 23) is not a one-to-one correspondence, nor should it be understood as an exact relation since, for example, λ 6 is not well separated from the larger exponents.Its physical meaning is that the QIs are among the slowest d.o.f. of the ISS dynamics.Such a claim is one of the core points of this work.In Sect.V, we show its statistical validity in the geometric framework established by a principal component analysis of the orbital solutions.Moreover, Sect.IV D shows that Eq. ( 23) can be stated more precisely in the case of a simplified dynamics that underlies H 2n .We remark that E 2n , C inc , C 2 constitute a set of three QIs that are independent and nearly in involution, and it is thus meaningful to associate three different Lyapunov exponents with them.On the one hand, the independence is easily checked at degree 2 as the vectors γ 1 , γ 2 , γ 3 are linearly independent.

D. New truncation of the Hamiltonian
The fundamental role of the external modes g 6 , g 7 , s 6 , s 7 in Table III IV and related to C 2 are resonant for a few timesteps and their time statistics is very tentative.More precise estimations of the half-widths should be obtained over an ensemble of different orbital solutions, possibly spanning more than 5 Gyr.In any case, the fundamental point here is the drastic reduction in the size of the uppermost resonances with respect to Table III, and this is a robust result.We remark that resonances of order 12 and higher may also carry an important contribution at these scales, but they are excluded by the truncation at degree 10 adopted in Ref. [19] to establish the resonant harmonics, so they do not appear in the tables of this work.
Hamiltonian H • 2n .The implications of Table IV suggest to introduce an additional truncation in the Hamiltonian H 2n .This consists in dropping the harmonics of Eq. ( 6) that involve external modes other than g 5 : where φ 1 (t) = −g 5 t and • = ( 1 , 0, . . ., 0), with 1 ∈ Z. Consistently with the absence of symmetry-breaking resonances related to E 2n in Table IV, the corresponding dynamics admits the exact integral of motion which represents the transformed Hamiltonian under the canonical change of variables that eliminates the explicit time dependence in Eq. ( 24).We point out that, as the additional truncation is applied to the action-angle formulation of Eq. ( 6), the external modes other than g 5 still enter the definition of the proper modes of the forced Laplace-Lagrange dynamics [4].The orbital solution arising from H • 2n is initially very close to that of H 2n .A frequency analysis over the first 20 Myr shows that the differences in the fundamental frequencies of the motion between H • 2n and H 2n are of the order of 10 −3 arcsec yr −1 , an order of magnitude smaller than the typical frequency differences between H 4 and H 6 [4, Table 3].Therefore, even though H • 2n constitutes a simplification of H 2n , it should not be regarded as a toy model.Its dynamics, in particular, still possesses 8 d.o.f.We compute the Lyapunov spectrum of the Hamiltonian H • 4 in the same way as described in Sect.III in the case of H 2n .Since its dynamics turns out to be much more stable than that of H 4 (see Sect.VI, Fig. 7), we extend the computation to a time span of 100 Gyr.The marginal ensemble PDFs of the positive FT-LCEs are shown in Fig 1b .Comparing to the Lyapunov spectrum of H 4 , one notices that the distributions of the leading exponents turn out to be quite similar, apart from being more spaced and except for a slight decrease in their median values.However, such a decrease is more pronounced for smaller exponents, and the drop in the smallest exponents is drastic.The smallest one, λ 8 , decreases monotonically, consistently with the fact that E • 4 from Eq. ( 25) is an exact integral of motion.The exponent λ 7 drops by more than an order of magnitude, and apparently begins to stabilize around a few 10 −4 arcsec yr −1 , while λ 6 also reduces significantly, by a factor of three, to about 0.005 yr −1 .The drop in the smallest exponents agrees remarkably well with that of the half-width of the leading symmetry-breaking resonances when switching from Table III to Table IV.One can indeed write where ∆ω R3,Q stands for the half-width of the uppermost resonances of ranking R 3 related to the quasi-integral Q.The hierarchy of the three smallest exponents in the spectrum of Fig. 1b These one-to-one correspondences are a particular case of Eq. ( 23) and support the physical intuition behind it.In Sect.V, we prove the validity of Eq. ( 27) in the geometric framework established by a principal component analysis of the orbital solutions of H • 2n .Numerical integrations.We compute ensembles of 1080 orbital solutions of the dynamical models H • 4 and H • 6 , with initial conditions very close to the nominal ones of Gauss's dynamics and spanning 100 Gyr in the future.This closely follows what we did in Ref. [7] in the case of the models H 2n .The bottom row of Fig. 2 shows the filtered dimensionless QIs along the nominal solutions of the two models over the first 5 Gyr.The hierarchy of the QIs stated in Eq. ( 27) is manifest.The quantity C 2 has secular variations much slower than C inc , while the latter is itself slower with respect to its counterpart in the orbital solutions of H 2n .We remark that, as E • 2n is an exact integral of motion for the model H • 2n , we do not plot it.From Fig. 2 it is also evident how difficult can be the retrieval of the short-lasting resonances affecting C 2 from a solution of H • 2n spanning only a few billion years.The hierarchy of the QIs is confirmed by a statistical analysis in Appendix C. Figure 10 shows the entire time evolution of the distributions of the filtered dimensionless QIs over the stable orbital solutions of the ensembles of 1080 numerical integrations.Figure 11 details the growth of the QI dispersion over time.As suggested by Table IV, the drop in the diffusion rates of the QIs when switching from H 2n to H • 2n is manifest.

V. STATISTICAL DETECTION OF SLOW VARIABLES
Section IV shows how the slow-fast nature of the ISS dynamics, indicated by the Lyapunov spectrum, emerges from the quasi-symmetries of the resonant harmonics of the Hamiltonian.QIs of motion can be introduced semianalytically and they constitute slow quantities when evaluated along stable orbital solutions.In this section, we consider the slow variables that can be systematically retrieved from a numerically integrated orbital solution by means of a statistical technique, the principal component analysis.We show that, in the case of the forced secular ISS, the slowest variables are remarkably close to the QIs, and this can be established in a precise geometric framework.

A. Principal component analysis
PCA is a widely used classical technique for multivariate analysis [50,51].For a given dataset, PCA aims to find an orthogonal linear transformation of the variables such that the new coordinates offer a more condensed and representative view of the data.The new variables are called principal components (PCs).They are uncorrelated and ordered according to decreasing variance: the first PC and last one have, respectively, the largest and the smallest variance of any linear combination of the original variables.While one is typically interested in the PCs of largest variance, in this work we employ the variance of the time series of a dynamical quantity to assess its slowness when compared to the typical variations of the action variables (see Sect.IV B).We thus perform a PCA of the action variables I and focus on the last PCs, as they give a pertinent statistical definition of slow variables.We stress that, when coupled to a lowpass filtering of the time series, the statistical variance provides a measure of chaotic diffusion.
Implementation.Our procedure for the PCA is described briefly as follows [for general details see, e.g., 52,53].Let I(t) = (X(t), Ψ(t)) be the 8-dimensional time series of the action variables evaluated along a numerical solution of the equations of motion.As in Sect.IV B, we apply the KZ low-pass filter with three iterations of the moving average and a cutoff frequency of 1 Myr −1 to obtain the filtered time series Î(t) [4,48].In this way, the short-term quasi-periodic oscillations are mostly suppressed, which better reveals the chaotic diffusion over longer timescales.We finally define the meansubtracted filtered action variables over the time interval [t 0 , t 0 + T ] as Ĩ(t) = Î(t) − n −1 n−1 i=0 Î(t 0 + i∆t), where the mean is estimated by discretization of the time series with a sampling step ∆t such that T = (n − 1)∆t.The discretized time series over the given interval is stored in an 8 × n matrix: The PCA of the data matrix D consists in a linear transformation P = A T D, where A is an 8 × 8 orthogonal matrix (i.e.A −1 = A T ) defined as follows.By writing A = [a 1 , . . ., a 8 ], the column vectors a i ∈ R 8 are chosen to be the normalized eigenvectors of the sample covariance matrix, in order of decreasing eigenvalues: (n − 1) −1 DD T = AΣA T , where Σ = diag(σ 1 , . . ., σ 8 ) and σ 1 ≥ • • • ≥ σ 8 .The PCs are defined as the new variables after the transformation, that is, PC i = a i • I with i ∈ {1, . . ., 8}.The uncorrelatedness and the ordering of the PCs can be easily seen from the diagonal form of their sample covariance matrix, (n − 1) −1 PP T = Σ, from which it follows that the variance of PC i is σ i .
Among all the linear combinations in the action variables I, the last PC, i.e., PC 8 , has the smallest variance over the time interval [t 0 , t 0 + T ] of a given orbital solution.The second last PC, i.e., PC 7 , has the second smallest variance and is uncorrelated with PC 8 , and so on.It follows that the linear subspace spanned by the last k PCs is the k-dimensional subspace of minimum variance: the variance of the sample projection onto this subspace is the minimum among all the subspaces of the same dimension.These properties indicate that the last PCs provide a pertinent statistical definition of slow variables along numerically integrated solutions of a dynamical system.The linear structure of the PCA, in particular, seems adapted to quasi-integrable systems close to a quadratic Hamiltonian, like the ISS.In such a case, one may reasonably expect that the slow variables are, to a first approximation, linear combinations of the action variables.We remark that the mutual orthogonality allows us to associate a linear d.o.f. to each PC.
Aggregated sample.Instead of considering a specific solution, it is also possible to take the same time interval from m different solutions, and stack them together to form an aggregated sample: where D i is the data matrix of Eq. ( 28) for the ith solution.Since this work deals with a non-stationary dynamics, as the ISS ceaselessly diffuses in the phase space [7], we always consider the same time interval for each of the m solutions.The aggregated sample is useful in capturing globally the behavior of the dynamics, because it averages out temporary and rare episodes arising along specific solutions.

B. Principal components and quasi-integrals
Both the QIs and the last PCs represent slow variables, but are established through two different methods.Equations ( 23) and ( 27) claim that the QIs found semianalytically in Sect.IV are among the slowest d.o.f. of the ISS dynamics.This naturally suggests to compare the three QIs with the three last PCs retrieved from numerically integrated orbital solutions.In this part, we first introduce the procedure that we implement to establish a consistent and systematic correspondence between QIs and PCs.We then present both a visual and a quantitative geometric comparison between them.

Tweaking the QIs
The three last components PC 8 , PC 7 , PC 6 are represented by the set of vectors S PCs = {a 8 , a 7 , a 6 } belonging to R 8 .By construction, these PCs have a linear, hierarchical, and orthogonal structure.In other words: the PCs are linear combinations of the action variables I; denoting by the order of statistical variance, one has PC 8 PC 7 PC 6 ; the unit vectors (a i ) 8 i=6 are orthogonal to each other.On the other hand, the QIs of motion C inc , C 2 , E 2n do not possess these properties.Therefore, we adjust them in such a way to reproduce the same structure.
a. Linearity.While C inc and C 2 are linear functions of the action variables, E 2n is not when n > 1.Nevertheless, as we explain in Sect.IV A, as far as one considers stable orbital solutions, the linear LL approximation E 2 = γ 3 •I reproduces E 2n reasonably well.Therefore, we consider the three linear QIs of motion C inc , C 2 , E 2 , which are respectively represented by the set of R 8 -vectors S QIs = {γ 1 , γ 2 , γ 3 }.In this way, the 3-dimensional linear subspaces of the action space spanned by the sets S QIs and S PCs can be compared.b.Ordering.We define a set of QIs that are ordered by statistical variance, as it is the case for the PCs.We follow two different approaches according to model H • 2n in Eq. ( 24) or H 2n in Eq. ( 6) (clearly n > 1).c. Orthogonality.We apply the Gram-Schmidt process to the ordered set S QIs to obtain the orthonormal basis S QIs = {α 1 , α 2 , α 3 }.The set S QIs clearly spans the same subspace as S QIs .Moreover, the Gram-Schmidt process preserves the hierarchical structure, that is, the two m-dimensional subspaces spanned by the first m ≤ 3 vectors of S QIs and S QIs , respectively, are identical.
In the end, we obtain a linear, ordered, and orthogonal set of modified QIs of motion {QI 1 , QI 2 , QI 3 }, where QI i = α i • I.

Visual comparison
We now visually compare the vectors α 1,2,3 of the modified QIs with the corresponding vectors a 8,7,6 of the last three PCs.We use the ensembles of 1080 numerically integrated orbital solutions of the models H 4 and H 4 considered in Sects.IV B and IV D, respectively.The nominal solution of each set is denoted as sol.#1 from now on.For the model H 4 , we also consider two other solutions: sol.#2 that represents a typical evolution among the 1080 solutions, and sol.#3 representing a rarer one.The particular choice of these two solutions is detailed in Sect.V B 3.
Hamiltonian H • 4 .The modified QIs can be explicitly derived in this case and comprise interpretable physical Here, QI 1 is proportional to C2 and QI 2 is proportional to C ⊥ 2 ; see Eq. (20).
quantities.One has QI 1 proportional to C 2 and QI 2 proportional to C ⊥ 2 .Moreover, QI 3 is the component of E 2 that is orthogonal to both C 2 and C ⊥ 2 .Figure 3 shows the comparison between the modified QIs and the corresponding PCs for three different time intervals along sol.#1 of H • 4 (see Fig. 2 bottom left for its time evolution).The agreement of the pairs (QI 1 , PC 8 ), (QI 2 , PC 7 ), and (QI 3 , PC 6 ) across different intervals is manifest and even impressive.One can appreciate that the "slower" the PC, the more similar it is to its corresponding QI.The overlap between the modified QIs and the three last PCs means that the QIs of motion span the slowest 3dimensional linear subspace of the action space.Therefore, to a linear approximation, they represent the three slowest d.o.f. of the H • 4 dynamics.The quasi-integral C 2 represents the slowest linear d.o.f.: it coincides with the last principal component PC 8 , which has the smallest variance among all the linear combinations of the action variables.C inc and E 2 represent the second and the third slowest linear d.o.f., respectively: the component of C inc orthogonal to C 2 , i.e., C ⊥ 2 , matches the second last principal component PC 7 ; the component of E 2 orthogonal to the subspace generated by (C 2 , C inc ) matches the third last principal component PC 6 .The strong hierarchical structure of the slow variables for the simplified dynamics H • 4 is clearly confirmed by the almost frozen basis vectors of the PCs.
Hamiltonian H 4 .In this case, the QIs of motion C inc , C 2 , E 2 do not show a clear hierarchical structure in terms of statistical variance.Therefore, we consider the whole subspace spanned by the three QIs with respect to that spanned by the three last PCs.Since it is not easy to visually compare two 3-dimensional subspaces of R 8 , we compare their basis vectors instead.The basis α 1,2,3 of modified quasi-integrals QI 1,2,3 is computed according to the algorithm presented in Sect.V B 1.
Figure 4 presents the comparison between the modified QIs and the corresponding PCs across three different time intervals of three solutions of H 4 (see Fig. 5 for their time evolution).The first two, sols.#1 and #2, show thorough agreement between the pairs of QIs and PCs across all intervals, which indicates close proximity between the two subspaces V QIs = span(S QIs ) and V PCs = span(S PCs ).One can appreciate that the directions of the basis vectors are quite stable.The last component PC 8 , in particular, remains close to C inc .The slowest linear d.o.f. of H 4 can thus be deduced to be close to C inc , in line with the discussion in Sect.IV C. Such a result shows how interesting physical insight can be gained through the PCA.Some changes in the basis vectors can arise, however, as for the first time interval of sol.#2.This may be expected from a dynamical point of view.Differently from H • 4 , there is no pronounced separation between the slowest d.o.f. at the bottom of the Lyapunov spectrum in Fig. 1a: the marginal distributions of consecutive exponents can indeed touch or overlap each other.Therefore, the hierarchy of slow variables is not as frozen as in H • 4 and it can change along a given orbital solution.
Solutions #1 and #2 represent typical orbital evolutions.If the same time intervals of all the 1080 solutions are stacked together to form an aggregated sample on which the PCA is applied, the features mentioned above persist: the agreement between QIs and PCs, the stability of the basis vectors, and the similarity between PC 8 and C inc (see Fig. 4).Once again, the PCA confirms that the subspace spanned by the three QIs is overall close to the slowest 3-dimensional linear subspace of the action space.Therefore, to a linear approximation, they represent the three slowest d.o.f. of the H 4 dynamics.We remark that the slowness of the 3-dimensional subspace spanned by the QIs is a much stronger constraint than the observation that each QI is a slow variable.To give an example, let Q = q • I be a slow variable with unit vector q.If is an arbitrary small vector, i.e.
1, then Q = ( q + ) • I can also be considered as a slow variable, whereas the normalized difference of two quantities, •I, is generally not.Therefore, the linear subspace spanned by Q and Q , that is, by q and , is not a slow 2-dimensional subspace.
Solution #3 in Fig. 4 represents an edge case (see Fig.  sol.#1 Time evolution over 5 Gyr of the dimensionless QIs of motions ( Cinc, C2, E) and of two representatives of the dimensionless action variables ( X1, Ψ3) for three solutions of H4, that is, sol.#1 (top), sol.#2 (middle), and sol.#3 (bottom).E stands for E4.The time series are low-pass filtered with a cutoff frequency of 1 Myr −1 and the mean over 5 Gyr is subtracted.
for its time evolution).Typically, the variances of the QIs are at least one order of magnitude smaller than those of the action variables, which allows a clear separation.Nevertheless, the distinction between the QIs and faster d.o.f.can be more difficult in two rare possibilities.Firstly, if the change in a QI accumulates continually in one direction, its variance can inflate over a long time interval.This is the case for the interval [0, 5] Gyr of sol.#3.Secondly, the variance of a variable that is typically fast can suddenly dwindle during a certain period of time, for example, Ψ 3 over the interval [1,2] Gyr of sol.#3.In both cases, the slow subspace defined by the three last PCs can move away from the QI subspace due to the contamination by d.o.f. that are typically faster.This is reflected in the mismatch of QI 3 and PC 6 on the last two time intervals of sol.#3 in Fig. 4. We remark that PC 8,7 are still relatively close to QI 1,2 , which indicates that the slowest 2-dimensional subspace spanned by PC 8,7 still resides inside the QI subspace.It should be stressed that this disagreement between QIs and PCs does not mean that the QIs are not slow variables in this case.The mismatch has a clear dynamical origin instead.The resonance tables of this work have been retrieved from a single, very long orbital solution, with the idea that its time statistics is representative of the ensemble statistics over a set of initially very close solutions [19].Therefore, the QIs derived from these tables characterize the dynamics in a global sense.The network of resonances can temporarily change in an appreciable way along specific solution, or be very particular along rare orbital solutions.In these cases, a mismatch between the last PCs and the present QIs may naturally arise.Moreover, the contamination of the QIs by d.o.f. that are typically faster may also be expected from the previously mentioned lack of a strong hierarchical structure of the slow variables.The Lyapunov spectrum in Fig. 1a shows that the marginal distributions of the exponents λ 5 and λ 6 , for example, are not separate but overlap each other.

Distance between the subspaces of PCs and QIs
The closeness of the two 3-dimensional linear subspaces V PCs , V QIs ⊂ R 8 spanned by the sets of vectors S PCs and S QIs , respectively, can be quantitatively measured in terms of a geometric distance.This can be formulated using the principal (canonical) angles [55][56][57].
Let A and B be two sets of m ≤ n independent vectors in R n .The principal vectors (p k , q k ) m k=1 are defined recursively as solutions to the optimization problem: between the two subspaces span(A) and span(B) are then defined by The principal angle θ 1 is the smallest angle between all pairs of unit vectors in span(A) and span(B); the principal angle θ 2 is the smallest angle between all pairs of unit vectors that are orthogonal to the first pair; and so on.Given the matrices defining the two subspaces, the principal angles can be computed from the singular value decomposition of their correlation matrix.The result is the canonical correlation matrix diag(cos θ 1 , . . ., cos θ m ).This cosine-based method is often ill-conditioned for FIG. 6. PDF of the distance between two random 3dimensional linear subspaces of R 8 (blue, 10 5 draws) compared with the PDF of the distance between the two subspaces VPCs (PC8,7,6) and VQIs (QI1,2,3) arising from the time interval [0, 5] Gyr of 1080 solutions of H • 4 (top) and 10 800 solutions of H4 (bottom) (green).For each model, the subspace distance from the same time interval of representative solutions (vertical red lines) and of the aggregated sample of all the solutions (vertical dark green line) are indicated.The subspace distance is given by Eq. (31).
small angles.In such case, a sine-based algorithm can be employed [58].In this work, we use the combined technique detailed in Ref. [59].
Once the principal angles have been introduced, different metrics can be defined to measure the distance between two subspaces.In this work, we choose the normalized chordal distance [57]: The distance is null if A and B are the same subspace and equal to 1 when they are orthogonal.We use this metric to show that the subspace closeness suggested by Figs. 3 and 4 is indeed statistically significantly.More precisely, we provide evidence against the null hypothesis that the distribution of distances between V PCs and V PCs , arising from the H • 4 and H 4 dynamics, coincides with that of randomly drawn 3-dimensional subspaces of R 8 .The PDF of the distance between two random 3-dimensional subspaces of R 8 is shown in Fig. 6 in blue (such random subspaces can be easily generated by sampling sets of 3 vectors uniformly on the unit 7-sphere [60]).While the range of possible distances is [0,1], the distribution concentrates on the right-hand side of the interval, with a probability of approximately 99.3% that the distance is larger than 0.6.In this regard, we remark that the notion of distance in high-dimensional spaces is very different from our intuition in a 3-dimensional world.If we draw randomly two vectors in a very high-dimensional space, it is extremely likely that they will be close to mutual orthogonality.
The upper panel of Fig. 6 shows in green the PDF of the distance between V PCs and V QIs arising from the time interval [0, 5] Gyr of the 1080 orbital solutions of model In the lower panel, we consider a larger ensemble of 10 800 solutions of model H 4 spanning the same time interval [7], and plot the corresponding PDF of the distance between V PCs and V QIs .In both cases, the distance stemming from the aggregated sample of all the solutions is indicated by a vertical dark green line.We also report the distances from the specific solutions considered in Figs. 3 and 4 as vertical red lines.As the PDFs of both models peak at small distances, there results a strong evidence that the distribution of distances between the subspaces spanned by the PCs and the QIs is not that of random subspaces.In this sense, the closeness of the subspaces V PCs and V QIs is a statistically robust result.In the case of the simplified dynamics H • 4 , the PDF peaks around a median of roughly 0.08 and has small variance.Switching to model H 4 , the median increases to about 0.26 and the PDF is more spread out, with a long tail toward larger distances.The differences between the PDFs of the two models follow quite naturally the discussion in Sect.V B 2: a quasi frozen hierarchy of the slowest variables for H • 4 ; a larger variance for H 4 related to contamination by d.o.f. that are typically faster and to variations of the resonant network with respect to the nominal solution of Gauss's dynamics which is used to infer the QIs.Solution #3 in Fig. 4 represents in this sense an edge case of the distance distribution, while sol.#2 is a typical solution close to the PDF median.

VI. IMPLICATIONS ON LONG-TERM STABILITY
The existence of slow variables can have crucial implications on the stability of the ISS.The QIs of motion can effectively constraint in an adiabatic way the chaotic diffusion of the planet orbits over long timescales, forbidding in general a dynamical instability over a limited time span, e.g., several billions of years.Here we give compelling arguments for such a mechanism.
Figure 7 shows the cumulative distribution function (CDF) of the first time that Mercury eccentricity reaches a value of 0.7, from the ensembles of 1080 orbital solutions of H that such a high eccentricity is a precursor of the dynamical instability (i.e., close encounters, collisions, or ejections of planets) of the ISS [6].We also report the same CDF for the models H 4 and H 6 , which we recently computed in Ref. [7].One can appreciate that the time corresponding to a probability of instability of 1% is greater than 100 Gyr for the H 2n and H 2n relates to the smallest Lyapunov exponents (Fig. 1), and this is accompanied by a much slower diffusion of the QIs for H • 2n (Figs. 2, 10, and 11), Fig. 7 indicates that the dynamical half-life of the ISS is linked to the speed of diffusion of these slow quantities in a critical way.We stress that the slower diffusion toward the dynamical instability in the H • 2n model derives from neglecting the external forcing mainly exerted by Saturn, Uranus, and Neptune.
We also observe that, to a linear approximation, the knowledge of C inc and E 2 allows us to bound the variations of the action variables X, Ψ. Recalling that the actions are positive quantities, from Eq. ( 19) one sees that fixing a value of C inc puts an upper bound to the variations of the inclination actions Ψ.As a consequence, at degree 2 in eccentricities and inclinations, fixing a value of the QI with γ 3 = (γ ecc 3 , γ inc 3 ), also bounds the upper variations of the eccentricity actions X, since the components of γ ecc 3 have all the same sign, as those of γ inc 3 (see Appendix B).This is an important point, as we state in Sect.I that the lack of any bound on the chaotic variations of the planet orbits is one of the reasons that complicate the understanding of their long-term stability.We remark that the secular planetary phase space can be bound by fixing the value of the total AMD, that is, C ecc +C inc [47].A statistical study of the density of states that are a priori accessible can then be realized [61].It is not, however, fully satisfying to consider a fixed value of total AMD of the ISS, as we show that C ecc is changed by some of the leading resonances of the Hamiltonian, as a result of the eccentricity forcing mainly exerted by Jupiter through the mode g 5 .Moreover, the destabilization of the ISS consists indeed in a large transfer of eccentricity AMD, C ecc , from the outer system to the inner planets through the resonance g 1 − g 5 [5,6,36,62].It should be noted that C ecc can still be consider as a slow quantity with respect to an arbitrary function of the action variables, as it is only changed by the subset of the leading resonances involving the external mode g 5 .This slowness has indeed been observed on stable orbital solutions of the Solar System [47] and supports the statistical hypothesis in Ref. [61] that allows one to obtain a very reasonable first guess of the long-term PDFs of the eccentricities and inclinations of the inner planets.
The emerging picture explains the statistical stability of the ISS over billions of years in a physically intuitive way.The chaotic behavior of the planet orbits arises from the interaction of a number of leading resonant harmonics of the Hamiltonian, which determine the Lyapunov time.The strongest resonances are characterized by some exact symmetries, which are only broken by weak resonant interactions.These quasi-symmetries naturally give birth to QIs of motion, quantities that diffuse much more slowly than the LL action variables, constraining the variations of the orbits.The long dynamical half-life of the ISS is connected to the speed of this diffusion, which eventually drives the system to the instability.It should be stressed that, besides the speed of diffusion, the lifetime of the inner orbits also depends on the initial distance of the system from the instability boundary defined by the resonance g 1 − g 5 .This geometric aspect includes the stabilizing role of general relativity [5,6], which moves the system away from the instability boundary by 0.43 yr −1 , and the destabilizing effect of terms of degree 6 in eccentricities and inclinations of the planets [7].

VII. DISCUSSION
This work introduces a framework that naturally justifies the statistical stability shown by the ISS over a timescale comparable to its age.Considering a forced secular model of the inner planet orbits, the computation of the Lyapunov spectrum indicates the existence of very different dynamical timescales.Using the computer algebra system TRIP, we systematically analyze the Fourier harmonics of the Hamiltonian that become resonant along a numerically integrated orbital solution spanning 5 Gyr.We uncover three symmetries that characterize the strongest resonances and that are broken by weak resonant interactions.These quasi-symmetries generate three QIs of motion that represent slow variables of the secular dynamics.The size of the leading symmetrybreaking resonances suggests that the QIs are related to the smallest Lyapunov exponents.The claim that the QIs are among the slowest d.o.f. of the dynamics constitutes the central point of this work.On the one hand, it is supported by the analysis of the underlying Hamiltonian H • 2n , in which one neglects the forcing mainly exerted by Saturn, Uranus, and Neptune and, as a consequence, the diffusion of the QIs is greatly reduced.On the other hand, the geometric framework established by the PCA of the orbital solutions independently confirms that the QIs are statistically the slowest linear variables of the dynamics.We give strong evidence that the QIs of motion play a critical role in the statistical stability of the ISS over the Solar System lifetime, by adiabatically constraining the long-term chaotic diffusion of the orbits.

A. Inner Solar System among classical quasi-integrable systems
It is valuable to contextualize the dynamics of the ISS in the class of classical quasi-integrable systems.A comparison with the Fermi-Pasta-Ulam-Tsingou problem, in particular, deserves to be made.This concerns the dynamics of a one-dimensional chain of identical masses coupled by nonlinear springs.For weak nonlinearity, the normal modes of oscillation remain far from the energy equipartition expected from statistical mechanics for a very long time [13].One way to explain the lack of energy equipartition reported by Fermi and collaborators is through the closeness of the FPUT problem to the integrable Toda dynamics [63][64][65].This translates in a very slow thermalization of the action variables of the Toda problem and of the corresponding integrals of motion along the FPUT flow [15,[65][66][67][68][69].In the framework of the present study, the very long dynamical half-life of the ISS is also likely to be the result of the slow diffusion of some dynamical quantities, the QIs of motion.We find, in particular, an underlying Hamiltonian H • 2n for which this diffusion is greatly reduced, as a consequence of neglecting the forcing mainly exerted by Saturn, Uranus, and Neptune.This results in a dynamics that can be considered as stable in an astronomical sense.We stress that, differently from the FPUT problem, H • 2n is not integrable as the Toda Hamiltonian.It is indeed chaotic and shares with the original Hamiltonian H 2n the leading Lyapunov exponents.The QIs that we find in this work are only a small number of functions of the action-angle variables of the integrable LL dynamics, and are related to the smallest Lyapunov exponents of the dynamics.Our study suggests that in the FPUT problem the very slow thermalization occurring beyond the Lyapunov time might be understood in terms of combinations of the Toda integrals of motion diffusing over very different timescales.
The long-term diffusion in chaotic quasi-integrable systems should be generally characterized by a broad range of timescales that results from the progressive, hierarchical breaking of the symmetries of the underlying integrable problem by resonant interactions [70][71][72].A hierarchy of Lyapunov exponents spanning several orders of magnitude, in particular, should be common among this class of systems [e.g., 73].

B. Methods
The long-term dynamics of the ISS is described by a moderate but not small number of d.o.f., which places it far from the typical application fields of celestial mechanics and statistical physics.The first discipline often studies dynamical models with very few degrees of freedom, while the second one deals with the limit of a very large number of bodies.Chaos also requires a statistical description of the inner planet orbits.But the lack of a statistical equilibrium, resulting from a slow but ceaseless diffusion of the system, places the ISS outside the standard framework of ergodic theory.The kind of approach we develop in this work is heavily based on computer algebra, in terms of systematic series expansion of the Hamiltonian, manipulation of the truncated equations of motion, extraction of given Fourier harmonics, retrieval of polynomial roots, etc. [4,19].This allows us to introduce QIs of motion in a 16-dimensional dynamics by analyzing how action-space symmetries are progressively broken by resonant interactions.Our effective method based on the time statistics of resonances arising along a single, very long numerical integration is alternative to formal approaches that define QIs via series expansions [e.g.74,75].The practical usefulness of these formal expansions for a dynamics that covers an intricate, high-dimensional network of resonances seems indeed doubtful.Through the retrieval of the half-widths of the symmetry-breaking resonances, computer algebra also permits us to extend the correspondence between the Lyapunov spectrum and the spectrum of resonances well beyond the standard relation linking the Lyapunov time to the strongest resonances [76].
In the context of dynamical systems with a number of d.o.f. that is not small, this work also considers an approach based on PCA.The role of this statistical technique can be twofold.We use PCA as an independent test to systematically validate the slowness of the QIs.While being introduced semi-analytically as dynamical quantities that are not affected by the leading resonances, they can indeed be related to the last PCs.By extension, the first PCs should probe the directions of the main resonances.This leads to a second potential application of the PCA, which should offer a way to retrieve the principal resonant structure of a dynamical system.In this sense, PCA represents a tool to systematically probe numerical integrations of a complex dynamics and distill important hidden insights.We emphasize that PCA is the most basic linear technique of dimensionality reduction and belongs to the more general class of the unsupervised learning algorithms.There are more sophisticated methods of feature extraction that can be more robust [e.g.77,78] and can incorporate nonlinearity [79].These methods are often less intuitive to understand, less straightforward to apply, and harder to interpret than PCA.Yet, they might be more effective and worth pursuing for future works.
With long-term numerical integration and a computer algebra system at one's disposal, the entire strategy we develop in this work can in principle be applied to other planetary systems and quasi-integrable Hamiltonian dynamics with a moderate number of d.o.f.To quantitatively estimate the numerical precision on the computed FT-LCEs, we exploit the symmetry of the spectrum stated in Eq. (9).For a single orbital solution, the relative numerical error on each exponent λ i can be estimated as We plot in Fig. 8b the medians of i for the ensemble of 150 orbital solutions of Fig. 1a.The relative errors decrease asymptotically with time, as expected.Even in the case of the smallest exponent, λ 8 , the median error is less than 10% at 10 Gyr.Hamiltonian H 6 .We compute for comparison the FT-LCEs of the forced ISS truncated at degree 6 in eccentricities and inclinations, that is, H 6 .We consider 150 stable orbital solutions with initial conditions very close to the nominal values of Gauss's dynamics and random sets of initial tangent vectors, as we do for the truncation at degree 4. Figure 9 shows the [5th, 95th] percentile range of the marginal PDF of each FT-LCE estimated from the ensemble of solutions.Apart from being somewhat larger, the asymptotic distributions of the exponents are very similar to those of H 4 shown in Fig. 1a.

Appendix C: Ensemble distributions of the quasi-integrals over time
To retrieve the long-term statistical behavior of the QIs, we consider the ensembles of 1080 numerical integrations of the dynamical models H 4 and H 6 , with very close initial conditions and spanning 100 Gyr in the future, that have been presented in Ref. [7].We also consider the similar ensembles of solutions for the simplified Hamiltonians H • 4 and H • 6 that we introduce in Sect.IV D. We report in Fig. 10 the time evolution of the ensemble PDFs of the low-pass filtered dimensionless QIs and dimensionless actions X 1 , Ψ 3 for the different models (the cutoff frequency of the time filter is set to 1 Myr −1 , as in Sec.IV B).More precisely, to highlight the growth of the statistical dispersion, we consider at each time the PDF of the signed deviation from the ensemble mean, so that all the plotted distributions have a null mean.At each time, the PDF estimation takes into account only the stable orbital solutions, that is, those solutions whose running maximum of Mercury eccentricity is smaller than 0.7 [7]. Figure 10 shows that the QIs are indeed slow quantities when compared to the LL action variables.The growth of the QI dispersion is detailed in Fig. 11, where we report the time evolution of the interquartile range (IQR) of their distributions.After a transient phase lasting about 100 Myr and characterized by the exponential separation of close trajectories, the time growth of the IQR follows a power law typical of diffusion processes.and its PDF has null dispersion.

8 (b) H • 4 FIG. 1 .
FIG. 1. Positive FT-LCEs λi of the forced secular ISS from Hamiltonians H4 (a) and H • 4 [Eq.(24)] (b), and corresponding characteristic timescales λ −1 i .The bands represent the [5th, 95th] percentile range of the marginal PDFs estimated from an ensemble of 150 stable orbital solutions with very close initial conditions.The lines denote the distribution medians.The H • 4
most resonances drops by two orders of magnitude.We stress that such harmonics are resonant for very short periods of time along the 5 Gyr spanned by the nominal solution of Gauss's dynamics.To retrieve the time statistics of the resonances affecting C 2 , we indeed choose to repeat the computations of Ref.[19] by increasing the cutoff frequency of the low-pass filter applied to time series of the action-angle variables from (5 Myr) −1 to 1 Myr −1 [19, Appendices F.2 and G.5].The filtered time series have then been resampled with a timestep of 50 kyr.Many harmonics we show in Table

5
FIG. 5.Time evolution over 5 Gyr of the dimensionless QIs of motions ( Cinc, C2, E) and of two representatives of the dimensionless action variables ( X1, Ψ3) for three solutions of H4, that is, sol.#1 (top), sol.#2 (middle), and sol.#3 (bottom).E stands for E4.The time series are low-pass filtered with a cutoff frequency of 1 Myr −1 and the mean over 5 Gyr is subtracted.

H 4 :
FIG. 8. (a) Positive FT-LCEs of Hamiltonian H4 and corresponding characteristic timescales for a single initial condition and an ensemble of 150 random sets of initial tangent vectors.The bands represent the [5th, 95th] percentile range of the marginal PDFs.The lines denote the distribution medians.(b) Medians of the relative numerical errors i on the FT-LCEs λi, as defined in Eq. (A1), for the ensemble of 150 orbital solution of Fig. 1a.

8 FIG. 9 .
FIG. 9. Positive FT-LCEs λi of Hamiltonian H6 and corresponding characteristic timescales λ −1 i .The bands represent the [5th, 95th] percentile range of the marginal PDFs estimated from an ensemble of 150 stable orbital solutions with very close initial conditions.The lines denote the distribution medians.

FIG. 10 .
FIG. 10.Time evolution over 100 Gyr of the PDF of the signed deviation from the mean of the low-pass filtered dimensionless QIs and dimensionless actions X1, Ψ3.Estimation from an ensemble of 1080 numerical orbital solutions for different models (H4, H6, H • 4 , and H • 6 ).First row : C inc .Second row : C2.Third row : E4 (H4) and E6 (H6).Fourth row : X1.Fifth row : Ψ3.The time of each curve is color coded.At each time, the estimation only takes into account stable solutions, that are those with a running maximum of Mercury eccentricity smaller than 0.7.The quantity E • 2n is an exact integral of motion for the model H • 2n

FIG. 11 .
FIG. 11.Time evolution of the interquartile range (IQR) of the ensemble PDFs of the QIs shown in Fig. 10.Left: C inc .Middle: C2.Right: E4 (H4) and E6 (H6).The quantity E • 2n is an exact integral of motion for the model H • 2n and its PDF has a null IQR.

TABLE I .
Summary of the different models of forced secular ISS considered in this work.Gauss's dynamics results from first-order averaging of the N -body Hamiltonian over the mean longitudes of the planets.The dynamics generated by H2n and H2n are practically equivalent and treated as such.The H • 2n models are introduced and discussed in Sec.IV D. ×Z 7 .At degree 2, one has H 2 = −ω LL •I, where ω LL = (g LL , s LL ) ∈ R 4 ×R 4 are the LL fundamental precession frequencies of the inner planet perihelia and nodes.Hamiltonian H 2n is in quasi-integrable form.
, with unchanged action variables, allows us to remove the explicit time dependence in these harmonics.Quantity E 2n coincides with the transformed Hamiltonian and the harmonics in TableIIdo not contribute to its time derivative.
is therefore unaffected by the resonances listed in TableII.In an equivalent way, the time-dependent canonical transformation θ → θ +g 5 t 1 8 Table II, naturally repre-sent quasi-symmetries when considering the entire spectrum of resonances R 1 .They are indeed broken at some point by weak resonances (see Sect.IV C).Quantities E 2n , C inc , and C 2 are the corresponding QIs of motion.The persistence of the three symmetries under the 30 leading resonances is somewhat surprising.Concerning C inc and C 2 , for example, one might reasonably expect that, since the ISS has 8 d.o.f., the subspace spanned by the wave vectors of just a dozen of harmonics should already have maximal dimension, destroying all possible symmetries.
• 2n is exactly conserved and not shown).The time series are low-pass filtered with a cutoff frequency of 1 Myr −1 and the mean over 5 Gyr is subtracted.The variations of the QIs are enlarged in the insets.The H • 2n models are introduced and discussed in Sec.IV D.
FIG.2.Time evolution over 5 Gyr of the dimensionless QIs ( C inc , C2, E) and of two representatives of the dimensionless action variables ( X1, Ψ3) along the nominal orbital solutions of different models.Top row : H4 and H6 ( E stands for E4 and E6, respectively).Bottom row : H • 4 and H • 6 from Eq. (24) ( E reports the 10 strongest symmetry-breaking resonances that change E 2n , C inc , C 2 , respectively.As in Table

TABLE III .
Top of ranking R2.First 10 symmetry-breaking resonances of H10 along the 5-Gyr nominal solution of Gauss's dynamics, that change E2n, C inc , and C2, respectively (see TableIIfor details).

TABLE IV .
Top of ranking R3.First 10 symmetry-breaking resonances of H10 along the 5-Gyr nominal solution of Gauss's dynamics, that only involve g5 among the external modes and change E2n, C inc , and C2, respectively.
raises the question of which symmetry-breaking resonances persist if one excludes all the Fourier harmonics that involve external modes other than g 5 .Therefore, we define a new ranking R 3 by extracting such resonances from ranking R 2 .TableIVreports the 10 strongest resonances per each broken symmetry.The difference with respect to Table III is manifest.As g 5 is the only external mode remaining, there are no resonances left that can contribute to the time evolution of E 2n .For the remaining two QIs, the only harmonics that appear in Table IV are of order 8 or higher, and this is accompanied by a significant drop in the half-width of the leading resonances.In the case of C inc , the half-width of the uppermost resonances is now around 0.005 yr −1 .One can appreciate that the activation times τ res of the resonances do not exceed a few percent, differently from Table III.The most impressive change is, however, related to C 2 : only harmonics of order 10 appear in Table IV and the half-width of the upper consistently follows that of the QIs suggested in Table IV by the very different sizes of the leading resonances.In other words, one can state: H • 2n : A strong hierarchy of statistical variances among the QIs emerges from the size of the leading symmetry-breaking resonances in Table IV and from the orbital solutions in Figs. 2, 10, and 11.One has E • 2n ≺ C 2 ≺ C inc .While E • 2n is an exact non-linear integral of motion, we expect that its linear truncation E • 2 = E 2 varies more than C 2 and C inc .Therefore, we consider the ordered set of QIs of motion {C 2 , C inc , E 2 } represented by the ordered set of vectors S QIs = {γ 2 , γ 1 , γ 3 }.H 2n : Since the leading resonances affecting the QIs in Table III have comparable sizes, there is no clear order of statistical variances that can be inferred.We then implement a systematic approach that orders the QIs by simply inheriting the ordering of the PCs.More precisely, we define a set of ordered vectors S QIs through the projections of the three last PCs onto the linear subspace generated by the QIs: S QIs = {proj SQIs (a 8 ), proj SQIs (a 7 ), proj SQIs (a 6 )}[54].As a result, the new set of QIs mirrors the hierarchical structure of the PCs.We stress that S QI spans the same subspace of R 8 as S QI , since the ordered QIs are just linear combinations of the original ones.
• 4 model, while it is about 15 Gyr for H 4 .At degree 6, this time still increases from 5 Gyr for H 6 to about 20 Gyr in H •