Long-term dynamics of the inner planets in the Solar System

Although the discovery of the chaotic motion of the inner planets in the Solar System (Mercury to Mars) was made more than 30years ago, the secular chaos of their orbits still demands more analytical analyses. In addition to the high-dimensional structure of the motion, this is probably related to the lack of an adequately simple dynamical model. We consider a new secular dynamics for the inner planets, with the aim of retaining a fundamental set of interactions that explains their chaotic behaviour and at the same time is consistent with the predictions of the most precise orbital solutions currently available. We exploit the regularity in the secular motion of the outer planets (Jupiter to Neptune) to predetermine a quasi-periodic solution for their orbits. This reduces the secular phase space to the degrees of freedom dominated by the inner planets. In addition, the low masses of the inner planets and the absence of strong mean- motion resonances permits us to restrict ourselves to ﬁrst-order secular averaging. The resulting dynamics can be integrated numerically in a very efﬁcient way through Gauss’s method, while computer algebra allows an analytical inspection of planet interactions when the Hamiltonian is truncated at a given total degree in eccentricities and inclinations. The new model matches reference orbital solutions of the Solar System over timescales shorter than or comparable to the Lyapunov time very satisfactorily. It correctly reproduces the maximum Lyapunov exponent of the inner system and the statistics of the high eccentricities of Mercury over the next ﬁve billion years. The destabilizing role of the g 1 − g secular resonance also arises. A numerical experiment, consisting of a thousand orbital solutions over one hundred billion years, reveals the essential properties of the stochastic process driving the destabilization of the inner Solar System and clariﬁes its current metastable state.


Introduction
The inner solar system, with the secular chaos of its planetary orbits (Laskar 1989(Laskar , 1990cSussman & Wisdom 1992), holds a special status among the dynamical systems of celestial mechanics. Even if its chaotic behaviour has been demonstrated numerically three decades ago, still no analytical study corroborates or rules out the role in chaos generation of the secular resonances 2(g 3 − g 4 ) − (s 3 − s 4 ) and (g 3 − g 4 ) − (s 3 − s 4 ), between the Earth and Mars fundamental precession frequencies, as proposed in (Laskar 1990c(Laskar , 1992 and supported by (Laskar et al. 2004). Such a fact probably relates to the complex network of dynamical interactions among the inner planets: the fundamental precession frequencies of their orbits vary in an intricate way over a 0.1 yr −1 scale in a few tens of million years (Myr), with the exception of the Venus-dominated eccentricity mode g 2 which has somewhat smaller variations (Laskar 1990c;Laskar et al. 2004). This behaviour reveals the essential high-dimensional structure of the inner solar system, which limits the possibility of faithfully modelling its dynamics with few effective degrees of freedom (e.g., Lithwick & Wu 2011;Batygin et al. 2015).
Getting analytical insight into the motion of the inner planets requires an appropriate dynamical modelling of the long-term evolution of their orbits. On the one hand, such a model must be consistent with the predictions of the reference numerical integrations available in literature (Laskar 1990c;Laskar et al. 2004;Laskar 2008;Laskar & Gastineau 2009), to ensure it reproduces the dynamical features of the inner system with suf-ficient precision. On the other hand, the corresponding Hamiltonian should be put in a form suitable for the systematic application of canonical perturbation techniques (Hori 1966;Deprit 1969), which is essential to an unbiased analysis of such a high-dimensional dynamics. Moreover, the possibility to numerically integrate the equations of motion in an efficient way is fundamental to study the chaotic evolution of the orbital solutions in a statistical way. Unfortunately, the construction of such a model turns out to be a delicate task. In principle, one could just consider the full N-body Hamiltonian of the Newtonian gravitational interactions among the solar system planets, with the addition of the leading corrections coming from general relativity and the Earth-Moon interaction. This already reproduces the precession frequencies of the inner orbits with a precision better than 0.01 yr −1 (Laskar 1999). However, such a Hamiltonian is unnecessarily complicated, as it includes short-time harmonics, with periods of less than 5 000 yr (e.g., Carpino et al. 1987), which are known to generate small quasi-periodic oscillations in the inner orbits, without being implied in chaos generation. Indeed, the inner planets are not involved in any relevant meanmotion resonance. At the same time, the long-term numerical integration of the corresponding equations of motion is very time consuming, the resulting solutions needing to be filtered to extract the secular trend of the orbits (Carpino et al. 1987;Nobili et al. 1989). These facts suggest to consider a secular Hamiltonian to directly describe the slow movement of the planet perihelia and nodes, after proper averaging over the short-time orbital motion (Laskar 1984(Laskar , 1985. Secular dynamics includes the es-Article number, page 1 of 25 arXiv:2105.14976v1 [astro-ph.EP] 31 May 2021 sential planet interactions responsible for chaos in the inner solar system and allows to perform the fastest long-term numerical integrations (Laskar 1988(Laskar , 1989(Laskar , 1994(Laskar , 2008. Indeed, this made the discovery of chaos in the inner system possible before the use of symplectic integration schemes (Sussman & Wisdom 1992). Unfortunately, an effective secular model for the entire solar system has to be of high order in planet masses, principally because of the 5:2 near mean-motion resonance between Jupiter and Saturn, the so-called Great Inequality (Laplace 1785;Laskar 1996). A simple averaging of the N-body Hamiltonian over the planet mean longitudes, resulting in a first-order secular dynamics in planet masses, would reproduce very poorly the fundamental frequencies g 5 and g 6 , which dominate the perihelion precession of Jupiter and Saturn, respectively (Laskar 1988) 1 . The construction of such higher-order models requires the manipulation of large Poisson series and the use of sophisticated computer algebra systems. This is probably the reason why they are still not widely used, at least as a basis of extensive research.
This paper introduces a new secular dynamics for the orbits of the inner planets. It is based on the practical long-term regularity of the outer planet trajectories (Laskar 1990c;Laskar et al. 2004;Hoang et al. 2021), the smallness of the inner planet masses and the absence of relevant mean-motion resonances in the inner system (Sect. 2). We show that the present model can be numerically integrated via the so-called Gauss (1818)'s method (Sect. 3), while its Hamiltonian is suitable for a systematic expansion in planet eccentricities and inclinations, once one employs the algorithm of (Laskar & Robutel 1995) and a computer algebra system like TRIP  (Sect. 4). We compare its orbital solution to a reference N-body integration over short (secular) timescales (Sect. 5) and compute the corresponding maximal Lyapunov exponent in a statistical fashion (Sect. 6). We then determine the percentages of the high Mercury eccentricities over the next 5 billion years and highlight the destabilizing role of the g 1 − g 5 secular resonance (Sect. 7). Finally, we perform a new numerical simulation involving a thousand orbital solutions over one hundred billion years, to characterize the effective stochastic process which drives the destabilization of the inner solar system, and discuss its implications on the conjecture of marginal stability of a secularly evolving planetary system formulated in (Laskar 1996) (Sect. 8).

Dynamical model
We model the dynamics of the largest bodies in the solar system, by considering the Sun and the N = 8 planets as point masses m 0 , (m k ) k=1,N , indexed in order of increasing semi-major axis. The barycentric coordinates of the bodies and the corresponding momenta are denoted by (R 0 ,R 0 = m 0Ṙ0 ), (R k ,R k = m kṘk ) k=1,N . By employing the canonical heliocentric variables of (Poincaré 1896), the Hamiltonian of the Newtonian gravitational interactions among the bodies reads where r k = R k − R 0 are the planet heliocentric coordinates and r k =R k their conjugated momenta, µ k = m 0 m k /(m 0 + m k ) are the reduced masses of the planets and G is the gravitational constant (Laskar 1991;Laskar & Robutel 1995). The Hamiltonian H is a perturbation to the union of disjoint Kepler problems, so that it is useful to introduce a set of canonical variables which trivially integrates the unperturbed problems (e.g., Morbidelli 2002). By adopting from now on the momentumcoordinate ordering of conjugate pairs, appropriate variables are (Λ k , λ k ; x k , − jx k ; y k , − jȳ k ) k=1,N , defined as Λ k = µ k G(m 0 + m k )a k , where a k are the planet semi-major axes, λ k the mean longitudes, e k the eccentricities, i k the inclinations, k the longitudes of the perihelia and Ω k the longitudes of the nodes (Laskar 1991;Laskar & Robutel 1995). Throughout the paper, j = √ −1 stands for the imaginary unit, E represents the exponential operator and the overbar denotes the conjugate of a complex variable. The variables x k and y k are the Poincaré's rectangular coordinates in complex form; we shall refer to them throughout the paper as Poincaré's complex variables, or simply Poincaré's variables. With such a choice of canonical variables, the integrable part of the Hamiltonian (1) reads so that the Poincaré's complex variables are constants of motion for the Kepler problem.
In the regime of small orbital eccentricities and inclinations, which characterizes the solar system, Eqs. (2) give x k = √ Λ k /2e k E j k + O(e 3 k ) and y k = √ Λ k /2i k E jΩ k + O(e 2 k i k , i 3 k ), and the Poincaré's variables are also small. The principal part of the two-body perturbation in Eq. (1) can be thus expanded as a Fourier series in the planet mean longitudes, with polynomial coefficients depending on the Poincaré's variables (e.g., Laskar 1990bLaskar , 1991Laskar & Robutel 1995), where N = (n, n ,n,n , m, m ,m,m ) is a tuple of non-negative integers 2 . Following (Laskar & Robutel 1995), we have defined the dimensionless Poincaré's variables 2Λ and the semi-major axis ratio α = a/a , with a < a . The analytical expression of the coefficients Γ , N (α), only depending on the semi-major axis ratio, is given in (Laskar & Robutel 1995) in terms of Laplace coefficients. The indirect part of the twobody perturbation can also be expanded as a Fourier series in the mean longitudes, The computation of the Fourier coefficients T , is outlined in Appendix A. Equations (4) and (5) allow to explicitly compute the Fourier expansion of the Hamiltonian perturbing function, where λ stands for the vector of the planet mean longitudes, λ = (λ 1 , . . . , λ 8 ), and the coefficients H depend on all the remaining canonical variables.

Secular dynamics
The long-term dynamics of the solar system planets, in particular that of the inner ones, essentially consists of the slow precession of their perihelia and nodes, driven by secular, i.e. orbitaveraged, gravitational interactions (Laskar 1990c;Laskar et al. 2004). A secular Hamiltonian, describing such long-term motion, can be introduced in its simplest form by searching for a change of variables that eliminate, at first order in the planet masses, all the harmonics with non-null wavevectors in the Fourier expansion of the perturbation (6). In canonical perturbation theory (Hori 1966;Deprit 1969;Morbidelli 2002), this elimination is achieved through a canonical transformation defined as the time-1 flow of a generating Hamiltonian S satisfying the homologic equation where the braces represent the Poisson bracket. The anglebracket operator stands for averaging over the mean longitudes, with the integration defined over the hypertorus T N , at fixed values of all the remaining canonical variables. This means that H 1 is the Fourier coefficient H 0 corresponding to the null harmonic in the expansion (6). The homologic equation (7) gives the generating Hamiltonian S as a formal Fourier series, where n = ∂H 0 /∂Λ is the vector of the planet mean motions. The secular Hamiltonian H is formally given by the Lie transform generated by the function S and applied to the Hamiltonian H, where L S · = {S , ·} is the Lie derivative associated to the generating Hamiltonian, L 0 S is defined as the identity operator and L n S · = L S L n−1 S · for n ≥ 1. The Hamiltonian H in Eq. (10) is expressed in the new canonical variables (Λ k ,λ k ;x k , − jx k ; y k , − jŷ k ) k=1,N , which we shall call the secular variables. They are related to the original ones via the Lie transforms The original variables are therefore the superposition of the secular ones and short-time oscillations generated by the n ≥ 1 terms of the Lie transforms. Since in the present study we are interested in the long-term dynamics of the planets, we shall only focus on the secular variables. Therefore, to keep simpler notation, we shall omit the hat on the secular variables from now on. Differently from the outer planets, the inner ones are not currently involved in strong mean-motion resonances. More precisely, a maximal contribution of only 0.07 yr −1 to the fundamental precession frequencies of the inner orbits in the Laplace-Lagrange solution arises from the second order in the planet masses (Laskar 1985, Table 8). This same contribution is 0.9 yr −1 for the outer planets, mainly due to the Great Inequality. Building on this fact, in the present work we choose to truncate the series (10) at first order in the planet masses, i.e. we neglect quadratic and higher-order terms with respect to the Fourier coefficients H . Indeed, we expect the main contribution to the precession frequencies of the inner orbits to come from the linear terms, given the small masses of the inner planets. In absence of strong mean-motion resonances, in Eq. (9) the denominators · n involving at least one inner planet are sufficiently far from zero. Under these assumptions, the higher-order terms only produce small corrections to the dynamics generated by the leading ones. Using the homologic equation (7), one thus obtains The resulting secular Hamiltonian H is simply the average of the N-body Hamiltonian (1) over the planet mean longitudes. The averaging process is mathematically equivalent to replacing each planet by its instantaneous Keplerian orbit, with the corresponding mass distributed along it in a way that is inversely proportional to the local orbital speed of the planet. The secular dynamics is thus the slow gravitational interaction of such Keplerian rings. This equivalence was pointed out by Gauss (1818) and it is thus referred to as Gauss' dynamics. As a result of the averaging over the mean longitudes, the Λ k variables are constant of motion in the secular dynamics, and so are the semi-major axes of the planets. The resulting Hamiltonian system consists of two degrees of freedom for each planet 3 , corresponding to the two Poincaré's complex variables x k and y k .
General relativity and minor effects. The dynamical interactions described by the Hamiltonian (1) are not sufficient to finely reproduce the precession frequencies of the inner orbits; to this end, additional physical effects must be taken into account (Laskar 1999). Indeed, it is well-known that general relativity contributes with 0.430 yr −1 to the secular precession of Mercury perihelion. Such a correction is critical for the statistics of the long-term destabilization of the inner orbits (Laskar 2008;Laskar & Gastineau 2009), since it moves the system away from the g 1 − g 5 secular resonance, which is responsible for the very high eccentricities of Mercury (Laskar 2008;Batygin & Laughlin 2008;Boué et al. 2012). We shall therefore include in the Hamiltonian (12) the leading secular contribution of general relativity, which reads where c is the speed of light (e.g., Saha & Tremaine 1992;Mogavero 2017). We point out that the first term in the summation only depends on the semi-major axis a k and is thus a constant quantity in the secular dynamics. The next largest effect to be taken into account would be the Earth-Moon gravitational interaction, accounting for a 0.077 yr −1 contribution to the secular perihelion precession of the Earth (Laskar 1999). However, this is the order of magnitude of the contribution coming from the Hamiltonian terms at second order in the planet masses, which are neglected in the present model. Indeed, at degree 2 in planet eccentricities and inclinations, (Laskar 1985, Table 8) reported contributions of 0.073 yr −1 and 0.062 yr −1 to the Earthdominated and Mars-dominated eccentricity modes g 3 and g 4 , respectively, the other frequency corrections being at least ten times smaller. Since the denominators · n appearing in Eq. (9) and involving the inner planets are not close to zero, we expect the contributions from higher degrees to be generally smaller than those from degree 2. We thus choose to not include the Earth-Moon interaction in the present model, preferring to deal with a Hamiltonian whose dependence on eccentricities and inclinations is exact, as compared to a non-decisive increase of precision on the orbit precession frequencies. Moreover, we know that the N-body Hamiltonian (1), corrected for general relativity, reproduces the maximum Lyapunov exponent of the inner solar system (Rein & Tamayo 2018). Therefore, apart from irrelevant constant terms only depending on the semi-major axes in Eqs. (3) and (13), the secular Hamiltonian of the 8 solar system planets considered in the present study reads

Forced inner planets
As discussed in the Introduction, the secular Hamiltonian (14) would not correctly reproduce the frequencies of the Jupiter and Saturn dominated eccentricity modes g 5 and g 6 , respectively. Fortunately enough, the very small variations over billions of years of the precession frequencies of the outer orbits, compared to those of the inner one (Laskar 1990c;Laskar et al. 2004;Hoang et al. 2021), naturally suggest a way to make out of the Hamiltonian (14) a very effective model for the inner system. This is achieved by choosing, once for all, an explicit quasiperiodic time dependence for the orbits of the giant planets, i.e. by expressing the corresponding Poincaré's complex variables as finite Fourier series, where t denotes the time,x k andỹ k are complex amplitudes, m k and n k integer vectors, and φ(t) = ω o t, with ω o the vector of the (constant) fundamental precession frequencies of the outer orbits, denoted as ω o = (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ) 4 (Laskar 1990c). The number of harmonics M k , N k appearing in the decompositions depends on the planet. When such predetermined 4 The frequency of the Jupiter-dominated inclination mode s 5 is null because the total angular momentum of the system is conserved. time dependence is injected in Eq. (14), one obtains the Hamiltonian H of a forced secular inner solar system, The explicit time dependence in the Hamiltonian H physically means that the inner planets interact with each other while moving in an external gravitational potential generated by the giant planets 5 . The inner planets thus constitute an open system and the corresponding dynamics does not possess any fundamental integral of motion, such as the energy or angular momentum. As a result of the predetermination of the outer orbits, in addition to the explicit time dependence, the Hamiltonian H possesses 8 degrees of freedom. As already stated, compared to the secondorder secular system of (Laskar 1984(Laskar , 1985, the present model neglects terms of order higher than one in the planet masses. Nevertheless, the precession frequencies of the outer orbits can be set to very precise values in Eq. (15). Moreover, the Hamiltonian (16) is not truncated neither in eccentricities nor in inclinations, so that the range of validity of its dynamics extends to very excited states, provided that the Keplerian orbits of the planets do not cross each other. This makes the present model perfectly suited to explore the very long-term evolution of the inner system, even when high eccentric and inclined orbits become statistically recurrent.

Construction of the outer planet solution. Initial conditions
The quasi-periodic form of the outer orbits in Eq. (15) is established numerically and explicitly reported in Appendix D. The full equations of motion of the main bodies of the solar system are numerically integrated over 30 Myr in the future, following the comprehensive model of ). The initial conditions of the integration were adjusted through least squares to the high precision planetary ephemeris INPOP13b Fienga et al. 2014) extending over 1 Myr. Throughout the paper, we shall refer to this direct numerical integration as LaX13b. Through frequency analysis (Laskar 1988;Laskar et al. 1992;Laskar 1993Laskar , 2005 of the orbital solution, the leading secular harmonics of the dimensionless Poincaré's variables of the outer planets (X k , Y k ) k=5,8 are extracted to construct the Fourier series (15), as in (Laskar 1988(Laskar , 1990c. This is performed up to a numerical precision that does not allow anymore to recognize in an unambiguous way new harmonics as linear combination of the fundamental frequencies ω o . The precision of the outer planet solution is shown in Table D.4 of Appendix D. We report there the root mean square of the dimensionless Poincaré's variables of the outer planets in the LaX13b solution, and that of the corresponding residuals after subtraction of their quasi-periodic decomposition from Eq. (15) and filtering out of the non-secular Fourier components. The effectiveness of the quasi-periodic approximation is illustrated in Figs. D.1 and D.2. The choice of the secular semi-major axes and that of the initial conditions for the secular Poincaré's variables of the inner planets deserve a discussion. In principle, the initial secular variables have to be computed by inverting the Lie transforms in Eq. (11). They should not be set to the initial value of the corresponding original variable, as this could cause an offset in the secular frequencies of the motion, due to the short-time (highfrequency) oscillations generated by the Lie transform. This is Table 1: Average wall-clock time on an Intel(R) Core(TM) i7-7700T CPU @ 2.90GHz for the numerical integration of the forced inner system according to the total degree of truncation of the Hamiltonian in eccentricities and inclinations (H 2n is given in Eq. (25), while H in Eq. (16)).
H 2n 4 6 8 10 H 100-Myr orbital solution 1.3 s 7.5 s 35 s 4 m 01 s 2 m 48 s a general rule when constructing averaged dynamical systems (e.g., Laskar & Simon 1988). To avoid the explicit computation of the generating function S in Eq. (9) and the related Lie transforms, appropriate initial conditions for the secular solar system can be effectively computed by filtering out the short-time components of the numerical solution of the non-averaged system, as it is done in the present study. The constant term in the frequency analysis of the variables Λ k , performed on the LaX13b solution, provides the value of the secular semi-major axes of the planets (

Numerical integration. Gauss's method in Hamiltonian formalism
The dimensionless Hamilton's equations for the forced inner system (16) reaḋ for k ∈ {1, 2, 3, 4}. Because of the averaging over the mean longitudes in Eq. (14), the above equations include double integrals of the form where Z stands for X or Y, alternatively. The equations of motion (17) are thus non-algebraic, and the numerical computation of the derivatives is in principle much more complex than in Nbody dynamics or in polynomial secular systems as in (Laskar 1985(Laskar , 1990c. Nevertheless, it is well-known that the Eqs. (17) can be integrated numerically in a very efficient way through the so-called Gauss's method (Gauss 1818; Bour 1855; Hill 1882).
Indeed, the vector f in Eq. (19) is proportional to the gravitational force exerted on a test particle by a Keplerian ring, and can be analytically expressed in terms of complete Legendre's elliptic integrals (e.g., Olver et al. 2020, Chapter 19). This was first shown by Gauss (1818), who introduced at same time the arithmetic-geometric mean to numerically evaluate such special functions in a few elementary iterations and with high precision. Building on its modern derivation in (Musen 1970), we implemented Gauss's method into the present complex Hamiltonian formalism (17). This allows, in particular, to eliminate the degeneracy typically appearing at e = 0 (circular orbits) and i = 0 (equatorial orbits) (e.g., Fouvry et al. 2020), which is fundamental to guarantee high numerical precision. We employ a state-of-the-art algorithm, based on piecewise minimax rational function approximation, to numerically compute the complete elliptic integrals in double floating-point precision at the cost of elementary functions (Fukushima 2015). The remaining simple integral in Eq. (18) is then effectively evaluated via the trapezoidal rule, which converges exponentially fast with the number of function evaluations, because of the periodicity of the integrand (e.g., Touma et al. 2009). We apply the trapezoidal rule in an adaptive way, by doubling the number of function evaluations until the estimated relative error on the integral is smaller than 10 −12 . It should be noted that, as the convergence of the numerical integral is exponential, the resulting error is often orders of magnitude smaller than this tolerance. We integrate the Eqs. (17) using an Adams PECE method of order 12 and a conservative timestep of 250 years, as in (Laskar 1994). In absence of any fundamental integral of motion (see Sect. 2.2), we estimate the integration error after a time T following (Laskar 1994). Starting from the nominal initial conditions, we integrate the dynamics over the time interval [0, T /2] and then backwards to the initial time. By denoting the deviations from the initial coordinates of the system in the phase space as (δX k (T ), δY k (T )) k=1,4 , we define the relative integration error as where (X k (0), Y k (0)) k=1,4 are the initial conditions of the system. Over the time interval T ∈ [0, 20] Myr, the average integration error is well described by the power law δ(T ) = δ 0 [1 + (T/T 0 ) α ], with δ 0 = 9·10 −15 , T 0 = 0.6 Myr and α = 1.2. Therefore, one has δ(10 Myr) ≈ 3 · 10 −13 , and δ(10 Gyr) ≈ 10 −9 if one extrapolates by ignoring the chaotic behaviour of the solution. The integration error thus turns out to be similar to that of the secular system in (Laskar 1994). We emphasize that the present fixed-timestep integration scheme allows to reach very highly excited orbital states. Indeed, violations of the precision goal on the integral (18) typically occurs only when the system is already close to the first intersection of the instantaneous Keplerian orbits of two planets (see Sect. 7), and are related to an intrinsic discontinuity in the equations of motion (17) at orbit crossing (e.g., Touma et al. 2009). From a computational perspective, the choice of a multistep method is well adapted to the present dynamical system, as the derivative evaluation constitutes by far the most expensive step of the integration scheme. Once implemented in the C programming language, Gauss's method allows to compute 100-Myr orbital solutions of the forced inner planets in a few minutes on a PC, as shown in Table 1. Therefore, the computational cost turns out to be comparable to that of the secular system of Laskar (1990c), in spite of the different complexity of the corresponding Hamiltonians.
The rather involved analytical derivations adapting Gauss's method to the present complex Hamiltonian formalism shall be the subject of a forthcoming paper (Mogavero 2021), along with the release of the C program implementing it.

Analytical expansion. Computer algebra
In addition to its fast numerical integration, the present dynamical model allows for a systematic analytical development of its Hamiltonian by means of computer algebra. In this study we employ TRIP, a computer algebra system dedicated to perturbation series, specially those of celestial mechanics (Laskar 1990a;. The main objects of its symbolic kernel are the Poisson series, i.e. multivariate Fourier series whose coefficients are multivariate Laurent series, where (z p ) p=1,n and (ϕ p ) p=1,m are complex and real variables, respectively, k = (k p ) p=1,n ∈ Z n , = ( p ) p=1,m ∈ Z m and C k, are complex coefficients. Following (Laskar & Robutel 1995), the smallness of eccentricities and inclinations in the solar system can be exploited to develop the averages appearing in Eq. (14) as formal series in The rotational invariance of the secular Hamiltonian (14) requires n + n + m + m −n −n −m −m = 0, while from the planar symmetry it follows that m + m +m +m is an even integer. These relations imply that the monomials in Eq. (22) are even with respect to both (x, x ,x,x ) and (y, y ,ȳ,ȳ ) variables. The relativistic terms in Eq. (14) are readily expanded, By truncating the series (22) and (23) at a given total degree 2n (n ∈ N 0 ) in Poincaré's variables (x k ,x k , y k ,ȳ k ) k=1,8 , and after substitution in Eq. (14), one obtains a polynomial secular Hamiltonian for the ensemble of the solar system planets, where H (2p) groups all the monomials of same total degree 2p. Practically, the expansion (24) is readily obtained in TRIP, which implements the algorithm of (Laskar & Robutel 1995) to compute a truncation of the series (22). The truncated Hamiltonian (24) is in the form of a Poisson series (21) in Poincaré's complex variables, with no angular dependencies. The crucial point here is that its coefficients only depend on the planet masses and secular semi-major axes, which are constant parameters in the present dynamics. They can be thus numerically evaluated once for all, to obtain very compact series. This is shown in Table 2, which gives the total number of monomials in (24) with respect to the degree of truncation, along with the disk usage of the corresponding series when stored in a plain-text file, and their typical computation time in TRIP. The expansion (24) straightforwardly provides a truncated Hamiltonian for the forced inner system (16), The computation of the equations of motion for the truncated Hamiltonian (25) only requires to take the derivatives of a multivariate polynomial, and can be thus systematically performed in TRIP. We put the resulting polynomials in Horner form, to ensure speed and stability of their numerical evaluation, and wrap them into a C code to achieve the best computational performance 6 . The equations of motion are then integrated through the Adams PECE method of order 12 employed for Gauss's dynamics, with the same timestep of 250 years. The typical computational cost of a 100-Myr orbital solution is shown in Table 1, according to the degree of truncation. In light of the pertinence of the forced inner system that will be shown in Sect. 5, the truncated Hamiltonian (25), at degrees 4 and 6 in particular, constitutes the state of the art of very fast, still realistic, dynamical models of the inner solar system.

Forced Laplace-Lagrange dynamics
For small eccentricities and inclinations, the secular Hamiltonian (14) is a perturbation to the integrable Laplace-Lagrange (LL) problem. Indeed, its truncation at degree 2 in Poincaré's variables gives the quadratic form where we have defined the column vectors x = (x 1 , . . . , x 8 ) T and y = (y 1 , . . . , y 8 ) T , T is the transposition operator and the dagger stands for Hermitian transposition, i.e. x † =x T . The matrices M and N are real and symmetric. Since the variables x and y are not coupled in the LL Hamiltonian, the following derivations focus on the degrees of freedom related to x. A similar treatment holds for the y variables. After substitution of the time dependence of the outer planet variables given in Eq. (15), one obtains where we have defined the inner and outer planet column vectors, x 6 (t), x 7 (t), x 8 (t)) T , respectively, and the real matrix M has been written as a block matrix, with M ii and M oo symmetric 4 × 4 matrices. By discarding terms only depending on time, the Hamiltonian reads where Re stands for the real part of a complex quantity. The Hamiltonian (28) corresponds to a forced Laplace-Lagrange system. Since it is quadratic in the inner vector x i , one can analytically solve the corresponding equations of motion and introduce appropriate action-angle variables.

Solution of the equations of motion
The matrix M ii is real and symmetric and can thus be diagonalized through an orthogonal matrix O M , The columns of O M are the eigenvectors of M ii , while the diagonal entries of D M = −diag( g LL ) are the corresponding real  (Eqs. 42,44,48) The transformed Hamiltonian reads The corresponding Hamilton's equations are given bẏ and constitute a first-order inhomogeneous matrix ordinary differential equation. The general solution can be written as The free solution x i,F (t) is the general integral to the associated homogeneous equationẋ i = − j D M x i , representing the autonomous perihelia precession of the inner orbits, while the forced solution x i,f (t) is a particular integral of the complete Eq. (32), which arises from the gravitational forcing of the outer planets. We define, once for all, the forced solution to be: By employing the decomposition given in Eq. (15), the components of the forced solution are The (constant) denominators appearing in Eq. (35) are far from zero, as the inner planets are not involved in the corresponding secular resonances. The forced solution is thus well defined. 7 The minus sign in the definition of the matrix D M is such that the precession frequencies of the perihelia˙ have the same sign as g LL . 8 In the present notation, the canonical momenta constitute column vectors, while the coordinates are wrapped in row vectors.

Proper modes
The following derivation shows that there exists a canonical transformation depending on time, As in the case of an autonomous Laplace-Lagrange system, the new canonical variables u will physically correspond to the free part x i,F (t) of the solution in Eq. (33) (e.g., Morbidelli 2002). We thus begin by defining a new set of variables u such that By using the fact that x i,f (t) is a solution to Eq. (32), and discarding terms only depending on time, the forced LL Hamiltonian in Eq. (31) can be written as We now ask that the change of variables ( so that the last two terms in Eq. (37) shall be killed in the transformed Hamiltonian. By integrating the above equation with respect to time, one obtains where f (x i , u † ) is an unknown function, which has to be determined. Since the generating function F depends on the old momenta x i and the new coordinates − ju † , it must verify the relations (e.g., Landau & Lifshitz 1969) which imply With the choice f (x i , u † ) = ju † x i , these relations are both equivalent to Eq. (36).
Article number, page 7 of 25 A&A proofs: manuscript no. iss The previous derivations imply that the variable transforma- We call the complex variables u the proper modes of the Poincaré's variables x i . In the forced Laplace-Lagrange dynamics, they simply rotate in the complex plane, with constant angular frequencies given by g LL . According to Eq. (42), the corresponding Poincaré's variables thus result in a superposition of such independent harmonic oscillations and those arising from the forcing of the outer planets. When higher-degree terms of the Hamiltonian (25) are taken into account, the dynamics of the proper modes becomes coupled, and their frequency spectrum is no more monochromatic. Physically, the proper modes u correspond to the variables (z * k ) k=1,4 defined in (Laskar 1990c). Nevertheless, it is important to emphasize that their mathematical definitions do not coincide. Indeed, the transformation matrix O M do not include contributions at the second order in planet masses, as that considered in (Laskar 1990c) does. Moreover, the transformation (42) is nonlinear in the forcing of the outer planets, in the sense that harmonics of order 9 |m| greater than one are present in the decomposition (15) (see Table D.5) and thus in the forced solution (35).
Similar derivations allow to define the proper modes v of the Poincaré's variables y i , by means of a time-dependent canonical change of variables (y We have defined where N io = O T N N io . Once the proper modes u, v have been defined, it is very useful to establish a corresponding truncated Hamiltonian in the form where To define higher-degree truncations, we first consider the non-quadratic part δ H of the secular Hamiltonian (14), truncated at a given total degree in Poincaré's variables ( Then, the proper modes u, v are injected in (49) by means of Eqs. (42) and (44), and the Poincaré's variables of the outer planets replaced by their quasi-periodic decompositions given in Eq. (15). However, such substitutions do not conserve the degree of the terms in the expansion, since non-linear harmonics are present in Eqs. (15), (35) and (47), i.e. harmonics of order |m| or |n| greater than one. Such harmonics are in principle much smaller than the linear ones, and generate higher-degree terms in the substitution process. To define the terms H (2p) , p ≥ 2 in a consistent way, we thus introduce a fictitious real variable to redefine the quasiperiodic decompositions (15) as meaning that each harmonic is treated as of degree |m k | or |n k | for the purpose of the expansion 10 . The forced solutions (35) and (47) are modified accordingly. The series resulting from the substitution of variables in (49) can be thus truncated at total degree 2n with respect to Finally, addition of the quadratic terms H 2 gives, by definition, the truncated Hamiltonian (48). Truncation is effectively performed in TRIP in parallel with variable substitutions, thanks to a dedicated monomial truncated product , which allows to minimize the computational cost of the expansion. It is worthwhile to note that, at the end of the expansion process, the Hamiltonian terms only depending on time can be removed, since they are dynamically irrelevant. When numerical evaluations involving the terms of the Hamiltonian have to be performed, the fictitious variable is simply set to 1. In Table 2 we show the total number of monomials of variables {(u k ,ū k , v k ,v k ) k=1,4 , (E jφ k (t) ) k=1,7 } in the Hamiltonian H 2n according the degree of truncation. We also report the typical computation time of the series in TRIP, related to the transformation , (E jφ k (t) ) k=1,7 }. One may note the dramatic increase of the number of terms due to such change of variables. We point out that the precise number of monomials in H 2n depends on the particular quasi-periodic decomposition (15) we used in this work (see Appendix D).

Action-angle variables
Action-angle variables corresponding to the proper modes u, v are introduced by standard canonical transformations for k ∈ {1, 2, 3, 4}. The Laplace-Lagrange Hamiltonian thus reads H 2 = − g LL · X − s LL · Ψ and it is trivially integrable. Such transformations allow to expand the Hamiltonian H 2n as a Fourier series of the angle variables χ, ψ and φ(t) = ω o t,  Table 3: Frequency in arcsec yr −1 of the largest-amplitude Fourier harmonic for each proper mode of the Poincaré's variables (x k , y k ) k=1,4 over the time interval [0, 20] Myr. For the solutions LaX13b and La90 the proper modes (z k , ζ k ) k=1,4 defined in (Laskar 1990c) are employed, while for the present solutions we use the proper modes (u k , v k ) k=1,4 . where H k, 2n are complex amplitudes and we have employed a compact notation for the action-angle variables, Only a finite number of harmonics have non-zero amplitude in Eq. (52), since we deal with a truncated Hamiltonian. We also recall that, the Hamiltonian being a real function, one has 2n for all k, . In Table 2 we show the number of harmonics in Eq. (52) according the truncation degree of the Hamiltonian, i.e. the number of different wavevectors (k, ) up to a global minus sign.
It is worthwhile to mention that the explicit time dependency of the Hamiltonian H 2n can be easily absorbed in a phase-space extension, through the definition of action-angle variables for the trivial degrees of freedom of the outer orbits. Indeed, one introduces (Φ, φ) such that the new Hamiltonian reads with Φ = (Φ k ) k=1,7 . The dynamics of the additional angles is thus consistently given byφ = ω o , while that of the actions Φ is irrelevant. Such an autonomous formulation is useful in the context of canonical perturbation theory.

Comparison with reference dynamical models
It is essential, when constructing averaged models, to compare the resulting trajectories to those of the nominal Hamiltonian, or of a more comprehensive dynamical model. This allows to validate both the underlying averaging approximations and the choice of the initial conditions for the secular variables (see Sect. 2.3 and Appendix D). In this Section, we compare the orbital solution of the forced inner system (Eq. 16) with the direct numerical integration of the full equations of motion LaX13b, which has been used to predetermine the secular motion of the outer planets (Sect. 2.3 and Appendix D). We shall also employ in the comparison the orbital solution of (Laskar 1990c, La90 from now on), as it includes all the dynamical interactions taken into account in the present work and constitutes the most precise secular system to this day. The comparison between the different models begins at short secular timescales in Fig. 1, where we show the inner planet eccentricities and inclinations over 250 000 years in the future. The plots indicate that the present model correctly reproduces the secular behaviour of the LaX13b solution. This is particularly manifest in the eccentricity plots, where the forced model constitutes the average of the short-time (orbital) oscillations of the direct integration, as the La90 solution also does. In particular, the behaviour of the solutions at the origin of time shows that the initial conditions for the secular variables of this work have been correctly determined. Some small periodic differences between the two secular models appear on some of the plots, with the La90 solution generally being more precise in reproducing the long-term average of the direct integration. However, such deviations are practically irrelevant at these scales.
In Fig. 2 the comparison is extended to a longer time interval of 10 Myr, to evaluate how long the present solution remains close to the direct integration. The solution La90 practically coincides with LaX13b on these plots (see Laskar et al. 2004) and it is thus not shown. The general agreement shown over the first few million years by the curves in Fig. 2 is still very satisfactory. For longer times, the divergence of the two models becomes important. First of all, the orbital oscillations show sometime different amplitudes. Nevertheless, as we are not interested in the construction of a secular ephemeris (the solution La90 is clearly more adapted to this end), this is not the most relevant point. When considering that the inner solar system has a Lyapunov time of about 5 Myr (see Sect. 6), the crucial aspect is rather the slow dephasing of the solutions appearing from this comparison, as it translates a difference in the fundamental frequencies of the motion. To ensure that the forced model faithfully reproduces the resonant structure of the inner system, it is essential to quantify such deviations. To this end, we perform a frequency analysis of the proper modes of the Poincaré's variables (x k , y k ) k=1,4 for the different orbital solutions over the [0, 20] Myr time interval. For the solutions LaX13b and La90 we employ the proper modes (z k , ζ k ) k=1,4 defined in (Laskar 1990c), while for the present solution the proper modes (u k , v k ) k=1,4 are used (the employ of (z k , ζ k ) k=1,4 gives the same results at the numerical precision used in the comparison). Table 3 shows the dominant frequency in the Fourier spectrum for each proper mode, according the notation of (Laskar 1990c). The average absolute difference between the LaX13b solution and the present one amounts to a few hundredths of arc second per year, with a maximum of 0.06 yr −1 in the case of s 3 and s 4 . This is in agreement with the behaviour shown in Fig. 2 and the general expectation from the discussion on the model precision in Sect. 2.1. It is important to note that, when considering combinations of the fundamental frequencies (which determine the resonant structure), the differences between the two models can be even smaller. Indeed, in the case of the two leading resonances (g 1 − g 5 ) − (s 1 − s 2 ) and 2(g 3 − g 4 ) − (s 3 − s 4 ) described in (Laskar 1990c), the deviations amount to just 0.01 yr −1 and 0.004 yr −1 , respectively. Table 3 also reports the fundamental frequencies of the truncated model H 2n (Eq. 25) for different total degrees of truncation. When considering the combinations of frequencies (g 1 − g 5 ) − (s 1 − s 2 ) and 2(g 3 − g 4 ) − (s 3 − s 4 ), the models H 4 and H 6 practically show the same deviations as the full forced model H with respect to the solution LaX13b. This suggests that the truncated model should already realistically reproduce the resonant structure of the inner system at the lowest degrees. We also note that the orbital solutions of the models H 8 and H 10 practically coincide with that of the full forced model H over the first 20 Myr.

Finite-time maximum Lyapunov exponent
Once the adequacy of the present dynamical model has been established over a few million years, one has to assess if it is able to correctly reproduce the resonant structure of the inner solar system over longer times. Even if the fundamental frequencies in Table 3 suggest this is actually the case, the computation of the maximum Lyapunov exponent (MLE) remains an essential test, as its value is related to the width of the leading resonant harmonics of the Hamiltonian, which are the dynamical sources of stochasticity (Chirikov 1979). At this point, it is important to realize that the non-null probability of an unstable evolution of the inner planets (Laskar 1994(Laskar , 2008Batygin & Laughlin 2008;Laskar & Gastineau 2009) prevents the existence of the MLE as an infinite-time limit, as in its usual mathematical definition (Oseledec 1968); pertinent to the present case is the consideration of a finite-time MLE (FT-MLE). As we do not have at our disposal an efficient numerical algorithm to evaluate the variational equations of Gauss's dynamics, we consider the FT-MLE given by the two-particle algorithm proposed in (Benettin et al. 1976), which has been used in (Laskar 1989). This method computes the divergence of close trajectories of a dynamical system, by considering the motion of reference and shadow particles, initially separated by a tiny vector d 0 in the phase space. At time intervals τ, a renormalization procedure applied to the trajectory separation d(t) brings again the shadow particle at a distance ||d 0 || from the reference one. The resulting FT-MLE depends in a intricate way on the initial position of the reference particle in the phase space, so that its asymptotic evolution is chaotic (and independent of the choice of the vector d 0 (Benettin et al. 1976)). Therefore, its computation acquires full physical significance only for an ensemble of trajectories. In Fig. 3 we show the computation of the FT-MLE of the forced inner solar system over 5 Gyr in the future, for an en- On the left vertical axis, we report the FT-MLE expressed as an angular frequency in arcsec yr −1 , while on the right one the corresponding Lyapunov time, defined as FT-MLE −1 , is given in million years. The dark-grey region represents, at each renormalization time, the interquantile range of the observed probability density function (PDF) of the FT-MLE, which encloses by definition 50% of the solutions around the median (shown by the red line). In the same manner, the light-grey region corresponds to the [1 th , 99 th ]-percentile range of the PDF, enclosing the 98% of the probability. The two FT-MLEs reported in (Laskar 1989) are also shown with black lines.
As stated above, the initial conditions of each solution in the ensemble are very close to the nominal ones. They are chosen by taking the relative variation of each coordinate of the nominal phase-space vector as a normal random variable with zero mean and a standard deviation of 10 −9 . The initial conditions thus follow a multivariate Gaussian distribution centred at the nominal phase-space position of the system, with a diagonal covariance matrix. The width of such distribution is of the same order of magnitude of that in (Laskar 2008, Table 1). The initial position of each shadow particle, around the corresponding reference one, is then chosen according the same kind of Gaussian distribution, with a relative standard deviation of 10 −8 . Indeed, such a value should be close to minimize the accumulation over long timescales of numerical errors, due to the floating-point implementation of the algorithm of (Benettin et al. 1976), when working in double precision (Tancredi et al. 2001;Mei & Huang 2018). In Fig. 3 the computation of the FT-MLE for the nominal initial conditions is shown by dotted blue lines for a set of 16 different initial tangent vectors d 0 . This manifestly shows that the present FT-MLEs are asymptotically independent of d 0 , the asymptotic regime being practically reached in a few hundred million years. In the same manner, our computation has been tested to be asymptotically independent of the renormalization time τ, set to 5 Myr in our computation, and of the norm chosen for the phase-space vectors, taken here to be Euclidean as usual.
Article number, page 12 of 25 F. Mogavero & J. Laskar: Long-term dynamics of the inner solar system For times shorter than 100 Myr, Fig. 3 shows that the distribution of the FT-MLE only reflects the choice of different initial tangent vectors d 0 , around essentially the same reference trajectory, i.e. the nominal one. Therefore, this part of the plot can be discarded because of its non-asymptotic character. We also note that, for t < 100 Myr, the PDF of the FT-MLE tends to shrinks with increasing time, as expected by the asymptotic behaviour of the method employed (see the blue lines and two curves from (Laskar 1989)). For longer times, the PDF begins to broaden, because each orbital solution acquires a macroscopically different FT-MLE, which chaotically depends on its initial conditions. The time-asymptotic regions shown in Fig. 3 thus represent the probabilistic knowledge of the FT-MLE of the inner solar system, which arises from its chaotic behaviour and our determination of the current planet positions and velocities.
The dark-blue line in Fig. 4 shows the kernel density estimation (Rosenblatt 1956;Parzen 1962) of the PDF of the Lyapunov time FT-MLE −1 at 5 Gyr, along with the corresponding cumulative distribution function (CDF) on the upper horizontal axis. We employed the standard Gaussian kernel and Silverman's rule of thumb for the selection of the optimal bandwidth (Silverman 1986). We also show, through the light-blue region, the pointwise confidence interval at 98% level on the estimated PDF, obtained via nonparametric bootstrap, i.e. by resampling the original data with replacement (Efron 1979). One can appreciate that the value of about 5 Myr typically reported in literature (Laskar 1989;Sussman & Wisdom 1992) is right in the bulk of the present PDF, even though the distribution rather peaks at 3.7 Myr and the average value is 4.3 Myr. The value of 6.5 Myr found in (Rein & Tamayo 2018) can also be accounted for by our computation. The curves of (Laskar 1989) suggest that the asymptotic PDFs of the two secular models should largely overlap. Finally, we note that numerical value of 1.1 Myr reported in (Batygin et al. 2015), as the Lyapunov time of a simplified Mercury dynamics, does not agree with the present findings. Figures  3 and 4 show that, even considering the diffusion of the planet orbits over 5 billion years, such a high FT-MLE is practically never reached by the present dynamical model, even though the latter includes all the Hamiltonian harmonics taken into account in (Batygin et al. 2015).

Orbit excitation over 5 Gyr
In this Section, we perform a statistical study of the orbit excitation in the forced inner system over 5 Gyr in the future, the aim being to assess if the present model is able to reproduce the rate of the high Mercury eccentricities observed in (Laskar 2008;Laskar & Gastineau 2009). Following (Laskar 1994), we characterize the unstable evolutions of the system by defining a secular collision as the intersection of the instantaneous Keplerian orbits of a pair of planets. Numerically, such an event is easily detected by tracking the relative positions of their mutual nodes, i.e. the intersections of each orbit with the orbital plane of the other planet. When two orbits crosses, a pair of mutual nodes exchanges their positions along the line of nodes (i.e. the intersection of the two orbital planes). Such a computation is inexpensive in the context of Gauss's dynamics and can be performed at each timestep of the integration scheme. When a secular collision occurs, the present averaged model is no longer a good long-term approximation of the original N-body dynamics, in which the involved planets can experience a close encounter, possibly leading to a physical collision or to the escape of one of two planets from the solar system (Laskar 1994). We thus denote by τ c the time of the first secular collision for a given orbital so-  LG09 represents the 2 501 orbital solutions in (Laskar & Gastineau 2009), while La08 stands for the 478 integrations in (Laskar 2008). The statistics of this work derives from an ensemble of 10 560 orbital solutions.
lution. The statistics of the first-collision time τ c shall provide an estimate of the rate of catastrophic events in the real inner system. Moreover, by stopping the numerical integration at the first secular collision, we will not be involved in extrapolating the present secular dynamics out of its validity range, and the use of a fixed-timestep integration scheme (see Sect. 3) turns out to be sufficient end even suitable. We note that a secular collision between a planet and the Sun can be simply defined to occur whenever its pericenter distance a(1 − e) is smaller than the Sun radius.
We perform an ensemble of 10 560 numerical integrations of the forced inner system over 5 Gyr. The initial conditions are chosen according the same multivariate Gaussian distribution centred at the nominal values as in Sect. 6. The relative standard deviation of each phase-space coordinate is 10 −9 , of the same order of magnitude of that in (Laskar 2008, Table 1). Table 4 shows the observed probability of having, over the entire time span, a maximum Mercury eccentricity greater than a given value. The statistics of the present dynamical model is compared to those resulting from (Laskar 2008 , Table 3) and (Laskar & Gastineau 2009, see Table C.1 in Appendix C of the present paper). The statistical bounds represent the Wilson (1927) score interval corresponding to a 98% confidence level. At the present statistical precision 11 , the probabilities of the high Mercury eccentricities in the forced inner system are fully compatible with those arising in the secular model of (Laskar 1990c). They also represent a very good estimate of the probabilities resulting from N-body dynamics, even if somewhat lower. This is in agreement with the general expectations from an averaged model, which has less degrees of freedom (Laskar 1994(Laskar , 2008. When comparing the midpoints of the statistical intervals, Table 4 could also suggest that the true probabilities of the present model are slightly lower than those in (Laskar 2008). If true, such a fact would agree with the expectation of higher regularity of the forced model, as the degrees of freedom dominated by the outer planets (their fundamental frequencies in particular) are frozen into a quasi-periodic time dependency. Moreover, the Hamiltonian terms at second order in planet masses could also contribute to somewhat larger chaotic excursions of the inner orbits. 11 As there is no reason for the true probabilities of the three different dynamical models to be exactly the same, for a very large number of orbital solutions their statistical intervals are not likely to overlap.
Article number, page 13 of 25 A&A proofs: manuscript no. iss

Secular collisions
Among the 10 560 orbital solutions of our ensemble, we found 42 secular collisions, corresponding to a 98% confidence interval (0.28%, 0.57%) for the probability P(τ c ≤ 5 Gyr). Out of the 2501 integrated solutions of (Laskar & Gastineau 2009), there result 6 collisions between Mercury and the Sun, 9 collisions between Mercury and Venus, and 1 collision between Mars and the Sun (see Appendix C). These 16 events corresponds to a 98% confidence interval (0.36%, 1.13%) of probability. Our definition of secular collision thus provides a very good estimate of the probability of a physical collision in the non-averaged inner solar system over the next 5 Gyr, probably somewhat lower in agreement with the previous discussion. All the secular collisions in our ensemble of orbital solutions involve the Mercury-Venus couple, with a maximum Mercury eccentricity ranging from 0.797 to 0.966. In particular, in 41 out of 42 collisional solutions it is larger than 0.8, while in 15 solutions it exceeds 0.9. This implies that the statistics corresponding to the cases e max = 0.8 and 0.9 in Table 4 is affected by the choice of stopping the numerical integration when the first secular collisions is detected, independently of the fact that this could not physically correspond to any unstable further evolution. We emphasize that, as a consequence, the statistics of the secular collisions other than the Mercury-Venus one cannot be correctly reproduced by the present numerical computation, since such events, as observed in (Laskar & Gastineau 2009), always require a very large Mercury eccentricity at some previous time.
7.2. Secular resonance g 1 − g 5 The destabilizing role of the secular resonance g 1 − g 5 , able to drive Mercury orbit to very high eccentricities, has been first established numerically (Laskar 2008;Batygin & Laughlin 2008) and then confirmed by an analytical study (Boué et al. 2012). We conclude this Section by showing that the g 1 − g 5 resonance is indeed reached by the collisional orbital solutions of our ensemble. While g 5 ≈ 4.257 yr −1 is constant in the present dynamical model, the frequency g 1 , which dominates the spectrum of the Poincaré's variable x 1 and thus the precession of Mercury perihelion, is numerically computed in the following way. Given an orbital solution, we sample the corresponding Poincaré's variables x with a timestep ∆t = 1 000 yr, i.e. we consider the time series x(t n ) with t n = n∆t and n ∈ N 0 . The proper mode u 1 (t n ) is thus computed through Eq. (42) and the corresponding angle variable χ 1 (t n ), following the definition in Eq. (51), retrieved as an unwrapped phase, i.e. as a continuous function of time. The instantaneous (angular) frequency ω 1 of the proper mode u 1 is given by (Cohen 1995) We could analytically compute the derivative of Eq. (42) to obtain the time series ω 1 (t n ). In practice, we employ numerical differentiation via the Lagrange five-point formula (Olver et al. 2020, Chapter 3), as very high numerical precision is unnecessary here. As it is, the frequency ω 1 (t n ) presents short-time oscillations which hide its long-term evolution. To suppress such high-frequency fluctuations, we thus define the time series g 1 (t n ) as the output of the low-pass Kolmogorov-Zurbenko (KZ) filter (Yang & Zurbenko 2010, and references therein) applied to ω 1 (t n ), The KZ filter is defined as an iteration of the common moving average and has been comprehensively characterized in Appendix B. For the present application, we use three iterations of the moving average and a cutoff frequency of 1 Myr −1 (see Appendix B). This particular choice is motivated by the typical duration of the libration episodes of the g 1 − g 5 resonance, which only last a few million years (Laskar 2008;Batygin & Laughlin 2008).
In the top panel of Fig. 5 we present the time evolution of the fundamental frequency g 1 along the 42 collisional integrations. We isolate in the bottom panel six solutions to show the typical evolution of a single curve more clearly. As previously discussed, the curves stop at the first secular collision. The horizontal line represents the constant fundamental frequency g 5 and thus the location of the g 1 − g 5 resonance. Within the first 2 Gyr, the chaotic diffusion of g 1 along the solutions of Fig.  5 is limited to no more than 0.25 yr −1 with respect to its initial value. Nevertheless, when lower values of about 5 yr −1 are eventually reached, the dynamics of g 1 starts to be characterized by large random variations. This seems to reproduce the large chaotic zone related to the resonance g 1 − g 5 , with a half-width of more than 1 yr −1 , which has been reported in (Laskar 2008, Fig . 6(a2)) 12 . Along all the curves, the fluctuations eventually lead to a crossing of the resonance (corresponding locally to a libration of the related argument), shortly before the first secular collision.

A numerical experiment over 100 Gyr
In this last Section, we present the results of a numerical experiment consisting in integrating an ensemble of orbital solutions of the forced inner system over 100 Gyr in the future. Through such simulations, we aim to characterize the kind of stochastic process which arises from the Hamiltonian (16) and drives the destabilization of the inner system, in a regime where the highly excited orbits do not represent rare events anymore, as they do in Sect. 7. Obviously, this kind of simulation does not correspond in anyway to the real future evolution of the orbits of the inner planets over much more than 5 Gyr. Indeed, it is well known that over the next 7 Gyr, the Sun will experience a significant mass loss, down to 0.54 M , as it leaves the main sequence (Sackmann et al. 1993). As a consequence, the semi-major axes of the planet will adiabatically inflate by a factor of 1.85. Even more critically for the fate of the inner system, it is likely that the innermost planets Mercury, Venus and the Earth, will be engulfed by the Sun as its radius expands along the red giant branch, marking the end of their existence (Rybicki & Denis 2001;Schröder & Smith 2008). Therefore, this Section intends to investigate the very long-term stochastic behaviour of the dynamical system (16), without inferring any conclusion on the real evolution of the inner system over such timescales.
We numerically integrate an ensemble of 1 042 orbital solutions of the forced inner system over 100 Gyr. The initial conditions are chosen according the same Gaussian distribution as in Sect. 7, but with different realizations, to obtain different evolutions over the first 5 Gyr. As in the previous Section, the numerical integration is stopped at the first secular collision. In Fig. 6 we show the probability density function of the time τ c of the first secular collision. In the context of stochastic processes, this represents the hitting time (or first hit time) corresponding to the subspace of the phase space realizing a secular collisions between planets. The blue curve in the lower panel of Fig. 6 represents the PDF estimated via the kernel density estimation (KDE) method (Rosenblatt 1956;Parzen 1962). We employed a standard Gaussian kernel and Silverman's rule of thumb for the selection of the optimal bandwidth (Silverman 1986). The PDF is normalized in such a way that its integral over 100 Gyr is equal to the overall percentage of secular collisions obtained in the experiment, i.e. 890/1042 ≈ 85.4%. The blue region in the lower panel of Fig. 6 represents the pointwise confidence interval for the PDF at 98% level, obtained via nonparametric bootstrap, i.e. by resampling the original data with replacement (Efron 1979). This region thus corresponds pointwisely to the [1 th ,99 th ]-percentile range of the bootstrapped estimated PDF. We emphasize that, for simplicity, we do not correct for the bias arising in the kernel density estimation, so that it constitutes, strictly speaking, a confidence region for the expectation of the KDE and not for the true subjacent PDF (Chen 2017). When referred to the true PDF, such a confidence region typically represents an undercoverage at the points where the function bends, i.e. where its second derivative is large. Even if we could consider simple corrections for the KDE bias (Chen 2017), this would not be completely satisfactory given the very small number of events which we deal with over the first few Gyr. Anyway, such corrections on the confidence intervals turn out to be sufficiently small to not affect at all the kind of conclusions we infer in this Section. To give a further idea of the statistical variance, in Fig. 6 we also show in red colour the PDF of τ c over the first 5 Gyr estimated from the ensemble of 10 560 orbital solutions presented in Sect. 7. We finally report in the upper panel of Fig. 6 the estimated cumulative distribution function (CDF) as a function of time, which gives the probability of a collisional evolution within a given time, along with its pointwise confidence interval at 98% obtained through bootstrap.

Marginal stability of the inner solar system
The PDF of the time τ c of the first secular collision peaks at 17.6 Gyr, with a 98% confidence interval (15.4, 20.1) Gyr, while its median is 40.8 Gyr, with a 98% confidence interval (38.2, 43.5) Gyr. The percentage of collisional evolutions over 100 Gyr has a 98% confidence interval (82.8%, 87.9%). These findings clarify the state of matastability which characterizes the inner solar system, according to the definition of (Laskar 1996). Over a timescale comparable to its age, i.e. over the next 5 Gyr, the inner system is statistically very stable, with a probability of only about 0.5% of an unstable evolution (see Sect. 7.1). Nevertheless, here we find that such probability increases very fast with time, so that the system rapidly becomes statistically unstable over longer timescales, meaning that unstable evolutions do not represent rare events anymore. In particular, the probability of an instability is about 20% at 20 Gyr already, and 50% at 40 Gyr.
In the conjecture of (Laskar 1996), this picture should characterize any secularly evolving planetary system: at each moment of its formation history, a planetary system should be in a state of marginal stability, that is practically stable over a timescale comparable to its age, while strong instabilities arise very fast over longer times, resulting in more stable configurations of the surviving planets.
Article number, page 15 of 25 A&A proofs: manuscript no. iss 8.2. Destabilization of the inner solar system as a stochastic process Recently, in the framework of the simplified model introduced in (Batygin et al. 2015), Woillez & Bouchet (2020) applied the theory of the white noise limit for slow-fast dynamical systems (e.g., Gardiner 1985), to describe the long-term dynamics of Mercury through an effective stochastic process. They show that the instability timescale of Mercury orbit is reproduced by the hitting time of a one-dimensional Wiener process with a reflecting barrier. The corresponding PDF peaks at 1.3 Gyr (Woillez & Bouchet 2020, Figure 4), in agreement with the analysis in (Batygin et al. 2015). Such a timescale is an order of magnitude smaller than that in Fig. 6, and the simplified dynamics of Mercury predicts a large probability of instability over the solar system age, in contrast with the findings of realistic models.
In this Section, we reconsider the stochastic process discussed in (Woillez & Bouchet 2020) to characterize the overall structure of the PDF in Fig. 6. To fix ideas, it could describe in a crude way the long-term wanderings of the fundamental frequency g 1 towards the destabilizing secular resonance g 1 − g 5 (see Sect. 7.2). Indeed, numerically, the PDF of g 1 turns out to be asymmetric, with a longer tail diffusing towards low values, while its upper end is practically fixed at about 6 yr −1 (Hoang et al. 2021, see also Fig. 5). Therefore, let suppose that g 1 performs a random walk described by a Wiener process starting at g 1,0 = g 1 (t = 0), and consider an upper reflecting barrier at g 1,max > g 1,0 . We assume that a secular collisions rapidly occurs whenever the secular resonance g 1 − g 5 is reached, as shown in Fig. 5; this is taken into account by considering an absorbing barrier at g 1 = g 5 . The PDF ρ(τ) of the hitting time τ = inf t {g 1 (t) ≤ g 5 } is thus given by (e.g., Schwarz 1992; Woillez & Bouchet 2020) where α = (g 1,0 − g 5 )/(g 1,max − g 5 ) and T 1 = (g 1,max − g 5 ) 2 /4D, with D being the diffusion coefficient of the Wiener process. After a least-squares search, we plot in Fig. 6 as a dashed black line the curve (57) with parameters T 1 = 27.6 Gyr and α = 0.9. Even if the real random walk performed by the frequency g 1 , and the destabilization mechanism in general, are likely to be much more complex, a Wiener process with a reflecting barrier is able to reproduce the overall behaviour of the observed PDF, at the level of the present statistical precision and apart from a certain excess of events in the tail of the distribution 13 . On the one hand, this shows that at short times the PDF behaves as where T 0 = α 2 T 1 . In this regime, the PDF coincides with that of the hitting time of a standard Wiener process (i.e., without reflecting barrier). The very fast decaying of the density probability for τ → 0 in Eq. (58) characterizes the rare destabilizations of Mercury orbit over the first few billions of years. In particular, the instanton phenomenology described in (Woillez & Bouchet 2020) could indeed apply to the real solar system (such a hint still needs to be tested through integration of realistic models, as the present one.). On the other hand, Eq. (57) shows that at large 13 We stress that such a result is independent of the hypothesis of which precise dynamical quantity (if any) actually undergoes the stochastic process considered here.
times the decay of the observed PDF is at least exponentially fast, with a characteristic time equal to 4T 1 . This behaviour results from the reflecting barrier taken into account in the stochastic process; it would not be reproduced by a standard Wiener process, as the corresponding PDF has a polynomial heavy tail proportional to τ −3/2 , as shown by Eq. (58). In particular, the average of the hitting time as given by Eq. (57) is finite, differently from Eq. (58). We note that (Woillez & Bouchet 2020) reported T 0 = 1.56 Gyr, a value which is more than ten times smaller than that in Fig. 6, i.e. T 0 = 22.4 Gyr, in agreement with the previous discussion about the validity of the simplified dynamics of Mercury in (Batygin et al. 2015). Here we emphasize that the incompatibility with realistic models seems to be related to the limitations of some of the simplifying assumptions 14 . Indeed, as discussed in Sect. 7, averaged models generally tends to produce slower instabilities with respect to the full dynamics, because the number of degrees of freedom is smaller. Nevertheless, the model of (Batygin et al. 2015) produces instabilities which are ten times faster than in the present dynamics, even though the latter includes all the dynamical interactions considered in the simplified one. Such a lack of continuity in the destabilization time with respect to the complexity of the dynamical model seems rather severe.
We conclude this Section by noting that a lower bound for the typical instability time of Mercury orbit can be estimated from the maximum Lyapunov exponent (Fig. 3), if one assumes that no dynamical constraints prevent g 1 to diffuse at a rate determined by the leading secular resonances. In particular, the resonance (g 1 − g 5 ) − (s 1 − s 2 ), with a libration frequency of 0.12 yr −1 , is among the leading ones affecting the frequency g 1 (Laskar 1990c). We may estimate an upper bound on the diffusion coefficient D as with 2π MLE = 0.1 yr −1 (this is also the lower limit for the MLE in Fig. 3). A lower bound on the typical destabilization time of Mercury orbit is thus given by with g 1,0 = 5.577 yr −1 (Table 3) and g 5 = 4.257 yr −1 . Eq. (60) gives essentially the same diffusion timescale as in (Batygin et al. 2015;Woillez & Bouchet 2020). This is already two order of magnitude greater than the Lyapunov time of the system, but still insufficient to explain the findings of Fig. 6. This means that some dynamical constraints effectively makes the diffusion towards the g 1 − g 5 resonance slower than in the estimates (59), (60). Certain secondary resonant harmonics, in particular, must play a determinant role in the destabilization of the inner solar system.

Conclusions
This work introduces a new secular model for the dynamics of the inner planets of the solar system. It exploits the practical constancy of the fundamental precession frequencies of the outer planet orbits over timescales of billion years, the smallness of the inner planet masses and the absence of strong mean-motion resonances in the inner system. The range of validity of the resulting dynamics extends to very high eccentricities and inclinations, up to a secular planetary collision, i.e. the geometric intersection of the instantaneous Keplerian orbits of two planets. The new model can be studied analytically by truncating the Hamiltonian at a given total degree in planet eccentricities and inclinations, with the aid of a computer algebra system like TRIP. It can also be integrated numerically in a very efficient way through Gauss's method. The orbital solution matches the predictions of a comprehensive model of the solar system at a very satisfactory level over timescales shorter than or comparable to the Lyapunov time. The new model also correctly reproduces the maximal Lyapunov exponent of the inner system and the statistics of the highly eccentric Mercury orbits over the next 5 billion years. Moreover, the destabilizing role of the secular resonance g 1 − g 5 clearly stands out. We performed a numerical experiment consisting of a thousand orbital solutions of the inner solar system over one hundred billion years, to explore a regime in which unstable orbits are statistically common. We first pointed out the fast growth of the rate of orbit instabilities in the framework of planetary system formation through successive metastable states. We then showed that the PDF of the time of the first secular collision is reasonably well reproduced by a Wiener process with a reflecting barrier, which could be performed, for example, by the fundamental precession frequency g 1 . Given the robustness of the statistical predictions of the present dynamical model over 5 Gyr, the main properties of this PDF, that is the characteristic time of its peak and the behaviours at short and long times, are likely to represent what would arise from the full dynamics of the planets. We finally argued that a dynamical mechanics is needed to explain the rarity of the large excursions of g 1 up to the g 1 − g 5 secular resonance within the next 5 Gyr. We emphasize that the new dynamical model can be straightforwardly implemented once truncated at degree 4 in planet eccentricities and inclinations, by using the corresponding expression of the two-body secular Hamiltonian reported in (Laskar & Robutel 1995, Appendix) and the quasi-periodic secular solution for the outer planets given in Appendix D.
where =r/µ denotes the velocity vector tangent to the instantaneous Keplerian ellipse of a given planet. In the reference frame of the Keplerian orbit, with origin at the Sun, x-axis directed towards the pericenter and z-axis normal to the orbital plane, the components of the vector read where E is the eccentric anomaly, n = G(m 0 + m)/a 3 the mean motion and r = ||r|| the heliocentric orbital distance. After transformation to the fixed reference frame, the velocity components are given by where the rotation matrices R(i, Ω) = R 3 (Ω)R 1 (i)R 3 (−Ω) and R 3 ( ) are given in (Laskar 1991;Laskar & Robutel 1995). We introduce the eccentric longitude F = E+ , which obeys the modified Kepler's equation where z = eE j and γ = 1 − XX/4. By using the fact that 1 − r/a = γ Re X E − jF , and after some algebra, one finds By employing the complex formalism of (Laskar & Robutel 1995), one can thus write The application of matrix R(i, Ω) in Eq. (A.3) then gives with δ = 1 − XX/2 and η = 1 − XX/2 − YȲ. The scalar product in Eq. (A.1) thus reads The coefficients T , are readily derived once the quantities a/r and (a/r)E jF , appearing in Eq. (A.6), are expanded in Fourier series of the mean longitude λ. One has By using Eq. (A.4) and the fact that dλ/dF = r/a, a classical calculation gives c = J ( e) E − j , where the J are the Bessel functions of the first kind. The following alternative expression of the Fourier coefficients, can be straightforwardly expanded in series of X andX, and is well-defined at zero eccentricity. In a similar manner, one obtains

Appendix B: Kolmogorov-Zurbenko filter
Low-pass filters are employed in time-series analysis to extract the long-term (low-frequency) components of a given signal. A dedicated low-pass filter has been constructed in (Carpino et al. 1987), to recover the secular changes of the giant planet orbits from the numerical output of a N-body integration. Here we propose the use of the Kolmogorov-Zurbenko filter (Yang & Zurbenko 2010) as an out-of-the-box, computationally advantageous and still effective choice. Let consider the real-valued time series ξ n = ξ(n∆t), n ∈ Z, (B.1) resulting from the discrete sampling of a continuous signal ξ(t), with constant sampling rate ω s = 2π/∆t. In the frequency domain, the time series is characterized by its discrete-time Fourier transform which is a periodic function with period ω s and can be thus restricted to the interval |ω| ≤ ω s /2. The action of a finite-impulseresponse filter F on the time series is defined as the discrete convolution of the signal ξ n with a given finite sequence d n , M being a non-negative integer. The sequence d n is the impulse response function of the filter, and L M = 2M + 1 is the filter length. A liner filter is characterized by the discrete-time Fourier transformd(ω) of its impulse response, which is called the frequency response of the filter and gives the spectrum of the output signal as while the response should be unitary at lower frequencies, i.e. d(ω) = 1 for ω ≤ ω c . Moreover, the filter must not change the phase of the Fourier components of the signal, i.e. Im(d(ω)) = 0. This is achieved by requiring d n = d −n ∈ R (Carpino et al. 1987). Nevertheless, a real filter is rather characterized by: a passband, |d(ω) − 1| ≤ ρ for |ω| ≤ ω p , ρ being a given maximum loss; a stopband, |d(ω)| ≤ α for |ω| ≥ ω c , α being a given maximum gain; a transition band between ω p and ω c , where the response function smoothly decreases from ≈1 to ≈0.
One of the simplest low-pass filters is the moving average (MA), which simply constitutes the local unweighted average of the signal over a time window L M ∆t. Computationally, the moving average is very advantageous; however, it is well-known that it provides a poor attenuation of the Fourier components of the signal in the stopband (see Fig. B.2). Kolmogorov proposed to bypass such a problem by applying the moving average iteratively (Yang & Zurbenko 2010, and reference therein). The Kolmogorov-Zurbenko filter is thus defined as The output signal of such a filter is given by   (Laskar & Gastineau 2009). In the 3 bottom rows of the last column (5 Gyr), the second number is the number reached after full completion of all the integrations, in June 2010.
In the last column of Table C.1, the number of solutions for which the eccentricity of Mercury extends beyond 0.9 increased from 19 to 21. Among these 21 solutions, 6 present a collision of Mercury with Venus (i.e. their center of mass distance is smaller than the sum of their radii), 9 a collision of Mercury with the Sun, 5 reach 5 Gyr without collision, but with very close encounters (for example, one solution presents an encounter of Mercury with Venus with less than 1800 km between the two surfaces), and one solution has a very close encounter of Mars with the Earth (surface distance of less than 794 km) followed by a collision of Mars with the Sun. In this work, we use a quasi-periodic solution for the secular evolution of the orbits of the outer planets (Jupiter to Uranus), denoted hereafter as QPSO. In order to achieve a realistic modelling, this solution is deduced from a numerically integrated full solution of the solar system. The solution of reference (LaX13b) uses an identical model as in ), but which has been initially fitted to the improved INPOP13b high precision planetary ephemeris Fienga et al. 2014), extended over 1 Myr. The model is then integrated over 30 Myr to derive a quasi-periodic approximation through frequency analysis (Laskar 1988(Laskar , 2005. The variables considered here are the Poincaré complex canonical variables in Eq. (2), scaled to suppress the semi-major axis (or Λ) dependence, that is It should be noted that X = z + O(e 3 ) and Y = ζ + O(e 2 sin(i/2)) where z = eE j and ζ = sin(i/2)E jΩ are the classical, non canonical, complex elliptic elements. The derivation of the QPSO model needs some care. Indeed, we want to obtain a model that contains only the secular frequencies pertaining to the outer planet system. The secular terms related to the inner planets are thus discarded. Because of the presence of the 5:2 close mean-motion resonance among Jupiter and Saturn, and the proximity of the 2:1 mean motion resonance in the Uranus-Neptune system, several short period terms (with arguments involving the planetary mean motions) have also to be discarded. The solution comprises a small number of terms and is fully given in Tables D.5 and D.6. It should be noted that only the terms for which the angular argument is recognized in an unambiguous way as a combination of the fundamental frequencies (Table D.3) are selected, as we want to derive an analytic model for the outer planets variables. These are expressed in the form where Z stands for X or Y, alternatively, ω o = (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ) is the vector of the fundamental secular frequencies of the outer planet system, given in Table  D.3, and k n ∈ Z 7 is a 7-uple of integers.
In order to test this model, we have compared it to the full solution LaX13b (Figs. D.1 and D.2). Only the real parts of X and Y are represented, as the imaginary parts lead to very similar plots. The full solution from LaX13b is plotted in purple, while the residuals, once the analytical model QPSO removed, are plotted in green. For the inclination variables (Y k ), these residuals are very small and appear as a straight line in the plots (Fig. D.2). This is not the case for the eccentricity variables (X k ) where a significant band of residuals appears in green (Fig. D.1). This is mostly due to the short-period terms that are present in the full LaX13b solution. When the largest of these terms are removed, the residuals become much smaller, as shown by the black curves (see also Table D.4).

Appendix D.1: Secular initial conditions
To complete the QPSO solution, or to integrate the secular equations by Gauss's averaging method (Sec. 3), one needs the val-ues of the secular semi-major axesâ k , or more precisely those of the secular canonical variablesΛ k ∝ √â k , which are constant in the secular system. In this work, the secular semi-major axes are defined as the square of the secular averages √ a k along the LaX13b solution (Table D.1), i.e.â k = √ a k 2 . The mean mean motions N k (secular averages of the mean motions) and the values of the mean longitudes at the origin λ 0k are derived as well from the LaX13b solution and provided in Table D.1. The remaining initial conditions of the secular system, corresponding to the complex variables (X k , Y k ), are obtained by a least-square fit of a polynomial of degree 5 in time to the full LaX13b solution over the first few thousand years, after removal of the shortperiod component by Fourier filtering (Table D.2).
arcsec yr −1 g 5 + 4.2574706495769208e+0 g 6 + 2.8245402509991674e+1 g 7 + 3.0879599203482901e+0 g 8 + 6.7303498995492750e−1 s 6 − 2.6347830405033751e+1 s 7 − 2.9925307659382381e+0 s 8 − 6.9173578620513787e−1 Table D.3: Values of the fundamental secular frequencies of the outer solar system derived by frequency analysis of the LaX13b solution over 30 Myr. These are the values actually used in the QPSO model. They may slightly differ from the corresponding values given in (Laskar et al. 2004. 0 0 − 1 0 1 0 11 − 5.1451015466983096e−6 − 6.1085028343284465e−7 1 0 1 0 0 − 1 0 Table D.6: Quasi-periodic decomposition (Eq. D.2) of the secular solution for the outer planets for the eccentricity variable Y. The first column is the index n of the terms ranked by decreasing amplitude. Columns 2 and 3 are the real and imaginary parts of the complex amplitude Y n . The last 7 columns are the integer coefficients k n of the secular fundamental frequencies (g 5 , g 6 , g 7 , g 8 , s 6 , s 7 , s 8 ). Article