Clément Ranc

Léo Bernus

Léo Goossens

Abdelkarim Boulfoul

Frédéric Dauvergne

Abed Bounemoura

Alain Albouy

Hervé Manche

Nathan Hara

Antoine Petite

Pierre Auclair-Desrotour

Jeremy Couturier

Mohamad Farhat

Ariane Courtot

Julie Vermersch

Daniel Baguet

Daniel Villegas Pinto

Matthias Sinesael

Margriet Lantink

Chaotic diffusion of the fundamental frequencies of the Solar System

for the relaxing soirées that we had. Many thanks to my

Chapter 1 Introduction

Historical background A hope of everlasting stability

The orbital motion of the Solar System has always been an important subject of inquiry in Astronomy. But it was not until Newton's discovery of the gravitational law in 1687 that mathematical descriptions of orbital dynamics could be formulated. Newton with his law proved that the orbit of a lone planet around its star is an ellipse. This is an elegant demonstration of the laws of planetary motion proposed earlier by Kepler. But the Solar System has more than one planet, each of which can even have satellites, like the Moon of the Earth. Even with one extra body, solving the dynamical equations becomes a much harder problem, since the gradual accumulation of the perturbation from the extra body can disturb and even destabilize the system over a long timescale. Newton himself was greatly troubled by the 3-body problem of Earth-Moon-Sun that he considered it to be his greatest failure. The long term stability of the Solar System seemed to be impossible that Newton thought a divine intervention is needed to restore the system once it goes off-course (e.g. [START_REF] Trinh | Chaos and Harmony: perspectives on scientific revolutions of the twentieth century[END_REF].

The first advancement to this n-body problem was laid out during the late 18 th century by the mutual effort of [START_REF] Lagrange | Recherches sur les équations séculaires des mouvemens des noeuds, et des inclinaisons des orbites des planètes[END_REF] and Laplace (1776) (see [START_REF] Laskar | Is the solar system stable?[END_REF] for historical details). In these works, the two astronomers simplified the problem with a secular approach. First, they averaged the fast orbital timescale and only considered the long-term variation of the orbits. Then, they retained the dominant interaction terms by treating the Solar System perturbatively due to its two properties: the masses of the planets are much smaller compared to the Sun, and the planetary orbits are nearly circles residing close to a common plane. With these simplifications, the long-term evolution of planetary orbits are determined by a system of linear differential equations with constant coefficients. Accuracy of the orbit prediction could be increased by including additional perturbing terms, which induce corrections to the solutions. A notable example was the great 5:2 inequality of Jupiter and Saturn orbital periods (see Wilson, 1985, for historical details). The prediction was consistent with the astronomical observations at the time.

In fact, every apparent deviation in planetary observation seemed to be resolved by a more accurate theory based on Newtonian gravitation. The success of his theory propelled [START_REF] Laplace | Essai Philosophique sur les Probabilites[END_REF] to articulate his optimistic belief of determinism, that is, any prediction of future is guaranteed with sufficient information. For the Solar System, the theory reinforced the belief of its regularity, and by extension, its perpetual stability.

Shift of paradigm

The downfall of the deterministic long-term prediction of the regular planetary orbits was foreshadowed by [START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF], who initiated the shift of paradigm to the modern celestial mechanics that we know today (see [START_REF] Laskar | Is the solar system stable?[END_REF]. Poincaré developed the concept and applications of phase space, where all the states of a dynamical system reside (see [START_REF] Nolte | The tangled tale of phase space[END_REF]; he also re-examined the problem via Hamiltonian formalism. With these tools, he showed that the perturbation series that astronomers used to calculate the planetary motions beyond Lagrange-Laplace theory is generally not convergent. It is not possible to keep improving accuracy of the prediction forever by continual addition of extra terms. More importantly, he introduced the notion of what is called now chaos by pointing out that solutions of many deterministic systems display sensitive dependence to the initial conditions. This revolutionary idea directly refuted the applicability of Laplace's determinism, which permeated the intellectual atmosphere at the time.

Chaos began picking up steam after the birth of electronic computers. It was first officially observed by [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] when he studied the impact of rounding error in numerical integrations of his model for the Rayleigh-Benard convection. The popular butterfly effect is later coined by [START_REF] Lorenz | Predictability: does the flap of a butterfly's wing in Brazil set off a tornado in Texas?[END_REF]'s talk "Does the flap of a butterfly's wings in Brazil set off a tornado in Texas?", to conceptualize the idea that small changes can result in large differences in some situation.

Chaotic dynamics subsequently emerged in many diverse fields: Astrophysics [START_REF] Hénon | The applicability of the third integral of motion: some numerical experiments[END_REF], hydrodynamics [START_REF] Ruelle | On the nature of turbulence[END_REF], biology, of which cardiovascular system [START_REF] Guevara | Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells[END_REF], brain [START_REF] Skarda | Chaos and the new science of the brain[END_REF] and respiratory system [START_REF] Winkler | A chaotic systems analysis of the nasal cycle[END_REF] were modelled, ecology [START_REF] Smale | On the differential equations of species in competition[END_REF], sociology [START_REF] Baker | Chaos, order, and sociological theoy[END_REF], chemistry [START_REF] Hudson | Chaos in the Belousov-Zhabotinskii reaction[END_REF], engineering [START_REF] Matsumoto | A chaotic attractor from Chua's circuit[END_REF], finance and economics [START_REF] Grandmont | Nonlinear Difference Equations, Bifurcations and Chaos[END_REF], etc. It was apparent that in non-linear modelling, chaos becomes the norm, not the exception. The question now pointed towards the most ancient n-body problem, and arguably the most important non-linear dynamics: the Solar System.

The route to chaotic Solar System

Paradoxically in celestial mechanics, the domain where the first notion of chaos emerged, the search for the perpetual regularity of the planetary system was still carried out for most of the 20 th century. The quest was bolstered by the formulation of the [START_REF] Kolmogorov | On conservation of conditionally periodic motions for a small change in Hamilton's function[END_REF][START_REF] Arnol'd | Small Denominators and Problems of Stability of Motion in Classical and Celestial Mechanics[END_REF][START_REF] Möser | On invariant curves of area-preserving mappings of an annulus[END_REF] (KAM) theorem. The theorem states that if the non-linear perturbation of a Hamiltonian dynamical system is weak enough, then quasiperiodic solutions can still exist. This theory can be applied to the planetary motion, albeit with a big caveat first remarked by Henon (see [START_REF] Laskar | Michel Hénon and the stability of the Solar System[END_REF], for historical details), that the planetary masses required for the theory to work need to be smaller than the 1.1. Historical background mass of an electron. Notwithstanding, the KAM theory was still considered as an informal "proof" of the Solar System stability.

A direct way to tell if a complex dynamics, like the Solar System, is chaotic or not, is via numerical integration. Yet, due to the limit of computational power at the time, it was not easy to integrate an n-body problem from Newtonian equations with all considerations of rounding error, relativistic corrections, etc., and it was certainly much more difficult to extend the integration for millions of years.

The limitation was partially due to the short integrating timestep, which has to be much smaller than the orbital period of the fastest planet, Mercury, of roughly 80 days. If the timestep is taken to be 0.3 day for example, integrating the system for 1 million years (Myr) will require a staggering 10 9 iterations of calculation. Moreover, integration of a system is serial computing by nature, and is thus difficult to parallelize. Physicist in early days had to either make strong assumption of the system and often excluded important part of the system itself (e.g. [START_REF] Carpino | Long-term numerical integrations and synthetic theories for the motion of the outer planets[END_REF], only considered the outer planets), or be limited by a short interval of integration (e.g. Quinn et al., 1991, only integrated for 3 Myr). With these difficulties presented, the first person to answer the question "Is the Solar System chaotic?", needed to have good computation power and a good approach.

Jacques Laskar, a researcher from Observatory of Paris, had both. With the assistance of the quickly-developing modern computational power on one hand and the historic analytical methods of averaging in the tradition of Lagrange, Laplace and Le Verrier on the other, he was able to finish the work of Poincare. By numerically integrating the secular model of the Solar System over 200 Myr, [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF] proved that the Solar System is chaotic with the Lyapunov time of 5 Myr1 . This was possible due to a very large timestep of 500 years because the fast orbital motions were averaged. The system of averaged equation, which was developed based on previous works (e.g. [START_REF] Brumberg | Construction of a general planetary theory of the first order[END_REF] with dedicated computer algebra programs, consists of 153 825 terms and is accurate up the second order in planetary masses and 5 th order in eccentricity and inclination.

Only three years later, the MIT researchers, [START_REF] Sussman | Chaotic evolution of the solar system[END_REF] presented the first long direct integration of the whole Solar System over 100 Myr. This was in itself a remarkable achievement and only feasible due to a new efficient symplectic integration technique and a new more powerful computer. This direct integration confirmed that the Solar System is chaotic with roughly the same Lyapunov time of 5 Myr. A tiny uncertainty in the initial conditions, which is bound to happen due to imperfect ephemerides and modelling, will blow up exponentially by increasing 10-fold after approximately every 10 Myr.

The route to unstable Solar System

The hope of proving a regular Solar System was finally shattered, yet one could still have faith that Solar System is stable, at least during its lifetime of 5 billion years (Gyr). This hope was short-lived as [START_REF] Laskar | Large-scale chaos in the solar system[END_REF] showed the possibility of an unstable trajectory where Mercury could either collide with Venus or escape from the system via a rare event sampling scheme. This result showed that there is no dynamical constraint, at least for the terrestrial planets, to prevent the eventual destabilization. A catastrophic event, which is either planetary collision or planetary ejection or both, could potentially occur in the next 5 Gyr. This existence of instability invited even more questions about its quantification, characterization and description.

It should be noted that there are two notions of stability in celestial mechanics. The first one characterizes the regularity of a trajectory, or more precisely, the lack of exponential divergence between two close solutions. The solution of a stable dynamics in this sense can be analytically computed from a convergent asymptotic series. The second notion concerns with the stability of the architecture of the dynamical system. Stability in this sense is defined as the lack of instability, which arises when a catastrophic event, such as a planetary collision, occurs. For the sake of clarity, for the rest of the manuscript we use regularity for the notion of stability in the first sense, and retain the original meaning for the second.

Modern understanding of the Solar System

We came into the 21 st century with a new belief of the chaotic Solar System. It is apparent that trying to predict the Solar System for its remaining future (5 Gyr), which is 1000 times larger than the its Lyapunov time (5 Myr), is a hopeless venture. Yet, the new century also equips us with new updated technology to face new challenges, we have better computers for more intensive integrations, and new space technology for more accurate planetary measurements. In the spirit of keep pushing forwards to obtain the most precise solution, long-term integrations of realistic models of the Solar System were realized [START_REF] Ito | Long-term integrations and stability of planetary orbits in our Solar system[END_REF][START_REF] Varadi | Successive refinements in long-term integrations of planetary orbits[END_REF][START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]Laskar et al., , 2011a;;[START_REF] Zeebe | Numerical solutions for the orbital motion of the Solar System over the past 100 Myr: limits and new results[END_REF] with new planetary ephemerides (e.g. [START_REF] Fienga | INPOP08, a 4-D planetary ephemeris: from asteroid and time-scale computations to ESA Mars Express and Venus Express contributions[END_REF]. The limit of validity for a long-term solution is pushed towards about 60 Myr. Beyond this threshold, a deterministic solution of the Solar System is no longer reliable.

The limit imposed by chaos raises interesting and practical problems concerning the two time directions: future and past. Before addressing these two problems in section 1.3 and section 1.4 respectively, I shall introduce in this section important and relevant characteristics of the Solar System, as well as recent advancements of the field during these recent years.

Overview

The fundamental frequencies Figure 1.1 shows the eccentricity of the Earth over 10 million years ago in the past (Ma, Laskar et al., 2004). Over this period, the evolution of eccentricity, and more generally all orbital components, can be well approximated by a quasi-periodic series. In fact, the underlying frequencies of the series can be considered as combinations of the fundamental secular frequencies of the Solar System.

Each planet is associated with two characteristic frequencies, g j and s j , where j denotes the order of the planet in terms of distance from closest to furthest to the Sun. The frequencies (g j ) j=1,8 portray the precession of perihelion and evolution of orbital eccentricity; the frequencies (s j ) j=1,8 characterize the precession of the orbital plan in precise long-term prediction. Therefore one needs to adopt a probabilistic approach to probe beyond the chaotic limitation. Instead of integrating one solution from one initial condition, one needs to propagate in time a probability density function (PDF), which should mimic the uncertainty in the initial conditions as well as in the model.

Formally, evolution of a density function in phase space of a dynamical system could be captured exactly by the Liouville equation [START_REF] Liouville | Note sur la Théorie de la Variation des constantes arbitraires[END_REF]. For a low-dimensional chaotic map, the Liouville equation can be solved directly to generate a map for the density function (e.g. [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF] for the Chirikov standard map, [START_REF] Ford | The Arnol'd cat: Failure of the correspondence principle[END_REF] for Arnold map ). One can also derive from the Liouville equation (e.g. [START_REF] Landau | On the theory of phase transitions[END_REF]; see [START_REF] Escande | When can the Fokker-Planck equation describe anomalous or chaotic transport[END_REF], for further discussion) a Fokker-Plank equation to describe the transport of the PDF. The expression of the diffusion coefficient of the transport can be accordingly obtained (e.g. [START_REF] Rechester | Calculation of turbulent diffusion for the Chirikov-Taylor model[END_REF] for standard map; Lieberman and Godyak, 1998 for Fermi map). However, for a complex system involving many degrees of freedom like the Solar System, both integrating numerically Liouville equation or finding a Fokker-Plank approximation are extremely difficult problems.

Therefore, in practice, a Monte Carlo method is often used as follows. An ensemble of initial conditions, which are randomly sampled according to an initial PDF, is numerically integrated to obtain a sample of time-evolving solutions. The solutions altogether represent approximately the time-evolving PDF. The solutions are equiprobable, each depicts a possible future (or past) pathway of the system. They are instrumental to the statistics of the system, from which we can draw probabilistic long-term prediction.

One can also see each solution like a mobile scanner probing the phase space. Imagine dropping a droplet of ink in a cup of water. Each of the ink particle will eventually part away from the clump in its turbulent motion to discover the surrounding accessible environment. After some time, the concentrated ink droplet forgets its initial state and dilute itself in the water. Imagine again that we cannot see neither the water nor the cup holding it, but we want to measure them. Just by following the ink, we can have a good picture of the cup (what it looks like, its dimension, etc.). In the same manner, we can understand the constraints in the phase space imposed by the dynamics just by looking at the ensemble of its chaotic solutions. [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF] performed a statistical analysis over 1001 different integrations of secular equations over 5 Gyr starting from closely-spaced initial conditions. The PDF of the eccentricity and inclination of the Solar System planets were estimated from the sample (Fig. 1.5). From the evolution of the PDFs through time, the distinction between the ISS and OSS is once again clear. The PDFs of the OSS planets remain unchanged for the whole 5 Gyr, while those of the ISS slowly diffuse overtime. Moreover, it was shown that the outer planets' PDFs are in fact the PDFs of quasi-periodic motions. For the ISS, the PDFs can be reasonably fitted by Rice distributions with parameters characterizing the diffusion. The eccentricity can range up to 0.45 for Mercury in 5 Gyr and 0.15 for Mars, while Venus and Earth's variations are smaller with values being below 0.08.

Density estimation of orbital elements

It is natural and tempting at the same time to make sense of the descriptive statistics of a chaotic system from a perspective of statistical mechanics. For example, to explore the distribution of extrasolar planet orbits, [START_REF] Tremaine | The statistical mechanics of planet orbits[END_REF] makes use of certain stability conditions with the shear-sheet approximation of a planetary system, and the equiprobability hypothesis, namely planetary configurations are uniformly distributed conforming to the stability constraints. What makes such a statistical model attractive is that, from a few simple hypotheses, one can have a precise prediction for the distribution of any pertinent orbital parameters, such as eccentricity and semi-major axis relations.

Inspired by [START_REF] Tremaine | The statistical mechanics of planet orbits[END_REF]'s work, [START_REF] Mogavero | Addressing the statistical mechanics of planet orbits in the solar system[END_REF] showed that the [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF]'s distribution of ISS planets can be understood via the constraint of its total angular momentum deficit (AMD) [START_REF] Laskar | Large scale chaos and the spacing of the inner planets[END_REF][START_REF] Laskar | On the spacing of planetary systems[END_REF]. AMD is the difference of the angular momentum between a circular planar orbit and its elliptic inclined counterpart. AMD measures the non-linearity of the secular system, and is thus used in the stability analysis in the exoplanetary system [START_REF] Laskar | AMD-stability and the classification of planetary systems[END_REF]. For a secular isolated system, its total AMD is conserved. [START_REF] Mogavero | Addressing the statistical mechanics of planet orbits in the solar system[END_REF]'s model assumes the conservation of the ISS's AMD and the hypothesis of equiprobability for all states preserving this constraint in phase space. The correspondence between the PDFs of eccentricity and inclination of this statistical theory and numerical simulation is excellent for Earth and Venus, acceptable for Mars but poor for Mercury. Indeed, the hypothesis is problematic since Mercury seems to still remember its initial state despite 5 Gyr of chaotic evolution. In addition, the AMD of inner planets is not conserved, because the ISS is not isolated from the OSS due to the coupling secular resonances (e.g. g 1 -g 5 -s 1 + s 2 ) between the two systems. Nevertheless, the result is promising given the simple hypothesis of the model. For a better comparison, once should look for better dynamical constraints, and prolong the integration so that the system has sufficient relaxation time to forget its initial state.

Destabilization of the ISS 1.3.1 Probability of instability

The next logical question after the confirmation of a possibly unstable future of the ISS [START_REF] Laskar | Large-scale chaos in the solar system[END_REF] is: What is the probability of destabilization for its remaining life time of roughly 5 Gyr? From the 1001 solutions of a secular model of the Solar System at degree 6 in eccentricity and inclination, [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF] showed that this probability is about 1%. Nevertheless, one can be skeptical about the veracity of this prediction. Since the model is equivalent to a secular Hamiltonian developed in eccentricity and inclination as small parameters, in the event of orbital excitation to high eccentricity, the expansion of the secular Hamiltonian is generally divergent [START_REF] Wintner | The analytical foundations of celestial mechanics[END_REF]. Furthermore, such an event can give rise to the crossing of two planetary orbits, during which emerge important dynamical factors such as mean motion resonances and variation of planetary semi-major axis; these factors cannot be captured by the averaged model.

A year later, with a huge computational power in disposal, [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF] were able to obtain 2501 direct numerical simulations of the Solar System with an N -body code. This more complete and non-averaged model confirmed the result of the secular model, the probability of a catastrophe in the ISS is roughly 1%. This estimation is subsequently reproduced in later studies [START_REF] Zeebe | Highly stable evolution of Earth's future orbit despite chaotic behavior of the Solar System[END_REF][START_REF] Abbot | Rare event sampling improves mercury instability statistics[END_REF]. The lower panel in the Fig. 1.6 shows the maximum eccentricity of Mercury of 2501 integrations of [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF]. Among these, there were 20 solutions with eccentricity of Mercury higher than 0.9; most of them became unstable in the late stage of the integration, which displays the evolution of the precession frequency of Mercury g 1 for collisional secular solutions of the ISS.

When the two frequencies are comparable, that is g 1 -g 5 ≈ 0, the resonance between Mercury and Jupiter is activated as the solution enters and explores its chaotic zone. The system becomes much more chaotic during this activation with much larger Lyapunov exponent. As a result, the typical dynamical timescale of this stage, which only last for a few Myr, is short [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF]. This is indicated in the rapid bounce back once g 1 reaches g 5 in the unstable solution (Fig. 1.7). Despite its short duration, the resonance is strong and very efficient. The activation of the resonance creates a pathway to significantly accelerate the transfer of energy and AMD between the smallest planet and the biggest planet in the Solar System, essentially setting up a highway channel between the ISS and OSS. As a consequence, the Mercury eccentricity is pumped to a high value, which is often considered as the first sign of instability. Moreover, since the ISS is tightly connected via a web of secular resonances [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF], the excitation from Mercury can be readily transferred to other planets, thereby potentially destabilizing the whole system [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF].

To model the excitation of Mercury eccentricity, [START_REF] Boué | A simple model of the chaotic eccentricity of Mercury[END_REF] construct an integrable 1 degree-of-freedom Hamiltonian of a massless Mercury in vicinity of the resonance g 1 -g 5 . The Hamiltonian is expanded in ratio of semi-major axis between two planets instead of eccentricity as usual, so that the model is valid for the full range of eccentricity value. With this simple dynamics, a rapid excitation of eccentricity over several Myrs, which is observed in unstable numerical solutions, is captured.

Contribution of General Relativity

Historically, the theory of general relativity (GR) plays an important role in the dynamics of the Solar System. The prediction from the Newtonian equations for the perihelion precession of Mercury deviates from the astronomical observations. This problematic anomaly, which was first remarked by Le [START_REF] Verrier | Lettre de M. Le Verrier à M. Faye sur la théorie de Mercure et sur le mouvement du périhélie de cette planète[END_REF], prompted the astronomer to formulate an even more problematic hypothesis of a missing planet. The anomaly is actually out of scope of the classical Newtonian framework. It was not until the establishment of GR, that [START_REF] Einstein | Die Grundlage der allgemeinen Relativitätstheorie[END_REF] were able to explain the anomalous precession of Mercury. In the framework of classical equations, the effect of GR is often included by a post-Newtonian expansion in 1/c 2 as the small parameter (see Will, 2018).

The effect is most prominent for Mercury because it is the closest planet to the Sun. The perihelion frequency of Mercury g 1 = 5.59 arcsec yr -1 , out of which roughly 0.4 arcsec yr -1 can be attributed to relativistic effects. For other planets, the impacts are smaller but can still be detected in recent observations [START_REF] Anderson | Recent developments in solar-system tests of general relativity[END_REF]. The change in precession of Venus is 0.086 arcsec yr -1 , 0.038 arcsec yr -1 for Earth, and 0.01 arcsec yr -1 for Mars (e.g. Laskar, 1999, table 4) For outer planets, the effect is even smaller (e.g. 0.0006 arcsec yr -1 for Jupiter).

This shift in the frequency g 1 of Mercury is substantial if we compare it to its chaotic variation, which is 0.13 arcsec yr -1 after 100 Myr (table 1.1). Dynamically, this GR shift is comparable to the half-width of the strongest secular resonance [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. If we consider the instability of the ISS as the chaotic diffusion of g 1 towards g 5 , the relativistic shift widen the gap between two frequencies by almost 50% from 0.9 arcsec yr -1 (without GR) to 1.33 arcsec yr -1 (with GR). It is therefore not unexpected that a purely Newtonian model of the Solar System is more unstable than the one with GR [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Batygin | On the Dynamical Stability of the Solar System[END_REF][START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF] What is surprising is the degree of change brought about by the relativistic effect. It utterly transforms the nature of the Solar System evolution in its lifetime. The contrasting difference is apparent in the comparison between two scenarios shown in Fig. 1.6. Without GR, the ISS is no longer marginally stable with small chance of being destabilized. In fact an unstable future, with the probability of 60% in the next 5 Gyr, is expected to be a normal occurrence. This is also true for the nearer future. In the next 500 Myr, there is more than 1-out-of-10 chance that the system is destabilized [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF].

The BMH model

To our knowledge, the only attempt to model the instability of the ISS in its full scope, that is to give prediction of instability time within a chaotic model, is the work of (Batygin, Morbidelli and Holman, 2015, BMH). Following the work of [START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF][START_REF] Boué | A simple model of the chaotic eccentricity of Mercury[END_REF], [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF] constructs a model of massless Mercury forced by several eigenmodes, which are composed of secular frequencies of Jupiter (g 5 ) and Venus (g 2 , s 2 ). The BMH model has 2.5 degrees of freedom. Its Hamiltonian is expanded up to degree 4 in eccentricity and inclination of Mercury and only retains several supposedly relevant harmonics while discarding the rest. In short, the Hamiltonian, which is made up of about a dozen of monomials and four harmonics, is a small subset of the Hamiltonian of the forced ISS at degree 4 [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF], which has 6304 monomials and 2748 harmonics.

Despite its simplicity, the BMH model possesses two important aspects of Mercury: chaos and instability. Moreover, the model also captures the qualitative disparity between the chaotic timescale and the instability timescale: its expected instability time of roughly 1 Gyr is three orders of magnitude larger than the Lyapunov time of about 1 Myr. Nevertheless, the simple model clearly falls short in terms of quantitative aspects: (1) It is too chaotic, with Lyapunov time 5 times smaller than the expected value, and (2) the model is too unstable, it predicts Mercury will be almost certainly unstable after 5 Gyr.

Nevertheless, the physics of this simple model is interesting. Its Hamiltonian, which varies slowly, defines in the phase space a bounded region near the initial condition, and an unbounded region further away. Because of the bounded constraint, the dynamics is bifurcated into two different timescales: the fast Lyapunov time, with which solutions wander inside the bounded topological subspace determined by the Hamiltonian; the slow diffusion timescale of the Hamiltonian itself, with which the bounded subspaces move and bring the solutions with it. The dynamics is stable as long as its Hamiltonian bounds the motion, however, as the Hamiltonian slowly diffuses towards certain values where it can no longer provide the same restriction, the solution escapes into the open phase space and becomes unstable.

Based on this model, [START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF] show that the instability statistics of the BMH model can be reproduced by a simple stochastic model. In this model, which the Hamiltonian randomly diffuses according to the 1-dimensional standard Brownian motion with reflective barrier in one side and the critical value of the topological transition leading to instability in the other. With this statistical model, the analytical probability function of the hitting time can be derived exactly. The rare events will concentrate around the instanton path, which is just a straight line from the initial point to the critical value. In fact, since the BMH model is too unstable and will destabilize after 5 Gyr with high probability, so the instanton, defined from the rare event characterization, is no longer valid at 5 Gyr. Therefore, [START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF] has to reduce the time span to 445 Myr in order to study the instanton.

Problem of instability

Although several descriptive aspects of the destabilization have been revealed, as its triggering resonance is known and its probability was quantified, a global picture still remains obscure. The problem appears even more puzzling upon close examination. The ISS, which is forced by the outer planets, is an open system. It has no apparent integrals of motions to constrain its chaotic motion, thus a sizeable unstable region in phase space is accessible for the dynamics. Moreover, the system has plenty of time to discover the phase space including unstable regions, since the remaining lifespan of the Sun is roughly 5 Gyr, that is 3 orders of magnitudes longer than the Lyapunov time of the dynamics. All these aspects beg the question of why the ISS is so resistant to calamity? why is there only 1% chance of destabilization over 5 Gyr? given that the dynamics of the ISS is in a vigorous and unbounded chaotic motion for so long. Imagine a drunk man walking on a sidewalk next to a busy highway without any visible barrier all night, it would be puzzling if there are no accidents in the morning.

This paradox has been recognized for a long time (e.g. [START_REF] Lecar | Chaos in the solar system[END_REF], but it was not until recently that it was articulated and properly addressed by the BMH model [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF]. The quantitative predictions of the model are inconsistent with our actual knowledge of ISS. The model is more chaotic and more unstable, as its expected destabilization time (∼ 1 Gyr) is at least one order of magnitude smaller than a more realistic model (∼ 40 Gyr, c.f. [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. Nevertheless, the BMH dynamics succeeds in producing two well-separated timescales, Lyapunov time and destabilization time. In this slow-fast dynamics, the destabilization corresponds to the topological transition from bounded to unbounded constraint determined by its Hamiltonian.

The attempt of [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF] to model the instability of the ISS raises two important questions concerning its generalization. First, is it possible to build a simple model that reproduces correctly the statistics of the stability? Indeed, BMH model starts from the linear Lagrange-Laplace theory and then "sequentially enhance the complexity until the desired behaviour is adequately represented". One can thus practically enhance the complexity by keep adding more terms to the Hamiltonian until the model is complex enough and the obtained statistics is acceptable. To this end, since there are an infinite number of terms one can choose, a better question should be what are the important terms, which play essential roles in the stability of the ISS, and cannot be neglected, that we need to add to the Hamiltonian? Answering this question would reveal the intricate dynamical interactions, which are represented by corresponding harmonics, between planets and their contributions to the system stability.

The second question rises from the slow-fast nature of the BMH model. Does the ISS have a similar slow-fast dynamics? If yes, how can we tell? It should be remarked that the approach of BMH cannot be simply generalized to the whole ISS. Since its bounding constraint emerges partly from the low dimensionality of the model. If more planets are added and dimensionality is accordingly increased, the Hamiltonian no longer defines a bounded region in phase space, there will be no topological transition as the solution can a priori escape from the get go. Moreover, apart from the Hamiltonian, there are no other visible constraints.

In the chapter 4 and the chapter 5 of this thesis, we attempt to resolve these two questions.

Geology and Astronomy

The question of destabilization is no longer pertinent for the past evolution of the Solar System, for Mercury exists now, it must have existed before. Indeed, the orbital evolution of the Solar System since the last 5 billions of years ago, by definition, has already occurred. Although there can be many possible futures, there is only one true past of the Solar System, recovery of which is an important problem. Because of the chaos, integration from a realistic model with precise ephemerides only takes us as far as 60 million years ago (Laskar et al., 2011a). Beyond this threshold, we need to look for the traces that the dynamics left on its path, which can be found in the geological records on Earth. Therefore, the recovery of the history of the Solar System becomes a problem of recovery of its imprints from the geological data.

Overview

The connection between these two seemingly disconnected fields, namely astronomy and geology, has been hypothesized and demonstrated for the last 150 years. The connection officially started with the work of [START_REF] Agassiz | Etudes sur les glaciers. Aux frais de l[END_REF], who observed and argued for the cyclic nature of the ice age, and prompted the search for the origin of these glacial cycles. [START_REF] Adhémar | Révolutions de la mer, déluges périodiques[END_REF] and later [START_REF] Croll | Climate and time[END_REF] with the solution of Le Verrier (1856) laid out the basic principal of astronomical forcing as the explanation for the long-term climate swings and the resulted periodic glaciation. The theory with all necessary ingredients was pieced together by [START_REF] Milankovitch | Canon of insolation and the iceage problem[END_REF] (see [START_REF] Imbrie | Astronomical theory of the Pleistocene ice ages: A brief historical review[END_REF], for more historical details).

The central element of the theory that links the long term climate cycles and the astronomical variation of Earth's orbital element is the amount of sun light cast on Earth, termed insolation (incoming solar radiation). This variation of insolation affects atmospheric temperature and circulation, which in turn influence oceanic temperature, oceanic circulation, rainfall pattern and biological productivity, etc. All these factors control the deposition, transport and composition of sedimentation, leaving on its trail physical, biological and chemical markers, and forming stratigraphic succession. These markers -measured on sedimentary records, are thus used as proxy for the climatic variations; notable examples are rock magnetism, sediment compositions, isotopic ratios, and even qualitative features like color and facies variation (see e.g., [START_REF] Strasser | Cyclostratigraphy concepts, definitions, and applications[END_REF][START_REF] Hinnov | Cyclostratigraphy and astrochronology[END_REF]. Many of these stratigraphic successions, such as layered sedimentary rocks on mountain sides (i.e. strata), are visible to the naked eye and can be visually identified with orbital forcings (e.g., [START_REF] Hilgen | Integrated stratigraphy and astronomical calibration of the Serravallian/Tortonian boundary section at Monte Gibliscemi (Sicily, Italy)[END_REF][START_REF] Kuiper | Synchronizing rock clocks of Earth history[END_REF][START_REF] Wu | Time-calibrated Milankovitch cycles for the late Permian[END_REF].

The evidence for Milankovitch theory was established by the landmark study of [START_REF] Hays | Variations in the Earth's orbit: pacemaker of the ice ages[END_REF] entitled "Variations in the Earth's orbit: Pacemaker of the Ice Ages". The study integrated necessary ingredients for an unambiguous demonstration of astronomical forcings on paleoclimate: A long and continuous record (deep-sea core documenting 500 000 years of history), on which robust proxy for climate variation (marine δ 18 O3 as a proxy for ocean temperature), and controlled by independent timing method (radiocarbon dating). This study set up basic rationale for the new field cyclostratigraphy, and inspired geologists to begin the extensive search for astronomical signals in stratigraphic records (see [START_REF] Fischer | Climatic rhythms recorded in strata[END_REF].

Since then, interest in the field surged as many geological records of climate responding to insolation forcing have been identified. Cyclostratigraphy is said to have entered a "golden age" [START_REF] Hinnov | Cyclostratigraphy and astrochronology in 2018[END_REF], with multi-million year-long data ranging from the deep-sea sequence during the recent ice age several millions of years ago [START_REF] Lisiecki | A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records[END_REF] to the deep-time strata in Australia during the Precambrian 2.5 billions of years ago [START_REF] Lantink | Milankovitch cycles in banded iron formations constrain the Earth-Moon system 2.46 billion years ago[END_REF].

The astronomical solutions, when used to calibrate a stratigraphic sequence, provide unprecedented high-resolution temporal constraints to the data. This practice, known as astrochronology, is about tuning the geological record to astronomical forcing (see Meyers, 2019, for definitions), and convert the record from depth domain into time domain. By anchoring the geological archive to an absolute astronomical reference, the standard geological timescales can be robustly constructed (e.g., [START_REF] Gradstein | A new geologic time scale, with special reference to Precambrian and Neogene[END_REF][START_REF] Gradstein | The geologic time scale[END_REF][START_REF] Gradstein | The chronostratigraphic scale[END_REF]. It should be remarked that the value of absolute geological ages goes far beyond archiving Earth's geological history. For instance, in order to understand the impact of the carbon cycle on the historical climate change and thereby estimate its future projection, it is essential to know the precise timing of forcings like volcanism, impacts and orbital forcing, etc., as well as the timing of the climate indicators (e.g., [START_REF] Cramwinckel | Synchronous tropical and polar temperature evolution in the Eocene[END_REF].

Milankovitch cycles

The long-term variations of orientation and orbits of the Earth, which have direct impacts on its long-term climate change, are called Milankovitch cycles. The Milankovitch cycles can be divided into three groups: (1) fluctuation of orbit shape (eccentricity), (2) wobbling of the rotational axis (precession), (3) oscillation of the angle of the axial tilt with respect to the orbital plane (obliquity). It should be noticed that although the orbital variations are often referred to as "cycle", they are much closer to quasi-periodic motion and should thus be regarded as such. The Milankovitch cycles are obtained via astronomical solutions of the orbital dynamics of the Solar System, and the rotational dynamics of the Earth-Moon system. They are the primary inputs to compute the insolation that Earth receives on the top of its atmosphere. The insolation equation was first introduce by [START_REF] Meech | On the relative intensity of the heat and light of the Sun upon different latitudes of the Earth[END_REF], and refined by [START_REF] Milankovitch | Canon of insolation and the iceage problem[END_REF]; modern calculation of insolation is often referred to [START_REF] Berger | Long-term variations of daily insolation and Quaternary climatic changes[END_REF]Laskar et al., 1993a;[START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF], see Hinnov, 2018, for a recent review). In the followings, I will discuss about how the Milankovitch cycles qualitatively affect the insolation, and their quantitative expression in insolation variations, which were imprinted and can be recovered in stratigraphic sequences.

Precession and Obliquity cycle

The rotational axis of Earth is currently tilted with respect to the perpendicular axis of its orbital plane by an angle of 23.5 • . This angle, which is called obliquity or axial tilt, presently oscillates from 22.1 • to 24.6 • (Laskar et al., 1993b). In addition to being tilted, the rotational axis of the Earth also precesses, similar to the precession of a spining top.

Effects on insolation

Both precession and obliquity cycles have direct impact on the latitudinal and seasonal distribution of insolation that Earth receives. Here I briefly demonstrate this influence, more details are referred to [START_REF] Ruddiman | Earth's Climate: past and future[END_REF]. Obliquity breaks the symmetry of the incidence angle of insolation with respect to the postion on Earth's orbit. As one hemisphere always points towards the Sun half of the year and spend the other half tilt away from the Sun, obliquity is the reason why the Earth has seasons. Therefore, variation of obliquity can make seasons either more extreme (increasing obliquity) or more moderate (decreasing obliquity). Due to geometrical effect, obliquity's cycles have greater influence on regions of higher latitude, and in the winter season.

To understand the impact of the precession on insolation of each hemisphere, one needs to put the inclined rotational axis on an eccentric orbit. The orientation of the rotational axis determines the position of winter (when it directs away from the Sun) and summer (when it directs towards the Sun) of a hemisphere on the orbit. Take the northern hemisphere for example, winter currently occurs near the perihelion, where the Earth is closest to the Sun, and summer is around the aphelion, where the Earth is furthest from the Sun. This moderates the contrast between the two seasons, as Northern hemisphere in winter is relatively warmer due to closer distance to the Sun and inversely, it is cooler in summer. On the contrary, the current orientation of the rotational axis widens seasonal difference for the Southern hemisphere.

This effect is only possible on an eccentric orbit, where the planet can approach or recede from the star along its orbit. Therefore, the intensity of insolation variation is controlled by how eccentric the orbit is, hence the eccentricity cycle modulate the precession cycle via its amplitude. As the rotational axis precesses, the positions of winter and summer of a hemisphere also precess along the eccentric orbit. Moreover, the eccentric orbit of Earth also rotates on its orbital plane, this motion is called apsidal precession of perihelion. The combined effect of this apsidal precession and the precession of the rotational axis (precession for short) is called climatic precession.

Precession frequency. On a short timescale, the rotational axis of the Earth can be considered to precess with a constant rate p, known as precession frequency. The period of one revolution with the current value for p is roughly 25.7 kyr [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]. Combining with the apsidal precession of perihelion of the Earth, marked by the 5 inner fundamental frequencies (g i ) i=1,5 , the climatic precession can be analytically approximated in the short term by quasi-periodic series of p + g i ; the period of climatic precession on average is around 23 kyr. The amplitude of this climatic precession is determined by the eccentricity cycles. Similarly, the evolution of Earth's obliquity in the near future can be characterized by the combination of the secular inclination frequencies s i and the precession frequency p; the periods of the three dominant terms p + s 3 , p + s 4 , p + s 6 are is small, and scales as e 2 . On the other hand, eccentricity is very important for seasonal insolation variation. The difference in distance between the closest point (perihelion) and the furthest point (aphelion) of Earth's orbit to the Sun scales with 2e. Considering also that the solar radiation is proportional to inverse squared distance, the difference in insolation on Earth between perihelion and aphelion can be accordingly shown to vary as 4e. This difference shown in number is significant: the current value is roughly 6%, and was more than 20% 1 Myr ago as Earth's orbit was more eccentric.

The eccentricity influences the seasonal contrast of insolation: it can either amplify or diminish the difference between winter and summer of a hemisphere depending on the orientation of the rotational axis. Therefore, the major role of eccentricity cycle is to modulate the precession cycle, and control the amplitude of expression for the precession in insolation.

The grand cycles. The evolution of the eccentricity of the Earth in the past 10 Myr is shown in the Fig. 1.1. The value can go up and down from 0 to 0.06 with many overlapping quasi-periodic cycles. Currently, we are at a trough of these cycles, with eccentricity near that lower range at 0.0167. The Fig. 1.8 shows the periodogram of the last 10 Myr of the Earth's eccentricity. Most if not all of the peaks can be explained as the combinations of the five secular fundamental frequencies of the Solar System (g i ) i=1,5 . They are the eccentricity cycles. The most prominent peak with period of 405 kyr corresponds to the difference between Venus and Jupiter frequencies g 2 -g 5 . On the right side of this peak, there are four main shorter eccentricity cycles with periods around 100 kyr: g 3 -g 2 , g 4 -g 2 , g 3 -g 5 and g 4 -g 5 . On the left side are the longer cycles: g 2 -g 1 , g 1 -g 5 and g 4 -g 3 . The longest eccentricity cycle g 4 -g 3 has period of around 2.4 Myr. Since the eccentricity cycles are relatively long, and they determine the amplitude of climate precession cycles, they are often called "grand cycles" or "modulating cycles" [START_REF] Hinnov | Cyclostratigraphy and its revolutionizing applications in the earth and planetary sciences[END_REF].

The metronomes. As time progresses, the fundamental frequencies change along the dynamical solution, and these cycle change too. Due to the chaotic motion of the Solar System, the motion can only be traced back as far as about 60 Myr. Nevertheless, some of the astronomical forcing are stable and prominent enough to be used to calibrate geological records beyond the chaotic threshold (see [START_REF] Laskar | Astrochronology[END_REF], for a review). These stable cycles are called metronomes, and they are the backbone of modern cyclostratigraphy. The most commonly-used metronome is the 405-kyr eccentricity cycle g 2 -g 5 . The cycle is stable because both g 2 and g 5 do not vary significantly (see table 1.1). The g 2 -g 5 metronome is also the strongest eccentricity cycle (Fig. 1.8), and is thus often visible in astronomically-forced stratigraphic successions. Recently, the inclination cycle s 3 -s 6 with period of 174 kyr has also been suggested as a metronome for the obliquity-dominated stratigraphic sequences [START_REF] Boulila | Towards a robust and consistent middle Eocene astronomical timescale[END_REF][START_REF] Charbonnier | Obliquity pacing of the hydrological cycle during the Oceanic Anoxic Event 2[END_REF]. Unlike g 2 -g 5 , the s 3 -s 6 cycle is not the strongest component in the obliquity spectrum, but it is far from other cycles and can be relatively easy identified.

Resonance transition. The transitions between two overlapping resonances, which is a strong evidence for chaos in the Solar System (Laskar, 1990) and is often found in the orbital solution (e.g., Laskar et al., 2011a), can also be detected in the geological records of paleoclimate. The detected transition is between the two secular resonances of perihelion precession of Mars and Earth: θ 2:1 = 2(g 4 -g 3 ) -(s 4 -s 3 ) and θ 1:1 = (g 4 -g 3 ) -(s 4 -s 3 ). As previously discussed in Sect. 1.2.2, the Solar System is currently in the θ 2:1 resonance, which means that this combination of resonance is currently close to 0. At present, the period of the grand eccentricity cycle g 4 -g 3 is around 2.4 Myr, and the period of the grand obliquity cycle s 4 -s 3 = 2(g 4 -g 3 ) is roughly 1.2 Myr. Along the chaotic evolution of the Solar System, if a transition from θ 2:1 to θ 1:1 occurs, the periods of these two grand cycles will be equal. Therefore the resonance transition, if it occurred in the distant past of the Solar System, can be detected by identifying these two cycles in the geological records. [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF] has found the evidence for this transition in the Libsack core in North America during the interval 83-90 Ma. Since then, many other geological evidences in different places from different intervals has been proposed (e.g., [START_REF] Westerhold | Astronomical calibration of the Ypresian timescale: implications for seafloor spreading rates and the chaotic behavior of the solar system[END_REF][START_REF] Zeebe | Solar System chaos and the Paleocene-Eocene boundary age constrained by geology and astronomy[END_REF][START_REF] Wu | Continental geological evidence for Solar System chaotic behavior in the Late Cretaceous[END_REF].

Problems of cyclostratigraphy

Despite its ongoing role in revolutionizing the geological dating and setting up important inputs and constraints to the general framework for geoscience, there is still an underlying skepticism of the reliability of cyclostratigraphy. The challenges of cyclostratigraphy stem from two different fields, namely geology and astronomy, that it attempt to connect. The first problem lies on the quality of the geological data and its ability to record faithfully astronomical signals. The intermediate processes which propagate the astronomicalforcing driven insolation to the geological record are prohibitively complex, with numerous different other factors that can contaminate the signal (see [START_REF] Meyers | Cyclostratigraphy and the problem of astrochronologic testing[END_REF], for more indepth discussion). As a result, the stratigraphic sequences can be highly irregular (due to the variable formation rate, e.g. [START_REF] Olsen | Long-period Milankovitch cycles from the Late Triassic and Early Jurassic of eastern North America and their implications for the calibration of the Early Mesozoic time-scale and the long-term behaviour of the planets[END_REF][START_REF] Haerdle | Asymptitic optimal data-driven bandwidths for regression under dependence[END_REF], riddled with missing gaps (known as hiatus, e.g. [START_REF] Gale | Orbital tuning of Cenomanian marly chalk successions: towards a Milankovitch time-scale for the Late Cretaceous[END_REF], distorted by physio-chemical processes after sedimentary deposition (e.g. diagenesis, [START_REF] Hallam | Origin of minor limestone-shale cycles: Climatically induced or diagenetic?[END_REF], contaminated by different sources other than the Milankovitch cycles (e.g. stochastic climate variability, [START_REF] Hasselmann | Stochastic climate models part I[END_REF]. Moreover, even for the analysis of the same stratigraphic sequence, quantitative outcomes can be different according to different investigators [START_REF] Sinnesael | The Cyclostratigraphy Intercomparison Project (CIP): consistency, merits and pitfalls[END_REF]. The problems of false negatives and false positives in existence of detectable astronomical signal in geological records are also brought into discussion [START_REF] Vaughan | Detecting cycles in stratigraphic data: Spectral analysis in the presence of red noise[END_REF][START_REF] Hilgen | Stratigraphic continuity and fragmentary sedimentation: the success of cyclostratigraphy as part of integrated stratigraphy[END_REF][START_REF] Weedon | Cyclostratigraphy, stratigraphic gaps and the duration of the Hettangian Stage (Jurassic): insights from the Blue Lias Formation of southern BritainGP Weedon et al. Cyclostratigraphy of the Blue Lias Formation[END_REF]. Nevertheless, with the remarkable progress that cyclostratigraphy as a new field has seen in the recent years, many of these problems have been acknowledged and addressed (e.g. [START_REF] Hinnov | Cyclostratigraphy and astrochronology in 2018[END_REF][START_REF] Meyers | Proterozoic Milankovitch cycles and the history of the solar system[END_REF][START_REF] Sinnesael | The Cyclostratigraphy Intercomparison Project (CIP): consistency, merits and pitfalls[END_REF]. Some stratigraphic records with well-preserved orbital signal and analyzed with well-defined techniques can provide consistent and reproducible estimations for the fundamental secular frequencies (e.g. [START_REF] Olsen | Mapping Solar System chaos with the Geological Orrery[END_REF], precession frequency (e.g. [START_REF] Meyers | Cyclostratigraphy and the problem of astrochronologic testing[END_REF]Malinverno, 2018, Lantink et al., 2022), and offer evidence for the chaotic transition (e.g. [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF].

The second challenge of cyclostratigraphy is also the challenge of celestial mechanics of the Solar System. The time limit of a precise orbital solution is about 60 Myr (Laskar et al., 2011a), this limit is even shorter for obliquity and precession solution [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]. Although the exact astronomical reference is infeasible beyond 60 Myr, some stable components of the astronomical forcing such as the 405-kyr g 2 -g 5 cycle can still be applicable. These metronomes are often used in tandem with independent radioisotopic dating to calibrate the stratigraphic sequences beyond the Cenozoic [START_REF] Hinnov | Cyclostratigraphy and its revolutionizing applications in the earth and planetary sciences[END_REF]. However, the chaotic uncertainty of the orbital forcing, even for the metronomes, in the geological data beyond 60 Myr always persists and gets worse with time due to chaotic diffusion. Thus, the second problem of cyclostratigraphy concerning the astronomical signal is how to quantify properly this chaotic uncertainty. This is also the problem that this thesis attempts to address.

In the resolution of these challenges, we can hope to trace back the orbital variations in the geological records and use them to constrain the astronomical solutions beyond the predictability horizon.

Objectives and organization of this thesis

Because of the chaotic nature of the Solar System [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF], the study of its evolution in the long term has to be set in a statistical paradigm. In this paradigm, more is better. The question of obtaining the most accurate parameters and initial conditions from the most complete model is no longer pertinent, one is never enough. We need thousands, and possibly millions, to quantify the long-term statistics of the Solar System (e.g. [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF]. In this context, there are two important non-resolved problems concerning the two time directions of its evolution.

The problem in the future consists of questions regarding the destabilization of the inner planets: its quantification, characterization, description. I will address this problem in the chapter 4 and chapter 5. They are the first two pieces of the forthcoming series of papers aiming to provide a thorough understanding of the destabilization of the Solar System. The chapter 4 (published in Hoang, Mogavero and Laskar 2022) examines the problem with extensive numerical simulations from a hierarchy of comprehensive secular models of the forced ISS (Laskar, 1990;[START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. From these simulations, we (1) present the descriptive statistics of instability from a wide range of models, and ( 2) show what it takes to produce a consistent estimation by pointing out the surprisingly important impacts of the higher order and non-resonant interaction. In chapter 5, we demonstrate the hidden slow-fast nature from the chaotic dynamics of the ISS, in which there exist a wide range of dynamical timescale, implicitly implied in the spectrum of the Lyapunov exponents. We provide concrete evidences by identifying the three quasi-integrals of motion via a semi-analytical and a statistical approach. These quasi-integrals are the dynamical constraints that stabilize and prevent the ISS from destabilizing quickly within its remaining timespan.

In the reverse direction -the past, the recovery of the history of the Solar System is a more important and relevant question. Since a precise orbital solution from integrating a Solar System model is only valid up to 60 Myr (Laskar et al., 2011a), it is required to look for the traces in the geological records to overcome this limit. This problem is, in its core, multi-disciplinary. In the chapter 3 (published in Hoang, Mogavero and Laskar 2021), we address the problem from the astronomical perspective. We provide a large ensemble of solutions to quantify the chaotic diffusion of the fundamental frequencies of the Solar System, thereby offering a global view onto the astronomical forcing, particularly for the metronomes, the eccentricity grand cycles and the resonance transition.

The chapter 2 introduces basic notations and description of secular dynamics of the Solar System.

Chapter 2

Secular dynamics of planetary system

This chapter introduces basic notation and classical results for the secular Hamiltonian of the Solar System. A detailed development and discussion of the secular Hamiltonian is referred to (Laskar, 1991;[START_REF] Laskar | Stability of the planetary three-body problem[END_REF][START_REF] Morbidelli | Modern celestial mechanics: aspects of solar system dynamics[END_REF], see also [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF] for the model of forced inner Solar System.

General description

The Hamiltonian H of the planets of the Solar System can be written as:

H = N k=1 p k 2 2µ k -G m 0 m k r k + N k=1 N l=k+1 p k • p l m 0 -G m k m l |r k -r l | , (2.1)
where r k 's are the Heliocentric coordinates, p k 's their conjugated momenta, m 0 the Sun mass, m k the planet masses, µ k = m 0 m k /(m 0 + m k ) the reduced masses, and G the gravitational constant, N is the number of planets (N =8 for all the planets of the Solar System) ; index k indicates the planets in radially ascending order, from Mercury to Neptune. The two first terms on the right side of Eq. (2.1) make up the Keplerian Hamiltonian -H 0 ; the third term compensates the coordinate change from inertial coordinates to heliocentric coordinates, and the last one represents the gravitational interaction between any two planets. Because the Keplerian Hamiltonian is much larger than the interaction, the problem should be treated as the perturbation of integrable two-body Hamiltonians. It is then natural to use a set of canonical action-angle variables. An adapted set of such coordinates for planetary motions are Poincare's rectangular coordinates in complex form (λ k , Λ k , x k , -ix k , y k , -iȳ k ), where λ k is the mean longitude, Λ k = µ k (m 0 + m k )a k (a k is the semi-major axis), the bar over head represent complex conjugate and

x k = Λ k 1 -1 -e 2 k E j̟ k , (2.2)
y k = Λ k 1 -e 2 k (1 -cos I k )E jΩ k , (2.3)
where e k is the eccentricity, I k the inclination, ̟ k the longitude of the perihelion, Ω k the longitudes of the nodes, and E represents the exponential operator, j stands for the imaginary unit. The two-body Keplerian Hamiltonian H 0 can be expressed neatly:

H 0 = N k=1 -G(m 0 + m k ) 2 µ 3 k 2Λ 2 k . (2.4)
With only H 0 , we recover the two-body problem, where all variables are constant, except for the mean longitude λ k , which increases linearly in time with the rate n k = ∂H 0 ∂Λ k . It is also this variable that will be averaged in the secular approach at the first order in planetary masses,

H = 1 (2π) N 2π 0 dλ 1 . . . 2π 0 dλ N H = H 0 + H, (2.5)
where the secular Hamiltonian H in a Newtonian setting is given by:

H = - N k=1 N l=k+1 Gm k m l (2π) 2 2π 0 2π 0 dλ k dλ l |r k -r l | , (2.6)
which represents the total orbitally averaged interaction between two planets. The function |r kr l | -1 in Eq. (2.6) should then be developed in Poincare's variables. It should be noticed that the averaged over closed orbits of the third term in Eq. (2.1) amounts to zero. Because H 0 is constant, it could be left out and H is thus sufficient to describe the secular evolution of the planetary system. However, for a long-term dynamics, relativistic correction could be important, especially for Mercury; therefore they must be taken into account in the secular Hamiltonian (e.g., Will, 2018):

H = - N k=1 N l=k+1 Gm k m l (2π) 2 2π 0 2π 0 dλ k dλ l |r k -r l | - N k=1 3G 2 m 2 0 m k c 2 a 2 k 1 -e 2 k ,
(2.7)

Truncation of secular Hamiltonian

The secular Hamiltonian H of the entire Solar System, at first order in planetary masses, can be expanded in series of the complex Poincaré variables of the planets, i.e. (x i , xi , y i , ȳi ) N i=1 [START_REF] Laskar | Stability of the planetary three-body problem[END_REF]. The secular Hamiltonian can be then truncated in eccentricity and inclination at different degrees. Truncation at total degree 2n results in a polynomial Hamiltonian H 2n . At the lowest degree, H 2 describes an integrable Laplace-Lagrange dynamics. Its analytical solution can be obtained by a canonical transformation to the complex proper modes variables (u, v), with corresponding action-angle variables (X, χ; Ψ, ψ) such that (u = √ XE -jχ ; v = √ ΨE -jψ ). Their vector forms can be written as: X = (X 1 , . . . , X N ), Ψ = (Ψ 1 , . . . , Ψ N ), χ = (χ 1 , . . . , χ N ), and ψ = (ψ 1 , . . . , ψ N ). When expressed in these action-angle variables, the truncated Hamiltonian is a finite Fourier series: (2.8) where I = (X, Ψ) and θ = (χ, ψ) are the action and angle variables, respectively and k is the wave vector of a given harmonic.

H 2n (I, θ) = k Hk 2n (I)E jk•θ ,

Equations of motion

The Hamiltonian H 2n can be expressed in three different sets of canonical variables. In the complex Poincare's rectangular coordinates, the equations of motion for H 2n (x, x, y, ȳ):

dx dt = -j ∂ H 2n ∂ x ,
(2.9)

dy dt = -j ∂ H 2n ∂ ȳ , (2.10)
where x = (x 1 , . . . , x N ) and y = (y 1 , . . . , y N ). Similarly for the equations of motions of the secular Hamiltonian H 2n (u, ū, v, v) in complex proper mode variables, where :

du dt = -j ∂ H 2n ∂ ū , (2.11) dv dt = -j ∂ H 2n ∂ v , (2.12)
where u = (u 1 , . . . , u N ) and v = (v 1 , . . . , v N ). Finally, the equations of motion for the secular Hamiltonian in the action angle variables of the proper modes H 2n (I, θ):

dI dt = - ∂ H 2n ∂θ , (2.13) dθ dt = ∂ H 2n ∂I . (2.14)
In practice, the integration are often performed in the first two sets of the complex variables depending on the objective. If the objective is to simply integrate the secular dynamics at degree 2n, namely the Hamiltonian H 2n , then it is best to use Eqs. (2.9 -2.10), since the representation of H 2n is most compact in (x, x, y, ȳ). On the other hand, if the objective is to construct a particular secular Hamiltonian according to specific choice of harmonics, for the purpose of either simplification or investigation of the their impacts, the Eqs. (2.11 -2.12), in the complex proper mode variables (u, ū, v, v) should be used. Notwithstanding, the Hamiltonian in action angle variables H 2n (I, θ) provides better intuition to understand the secular dynamics. For example, the Lagrange-Laplace secular Hamiltonian at degree 2 in eccentricity and inclination is simply a linear combination of the action variables:

H 2 (I, θ) = -ω LL • I, (2.15)
where the coefficient vector ω LL = (g LL 1 , . . . , g LL N , s LL 1 , . . . , s LL N ) contains the dominant constant contributions of the fundamental secular frequencies of the Solar System. Thus at degree 2, the action I are constant and the angles θ rotate with constant rates ω LL . One can also see from Eq. (2.13) with the Fourier expansion of Hamiltonian (Eq. 2.8), that it is convenient for separating the variation of the action vector I due to each harmonics of wave vector k: (2.16) where c.c. denotes complex conjugate. As a result, the direction of variation is also k.

dI dt k = k -j Hk 2n (I)E jk•θ + c.c. ,
The cumulative variation of a harmonic, that is the integral of Eq. (2.16) with respect to time, depends on the amplitude of the harmonic as well as its proximity to resonance.

Forced secular inner Solar System

As previously discussed in the Sect. 1.2.1, the ISS and the OSS are markedly different.

The planets in the OSS, due to being substantially more massive than the inner planets, are not much perturbed by the ISS. The trajectories of the OSS are also very stable, they can be effectively modeled by quasi-periodic series [START_REF] Laskar | Large-scale chaos in the solar system[END_REF] with constant outer secular frequencies [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF](Laskar et al., , 2011a)). Besides, the dynamics of the ISS is absent of strong mean motion resonances. Basing on these observations, [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF] propose a forced secular model of ISS. In this model, trajectories of the outer planets pre-determined quasi-periodic finite series:

x k (t) = M k l=1 xkl E jm kl •ωot
(2.17)

y k (t) = N k l=1 ỹkl E jn kl •ωot (2.18)
for k ∈ {5, 6, 7, 8}, where t is time, ω o is the vector of constant outer secular frequencies ω o = (g 5 , g 6 , g 7 , g 8 , s 5 , s 6 , s 7 , s 8 ), xkl and xkl are constant amplitudes, M k and N k are the number of components in the approximating series, m kl and n kl are wave vectors. Substituting the Eqs. (2.17 -2.18) to the Eq. (2.7), the Hamiltonian of the forced secular inner Solar System is obtained:

H[(x k , y k ) k=1,4 , t] = H[(x k , y k ) k=1,4 , (x k = x k (t), y k = y k (t)) k=5,8 ], (2.19)
with which the total number of degrees of freedom of the dynamics is reduced to 8. The truncation process of this secular Hamiltonian is also similar to the Sec. 2.1.1 .

Lagrange-Laplace secular dynamics

Because the eccentricity and inclinations of the planets in the Solar System are relatively small, H can be developed as a power series of these variables. By keeping only the second degree of the secular Hamiltonian, we retrieve the linear Lagrange-Laplace (LL) Hamiltonian:

H LL = -(x ⊤ Ax + ȳ⊤ By), (2.20)
where two matrix A and B are given in the appendix (A). The equations of motion read:

dx dt = -j ∂ H LL ∂ x = jAx, (2.21) dy dt = -j ∂ H LL ∂ ȳ = jBy. (2.22)
Resolving the Lagrange-Laplace dynamics is important, since its integrability allows us to define the proper variables (i.e. u, v), and their corresponding action-angle variable (I, θ), which were both discussed in the previous section. The Eqs. (2.21 -2.22) can be solved by diagonalizing the matrix through a linear transformation:

x = S a u, (2.23) y = S b v, (2.24)
where S a and S b are the real matrices whose columns are eigenvectors of A and B respectively. With these new variable, the equations of motion read:

du dt = jD a u, D a = S a -1 AS a (2.25) dv dt = jD b v, D b = S b -1 BS b . (2.26) (2.
du k dt = jg k u k , dv k dt = js k v k (2.28)
will have solutions:

u k = u k (0)E jg k t , v k = v k (0)E js k t , (2.29)
where u k and v k are the proper modes of the Solar System. The complete solutions will then be superpositions of such proper modes.

x i (t) = N k=1 S a ik u k (t) = N k=1 S a ik u k (0) exp(jg k t),
(2.30)

y i (t) = N k=1 S b ik v k (t) = N k=1 S b ik v k (0) exp(js k t).
(2.31) One of the mode, which is historically chosen as s 5 , is put to zero as a consequence of the conservation of total angular momentum of the system. From Eqs. (2.30 -2.31) and the definition of the Poincaré coordinates in Eqs. (2.2 -2.3), it is clear that eccentricity cycles are made up of g k , and the obliquity cycles of s k . These frequencies are constant in Laplace-Lagrange dynamics, however, in models of higher degree in eccentricity and inclination, the non-linearity, despite being small, can induce chaotic behavior. The fundamental frequencies are no longer constant in the chaotic model, but change with time.

Lagrange-Laplace dynamics of the forced inner Solar System

In the forced secular ISS, the first truncation of Hamiltonian (Eq. 2.19) at degree two in eccentricity and inclination constitutes the integrable Lagrange-Laplace dynamics. For simplicity the inclination part will be omitted, the LL Hamiltonian reads:

H LL = -x ⊤ Ax = -[x ⊤ in , x⊤ out ] A in A int A ⊤ int A out x in x out , = -x ⊤ in A in x in -x⊤ in A int x out -x⊤ out A ⊤ int x in -x⊤ out A out x out (2.32)
where

x in = [x 1 , . . . , x 4 ], x out = [x 5 , . . . , x 8 ],
and A is specified in the appendix A. The Hamiltonian of the heavy outer Solar System could then be uncoupled, and reads

H out = -x ⊤ out A out x out , (2.33)
which then have solutions of Eq. (2.30). More generally, we can assume that the outer planets follows a pre-constructed quasi-periodic motion (Eq. 2.17 -2.18). The Hamiltonian of the inner Solar System now reads: (2.34) where the last two terms represent the forcing of the giant planets. The equation of motion of the inner Solar System is

H in = -x ⊤ in A in x in -x⊤ in A int x out (t) -x⊤ out (t)A int x in ,
dx in dt = -j ∂ H in ∂ xin = jA in x in + jA int x out (t). (2.35)
This equation is solved similarly as the original LL Eq. (2.21). The matrix A in is diagonalized with the new variables u in = S ⊤ x in , where S is a transformation matrix made of the eigenvectors of A in . The Eq. (2.35) now reads:

du in dt = iD in u in + jA ′ int x out (t), (2.36)
where D in is a diagonal matrix of the eigenvalues of A in , that are (g LL 1 , g LL 2 , g LL 3 , g LL 4 ), and A ′ int = S ⊤ A int is the new interaction matrix. The Eq. (2.36) have the solutions:

u in (t) = u free (t) + u forced (t), (2.37)
where u free (t) = exp(jtD in )u(0), (2.38)

and u forced (t) = exp(jtD in ) t 0 ds exp(-jsD in )A ′ int x out (t).
(2.39)

Finally, the complete LL solution of the forced inner Solar System is:

x in = S(u free (t) + u forced (t)).
(2.40)

Introduction

Milankovitch (1941) hypothesized that some of the past large climate changes on the Earth originated from the long-term variations in its orbital and rotational elements. These variations are imprinted along the stratigraphic sequences of sediments. Using their correlations with an orbital solution [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF], 2011a[START_REF] Laskar | Astrochronology[END_REF], some of the geological records can be dated with precision. This method, named astrochronology, has become a standard practice in the stratigraphic community and has proven to be a powerful tool for reconstructing the geological timescale (e.g., [START_REF] Gradstein | A new geologic time scale, with special reference to Precambrian and Neogene[END_REF][START_REF] Gradstein | The geologic time scale[END_REF][START_REF] Gradstein | The chronostratigraphic scale[END_REF]. The climate rhythms found in the geological records are directly related to the Earth's precession constant and to the fundamental secular frequencies of the Solar System: the precession frequencies (g i ) i=1,8 of the planet perihelia and the precession frequencies (s i ) i=1,8 of their ascending nodes. The evolution of these fundamental frequencies is accurately determined up to 60 Myr [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF], 2011a[START_REF] Laskar | Astrochronology[END_REF]. Beyond this limit, even with the current highest precision ephemerides, it is hopless to obtain a precise past history of the Solar System simply via numerical integration. This limit does not lie in the precision of the determination of the initial conditions but originates in the chaotic nature of the Solar System [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF], 1990, Laskar et al., 2011a). However, because the astronomical signal is recorded in the geological data, it appears to be possible to trace back the orbital variations in the geological record and thus to constrain the astronomical solutions beyond their predictability horizon [START_REF] Olsen | Long-period Milankovitch cycles from the Late Triassic and Early Jurassic of eastern North America and their implications for the calibration of the Early Mesozoic time-scale and the long-term behaviour of the planets[END_REF][START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF][START_REF] Olsen | Mapping Solar System chaos with the Geological Orrery[END_REF]. Nevertheless, a deterministic view of the Solar System is no longer reliable beyond 60 Myr, and a statistical approach should be adopted. Geological constraints should likewise be retrieved in a statistical setting. In this spirit, a recent Bayesian Markov chain Monte Carlo (MCMC) approach has been proposed to provide geological constraints on the fundamental frequencies [START_REF] Meyers | Proterozoic Milankovitch cycles and the history of the solar system[END_REF]. For such a Bayesian approach to give any meaningful constraint, proper prior distributions of the fundamental secular frequencies are required, and therefore a statistical study of the orbital motion of the Solar System planets is needed. This constitutes the motivation for the present study. [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF] performed the first statistical analysis of the long-term chaotic behavior of the planetary eccentricities and inclinations in the Solar System. [START_REF] Mogavero | Addressing the statistical mechanics of planet orbits in the solar system[END_REF] reconsidered the problem from the perspective of statistical mechanics. Our study is a follow-up of [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF]. We study fundamental frequencies instead of orbital elements because they are more robust and are closer to the proxies that can be traced in the geological records. This study is based on the numerical integrations of 120 000 different solutions of the averaged equations of the Solar System over 500 Myr, 40 000 of which were integrated up to 5 Gyr. The initial conditions of the solutions are sampled closely around a reference value that is compatible with our present knowledge of planetary motion.

Dynamical model

Secular equations

We used the secular equations of motions of [START_REF] Laskar | Accurate methods in general planetary theory[END_REF](Laskar, , 1990(Laskar, , 2008, and references therein), and references therein). They were obtained via series expansions in planetary masses, eccentricities, and inclinations as well as through second-order analytical averaging over the rapidly changing mean longitudes of the planets. The expansion was truncated at the second order with respect to the masses and to degree 5 in eccentricities and inclinations. The equations include corrections from general relativity and Earth-Moon gravitational interaction. This leads to the following system of ordinary differential equations:

dω dt = √ -1{Γ + Φ 3 (ω, ω) + Φ 5 (ω, ω)}, (3.1)
where

ω = (z 1 , . . . , z 8 , ζ 1 , . . . , ζ 8 ) with z k = e k exp(̟ k ) and ζ k = sin(i k /2) exp(Ω k ).
The variable ̟ k is the longitude of the perihelion, Ω k is the longitude of the ascending node, e k is eccentricity, and i k is inclination. The function Φ 3 (ω, ω) and Φ 5 (ω, ω) are the terms of degree 3 and 5. The 16 × 16 matrix Γ is the linear Laplace-Lagrange system, which is slightly modified to make up for the higher-order terms in the outer Solar System. With an adapted initial condition, the secular solution is very close to the solution of direct integration over 35 Myr [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]. The major advantage of the secular system over direct integration is speed. Numerically integrating averaged equations is 2000 times faster than non-averaged ones due to the much larger step size: 250 years instead of 1.8265 days. It is thus desirable to employ the secular equations to study the statistics of the Solar System. However, we also compare their predictions to those of a non-averaged comprehensive dynamical model in Sect. 3.4.

Frequency analysis

We employed the frequency analysis (FA) technique proposed by [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency analysis for multi-dimensional systems. Global dynamics and diffusion[END_REF] to extract the fundamental secular frequencies from the integrated solutions. The method finds a quasi-periodic approximation f ′ (t) = N k=1 a k e iν k t of a function f (t) over a time span interval [0, T ]. It first finds the strongest mode, which corresponds to the maximum of the function:

φ(σ) = f (t)|e iσt = 1 T T 0 χ(t)f (t)e -iσt dt, (3.2)
where χ(t) is a weight function that improves the precision of the maximum determination; it was chosen to be the Hanning window filter, that is, χ(t) = 1 + cos(πt/T ). The next step is the Gram-Schmidt orthogonalization. The complex amplitude a 1 of the first frequency ν 1 is calculated via the orthogonal projection of the function f (t) on e iν 1 t . This mode is then subtracted from the function

f (t) to get a new function, f 1 (t) = f (t) -a 1 e νt .
The process is then repeated with this newly obtained function until N desired strongest modes are obtained. This technique works very well for weakly chaotic systems such as the Solar System when variables can be decomposed into quasi-periodic modes over a sufficiently short period of time. It has been proven that this algorithm converges toward the true frequencies much faster than the classical fast Fourier transform [START_REF] Laskar | Frequency map analysis and quasiperiodic decompositions[END_REF]. Therefore, it is a good tool for studying the chaotic diffusion of the fundamental frequencies. In this work we used a routine naftab written in the publicly available computer algebra software TRIP [START_REF] Gastineau | Trip: a computer algebra system dedicated to celestial mechanics and perturbation series[END_REF], developed at IMCCE, to directly apply the frequency analysis.

To extract the fundamental secular frequencies of the Solar System, we applied the FA to the proper variables (z

• i , ζ • i ) i=1,8
of the secular equations (Laskar, 1990). Each fundamental frequency is obtained as the frequency of the strongest Fourier component of the corresponding proper variable. To track the time evolution of the frequencies, the FA was applied with a moving interval whose sliding step was 1 Myr. The interval sizes were 10 Myr, 20 Myr, and 50 Myr.

Estimation of probability density functions

The samples of this study consist of the secular frequencies of the astronomical solutions that were obtained by integrating the secular equations (Eq. 3.1) from very close initial conditions. Due to this initial proximity, the correlation of the solutions in the samples lasts for a long period of time but will eventually diminish. Our objective is to obtain a robust estimation of the marginal probability density functions (PDFs) from these correlated samples. In fact, this correlation is the main motivation for our use of the estimation methods in this section. Details of our samples are described in the first part of Sect. 3.4. We used kernel density estimation (KDE) to estimate the time-evolving marginal PDFs of the fundamental frequencies of the Solar System. In addition, the statistical uncertainty of our density estimations (i.e., PDF estimations) was measured by the moving block bootstrap (MBB) method. To our knowledge, this application of MBB for the KDE of a

Moving block bootstrap

Since the seminal paper by [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF], bootstrap has become a standard resampling technique for evaluating the uncertainty of a statistical estimation. Bias and variance errors of KDE with the choice of bandwidth BW 1/5 (Eq. 3.4) are on the same order of magnitude, and hence one should either under-smooth (i.e., choose a smaller bandwidth) to minimize the bias (Hall 1995) or use an explicit bias correction with appropriate Studentizing [START_REF] Cheng | Nonparametric inference via bootstrapping the debiased estimator[END_REF][START_REF] Calonico | On the effect of bias estimation on coverage accuracy in nonparametric inference[END_REF]. However, naively applying the i.i.d. bootstrap procedure on dependent data could underestimate the true uncertainty because all the dependence structure would be lost just by "scrambling" data together [START_REF] Hart | Some automated methods of smoothing time-dependent data[END_REF]. To remedy the problem, [START_REF] Kunsch | The jackknife and the bootstrap for general stationary observations[END_REF][START_REF] Liu | Moving blocks jackknife and bootstrap capture weak dependence[END_REF] independently devised the MBB, which later became standard practice for evaluating the uncertainty in dependent data (see Kreiss and Lahiri 2012 for a review). Although the MBB for smooth functional has been intensively studied, the literature on MBB for the KDE of dependent data is very limited. Recently, [START_REF] Kuffner | Block bootstrap optimality for density estimation with dependent data[END_REF] formulated an optimality theory of the block bootstrap for KDE under an assumption of weak dependence. They proposed both under-smoothing and an explicit bias correction scheme to obtain the sampling distribution of the KDE. However, good tuning parameters, which are generally difficult to find if the data are from an unknown distribution, are required to provide a decent result. In this paper we propose overcoming this problem with an inductive approach: The optimal tuning parameters obtained in a known model are tested on different models and then extrapolated to the subject of our study, the Solar System.

Procedure of MBB

We briefly describe the under-smooth MBB for the KDE method (a more detailed description can be found in [START_REF] Kuffner | Block bootstrap optimality for density estimation with dependent data[END_REF]. We suppose that X = {X 1 , X 2 , . . . , X n } is a dependent sample of a mixing process with an underlying density function p(x). The KDE of the sample is ph = ph (x|X); the hat above a given quantity denotes its estimated value. We used MBB to estimate the distribution of δ(x) = ph (x) -p(x), where x ∈ Ω and Ω is the domain of interest. Let l be an integer satisfying 1 ≤ l ≤ n. Then B i,l = {X i , X i+1 , . . . , X i+l-1 } with i ∈ {1, . . . , n -l + 1} denotes all the possible overlapping blocks of size l of the sample. Supposing, for the sake of simplicity, that l divides n, then b = n/l. The MBB samples are obtained by selecting b blocks randomly with replacement from {B 1,l , . . . , B n-l+1,l }. Serial concatenation of b blocks will give n bootstrap observations:

X * l = {B * 1,l , . . . , B * b,l }.
By choosing sufficiently large values of l (preferably larger than the correlation length), the MBB sample can retain the structure of the sample dependence. The schematic diagram of the MBB for a sample of n data point and block size l = 3 is shown in Fig. 3.2. For k > 0, the KDE of the bootstrap sample is p * k,l = pk (x|X * l ) and its expectation is

E[p * k,l ] = pk (x|B 1,l , . . . , B n-l+1,l ). We define δ * k,l (x) = k h (p * k,l -E[p * k,l ]) (3.5)
such that if h is chosen properly to reduce the bias to be asymptotically negligible with respect to the stochastic variation, then the MBB distribution P (δ * k,l (x)|X) is a consistent estimator of the error distribution P (δ(x)) when h → 0, nh → ∞, k → 0, lk → ∞, and n/l → ∞. We note that if l = 1, then k = h and MBB reverts to the under-smoothing procedure for the i.i.d sample studied by Hall et al., 1995. The efficiency of this estimator depends sensitively on two tuning parameters, l and k. We are interested in the uncertainty of the KDE, which is characterized by the confidence interval CI 1-α (x) and the confidence band CB 1-α , which are defined as:

P (|δ(x)| < CI 1-α (x)) = 1 -α, (3.6) P (|δ(x)| < CB 1-α ∀x ∈ Ω) = 1 -α, (3.7)
where α denote the level of uncertainty; for example, α = 0.05 denotes 95% CI. They can be estimated by the MBB distribution P (δ * k,l (x)|X) as:

P (|δ * k,l (x)| < CI 1-α (x)) = 1 -α, (3.8) P (|δ * k,l (x)| < CB 1-α ∀x ∈ Ω) = 1 -α. (3.9)
In this paper we also use CB and CI without the hat overhead to denote estimated values.

Our choice of the parameters of the MBB procedure, l and k, is based on the effective sample size n eff , defined as [START_REF] Kass | Markov chain Monte Carlo in practice: a roundtable discussion[END_REF]:

n eff = n 1 + 2 ∞ k=1 ρ(k) , (3.10)
where ρ(k) is the sample autocorrelation of lag k. The block length l is chosen by the sample correlation size l corr as l = l corr := n n eff , (3.11) and the bootstrap bandwidth is parametrized as

k = h(c 0 + (1 -c 0 )l -γ corr ), (3.12)
where γ and c 0 are two optimizing constants. The reason for this choice of parametrization is twofold. First, when l corr → 1, the sample become independent, and then k → h. Secondly, the rate of change of k with respect to l should be greater when l corr is small than when it is large. Therefore, when l corr ≫ 1, the optimal value of k should be quite stable. We also observe experimentally that the optimal value of k is indeed relatively stable at around 2h as long as l = l corr ≫ 1. So we simply chose γ = 1 and c 0 = 2. This choice of parameters turns out to be quite robust, as demonstrated by the two numerical experiments in Sect. 3.3.3. The literature on KDE and MBB focuses on stationary, weakly dependent sequences. The data in our case, however, are different: They are not strictly stationary; the formulation of the mixing condition might be different [START_REF] Hang | Kernel density estimation for dynamical systems[END_REF]; the correlation in the sample is not constant but evolving with time; finally, and most importantly, the data structure is different. Our sample units, which are the orbital solutions, are ordered based on their initial distances in phase space. The solutions evolve over time but their order remains unchanged. Therefore, statistical notions such as correlation and stationarity should be considered within this framework for the sample at fixed values of time. Because of the differences presented, an optimality theory, which is not currently available, might be needed for this case. However, we assume that it would not differ significantly from the orthodox analysis and that a decent working MBB procedure could be obtained with some good choice of parameters. This is tested in the section below.

Application to the Solar System

The goal of this paper is to have a consistent statistical description of the propagation of the dynamical uncertainty on the fundamental secular frequencies of the Solar System induced by its chaotic behavior: that is, simply put, to obtain their time-evolving marginal PDF. We first sampled the initial orbital elements of the Solar System planets that were close to the reference value, and then numerically integrated the secular equations (Eqs. 3.1) from those initial conditions to obtain a sample of orbital solutions. Kernel density estimation was then used to estimate the marginal PDF of the frequencies of the sample at a fixed value of time, and finally the MBB method was applied to estimate the uncertainty of the density estimation.

The evolution of our sample can be divided into two stages: Lyapunov divergence and chaotic diffusion. In the first stage, because the initial density is extremely localized around the reference value, all solutions essentially follow the reference trajectory; the difference between the solutions is very small but diverges exponentially with a characteristic Lyapunov exponent of ∼ 1/5 Myr -1 [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF]. The solutions in the first stage are almost indistinguishable and the correlation between them is so great that regardless of how many solutions in the sample we integrated, the effective size of the sample is close to one. The second stage begins when the differences between the solutions are large enough to become macroscopically visible. The Lyapunov divergence saturates and gives place to chaotic diffusion. The correlation between the solutions starts to decrease, the distribution of the sample settles, and the memory of the initial conditions fades. It can take several hundred million years for the sample to forget its initial configuration. Contrary to the exponential growth in the first stage, the dispersion of the samples expands slowly with a power law in time (see Fig. 3.12). The time boundary between the two stages depends on the dispersion of the initial conditions:The wider they are, the faster the second stages come, and vice versa. If they are chosen to represent the uncertainty of the current ephemeris, the second stage should take place around 60 Myr in the complete model of the Solar System (Laskar et al., 2011a;[START_REF] Laskar | Strong chaos induced by close encounters with Ceres and Vesta[END_REF].

In this section we focus on the statistical description of the fundamental frequencies of the Solar System in the second stage. The aim is to obtain a valid estimation of time-evolving PDFs of the frequencies beyond 60 Myrs. However, the PDF evolution generally depends on the choice of initial conditions. Moreover, the simplification of the secular equations compared to the complete model of the Solar System could, a priori, provide results that are not sufficiently accurate.

Choice of initial conditions

For a complete model of the Solar System, the initial conditions should be sampled in such a way that they are representative of the current planetary ephemeris uncertainty. Nevertheless, for a simplified secular model (Eq. 3.1), the difference from the complete model is greater than the ephemeris uncertainty. An optimized secular solution follows the complete solution initially but departs from it long before 60 Myr [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]. A direct adaptation of the current planetary ephemeris uncertainty to the initial conditions of the secular model can thus be misleading. Therefore, we adopted a more cautious approach, that is, to study first the effect of sampling the initial conditions on

Solutions

Offsets ǫ X i -5000 ǫ to 5000 ǫ 10 -10 Y i -5000 ǫ to 5000 ǫ 10 -8 Z i -5000 ǫ to 5000 ǫ 10 -11

Table 3.1 -Offsets of the initial eccentricities of the four planets: {Mercury, Venus, Earth, Mars}, which corresponds to i = {1, 2, 3, 4}. Different integrations correspond to offsets of N ǫ in eccentricity of a single planet for N = -5 × 10 3 , . . . , +5 × 10 3 , while other variables are kept to their nominal values.

the PDF estimation.

Initial conditions can be sampled in many ways, especially in a high-dimensional system such as the Solar System. Our choice of initial conditions is quite particular, but they encompass different possible ways of sampling initial conditions (Table 3.1). There were three batches, {X i }, {Y i }, and {Z i }, which correspond to three different variation sizes ǫ of initial conditions. Each batch was composed of four different sets of samples. Each set contained 10 000 initial conditions, where the eccentricity of the associated planet is linearly spaced from the reference value, with the spacing ǫ corresponding to the batch it belongs to. We then integrated the secular equations (Eqs. 3.1) from these initial conditions 500 Myr into the past and 500 Myr into the future. For the batch {Z i }, the integration time is 5 billion years in both directions. The frequencies were then extracted using FA (Sect. 3.2.2). It should be noted that we do not aim to obtain the joint probability distribution of all the fundamental frequencies, but rather their individual marginal PDFs (i.e., the PDF of one frequency at a time). The marginal PDFs of the frequencies were estimated by the KDE with the rule-of-thumb bandwidth (BW 1/5 ); upper bounds of their 95% CIs are measured by the MBB method with the bandwidth h = BW 1/3 and the optimized parameters (Eqs. 3.11-3.12) from 1 000 MBB samples.

We first compare the evolution of the density of the four sets in each batch in Sect. 3.4.2, and the second test is performed to compare the statistics between the batches in Sect. 3.4.3. The robustness of the secular statistics is assessed by these two tests, which additionally shed light on the initial-condition-dependence aspect of the statistics. All the density estimations from these sets are compared with those of the 2 500 complete solutions obtained in the previous work of [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF] to test the accuracy of the secular statistics. It should be recalled that this numerical experiment needed 8 million hours of CPU time, the output of which was saved and could thus be used in the present study.

When comparing two sets of different sizes, because the rates of divergence in the first stage are similar, the wider set reaches the chaotic diffusion phase faster than the more compact set; hence, it is essentially diffusing ahead for a certain time in the second stage. Therefore, in order to have a relevant comparison, a proper time shift was introduced to compensate for this effect. We shifted {X i } and {Y i } ahead by 30 Myr and 20 Myr, respectively, while keeping the time of {Z i } as reference. This choice was motivated by the fact that the transition to the chaotic diffusion of {Z i } is around 50 -60 Myr, which is indicated by the direct integration of the Solar System (Laskar et al., 2011a;[START_REF] Laskar | Strong chaos induced by close encounters with Ceres and Vesta[END_REF].

First test: Different samples of the same variation size ǫ

Comparing the density estimation of different sets in the same batch was the first test of prediction robustness from our secular model. The evolution of the density estimations are sensitive to how initial conditions are sampled, and therefore this initial-condition sensitivity must be quantified for a valid prediction. In this first test we compared the time-evolving PDFs whose initial conditions are sampled with the same variation size ǫ but in different variables.

The result is quite clear. The different sets of the same batch slowly lose the memory of their initial differences due to chaos and then converge toward the same distribution. This convergence is illustrated by Fig. 3.7, which shows that the density estimation of (g) i=1,4 of the four sets of the batch {Z i } nearly overlap with one another at 150 Myr in the future. The rates of convergence of different batches are different. So although {X i } and {Y i } exhibit the same behavior as {Z i }, they converge differently with disparate rates: At around 100-150 Myr in the future, the density estimations of the frequencies of {X i } nearly overlap with one another; this occurs at 150-200 Myr for {Y i }, depending on the frequency. Interestingly, for the samples that are integrated in the past, the rates of convergence are higher and the overlap generally happens at around -100 Myr (see Fig.

Second test: Different samples of different variation sizes

Comparing the density estimation of the three batches, {X i }, {Y i }, and {Z i }, was our second test of robustness. Although the initial conditions of the three batches were varied around the same reference values, the ways they were sampled were different since the variation sizes were different. Differences in the initial variation sizes mean that the batches enter the diffusion stage at different times and also at different points in the phase space, so that the convergence between batches, if it exists, takes longer. The result of our test is summarized by the density estimation of the frequencies at two times in the past, -100 Myr and -200 Myr (Fig. 3.8). At -100 Myr, the density estimations from different sets of each batch cluster around one another as described in the previous test. Each batch forms a cluster of density estimations, and the differences between the three clusters are noticeable. Moreover, the estimation uncertainty, depicted by the colored band, is quite large. Fast-forward 100 Myrs of chaotic mixing: at -200 Myr, density estimations of the frequencies spread out, estimation uncertainty shrinks, and, most importantly, differences between the three batches are much smaller and continue to diminish even further with time. In the opposite time direction, the same phenomenon is observed but the rate of convergence between the batches is slower (Fig. 3.9). At 150 Myr in the future, the density estimations of the frequencies of {X i } and {Z i } have practically converged and those of {Y i } are still trying to, and yet differences between the three batches are noticeable. However, the estimations of all 12 sets from the three batches nearly overlap with one another at 350 Myr, which demonstrates that the effect of different initial samplings vanishes via chaotic mixing. It should be noted that when looking at some specific properties of the PDF, such as means and variances, the differences between the sets are small. For example, the differences between the means of the PDF estimation of the 12 sets are generally smaller than 0.1 ′′ /yr for most of the fundamental frequencies at -100 Myr; at -200 Myr, the differences diminish to twice as small. order as the variability of the results from the same secular model. These differences in the means of the PDFs are generally smaller than 0.05 "/yr at 350 Myr. Although some minor differences between the two models are still visible, especially in s 1 for example, these differences diminish with time. This convergence between the two models strongly suggests the compatibility of the secular system with the direct integration of the realistic model of the Solar System used in [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF].

A complementary test on frequency analysis

The fundamental frequencies of the Solar System are central in our work, and the method to obtain them is thus essential. In this section we briefly examine the FA method (Sect. 3.2.2). The FA method searches a quasi-periodic approximation of the solution of the Solar System with constant frequencies over a time window ∆t. So unique frequencies at time t are extracted from an oscillating sequence in the time interval [t -∆t/2, t + ∆t/2]. For a quasi-periodic solution, the longer we choose the time interval ∆t to be, the more accurate the extracted frequencies are. Nevertheless, the fundamental frequencies of the Solar System are expected to vary over a few Lyapunov times, that is, over 5 Myr. Therefore, we have a trade-off between the extraction accuracy and the time variation of the frequencies when choosing ∆t: When ∆t is too large, the obtained frequency will tend to be the average of its variation over the same period. We chose ∆t = 20 Myr as the standard FA interval. In some circumstances that require the detection of rapid changes in frequencies, such as the resonance transition, a smaller ∆t is more favorable (see Sect. 3.6.3).

The extracted frequencies are sensitive to the choice of ∆t, and yet their density estimation is relatively robust. Figure 3.10 compares the density estimation of the eccentricity frequencies (g i ) 1,4 , which are extracted via FA with three different ∆t at two different times. The differences are generally small but still notable for g 2 and g 4 at 150 Myr, and they diminish with time.

Parametric fitting

From Sect. 3.4 we see that the long-term PDFs of the secular frequencies possess a distinct Gaussian-like shape that flattens as time passes. It is interesting to approximate these densities by a simple parametric model, such that its parameters can characterize the shape of the density and summarize its evolution. The model with its fitting parameters can also be used as an approximation of the numerical densities for later application. For this purpose, we used the density estimation of the secular frequencies of the Solar System from the batch {Z i }, which is composed of 40 000 different orbits over 5 Gyr. The inner fundamental frequencies (i.e., the frequencies of the inner planets, (g i , s i ) i=1,4 ) are obtained by FA over an interval of 20 Myr. For the outer fundamental frequencies (i.e., the frequencies of the outer planets, (g i ) i=5,8 , (s i ) i=6,8 ), the FA interval is 50 Myr. The frequency s 5 is zero due to the conservation of the total angular momentum of the Solar System. Table 3.2 -Linear and power-law fits for the time evolution of the parameters (Fig. 3.12)

µ 0 ["/yr] a ["/yr] 2 b α 0 µ 1 ["/yr] σ 2 1 ["/yr] 2 A
of the skew Gaussian mixture model (Eq. 3.19) for the fundamental frequencies of the Solar System. Column 1 contains the considered secular frequencies. In Col. 2 we show the linear fits of µ 0 that represent the center of the distribution. The power-law fit of σ 2 0 has the form σ 2 0 (t) = aT b (Eq. 3.20), where T = t/(1 Gyr) and a and b are given in Cols. 3 and 4, respectively. Linear fits of the skewness parameter α 0 are given in Col. 5. The last three columns show linear fits of the secondary mode of g 4 and s 3 (Eq. 3.19). The parameter σ 2 0 is fitted from 200 Myr to 5 Gyr in the past, while all the others are fitted from 500 Myr. [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF] found that the 250 Myr-averaged marginal PDFs of the eccentricity and inclination of the inner Solar System planets are described quite accurately by the Rice distribution, which is essentially the distribution of the length of a 2D vector when its individual components follow independent Gaussian distributions. In in our case, the density estimations of the fundamental frequencies of the Solar System resemble Gaussian distributions, but many of them get skewed as time passes; this is especially true for the inner frequencies. To account for this skewness, we propose the skew normal distribution as the fitting distribution to the density estimation of the frequencies:

Skew Gaussian mixture model

f µ 0 ,σ 0 ,α 0 (x) = 2 σ 0 φ x -µ 0 σ 0 Φ α x -µ 0 σ 0 , (3.17)
where α is the parameter characterizing the skewness, φ(x) denotes the standard normal probability density distribution with mean µ 0 and standard deviation σ 0 , and Φ(x) is its cumulative distribution function given by .18) where erf denotes the error function. Some of the frequencies, interestingly, have several secondary modes in their density estimation apart from their primary one. Most of the secondary modes, if they exist, are quite small compared to the primary one. They are also often short lived; most of them emerge at the beginning of the diffusion stage and disappear quickly thereafter. much smaller than the primary mode: A 0 ≫ m i=1 A i . In our case here, m = 1 for g 4 and s 3 , while the asymptotic secondary modes for the other frequencies can be considered negligible.

Φ(x) = x -∞ φ(t) dt = 1 2 1 + erf x √ 2 , ( 3
For the outer fundamental frequencies, we do not observe significant skewness or secondary modes in their density estimations. Therefore, the density estimations of the outer fundamental frequencies can be approximated by a simple Gaussian distribution.

The density estimations of the frequencies are shown at three well-separated times in Fig. 3.11. The diffusion is clearly visible as the density estimations get more and more disperse over time. Moreover, the density estimations get more skewed as time goes by. For g 3 and g 4 , the skewness of the density estimations is small, so we assume α 0 = 0 (Eq. 3.17) for these frequencies for the sake of simplicity. The fundamental frequencies of the outer planets are very stable. Their variations are much smaller than those of the inner frequencies. Taking as an example the most unstable of the outer frequencies, g 6 , its standard deviation at -4 Gyr is only about 2 • 10 -3 "/yr, while that of the most stable inner frequency, g 2 , is 15 times larger. Therefore, when considering a combination involving an inner fundamental frequency and an outer one, the latter can be effectively regarded as a constant.

The result of our fitting model is also plotted in Fig. 3.11. It is remarkable that the density estimations of the frequencies are well approximated by the fitting curve over three different epochs. It should be noticed that the base of our fitting model -the skew Gaussian distribution -only has three parameters. Additional parameters are only needed for some frequencies, for example g 4 and s 3 . Nevertheless, such additional parameters only account for the minor features, and three parameters are sufficient to represent the bulk of the density estimations over a long timescale.

Evolution of the parameters

The parameters of our fitting models are extracted by the method of least squares, implemented by the routine curve_fit in the scipy package in Python. To retrieve the statistical distribution (mean and standard deviation) of the parameters of a given model, we implemented a bootstrap approach based on Eq. (3.16), with the assumption that pointwise standard errors of the KDE estimated by MBB are independent. We remark that, with such an assumption, the variance of the fitting parameters tends to be underestimated.

The time evolution of the mean of the parameters of the fitting models is shown in Fig. 3.12, along with their ±2 standard error, for both the past and the future. It turns out that the evolution of σ 2 0 is robustly fitted from 200 Myr to 5 Gyr in the past by the power-law function σ 2 0 (t) = a T b , (3.20) where T = t/(1 Gyr), as shown in Fig. 3.12. For all the other parameters, we performed a linear fit. All these fits are summarized in Table 3.2.

The differences between past and future evolutions are small and generally tend to decrease with time. Therefore, the fit for the parameters in the past given in Table 3.2 should also be representative of their future evolution. In general, the parameters follow relatively smooth curves with distinct tendencies. The skewness parameters α 0 , increasing in absolute value, show that the PDFs of the inner fundamental frequencies get more and more skewed over time. The center of the distributions, indicated by µ 0 , does not change significantly compared to α 0 and σ 2 0 (Fig. 3.12). The secondary modes of g 4 and s 3 are also quite stable.

The diffusion of the frequencies is quantified by the increasing σ 2 0 , which is closely linked to their distribution variance. As the exponents b of the power laws in Table 3.2 are different from unity, the chaotic diffusion of the fundamental frequencies turns out to be an anomalous diffusion process. Interestingly, all the inner frequencies clearly undergo subdiffusion, that is, the exponent of the power law b is smaller than 1 (b = 1 corresponding to Brownian diffusion). Therefore, an extrapolation of the variance of the inner frequencies based on the assumption of linear diffusion over a short time interval would generally lead to its overestimation over longer times. On the contrary, the exponents b are either smaller or larger than unity for the outer frequencies. It should be noted that, because the variations in the outer frequencies are very small, the value of the corresponding exponents b might be overestimated due to the finite precision of FA.

Geological application

The aim of this project is to have a reliable statistical picture of the secular frequencies of the Solar System beyond 60 Myr. It is interesting to put recent geological results in this astronomical framework. First, it can be used as a geological test of our study, and secondly the application provides a glimpse of how astronomical data could be used in a cyclostratigraphy study. We first show how the uncertainty of a widely used dating tool in astrochronology can be quantified, then we apply our work to two recent geological results, which are from the Newark-Hartford data [START_REF] Olsen | Mapping Solar System chaos with the Geological Orrery[END_REF] and the Libsack core [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF].

Astronomical metronomes

Although it is not possible to recover the precise planetary orbital motion beyond 60 Myr, some astronomical forcing components are stable and prominent enough such that they can be used to calibrate geological records in the Mesozoic Era or beyond (see Laskar 2020 for a review). The most widely used is the 405 kyr eccentricity cycle g 2 -g 5 , which is the strongest component of the evolution of Earth's eccentricity. The inclination cycle s 3 -s 6 has also recently been suggested for the time calibration of stratigraphic sequences [START_REF] Boulila | Towards a robust and consistent middle Eocene astronomical timescale[END_REF][START_REF] Charbonnier | Obliquity pacing of the hydrological cycle during the Oceanic Anoxic Event 2[END_REF]. Although s 3 -s 6 is not the strongest among the obliquity cycles, it is quite isolated from other cycles and thus easy to identify [START_REF] Laskar | Astrochronology[END_REF]. The main reason that g 2 -g 5 and, possibly, s 3 -s 6 can be used as astronomical metronomes is their stability. Indeed, their uncertainty has to be small to be reliably used for practical application.

In previous work, the uncertainty of the frequency combination was derived from the analysis of a few solutions [START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF](Laskar et al., , 2011a;;[START_REF] Laskar | Astrochronology[END_REF]. Here we used a much larger number of solutions and the KDE-MBB method to derive both the PDF of the frequencies and their statistical errors. Starting from the results of all the sets of orbital solutions (Sect. 3.4), we produced the compound density estimation of the fundamental frequencies and their relevant combinations following the conservative approach in Eq. (3.16). These data are archived with the paper. are all consistent with our density estimation as there is a non-negligible possibility of finding a secular solution that agrees with the geological data. It should be noted that certain frequencies are significantly correlated, g 3 and g 4 for example. We cannot assume that they are independent; therefore, we have to calculate the density estimation of their combination directly.

Given the unavailability of the uncertainty in geological frequencies, the probability of finding the geological frequencies in a numerical orbital solution cannot be obtained directly. However, we can use La2010d, which is the solution from the complete model of the Solar System that matches best with the Newark-Hartford data, as the benchmark for our secular statistics. There are several criteria for determining how good a solution is (i.e., how well it could match with the geological data). A simple and rather straightforward criterion that we used is δ

= 1 4 4 i=1 (g i -g * i ) 2
, where g i and g * i are the frequencies from the astronomical solution and the geological data, respectively. A better-suited solution will have smaller δ and vice versa. We found that in the range from -200 Myr to -220 Myr, out of the 120 000 solutions, there are around 5000 (roughly 4.2% of the total number) that have smaller s than that of La2010d at -210 Myr, which is the value originally used to compare with the geological data. It should be noted that La2010d is one of 13 available complete solutions. The 95% CI of the probability of obtaining such a good matching solution from the complete model of the Solar System is thus (1.37%, 33.31%) [START_REF] Wilson | Probable inference, the law of succession, and statistical inference[END_REF]. Therefore, with the criterion δ, our result is statistically compatible with that of [START_REF] Olsen | Mapping Solar System chaos with the Geological Orrery[END_REF]. Laskar (1990) and [START_REF] Laskar | Confirmation of resonant structure in the solar system[END_REF] presented several secular resonances to explain the origin of chaos in the Solar System. In particular, the argument of the resonance (s 4 -s 3 ) -2(g 4 -g 3 ) is currently in a librational state, that is, (s 4 -s 3 ) -2(g 4 -g 3 ) = 0, (3.21) and moves out to the rotational state around -50 Myr. The dynamics can even switch to the librational state of a new resonance: This change was also visible in the La2004 astronomical solution, and the long cycle was attributed to the frequency combination g 4 -g 3 , which is visible from the spectrum of the eccentricity of the Earth. Although the exact value before and especially after the transition is not clear, the change in the period is visible from the band power of the core (Fig. 1 of [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF].

Libsack core

(s 4 -s 3 ) -(g 4 -g 3 ) = 0. ( 3
This change in g 4 -g 3 observed in the Libsack core corresponds to a transition from the resonance (s 4 -s 3 ) -2(g 4 -g 3 ), which is the resonance that the Solar System is currently at, to the resonance (s 4 -s 3 ) -(g 4 -g 3 ). With a large number of astronomical solutions, we could better understand this phenomenon. The transition is usually very fast: The frequency changes quickly to another value and then reverts back just as quickly. Therefore, to study the transition we used a smaller window for the FA, 10 Myr instead of 20 Myr. Figure 3.15 shows the density estimation of g 4 -g 3 at 90 Myr as well as at 300 Myr in the past. Both have a principle population in the range [0.4, 0.8] ′′ /yr and a small but not insignificant one centered around 1.0 ′′ /yr, which corresponds to the small chunk of g 4 centered at 18.2 ′′ /yr. The transition observed is a sudden jump in frequency from the main population to the secondary one of g 4 , and therefore g 4 -g 3 as well.

The size of the secondary population is defined as the proportion of the solutions whose g 4 -g 3 > 0.9 "/yr and is denoted as P(g 4 -g 3 > 0.9 ′′ /yr). The rate of transition is defined as the proportion of the solutions whose g 4 -g 3 > 0.9 ′′ /yr over 10 Myr and is denoted as P 10M yr (g 4 -g 3 > 0.9 "/yr). Both are shown in Fig. 3.16. During the predictable period, that is, from now until -50 Myr, no transition is observed. After -50 Myr, a transition can occur; its rate rises until 100 Myr, when the percentage of a secondary population stays relatively stable at around 2.1% ± 0.5% at a time, and the rate of the transition could also be determined to be 8% ± 1% every 10 Myr during this period. At 80-90 Myr, when the transition was detected in the Libsack core, that rate of transition during this period is found with our numerical solutions to be 7.7% ± 4.5%.

Conclusion

In this work we give a statistical description of the evolution of the fundamental frequencies of the Solar System beyond 60 Myr, that is, beyond the predictability horizon of the planetary motion, with the aim to quantify the uncertainty induced by its chaotic behavior. The base of our analysis is 120 000 orbital solutions of the secular model of the Solar System. The PDF of the frequencies is estimated by the KDE, and its uncertainty is evaluated by the MBB; both methods are tested via numerical experiments.

We benchmarked the secular model by sampling the initial conditions in different ways and then compared the density estimation of their solutions with one another and finally with the complete model. The results are twofold. First, regardless of how initial conditions were sampled, their density estimation will converge toward a single PDF; after this overlap, a robust estimation is guaranteed. Secondly, the density estimation of the secular model is compatible with that of the complete model of the Solar System. This agreement means that the results of the secular model, with superior computational simplicity, can be used for application to geological data.

We observe that the density estimations of the fundamental frequencies can be well fitted by skew Gaussian mixture models. The time evolution of the parameters σ 2 0 , related to the frequency variances, follows power-law functions. Interestingly for the inner fundamental frequencies, the exponents of such power laws are all smaller than 1, which indicates that they undergo sub-diffusion processes.

We show several examples of how this result can be used for geological applications. First, the uncertainty of any astronomical frequency signal is fully quantified, so that, for example, a proper quantitative response can be given to the question of how stable the astronomical metronomes are. With this statistical framework, previous results from geological records beyond 60 Myr can also be interpreted with a more comprehensive approach. A more quantitative answer, not only about the possibility but also about the probability of the occurrence of an astronomical signal in geological data, can be made. Apart from these direct applications, a more systematic approach could make full use of the density estimation of frequencies. The method TimeOpt from Meyers and Malinverno, 2018, for example, shows that it is possible to combine the uncertainty from astronomical signals with geological records to derive an effective constraint for both. In fact, any similar Bayesian method could use the density estimation of frequencies as proper priors.

Chapter 4

Long-term instability of the inner Solar System: numerical experiments

This chapter is published in [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF].

Introduction

Even though the planet orbits in the inner Solar System (ISS) are chaotic with a Lyapunov time of about 5 million years [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF](Laskar, , 1990;;[START_REF] Sussman | Chaotic evolution of the solar system[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF], they are still statistically very stable over a timescale that is a thousand times longer. The probability of a Mercury eccentricity higher than 0.7 over the next 5 billion years, for example, is about 1% from direct integrations of the Solar System [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF][START_REF] Abbot | Rare event sampling improves mercury instability statistics[END_REF]. This percentage agrees with the statistics of a dynamical instability observed in secular models where the dynamics is averaged over the planet mean longitudes [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. The statistical stability of the ISS over the remaining lifetime of the Sun as a main sequence star is intriguing, if one considers that it represents an open system, as it is forced by the very regular motion of the outer planets (Laskar, 1990;[START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. No exactly conserved quantities, such as the energy or angular momentum, can bound a priori the chaotic wanderings of the system in its high-dimensional phase space.

The disproportion between the Lyapunov time and the destabilisation timescale of the ISS has been addressed by [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF], building on previous works by [START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF] and [START_REF] Boué | A simple model of the chaotic eccentricity of Mercury[END_REF]. [START_REF] Boué | A simple model of the chaotic eccentricity of Mercury[END_REF] consider the firstorder secular dynamics of a mass-less Mercury in the gravitational field of all the other planets, whose orbits are predetermined to a quasi-periodic form. They use a multipolar expansion of the Hamiltonian to show that very high Mercury eccentricities appear in the reduced phase space of the resonance g 1 -g 5 (involving the fundamental precession frequencies of the Mercury and Jupiter perihelia), which confirms the role of this harmonic in the destabilisation of the ISS [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Batygin | On the Dynamical Stability of the Solar System[END_REF][START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF]. [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF] eccentricities and inclinations of the planets, and study a few of its Fourier harmonics related to the fundamental frequencies g 1 , g 2 , g 5 , s 1 , and s 2 . Their simplified dynamics is however much more unstable than realistic models, the typical time for the destabilisation of Mercury orbit being around 1 Gyr [START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF]. Recently, Mogavero and Laskar (2021, ML21 from now on) have proposed the model of a forced secular ISS, in which the outer planets only are frozen to quasi-periodic orbits. With a numerical experiment over 100 Gyr, they estimate the dynamical half-life of Mercury at 40 Gyr, consistently with the small probability of an instability over 5 Gyr.

Here we employ the computer algebra software TRIP (Gastineau andLaskar, 2011, 2021) to perform truncation of the forced secular ISS at different degrees in eccentricities and inclinations. Through ensembles of 10 3 to 10 5 numerical integrations spanning 5 to 100 Gyr, we show how dynamical contributions usually deemed as unimportant, that is, high-degree terms of the Hamiltonian and non-resonant harmonics, strongly affect the probability of an instability over 5 Gyr.

Dynamical models

In the forced secular model of the ISS (detailed presentation in ML21), the orbits of the outer planets are predetermined to a quasi-periodic form, whose frequencies and amplitudes are inferred from frequency analysis [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency map analysis and quasiperiodic decompositions[END_REF]) of a comprehensive model of the Solar System (Laskar et al., 2011a). The secular gravitational interactions are considered at first order in planetary masses, which corresponds to Gauss' dynamics of Keplerian rings [START_REF] Gauss | Determinatio attractionis quam in punctum quodvis positionis datae exerceret [etc[END_REF], and the leading contribution of general relativity (GR) is included. We remark that the shift of the frequency g 1 due to GR is comparable to the half-width of the principal secular resonances [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF], and cannot be considered as a small correction.

With the aid of TRIP, the secular Hamiltonian H of the entire Solar System, at first order in planetary masses, can be expanded in series of the complex Poincaré variables of the planets, i.e. (x i , xi , y i , ȳi ) 8 i=1 [START_REF] Laskar | Stability of the planetary three-body problem[END_REF]. The planets are indexed in order of increasing semi-major axis, as usual. Truncation at total degree 2n results in a polynomial Hamiltonian H 2n . When the predetermined orbits of the outer planets (x i (t), y i (t)) 8 i=5 are substituted, one obtains the Hamiltonian of the forced ISS truncated at degree 2n, i.e. H 2n ((

x i , y i ) 4 i=1 , t) = H 2n ((x i , y i ) 4 i=1 , (x i = x i (t), y i = y i (t)) 8 i=5
). The non truncated Hamiltonian, formally H = H ∞ , represents Gauss' dynamics of the forced ISS.

At the lowest degree, H 2 describes an integrable forced Laplace-Lagrange dynamics. Its analytical solution can be obtained by a canonical transformation to the complex proper modes variables (u i , v i ) 4 i=1 , with corresponding action-angle variables (X i , χ

i ; Ψ i , ψ i ) such that (u i = √ X i E -jχ i ; v i = √ Ψ i E -jψ i ) 1 .
When expressed in these action-angle variables, the truncated Hamiltonian is a finite Fourier series:

H 2n (I, θ, t) = k,ℓ H k,ℓ 2n (I)E j(k•θ+ℓ•ωot) , H k,ℓ 2n = n p=1 H k,ℓ (2p) , ( 4.1)
where I = (X, Ψ) and θ = (χ, ψ) are the eight-dimensional vectors of the action and angle variables, respectively, t is the time, ω o = (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ) is the septuple of the constant fundamental frequencies of the outer orbits (Laskar, 1990), and (k, ℓ) ∈ Z 8 × Z 7 is the wave vector of a given harmonic. The amplitude of a harmonic H k,ℓ 2n consists of partial contributions H k,ℓ (2p) from terms of the same degree 2p ≤ 2n. To identify these partial contributions we define

F k,ℓ (2p) = H k,ℓ (2p) E j(k•θ+ℓ•ωot) + c.c. for k = 0, l = 0, F 0,0 (2p) = H 0,0 (2p) , (4.2)
where c.c. refers to the complex conjugate. The order of a harmonic is defined as the even integer (k, ℓ) 1 ≤ 2n, where • 1 denotes the 1-norm. Since the quasi-periodic form of the outer orbits contains harmonics of order higher than one, the dynamics of H 2n and H 2n are not exactly the same. Yet, the difference is unimportant for the results of this work, so we shall treat the two Hamiltonians as equivalent from now on.

Second order in planetary masses. To investigate the effect of the order of the secular averaging on the long-term statistics, we employ the autonomous polynomial equations of motion of [START_REF] Laskar | Accurate methods in general planetary theory[END_REF]Laskar ( , 1990) ) for the ensemble of the Solar System planets. These equations formally derive from a Hamiltonian of order two in masses and degree 6 in eccentricities and inclinations, and will be denoted as L 6 throughout the paper. In this work, we also implement a variant of this dynamics, in which the equations for the inner planets are truncated at total degree 3 in eccentricities and inclinations, while those of the outer planets are kept at degree 5 (Appendix B.1). This new model, denoted as L 4 , is meant as an analogue of H 4 at second order in masses.

Numerical Experiments

We systematically derive the equations of motion for the truncated Hamiltonians H 2n in TRIP. They are numerically integrated via an Adams PECE scheme of order 12, with a time step of 250 years. Typical integration times are given in ML21 (table 1). All the the leading contribution coming from the Hamiltonian terms of degree 6. The practical stability of H 4 over 5 Gyr is unexpected, since it still reproduces the chaotic dynamics of the ISS with the same long-term statistical distribution of the maximum Lyapunov exponent as in H 6 or Gauss' dynamics [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. It also shows the same destabilisation mechanism, that is, the activation of the resonance g 1 -g 5 . Previous works on the instability of Mercury orbit studied a simplified dynamics in which only a few Fourier harmonics of H 4 are considered [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF][START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF]. This latter model is actually considerably more unstable than the reference model LG09, in deep contrast with the practical stable dynamics of H 4 over 5 Gyr. The secular models of degree higher than 4 give predictions that generally agree with the N -body integrations LG09. While Gauss's dynamics tends to underestimate, and H 8 overestimates the instability probability, H 10 gives accurate predictions (this is probably a coincidence related to the behaviour of the series expansion). At times shorter than 5 Gyr and for a threshold of 0.7, the estimations of H 6 , H 10 , and L 6 generally agrees with the results of a refined method of rare event detection applied to direct integrations [START_REF] Abbot | Rare event sampling improves mercury instability statistics[END_REF] 3 : for L 6 , the probability of Mercury orbit having eccentricity larger than 0.7 in the next 2 Gyr is 0.025% with a 90% confidence interval (0.019%, 0.034%), while it is 0.2% with a 90% confidence interval (0.18%, 0.22%) over 3 Gyr. The Hamiltonians of degree 4 and 6 show a relatively pronounced decay of the probability from 0.7 to 0.9 (see also Fig. B.1). This interesting phenomenon is inherent to degree 6 or lower, because the probability of the models of higher degree is fairly constant across high values of Mercury's eccentricity: if e 1 reaches 0.7 along a solution, it also probably goes beyond 0.9.

Ranking of harmonics according to their contributions to g 1

In order to explain the difference between the statistics of H 4 and H 6 , we shall switch to the proper mode variables and the Fourier representation of Eq. (4.1). The Hamiltonian H 6 contains substantially more harmonics than H 4 , 69 339 compared to 2 748. For each harmonic of H 4 , H 6 includes additional terms of degree 6 in its amplitude. Despite the large difference in the number of terms, most of the contributions of H 6 are negligible. We aim to identify here the Fourier harmonics that have an important impact on the destabilisation mechanism, that is, the activation of the resonance g 1 -g 5 . Because g 5 is constant in the forced dynamics, we shall focus on the fundamental precession frequency of Mercury perihelion g 1 . Following ML21, the instantaneous value of the frequency g 1 for the Hamiltonian H 2n is defined as:

ĝ(2n) 1 = -θ1 = - ∂H 2n ∂I 1 = k,ℓ n p=1 ĝk,ℓ 1(2p) , ( 4.4)
where the partial contributions at degree 2p of each harmonic are ĝk,ℓ 4.5), (4.6)). For each partial degree, the harmonics are ranked according to |C k,ℓ (2p) | with T = 2Gyr, which is shortly after the first activation of the resonance g 1 -g 5 .

1(2p) = - ∂F k,ℓ (2p) (I) ∂I 1 = - ∂ H k,ℓ (2p) (I) ∂I 1 E j(k•θ+ℓ•ωot) + c.c. for k = 0, l = 0, ĝ0,0 1(2p) = - ∂F 0,0 (2p) (I) ∂I 1 = - ∂ H 0,0 (2p) (I) ∂I 1 . (4.5) i Harmonic F i (4) C k,ℓ (4) Harmonic F i (6) C k,ℓ (6) 1 0 -0.836 0 -0.614 2 s 1 -s 2 -0.790 s 1 -s 2 -0.573 3 2g 1 -s 1 -s 2 0.413 2g 1 -s 1 -s 2 0.268 4 2g 1 -2s 1 0.366 2g 1 -2s 1 0.237 5 g 1 -g 5 -0.126 g 1 -g 2 -0.071 6 2g 1 -2s 2 0.117 2g 1 -2s 2 0.054 7 g 1 -g 2 + s 1 -s 2 -0.054 g 1 -g 5 -0.053 8 g 1 -g 2 0.035 g 1 -g 2 + s 1 -s 2 -0.044 9 g 1 -g 2 -s 1 + s 2 -0.027 2s 1 -2s 2 0.044 10 s 1 -s 3 0.024 g 1 -g 2 -s 1 + s 2 -0.039
In this form, each harmonic manifests its importance via its direct contribution to g 1 , which varies along an orbital solution according to the position in the phase space, i.e. ĝk,ℓ 1(2p) (t) = ĝk,ℓ 1(2p) (I(t), θ(t), t). To identify the main harmonics involved in the destabilisation of the dynamics, Eqs. (4.4) and (4.5) are evaluated along unstable solutions. Short-term oscillations are suppressed by the low-pass Kolmogorov-Zurbenko (KZ) filter (Yang and Zurbenko, 2010), which is applied to the instantaneous frequency g 1 and its harmonic contributions. We use the KZ filter with 3 iterations of the moving average and a cutoff frequency of (1 Myr) -1 (ML21, appendix B) to obtain the filtered values g

(2n) 1 = KZ(ĝ (2n) 1 ), g k,ℓ 1(2p) = KZ(ĝ k,ℓ 1(2p) ) (4.6)
The harmonics can then be ranked according to the maximum value of their absolute filtered contribution over the time interval [0, T ]. The timespan T is chosen to be slightly larger than the time of the first activation4 of the resonance g 1 -g 5 . After this point, the system either exhibits a secular collision right away or enters a period of excited dynamics before an eventual collision. This unstable state typically lasts longer for a solution of H 6 than for a Hamiltonian of higher degree.

We establish the harmonic ranking on an unstable solution of H 6 , whose Mercury's eccentricity over time is shown in Fig. 4.2 (the ranking of the leading harmonics is quite robust when we switch to other unstable solutions). The first entrance into the chaotic zone of the resonance g 1 -g 5 occurs just after 1.97 Gyr (see Fig. 4.3a), during which the eccentricity of Mercury is pumped to 0.65 and the harmonic contributions generally reach their maximum values (see Fig. 4.3b). The ranking is computed over the first 2 Gyr to capture the contributions of the harmonics at the resonance. Table 4.2 shows two harmonic rankings based on the partial contributions at degree 4 and 6, respectively. It is surprising to find that the contributions to g 1 at degree 6 are slightly less, but still roughly the same amount as those at degree 4. Because the principal contributions at degree 6 come from harmonics of order 2 and 4, what H 6 mainly offers is not new resonances, but rather corrections to the existing harmonics of H 4 . The corrections at degree 6 help to push g 1 toward g 5 and bring the solution closer to the destabilizing resonance. Geometrically speaking, in the phase space the resonance g 1 -g 5 defined by H 6 is closer to the current ISS than that of H 4 .

Figure 4.3b gives a closer look at the time evolution of the leading harmonic contributions to g 1 at degree 6. They are small at the beginning when the solution is stable, but get much bigger when the eccentricity of Mercury becomes higher, that is during and after the first activation of the resonance g 1 -g 5 at 1.97 Gyr. During this period, which is shown in the lower panel of Fig. 4.3b, the strongest terms are the null-frequency harmonic, i.e. the integrable part of the Hamiltonian H 0,0 (6) , and the harmonic s 1 -s 2 , which also enters resonance. These two terms tend to destabilize the system by decreasing g 1 by substantial amounts, which are even greater than the leading GR correction of 0.4 ′′ yr -1 at degree 2 at some point. In the opposite direction, the two harmonics 2g 1 -(s 1 + s 2 ) and 2(g 1 -s 1 ) raise g 1 , moving it away from g 5 . Although these terms are non resonant, they are extremely crucial for the stability of Mercury orbit (see Section 4.4.3). Other harmonics also contribute to g 1 at degree 6 in an alternating pattern, but to a lesser extent.

To confirm the crucial role of the terms of degree 6, we add them to H 4 to construct partial Hamiltonians [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]:

H 4,m = H 4 + m i=1 F i (6) , (4.7)
where

F i (6) = F k i ,ℓ i (6)
is the ith harmonic from the ranking at degree 6 of Table 4.2, and m is the total number of such harmonics that are considered. Figure 4.3a shows the filtered g 1 computed from different Hamiltonians along the same unstable trajectory of H 6 of Fig. 4.2. Initially, when the solution is stable and Mercury's eccentricity is relatively low, the frequency g (4)

1 of H 4 is almost indistinguishable from the corresponding g (6)

1 of H 6 . Across the activation of the resonance g 1 -g 5 , the difference between the two frequencies becomes considerable: g (6)

1 almost reaches g 5 , while g (4)

1 does not. The difference is mainly due to the integrable term H 0,0 (6) , which is included in H 4,1 , and to the first leading harmonics contained in H 4,4 .

The statistics of the high Mercury eccentricities from H 4,m should approximate that of H 6 better than H 4 . In order to test this expectation, we integrate the dynamics of H 4,1 , H 4,4 and H 4,51 from 10 800 and 1 080 initial conditions over 5 Gyr and 100 Gyr, respectively. The initial conditions are taken from the same ensembles employed for H 2n . The CDFs of the first time that Mercury eccentricity reaches 0.7 for H 4,m are shown in Fig. 4.1. The wide discrepancy between H 4 and H 6 is first bridged by adding the integrable term H 0,0 (6) , with which the curve of H 4,1 attains a probability of 0.2% at 5 Gyr. Including the next three leading harmonics brings the curve to the same level as Gauss' dynamics. Adding additional terms makes the statistics oscillate around that of H 6 .

The impact of the choice of the initial conditions on the present analysis deserves a discussion. As stated in Section 4.3, the nominal initial conditions of the truncated forced dynamics H 2n are chosen to be the same as those of Gauss' dynamics H. In principle, they should be adapted to each model according to the harmonics that are dropped from the to be the same as in Section 4.4.2, that is, m ∈ {1, 4, 51}. We integrate the equations of motion defined by H * 4 and H * 4,m over 5 Gyr from the same ensemble of initial conditions defined in Section 4.3, to obtain 10 800 solutions.

Figure 4.4 shows the comparison between H 4 , H 4,m and H * 4 , H * 4,m for the CDF of the first time that Mercury's eccentricity reaches 0.7 over 5 Gyr. For all the models, removing the three non-resonant harmonics makes the dynamics significantly more unstable, with at least one order of magnitude of difference. For comparison, the dynamics of H 4 is a thousand times more stable than H 6 over 5 Gyr, but taking away the three harmonics brings the model H * 4 basically to the same level of instability of H 6 . If we consider the Hamiltonian H * 4,51 , based on H 4,51 which is the closest dynamics to H 6 among the presented partial Hamiltonians, its probability of instability is around 20% at 5 Gyr, that is, twenty times more than the instability rate of H 6 . These numerical experiments show the crucial role of these non-resonant harmonics in stabilizing the ISS. Interestingly enough, all the three harmonics permits the exchange of angular momentum deficit between the eccentricity and inclination degrees of freedom, that is, between the proper modes (u i ) and (v i ). These results also show the sensitivity of the destabilisation probability to the details of the dynamics, and may explain, at least partially, the great instability shown by the simplified models considered in literature [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF][START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF].

Statistics over 100 Gyr

To explore the dynamics in a regime where highly excited orbits no longer represent rare events, we follow ML21 and prolong 1080 integrations of the different Hamiltonians previously considered to 100 Gyr. Figure 4.5 shows the CDFs of the first time that Mercury eccentricity reaches 0.7 and the corresponding probability density functions (PDFs). The PDFs are estimated by the debiased kernel density estimation (KDE) method [START_REF] Cheng | Nonparametric inference via bootstrapping the debiased estimator[END_REF], with Gaussian kernel and Silverman (1986)'s rule-of-thumb bandwidth (Appendix B.4). We use the log transformation and the pseudo-data method [START_REF] Cowling | On pseudodata methods for removing boundary effects in kernel density estimation[END_REF] to remove the boundary effects induced by the KDE at 0 and 100 Gyr, respectively. The confidence intervals of the PDFs are estimated by bootstrap [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF] of the debiased KDEs; for the CDFs, we use Wilson's score interval.

The CDFs of H 6 and H 8 are close to each other, with medians of 30 Gyr, while that of H 10 is around 35 Gyr. The increasing values of the medians may suggest a convergence toward the value of 40 Gyr of Gauss' dynamics. On the other hand, the difference between H 4 and the other truncated forced dynamics is still considerable. The median time for H 4 is 75 Gyr, roughly doubling the value of H 6 . If we assume that the PDFs follow a Levy distribution ρ(τ ) = (T 0 /πτ 3 ) 1/2 E -T 0 /τ over short times (ML21), with T 0 proportional to the median of the distribution, one easily understands how a difference by a factor of two in the medians of the PDFs results in very different probabilities over 5 Gyr. Indeed, the fact that the destabilisation over 5 Gyr is a rare event greatly amplifies the disparity between H 4 and the models of higher degree.

There is practically no difference between the statistics of H 6 and L 6 over this timescale, which confirms the secondary effect of the second order in masses for the forced ISS and the statistics of the high Mercury eccentricities in particular. However, this effect is magnified for H 4 , the CDF of L 4 approaching halfway the curves of higher degrees, with of its action variables [START_REF] Hénon | Integrals of the Toda lattice[END_REF][START_REF] Flaschka | The Toda lattice. II. Existence of integrals[END_REF][START_REF] Ferguson | Nonlinear Normal Modes for the Toda Chain[END_REF][START_REF] Benettin | The Fermi-Pasta-Ulam Problem and Its Underlying Integrable Dynamics[END_REF]. Although not integrable, and indeed chaotic, the Hamiltonian H 4 plays a role similar to the Toda Hamiltonian, as it does not allow essentially any dynamical instability over 5 Gyr. The main question at this point is why the dynamics of H 4 is practically stable over 5 Gyr. Once this is assessed, the small 1% probability of an instability of the ISS may be conceived as a natural perturbative effect of terms of degree 6 and higher.

Chapter 5

Timescales of chaos in the Inner Solar System

This chapter was submitted and is currently under review for publication.

Introduction

The planetary orbits in the inner Solar System (ISS) are chaotic, with a Lyapunov time distributed around 5 million years (Myr) [START_REF] Laskar | A numerical experiment on the chaotic behaviour of the solar system[END_REF](Laskar, , 1990;;[START_REF] Sussman | Chaotic evolution of the solar system[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. 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) [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF][START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. The dynamical half-life of Mercury orbit has recently been estimated at 30 to 40 billion years [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF][START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. A disparity of nearly four orders of magnitude between the Lyapunov time and the timescale of dynamical instability is intriguing, since the chaotic variations of the orbits of the inner planets cannot be constrained a priori. While the total energy and angular momentum of the Solar System are conserved, the disproportion of masses between the outer and inner planets implies that unstable states of the ISS are in principle easily realizable through exchanges of these quantities. The surprising stability of the ISS deserves a global picture in which it can emerge more naturally.

To our knowledge, the only study addressing the timescale separation in the long-term dynamics of the ISS is based on the simplified secular dynamics of a mass-less Mercury [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF]: All the other planets are frozen on regular quasi-periodic orbits; secular interactions are expanded to first order in masses and degree 4 in eccentricities and inclinations; an a priori choice of the relevant terms of the Hamiltonian is made. The typical instability time of about 1 Gyr [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF][START_REF] Woillez | Instantons for the destabilization of the inner Solar System[END_REF] 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 [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. We have indeed 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 system1 . From the perspective of these latest findings, the small probability of 1% of an instability over the age of the Solar System may be naturally regarded as a perturbative effect of terms of degree 6 and higher. Clearly, the striking stability of the dynamics at degree 4 is even more impressive in the present context, and remains to be explained.

A strong separation in dynamical timescales is not uncommon among classical quasiintegrable systems (e.g. [START_REF] Milani | An example of stable chaos in the Solar System[END_REF][START_REF] Morbidelli | On the Relationship Between Lyapunov Times and Macroscopic Instability Times[END_REF]. This is notably evinced by the Fermi-Pasta-Ulam-Tsingou (FPUT) problem, which deals with a chain of coupled weakly-anharmonic oscillators [START_REF] Fermi | Studies of Nonlinear Problems[END_REF]. Far from KAM and Nekhoroshev regimes (as is likely to be pertinent to the ISS, see Sect. 5.3), 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 [START_REF] Lam | Stochastic Perturbation of Integrable Systems: A Window to Weakly Chaotic Systems[END_REF][START_REF] Goldfriend | Equilibration of quasi-integrable systems[END_REF]. In other words, the Lyapunov time and the diffusion/instability time scale differently with the size of the terms that break integrability, and this can result in very different timescales (e.g. [START_REF] Morbidelli | On the Relationship Between Lyapunov Times and Macroscopic Instability Times[END_REF]. However, this argument is as general as poorly satisfactory in addressing quantitatively the timescale separation in a complex problem as the present one. Moreover, even though order-of-magnitude estimates of the chaotic diffusion in the ISS suggest that it may take hundreds of million years to reach the destabilizing secular resonance g 1 -g 5 2 , the low probability of an instability over 5 Gyr still remains unexplained [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. 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 [START_REF] Laskar | Large Scale Chaos and Marginal Stability in the Solar System[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF] 3 . With its 8 secular degrees of freedom (DOFs), this system also constitutes a peculiar bridge between the low-dimensional dynamics often addressed in celestial mechanics and the systems with a large number of bodies studied in statistical mechanics: It cannot benefit from the straightforward application of standard methods of the two fields (e.g. Mogavero and Laskar, 2022, Appendix A).

This work aims to open a window on the long-term statistical behavior of the inner planet orbits. Section 5.2 briefly recalls the dynamical model of forced secular ISS introduced in [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. Section 5.3 presents the numerical computation of its Lyapunov spectrum. Section 5.4 introduces the quasi-symmetries of the resonant harmonics of the Hamiltonian and the corresponding quasi-integrals of motion. Section 5.5 establishes a geometric connection between the quasi-integrals and the slowest DOFs of the dynamics via a principal component analysis of the orbital solutions. Section 5.6 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.

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 (Laskar, 1990;[START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF]. At first order in planetary masses, the secular Hamiltonian corrected for the leading contribution of general relativity reads (e.g. [START_REF] Morbidelli | Modern celestial mechanics: aspects of solar system dynamics[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF])

H = - 8 i=1   i-1 l=1 Gm i m l r i -r l + 3G 2 m 2 0 m i c 2 a 2 i 1 -e 2 i   .
(5.1)

The planets are indexed in order of increasing semi-major axes (a i ) 8 i=1 , m 0 and m i are the Sun and planet masses, respectively, e i the eccentricities, G the gravitational constant and c the speed of light. The vectors r i are the heliocentric positions of the planets, and the bracket operator represents the averaging over the mean longitudes resulting from the elimination of the non-resonant Fourier harmonics of the N -body Hamiltonian [START_REF] Morbidelli | Modern celestial mechanics: aspects of solar system dynamics[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. Hamiltonian (5.1) generates Gauss's dynamics of Keplerian rings [START_REF] Gauss | Determinatio attractionis quam in punctum quodvis positionis datae exerceret [etc[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF], whose semi-major axes a i are constants of motion of the secular dynamics.

By developing the 2-body perturbing function (Laskar, 1991;[START_REF] Laskar | Stability of the planetary three-body problem[END_REF] in the computer algebra system TRIP (Gastineau andLaskar, 2011, 2021), the secular Hamiltonian can be systematically expanded in series of the Poincaré rectangular coordinates in complex form,

x i = Λ i 1 -1 -e 2 i E j̟ i , y i = 2Λ i 1 -e 2 i 1 4 sin(I i /2)E jΩ i , (5.2) where Λ i = µ i [G(m 0 + m i )a i ] 1/2 , µ i = m 0 m i /(m 0 + m i )
being the reduced masses of the planets, I i the inclinations, ̟ i the longitudes of the perihelia and Ω i the longitudes of the nodes 4 . Pairs (x i , -jx i ) and (y i , -jy i ) are canonically conjugate momentum-coordinate variables. When truncating at a given total degree 2n in eccentricities and inclinations, the expansion provides Hamiltonians

H 2n = H 2n [(x i , xi , y i , ȳi ) 8 i=1 ] that are multivariate polynomials.
Valuable insight into the dynamics of the inner planets is provided by the model of a forced ISS recently proposed [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. It exploits the great regularity of the long-term motion of the outer planets (Laskar, 1990;[START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF][START_REF] Hoang | Chaotic diffusion of the fundamental frequencies in the Solar System[END_REF] to predetermine their orbits to a quasi-periodic form: (5.3) for i ∈ {5, 6, 7, 8}, where t denotes time, xil and ỹil are complex amplitudes, m il and n il integer vectors, and ω o = (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ) represents the septuple of the constant fundamental frequencies of the outer orbits. Frequencies and amplitudes of this Fourier decomposition are established numerically by frequency analysis [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency map analysis and quasiperiodic decompositions[END_REF] of a comprehensive orbital solution of the Solar System (Mogavero and Laskar, 2021, Appendix D). Gauss's dynamics of the forced ISS is obtained by substituting the predetermined time dependence in Eq. (5.1),

x i (t) = M i l=1 xil E jm il •ωot , y i (t) = N i l=1 ỹil E jn il •ωot ,
H = H[(x i , y i ) 4 i=1 , (x i = x i (t), y i = y i (t)) 8 i=5 ],
(5.4)

so that H = H[(x i , y i ) 4 i=1 , t].
The resulting dynamics consists of two DOFs for each inner planet, corresponding to the x i and y i variables, respectively. Therefore, the forced secular ISS is described by 8 DOFs and an explicit time dependence. As a result of the forcing from the outer planets, no trivial integrals of motion exist and its orbital solutions live in a 16-dimensional phase space.

A truncated Hamiltonian H 2n for the forced ISS is readily obtained by substituting Eq. ( 5.3) in the truncated Hamiltonian H 2n of the entire Solar System. At the lowest degree, H 2 generates a linear, forced Laplace-Lagrange (LL) dynamics. This can be analytically integrated by introducing complex proper mode variables (u i , v i ) 4

i=1 via a time-dependent canonical transformation (x i , -jx i ) → (u i , -ju i ), (y i , -jy i ) → (v i , -jv i ) [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. Action-angle pairs (X i , χ i ), (Ψ i , ψ i ) are introduced as

u i = X i E -jχ i , v i = Ψ i E -jψ i .
(5.5)

When expressed in the proper modes, the truncated Hamiltonian can be expanded as a finite Fourier series: (5.6) where I = (X, Ψ) and θ = (χ, ψ) are the 8-dimensional vectors of the action and angle variables, respectively, and we introduced the external angles φ(t) = -ω o t. The wave vectors (k, ℓ) belong to a finite subset of Z 8 × Z 7 . At degree two, one has

H 2n (I, θ, t) = k,ℓ H k,ℓ 2n (I)E j(k•θ+ℓ•φ(t)) ,
H 2 = -ω LL • I,
where ω LL = (g LL , s LL ) ∈ R 4 × R 4 are the LL fundamental precession frequencies of the inner planet perihelia and nodes. Hamiltonian H 2n is in quasi-integrable form.

The quasi-periodic form of the outer orbits in Eq. ( 5.3) contains harmonics of order higher than one, that is, m il 1 > 1 and n il 1 > 1 for some i and l, where • 1 denotes the 1-norm. Therefore, the dynamics of H 2n and H 2n are not exactly the same [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. Still, the difference is irrelevant for the results of this work, so we shall treat the two Hamiltonians as equivalent from now on. Despite the simplifications behind Eqs. (5.1) and (5.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 (Laskar, 1990;[START_REF] Laskar | A long-term numerical solution for the insolation quantities of the Earth[END_REF][START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF]. It correctly reproduces the finite-time maximum Lyapunov exponent (FT-MLE) and the statistics of the high eccentricities of Mercury over 5 Gyr [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF].

Lyapunov spectrum

Ergodic theory provides a way, through the Lyapunov characteristic exponents (LCEs), to introduce a fundamental set of timescales for any differentiable dynamical system ż = F (z, t) defined on a phase space P ⊆ R P [START_REF] Oseledec | A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems[END_REF][START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF][START_REF] Gaspard | Chaos, Scattering and Statistical Mechanics[END_REF][START_REF] Skokos | The Lyapunov Characteristic Exponents and Their Computation[END_REF]. If Φ(z, t) denotes the associated flow and z(t) = Φ(z 0 , t) the orbit that emanates from the initial condition z 0 , the LCEs λ 1 ≥ λ 2 ≥ • • • ≥ λ P are the logarithms of the eigenvalues of the matrix Λ(z 0 ) defined as lim t→∞ M(z 0 , t) T M(z 0 , t)

1/2t = Λ(z 0 ), (5.7)

where M(z 0 , t) = ∂Φ/∂z 0 is the fundamental matrix and T stands for transposition [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF][START_REF] Gaspard | Chaos, Scattering and Statistical Mechanics[END_REF] 0). The multiplicative ergodic theorem of [START_REF] Oseledec | A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems[END_REF] states that if ρ is an ergodic (i.e. invariant and indecomposable) measure for the time evolution and has compact support, then the limit in Eq. (5.7) exists for ρ-almost all z 0 , and the LCEs are ρ-almost everywhere constant and only depend on ρ [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]. Moreover, one has

lim t→∞ 1 t log M(z 0 , t)ζ 0 = λ (i) if ζ 0 ∈ E (i) z 0 \ E (i+1) z 0 , (5.8)
for ρ-almost all z 0 , where λ (1) > λ (2) > . . . are the LCEs without repetition by multiplicity, and E (i) z 0 is the subspace of R P corresponding to the eigenvalues of Λ(z 0 ) that are smaller than or equal to exp λ (i) , with T z 0 P = E (1) [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF][START_REF] Skokos | The Lyapunov Characteristic Exponents and Their Computation[END_REF]. Once the LCEs have been introduced, a characteristic timescale can be defined from each positive exponent as λ -1 i . In the case of the maximum Lyapunov exponent, λ 1 , the corresponding timescale is commonly called the Lyapunov time.

z 0 ⊃ E (2) z 0 ⊃ • • • . The specific choice of the R P -norm • in Eq. (5.8) is irrelevant
For a Hamiltonian system with p DOFs (i.e. P = 2p), the fundamental matrix is symplectic and the set of LCEs is symmetric with respect to zero, that is,

∆λ i := λ i + λ 2p-i+1 = 0 for all 1 ≤ i ≤ p.
(5.9)

If the Hamiltonian is time independent, a pair of exponents vanishes. In general, the existence of an integral of motion C = C(z) implies a pair of null exponents, one of them being associated with the direction of the tangent space that is normal to the surface of constant C (e.g. [START_REF] Gaspard | Chaos, Scattering and Statistical Mechanics[END_REF].

The ISS is a clear example of a dynamical system that is out of equilibrium. Its phasespace density diffuses seamlessly over any meaningful timescale [START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Hoang | Chaotic diffusion of the fundamental frequencies in the Solar System[END_REF]. Therefore, the infinite time limit in Eq. (5.7) is not physically relevant. The non-null probability of a collisional evolution of the inner planets [START_REF] Laskar | Large-scale chaos in the solar system[END_REF][START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Batygin | On the Dynamical Stability of the Solar System[END_REF][START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF] 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 [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF][START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. Physically relevant quantities are the finite-time LCEs (FT-LCEs), λ i (z 0 , t), defined from the eigenvalues m 1 ≥ m 2 ≥ • • • ≥ m P of the time-dependent symmetric positive-defined matrix M(z 0 , t) T M(z 0 , t) as

λ i (z 0 , t) = 1 2t log m i (z 0 , t).
(5.10)

On the top of that, a computation of the exponents for an ensemble of trajectories is advisable for a non-ergodic dynamics [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. These considerations show how demanding the computation of the Lyapunov spectrum of the Solar System planets is. Nevertheless, a 5-Gyr integration of the forced ISS takes only a couple of hours for Gauss's dynamics (H) and a few minutes at degree 4 (H 4 ). This dynamical model thus provides an 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 [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF], based on Gram-Schmidt orthogonalization. Manipulation of the truncated Hamiltonian H 2n in TRIP allows us to systematically derive the equations of motion and the corresponding variational equations, which we integrate through an Adams PECE method of order 12 and a timestep of 250 years. Parallelization of the time evolution of the 16 tangent vectors, between two consecutive reorthonormalization steps of the [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF] algorithm, significantly reduces the computation time. Figure 5.1a shows the positive FT-LCEs expressed as angular frequencies over the next 10 Gyr for the Hamiltonian truncated at degree 4. The FT-LCEs are computed for 150 stable solutions, with initial conditions very close to the nominal values of Gauss's dynamics and random sets of initial tangent vectors (Mogavero and Laskar, 2022, 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 (Appendix C. [START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF]Fig. C.1a). In this asymptotic regime, the [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF] algorithm purely retrieves the FT-LCEs as defined in Eq. (5.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. (5.9) (Appendix C. [START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF]Fig. C.1b).

The spectrum in Fig. 5.1a has distinctive features. A set of intermediate exponents follow the MLE, ranging from 0.1 to 0.01 ′′ yr -1 , while the smallest ones fall below 0.01 ′′ yr -1 . Figure 5.1a reveals the existence of a hierarchy of exponents and corresponding timescales that spans two orders of magnitude, down to a median value of λ -1 8 ≈ 500 Myr. The number of positive exponents confirms that no integral of motion exists, as one may expect from the forcing of the outer planets. We also compute the spectrum for the Hamiltonian truncated at degree 6. As shown in Appendix C. 1 (Fig. C.2), the asymptotic distributions of the exponents are very similar to those at degree 4. This result suggests that long-term diffusion of the phase-space density is very close in the two cases. The different instability rates of the two truncated dynamics mainly relates to the geometry of the instability boundary, which is closer to the initial position of the system for H 6 than for H 4 [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF].

The relevance of the Lyapunov spectrum in Fig. 5.1a emerges from the fact that the existence of an integral of motion implies a pair of vanishing exponents. This is a pivotal point: By a continuity argument, the presence of positive exponents much smaller than the leading one constitutes a compelling indication that there are dynamical quantities whose chaotic decoherence over initially very close trajectories takes place over timescales much longer than the Lyapunov time. In the long term, such quantities should diffuse much more slowly than any LL action variable. Therefore, Fig. 5.1a 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 towards the unstable states. Next section addresses the emergence of these slow quantities from the symmetries of the Fourier harmonics that compose the Hamiltonian.

Quasi-integrals of motion

The emergence of a chaotic behavior of the planetary orbits can be explained in terms of the pendulum-like dynamics generated by each Fourier harmonic that composes the Hamiltonian in Eq. (5.6) (e.g. [START_REF] Chirikov | A universal instability of many-dimensional oscillator systems[END_REF]. One can write H 2n (I, θ, t) = H 0,0 2n (I) + Mn i=1 F i (I, θ, t), with (5.11) where (k i , ℓ i ) = (0, 0), M n is the number of harmonics in 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 F i whose frequency combination

F i (I, θ, t) = H k i ,ℓ i 2n (I) E j(k i •θ+ℓ i •φ(t)) + c.c.,
k i • θ + ℓ i • φ(t)
vanishes at some point along the motion. Using the computer algebra system TRIP, the harmonics of H 10 that enter into resonance along the 5-Gyr nominal solution of Gauss's dynamics have been systematically retrieved, together with the corresponding time statistic of the resonance half-widths ∆ω [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. The resonances have then been ordered by decreasing time median of their half-widths. The resulting ranking of resonances will be denoted as R 1 from now on. Table 5.1 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 5.1 also shows the order of each harmonic, defined as the even integer O = (k, ℓ) 1 . The support of the asymptotic ensemble distribution of the FT-MLE shown in Fig. 5.1a overlaps in a robust way with that of the time distribution of the half-width of the strongest resonances. In other words, 2πλ 1 ≈ ∆ω R 1 , (5.12)

(g i , s i ) 8 i=1 , that is, k • ω i + ℓ • ω o , with ω i = (g 1 , g 2 , g 3 , g 4 , s 1 , s 2 , s 3 , s 4 ). Table
where ∆ω R 1 stands for the half-width of the uppermost resonances of ranking R 1 . Equation (5.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 assumption that the time distribution of physical observables along the 5-Gyr nominal solution of Gauss's dynamics should be representative of their ensemble distribution (defined by a set of stable orbital solutions with very close initial conditions) at some large time of the order of billions of years [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. We remark that, strictly speaking, ranking R 1 is established on the Fourier harmonics of the Lie-transformed Hamiltonian H ′ 2n (Mogavero and Laskar, 2022, Appendix G). New canonical variables are indeed defined to transform H 2n in a Birkhoff normal form to Table 5.1 -Top of ranking R 1 . First 30 resonant harmonics of H 10 along the 5-Gyr nominal solution of Gauss's dynamics, in order of decreasing time median of the resonance half-width ∆ω (arcsec yr -1 ). O is the order of the harmonic. τ res is the fraction of time the harmonic is resonant. Only harmonics with τ res > 1% are shown. 5th and 95th percentiles of the time distribution of ∆ω as subscripts and superscripts, respectively. Adapted from [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF] degree 4. The goal is to let the interactions of the terms of degree 4 in H 2n to appear more explicitly in the amplitudes of the harmonics of higher degrees in H ′ 2n , the physical motivation being that the non-linear interaction of the harmonics at degree 4 constitutes the primary source of chaos [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. Keeping in mind the quasiidentity nature of the Lie transform, here we shall drop for simplicity the difference between the two Hamiltonians. Moreover, all the new analyses of this work will involve the original variables of Eq. (5.5).

Quasi-symmetries of the resonant harmonics

In addition to the dynamical interactions responsible for the chaotic behavior of the orbits, Table 5.1 provides information on the geometry of the dynamics in the action variable space. Ranking the Fourier harmonics allows to consider partial Hamiltonians constructed from a limited number m of leading terms [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF][START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF], that is,

H 2n,m = H 0,0 2n + m i=1 F i .
(5.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 shall be interested in how these symmetries are progressively destroyed when one increases the number of terms taken into account in Eq. (5.13).

Consider a set of m harmonics of H 2n and a dynamical quantity that is a linear combination of the action variables, that is,

C γ = γ • I,
(5.14) γ ∈ R 8 being a parameter vector. From Eq. (5.11), the partial contribution of the m harmonics to the time derivative of C γ along the flow of

H 2n is Ċγ,m = 2 m i=1 γ • k i Im{ H k i ,ℓ i 2n (I) E j(k i •θ+ℓ i •φ(t)) }, (5.15)
and Ċγ = Ċγ,Mn , where M n is the total number of harmonics with a non-null wave vector that appear in H 2n . Any quantity C γ with γ • k i = 0 is conserved by the one-DOF dynamics generated by the single harmonic F i . In other words, such a quantity would be an integral of motion if F i were the only harmonic to appear in the Hamiltonian.

Considering now m different harmonics, these do not contribute to the change of the quantity

C γ if γ ⊥ span(k 1 , k 2 , . . . , k m )
, that is, if the vector γ belongs to the orthogonal complement of the linear subspace of R 8 spanned by the wave vectors (k i ) m i=1 . We shall also consider the quantity

C ′ γ = H 2n + γ • I.
(5.16)

Due to the explicit time dependence in the Hamiltonian, the partial contribution of a set of m harmonics to the time derivative of

C ′ γ along the flow of H 2n is Ċ′ γ,m = 2 m i=1 (γ • k i + ℓ i • ω o ) Im{ H k i ,ℓ i 2n (I) E j(k i •θ+ℓ i •φ(t)) }, (5.17)
the harmonics in Table 5.1 verify the relation 4 i=1 k inc i = 0, where k inc = (k inc 1 , . . . , k inc 4 ). Therefore, denoting γ 1 = (0 4 , 1 4 ), the quantity

C inc := C γ 1 = Ψ 1 + Ψ 2 + Ψ 3 + Ψ 4 (5.19)
is conserved by these resonances. C inc is the angular momentum deficit (AMD) [START_REF] Laskar | Large scale chaos and the spacing of the inner planets[END_REF] contained in the inclination DOFs. This symmetry can be then 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 DOFs, C ecc = 4 i=1 X i , is not invariant under the leading resonances because of the eccentricity forcing mainly exerted by Jupiter through the mode g 5 . The conservation of C inc depends on two facts: the inclination modes s 6 , s 7 , s 8 of the external forcing do not appear in Table 5.1; low-order harmonics like 2g 1 -s 1 -s 2 , 2g 1 -2s 1 and 2g 1 -2s 2 are never resonant (even if they can raise large quasi-periodic contributions), so that two AMD reservoirs C ecc and C inc are decoupled in Table 5.1. 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 (Laskar, 1990;[START_REF] Morbidelli | Modern celestial mechanics: aspects of solar system dynamics[END_REF][START_REF] Boué | Spin axis evolution of two interacting bodies[END_REF].

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 5.1. Consider the 30 × 8 matrix K 30 whose rows are the wave vectors (k i ) 30 i=1 of the listed resonances. A singular value decomposition shows that the rank of K 30 is equal to 6. Therefore, the linear subspace V 30 = span(k 1 , k 2 , . . . , k 30 ) spanned by the wave vectors has dimension 6. A Gram-Schmidt orthogonalization allows to determine two linearly independent vectors that span its orthogonal complement V ⊥ 30 . One choice consists in V ⊥ 30 = span(γ 2 , γ ⊥ 2 ), with

C 2 := C γ 2 = -X 3 -X 4 + Ψ 1 + Ψ 2 + 2Ψ 3 + 2Ψ 4 , C ⊥ 2 := C γ ⊥ 2 = X 3 + X 4 + Ψ 1 + Ψ 2 .
(5.20)

Since the second symmetry clearly requires that γ 1 ∈ V ⊥ 30 , the three quantities

C inc , C 2 , C ⊥ 2
are not independent and one has indeed

C inc = (C 2 + C ⊥ 2 )/2. We remark that (C 2 -C ⊥ 2 )/2 = -X 3 -X 4 + Ψ 3 + Ψ 4 .
The additional symmetry can be thus interpreted in terms of a certain decoupling between the DOFs 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 5.1, naturally represent quasi-symmetries when considering the entire spectrum of resonances R 1 . They are indeed broken at some point by weak resonances (see Sect. 5.4.3). Quantities E 2n , C inc , and C 2 are the corresponding QIs of motion. The persistence of the three symmetries under the 30 leading resonances is somewhat surprising. Concerning C inc and C 2 , for example, one might reasonably expect that, since the ISS has 8 DOFs, 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 inc and C 2 , the quantity E 2n is a non-linear function of the action-angle variables. However, as far as stable orbital evolutions are concerned, the convergence of the series expansion of the Hamiltonian is sufficiently fast that the linear LL approximation

E 2 = H 2 + g 5 1 8 • I = C γ 3 , with γ 3 = -ω LL + g 5 1 8 , reproduces
reasonably well E 2n along the flow of H 2n for n > 1. The vector γ 3 will be used in Sect. 5.5, together with γ 1 and γ 2 , to deal with the geometry of the linear action subspace spanned by the QIs. The explicit expressions of these vectors are given in Appendix C.2. We mention that, differently from γ 1 and γ 2 , the components of γ 3 are not integer and they have the dimension of a frequency.

Slow variables

The QIs of motion E 2n , C inc , C 2 are clearly strong candidates for slow variables once evaluated along the orbital solutions. In what follows, to assess the slowness of a dynamical quantity when compared to the typical variations of the action variables, we shall consider the variance of its time series along a numerical solution.

We define the dimensionless QIs

C inc = C inc γ 1 C 0 , C 2 = C 2 γ 2 C 0 , E 2n = E 2n γ 3 C 0 , ( 5.21)
where C 0 stands for the current total AMD of the inner planets, that is, the value of C ecc + C inc at time zero. We stress that, by introducing the unit vectors

γ i = γ i / γ i for i ∈ {1, 2, 3}, one has C inc = C γ 1 /C 0 and C 2 = C γ 2 /C 0 . At degree 2, one also has E 2 = C γ 3 /C 0 .
We then consider the ensembles of numerical integrations of H 4 and H 6 , with very close initial conditions and spanning 100 Gyr in the future, that have been presented in ref. [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. The top row of Fig. 5.2 shows the time evolution over 5 Gyr of the dimensionless QIs and of two components of the dimensionless action vector I = I/C 0 along the nominal orbital solutions of the two ensembles. We subtract from each time series its mean over the plotted time span. The time series are low-pass filtered by employing the Kolmogorov-Zurbenko (KZ) filter with three iterations of the moving average [START_REF] Zurbenko | Kolmogorov-Zurbenko filters in spatiotemporal analysis[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. A cutoff frequency of 1 Myr -1 is chosen to highlight the long-term diffusion that can be hidden by short-time quasi-periodic oscillations. This is in line with our definition of quasi-integrals based on contribution from resonant harmonics only. We remark that C 2 and E 2n show very similar time evolutions along stable orbital solutions, as can be seen in the top row Fig. 5.2. This is explained by the interesting observation that the components of the unit vectors γ 2 and γ 3 differ from each other by only a few percent, as shown in Appendix C.2. However, we stress that the two vectors are in fact linearly independent: C 2 does not depend on the actions X 1 and X 2 , while E 2n does. The two QIs move away from each other when high eccentricities of Mercury are reached, that is, for large excursions of the Mercury-dominated action X 1 .

Table 5.2 -Top of ranking R 2 . First 10 symmetry-breaking resonances of H 10 along the 5-Gyr nominal solution of Gauss's dynamics, that change E 2n , C inc , and C 2 , respectively (see Table 5.1 for details).

i Fourier harmonic

F i O i τ res i ∆ω i E 2n g 1 + g 2 -2g 5 + s 2 -s 7 6 1% 0.018 0.020 0.016 g 1 -2g 2 + g 6 -s 2 +
C inc g 1 + g 2 -2g 5 + s 2 -s 7 6 1% 0.018 0.020 0.016 g 1 -2g 2 + g 6 -s 2 + s 6
6 4% 0.017 0.024 0.008 g 2 -g 6 + s 1 -s 6 4 8% 0.011 0.017 0.003 g 4 -g 6 -2s 3 + 3s 4 -s 6 8 5% 0.010 0.012 0.002 2g 1 -2g 5 + s 1 -s 7 6 1% 0.007 0.008 0.005 4g 1 -3g 2 -g 5 -s 2 + s 7 10 7% 0.006 0.015 0.003 g 3 -g 6 -2s 3 + 3s 4 -s 6 8 2% 0.006 0.011 0.003 4g 1 -g 2 -g 3 -s 1 -2s 2 + s 4 10 2% 0.006 0.009 0.001 2g 1 -g 2 -g 5 + 3s 1 -2s 2 -s 7 10 3% 0.006 0.008 0.002 3g 1 -3g 2 + s 1 -2s 2 + s 7 10 19% 0.006 0.009 0.002

C 2 g 1 + g 2 -2g 5 + s 2 -s 7 6 1% 0.018 0.020 0.016 g 1 -2g 2 + g 6 -s 2 + s 6 6 4% 0.017 0.024 0.008 g 3 -g 6 + s 2 -s 4 4 10% 0.017 0.022 0.009 g 4 -g 6 + s 2 -s 4 4 4% 0.016 0.021 0.007 g 2 -2g 4 + g 6 4 12% 0.016 0.027 0.009 g 1 -g 3 -g 5 + g 6 -s 1 + s 4
6 6% 0.015 0.024 0.011 g 2 + g 3 -3g 4 + g 6 6 3% 0.014 0.017 0.010 g 2 -2g 3 + g 6 + s 3 -s 4 6 3% 0.011 0.018 0.003 g 1 -g 3 -g 4 + g 6 + s 1 -s 2 6 4% 0.011 0.013 0.008 g 2 -g 6 + s 1 -s 6 4 8% 0.011 0.017 0.003

Weak resonances and Lyapunov spectrum

A fundamental result from Table 5.1 is that the symmetries introduced in Sect. 5.4.1 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 5.2 reports the 10 strongest symmetry-breaking resonances that change E 2n , C inc , C 2 , respectively. As in Table 5.1, only harmonics that are resonant for more than 1% of the 5-Gyr time span of the nominal solution of Gauss's dynamics are shown. The leading symmetry-breaking resonances have half-widths of about 0.01 ′′ yr -1 . For each QI, the dominant contribution comes from harmonics involving Fourier modes of the outer planet forcing other than g 5 : the Saturn-dominated modes g 6 , s 6 and the modes g 7 , s 7 mainly associated to Uranus. In the case of C inc , there is also a contribution that starts at about 0.006 ′′ yr -1 with F 8 = 4g 1 -g 2 -g 3 -s 1 -2s 2 + s 4 and comes from high-order internal resonances, that is, resonances that only involve the DOFs of the inner planets.

We remark that the decrease of the resonance half-width with the index of the harmonic in Table 5.2 is steeper for C inc than for E 2n , C 2 , and is accompanied by a greater presence of high-order resonances. This may notably explain why the secular variations of C inc are somewhat smaller in the top row of Fig. 5.2. We finally point out the important symmetry-breaking role of the modes g 7 , s 7 , representing the forcing mainly exerted by Uranus. Differently from what one might suppose, these modes cannot be completely neglected when addressing the long-term diffusion of ISS. This recalls the role of the modes s 7 and s 8 in the spin dynamics of Venus [START_REF] Correia | Long-term evolution of the spin of Venus. II. numerical simulations[END_REF], and is basically a manifestation of the long-range nature of the gravitational interaction.

As stated in Sect. 5.3, a pair of Lyapunov exponents would vanish if there were an exact integral of motion. In 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 Equation (5.12). Comparison of Table 5.2 with the Lyapunov spectrum in Fig. 5.1a shows that the time statistics of the half-widths of the symmetry-breaking resonances of ranking R 2 overlaps with the ensemble distribution of the three smallest FT-LCEs, that is, λ 6 , λ 7 , λ 8 . One can indeed write (5.22) where ∆ω R 2 stands for the half-width of the uppermost resonances of ranking R 2 . Table 5.2 and Fig. 5.1a suggest a relation between the QIs and the smallest Lyapunov exponents:

2πλ 6 ≈ ∆ω R 2 ,
λ 6 , λ 7 , λ 8 ←→ E 2n , C inc , C 2 .
(5.23) Equation (5.23) is not a one-to-one correspondence, nor it should be understood as an exact relation since, for example, λ 6 is not well separated from the larger exponents. Its physical meaning is that the QIs are among the slowest DOFs of the ISS dynamics. Such claim is one of the core points of this work. In Sect. 5.5, we will show its statistical validity in the geometric framework established by a principal component analysis of the orbital solutions. Moreover, Sect. 5.4.4 shows that Eq. (5.23) can be stated more precisely in the case of a simplified dynamics that underlies H 2n . We remark that E 2n , C inc , C 2 constitute a set of three QIs that are independent and nearly in involution, and it is thus meaningful to associate three different Lyapunov exponents with them. On the one hand, the independence is easily checked at degree 2 as the vectors γ 1 , γ 2 , γ 3 are linearly independent. On the other hand, one has the Poisson bracket {C inc , C 2 } = 0, since the two quantities are functions of the action variables only. One also has {E 2n ,

C inc } = {H 2n , C inc } = Ċinc and {E 2n , C 2 } = Ċ2 .
Only weak resonances contribute to these Poisson brackets and the three QIs are therefore nearly in involution.

A new truncation of the Hamiltonian

The fundamental role of the external modes g 6 , g 7 , s 6 , s 7 in Table 5.2 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 5.3 reports the 10 strongest resonances per each broken symmetry. The difference with respect to Table 5.2 is manifest. As g 5 is the only external mode remaining, there are no resonances left that can contribute to the time evolution of E 2n . For the remaining two QIs, the only harmonics that appear in Table 5.3 are of order 8 or higher, and this is accompanied by a significant drop in the half-width of the leading resonances. In the case of C inc , the half-width of the uppermost resonances is now around 0.005 ′′ yr -1 . One can appreciate that the activation times τ res of the resonances do not exceed a few percent, differently from Table 5.2. The most impressive change is however related to C 2 : only harmonics of order 10 appear in Table 5.3 and the half-width of the uppermost resonances drops by two orders of magnitude. We stress that such harmonics are resonant for very short periods of time along the 5 Gyr spanned by the nominal solution of Gauss's dynamics. To retrieve the time statistic of the resonances affecting C 2 we indeed chose to repeat the computations of ref. [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF] by increasing the cutoff frequency of the low-pass filter applied to time series of the action-angle variables from (5 Myr) -1 to 1 Myr -1 (Mogavero and Laskar, 2022, Appendices F.2 and G.5). The filtered time series have then been resampled with a timestep of 50 kyr. Many harmonics shown in Table 5.3 and related to C 2 are resonant for a few timesteps and their time statistic 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. Anyway, the fundamental point here is the drastic reduction in the size of the uppermost resonances with respect to Table 5.2, 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. [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF] to establish the resonant harmonics, so that they do not appear in the tables of this work.

Hamiltonian H • 2n . The implications of Table 5.3 suggest to introduce an additional truncation in the Hamiltonian H 2n . This consists in dropping the harmonics of Eq. (5.6) that involve external modes other than g 5 : (5.24) where φ 1 (t) = -g 5 t and ℓ • = (ℓ 1 , 0, . . . , 0), with ℓ 1 ∈ Z. Consistently with the absence of symmetry-breaking resonances related to E 2n in Table 5.3, the corresponding dynamics admits the exact integral of motion (5.25) which represents the transformed Hamiltonian under the canonical change of variables that eliminates the explicit time dependence in Eq. (5.24). We point out that, as the additional truncation is applied to the action-angle formulation of Eq. (5.6), the external modes other than g 5 still enter the definition of the proper modes of the forced Laplace-Lagrange dynamics [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. The orbital solution arising from H • 2n is initially very close to that of H 2n . A frequency analysis over the first 20 Myr shows that the differences in the fundamental frequencies of the motion between H • 2n and H 2n are of the order of 10 -3 arcsec yr -1 , an order of magnitude smaller than the typical frequency differences between H 4 and H 6 (Mogavero and Laskar, 2021, Table 3). Therefore, even though H • 2n constitutes a simplification of H 2n , it should not be regarded as a toy model. Its dynamics, in particular, still possesses 8 DOFs.

H • 2n (I, θ, t) = k,ℓ 1 H k,ℓ • 2n (I)E j(k•θ+ℓ 1 φ 1 (t)) ,
E • 2n = H • 2n + g 5 4 i=1 (X i + Ψ i ),
We compute the Lyapunov spectrum of the Hamiltonian H • 4 in the same way as described in Sect. 5.3 in the case of H 2n . Since its dynamics turns out to be much more stable then that of H 4 (see Sect. 5.6,Fig. 5.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 5 .1b. Comparing to the Lyapunov spectrum of H 4 , one notices that the distributions of the leading exponents turn out to be quite similar, apart from being more spaced and except for a slight decrease in their median values. However, such a decrease is more pronounced for smaller exponents, and the drop in the smallest exponents is drastic. The smallest one, λ 8 , decreases monotonically, consistently with the fact that E • 4 from Eq. (5.25) is an exact integral of motion. The exponent λ 7 drops by more than an order of magnitude, and apparently begins to stabilize around a few 10 -4 arcsec yr -1 , while λ 6 also reduces significantly, by a factor of three, to about 0.005 ′′ yr -1 . The drop in the smallest exponents agrees remarkably well with that of the half-width of the leading symmetry-breaking resonances when switching from Table 5.2 to Table 5.3. One can indeed write

2πλ 6 ≈ ∆ω R 3 ,C inc , 2πλ 7 ≈ ∆ω R 3 ,C 2 , λ 8 = 0, (5.26)
where ∆ω R 3 ,Q stands for the half-width of the uppermost resonances of ranking R 3 related to the quasi-integral Q. The hierarchy of the three smallest exponents in the spectrum of Fig. 5.1b consistently follows that of the QIs suggested in Table 5.3 by the very different sizes of the leading resonances. In other words, one can state:

λ 6 ←→ C inc , λ 7 ←→ C 2 , λ 8 ←→ E • 2n .
(5.27)

These one-to-one correspondences are a particular case of Eq. (5.23) and support the physical intuition behind it. In Sect. 5.5, we will prove the validity of Eq. (5.27) in the geometric framework established by a principal component analysis of the orbital solutions of H • 2n .

Numerical integrations. We compute ensembles of 1080 orbital solutions of the dynamical models H • 4 and H • 6 , with initial conditions very close to the nominal ones of Gauss's dynamics and spanning 100 Gyr in the future. This closely follows what did in ref. [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF] in the case of the models H 2n . The bottom row of Fig. 5.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. (5.27) is manifest. The quantity C 2 has secular variations much slower than C inc , while the latter is itself slower with respect to its counterpart in the orbital solutions of H 2n . We remark that, as E • 2n is an exact integral of motion for the model H • 2n , we do not plot it. From Fig. 5.2 it is also evident how difficult can be the retrieval of the short-lasting resonances affecting C 2 from a solution of H • 2n spanning only a few Gyr.

The hierarchy of the QIs is confirmed by a statistical analysis in Appendix C. 3. Figure C.3 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 C.4 details the growth of the QI dispersion over time. As suggested by Table 5.3, the drop in the diffusion rates of the QIs when switching from H 2n to H • 2n is manifest.

Statistical detection of slow variables

Section 5.4 shows how the slow-fast nature of the ISS dynamics, indicated by the Lyapunov spectrum, emerges from the quasi-symmetries of the resonant harmonics of the Hamiltonian. QIs of motion can be introduced semi-analytically and they constitute slow quantities when evaluated along stable orbital solutions. In this section, we will 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 (PCA). We will 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.

Principal component analysis

PCA is a widely used classical technique for multivariate analysis [START_REF] Pearson | LIII. On lines and planes of closest fit to systems of points in space[END_REF][START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]. For a given dataset, PCA aims to find an orthogonal linear transformation of the variables such that the new coordinates offer a more condensed and representative view of the data. The new variables are called principal components (PCs). They are uncorrelated and ordered according to decreasing variance: the first PC and last one have respectively the largest and the smallest variance of any linear combination of the original variables. While one is typically interested in the PCs of largest variance, in this work we employ the variance of the time series of a dynamical quantity to assess its slowness when compared to the typical variations of the action variables (see Sect. 5.4.2). We will 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 e.g. [START_REF] Jolliffe | Principal component analysis[END_REF][START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]. Let I(t) = (X(t), Ψ(t)) be the 8-dimensional time series of the action variables evaluated along a numerical solution of the equations of motion. As in Sect. 5.4.2, we apply the KZ low-pass filter with three iterations of the moving average and a cutoff frequency of 1 Myr -1 to obtain the filtered time series Î(t) [START_REF] Zurbenko | Kolmogorov-Zurbenko filters in spatiotemporal analysis[END_REF][START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF]. In this way, the short-term quasi-periodic oscillations are mostly suppressed, which better reveals the chaotic diffusion over longer timescales. We finally define the mean-subtracted filtered action variables over the time interval [t 0 , t 0 + T ] as Ĩ(t) = Î(t) -n -1 n-1 i=0 Î(t 0 + i∆t), where the mean is estimated by discretization of the time series with a sampling step ∆t such that T = (n -1)∆t. The discretized time series over the given interval is stored in a 8 × n matrix: D = [ Ĩ(t 0 ), Ĩ(t 0 + ∆t), . . . , Ĩ(t 0 + (n -1)∆t)].

(5.28)

The PCA of the data matrix D consists in a linear transformation P = A T D, where A is an 8×8 orthogonal matrix (i.e. A -1 = A T ) defined as follows. By writing A = [a 1 , . . . , a 8 ], the column vectors a i ∈ R 8 are chosen to be the normalized eigenvectors of the sample covariance matrix, in order of decreasing eigenvalues: (n -1) -1 DD T = AΣA T , where Σ = diag(σ 1 , . . . , σ 8 ) and σ 1 ≥ • • • ≥ σ 8 . The PCs are defined as the new variables after the transformation, that is, PC i = a i • I with i ∈ {1, . . . , 8}. The uncorrelatedness and the ordering of the PCs can be easily seen from the diagonal form of their sample covariance matrix, (n -1) -1 PP T = Σ, from which it follows that the variance of PC i is σ i . Among all the linear combinations in the action variables I, the last PC, i.e. PC 8 , has the smallest variance over the time interval [t 0 , t 0 + T ] of a given orbital solution. The second last PC, i.e. PC 7 , has the second smallest variance and is uncorrelated with PC 8 , and so on. It follows that the linear subspace spanned by the last k PCs is the k-dimensional subspace of minimum variance: the variance of the sample projection onto this subspace is the minimum among all the subspaces of the same dimension. These properties indicate that the last PCs provide a pertinent statistical definition of slow variables along numerically integrated solutions of a dynamical system. The linear structure of the PCA, in particular, seems adapted to quasi-integrable systems close to a quadratic Hamiltonian, like the ISS. In such a case, one may reasonably expect that the slow variables are, to a first approximation, linear combinations of the action variables. We remark that the mutual orthogonality allows to associate a linear DOF 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:

D agg = [D 1 , D 2 , . . . , D m ],
where D i is the data matrix of Eq. (5.28) for the i-th solution. Since this work deals with a non-stationary dynamics, as the ISS ceaselessly diffuses in the phase space [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF], 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.

Principal components and quasi-integrals

Both the QIs and the last PCs represent slow variables, but are established through two different methods. Equations (5.23) and (5.27) claim that the QIs found semi-analytically in Sect. 5.4 are among the slowest DOFs of the ISS dynamics. This naturally suggests to compare the three QIs with the three last PCs retrieved from numerically integrated orbital solutions. In this part, we first introduce the procedure that we implement to establish a consistent and systematic correspondence between QIs and PCs. We then present both a visual and a quantitative geometric comparison between them.

Tweaking the QIs

The three last components PC 8 , PC 7 , PC 6 are represented by the set of vectors S PCs = {a 8 , a 7 , a 6 } belonging to R 8 . By construction, these PCs have a linear, hierarchical, and orthogonal structure. In other words: the PCs are linear combinations of the action variables I; denoting by the order of statistical variance, one has PC 8 PC 7 PC 6 ; the unit vectors (a i ) 8

i=6 are orthogonal to each other. On the other hand, the QIs of motion C inc , C 2 , E 2n do not possess these properties. Therefore, we will adjust them in such a way to reproduce the same structure.

Linearity. While C inc and C 2 are linear functions of the action variables, E 2n is not when n > 1. Nevertheless, as explained in Sect. 5.4.1, as far as one considers stable orbital solutions, the linear LL approximation E 2 = γ 3 • I reproduces E 2n reasonably well. Therefore, we shall consider the three linear QIs of motion C inc , C 2 , E 2 , which are respectively represented by the set of R 8 -vectors S QIs = {γ 1 , γ 2 , γ 3 }. In this way, the 3-dimensional linear subspaces of the action space spanned by the sets S QIs and S PCs can be compared.

Ordering. We shall define a set of QIs that are ordered by statistical variance, as it is the case for the PCs. We follow two different approaches according to model H • 2n in Eq. (5.24) or H 2n in Eq. (5.6) (clearly n > 1).

H 2n

Since the leading resonances affecting the QIs in Table 5.2 have comparable sizes, there is no clear order of statistical variances that can be inferred. We then implement a systematic approach that orders the QIs by simply inheriting the ordering of the PCs. More precisely, we define a set of ordered vectors S ′ QIs through the projections of the three last PCs onto the linear subspace generated by the QIs: S ′ QIs = {proj S QIs (a 8 ), proj S QIs (a 7 ), proj S QIs (a 6 )}6 . As a result, the new set of QIs mirrors the hierarchical structure of the PCs. We stress that S ′ QI spans the same subspace of R 8 as S QI , since the ordered QIs are just linear combinations of the original ones.

Orthogonality. We apply the Gram-Schmidt process to the ordered set S ′ QIs to obtain the orthonormal basis S ′′ QIs = {α 1 , α 2 , α 3 }. The set S ′′ QIs clearly spans the same subspace as S QIs . Moreover, the Gram-Schmidt process preserves the hierarchical structure, that is, the two m-dimensional subspaces spanned by the first m ≤ 3 vectors of S ′ QIs and S ′′ QIs , respectively, are identical.

In the end, we obtain a linear, ordered, and orthogonal set of modified QIs of motion {QI 1 , QI 2 , QI 3 }, where QI i = α i • I.

Hamiltonian H •

4 . The modified QIs can be explicitly derived in this case and comprise interpretable physical quantities. One has QI 1 proportional to C 2 and QI 2 proportional to C ⊥ 2 . Moreover, QI 3 is the component of E 2 that is orthogonal to both C 2 and C ⊥ 2 . Figure 5.3 shows the comparison between the modified QIs and the corresponding PCs for three different time intervals along sol. #1 of H • 4 (see Fig. 5.2 bottom left for its time evolution). The agreement of the pairs (QI 1 , PC 8 ), (QI 2 , PC 7 ), and (QI 3 , PC 6 ) across different intervals is manifest and even impressive. One can appreciate that, the "slower" the PC, the more similar it is to its corresponding QI. The overlap between the modified QIs and the three last PCs means that the QIs of motion span the slowest 3-dimensional linear subspace of the action space. Therefore, to a linear approximation, they represent the three slowest DOFs of the H Hamiltonian H 4 . In this case, the QIs of motion C inc , C 2 , E 2 do not show a clear hierarchical structure in terms of statistical variance. Therefore, we consider the whole subspace spanned by the three QIs with respect to that spanned by three last PCs. Since it is not easy to visually compare two 3-dimensional subspaces of R 8 , we compare their basis vectors instead. The basis α 1,2,3 of modified quasi-integrals QI 1,2,3 is computed according to the algorithm presented in Sect. 5.5.2.

Figure 5.4 presents the comparison between the modified QIs and the corresponding PCs across three different time intervals of three solutions of H 4 (see Fig. 5.5 for their time evolution). The first two, sols. #1 and #2, show thorough agreement between the pairs of QIs and PCs across all intervals, which indicates close proximity between the two subspaces V QIs = span(S QIs ) and V PCs = span(S PCs ). One can appreciate that the directions of the basis vectors are quite stable. The last component PC 8 , in particular, remains close to C inc . In line with the discussion in Sect. 5.4.3, the slowest DOF of H 4 can be thus deduced to be close to C inc , to a linear approximation at least. Such a result shows how interesting physical insight can be gained through the PCA. Some changes in the basis vectors can arise, however, as for the first time interval of sol. #2. This may be expected from a dynamical point of view. Differently from H • 4 , there is no pronounced separation between the slowest DOFs at the bottom of the Lyapunov spectrum in Fig. 5.1a: the marginal distributions of consecutive exponents can indeed touch or overlap each other. Therefore, the hierarchy of slow variables is not as frozen as in H • 4 and it can change along a given orbital solution.

Solutions #1 and #2 represent typical orbital evolutions. If the same time intervals of all the 1080 solutions are stacked together to form an aggregated sample on which the PCA is applied, the features mentioned above persist: the agreement between QIs and PCs, the stability of the basis vectors and the similarity between PC 8 and C inc (see Fig. 5.4). Once again, the PCA confirms that the subspace spanned by the three QIs is overall close to the slowest 3-dimensional linear subspace of the action space. Therefore, to a linear approximation, they represent the three slowest DOFs of the H 4 dynamics. We remark that the slowness of the 3-dimensional subspace spanned by the QIs is a much stronger constraint than the observation that each QI is a slow variable. To give an example, let Q = q • I be a slow variable with unit vector q. If ǫ is an arbitrary small vector, i.e. ǫ ≪ 1, then Q ′ = ( q + ǫ) • I can be also considered as a slow variable, whereas the normalized difference of two quantities, ǫ • I, is generally not. Therefore, the linear subspace spanned by Q and Q ′ , that is, by q and ǫ, is not a slow 2-dimensional subspace.

Solution #3 in Fig. 5.4 represents an edge case (see Fig. 5.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 DOFs can be more difficult in two rare possibilities. Firstly, if the change in a QI accumulates continually in one direction, its variance can inflate over a long time interval. This is the case for the interval [0, 5] Gyr of sol. #3. Secondly, the variance of a variable that is typically fast can suddenly dwindle during a certain period of time, for example Ψ 3 over the interval [START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF]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 DOFs that are typically faster. This is reflected in the mismatch of QI 3 and PC 6 on the last two time intervals of sol. #3 in Fig. 5.4. We remark that PC 8,7 are still relatively close to QI 1,2 , which indicates that the slowest 2-dimensional subspace spanned by PC 8,7 still resides inside the QI subspace. It should be stressed that this disagreement between QIs and PCs does not mean that the QIs are not slow variables in this case. The mismatch has a clear dynamical origin instead. The resonance tables of this work has 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 [START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. 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 DOFs 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. 5.1a shows that the marginal distributions of the exponents λ 5 and λ 6 , for example, are not separate but overlap each other.

Distance between the subspaces of PCs and QIs

The closeness of the two 3-dimensional linear subspaces V PCs , V QIs ⊂ R 8 spanned by the sets of vectors S PCs and S QIs , respectively, can be quantitatively measured in terms of a geometric distance. This can be formulated using the principal (canonical) angles [START_REF] Jordan | Essai sur la géométrie à n dimensions[END_REF][START_REF] Van Loan | Matrix computations[END_REF][START_REF] Ye | Schubert varieties and distances between subspaces of different dimensions[END_REF].

Let A and B be two sets of m ≤ n independent vectors in R n . The principal vectors (p k , q k ) m k=1 are defined recursively as solutions to the optimization problem maximize p • q subject to: p ∈ span(A), q ∈ span(B),

p = 1, q = 1, p • p i = 0, q • q i = 0, i = 1, . . . , k -1,
(5.29) for k = 1, . . . , m. The principal angles 0 ≤ θ 1 ≤ • • • ≤ θ m ≤ π/2 between the two subspaces span(A) and span(B) are then defined by

cos θ k = p k • q k , k = 1, . . . , m.
(5.30)

The principal angle θ 1 is the smallest angle between all pairs of unit vectors in span(A) and span(B); the principal angle θ 2 is the smallest angle between all pairs of unit vectors that are orthogonal to the first pair; and so on. Given the matrices defining the two subspaces, the principal angles can be computed from the singular value decomposition of their correlation matrix. The result is the canonical correlation matrix diag(cos θ 1 , . . . , cos θ m ). This cosine-based method is often ill-conditioned for small angles. In such case, a sinebased algorithm can be employed [START_REF] Björck | Numerical methods for computing angles between linear subspaces[END_REF]. In this work, we use the combined technique detailed in [START_REF] Knyazev | Principal angles between subspaces in an Abased scalar product: algorithms and perturbation estimates[END_REF].

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 [START_REF] Ye | Schubert varieties and distances between subspaces of different dimensions[END_REF]:

d(A, B) = 1 m m k=1 sin 2 θ k 1/2 .
(5.31)

The distance is null if A and B are the same subspace and equal to 1 when they are orthogonal. We will use this metric to show that the subspace closeness suggested by Figs. 5.3 and 5.4 is indeed statistically significantly. More precisely, we will provide evidence against the null hypothesis that the distribution of distances between V PCs and V PCs , arising from the H • 4 and H 4 dynamics, coincides with that of randomly drawn 3dimensional subspaces of R 8 . The PDF of the distance between two random 3-dimensional subspaces of R 8 is shown in Fig. 5.6 in blue color (such random subspaces can be easily generated by sampling sets of 3 vectors uniformly on the unit 7-sphere, [START_REF] Muller | A Note on a Method for Generating Points Uniformly on N-Dimensional Spheres[END_REF]. While the range of possible distances is [0,1], the distribution concentrates on the right side of the interval, with a probability of approximately 99.3% that the distance is larger than 0.6. In this regard, we remark that the notion of distance in high-dimensional spaces is very different from our intuition in a 3-dimensional world. If we draw randomly two vectors in a very high-dimensional space, it is extremely likely that they will be close to mutual orthogonality.

The upper panel of Fig. 5.6 shows in green color the PDF of the distance between V PCs and V QIs arising from the time interval [0, 5] Gyr of the 1080 orbital solutions of model H • 4 . In the lower panel, we consider a larger ensemble of 10 800 solutions of model H 4 spanning the same time interval [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF], and plot the corresponding PDF of the distance between V PCs and V QIs . In both cases, the distance stemming from the aggregated sample of all the solutions is indicated by a vertical dark green line. We also dynamics which is used to infer the QIs. Solution #3 in Fig. 5.4 represents in this sense an edge case of the distance distribution, while sol. #2 is a typical solution close to the PDF median.

Implications on long-term stability

The existence of slow variables can have crucial implications on the stability of the ISS. The QIs of motion can effectively constraint in an adiabatic way the chaotic diffusion of the planet orbits over long timescales, forbidding in general a dynamical instability over a limited time span, e.g. several billions of years. Here we give compelling arguments for such a mechanism.

Figure 5.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 Sect. 5.4.4. 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 [START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF]. We also report the same CDF for the models H 4 and H 6 , which we recently computed in ref. [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. One can appreciate that the time corresponding to a probability of instability of 1% is greater than 100 Gyr for the H We also observe that, to a linear approximation, the knowledge of C inc and E 2 allows to bound the variations of the action variables X, Ψ. Recalling that the actions are positive quantities, from Eq. (5. [START_REF] Lantink | Milankovitch cycles in banded iron formations constrain the Earth-Moon system 2.46 billion years ago[END_REF] one sees that fixing a value of C inc puts an upper bound to the variations of the inclination actions Ψ. As a consequence, at degree 2 in eccentricities and inclinations, fixing a value of the QI

E 2 = γ 3 • I = γ ecc 3 • X + γ inc 3 • Ψ, (5.32)
with γ 3 = (γ ecc 3 , γ inc 3 ), also bounds the upper variations of the eccentricity actions X, since the components of γ ecc 3 have all the same sign, as those of γ inc 3 (see Appendix C.2). This is an important point, as we stated in Sect. 5.1 that the lack of any bound on the chaotic variations of the planet orbits is one of the reasons that complicate the understanding of their long-term stability. We remark that the secular planetary phase space can be bound by fixing the value of the total AMD, that is, C ecc + C inc [START_REF] Laskar | Large scale chaos and the spacing of the inner planets[END_REF]. A statistical study of the density of states that are a priori accessible can then be realized [START_REF] Mogavero | Addressing the statistical mechanics of planet orbits in the solar system[END_REF]. It is not, however, fully satisfying to consider a fixed value of total AMD of the ISS, as we showed 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

Discussion

This work introduces a framework that naturally justifies the statistical stability shown by the ISS over a timescale comparable to its age. Considering a forced secular model of the inner planet orbits, the computation of the Lyapunov spectrum indicates the existence of very different dynamical timescales. Using the computer algebra system TRIP, we systematically analyze the Fourier harmonics of the Hamiltonian that become resonant along a numerically integrated orbital solution spanning 5 Gyr. We uncover three symmetries that characterize the strongest resonances and that are broken by weak resonant interactions. These quasi-symmetries generate three QIs of motion that represent slow variables of the secular dynamics. The size of the leading symmetry-breaking resonances suggests that the QIs are related to the smallest Lyapunov exponents. The claim that the QIs are among the slowest DOFs 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.

The inner Solar System among classical quasi-integrable systems

It is valuable to contextualize the dynamics of the ISS in the class of classical quasiintegrable systems. A comparison with the FPUT dynamics, in particular, deserves to be made. In the FPUT problem, the proper modes of the chain of oscillators remain far from the energy equipartition expected from statistical mechanics for a very long time [START_REF] Fermi | Studies of Nonlinear Problems[END_REF]. One way to explain the lack of energy equipartition reported by Fermi and collaborators is through the closeness of FPTU problem to the integrable Toda dynamics [START_REF] Hénon | Integrals of the Toda lattice[END_REF][START_REF] Flaschka | The Toda lattice. II. Existence of integrals[END_REF][START_REF] Manakov | Complete integrability and stochastization of discrete dynamical systems[END_REF]. This translates in a very slow thermalization of the action variables of the Toda problem and of the corresponding integrals of motion along the FPTU flow [START_REF] Manakov | Complete integrability and stochastization of discrete dynamical systems[END_REF][START_REF] Ferguson | Nonlinear Normal Modes for the Toda Chain[END_REF][START_REF] Benettin | The Fermi-Pasta-Ulam Problem and Its Underlying Integrable Dynamics[END_REF][START_REF] Goldfriend | Equilibration of quasi-integrable systems[END_REF][START_REF] Christodoulidi | Stages of dynamics in the Fermi-Pasta-Ulam system as probed by the first Toda integral[END_REF][START_REF] Grava | Adiabatic Invariants for the FPUT and Toda Chain in the Thermodynamic Limit[END_REF]. 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 found, 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 FPTU 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 found 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 [START_REF] Ford | Equipartition of Energy for Nonlinear Systems[END_REF][START_REF] Onorato | Route to thermalization in the α-Fermi-Pasta-Ulam system[END_REF][START_REF] Pistone | Universal route to thermalization in weakly-nonlinear one-dimensional chains[END_REF]. A hierarchy of Lyapunov exponents spanning several orders of magnitude, in particular, should be common among this class of systems (e.g. [START_REF] Malishava | Lyapunov Spectrum Scaling for Classical Many-Body Dynamics Close to Integrability[END_REF].

Methods

The long-term dynamics of the ISS is described by a moderate but not small number of DOFs, 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 developed 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. [START_REF] Mogavero | Long-term dynamics of the inner planets in the Solar System[END_REF][START_REF] Mogavero | The origin of chaos in the Solar System through computer algebra[END_REF]. This allows us to introduce QIs of motion in a 16-dimensional dynamics by analyzing how action-space symmetries are progressively broken by resonant interactions. Our effective method based on the time statistics of resonances arising along a single, very long numerical integration is alternative to formal approaches that define QIs via series expansions (e.g. [START_REF] Contopoulos | A third Integral of Motion in a Galaxy[END_REF][START_REF] Kruskal | Asymptotic Theory of Hamiltonian and other Systems with all Solutions Nearly Periodic[END_REF]. 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 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 resonances7 .

In the context of dynamical systems with a number of DOFs that is not small, this work also considers an approach based on PCA. The role of this statistical technique can be twofold. We used PCA as an independent test to systematically validate the slowness of the QIs. While being introduced semi-analytically as dynamical quantities that are not affected by the leading resonances, they can indeed be related to the last PCs. By extension, the first PCs should probe the directions of the main resonances. This leads to a second potential application of the PCA, which should offer a way to retrieve the principal resonant structure of a dynamical system. In this sense, PCA represents a tool to systematically probe numerical integrations of a complex dynamics and distil important hidden insights. We emphasise that PCA is the most basic linear technique of dimensionality reduction and belongs to the more general class of the unsupervised learning algorithms. There are more sophisticated methods of feature extraction, that can be more robust (e.g. [START_REF] Candès | Robust principal component analysis[END_REF][START_REF] Markopoulos | Optimal algorithms for L_{1}subspace signal processing[END_REF] and can incorporate nonlinearity [START_REF] Lee | Nonlinear dimensionality reduction[END_REF]. 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.

Chapter 6 Conclusion

In this thesis, I study the long-term evolution of the Solar System in the past and in the future. The two time directions represent two important problems: the recovery of the past history and a possible future destabilization. The two problems are both important and different in their own ways: one is multi-disciplinary and requires practical geology expertise, the other is a celestial mechanics problem that is often treated analytically in Hamiltonian formalism. Nevertheless, due to the chaotic nature of the Solar System, the approaches to both problems share the common denominator of statistical analysis, which is based on a large number of long-term orbital solutions.

The first part of this thesis, presented in chapter 3 and published in [START_REF] Hoang | Chaotic diffusion of the fundamental frequencies in the Solar System[END_REF], showcases the first step in the overarching goal of recovering the Solar System orbital history by geological constraints. We develop the estimation of the chaotic uncertainty of the fundamental frequencies of the Solar System beyond 60 Myr, via their probability density function (PDF) by the kernel density estimator. We benchmark and implement the statistical estimators, namely moving block bootstrap, to gauge the error of the density estimation of the frequencies. The significance of these PDFs are twofold. First, they provide a global view from an astronomical perspective to assess the geological signals. Instead of working with only one or several orbital integrations, cyclostratigrapher can now make use of hundreds of thousands astronomical solutions to statistically evaluate their results. We illustrate the applications of our solutions on the two geological data sets: the Newark-Hartford records [START_REF] Olsen | Mapping Solar System chaos with the Geological Orrery[END_REF] and the Libsack core [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF]. The constraint of chaos provided is essential to hypothesis testing not only for geological findings but also for promising but unproven theory of the Solar System. For example, [START_REF] Spalding | An orbital window into the ancient Sun's mass[END_REF] propose that the Sun lost several percents of its mass during its lifetime, to explain for the faint young Sun paradox 1 . Their simulations show that the astronomical forcing g 2 -g 5 roughly scale linearly with the Sun mass and can thus be used as its proxy. The chaotic constraint that we provided for the secular frequencies is necessary in order to test this hypothesis. These are two general routes for the future works in this direction. In particular, there is a recent surge on the cyclostratigraphy literature of new geological confirmation of the resonance transition (e.g. [START_REF] Ma | Theory of chaotic orbital variations confirmed by Cretaceous geological evidence[END_REF]. A robust and consistent assessment of the transition probability, which is currently missing, is undoubtedly necessary for the community. With a large number of orbital solutions that we have, the proper response to this question, which is illustrated in the Sect. 3.6.3, is straight-forward and actionable. The second direction is also our primary objective, that is to effectively combine the geological signal with the chaotic astronomical solutions to provide constraints that are beneficial to both. A potential and promising route is to apply the Bayesian method, pioneered by [START_REF] Meyers | Proterozoic Milankovitch cycles and the history of the solar system[END_REF], who use simple and non-realistic Gaussian functions as priors for secular frequencies. This Bayesian approach can already be improved to a large extent by simply adopting our PDFs of secular frequencies as priors.

The second part of this thesis, presented in chapter 4 (published in Hoang, Mogavero and Laskar 2022) and chapter 5 (submitted for publication), studies the destabilization of the inner Solar System. The chapter 4 details a comprehensive semi-analytical study of instability statistics on a hierarchy of models. This chapter also marks the transition in our work from the descriptive statistics towards the exploratory statistics of the ISS, where physical insight can be extracted from the data. By either intentionally including or omitting secular terms to construct new secular Hamiltonians and comparing the difference in their instability statistics, we can show the impact of specific terms on the destabilization of the ISS. For instance, we have surprisingly found that the non resonant secular harmonics at degree 6 in eccentricity and inclination, which are often overlooked in the literature, are essential to correctly reproduce the destabilization timescale. Without these terms, the secular Hamiltonian at degree 4 is overly stable, with the instability probability within 5 Gyr dropping by several orders of magnitude.

The chapter 5 asserts the slow-fast nature of the ISS on the chaotic timescale. The first sign of the disparate timescale separation can be seen from the discrete spectrum of Lyapunov exponents spanning two orders of magnitude, the largest of which is the maximal Lyapunov exponent of 5 Myr -1 , and the smallest exponent is around 500 Myr -1 . Moreover, we show concrete evidences for the slow-fast dynamics by identifying the slow variables. From the ranking of active resonances (see Mogavero and Laskar, 2022, table 2), we find three symmetries among the leading resonances, from which three quasi-integrals of motion can be defined. We use the principle component analysis to look for and confirm that these quasi-integrals represent the slowest degree of freedom of the dynamics.

The role of these quasi-integrals for the stability of the inner planets will be demonstrated in the future work. I summarize the key idea as follows. The quasi-integrals set up constraints for the dynamics and thereby crucially slow down system in its the pathway towards the destabilization, similar to the role of the slowly-varying Hamiltonian in the simplified dynamics of Mercury of [START_REF] Batygin | Chaotic disintegration of the inner solar system[END_REF]. The quasi-integrals roughly define a bounded region; in this region, the dynamics is stable with rapid chaotic motion over the timescale of Lyapunov time. Due to the slow diffusion of the quasi-integrals, the bounded region accordingly expands slowly, so the dynamics can be stable for billions of years. When the bounded region touches the resonance g 1 -g 5 , and opens the possibility of activating it, the system can be destabilized. In this picture, there are two important factors controlling the instability statistics: the diffusion rate of the quasi-integrals and the relative distance of the system to the destabilizing resonance g 1 -g 5 . While the former Appendix A

Lagrange-Laplace matrix

The matrices A and B in the Laplace-Lagrange Hamiltonian (2.20) are (N x N) real matrices. The components of A are [START_REF] Laskar | Stability of the planetary three-body problem[END_REF]:

A jj = j-1 k=1 n j m k m 0 C 3 a k a j + N k=j+1 n j m k a j m 0 a k C 3 a j a k , A jk =    2n j m k m 0 a j a k C 2 a j a k , if j < k. 2n j m k m 0 C 2 a k a j , if j > k. (A.1)
The components of B:

B jj = - j-1 k=1 n j m k m 0 C 3 a k a j - N k=j+1 n j m k a j m 0 a k C 3 a j a k , B jk = j-1 k=1 n j m k m 0 C 3 a k a j + N k=j+1 n j m k a j m 0 a k C 3 a j a k . (A.2)
The functions C 2 (α) and C 3 (α) are: where we introduced the diagonal matrix D = diag(0, 1, 0, 1), with 0 = (0, 0, 0, 0) and 1 = (1, 1, 1, 1). It should be noted that the truncations behind the models L 6 and L 4 are defined with respect to the classical variables z k , ζ k , differently from the models H 2n which result from the expansion of H in the complex Poincaré variables x k ∝ (1 -(1 -e 2 k ) 1/2 ) 1/2 E j̟ k and y k ∝ (1 -e 2 k ) 1/4 sin(i k /2) E jΩ k . We define ensembles of initial conditions by slightly varying a single variable of an inner planet at a time, while keeping other variables identical to their reference values, as shown in Table B.1. For the integrations over 100 Gyr, we use the initial conditions varied from the variables (k i = e i cos ̟ i ) i=1,4 of the four inner planets, except for the solutions of L 4 , where only those varied from k 1 are used. For the solutions computed over 5 Gyr, the variables (e i ) i=1,4 are varied to obtain the initial conditions. The solutions integrated up to 100 Gyr are included in the analysis of the statistics of the first 5 Gyr. Equations (B.1) and (B.2) are integrated from these ensembles of initial conditions to obtain the solutions of L 6 and L 4 .

B.2 Statistics with different thresholds of Mercury's eccentricity

We compute the CDFs of the first hitting time of Mercury's eccentricity at the three levels 0.7, 0.8, and 0.9, in order to test the dependency of the instability statistics on different thresholds. The results are shown in Fig. B.1. Up to 5 Gyr, when the instability constitutes a rare event, the models of degree higher than 6 show consistency across high values of eccentricity. The difference between the CDFs of the three thresholds is relatively significant for the models at degree 6 (H 6 , L 6 ), and even more so at degree 4 (H 4 , L 4 ). For H 6 and L 6 , only about half of the integrations exceeding 0.7 also goes beyond 0.9 in 5 Gyr. It should be noted that if Mercury's eccentricity goes beyond 0.9, it is likely that a catastrophic event will shortly ensue, whether it is a secular collision (ML21) or a numerical instability in the truncated dynamics. Therefore, the expected time that a solution of H 6 spends in an unstable state of high Mercury eccentricity is longer, which makes H 6 a prime model for the study of the unstable states of the ISS. Over a longer timescale of 100 Gyr, when the destabilisation is no longer a rare event, the difference of the CDFs with respect to the choice of the eccentricity threshold is small for the models at degree 4 and negligible for the rest.

| p h -p(x)|, from the sample of the B bootstrap KDEs we can compute an asymptotically valid estimation of the piecewise confidence interval CI 1-α (x), defined as:

P (| p h -p(x)| < CI 1-α (x)) = 1 -α. (B.5)
Boundary correction. Kernel density estimation of a PDF on a finite interval can be affected by non-negligible bias at the boundaries. In our work, the interval is defined by the total integration time, that is, [0, 100] Gyr in Fig. 4.5. The nature of the two boundaries is different, and they should be treated differently. At t = 0, the integrations start closely around a nominal value of e 1 ≈ 0.2, therefore the PDF of the first hitting time of e 1 = 0.7 should be 0 when t = 0. This constraint suggests the log-transformation of the sample before applying the KDE [START_REF] Charpentier | Log-transform kernel density estimation of income distribution[END_REF]. The boundary at 100 Gyr has no similar constraints, and we employ a pseudodata method to correct the bias [START_REF] Cowling | On pseudodata methods for removing boundary effects in kernel density estimation[END_REF]. The idea is to use the original dataset to generate fictitious data outside the interval of interest. Let X (1) < • • • < X (n) be the order statistics of the data X 1 , . . . , X n on the interval [0,[START_REF] Lithwick | Theory of secular chaos and Mercury's orbit[END_REF]. The extra data points generated in the range (-∞, 0) are defined by the three-point rule:

X (-i) = -6X (i) + 4X (2i) -3X (3i) .
(B.6)

To adapt the upper limit of the interval [0, 100] Gyr to this rule, we simply transform the data as X (i) → (100 -X (i) )/100. The pseudodata are then generated according to Eq. (B.6), and the ensemble is back-transformed at the end. The number of pseudodata points is taken to be about 10% of the sample size.

C.3 Ensemble distributions of the quasi-integrals over time

To retrieve the long-term statistical behavior of the QIs we consider the ensembles of 1080 numerical integrations of the dynamical models H 4 and H 6 , with very close initial conditions and spanning 100 Gyr in the future, that have been presented in ref. [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF]. We also consider the similar ensembles of solutions for the simplified Hamiltonians H • 4 and H • 6 that have been introduced in Sect. 5.4.4. We report in Fig. C.3 the time evolution of the ensemble PDFs of the low-pass filtered dimensionless QIs and dimensionless actions X 1 , Ψ 3 for the different models (the cutoff frequency of the time filter is set to 1 Myr -1 , as in Sec. 5.4.2). 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 only takes into account the stable orbital solutions, that is, those solutions whose running maximum of Mercury eccentricity is smaller than 0.7 [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF] 2n when compared to H 2n . We recall that E • 2n is an exact integral of motion for the model H • 2n (see Sect. 5.4.4) and its PDF has null dispersion.

  27) D a and D b are diagonal matrices diag(g LL 1 , . . . , g LL N ) and diag(s LL 1 , . . . , s LL N ), which are made up of eigenvalues of A and B. The constant g LL k and s LL k are known as the Lagrange-Laplace secular frequencies. The equations

  3.8).

  .22) This transition corresponds to a change from the 2:1 resonance to the 1:1 resonance of two secular terms, g 4 -g 3 and s 4 -s 3 . Ma et al., 2017 found a sudden change in the period of a long cycle from 2.4 Myr to 1.2 Myr in the Libsack core of the Cretaceous basin from around -90 Myr to -83 Myr.

  Figure 5.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. 5.2 is confirmed by a statistical analysis in Appendix C.3. Figure C.3 shows the time evolution of the distributions of the same quantities as Fig. 5.2 over the stable orbital solutions of the entire ensembles of 1080 numerical integrations of ref. (Hoang et al., 2022). Figure C.4 details the growth of the QI dispersion over time.

  2n and H 2n relates to the smallest Lyapunov exponents (Fig. 5.1), and this is accompanied by a much slower diffusion of the QIs for H • 2n (Figs. 5.2, C.3 and C.4), Fig. 5.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 towards the dynamical instability in the H • 2n model derives from neglecting the external forcing mainly exerted by Saturn, Uranus and Neptune.

  are Laplace coefficients, which arises from the expansion of the inverse of a term related to distance; that is:(1 + α 2 -2α cos θ) α 2 -2α cos θ) k dθ. (A.5)Higher order planetary Hamiltonian should be referred to[START_REF] Laskar | Stability of the planetary three-body problem[END_REF] to N ǫ 10 -11 5000 100 Gyr e i -N ǫ to N ǫ 10 -11 10000 5 Gyr Table B.1 -Offsets of the initial variables k i = e i cos ̟ i and eccentricity e i , with i ∈ {1, 2, 3, 4} corresponding to the inner planets {Mercury, Venus, Earth, Mars}. Different initial conditions correspond to offsets of nǫ in a single variable of a single planet for n = -N, . . . , N , while other variables are kept to their nominal values. Each initial condition is used to compute a solution over the time interval [0, T ].

  Figure C.3 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. C.4, where we report the time evolution of the interquartile range (IQR) of their distributions. After a transient phase lasting about 100 Myr and characterized by the exponential separation of close trajectories, the time growth of the IQR follows a power law typical of diffusion processes. Figures C.3 and C.4 clearly show the slower diffusion of C inc and C 2 in the model H •

Table 4 . 1 -

 41 expand the secular Hamiltonian to degree 4 in Probability P (sup t≤5 Gyr e 1 (t) ≥ e max ) in percent and its 90% confidence interval, where e 1 is Mercury's eccentricity, for the dynamical models H 2n , H, L 2n , and LG09. LG09 represents the 2 501 direct integrations of[START_REF] Laskar | Existence of collisional trajectories of Mercury, Mars and Venus with the Earth[END_REF]. H denotes the 10 560 orbital solutions of Gauss' dynamics in ML21, and the two values of the last two rows of H correspond to the lower and upper bounds of the estimations, as explained in the text.

	e max 0.35 38.793 39.037 H 4 38.549	H 6 58.91 59.16 58.66	H 8 56.90 57.68 56.11	H 10 52.95 53.74 52.16	H 49.67 50.47 48.87	L 4 40.289 40.558 40.020	L 6 47.20 47.44 46.96	LG09 49.22 51.54 46.90
	0.4 13.294 13.464 13.124	31.62 31.85 31.39	29.05 29.77 28.33	25.33 26.03 24.65	21.19 21.85 20.55	13.698 13.887 13.510	24.47 24.67 24.26	25.55 27.63 23.58
	0.5	0.483 0.519 0.450	5.12 5.24 5.02	5.04 5.39 4.70	2.95 3.23 2.70	2.04 2.27 1.82	0.730 0.778 0.685	3.85 3.94 3.76	3.96 4.97 3.15
	0.6	0.012 0.019 0.008	1.53 1.59 1.47	2.40 2.65 2.17	1.15 1.33 0.99	0.52 0.65 0.42	0.056 0.070 0.044	1.84 1.90 1.77	0.96 1.53 0.60
	0.7	0.001 0.004 0.000	1.06 1.11 1.01	2.07 2.31 1.86	1.08 1.26 0.93	0.45 0.57 0.35	0.010 0.017 0.006	1.41 1.47 1.36	0.92 1.48 0.57
	0.8	0.001 0.004 0.000	0.75 0.79 0.71	1.91 2.14 1.70	1.02 1.19 0.87	0.40 0.57 0.28	0.41 0.52 0.32	0.006 0.011 0.003	1.20 1.25 1.15	0.88 1.43 0.54
	0.9	0.000 0.003 0.000	0.42 0.45 0.38	1.78 2.00 1.58	1.00 1.17 0.85	0.14 0.26 0.08	0.40 0.51 0.31	0.000 0.003 0.000	0.85 0.89 0.80	0.84 1.38 0.51

Table 4 .

 4 2 -Rankings of Fourier harmonics. Partial contributions to g 1 (arcsec yr -1 ) from the

	harmonics at degree 4 (F i (4) ) and 6 (F i (6) ), along the unstable solution of H 6 of Fig. 4.2. The maximum filtered contributions are denoted by C k,ℓ (2p) = g k,ℓ 1(2p) (t

⋆ ), with t ⋆ = arg max t≤T |g k,ℓ 1(2p) (t)| (Eqs. (

  Table 2).

	i 1 2 3 4 5 6 7 g 1 -2g 2 + g 4 + s 1 -2s 2 + s 4 Fourier harmonic F i O i τ res i 4 12% 0.332 0.526 ∆ω i 0.093 g 3 -g 4 -s 3 + s 4 4 19% 0.302 0.611 0.154 g 1 -g 2 + s 1 -s 2 6 23% 0.105 0.223 0.041 g 2 -g 5 -2s 1 + 2s 2 6 70% 0.076 0.159 0.023 2g 3 -2g 4 -s 3 + s 4 4 10% 0.074 0.178 0.056 g 1 -g 5 -s 1 + s 2 4 6% 0.066 0.098 0.025 g 2 -g 4 + s 2 -s 4 8 5% 0.061 0.074 0.051 8 4 17% 0.056 0.090 0.028 g 1 -g 3 + s 2 -s 3 9 6 5% 0.053 0.061 0.037 g 1 + g 3 -2g 4 + s 2 -s 3 10 8 9% 0.052 0.140 0.007 3g 3 -3g 4 -s 3 + s 4 11 6 5% 0.038 0.047 0.028 g 2 -g 3 -s 1 + 2s 2 -s 3 12 6 36% 0.038 0.083 0.016 g 1 -2g 3 + g 4 + s 2 -s 4 13 6 5% 0.037 0.043 0.028 2g 1 -g 3 -g 5 + s 2 -s 4 14 6 2% 0.033 0.036 0.031 g 4 -g 5 -s 2 + 2s 3 -s 4 8 25% 0.033 0.045 0.014 15 g 1 -2g 3 + g 4 + s 1 + s 3 -2s 4 16 4 23% 0.032 0.054 0.017 g 1 -g 4 + s 1 -s 4 17 0.023 g 1 -2g 2 + g 5 + 3s 1 -3s 2 10 6% 0.032 0.039 18 4 18% 0.031 0.073 0.016 g 1 -g 4 + s 2 -s 3 8 2% 0.031 0.042 0.023 19 3g 1 -g 2 -g 4 -g 5 + s 1 -s 3 20 6 29% 0.028 0.051 0.016 2g 1 -g 2 -g 3 + s 1 -s 3 21 6 3% 0.026 0.028 0.021 2g 1 -g 2 -g 4 + s 1 -s 3 22 0.012 3g 3 -3g 4 -2s 3 + 2s 4 10 8% 0.025 0.055 8 3% 0.023 0.036 0.012 23 2g 1 -g 2 -2g 3 + g 4 + s 1 -s 4 24 6 16% 0.023 0.048 0.010 2g 3 -g 4 -g 5 -s 1 + s 4 25 8 7% 0.021 0.030 0.008 g 1 -3g 3 + 2g 4 + s 2 -s 4 26 6 6% 0.021 0.032 0.004 g 1 -g 2 -g 3 + g 4 + s 1 -s 2 27 6 3% 0.021 0.022 0.017 g 1 + g 3 -2g 4 + s 1 -s 4 28 0.006 g 1 + g 2 -2g 5 -3s 1 + 3s 2 10 4% 0.020 0.028 8 4% 0.020 0.027 0.008 29 3g 1 -g 2 -g 4 -g 5 + s 2 -s 3 30 2g 1 -g 4 -g 5 + s 2 -s 4 6 7% 0.020 0.029 0.007

Table 5 . 3 -

 53 Top of ranking R 3 . First 10 symmetry-breaking resonances of H 10 along the 5-Gyr nominal solution of Gauss's dynamics, that only involve g 5 among the external modes and change E 2n , C inc , and C 2 , respectively. For C inc , Only harmonics that are resonant for more than one percent of time are shown, i.e. τ res > 1% -3g 2 + 2g 5 -s 1 + 2s 2 -s 4 10 0.01 o / oo 1e-4 1e-4

	i	Fourier harmonic F i O i E 2n	τ res i	∆ω i
		No resonances
	1 4g 1 -g 2 -g 3 -s 1 -2s 2 + s 4 10 C inc 8 2g 1 + g 2 -g 4 -2s 1 -s 2 + s 3 3g 1 -g 2 -s 1 -3s 2 + s 3 + s 4 10 8 2g 1 -g 4 + g 5 -s 1 -2s 2 + s 4 4g 1 -g 2 -g 4 -s 1 -2s 2 + s 4 10 4g 1 -g 2 -g 4 -3s 2 + s 3 10 8 3g 1 -g 4 -2s 1 -s 2 + s 4 8 3g 1 -g 3 -2s 1 -s 2 + s 4 8 2g 1 -g 3 + g 5 -s 1 -2s 2 + s 4 8 10 g 1 + g 2 -g 3 + g 5 -2s 1 -s 2 + s 3 1 2 3 4 5 6 7 8 9 C 2 1 g 1 1e-4 2.0% 0.006 0.009 0.001 2.5% 0.003 0.008 0.002 1.5% 0.003 0.005 0.001 1.3% 0.003 0.006 0.001 1.4% 0.003 0.004 0.001 1.3% 0.003 0.005 0.001 1.8% 0.003 0.004 0.001 1.6% 0.002 0.006 0.000 1.0% 0.002 0.004 0.000 1.3% 0.002 0.003 0.001
	2 3 4 5 6 7 8 9 10	2e-5 2g 1 -4g 2 + 2g 5 + s 2 -s 4 10 0.08 o / oo 8e-5 1e-4 3e-5 3g 2 -3g 5 + s 1 -2s 2 + s 3 10 0.01 o / oo 3e-5 3e-5 2e-6 g 1 -4g 2 + g 3 + 2g 5 -s 1 + s 2 10 0.13 o / oo 3e-5 5e-5 2e-5 g 1 + 3g 2 -4g 5 -s 2 + s 4 10 0.56 o / oo 3e-5 4e-5 2e-5 g 1 -4g 2 + 3g 5 + s 2 -s 4 10 0.01 o / oo 2e-5 2e-5 7e-6 3g 2 -3g 5 + s 1 -2s 2 + s 4 10 0.13 o / oo 1e-5 3e-5 9e-6 4g 2 -g 3 -g 4 -2g 5 + s 1 -s 3 10 0.03 o / oo 1e-5 1e-5 9e-6 2g 1 -5g 2 + g 4 + 2g 5 10 0.02 o / oo 9e-6 9e-6 2g 1 -5g 2 + g 3 + 2g 5 10 0.14 o / oo 6e-6 8e-6 3e-6

  H • 2n A strong hierarchy of statistical variances among the QIs emerges from the size of the leading symmetry-breaking resonances in Table 5.3 and from the orbital solutions in Figs. 5.2, C.3 and C.4. One has E • 2n ≺ C 2 ≺ C inc . While E • 2n is an exact non-linear integral of motion, we expect that its linear truncation E • 2 = E 2 varies more than C 2 and C inc . Therefore, we consider the ordered set of QIs of motion {C 2 , C inc , E 2 } represented by the ordered set of vectors S ′ QIs = {γ 2 , γ 1 , γ 3 }.

  • 4 dynamics. The quasi-integral C 2 represents the slowest linear DOF: it coincides with the last principal component PC 8 , which has the smallest variance among all the linear combinations of the action variables. C inc and E 2 represent the second and the third slowest linear DOF, respectively: the component of C inc orthogonal to C 2 , i.e. C ⊥ 2 , matches the second last principal component PC 7 ; the component of E 2 orthogonal to the subspace generated by (C 2 , C inc ) matches the third last principal component PC 6 . The strong hierarchical structure of the slow variables for the simplified dynamics H • 4 is clearly confirmed by the almost frozen basis vectors of the PCs.

• 4

 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 •

Lyapunov time is the inverse maximal Lyapunov exponent, which is the coefficient characterizing the exponential separation of close trajectories.

δ 18 O is a measure of ratio between two oxygen isotope 18 O and 16 O.

E represents the exponential operator, j stands for the imaginary unit.

The definition of instability in[START_REF] Abbot | Rare event sampling improves mercury instability statistics[END_REF] is however different.

Throughout the paper, by activation we mean the exploration of the chaotic zone of the resonance, independently of the entrance in a libration state.

The dynamics truncated at degree 4 produces nevertheless the same chaos of the full system, as measured by the finite-time maximum Lyapunov exponent(Mogavero and Laskar,

2022).2 This resonance involves the fundamental precession frequencies of Mercury and Jupiter perihelia[START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF][START_REF] Batygin | On the Dynamical Stability of the Solar System[END_REF][START_REF] Boué | A simple model of the chaotic eccentricity of Mercury[END_REF].3 "At each stage of its evolution, the system should have a time of stability comparable with its age"[START_REF] Laskar | Large Scale Chaos and Marginal Stability in the Solar System[END_REF]

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

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.

In this regard, it should be noted that a relation between QIs and Lyapunov exponents has been already highlighted in simple systems by some authors[START_REF] Contopoulos | On the number of isolating integrals in Hamiltonian systems[END_REF][START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF]

Faint young Sun paradox is the contradiction between the geological evidence of warm, wet early Earth, and the low-luminosity Sun in the same epoch from astrophysical models (seeFeulner,

2012).

apart from the important relativistic corrections.

thank him for all the passionate discussions and the enjoyable companionship throughout all the ups and downs of this journey. I also want to thank the CFM Foundation for Research for financing my PhD scholarship.

where σ is the standard deviation of the sample, IQR is its interquartile range, and β is a constant of choice. The bandwidth with β = 1/5 corresponds to the Silverman's rule of thumb [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. The BW β=1/5 is a version of the optimal choice of bandwidth for Gaussian distributed data that is slightly modified for better adaption to non-Gaussian data. The bias error and variance error of the KDE with this bandwidth will be on the same order of magnitude. Under-smoothing, that is, choosing a smaller bandwidth, shrinks the bias so that the total error is dominated by the variance error, which can then be estimated by the bootstrap method (Hall et al., 1995); the common value of β for under-smoothing is 1/3. When the sample is identically distributed but correlated, KDE is still valid under some mixing conditions [START_REF] Robinson | Nonparametric estimators for time series[END_REF][START_REF] Hart | Some automated methods of smoothing time-dependent data[END_REF]. Indeed, in the case of observations that are not too highly correlated, the dependence among the data that fall in the support of the kernel function K can actually be much weaker than it is among the entire sample. This principle is known as "whitening by windowing" [START_REF] Hart | Some automated methods of smoothing time-dependent data[END_REF]. Therefore, the correlation in the sample does not invalidate the use of KDE and only impacts the variability of the estimation (see Sect. 3.3.2). With regard to the choice of bandwidth, Hall et al., 1995 suggested that using the asymptotically optimal bandwidth for the independent data is still a good choice, even for some strongly dependent data sequences. The samples generated by a chaotic measure-preserving dynamical system (as in the case of the numerical integration of the Solar System) resemble those of mixing stochastic processes; therefore, the theory of KDE should also be applicable for dynamical systems, although the formulation might be different [START_REF] Bosq | Nonparametric estimation of the chaotic function and the invariant measure of a dynamical system[END_REF][START_REF] Maume-Deschamps | Exponential inequalities and functional estimations for weak dependent data: applications to dynamical systems[END_REF][START_REF] Hang | Kernel density estimation for dynamical systems[END_REF]. is primarily governed by the long-range interactions characterized by harmonics involving Saturn (g 6 , s 6 ) and Uranus (g 7 , s 7 ) secular frequencies, the latter is strongly influenced by the harmonics at degree 6 2 . Both factors are determined by and thus sensitive to these interactions, which were weak and hidden before our work. Yet, they are crucial for the dynamical stability of the system.

Appendix B

Long-term instability of the inner Solar System: numerical experiments

This section assembles the appendices for the chapter 4 from [START_REF] Hoang | Long-term instability of the inner Solar System: numerical experiments[END_REF].

B.1 Secular dynamics at second order in planetary masses

We use the secular equations of motions of [START_REF] Laskar | Accurate methods in general planetary theory[END_REF](Laskar, , 1990[START_REF] Laskar | Chaotic diffusion in the Solar System[END_REF], and references therein). They were obtained via series expansions in planetary masses, eccentricities, and inclinations, as well as through second-order analytical averaging over the rapidly changing mean longitudes of the planets. The expansion was truncated at the second order with respect to the masses and to degree 5 in eccentricities and inclinations. The equations include corrections from general relativity and Earth-Moon gravitational interaction. This leads to the following system of ordinary differential equations, denoted by L 6 throughout this paper:

where

The planets are indexed in order of increasing semi-major axis, as usual. The variable ̟ k is the longitude of the perihelion, Ω k is the longitude of the ascending node, e k is eccentricity, and i k is inclination. The function Φ 3 (ω, ω) and Φ 5 (ω, ω) are the terms of degree 3 and 5, respectively. The 16 × 16 matrix Γ is the linear Laplace-Lagrange system, which is slightly modified to make up for the higher-order terms in the outer Solar System.

To mimic H 4 , we define the new model L 4 by dropping the terms of degree 5 from the equations of the inner planets, that is:

Appendix C

Timescales of chaos in the Inner Solar System

This section assembles the appendices for the chapter 5.

C.1 Lyapunov spectrum

Convergence. We perform two tests to address the convergence of our implementation of the [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF] 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 C.1a 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 [START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF] algorithm looses memory of the initial tangent vectors and purely retrieves the FT-LCEs as defined in Eq. (5.10). Therefore, Fig. 5.1a shows asymptotically the dependence of the FT-LCEs on the initial condition z 0 and represents their statistical distribution over the phase-space domain explored by the dynamics in a non-ergodic way. The convergence of the computation is clearly slower for smaller exponents, but a comparison with Fig. 5.1a indicates that, even in the case of λ 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. (5.9). For a single orbital solution, the relative numerical error on each exponent λ i can be estimated as

We plot in The relative errors decrease asymptotically with time, as expected. Even in the case of the smallest exponent, λ 8 , the median error is less than 10% at 10 Gyr.