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

Federico Mogavero, Nam H. Hoang, Jacques Laskar

► To cite this version:

Federico Mogavero, Nam H. Hoang, Jacques Laskar. Timescales of Chaos in the Inner Solar System: Lyapunov Spectrum and Quasi-integrals of Motion. Physical Review X, 2023, 13, 10.1103/Phys-RevX.13.021018. insu-04849174

HAL Id: insu-04849174

https://insu.hal.science/insu-04849174v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Timescales of Chaos in the Inner Solar System: Lyapunov Spectrum and Quasi-integrals of Motion

Federico Mogavero<sup>*</sup>, Nam H. Hoang<sup>†</sup>, and Jacques Laskar<sup>‡</sup>

IMCCE, CNRS UMR 8028, Observatoire de Paris, Université PSL, Sorbonne Université,
77 Avenue Denfert-Rochereau, 75014 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 18 November 2022; accepted 3 March 2023; published 3 May 2023)

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

DOI: 10.1103/PhysRevX.13.021018

Subject Areas: Gravitation, Nonlinear Dynamics, Statistical Physics

## I. INTRODUCTION

The planetary orbits in the inner Solar System (ISS) are chaotic, with a Lyapunov time distributed around 5 million years (Myr) [1–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–7]. The dynamical half-life of Mercury orbit has recently been estimated at 30–40 billion years [4,7]. A disparity of nearly 4 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 dis-proportion 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 time-scale 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 quasiperiodic 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

<footnote>* Corresponding author.
federico.mogavero@obspm.fr</footnote>

<footnote>† Corresponding author.
nam.hoang-hoai@obspm.fr</footnote>

<footnote>‡ jacques.laskar@obspm.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
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 (see, e.g., Refs. [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 Sec. 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 or 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 such 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 precisely 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 (see, e.g., Ref. [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 quasisymmetries of the resonant harmonics of the Hamiltonian and the corresponding quasi-integrals (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 [4,21]

$$ \hat{H} = - \sum_{i=1}^{8} \left[ \sum_{l=1}^{i-1} \left\langle \frac{Gm_i m_l}{\|\mathbf{r}_i - \mathbf{r}_l\|} \right\rangle + \frac{3G^2 m_0^2 m_i}{c^2 a_i^2 \sqrt{1 - e_i^2}} \right]. \quad (1) $$

The planets are indexed in order of increasing semimajor axes ($a_i$)$_{i=1}^8$, $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 $\mathbf{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 nonresonant Fourier harmonics of the N-body Hamiltonian [4,21]. Hamiltonian (1) generates Gauss's dynamics of Keplerian rings [4,22], whose semimajor axes $a_i$ are constants of motion of the secular dynamics.

By developing the two-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,

$$ \begin{aligned} x_i &= \sqrt{\Lambda_i} \sqrt{1 - \sqrt{1 - e_i^2}} E^{j\Omega_i}, \\ y_i &= \sqrt{2\Lambda_i} (1 - e_i^2)^{1/4} \sin(\mathcal{I}_i/2) E^{j\Omega_i}, \end{aligned} \quad (2) $$

where $\Lambda_i = \mu_i [G(m_0 + m_i)a_i]^{1/2}$, $\mu_i = m_0 m_i / (m_0 + m_i)$ being the reduced masses of the planets, $\mathcal{I}_i$ the inclinations, $\Omega_i$ the longitudes of the perihelia, and $\Omega_i$ the longitudes of the nodes [27]. Pairs $(x_i, -j\bar{x}_i)$ and $(y_i, -j\bar{y}_i)$ are canonically conjugate momentum-coordinate variables. When truncating at a given total degree $2n$ in eccentricities and inclinations, the expansion provides Hamiltonians $\hat{H}_{2n} = \hat{H}_{2n}[(x_i, \bar{x}_i, y_i, \bar{y}_i)_{i=1}^8]$ 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 quasiperiodic form:

$$ x_i(t) = \sum_{l=1}^{M_i} \tilde{x}_{il} E^{jm_{il} \cdot \omega_o t}, \quad y_i(t) = \sum_{l=1}^{N_i} \tilde{y}_{il} E^{jn_{il} \cdot \omega_o t}, \quad (3) $$

for $i \in \{5, 6, 7, 8\}$, where $t$ denotes time, $\tilde{x}_{il}$ and $\tilde{y}_{il}$ are complex amplitudes, $\mathbf{m}_{il}$ and $\mathbf{n}_{il}$ integer vectors, and $\omega_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 (see Ref. [4], Appendix D). Gauss's dynamics of the forced ISS is obtained by substituting the predetermined time dependence in Eq. (1),
<table><caption>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 $\mathcal{H}_{2n}$ and $\mathbb{H}_{2n}$ are practically equivalent and treated as such. The $\mathbb{H}_{2n}^*$ models are introduced and discussed in Sec. IV D.</caption><thead><tr><th>Hamiltonian</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td>$\mathcal{H}$</td><td>Gauss's dynamics in complex Poincaré variables</td><td>Eq. (4)</td></tr><tr><td>$\mathcal{H}_{2n}$</td><td>Truncation of Gauss's dynamics at total degree 2n in eccentricities and inclinations</td><td>Ref. [4]</td></tr><tr><td>$\mathbb{H}_{2n}$</td><td>Truncated dynamics expressed in the action-angle variables of the Laplace-Lagrange dynamics $\mathcal{H}_2$</td><td>Eq. (6) and Ref. [4]</td></tr><tr><td>$\mathbb{H}_{2n}^*$</td><td>Fourier harmonics that involve outer planet modes other than $g_5$ are dropped from $\mathbb{H}_{2n}$</td><td>Eq. (24)</td></tr></tbody></table>

$$ \mathcal{H} = \hat{H}[(x_i, y_i)_{i=1}^4, (x_i = x_i(t), y_i = y_i(t))_{i=5}^8], \quad (4) $$

so that $\mathcal{H} = \mathcal{H}[(x_i, y_i)_{i=1}^4, 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 $\mathcal{H}_{2n}$ for the forced ISS is readily obtained by substituting Eq. (3) in the truncated Hamiltonian $\hat{\mathcal{H}}_{2n}$ of the entire Solar System. At the lowest degree, $\mathcal{H}_2$ generates a linear, forced Laplace-Lagrange (LL) dynamics. This can be analytically integrated by introducing complex proper mode variables $(u_i, v_i)_{i=1}^4$ via a time-dependent canonical transformation $(x_i, -j\bar{x}_i) \rightarrow (u_i, -j\bar{u}_i)$, $(y_i, -j\bar{y}_i) \rightarrow (v_i, -j\bar{v}_i)$ [4]. Action-angle pairs $(\mathbf{x}_i, \chi_i)$, $(\Psi_i, \psi_i)$ are introduced as

$$ u_i = \sqrt{\mathbf{x}_i} E^{-j\chi_i}, \quad v_i = \sqrt{\Psi_i} E^{-j\psi_i}. \quad (5) $$

When expressed in the proper modes, the truncated Hamiltonian can be expanded as a finite Fourier series:

$$ \mathbb{H}_{2n}(\mathbf{I}, \boldsymbol{\theta}, t) = \sum_{k, \ell} \tilde{\mathbb{H}}_{2n}^{k, \ell}(\mathbf{I}) E^{j[k\boldsymbol{\theta} + \ell\boldsymbol{\phi}(t)]}, \quad (6) $$

where $\mathbf{I} = (\mathbf{x}, \mathbf{\Psi})$ and $\boldsymbol{\theta} = (\boldsymbol{\chi}, \boldsymbol{\psi})$ are the eight-dimensional vectors of the action and angle variables, respectively, and we introduce the external angles $\boldsymbol{\phi}(t) = -\boldsymbol{\omega}_0 t$. The wave vectors $(\mathbf{k}, \boldsymbol{\ell})$ belong to a finite subset of $\mathbb{Z}^8 \times \mathbb{Z}^7$. At degree 2, one has $\mathbb{H}_2 = -\boldsymbol{\omega}_{LL} \cdot \mathbf{I}$, where $\boldsymbol{\omega}_{LL} = (\mathbf{g}_{LL}, \mathbf{s}_{LL}) \in \mathbb{R}^4 \times \mathbb{R}^4$ are the LL fundamental precession frequencies of the inner planet perihelia and nodes. Hamiltonian $\mathbb{H}_{2n}$ is in quasi-integrable form.

The quasiperiodic form of the outer orbits in Eq. (3) contains harmonics of order higher than one, that is, $\|\mathbf{m}_{il}\|_1 > 1$ and $\|\mathbf{n}_{il}\|_1 > 1$ for some $i$ and $l$, where $\|\cdot\|_1$ denotes the 1-norm. Therefore, the dynamics of $\mathcal{H}_{2n}$ and $\mathbb{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 $\dot{z} = F(z, t)$ defined on a phase space $\mathcal{P} \subseteq \mathbb{R}^P$ [31–34]. If $\Phi(z, t)$ denotes the associated flow and $z(t) = \Phi(z_0, t)$ the orbit that emanates from the initial condition $z_0$, the LCEs $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p$ are the logarithms of the eigenvalues of the matrix $\Lambda(z_0)$ defined as

$$ \lim_{t \to \infty} [\mathbf{M}(z_0, t)^{\mathrm{T}} \mathbf{M}(z_0, t)]^{1/2t} = \mathbf{\Lambda}(z_0), \quad (7) $$

where $\mathbf{M}(z_0, t) = \partial\Phi/\partial z_0$ is the fundamental matrix and $\mathbf{T}$ stands for transposition [32,33]. Introducing the Jacobian $\mathbf{J} = \partial F/\partial z$, the fundamental matrix allows us to write the solution of the variational equations $\dot{\boldsymbol{\zeta}} = \mathbf{J}(z(t), t)\boldsymbol{\zeta}$ as $\boldsymbol{\zeta}(t) = \mathbf{M}(z_0, t)\boldsymbol{\zeta}_0$, where $\boldsymbol{\zeta}(t) \in \mathcal{T}_{z(t)}\mathcal{P}$ belongs to the tangent space of $\mathcal{P}$ at point $z(t)$ and $\boldsymbol{\zeta}_0 = \boldsymbol{\zeta}(0)$. The multiplicative ergodic theorem of Oseledec [31] states that if $\rho$ 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 $\rho$-almost all $z_0$, and the LCEs are $\rho$-almost everywhere constant and only depend on $\rho$ [32]. Moreover, one has

$$ \lim_{t \to \infty} \frac{1}{t} \log \|\mathbf{M}(z_0, t)\boldsymbol{\zeta}_0\| = \lambda^{(i)} \quad \text{if } \boldsymbol{\zeta}_0 \in E_{z_0}^{(i)} \setminus E_{z_0}^{(i+1)}, \quad (8) $$

for $\rho$-almost all $z_0$, where $\lambda^{(1)} > \lambda^{(2)} > \dots$ are the LCEs without repetition by multiplicity, and $E_{z_0}^{(i)}$ is the subspace of $\mathbb{R}^P$ corresponding to the eigenvalues of $\mathbf{\Lambda}(z_0)$ that are smaller than or equal to $\exp\lambda^{(i)}$, with $\mathcal{T}_{z_0}\mathcal{P} = E_{z_0}^{(1)} \supset E_{z_0}^{(2)} \supset \dots$. The specific choice of the $\mathbb{R}^P$-norm $\|\cdot\|$ in Eq. (8) is irrelevant [32,34]. Once the LCEs have been introduced, a characteristic timescale can be defined from
each positive exponent as $\lambda_i^{-1}$. In the case of the maximum Lyapunov exponent $\lambda_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,

$$ \Delta\lambda_i := \lambda_i + \lambda_{2p-i+1} = 0 \quad \text{for all } 1 \le i \le p. \quad (9) $$

If the Hamiltonian is time independent, a pair of exponents vanishes. In general, the existence of an integral of motion $\mathcal{G} = \mathcal{G}(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 $\mathcal{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), $\lambda_i(z_0, t)$, defined from the eigenvalues $\mathbf{m}_1 \ge \mathbf{m}_2 \ge \dots \ge \mathbf{m}_p$ of the time-dependent symmetric positive-defined matrix $\mathbf{M}(z_0, t)^T \mathbf{M}(z_0, t)$ as

$$ \lambda_i(z_0, t) = \frac{1}{2t} \log \mathbf{m}_i(z_0, t). \quad (10) $$

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 nontrivial 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 $\mathcal{H}$ with initial conditions very close to their nominal values [38]. Its long-term distribution is quite large and does not shrink over time (see Ref. [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 (see Ref. [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 nonergodic 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 30 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 [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 nonergodic 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 ($\mathcal{H}$) and a few minutes at degree 4 ($\mathcal{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 $\mathcal{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 time step of 250 yr. 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 1(a) 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 (Ref. [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. 8(a)]. 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. 8(b)].

The spectrum in Fig. 1(a) has distinctive features. A set of intermediate exponents follow the FT-MLE, ranging from 0.1 to 0.01" yr$^{-1}$, while the smallest ones fall below 0.01" yr$^{-1}$. Figure 1(a) reveals the existence of a hierarchy
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 1. Positive FT-LCEs $\lambda_i$ of the forced secular ISS from Hamiltonians $\mathcal{H}_4$ (a) and $\mathcal{H}_4^*$ [Eq. (24)] (b), and corresponding characteristic timescales $\lambda_i^{-1}$. 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 $\mathcal{H}_4^*$ model is introduced and discussed in Sec. IV D.</figcaption></figure>

of exponents and corresponding timescales that spans 2 orders of magnitude, down to a median value of $\lambda_8^{-1} \approx 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 $\mathcal{H}_6$ than for $\mathcal{H}_4$ [7].

The relevance of the Lyapunov spectrum in Fig. 1(a) 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. 1(a) suggests that the secular orbits of the inner planets are characterized by a slow-fast 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 phase-space 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 pendulumlike dynamics generated by each Fourier harmonic that composes the Hamiltonian in Eq. (6) [44]. One can write $\mathcal{H}_{2n}(\mathbf{I}, \boldsymbol{\theta}, t) = \tilde{\mathcal{H}}_{2n}^{0,0}(\mathbf{I}) + \sum_{i=1}^{M_{2n}} \mathcal{F}_i(\mathbf{I}, \boldsymbol{\theta}, t)$, with

$$ \mathcal{F}_i(\mathbf{I}, \boldsymbol{\theta}, t) = \tilde{\mathcal{H}}_{2n}^{k_i, \ell_i}(\mathbf{I}) \mathrm{E}^{j[k_i \boldsymbol{\theta} + \boldsymbol{\ell}_i \boldsymbol{\phi}(t)]} + \text{c.c.}, \quad (11) $$

where $(k^i, \ell^i) \neq (0, 0)$, $M_{2n}$ is the number of harmonics in $\mathcal{H}_{2n}$ with a non-null wave vector, and c.c. stands for complex conjugate. Chaos arises from the interaction of resonant harmonics, that is, those harmonics $\mathcal{F}_i$ whose frequency combination $k^i \cdot \dot{\boldsymbol{\theta}} + \ell^i \cdot \dot{\boldsymbol{\phi}}(t)$ vanishes at some point along the motion. Using the computer algebra system TRIP, the harmonics of $\mathcal{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 $\Delta\omega$ [19]. The resonances have then been ordered by decreasing time median of their half-widths. The resulting ranking of resonances is denoted as $\mathcal{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 $(g_i, s_i)_{i=1}^8$, that is, $k \cdot \boldsymbol{\omega}_i + \boldsymbol{\ell} \cdot \boldsymbol{\omega}_o$, with $\boldsymbol{\omega}_i = (g_1, g_2, g_3, g_4, s_1, s_2, s_3, s_4)$. Table II also shows the order of each harmonic, defined
<table><caption>TABLE II. Top of ranking $\mathcal{R}_1$. First 30 resonant harmonics of $\mathbb{H}_{10}$ along the 5 Gyr nominal solution of Gauss's dynamics, in order of decreasing time median of the resonance half-width $\Delta\omega$ (arcsec yr$^{-1}$). $\mathcal{O}$ is the order of the harmonic. $\tau^{\text{res}}$ is the fraction of time the harmonic is resonant. Only harmonics with $\tau^{\text{res}} > 1\%$ are shown. 5th and 95th percentiles of the time distribution of $\Delta\omega$ as subscripts and superscripts, respectively. Adapted from Table 2 of Ref. [19].</caption><thead><tr><th><i>i</i></th><th>Fourier harmonic <i>F</i><sub><i>i</i></sub></th><th><i>O</i><sub><i>i</i></sub></th><th>τ<sup>res</sup><sub><i>i</i></sub></th><th>Δω<sub><i>i</i></sub></th></tr></thead><tbody><tr><td>1</td><td>g<sub>3</sub> - g<sub>4</sub> - s<sub>3</sub> + s<sub>4</sub></td><td>4</td><td>12%</td><td>0.332<sup>0.526</sup><sub>0.093</sub></td></tr><tr><td>2</td><td>g<sub>1</sub> - g<sub>2</sub> + s<sub>1</sub> - s<sub>2</sub></td><td>4</td><td>19%</td><td>0.302<sup>0.641</sup><sub>0.034</sub></td></tr><tr><td>3</td><td>g<sub>2</sub> - g<sub>5</sub> - 2s<sub>1</sub> + 2s<sub>2</sub></td><td>6</td><td>23%</td><td>0.105<sup>0.223</sup><sub>0.041</sub></td></tr><tr><td>4</td><td>2g<sub>3</sub> - 2g<sub>4</sub> - s<sub>3</sub> + s<sub>4</sub></td><td>6</td><td>70%</td><td>0.076<sup>0.159</sup><sub>0.023</sub></td></tr><tr><td>5</td><td>g<sub>1</sub> - g<sub>5</sub> - s<sub>1</sub> + s<sub>2</sub></td><td>4</td><td>10%</td><td>0.074<sup>0.178</sup><sub>0.056</sub></td></tr><tr><td>6</td><td>g<sub>2</sub> - g<sub>4</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>4</td><td>6%</td><td>0.066<sup>0.098</sup><sub>0.025</sub></td></tr><tr><td>7</td><td>g<sub>1</sub> - 2g<sub>2</sub> + g<sub>4</sub> + s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>8</td><td>5%</td><td>0.061<sup>0.074</sup><sub>0.051</sub></td></tr><tr><td>8</td><td>g<sub>1</sub> - g<sub>3</sub> + s<sub>2</sub> - s<sub>3</sub></td><td>4</td><td>17%</td><td>0.056<sup>0.090</sup><sub>0.028</sub></td></tr><tr><td>9</td><td>g<sub>1</sub> + g<sub>3</sub> - 2g<sub>4</sub> + s<sub>2</sub> - s<sub>3</sub></td><td>6</td><td>5%</td><td>0.053<sup>0.061</sup><sub>0.037</sub></td></tr><tr><td>10</td><td>3g<sub>3</sub> - 3g<sub>4</sub> - s<sub>3</sub> + s<sub>4</sub></td><td>8</td><td>9%</td><td>0.052<sup>0.140</sup><sub>0.007</sub></td></tr><tr><td>11</td><td>g<sub>2</sub> - g<sub>3</sub> - s<sub>1</sub> + 2s<sub>2</sub> - s<sub>3</sub></td><td>6</td><td>5%</td><td>0.038<sup>0.047</sup><sub>0.028</sub></td></tr><tr><td>12</td><td>g<sub>1</sub> - 2g<sub>3</sub> + g<sub>4</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>6</td><td>36%</td><td>0.038<sup>0.083</sup><sub>0.016</sub></td></tr><tr><td>13</td><td>2g<sub>1</sub> - g<sub>3</sub> - g<sub>5</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>6</td><td>5%</td><td>0.037<sup>0.043</sup><sub>0.028</sub></td></tr><tr><td>14</td><td>g<sub>4</sub> - g<sub>5</sub> - s<sub>2</sub> + 2s<sub>3</sub> - s<sub>4</sub></td><td>6</td><td>2%</td><td>0.033<sup>0.036</sup><sub>0.031</sub></td></tr><tr><td>15</td><td>g<sub>1</sub> - 2g<sub>3</sub> + g<sub>4</sub> + s<sub>1</sub> + s<sub>3</sub> - 2s<sub>4</sub></td><td>8</td><td>25%</td><td>0.033<sup>0.045</sup><sub>0.014</sub></td></tr><tr><td>16</td><td>g<sub>1</sub> - g<sub>4</sub> + s<sub>1</sub> - s<sub>4</sub></td><td>4</td><td>23%</td><td>0.032<sup>0.054</sup><sub>0.017</sub></td></tr><tr><td>17</td><td>g<sub>1</sub> - 2g<sub>2</sub> + g<sub>5</sub> + 3s<sub>1</sub> - 3s<sub>2</sub></td><td>10</td><td>6%</td><td>0.032<sup>0.039</sup><sub>0.023</sub></td></tr><tr><td>18</td><td>g<sub>1</sub> - g<sub>4</sub> + s<sub>2</sub> - s<sub>3</sub></td><td>4</td><td>18%</td><td>0.031<sup>0.073</sup><sub>0.016</sub></td></tr><tr><td>19</td><td>3g<sub>1</sub> - g<sub>2</sub> - g<sub>4</sub> - g<sub>5</sub> + s<sub>1</sub> - s<sub>3</sub></td><td>8</td><td>2%</td><td>0.031<sup>0.042</sup><sub>0.023</sub></td></tr><tr><td>20</td><td>2g<sub>1</sub> - g<sub>2</sub> - g<sub>3</sub> + s<sub>1</sub> - s<sub>3</sub></td><td>6</td><td>29%</td><td>0.028<sup>0.051</sup><sub>0.016</sub></td></tr><tr><td>21</td><td>2g<sub>1</sub> - g<sub>2</sub> - g<sub>4</sub> + s<sub>1</sub> - s<sub>3</sub></td><td>6</td><td>3%</td><td>0.026<sup>0.028</sup><sub>0.021</sub></td></tr><tr><td>22</td><td>3g<sub>3</sub> - 3g<sub>4</sub> - 2s<sub>3</sub> + 2s<sub>4</sub></td><td>10</td><td>8%</td><td>0.025<sup>0.055</sup><sub>0.012</sub></td></tr><tr><td>23</td><td>2g<sub>1</sub> - g<sub>2</sub> - 2g<sub>3</sub> + g<sub>4</sub> + s<sub>1</sub> - s<sub>4</sub></td><td>8</td><td>3%</td><td>0.023<sup>0.036</sup><sub>0.012</sub></td></tr><tr><td>24</td><td>2g<sub>3</sub> - g<sub>4</sub> - g<sub>5</sub> - s<sub>1</sub> + s<sub>4</sub></td><td>6</td><td>16%</td><td>0.023<sup>0.048</sup><sub>0.010</sub></td></tr><tr><td>25</td><td>g<sub>1</sub> - 3g<sub>3</sub> + 2g<sub>4</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>8</td><td>7%</td><td>0.021<sup>0.030</sup><sub>0.008</sub></td></tr><tr><td>26</td><td>g<sub>1</sub> - g<sub>2</sub> - g<sub>3</sub> + g<sub>4</sub> + s<sub>1</sub> - s<sub>2</sub></td><td>6</td><td>6%</td><td>0.021<sup>0.032</sup><sub>0.004</sub></td></tr><tr><td>27</td><td>g<sub>1</sub> + g<sub>3</sub> - 2g<sub>4</sub> + s<sub>1</sub> - s<sub>4</sub></td><td>6</td><td>3%</td><td>0.021<sup>0.022</sup><sub>0.017</sub></td></tr><tr><td>28</td><td>g<sub>1</sub> + g<sub>2</sub> - 2g<sub>5</sub> - 3s<sub>1</sub> + 3s<sub>2</sub></td><td>10</td><td>4%</td><td>0.020<sup>0.028</sup><sub>0.006</sub></td></tr><tr><td>29</td><td>3g<sub>1</sub> - g<sub>2</sub> - g<sub>4</sub> - g<sub>5</sub> + s<sub>2</sub> - s<sub>3</sub></td><td>8</td><td>4%</td><td>0.020<sup>0.027</sup><sub>0.008</sub></td></tr><tr><td>30</td><td>2g<sub>1</sub> - g<sub>4</sub> - g<sub>5</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>6</td><td>7%</td><td>0.020<sup>0.029</sup><sub>0.007</sub></td></tr></tbody></table>

as the even integer $\mathcal{O} = \|(\mathbf{k}, \mathcal{L})\|_1$. The support of the asymptotic ensemble distribution of the FT-MLE shown in Fig. 1(a) overlaps in a robust way with that of the time distribution of the half-width of the strongest resonances. In other words,

$$2\pi\lambda_1 \approx \Delta\omega^{\mathcal{R}_1}, \qquad (12)$$

where $\Delta\omega^{\mathcal{R}_1}$ stands for the half-width of the uppermost resonances of ranking $\mathcal{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 FT-MLE beyond 1 Gyr [4,19], and that extends to the entire spectrum in Fig. 1(a).

We remark that, strictly speaking, ranking $\mathcal{R}_1$ is established on the Fourier harmonics of the Lie-transformed Hamiltonian $\mathbb{H}'_{2n}$ (see Ref. [19], Appendix G). New canonical variables are indeed defined to transform $\mathbb{H}_{2n}$ in a Birkhoff normal form to degree 4. The goal is to let the interactions of the terms of degree 4 in $\mathbb{H}_{2n}$ appear more explicitly in the amplitudes of the harmonics of higher degrees in $\mathbb{H}'_{2n}$, the physical motivation being that the nonlinear interaction of the harmonics at degree 4 constitutes the primary source of chaos [19]. Keeping in mind the quasi-identity 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. Quasisymmetries 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,

$$\mathbb{H}_{2n,m} = \tilde{\mathbb{H}}_{2n}^{0,0} + \sum_{i=1}^{m} \mathcal{F}_{i}. \qquad (13)$$

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 $\mathbb{H}_{2n}$ and a dynamical quantity that is a linear combination of the action variables, that is,

$$C_{\gamma} = \gamma \cdot I, \qquad (14)$$

$\gamma \in \mathbb{R}^8$ being a parameter vector. From Eq. (11), the partial contribution of the $m$ harmonics to the time derivative of $C_{\gamma}$ along the flow of $\mathbb{H}_{2n}$ is

$$\dot{C}_{\gamma,m} = 2 \sum_{i=1}^{m} \gamma \cdot k^i \mathrm{Im}\{\tilde{\mathbb{H}}_{2n}^{k^i, \mathcal{L}^i}(\mathbf{I}) \mathrm{E}^{j[k^i \cdot \theta + \mathcal{L}^i \cdot \phi(t)]}\}, \qquad (15)$$
and $\dot{C}_{\gamma} = \dot{C}_{\gamma, M_{2n}}$, where $M_{2n}$ is the total number of harmonics with a non-null wave vector that appear in $\mathbb{H}_{2n}$. Any quantity $C_{\gamma}$ with $\gamma \cdot k^i = 0$ is conserved by the one-d.o.f. dynamics generated by the single harmonic $\mathcal{F}_i$. In other words, such a quantity would be an integral of motion if $\mathcal{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 $C_{\gamma}$ if $\gamma \perp \text{span}(k^1, k^2, \dots, k^m)$, that is, if the vector $\gamma$ belongs to the orthogonal complement of the linear subspace of $\mathbb{R}^8$ spanned by the wave vectors $(k^i)_{i=1}^m$. We also consider the quantity

$$C'_{\gamma} = \mathbb{H}_{2n} + \gamma \cdot I. \qquad (16)$$

Because of the explicit time dependence in the Hamiltonian, the partial contribution of a set of $m$ harmonics to the time derivative of $C'_{\gamma}$ along the flow of $\mathbb{H}_{2n}$ is

$$\dot{C}'_{\gamma,m} = 2 \sum_{i=1}^{m} (\gamma \cdot k^i + \mathcal{L}^i \cdot \omega_o) \text{Im}\{\tilde{\mathbb{H}}_{2n}^{k^i, \mathcal{L}^i}(\mathbf{I}) \text{E}^{j[k^i \cdot \theta + \mathcal{L}^i \cdot \phi(t)]}\}, \qquad (17)$$

and one has $\dot{C}'_{\gamma} = \dot{C}'_{\gamma, M_{2n}}$. Quantity $C'_{\gamma}$ is unchanged by the $m$ harmonics if $\gamma \cdot k^i + \mathcal{L}^i \cdot \omega_o = 0$ for $i \in \{1, 2, \dots, m\}$. Dynamical quantities $C_{\gamma}$ or $C'_{\gamma}$ that are unaffected by a given set of leading harmonics, that is, with null partial contribution in Eq. (15) or (17), are denoted as quasi-integrals of motion from now on. More specifically, we build our analysis on ranking $\mathcal{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 $\sum_{i=1}^{8} k_i + \sum_{i=1}^{7} \ell_i = 0$, where $k = (k_1, k_2, \dots, k_8)$ and $\ell = (\ell_1, \ell_2, \dots, \ell_7)$ [21,24,45,46]. Moreover, the Jupiter-dominated eccentricity mode $g_5$ is the only fundamental Fourier mode of the outer planet forcing to appear in Table II. The quantity

$$\mathcal{E}_{2n} := C'_{g_5 1_8} = \mathbb{H}_{2n} + g_5 \sum_{i=1}^{4} (\mathcal{X}_i + \Psi_i), \qquad (18)$$

with $1_8 = (1, 1, \dots, 1) \in \mathbb{R}^8$, is therefore unaffected by the resonances listed in Table II. In an equivalent way, the time-dependent canonical transformation $\theta \to \theta + g_5 t 1_8$, with unchanged action variables, allows us to remove the explicit time dependence in these harmonics. Quantity $\mathcal{E}_{2n}$ coincides with the

transformed Hamiltonian, and the harmonics in Table II do not contribute to its time derivative.

(b) Second symmetry.—We write the eccentricity and inclination parts of the harmonic wave vectors explicitly, that is, $k = (k^{\text{ecc}}, k^{\text{inc}})$ with $k^{\text{ecc}}, k^{\text{inc}} \in \mathbb{R}^4$. One can visually check that the harmonics in Table II verify the relation $\sum_{i=1}^{4} k_i^{\text{inc}} = 0$, where $k^{\text{inc}} = (k_1^{\text{inc}}, \dots, k_4^{\text{inc}})$. Therefore, denoting $\gamma_1 = (\mathbf{0}_4, \mathbf{1}_4)$, the quantity

$$C_{\text{inc}} := C_{\gamma_1} = \Psi_1 + \Psi_2 + \Psi_3 + \Psi_4 \qquad (19)$$

is conserved by these resonances. $C_{\text{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_{\text{ecc}} = \sum_{i=1}^{4} \mathcal{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_{\text{inc}}$ depends on two facts: the inclination modes $s_6, s_7, s_8$ of the external forcing do not appear in Table II and 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 quasiperiodic contributions), so that two AMD reservoirs $C_{\text{ecc}}$ and $C_{\text{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 \times 8$ matrix $\mathbf{K}_{30}$ whose rows are the wave vectors $(k^i)_{i=1}^{30}$ of the listed resonances. A singular value decomposition shows that the rank of $\mathbf{K}_{30}$ is equal to 6. Therefore, the linear subspace $V_{30} = \text{span}(k^1, k^2, \dots, 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 $V_{30}^{\perp}$. One choice consists in $V_{30}^{\perp} = \text{span}(\gamma_2, \gamma_2^{\perp})$, with

$$\begin{align}
C_2 &:= C_{\gamma_2} = -\mathcal{X}_3 - \mathcal{X}_4 + \Psi_1 + \Psi_2 + 2\Psi_3 + 2\Psi_4, \nonumber \\
C_2^{\perp} &:= C_{\gamma_2^{\perp}} = \mathcal{X}_3 + \mathcal{X}_4 + \Psi_1 + \Psi_2. \tag{20}
\end{align}$$

Since the second symmetry clearly requires that $\gamma_1 \in V_{30}^{\perp}$, the three quantities $C_{\text{inc}}, C_2, C_2^{\perp}$ are not independent and one has indeed $C_{\text{inc}} = (C_2 + C_2^{\perp})/2$. We remark that $(C_2 - C_2^{\perp})/2 = -\mathcal{X}_3 - \mathcal{X}_4 + \Psi_3 + \Psi_4$. 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 Table II, naturally represent quasisymmetries when considering the entire spectrum of resonances $R_1$. They are indeed broken at some point by weak resonances (see Sec. IV C). Quantities $\mathcal{E}_{2n}$, $C_{\text{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_{\text{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.

We remark that, differently from $C_{\text{inc}}$ and $C_2$, the quantity $\mathcal{E}_{2n}$ is a nonlinear 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 $\mathcal{E}_2 = \mathbb{H}_2 + g_5 \mathbf{1}_8 \cdot \mathbf{I} = C_{\gamma_3}$, with $\gamma_3 = -\omega_{\text{LL}} + g_5 \mathbf{1}_8$, reproduces reasonably well $\mathcal{E}_{2n}$ along the flow of $\mathbb{H}_{2n}$ for $n > 1$. The vector $\gamma_3$ is used in Sec. V, together with $\gamma_1$ and $\gamma_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 $\gamma_1$ and $\gamma_2$, the components of $\gamma_3$ are not integer and they have the dimension of a frequency.

## B. Slow variables

The QIs of motion $\mathcal{E}_{2n}$, $C_{\text{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,

$$ \hat{C}_{\text{inc}} = \frac{C_{\text{inc}}}{\|\gamma_1\|C_0}, \quad \hat{C}_2 = \frac{C_2}{\|\gamma_2\|C_0}, \quad \hat{\mathcal{E}}_{2n} = \frac{\mathcal{E}_{2n}}{\|\gamma_3\|C_0}, \quad (21) $$

where $C_0$ stands for the current total AMD of the inner planets, that is, the value of $C_{\text{ecc}} + C_{\text{inc}}$ at time zero. We stress that, by introducing the unit vectors $\hat{\gamma}_i = \gamma_i / \|\gamma_i\|$ for $i \in \{1, 2, 3\}$, one has $\hat{C}_{\text{inc}} = C_{\hat{\gamma}_1}/C_0$ and $\hat{C}_2 = C_{\hat{\gamma}_2}/C_0$. At degree 2, one also has $\hat{\mathcal{E}}_2 = C_{\hat{\gamma}_3}/C_0$. We then consider the ensembles of numerical integrations of $\mathcal{H}_4$ and $\mathcal{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 $\hat{\mathbf{I}} = \mathbf{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 filter with three iterations of the moving average [4,48]. A cutoff frequency of $1 \text{ Myr}^{-1}$ is chosen to highlight the long-term diffusion that can be hidden by short-time quasiperiodic 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 $\mathcal{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 $\hat{\gamma}_2$ and $\hat{\gamma}_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 $\mathcal{X}_1$ and $\mathcal{X}_2$, while $\mathcal{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 $\mathcal{X}_1$.

## C. Weak resonances and Lyapunov spectrum

A fundamental result from Table II is that the symmetries introduced in Sec. IVA 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 resonances $R_2$ is defined in this way. Table III reports the ten strongest symmetry-breaking resonances that change $\mathcal{E}_{2n}$, $C_{\text{inc}}$, $C_2$, respectively. As in Table 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_{\text{inc}}$, there is also a contribution that starts at about 0.006" yr$^{-1}$ with $\mathcal{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_{\text{inc}}$ than for $\mathcal{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_{\text{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
<figure><img src="image_5.png" /><figcaption>FIG. 2. Time evolution over 5 Gyr of the dimensionless QIs ($\hat{C}_{\text{inc}}$, $\hat{C}_2$, $\hat{\varepsilon}$) and of two representatives of the dimensionless action variables ($\hat{x}_1$, $\hat{\Psi}_3$) along the nominal orbital solutions of different models. Top row: $H_4$ and $H_6$ ($\hat{\varepsilon}$ stands for $\hat{\varepsilon}_4$ and $\hat{\varepsilon}_6$, respectively). Bottom row: $H_4^*$ and $H_6^*$ from Eq. (24) ($\hat{\varepsilon}_{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.</figcaption></figure>

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 Sec. 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. 1(a) 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, $\lambda_6$, $\lambda_7$, $\lambda_8$. One can indeed write

$$2\pi\lambda_6 \approx \Delta\omega^{R_2}, \qquad (22)$$

where $\Delta\omega^{R_2}$ stands for the half-width of the uppermost resonances of ranking $R_2$. Table III and Fig. 1(a) suggest a relation between the QIs and the smallest Lyapunov exponents:

$$\lambda_6, \lambda_7, \lambda_8 \leftrightarrow \mathcal{E}_{2n}, C_{\text{inc}}, C_2. \qquad (23)$$

Equation (23) is not a one-to-one correspondence, nor should it be understood as an exact relation since, for example, $\lambda_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 Sec. V, we show its statistical validity in the geometric framework established by a principal component analysis of the orbital solutions. Moreover, Sec. 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 $\mathcal{E}_{2n}$, $C_{\text{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 $\gamma_1, \gamma_2, \gamma_3$ are linearly independent. On the other hand, one has the Poisson bracket $\{C_{\text{inc}}, C_2\} = 0$, since the two quantities are functions of the action variables only. One also has $\{\mathcal{E}_{2n}, C_{\text{inc}}\} = \{H_{2n}, C_{\text{inc}}\} = \dot{C}_{\text{inc}}$ and $\{\mathcal{E}_{2n}, C_2\} = \dot{C}_2$. Only weak resonances contribute to these Poisson brackets and the three QIs are therefore nearly in involution.

## D. New truncation of the Hamiltonian

The fundamental role of the external modes $g_6, g_7, s_6, s_7$ in Table III 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$. Table IV reports the ten 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 $\mathcal{E}_{2n}$. For the remaining two QIs, the only harmonics that appear in Table IV are of
<table><caption>TABLE III. Top of ranking $\mathcal{R}_2$. First ten symmetry-breaking resonances of $\mathbb{H}_{10}$ along the 5 Gyr nominal solution of Gauss's dynamics, that change $\mathcal{E}_{2n}$, $C_{\text{inc}}$, and $C_2$, respectively (see Table II for details).</caption><thead><tr><th>i</th><th>Fourier harmonic F<sub>i</sub></th><th>O<sub>i</sub></th><th>τ<sup>res</sup><sub>i</sub></th><th>Δω<sub>i</sub></th></tr></thead><tbody><tr><td colspan="5"><i>E</i><sub>2n</sub></td></tr><tr><td>1</td><td>g<sub>1</sub> + g<sub>2</sub> - 2g<sub>5</sub> + s<sub>2</sub> - s<sub>7</sub></td><td>6</td><td>1%</td><td>0.018<sup>0.020</sup><sub>0.016</sub></td></tr><tr><td>2</td><td>g<sub>1</sub> - 2g<sub>2</sub> + g<sub>6</sub> - s<sub>2</sub> + s<sub>6</sub></td><td>6</td><td>4%</td><td>0.017<sup>0.024</sup><sub>0.008</sub></td></tr><tr><td>3</td><td>g<sub>3</sub> - g<sub>6</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>4</td><td>10%</td><td>0.017<sup>0.022</sup><sub>0.009</sub></td></tr><tr><td>4</td><td>g<sub>5</sub> - g<sub>7</sub> + s<sub>3</sub> - s<sub>4</sub></td><td>4</td><td>5%</td><td>0.017<sup>0.024</sup><sub>0.010</sub></td></tr><tr><td>5</td><td>g<sub>4</sub> - g<sub>6</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>4</td><td>4%</td><td>0.016<sup>0.021</sup><sub>0.007</sub></td></tr><tr><td>6</td><td>g<sub>2</sub> - 2g<sub>4</sub> + g<sub>6</sub></td><td>4</td><td>12%</td><td>0.016<sup>0.027</sup><sub>0.009</sub></td></tr><tr><td>7</td><td>g<sub>1</sub> - g<sub>3</sub> - g<sub>5</sub> + g<sub>6</sub> - s<sub>1</sub> + s<sub>4</sub></td><td>6</td><td>6%</td><td>0.015<sup>0.024</sup><sub>0.011</sub></td></tr><tr><td>8</td><td>2g<sub>3</sub> - 2g<sub>4</sub> + g<sub>5</sub> - g<sub>7</sub></td><td>6</td><td>25%</td><td>0.014<sup>0.022</sup><sub>0.006</sub></td></tr><tr><td>9</td><td>g<sub>2</sub> + g<sub>3</sub> - 3g<sub>4</sub> + g<sub>6</sub></td><td>6</td><td>3%</td><td>0.014<sup>0.017</sup><sub>0.010</sub></td></tr><tr><td>10</td><td>g<sub>2</sub> - 2g<sub>3</sub> + g<sub>6</sub> + s<sub>3</sub> - s<sub>4</sub></td><td>6</td><td>3%</td><td>0.011<sup>0.018</sup><sub>0.003</sub></td></tr><tr><td colspan="5"><i>C</i><sub>inc</sub></td></tr><tr><td>1</td><td>g<sub>1</sub> + g<sub>2</sub> - 2g<sub>5</sub> + s<sub>2</sub> - s<sub>7</sub></td><td>6</td><td>1%</td><td>0.018<sup>0.020</sup><sub>0.016</sub></td></tr><tr><td>2</td><td>g<sub>1</sub> - 2g<sub>2</sub> + g<sub>6</sub> - s<sub>2</sub> + s<sub>6</sub></td><td>6</td><td>4%</td><td>0.017<sup>0.024</sup><sub>0.008</sub></td></tr><tr><td>3</td><td>g<sub>2</sub> - g<sub>6</sub> + s<sub>1</sub> - s<sub>6</sub></td><td>4</td><td>8%</td><td>0.011<sup>0.017</sup><sub>0.003</sub></td></tr><tr><td>4</td><td>g<sub>4</sub> - g<sub>6</sub> - 2s<sub>3</sub> + 3s<sub>4</sub> - s<sub>6</sub></td><td>8</td><td>5%</td><td>0.010<sup>0.012</sup><sub>0.002</sub></td></tr><tr><td>5</td><td>2g<sub>1</sub> - 2g<sub>5</sub> + s<sub>1</sub> - s<sub>7</sub></td><td>6</td><td>1%</td><td>0.007<sup>0.008</sup><sub>0.005</sub></td></tr><tr><td>6</td><td>4g<sub>1</sub> - 3g<sub>2</sub> - g<sub>5</sub> - s<sub>2</sub> + s<sub>7</sub></td><td>10</td><td>7%</td><td>0.006<sup>0.015</sup><sub>0.003</sub></td></tr><tr><td>7</td><td>g<sub>3</sub> - g<sub>6</sub> - 2s<sub>3</sub> + 3s<sub>4</sub> - s<sub>6</sub></td><td>8</td><td>2%</td><td>0.006<sup>0.011</sup><sub>0.003</sub></td></tr><tr><td>8</td><td>4g<sub>1</sub> - g<sub>2</sub> - g<sub>3</sub> - s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>10</td><td>2%</td><td>0.006<sup>0.009</sup><sub>0.001</sub></td></tr><tr><td>9</td><td>2g<sub>1</sub> - g<sub>2</sub> - g<sub>5</sub> + 3s<sub>1</sub> - 2s<sub>2</sub> - s<sub>7</sub></td><td>10</td><td>3%</td><td>0.006<sup>0.008</sup><sub>0.002</sub></td></tr><tr><td>10</td><td>3g<sub>1</sub> - 3g<sub>2</sub> + s<sub>1</sub> - 2s<sub>2</sub> + s<sub>7</sub></td><td>10</td><td>19%</td><td>0.006<sup>0.009</sup><sub>0.002</sub></td></tr><tr><td colspan="5"><i>C</i><sub>2</sub></td></tr><tr><td>1</td><td>g<sub>1</sub> + g<sub>2</sub> - 2g<sub>5</sub> + s<sub>2</sub> - s<sub>7</sub></td><td>6</td><td>1%</td><td>0.018<sup>0.020</sup><sub>0.016</sub></td></tr><tr><td>2</td><td>g<sub>1</sub> - 2g<sub>2</sub> + g<sub>6</sub> - s<sub>2</sub> + s<sub>6</sub></td><td>6</td><td>4%</td><td>0.017<sup>0.024</sup><sub>0.008</sub></td></tr><tr><td>3</td><td>g<sub>3</sub> - g<sub>6</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>4</td><td>10%</td><td>0.017<sup>0.022</sup><sub>0.009</sub></td></tr><tr><td>4</td><td>g<sub>4</sub> - g<sub>6</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>4</td><td>4%</td><td>0.016<sup>0.021</sup><sub>0.007</sub></td></tr><tr><td>5</td><td>g<sub>2</sub> - 2g<sub>4</sub> + g<sub>6</sub></td><td>4</td><td>12%</td><td>0.016<sup>0.027</sup><sub>0.009</sub></td></tr><tr><td>6</td><td>g<sub>1</sub> - g<sub>3</sub> - g<sub>5</sub> + g<sub>6</sub> - s<sub>1</sub> + s<sub>4</sub></td><td>6</td><td>6%</td><td>0.015<sup>0.024</sup><sub>0.011</sub></td></tr><tr><td>7</td><td>g<sub>2</sub> + g<sub>3</sub> - 3g<sub>4</sub> + g<sub>6</sub></td><td>6</td><td>3%</td><td>0.014<sup>0.017</sup><sub>0.010</sub></td></tr><tr><td>8</td><td>g<sub>2</sub> - 2g<sub>3</sub> + g<sub>6</sub> + s<sub>3</sub> - s<sub>4</sub></td><td>6</td><td>3%</td><td>0.011<sup>0.018</sup><sub>0.003</sub></td></tr><tr><td>9</td><td>g<sub>1</sub> - g<sub>3</sub> - g<sub>4</sub> + g<sub>6</sub> + s<sub>1</sub> - s<sub>2</sub></td><td>6</td><td>4%</td><td>0.011<sup>0.013</sup><sub>0.008</sub></td></tr><tr><td>10</td><td>g<sub>2</sub> - g<sub>6</sub> + s<sub>1</sub> - s<sub>6</sub></td><td>4</td><td>8%</td><td>0.011<sup>0.017</sup><sub>0.003</sub></td></tr></tbody></table>

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_{\text{inc}}$, the half-width of the uppermost resonances is now around 0.005" yr$^{-1}$. One can appreciate that the activation times $\tau^{\text{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 uppermost resonances drops by 2 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

<table><caption>TABLE IV. Top of ranking $\mathcal{R}_3$. First ten symmetry-breaking resonances of $\mathbb{H}_{10}$ along the 5 Gyr nominal solution of Gauss's dynamics, that only involve $g_5$ among the external modes and change $\mathcal{E}_{2n}$, $C_{\text{inc}}$, and $C_2$, respectively.</caption><thead><tr><th>i</th><th>Fourier harmonic F<sub>i</sub></th><th>O<sub>i</sub></th><th>τ<sup>res</sup><sub>i</sub></th><th>Δω<sub>i</sub></th></tr></thead><tbody><tr><td colspan="5"><i>E</i><sub>2n</sub></td></tr><tr><td colspan="5">No resonances</td></tr><tr><td colspan="5"><i>C</i><sub>inc</sub><sup>a</sup></td></tr><tr><td>1</td><td>4g<sub>1</sub> - g<sub>2</sub> - g<sub>3</sub> - s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>10</td><td>2.0%</td><td>0.006<sup>0.009</sup><sub>0.001</sub></td></tr><tr><td>2</td><td>2g<sub>1</sub> + g<sub>2</sub> - g<sub>4</sub> - 2s<sub>1</sub> - s<sub>2</sub> + s<sub>3</sub></td><td>8</td><td>2.5%</td><td>0.003<sup>0.008</sup><sub>0.002</sub></td></tr><tr><td>3</td><td>3g<sub>1</sub> - g<sub>2</sub> - s<sub>1</sub> - 3s<sub>2</sub> + s<sub>3</sub> + s<sub>4</sub></td><td>10</td><td>1.5%</td><td>0.003<sup>0.005</sup><sub>0.001</sub></td></tr><tr><td>4</td><td>2g<sub>1</sub> - g<sub>4</sub> + g<sub>5</sub> - s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>8</td><td>1.3%</td><td>0.003<sup>0.006</sup><sub>0.001</sub></td></tr><tr><td>5</td><td>4g<sub>1</sub> - g<sub>2</sub> - g<sub>4</sub> - s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>10</td><td>1.4%</td><td>0.003<sup>0.004</sup><sub>0.001</sub></td></tr><tr><td>6</td><td>4g<sub>1</sub> - g<sub>2</sub> - g<sub>4</sub> - 3s<sub>2</sub> + s<sub>3</sub></td><td>10</td><td>1.3%</td><td>0.003<sup>0.005</sup><sub>0.001</sub></td></tr><tr><td>7</td><td>3g<sub>1</sub> - g<sub>4</sub> - 2s<sub>1</sub> - s<sub>2</sub> + s<sub>4</sub></td><td>8</td><td>1.8%</td><td>0.003<sup>0.004</sup><sub>0.001</sub></td></tr><tr><td>8</td><td>3g<sub>1</sub> - g<sub>3</sub> - 2s<sub>1</sub> - s<sub>2</sub> + s<sub>4</sub></td><td>8</td><td>1.6%</td><td>0.002<sup>0.006</sup><sub>0.000</sub></td></tr><tr><td>9</td><td>2g<sub>1</sub> - g<sub>3</sub> + g<sub>5</sub> - s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>8</td><td>1.0%</td><td>0.002<sup>0.004</sup><sub>0.000</sub></td></tr><tr><td>10</td><td>g<sub>1</sub> + g<sub>2</sub> - g<sub>3</sub> + g<sub>5</sub> - 2s<sub>1</sub> - s<sub>2</sub> + s<sub>3</sub></td><td>8</td><td>1.3%</td><td>0.002<sup>0.003</sup><sub>0.001</sub></td></tr><tr><td colspan="5"><i>C</i><sub>2</sub></td></tr><tr><td>1</td><td>g<sub>1</sub> - 3g<sub>2</sub> + 2g<sub>5</sub> - s<sub>1</sub> + 2s<sub>2</sub> - s<sub>4</sub></td><td>10</td><td>0.01‰</td><td>1e-4<sup>1e-4</sup></td></tr><tr><td>2</td><td>2g<sub>1</sub> - 4g<sub>2</sub> + 2g<sub>5</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>10</td><td>0.08‰</td><td>8e-5<sup>1e-4</sup></td></tr><tr><td>3</td><td>3g<sub>2</sub> - 3g<sub>5</sub> + s<sub>1</sub> - 2s<sub>2</sub> + s<sub>3</sub></td><td>10</td><td>0.01‰</td><td>3e-5<sup>3e-5</sup></td></tr><tr><td>4</td><td>g<sub>1</sub> - 4g<sub>2</sub> + g<sub>3</sub> + 2g<sub>5</sub> - s<sub>1</sub> + s<sub>2</sub></td><td>10</td><td>0.13‰</td><td>3e-5<sup>5e-5</sup></td></tr><tr><td>5</td><td>g<sub>1</sub> + 3g<sub>2</sub> - 4g<sub>5</sub> - s<sub>2</sub> + s<sub>4</sub></td><td>10</td><td>0.56‰</td><td>3e-5<sup>4e-5</sup></td></tr><tr><td>6</td><td>g<sub>1</sub> - 4g<sub>2</sub> + 3g<sub>5</sub> + s<sub>2</sub> - s<sub>4</sub></td><td>10</td><td>0.01‰</td><td>2e-5<sup>2e-5</sup></td></tr><tr><td>7</td><td>3g<sub>2</sub> - 3g<sub>5</sub> + s<sub>1</sub> - 2s<sub>2</sub> + s<sub>4</sub></td><td>10</td><td>0.13‰</td><td>1e-5<sup>3e-5</sup></td></tr><tr><td>8</td><td>4g<sub>2</sub> - g<sub>3</sub> - g<sub>4</sub> - 2g<sub>5</sub> + s<sub>1</sub> - s<sub>3</sub></td><td>10</td><td>0.03‰</td><td>1e-5<sup>5e-5</sup></td></tr><tr><td>9</td><td>2g<sub>1</sub> - 5g<sub>2</sub> + g<sub>4</sub> + 2g<sub>5</sub></td><td>10</td><td>0.02‰</td><td>9e-6<sup>9e-6</sup></td></tr><tr><td>10</td><td>2g<sub>1</sub> - 5g<sub>2</sub> + g<sub>3</sub> + 2g<sub>5</sub></td><td>10</td><td>0.14‰</td><td>6e-6<sup>8e-6</sup></td></tr></tbody></table>

<sup>a</sup>Only harmonics that are resonant for more than 1% of time are shown, i.e., $\tau^{\text{res}} > 1\%$.

of Ref. [19] by increasing the cutoff frequency of the low-pass filter applied to time series of the action-angle variables from $(5 \text{ Myr})^{-1}$ to $1 \text{ Myr}^{-1}$ (see Ref. [19], Appendixes F.2 and G.5). The filtered time series have then been resampled with a time step of 50 kyr. Many harmonics we show in Table IV and related to $C_2$ are resonant for a few time steps 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 $\mathbb{H}_{2n}$. The implications of Table IV suggest to introduce an additional truncation in the
Hamiltonian $\mathbb{H}_{2n}$. This consists in dropping the harmonics of Eq. (6) that involve external modes other than $g_5$:

$$ \mathbb{H}_{2n}^{\bullet}(\mathbf{I}, \boldsymbol{\theta}, t) = \sum_{k, \ell_1} \tilde{\mathbb{H}}_{2n}^{k, \ell}(\mathbf{I}) \mathrm{E}^{j[\mathbf{k}\boldsymbol{\theta}+\ell_1\boldsymbol{\phi}_1(t)]}, \quad (24) $$

where $\phi_1(t) = -g_5t$ and $\ell^* = (\ell_1, 0, \dots, 0)$, with $\ell_1 \in \mathbb{Z}$. Consistently with the absence of symmetry-breaking resonances related to $\mathcal{E}_{2n}$ in Table IV, the corresponding dynamics admits the exact integral of motion

$$ \mathcal{E}_{2n}^* = \mathbb{H}_{2n}^* + g_5 \sum_{i=1}^{4} (\mathfrak{X}_i + \Psi_i), \quad (25) $$

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 $\mathbb{H}_{2n}^*$ is initially very close to that of $\mathbb{H}_{2n}$. A frequency analysis over the first 20 Myr shows that the differences in the fundamental frequencies of the motion between $\mathbb{H}_{2n}^*$ and $\mathbb{H}_{2n}$ are of the order of $10^{-3}$ arcsec yr$^{-1}$, an order of magnitude smaller than the typical frequency differences between $\mathcal{H}_4$ and $\mathcal{H}_6$ (see Ref. [4], Table 3). Therefore, even though $\mathbb{H}_{2n}^*$ constitutes a simplification of $\mathbb{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 $\mathbb{H}_4^*$ in the same way as described in Sec. III in the case of $\mathcal{H}_{2n}$. Since its dynamics turns out to be much more stable than that of $\mathcal{H}_4$ (see Sec. 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 1(b). Comparing to the Lyapunov spectrum of $\mathcal{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, $\lambda_8$, decreases monotonically, consistently with the fact that $\mathcal{E}_4^*$ from Eq. (25) is an exact integral of motion. The exponent $\lambda_7$ drops by more than an order of magnitude, and apparently begins to stabilize around a few $10^{-4}$ arcsec yr$^{-1}$, while $\lambda_6$ also reduces significantly, by a factor of 3, 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

$$ \begin{aligned} 2\pi\lambda_6 &\approx \Delta\omega^{\mathcal{R}_3, \mathcal{L}_{\text{inc}}}, \\ 2\pi\lambda_7 &\approx \Delta\omega^{\mathcal{R}_3, \mathcal{L}_2}, \\ \lambda_8 &= 0, \end{aligned} \quad (26) $$

where $\Delta\omega^{\mathcal{R}_3, \mathcal{L}_2}$ stands for the half-width of the uppermost resonances of ranking $\mathcal{R}_3$ related to the quasi-integral $\mathcal{L}_2$. The hierarchy of the three smallest exponents in the spectrum of Fig. 1(b) 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:

$$ \begin{aligned} \lambda_6 &\leftrightarrow C_{\text{inc}}, \\ \lambda_7 &\leftrightarrow C_2, \\ \lambda_8 &\leftrightarrow \mathcal{E}_{2n}^*. \end{aligned} \quad (27) $$

These one-to-one correspondences are a particular case of Eq. (23) and support the physical intuition behind it. In Sec. V, we prove the validity of Eq. (27) in the geometric framework established by a principal component analysis of the orbital solutions of $\mathbb{H}_{2n}^*$.

**Numerical integrations.**—We compute ensembles of 1080 orbital solutions of the dynamical models $\mathbb{H}_4^*$ and $\mathbb{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 $\mathcal{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_{\text{inc}}$, while the latter is itself slower with respect to its counterpart in the orbital solutions of $\mathcal{H}_{2n}$. We remark that, as $\mathcal{E}_{2n}^*$ is an exact integral of motion for the model $\mathbb{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 $\mathbb{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 $\mathcal{H}_{2n}$ to $\mathbb{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 quasisymmetries 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 Sec. 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 low-pass 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., Refs. [52,53]). Let $I(t) = (\mathcal{X}(t), \Psi(t))$ the eight-dimensional time series of the action variables evaluated along a numerical solution of the equations of motion. As in Sec. IV B, we apply the Kolmogorov-Zurbenko low-pass filter with three iterations of the moving average and a cutoff frequency of $1 \text{ Myr}^{-1}$ to obtain the filtered time series $\hat{I}(t)$ [4,48]. In this way, the short-term quasiperiodic oscillations are mostly suppressed, which better reveals the chaotic diffusion over longer timescales. We finally define the mean-subtracted filtered action variables over the time interval $[t_0, t_0 + T]$ as $\tilde{I}(t) = \hat{I}(t) - n^{-1} \sum_{i=0}^{n-1} \hat{I}(t_0 + i\Delta t)$, where the mean is estimated by discretization of the time series with a sampling step $\Delta t$ such that $T = (n - 1)\Delta t$. The discretized time series over the given interval is stored in an $8 \times n$ matrix:

$$ \mathbf{D} = [\tilde{I}(t_0), \tilde{I}(t_0 + \Delta t), \dots, \tilde{I}(t_0 + (n-1)\Delta t)]. \quad (28) $$

The PCA of the data matrix $\mathbf{D}$ consists in a linear transformation $\mathbf{P} = \mathbf{A}^\mathrm{T}\mathbf{D}$, where $\mathbf{A}$ is an $8 \times 8$ orthogonal matrix (i.e., $\mathbf{A}^{-1} = \mathbf{A}^\mathrm{T}$) defined as follows. By writing $\mathbf{A} = [\mathbf{a}_1, \dots, \mathbf{a}_8]$, the column vectors $\mathbf{a}_i \in \mathbb{R}^8$ are chosen to be the normalized eigenvectors of the sample covariance matrix, in order of decreasing eigenvalues: $(n-1)^{-1}\mathbf{D}\mathbf{D}^\mathrm{T} = \mathbf{A}\mathbf{\Sigma}\mathbf{A}^\mathrm{T}$, where $\mathbf{\Sigma} = \text{diag}(\sigma_1, \dots, \sigma_8)$ and $\sigma_1 \ge \dots \ge \sigma_8$.

The PCs are defined as the new variables after the transformation; that is, $\mathbf{PC}_i = \mathbf{a}_i \cdot \mathbf{I}$, with $i \in \{1, \dots, 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}\mathbf{PP}^\mathrm{T} = \mathbf{\Sigma}$, from which it follows that the variance of $\mathbf{PC}_i$ is $\sigma_i$.

Among all the linear combinations in the action variables $I$, the last PC, i.e., $\mathbf{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., $\mathbf{PC}_7$, has the second smallest variance and is uncorrelated with $\mathbf{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: $\mathbf{D}_{\text{agg}} = [\mathbf{D}_1, \mathbf{D}_2, \dots, \mathbf{D}_m]$, where $\mathbf{D}_i$ is the data matrix of Eq. (28) for the $i$th solution. Since this work deals with a nonstationary 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 semi-analytically in Sec. 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.

### 1. Tweaking the QIs

The three last components $\mathbf{PC}_8, \mathbf{PC}_7, \mathbf{PC}_6$ are represented by the set of vectors $S_{\mathbf{PC}_8} = \{\mathbf{a}_8, \mathbf{a}_7, \mathbf{a}_6\}$ belonging to $\mathbb{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 $\le$ the order of statistical variance, one has $PC_8 \le PC_7 \le PC_6$; the unit vectors $(\boldsymbol{a}_i)_{i=6}^8$ are orthogonal to each other. On the other hand, the QIs of motion $C_{\text{inc}}, C_2, \mathcal{E}_{2n}$ do not possess these properties. Therefore, we adjust them in such a way to reproduce the same structure.

(1) *Linearity*.—While $C_{\text{inc}}$ and $C_2$ are linear functions of the action variables, $\mathcal{E}_{2n}$ is not when $n > 1$. Nevertheless, as we explain in Sec. IVA, as far as one considers stable orbital solutions, the linear LL approximation $\mathcal{E}_2 = \gamma_3 \cdot I$ reproduces $\mathcal{E}_{2n}$ reasonably well. Therefore, we consider the three linear QIs of motion $C_{\text{inc}}, C_2, \mathcal{E}_2$, which are respectively represented by the set of $\mathbb{R}^8$ vectors $S_{\text{QIs}} = \{\gamma_1, \gamma_2, \gamma_3\}$. In this way, the three-dimensional linear subspaces of the action space spanned by the sets $S_{\text{QIs}}$ and $S_{\text{PCs}}$ can be compared.

(2) *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 $\mathbb{H}_{2n}^*$ in Eq. (24) or $\mathbb{H}_{2n}$ in Eq. (6) (clearly $n > 1$).

(a) [$\mathbb{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 $\mathcal{E}_{2n}^* < C_2 < C_{\text{inc}}$. While $\mathcal{E}_{2n}$ is an exact nonlinear integral of motion, we expect that its linear truncation $\mathcal{E}_2^* = \mathcal{E}_2$ varies more than $C_2$ and $C_{\text{inc}}$. Therefore, we consider the ordered set of QIs of motion $\{C_2, C_{\text{inc}}, \mathcal{E}_2\}$ represented by the ordered set of vectors $S'_{\text{QIs}} = \{\gamma_2, \gamma_1, \gamma_3\}$.

(b) [$\mathbb{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'_{\text{QIs}}$ through the projections of the three last PCs onto the linear subspace generated by the QIs: $S'_{\text{QIs}} = \{\text{proj}_{S_{\text{QIs}}}(\boldsymbol{a}_8), \text{proj}_{S_{\text{QIs}}}(\boldsymbol{a}_7), \text{proj}_{S_{\text{QIs}}}(\boldsymbol{a}_6)\}$ [54]. As a result, the new set of QIs mirrors the hierarchical structure of the PCs. We stress that $S'_{\text{QI}}$ spans the same subspace of $\mathbb{R}^8$ as $S_{\text{QI}}$, since the ordered QIs are just linear combinations of the original ones.

(3) *Orthogonality*.—We apply the Gram-Schmidt process to the ordered set $S'_{\text{QIs}}$ to obtain the orthonormal basis $S''_{\text{QIs}} = \{\boldsymbol{\alpha}_1, \boldsymbol{\alpha}_2, \boldsymbol{\alpha}_3\}$. The set $S''_{\text{QIs}}$ clearly spans the same subspace as $S_{\text{QIs}}$. Moreover, the Gram-Schmidt process preserves the hierarchical structure;

that is, the two $m$-dimensional subspaces spanned by the first $m \le 3$ vectors of $S'_{\text{QIs}}$ and $S''_{\text{QIs}}$, respectively, are identical.

In the end, we obtain a linear, ordered, and orthogonal set of modified QIs of motion $\{\text{QI}_1, \text{QI}_2, \text{QI}_3\}$, where $\text{QI}_i = \boldsymbol{\alpha}_i \cdot \mathbf{I}$.

## 2. Visual comparison

We now visually compare the vectors $\boldsymbol{\alpha}_{1,2,3}$ of the modified QIs with the corresponding vectors $\boldsymbol{a}_{8,7,6}$ of the last three PCs. We use the ensembles of 1080 numerically integrated orbital solutions of the models $\mathcal{H}_4$ and $\mathbb{H}_4$ considered in Secs. IV B and IV D, respectively. The nominal solution of each set is denoted as sol. 1 from now on. For the model $\mathcal{H}_4$, we also consider two other solutions: sol. 2, which 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 Sec. VB 3.

**Hamiltonian $\mathbb{H}_4^*$.—** The modified QIs can be explicitly derived in this case and comprise interpretable physical quantities. One has $\text{QI}_1$ proportional to $C_2$ and $\text{QI}_2$ proportional to $C_2^\perp$. Moreover, $\text{QI}_3$ is the component of $\mathcal{E}_2$ that is orthogonal to both $C_2$ and $C_2^\perp$. Figure 3 shows the comparison between the modified QIs and the corresponding PCs for three different time intervals along sol. 1 of $\mathbb{H}_4^*$ (see Fig. 2, bottom left for its time evolution). The agreement of the pairs $(\text{QI}_1, \text{PC}_8)$, $(\text{QI}_2, \text{PC}_7)$, and $(\text{QI}_3, \text{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 three-dimensional linear subspace of the action space. Therefore, to a linear approximation, they represent the three slowest d.o.f. of the $\mathbb{H}_4^*$ dynamics. The quasi-integral $C_2$ represents the slowest linear d.o.f.: it coincides with the last principal component $\text{PC}_8$, which has the smallest variance among all the linear combinations of the action variables. $C_{\text{inc}}$ and $\mathcal{E}_2$ represent the second and the third slowest linear d.o.f., respectively: the component of $C_{\text{inc}}$ orthogonal to $C_2$, i.e., $C_2^\perp$, matches the second last principal component $\text{PC}_7$; the component of $\mathcal{E}_2$ orthogonal to the subspace generated by $(C_2, C_{\text{inc}})$ matches the third last principal component $\text{PC}_6$. The strong hierarchical structure of the slow variables for the simplified dynamics $\mathbb{H}_4^*$ is clearly confirmed by the almost frozen basis vectors of the PCs.

**Hamiltonian $\mathcal{H}_4$.—** In this case, the QIs of motion $C_{\text{inc}}, C_2, \mathcal{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 three-dimensional subspaces of $\mathbb{R}^8$, we compare their basis vectors instead. The basis $\boldsymbol{\alpha}_{1,2,3}$ of
<figure><img src="image_6.png" /><figcaption>FIG. 3. Vectors $\alpha_{1,2,3}$ representing the three modified QIs ($QI_{1,2,3}$, black circles) compared to the corresponding vectors $\alpha_{8,7,6}$ of the three last PCs ($PC_{8,7,6}$, red dots), for the intervals [0, 500] Myr (left-hand column), [1000, 2000] Myr (middle column), and [0, 5000] Myr (right-hand column) of sol. 1 and of the aggregated sample of 1080 solutions of model $\mathcal{H}_4^*$. Here, $QI_1$ is proportional to $C_2$ and $QI_2$ is proportional to $C_2^\perp$; see Eq. (20).</figcaption></figure>

modified quasi-integrals $QI_{1,2,3}$ is computed according to the algorithm presented in Sec. 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 $\mathcal{H}_4$ (see Fig. 5 for their time evolution). The first two, sol. 1 and sol. 2, show thorough agreement between the pairs of QIs and PCs across all intervals, which indicates close proximity between the two subspaces $V_{QIs} = \text{span}(S_{QIs})$ and $V_{PCs} = \text{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 $\mathcal{H}_4$ can thus be deduced to be close to $C_{inc}$, in line with the discussion in Sec. 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 $\mathcal{H}_4^*$, there is no pronounced separation between the slowest d.o.f. at the bottom of the Lyapunov spectrum in Fig. 1(a): the marginal distributions of consecutive exponents can indeed touch or overlap each other. Therefore, the hierarchy of slow variables is not as

<figure><img src="image_7.png" /><figcaption>FIG. 4. Vectors $\alpha_{1,2,3}$ representing the three modified QIs ($QI_{1,2,3}$, black circles) compared to the corresponding vectors $\alpha_{8,7,6}$ of the three last PCs ($PC_{8,7,6}$, red dots), for the intervals [0, 500] Myr (left-hand column), [1000, 2000] Myr (middle column), and [0, 5000] Myr (right-hand column) of sol. 1, sol. 2, and sol. 3 and of the aggregated sample of 1080 solutions of model $\mathcal{H}_4$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 5. Time evolution over 5 Gyr of the dimensionless QIs of motions ($\hat{C}_{\text{inc}}$, $\hat{C}_2$, $\hat{\mathcal{E}}$) and of two representatives of the dimensionless action variables ($\hat{\mathbf{x}}_1$, $\hat{\mathbf{\Psi}}_3$) for three solutions of $\mathcal{H}_4$, that is, sol. 1 (top), sol. 2 (middle), and sol. 3 (bottom). $\hat{\mathcal{E}}$ stands for $\hat{\mathcal{E}}_4$. The time series are low-pass filtered with a cutoff frequency of $1 \text{ Myr}^{-1}$ and the mean over 5 Gyr is subtracted.</figcaption></figure>

frozen as in $\mathbb{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 $\text{PC}_8$ and $C_{\text{inc}}$ (see Fig. 4). Once again, the PCA confirms that the subspace spanned by the three QIs is overall close to the slowest three-dimensional linear subspace of the action space. Therefore, to a linear approximation, they represent the three slowest d.o.f. of the $\mathcal{H}_4$ dynamics. We remark that the slowness of the three-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 = \hat{\mathbf{q}} \cdot \mathbf{I}$ be a slow variable with unit vector $\hat{\mathbf{q}}$. If $\boldsymbol{\epsilon}$ is an arbitrary small vector, i.e., $\|\boldsymbol{\epsilon}\| \ll 1$, then $Q' = (\hat{\mathbf{q}} + \boldsymbol{\epsilon}) \cdot \mathbf{I}$ can also be considered as a slow variable, whereas the normalized difference of two quantities, $\hat{\boldsymbol{\epsilon}} \cdot \mathbf{I}$, is generally not. Therefore, the linear subspace spanned by $Q$ and $Q'$, that is, by $\hat{\mathbf{q}}$ and $\hat{\boldsymbol{\epsilon}}$, is not a slow two-dimensional subspace.

Solution 3 in Fig. 4 represents an edge case (see Fig. 5 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, $\Psi_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 $\text{QI}_3$ and $\text{PC}_6$ on the last two time intervals of sol. 3 in Fig. 4. We remark that $\text{PC}_{8,7}$ are still relatively close to $\text{QI}_{1,2}$, which indicates that the slowest two-dimensional subspace spanned by $\text{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. 1(a) shows that the marginal distributions of the exponents $\lambda_5$ and $\lambda_6$, for example, are not separate but overlap each other.

### 3. Distance between the subspaces of PCs and QIs

The closeness of the two three-dimensional linear subspaces $V_{\text{PCs}}$, $V_{\text{QIs}} \subset \mathbb{R}^8$ spanned by the sets of vectors $S_{\text{PCs}}$ and $S_{\text{QIs}}$, respectively, can be quantitatively measured in terms of a geometric distance. This can be formulated using the principal (canonical) angles [55–57].
Let A and B be two sets of $m \le n$ independent vectors in $\mathbb{R}^n$. The principal vectors $(p_k, q_k)_{k=1}^m$ are defined recursively as solutions to the optimization problem:

$$ \begin{array}{ll} \text{maximize} & p \cdot q \\ \text{subject to} & \begin{array}{l} p \in \operatorname{span}(A), q \in \operatorname{span}(B), \\ \|p\| = 1, \|q\| = 1, \\ p \cdot p_i = 0, q \cdot q_i = 0, \quad i = 1, \dots, k-1, \end{array} \end{array} \quad (29) $$

for $k = 1, \dots, m$. The principal angles $0 \le \theta_1 \le \dots \le \theta_m \le \pi/2$ between the two subspaces $\operatorname{span}(A)$ and $\operatorname{span}(B)$ are then defined by

$$ \cos \theta_k = p_k \cdot q_k, \quad k = 1, \dots, m. \quad (30) $$

The principal angle $\theta_1$ is the smallest angle between all pairs of unit vectors in $\operatorname{span}(A)$ and $\operatorname{span}(B)$; the principal angle $\theta_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 $\operatorname{diag}(\cos \theta_1, \dots, \cos \theta_m)$. This cosine-based method is often ill conditioned for 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]:

$$ d(A, B) = \left( \frac{1}{m} \sum_{k=1}^{m} \sin^2 \theta_k \right)^{1/2}. \quad (31) $$

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_{\text{PCs}}$ and $V_{\text{PCs}}$, arising from the $\mathbb{H}_4^*$ and $\mathcal{H}_4$ dynamics, coincides with that of randomly drawn three-dimensional subspaces of $\mathbb{R}^8$. The PDF of the distance between two random three-dimensional subspaces of $\mathbb{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 three-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.

<figure><img src="image_9.png" /><figcaption>FIG. 6. PDF of the distance between two random three-dimensional linear subspaces of $\mathbb{R}^8$ (blue, $10^5$ draws) compared with the PDF of the distance between the two subspaces $V_{\text{PCs}}$ ($PC_{8,7,6}$) and $V_{\text{QIs}}$ ($QI_{1,2,3}$) arising from the time interval [0, 5] Gyr of 1080 solutions of $\mathbb{H}_4^*$ (top) and 10 800 solutions of $\mathcal{H}_4$ (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).</figcaption></figure>

The upper panel of Fig. 6 shows in green the PDF of the distance between $V_{\text{PCs}}$ and $V_{\text{QIs}}$ arising from the time interval [0, 5] Gyr of the 1080 orbital solutions of model $\mathbb{H}_4^*$. In the lower panel, we consider a larger ensemble of 10 800 solutions of model $\mathcal{H}_4$ spanning the same time interval [7], and plot the corresponding PDF of the distance between $V_{\text{PCs}}$ and $V_{\text{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_{\text{PCs}}$ and $V_{\text{QIs}}$ is a statistically robust result. In the case of the simplified dynamics $\mathbb{H}_4^*$, the PDF peaks around a median of roughly 0.08 and has small variance. Switching to model $\mathcal{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 Sec. V B 2: a quasifrozen 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 constrain 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_4^*$ and $H_6^*$ introduced in Sec. IV D. We recall 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_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_6^*$. The dynamics arising from $H_4^*$ and $H_6^*$ can be considered as stable in an astronomical sense. Recalling that the main difference between $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 $\mathcal{E}_2$ allows us to bound the variations of the action variables $\mathcal{X}, \mathcal{P}$. 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 $\mathcal{P}$. As a consequence, at degree 2 in eccentricities and inclinations, fixing a value of

$$ \mathcal{E}_2 = \gamma_3 \cdot \mathbf{I} = \gamma_3^{\text{ecc}} \cdot \mathcal{X} + \gamma_3^{\text{inc}} \cdot \mathcal{P}, \quad (32) $$

with $\gamma_3 = (\gamma_3^{\text{ecc}}, \gamma_3^{\text{inc}})$, also bounds the upper variations of the eccentricity actions $\mathcal{X}$, since the components of $\gamma_3^{\text{ecc}}$ have all the same sign, as those of $\gamma_3^{\text{inc}}$ (see Appendix B). This is an important point, as we state in Sec. I that the lack

<figure><img src="image_10.png" /><figcaption>FIG. 7. Cumulative distribution function of the first time that Mercury eccentricity reaches a value of 0.7, from 1080 orbital solutions of different models over 100 Gyr. The shaded regions represent the 90% piecewise confidence intervals from bootstrap.</figcaption></figure>

of any bound on the chaotic variations of the planet orbits is one of the reasons that complicates 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 considered 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 us 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 quasisymmetries 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'' \text{ 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 quasisymmetries generate three QIs of motion that represent slow variables of the secular dynamics. The size of the leading symmetry-breaking 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–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–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–72]. A hierarchy of Lyapunov exponents spanning several orders of magnitude, in particular, should be common among this class of systems (see, e.g., Ref. [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 [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 semianalytically 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 [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.

## ACKNOWLEDGMENTS

The authors thank M. Gastineau for his assistance with TRIP. F. M. is supported by a grant of the French Agence Nationale de la Recherche (AstroMeso ANR-19-CE31-0002-01). N. H. H. is supported by a Ph.D. scholarship of the Capital Fund Management (CFM) Foundation for Research. 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 program Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

## APPENDIX A: LYAPUNOV SPECTRUM

*Convergence.*—We perform two tests to address the convergence of our implementation of the Benettin *et al.* [43] method. We first compute the FT-LCEs for a single initial condition of $H_4$ and an ensemble of 150 different random sets of initial tangent vectors. Figure 8(a) shows the [5th, 95th] percentile range of the resulting marginal distributions of the positive FT-LCEs over a time span of 10 Gyr. The distributions shrink with increasing time, eventually collapsing on single time-dependent values. In this asymptotic regime, the Benettin *et al.* [43] algorithm loses memory of the initial tangent vectors and purely retrieves the FT-LCEs as defined in Eq. (10). Therefore, Fig. 1(a) shows asymptotically the dependence of the FT-LCEs on the initial condition $z_0$ and represents their statistical distribution over the phase-space domain

<figure><img src="image_11.png" /><figcaption>FIG. 8. (a) Positive FT-LCEs of Hamiltonian $H_4$ 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 $\epsilon_i$ on the FT-LCEs $\lambda_i$, as defined in Eq. (A1), for the ensemble of 150 orbital solution of Fig. 1(a).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 9. Positive FT-LCEs $\lambda_i$ of Hamiltonian $\mathcal{H}_6$ and corresponding characteristic timescales $\lambda_i^{-1}$. 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.</figcaption></figure>

explored by the dynamics in a nonergodic way. The convergence of the computation is clearly slower for smaller exponents, but a comparison with Fig. 1(a) indicates that, even in the case of $\lambda_8$, the numerical uncertainty on the FT-LCEs of each orbital solution at 10 Gyr is negligible with respect to the width of their ensemble distributions.

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 $\lambda_i$ can be estimated as

$$ \epsilon_i = \left| \frac{\Delta\lambda_i}{\lambda_i} \right| \qquad (A1) $$

We plot in Fig. 8(b) the medians of $\epsilon_i$ for the ensemble of 150 orbital solutions of Fig. 1(a). The relative errors decrease asymptotically with time, as expected. Even in the case of the smallest exponent, $\lambda_8$, the median error is less than 10% at 10 Gyr.

Hamiltonian $\mathcal{H}_6$.—We compute for comparison the FT-LCEs of the forced ISS truncated at degree 6 in eccentricities and inclinations, that is, $\mathcal{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 $\mathcal{H}_4$ shown in Fig. 1(a).

## APPENDIX B: VECTORS $\gamma_1, \gamma_2, \gamma_3$

We report here the explicit expressions of the vectors $(\gamma_i)_{i=1}^3$. We first give the components of the vector $\omega_{LL}$ of the fundamental precession frequencies of the inner orbits in the forced Laplace-Lagrange dynamics (including the leading correction of general relativity) [4]:

$$ \begin{aligned} \omega_{LL} &= (g_{LL}, s_{LL}) \\ &\approx (5.87, 7.46, 17.4, 18.1, -5.21, -6.59, -18.8, -17.7), \end{aligned} \qquad (B1) $$

in units of arcsec yr$^{-1}$ (see Refs. [80–82] for comparison with the frequencies of the Laplace-Lagrange dynamics of the entire Solar System). One then has

$$ \begin{aligned} \gamma_1 &= (\mathbf{0}_4, \mathbf{1}_4) = (0, 0, 0, 0, 1, 1, 1, 1), \\ \gamma_2 &= (0, 0, -1, -1, 1, 1, 2, 2), \\ \gamma_3 &= -\omega_{LL} + g_5 \mathbf{1}_8 \\ &\approx (-1.61, -3.20, -13.2, -13.9, 9.47, 10.8, 23.0, 22.0), \end{aligned} \qquad (B2) $$

with the components of $\gamma_3$ in units of arcsec yr$^{-1}$. We recall that $g_5 \approx 4.257''$ yr$^{-1}$ is a constant in the forced model of the ISS. The corresponding unit vectors $(\hat{\gamma}_i)_{i=1}^3$ are given by

$$ \begin{aligned} \hat{\gamma}_1 &= (0, 0, 0, 0, 1, 1, 1, 1)/2, \\ \hat{\gamma}_2 &= (0, 0, -1, -1, 1, 1, 2, 2)/2\sqrt{3}, \\ \hat{\gamma}_3 &\approx (-0.04, -0.08, -0.33, -0.35, 0.24, 0.27, 0.58, 0.55). \end{aligned} \qquad (B3) $$

Since $1/2\sqrt{3} \approx 0.289$, the components of $\hat{\gamma}_3$ are only a few percent away from those of $\hat{\gamma}_2$. Therefore, along stable orbital solutions with typical bounded variations of the Mercury-dominated action variable $\mathcal{X}_1$, the two quantities $C_2$ and $\mathcal{E}_{2n}$ exhibit very similar time evolutions. This is not the case anymore when Mercury orbit reaches high eccentricities.

## 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 $\mathcal{H}_4$ and $\mathcal{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 $\mathcal{H}_4^*$ and $\mathcal{H}_6^*$ that we introduce in Sec. IV D. We report in Fig. 10 the time evolution of the ensemble PDFs of the low-pass filtered dimensionless QIs and dimensionless actions $\mathcal{X}_1$,
<figure><img src="image_13.png" /><figcaption>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 $\mathcal{X}_1$, $\Psi_3$. Estimation from an ensemble of 1080 numerical orbital solutions for different models ($\mathcal{H}_4$, $\mathcal{H}_6$, $\mathcal{H}_4^\bullet$, and $\mathcal{H}_6^\bullet$). First row: $\hat{C}_{\text{inc}}$. Second row: $\hat{C}_2$. Third row: $\hat{\mathcal{E}}_4(\mathcal{H}_4)$ and $\hat{\mathcal{E}}_6(\mathcal{H}_6)$. Fourth row: $\hat{\mathcal{X}}_1$. Fifth row: $\hat{\Psi}_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 $\hat{\mathcal{E}}_{2n}^\bullet$ is an exact integral of motion for the model $\mathcal{H}_{2n}^\bullet$ and its PDF has null dispersion.</figcaption></figure>

$\Psi_3$ for the different models (the cutoff frequency of the time filter is set to $1 \text{ 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

<figure><img src="image_14.png" /><figcaption>FIG. 11. Time evolution of the interquartile range (IQR) of the ensemble PDFs of the QIs shown in Fig. 10. Left: $\hat{C}_{\text{inc}}$. Middle: $\hat{C}_2$. Right: $\hat{\mathcal{E}}_4(\mathcal{H}_4)$ and $\hat{\mathcal{E}}_6(\mathcal{H}_6)$. The quantity $\hat{\mathcal{E}}_{2n}^\bullet$ is an exact integral of motion for the model $\mathcal{H}_{2n}^\bullet$ and its PDF has a null IQR.</figcaption></figure>
the exponential separation of close trajectories, the time growth of the IQR follows a power law typical of diffusion processes. Figures 10 and 11 clearly show the slower diffusion of $C_{\text{inc}}$ and $C_2$ in the model $H_{2n}^*$ when compared to $H_{2n}$. We recall that $\mathcal{E}_{2n}^*$ is an exact integral of motion for the model $H_{2n}^*$ (see Sec. IV D) and its PDF has null dispersion.

[1] J. Laskar, *A Numerical Experiment on the Chaotic Behaviour of the Solar System*, Nature (London) **338**, 237 (1989).

[2] J. Laskar, *The Chaotic Motion of the Solar System—A Numerical Estimate of the Size of the Chaotic Zones*, Icarus **88**, 266 (1990).

[3] G. J. Sussman and J. Wisdom, *Chaotic Evolution of the Solar System*, Science **257**, 56 (1992).

[4] F. Mogavero and J. Laskar, *Long-Term Dynamics of the Inner Planets in the Solar System*, Astron. Astrophys. **655**, A1 (2021).

[5] J. Laskar, *Chaotic Diffusion in the Solar System*, Icarus **196**, 1 (2008).

[6] J. Laskar and M. Gastineau, *Existence of Collisional Trajectories of Mercury, Mars and Venus with the Earth, Nature (London)* **459**, 817 (2009).

[7] N. H. Hoang, F. Mogavero, and J. Laskar, *Long-Term Instability of the Inner Solar System: Numerical Experiments*, Mon. Not. R. Astron. Soc. **514**, 1342 (2022).

[8] K. Batygin, A. Morbidelli, and M. J. Holman, *Chaotic Disintegration of the Inner Solar System*, Astrophys. J. **799**, 120 (2015).

[9] E. Woillez and F. Bouchet, *Instantons for the Destabilization of the Inner Solar System*, Phys. Rev. Lett. **125**, 021101 (2020).

[10] The dynamics truncated at degree 4 produces nevertheless the same chaos of the full system, as measured by the finite-time maximum Lyapunov exponent [19].

[11] A. Milani and A. M. Nobili, *An Example of Stable Chaos in the Solar System*, Nature (London) **357**, 569 (1992).

[12] A. Morbidelli and C. Froeschlé, *On the Relationship between Lyapunov Times and Macroscopic Instability Times*, Celest. Mech. Dyn. Astron. **63**, 227 (1996).

[13] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou, *Studies of Nonlinear Problems, Technical Report, Los Alamos National Laboratory*, 1955.

[14] K.-D. N. T. Lam and J. Kurchan, *Stochastic Perturbation of Integrable Systems: A Window to Weakly Chaotic Systems*, J. Stat. Phys. **156**, 619 (2014).

[15] T. Goldfriend and J. Kurchan, *Equilibration of Quasi-integrable Systems*, Phys. Rev. E **99**, 022146 (2019).

[16] This resonance involves the fundamental precession frequencies of Mercury and Jupiter perihelia [5,36,62].

[17] J. Laskar, *Large Scale Chaos and Marginal Stability in the Solar System*, Celest. Mech. Dyn. Astron. **64**, 115 (1996).

[18] “At each stage of its evolution, the system should have a time of stability comparable with its age” [17].

[19] F. Mogavero and J. Laskar, *The Origin of Chaos in the Solar System through Computer Algebra*, Astron. Astrophys. **662**, L3 (2022).

[20] J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A. C. M. Correia, and B. Levrard, *A Long-Term Numerical Solution for the Insolation Quantities of the Earth*, Astron. Astrophys. **428**, 261 (2004).

[21] A. Morbidelli, *Modern Celestial Mechanics: Aspects of Solar System Dynamics* (Taylor & Francis, London, 2002).

[22] C. F. Gauss, *Determinatio attrationis, quam in punctum quodvis positionis datae exerceret planeta, SI eius massa per tota orbitam, ratione temporis, quo singulae partes describuntur, uniformiter errat dispertita* (Heinrich Dieterich, Göttingen, 1818).

[23] J. Laskar, *Analytical Framework in Poincaré Variables for the Motion of the Solar System*, in *Predictability, Stability, and Chaos in N-Body Dynamical Systems*, edited by A. E. Roy, NATO Advanced Study Institutes, Series B, Vol. 272 (Plenum Press, New York, 1991), pp. 93–114, 10.1007/978-1-4684-5997-5_7.

[24] J. Laskar and P. Robutel, *Stability of the Planetary Three-Body Problem. I. Expansion of the Planetary Hamiltonian*, Celest. Mech. Dyn. Astron. **62**, 193 (1995).

[25] M. Gastineau and J. Laskar, *TRIP: A Computer Algebra System Dedicated to Celestial Mechanics and Perturbation Series*, ACM Commun. Comput. Algebr. **44**, 194 (2011).

[26] M. Gastineau and J. Laskar, *TRIP 1.4.120, TRIP Reference Manual* (IMCCE, Paris Observatory, 2021).

[27] E represents the exponential operator and $j$ stands for the imaginary unit. The overline on variables denotes complex conjugate.

[28] N. H. Hoang, F. Mogavero, and J. Laskar, *Chaotic Diffusion of the Fundamental Frequencies in the Solar System*, Astron. Astrophys. **654**, A156 (2021).

[29] J. Laskar, *Secular Evolution of the Solar System over 10 Million Years*, Astron. Astrophys. **198**, 341 (1988).

[30] J. Laskar, *Frequency Map Analysis and Quasiperiodic Decompositions*, in *Hamiltonian Systems and Fourier Analysis: New Prospects For Gravitational Dynamics*, edited by D. Benest, C. Froeschlé, and E. Lega (Cambridge Scientific Publishers Ltd., Cambridge, 2005), pp. 93–114.

[31] V. I. Oseledec, *A Multiplicative Ergodic Theorem. Lyapunov Characteristic Numbers for Dynamical Systems*, Trans. Moscow Math. Soc. **19**, 197 (1968).

[32] J. P. Eckmann and D. Ruelle, *Ergodic Theory of Chaos and Strange Attractors*, Rev. Mod. Phys. **57**, 617 (1985).

[33] P. Gaspard, *Chaos, Scattering and Statistical Mechanics*, Cambridge Nonlinear Science Series (Cambridge University Press, Cambridge, England, 1998).

[34] C. Skokos, *The Lyapunov Characteristic Exponents and Their Computation*, in *Dynamics of Small Solar System Bodies and Exoplanets*, edited by J. J. Souchay and R. Dvorak (Springer, Berlin, 2010), p. 63–135.

[35] J. Laskar, *Large-Scale Chaos in the Solar System*, Astron. Astrophys. **287**, L9 (1994).

[36] K. Batygin and G. Laughlin, *On the Dynamical Stability of the Solar System*, Astrophys. J. **683**, 1207 (2008).

[37] The FT-LCEs do not depend on any tangent vector $\xi_0$.

[38] All the initial conditions mentioned in this work are sampled from a multidimensional Gaussian distribution that is centered at the nominal initial conditions of Gauss's dynamics (see Ref. [4], Appendix D) and has a relative width of $10^{-9}$ (see Ref. [19], Appendix C, and Ref. [7], Sec. III, for details).
[39] A. Morbidelli and A. Giorgilli, *On a Connection between KAM and Nekhoroshev's Theorems*, Physica (Amsterdam) **86D**, 514 (1995).

[40] J. Laskar, *A Few Points on the Stability of the Solar System*, in *Chaos, Resonance, and Collective Dynamical Phenomena in the Solar System*, edited by S. Ferraz-Mello (Kluwer Academic Publishers, Dordrecht, 1992), Vol. 152, p. 1.

[41] G. Brown and H. Rein, *A Repository of Vanilla Long-Term Integrations of the Solar System*, Res. Notes Am. Astron. Soc. **4**, 221 (2020).

[42] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, *All Lyapunov Characteristic Numbers Are Effectively Computable*, C.R. Seances Acad. Sci., Ser. B **286**, 431 (1978).

[43] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, *Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems—A Method for Computing All of Them. I—Theory. II—Numerical Application*, Meccanica **15**, 9 (1980). **15**, 21 (1980).

[44] B. V. Chirikov, *A Universal Instability of Many-Dimensional Oscillator Systems*, Phys. Rep. **52**, 263 (1979).

[45] J. Laskar, *Systèmes de Variables et Éléments*, in *Modern Methods in Celestial Mechanics*, edited by D. Benest and C. Froeschle (Editions Frontières, Gif-Sur-Yvette, France, 1990), pp. 63–87.

[46] G. Boué and J. Laskar, *Spin Axis Evolution of Two Interacting Bodies*, Icarus **201**, 750 (2009).

[47] J. Laskar, *Large Scale Chaos and the Spacing of the Inner Planets*, Astron. Astrophys. **317**, L75 (1997).

[48] I. G. Zurbenko and D. Smith, *Kolmogorov–Zurbenko Filters in Spatiotemporal Analysis*, Wiley Interdiscip. Rev. Comput. Stat. **10**, e1419 (2018).

[49] A. C. M. Correia and J. Laskar, *Long-Term Evolution of the Spin of Venus. II. Numerical Simulations*, Icarus **163**, 24 (2003).

[50] K. Pearson, *LIII. On Lines and Planes of Closest Fit to Systems of Points in Space*, Phil. Trans. R. Soc. A **2**, 559 (1901).

[51] H. Hotelling, *Analysis of a Complex of Statistical Variables into Principal Components*, J. Educ. Psychol. **24**, 417 (1933).

[52] I. T. Jolliffe, *Principal Component Analysis* (Springer, New York, 2002).

[53] I. T. Jolliffe and J. Cadima, *Principal Component Analysis: A Review and Recent Developments*, Phil. Trans. R. Soc. A **374**, 20150202 (2016).

[54] The projection of a vector $q$ onto the subspace spanned by a set of vectors $S$ can be written in a vectorial form as $B(B^T B)^{-1} B^T q$, where the column space of the matrix $B$ is the subspace spanned by the set $S$.

[55] C. Jordan, *Essai sur la Géométrie à n Dimensions*, Bull. Soc. Math. Fr. **3**, 103 (1875).

[56] C. F. Van Loan and G. Golub, *Matrix Computations* (The Johns Hopkins University Press, London, 1996).

[57] K. Ye and L.-H. Lim, *Schubert Varieties and Distances between Subspaces of Different Dimensions*, SIAM J. Matrix Anal. Appl. **37**, 1176 (2016).

[58] A. Björck and G. H. Golub, *Numerical Methods for Computing Angles between Linear Subspaces*, Math. Comput. **27**, 579 (1973).

[59] A. V. Knyazev and M. E. Argentati, *Principal Angles between Subspaces in an A-Based Scalar Product: Algorithms and Perturbation Estimates*, SIAM J. Sci. Comput. **23**, 2008 (2002).

[60] M. E. Muller, *A Note on a Method for Generating Points Uniformly on n-Dimensional Spheres*, Commun. ACM **2**, 19 (1959).

[61] F. Mogavero, *Addressing the Statistical Mechanics of Planet Orbits in the Solar System*, Astron. Astrophys. **606**, A79 (2017).

[62] G. Boué, J. Laskar, and F. Farago, *A Simple Model of the Chaotic Eccentricity of Mercury*, Astron. Astrophys. **548**, A43 (2012).

[63] M. Hénon, *Integrals of the Toda Lattice*, Phys. Rev. B **9**, 1921 (1974).

[64] H. Flaschka, *The Toda Lattice. II. Existence of Integrals*, Phys. Rev. B **9**, 1924 (1974).

[65] S. V. Manakov, *Complete Integrability and Stochastization of Discrete Dynamical Systems*, Zh. Eksp. Teor. Fiz. **67**, 543 (1974) [Sov. J. Exp. Theor. Phys. **40**, 269 (1975)].

[66] J. Ferguson, W. E., H. Flaschka, and D. W. McLaughlin, *Nonlinear Normal Modes for the Toda Chain*, J. Comput. Phys. **45**, 157 (1982).

[67] G. Benettin, H. Christodoulidi, and A. Ponno, *The Fermi-Pasta-Ulam Problem and Its Underlying Integrable Dynamics*, J. Stat. Phys. **152**, 195 (2013).

[68] H. Christodoulidi and C. Efthymiopoulos, *Stages of Dynamics in the Fermi-Pasta-Ulam System as Probed by the first Toda Integral*, Math. Biosci. Eng. **1**, 359 (2019).

[69] T. Grava, A. Maspero, G. Mazzuca, and A. Ponno, *Adiabatic Invariants for the FPUT and Toda Chain in the Thermodynamic Limit*, Commun. Math. Phys. **380**, 811 (2020).

[70] J. Ford, *Equipartition of Energy for Nonlinear Systems*, J. Math. Phys. **2**, 387 (1961).

[71] M. Onorato, L. Vozella, D. Proment, and Y. V. Lvov, *Route to Thermalization in the α-Fermi-Pasta-Ulam System*, Proc. Natl. Acad. Sci. U.S.A. **112**, 4208 (2015).

[72] L. Pistone, S. Chibbaro, M. D. Bustamante, Y. V. Lvov, and M. Onorato, *Universal Route to Thermalization in Weakly-Nonlinear One-Dimensional Chains*, Math. Biosci. Eng. **1**, 672 (2019).

[73] M. Malishava and S. Flach, *Lyapunov Spectrum Scaling for Classical Many-Body Dynamics Close to Integrability*, Phys. Rev. Lett. **128**, 134102 (2022).

[74] G. Contopoulos, *A Third Integral of Motion in a Galaxy*, Z. Astrophys. **49**, 273 (1960).

[75] M. Kruskal, *Asymptotic Theory of Hamiltonian and Other Systems with All Solutions Nearly Periodic*, J. Math. Phys. (N.Y.) **3**, 806 (1962).

[76] In this regard, it should be noted that a relation between QIs and Lyapunov exponents has already been highlighted in simple systems [43,83].

[77] E. J. Candès, X. Li, Y. Ma, and J. Wright, *Robust Principal Component Analysis?*, J. Assoc. Comput. Mach. **58**, 1 (2011).

[78] P. P. Markopoulos, G. N. Karystinos, and D. A. Pados, *Optimal Algorithms for L₁-Subspace Signal Processing*, IEEE Trans. Signal Process. **62**, 5046 (2014).

[79] J. A. Lee and M. Verleysen, *Nonlinear Dimensionality Reduction*, Vol. 1 (Springer, New York, 2007).
[80] V. A. Brumberg and J. Chapront, *Construction of a General Planetary Theory of the First Order*, Celest. Mech. **8**, 335 (1973).

[81] P. Bretagnon, *Long-Period Terms in the Solar System*, Astron. Astrophys. **30**, 141 (1974).

[82] J. Laskar, *Accurate Methods in General Planetary Theory*, Astron. Astrophys. **144**, 133 (1985).

[83] G. Contopoulos, L. Galgani, and A. Giorgilli, *On the Number of Isolating Integrals in Hamiltonian Systems*, Phys. Rev. A **18**, 1183 (1978).