Issue 
A&A
Volume 655, November 2021



Article Number  A1  
Number of page(s)  25  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202141007  
Published online  28 October 2021 
Longterm dynamics of the inner planets in the Solar System
Astronomie et Systèmes Dynamiques, Institut de Mécanique Céleste et de Calcul des Éphémérides, CNRS UMR 8028, Observatoire de Paris, Université PSL, Sorbonne Université,
77 avenue DenfertRochereau,
75014
Paris,
France
email: federico.mogavero@obspm.fr
Received:
6
April
2021
Accepted:
30
May
2021
Although the discovery of the chaotic motion of the inner planets in the Solar System (Mercury to Mars) was made more than 30 years ago, the secular chaos of their orbits still demands more analytical analyses. In addition to the highdimensional 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 quasiperiodic 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 meanmotion resonances permits us to restrict ourselves to firstorder secular averaging. The resulting dynamics can be integrated numerically in a very efficient 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 five billion years. The destabilizing role of the g_{1} − g_{5} 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 clarifies its current metastable state.
Key words: chaos / celestial mechanics / planets and satellites: dynamical evolution and stability / methods: analytical / methods: numerical
© F. Mogavero and J. Laskar 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The inner Solar System, with the secular chaos of its planetary orbits (Laskar 1989, 1990c; Sussman & Wisdom 1992), holds a special status among the dynamical systems of celestial mechanics. Even though 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 fundamental precession frequencies of Earth and Mars, as proposed in Laskar (1990c, 1992) and supported by Laskar et al. (2004). This fact is probably related 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 dozen million years (Myr), with the exception of the Venusdominated eccentricity mode g_{2}, which has somewhat smaller variations (Laskar 1990c; Laskar et al. 2004). This behaviour reveals the essential highdimensional structure of the inner Solar System, which limits the possibility of faithfully modelling its dynamics with a few effective degrees of freedom (e.g. Lithwick & Wu 2011; Batygin et al. 2015).
Analytical insight into the motion of the inner planets requires an appropriate dynamical modelling of the longterm 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, 2008; Laskar et al. 2004; Laskar & Gastineau 2009) to ensure that it reproduces the dynamical features of the inner system with sufficient precision. On the other hand, the corresponding Hamiltonian should be set 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 highdimensional dynamics. Moreover, the possibility of numerically integrating the equations of motion in an efficient way is fundamental for studying the chaotic evolution of the orbital solutions in a statistical way. Unfortunately, the construction of such a model is a delicate task. In principle, we might just consider the full Nbody 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, this Hamiltonian is unnecessarily complicated because it includes shorttime harmonics, with periods shorter than 5000 yr (e.g. Carpino et al. 1987), which are known to generate small quasiperiodic oscillations in the inner orbits, without being implied in chaos generation. The inner planets are not, indeed, involved in any relevant meanmotion resonance. At the same time, the longterm 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 considering a secular Hamiltonian to directly describe the slow movement of the planet perihelia and nodes after proper averaging over the shorttime orbital motion (Laskar 1984, 1985). Secular dynamics includes the essential planet interactions responsible for chaos in the inner Solar System and allows performing the fastest longterm numerical integrations (Laskar 1988, 1989, 1994, 2008). This enabled the discovery of chaos in the inner system 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 meanmotion resonance between Jupiter and Saturn, the socalled great inequality (Laplace 1785; Laskar 1996). A simple averaging of the Nbody Hamiltonian over the planet mean longitudes, resulting in a firstorder secular dynamics in planet masses, would reproduce the fundamental frequencies g_{5} and g_{6} (which dominate the perihelion precession of Jupiter and Saturn, respectively) very poorly^{1} (Laskar 1988). The construction of higherorder 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 longterm regularity ofthe outer planet trajectories (Laskar 1990c; Laskar et al. 2004; Hoang et al. 2021), the low masses of the inner planet, and the absence of relevant meanmotion resonances in the inner system (Sect. 2). We show that our model can be numerically integrated with the socalled Gauss (1818) method (Sect. 3), while its Hamiltonian is suitable for a systematic expansion in planet eccentricities and inclinations when the algorithm of Laskar & Robutel (1995) and a computer algebra system such as TRIP (Gastineau & Laskar 2011, 2020) are employed (Sect. 4). We compare its orbital solution to a reference Nbody integration over short (secular) timescales (Sect. 5) and compute the corresponding maximum 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 that drives the destabilization of the inner Solar System. We also discuss the implications on the conjecture of a marginal stability of a secularly evolving planetary system formulated in Laskar (1996) (Sect. 8).
2 Dynamical model
We modeled the dynamics of the largest bodies in the Solar System by considering the Sun and the N = 8 planets as point masses m_{0}, , indexed in order of increasing semimajor axis. The barycentric coordinates of the bodies and the corresponding momenta are denoted by , . By employing the canonical heliocentric variables of Poincaré (1896), the Hamiltonian of the Newtonian gravitational interactions among the bodies reads (1)
where r_{k} = R_{k} −R_{0} are the planet heliocentric coordinates and 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 unionof disjoint Kepler problems, so that it is useful to introduce a set of canonical variables that trivially integrates the unperturbedproblems (e.g. Morbidelli 2002). From now, we adopt the momentumcoordinate ordering of conjugate pairs. Appropriate variables are(Λ_{k}, λ_{k}; y_{k}, −jȳ_{k}_{)k=1,N}, defined as (2)
where a_{k} are the planet semimajor 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, 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é rectangular coordinates in complex form. We refer to them throughout as Poincaré’s complex variables, or simply Poincaré’s variables. With this choice of canonical variables, the integrable part of the Hamiltonian (1) reads (3)
so that the Poincaré complex variables are constants of motion for the Kepler problem.
In the regime of low orbital eccentricities and inclinations, which characterizes the Solar System, Eqs. (2) give and , and the Poincaré variables are also small. The principal part of the twobody perturbation in Eq. (1) can therefore be expanded as a Fourier series in the planet mean longitudes, with polynomial coefficients depending on the Poincaré variables (e.g. Laskar 1990b, 1991; Laskar & Robutel 1995), (4)
where is a tuple of nonnegative integers^{2}. Following Laskar & Robutel (1995), we defined the dimensionless Poincaré variables , and the semimajor axis ratio α = a∕a′, with a < a′. The analytical expression of the coefficients , only depending on the semimajor 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, (5)
The computation of the Fourier coefficients is outlined in Appendix A. Equations (4) and (5) allow explicitly computing the Fourier expansion of the Hamiltonian perturbing function, (6)
where λ stands for the vector of the planet mean longitudes, λ = (λ_{1}, …, λ_{8}), and the coefficients depend on all the remaining canonical variables.
2.1 Secular dynamics
The longterm dynamics of the Solar System planets, in particular that of the inner planets, essentially consists of the slow precession of their perihelia and nodes, driven by secular, that is, orbitaveraged, gravitational interactions (Laskar 1990c; Laskar et al. 2004). A secular Hamiltonian describing this longterm 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 nonnull wave vectors ℓ 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 time1 flow of a generating Hamiltonian S satisfying the homologic equation (7)
where the braces represent the Poisson bracket. The anglebracket operator stands for averaging over the mean longitudes, (8)
with the integration defined over the hypertorus at fixed values of all the remaining canonical variables. This means that ⟨H_{1} ⟩ is the Fourier coefficient corresponding to the null harmonic in the expansion (6). The homologic Eq. (7) gives the generating Hamiltonian S as a formal Fourier series, (9)
where n = ∂H_{0}∕∂Λ is the vector of the planet mean motions. The secular Hamiltonian is formally given by the Lie transform generated by the function S and applied to the Hamiltonian H, (10)
where L_{S}⋅ = {S, ⋅} is the Lie derivative associated with the generating Hamiltonian, is defined as the identity operator, and for n ≥ 1. The Hamiltonian in Eq. (10) is expressed in the new canonical variables , which we call the secular variables. They are related to the original variables by the Lie transforms (11)
The original variables are therefore the superposition of the secular variables and shorttime oscillations generated by the n ≥ 1 terms of the Lie transforms. Because we are interested in the longterm dynamics of the planets, we only focus on the secular variables. To keep a simpler notation, we therefore omit the hat on the secular variables from now on.
Differently from the outer planets, the inner planets are not currently involved in strong meanmotion resonances. More precisely, a maximum contribution of only 0.07′′ yr^{−1} to the fundamental precession frequencies of the inner orbits in the LaplaceLagrange 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 of Jupiter and Saturn. Building on this fact, we chose to truncate the series (10) at first order in the planet masses, that is, we neglected quadratic and higherorder terms with respect to the Fourier coefficients . We expect the main contribution to the precession frequencies of the inner orbits to come from the linear terms, given the low masses of the inner planets. In absence of strong meanmotion resonances, in Eq. (9) the denominators ℓ ⋅ n involving at least one inner planet are sufficiently far from zero. Under these assumptions, the higherorder terms only produce small corrections to the dynamics generated by the leading terms. Using the homologic Eq. (7), we thus obtain (12)
The resulting secular Hamiltonian is simply the average of the Nbody 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 these 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 constants of motion in the secular dynamics, and so are the semimajor axes of the planets.The resulting Hamiltonian system consists of two degrees of freedom for each planet^{3}, corresponding to the two Poincaré 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). It is well known that general relativity contributes with 0.430′′ yr^{−1} to the secular precession of the Mercury perihelion. This correction is critical for the statistics of the longterm destabilization of the inner orbits (Laskar 2008; Laskar & Gastineau 2009) because 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 therefore included in the Hamiltonian (12) the leading secular contribution of general relativity, which reads (13)
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 semimajor 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 EarthMoon 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 is neglected in our model. At degree 2 in planet eccentricities and inclinations, Laskar (1985, Table 8) reported contributions from the second order of 0.073′′ yr^{−1} and 0.062′′ yr^{−1} to the Earthdominated and Marsdominated eccentricity modes g_{3} and g_{4}, respectively,the other frequency corrections being at least ten times smaller. Because 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 therefore chose to exclude the EarthMoon interaction in our model, preferring to work with a Hamiltonian whose dependence on eccentricities and inclinations is exact, as compared to a nondecisive increase in precision on the orbit precession frequencies. Moreover, we know that the Nbody 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 semimajor axes in Eqs. (3) and (13), the secular Hamiltonianof the eight Solar System planets considered in our study reads (14)
2.2 Forced inner planets
As discussed in the Introduction, the secular Hamiltonian (14) would not correctly reproduce the frequencies of the Jupiter and Saturndominated 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 orbits (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 and for all, an explicit quasiperiodic time dependence for the orbits of the giant planets, that is, by expressing the corresponding Poincaré complex variables as finite Fourier series, (15)
where t denotes the time, and ỹ_{kℓ} are complex amplitudes, m_{kℓ} and n_{kℓ} are 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 inthe decompositions depends on the planet. When this predetermined time dependence is injected into Eq. (14), we obtain the Hamiltonian of a forced secular inner Solar System, (16)
The explicit time dependence in the Hamiltonian 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 possesses eight degrees of freedom. As already stated, compared to the secondorder secular system of Laskar (1984, 1985), our 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 in eccentricities or 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 our model perfectly suited to exploring the very longterm evolution of the inner system, even when highly eccentric and inclined orbits become statistically recurrent.
2.3 Construction of the outer planet solution. Initial conditions
The quasiperiodic form of the outer orbits in Eq. (15) was established numerically and is explicitly reported in Appendix D. The full equations of motion of the main bodies of the Solar System were numerically integrated over 30 Myr in the future, following the comprehensive model of Laskar et al. (2011). The initial conditions of the integration were adjusted through least squares to the highprecision planetary ephemeris INPOP13b (Verma et al. 2014; Fienga et al. 2014) that extends over 1 Myr. Throughout the paper, we refer to this direct numerical integration as LaX13b. Through frequency analysis (Laskar 1988, 1993, 2005; Laskar et al. 1992) of the orbital solution, the leading secular harmonics of the dimensionless Poincaré variables of the outer planets were extracted to construct the Fourier series (15) as in Laskar (1988, 1990c). This was performed up to a numerical precision that no longer allows recognizing new harmonics as linear combination of the fundamental frequencies ω_{o} in an unambiguous way. The precision of the outer planet solution is shown in Table D.4. We report there the root mean square of the dimensionless Poincaré variables of the outer planets in the LaX13b solution, and that of the corresponding residuals after subtraction of their quasiperiodic decomposition from Eq. (15) and filtering out of the nonsecular Fourier components. The effectiveness of the quasiperiodic approximation is illustrated in Figs. D.1 and D.2.
The choice of the secular semimajor axes and that of the initial conditions for the secular Poincaré variables of the inner planets deserves 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 because this might cause an offset in the secular frequencies of the motion as a result of the shorttime (highfrequency) oscillations generated by the Lie transform. This is a general rule when averaged dynamical systems are constructed (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 shorttime components of the numerical solution of the nonaveraged system, as we did here. The constant term in the frequency analysis of the variables Λ_{k}, performed on the LaX13b solution, provides the value of the secular semimajor axes of the planets (Table D.1). Moreover, a polynomial expansion in time over the first few thousand years of the same solution, once filtered, provides the nominal initial conditions for the secular Poincaré complex variables of the inner planets (Table D.2).
3 Numerical integration. Gauss’s method in Hamiltonian formalism
The dimensionless Hamilton equations for the forced inner system (16) read (17)
for k ∈{1, 2, 3, 4}. Because of the averaging over the mean longitudes in Eq. (14), these equations include double integrals of the form
where stands for or , alternatively. The equations of motion (17) are thus nonalgebraic, 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, 1990c). Nevertheless, it is well known that the Eqs. (17) can be integrated numerically in a very efficient way through the socalled Gauss method (Gauss 1818; Bour 1855; Hill 1882). 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 elliptic integrals (e.g. Olver et al. 2020, Chapter 19). This was first shown by Gauss (1818), who at same time introduced the arithmeticgeometric mean to numerically evaluate these 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 our 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 guaranteeing high numerical precision. We employed a dedicated algorithm based on piecewise minimax rational function approximation to numerically compute the complete elliptic integrals in double floatingpoint precision at the cost of elementary functions (Fukushima 2015). The remaining simple integral in Eq. (18) was then effectively evaluated with 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 applied the trapezoidal rule in an adaptive way by doubling the number of function evaluations until the estimated relative error on the integral was smaller than 10^{−12}. Because the convergence of the numerical integral is exponential, the resulting error is often orders of magnitude smaller than this tolerance.
We integrated Eqs. (17) using an Adams PECE method of order 12 and a conservative time step of 250 yr, as in Laskar (1994). In absence of any fundamental integral of motion (see Sect. 2.2), we estimated the integration error after a time T following Laskar (1994). Starting from the nominal initial conditions, we integrated 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 , we defined the relative integration error as (20)
where 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 , with δ_{0} = 9 × 10^{−15}, T_{0} = 0.6 Myr and α = 1.2. Therefore we have δ(10 Myr) ≈ 3 × 10^{−13}, and δ(10 Gyr) ≈ 10^{−9} when we extrapolate by ignoring the chaotic behaviour of the solution. The integration error is therefore similar to that of the secular system in Laskar (1994). We emphasize that our fixedtimestep integration scheme allows us to reach very highly excited orbital states. Violations of the precision goal on the integral (18) typically occur only when the system is already close to the first intersection of the instantaneous Keplerian orbits of two planets (see Sect. 7), and they 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 our dynamical system because the derivative evaluation constitutes by far the most expensive step of the integration scheme. When it is implemented in the C programming language, Gauss’s method allows us 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 is 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 to adapt Gauss’s method to our complex Hamiltonian formalism is the subject of a forthcoming paper (Mogavero, in prep.), along with the release of the C program implementing it.
4 Analytical expansion. Computer algebra
In addition to its fast numerical integration, our dynamical model allows a systematic analytical development of its Hamiltonian by means of computer algebra. We employed TRIP, a computer algebra system dedicated to perturbation series, especially those of celestial mechanics (Laskar 1990a; Gastineau & Laskar 2011, 2020). The main objects of its symbolic kernel are the Poisson series, that is, multivariate Fourier series whose coefficients are multivariate Laurent series, (21)
where and are complex and real variables, respectively, , and C_{k,ℓ} are complex coefficients.
Following Laskar & Robutel (1995), the low eccentricities and inclinations in the Solar System can be exploited to develop the averages appearing in Eq. (14) as formal series in Poincaré’s variables . Eq. (4) gives (22)
The rotational invariance of the secular Hamiltonian (14) requires , while from planar symmetry, itfollows that is an even integer. These relations imply that the monomials in Eq. (22) are even with respect to the and (y, y′, ȳ, ȳ′) variables separately. The relativistic terms in Eq. (14) are readily expanded, (23)
By truncating the series (22) and (23) at a given total degree 2n () in Poincaré’s variables and after substitution in Eq. (14), we obtain a polynomial secular Hamiltonian for the ensemble of the Solar System planets, (24)
where groups all the monomials of same total degree 2p. 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 semimajor axes, which are constant parameters in our dynamics. They can therefore be numerically evaluated once and for all, to obtain very compact series. This is shown in Table 2, which gives the total number of monomials in Eq. (24) with respect to the degree of truncation, along with the disk usage of the corresponding series when stored in a plaintext file, and their typical computation time in TRIP.
The expansion (24) straightforwardly provides a truncated Hamiltonian for the forced inner system (16), (25)
The computation of the equations of motion for the truncated Hamiltonian (25) only requires taking the derivatives of a multivariate polynomial, and can thus be systematically performed in TRIP. We set the resulting polynomials in Horner form to ensure speed and stability of their numerical evaluation, and wrapped them into a C code to achieve the best computational performance^{6}. The equations of motion were then integrated through the Adams PECE method of order 12 employed for Gauss’s dynamics, with the same time step of 250 years. The typical computational cost of a 100 Myr orbital solution is shown in Table 1, listed according to the degree of truncation. In light of the pertinence of the forced inner system that we shown in Sect. 5, the truncated Hamiltonian (25), at degrees 4and 6 in particular, is the best of very fast and still realistic dynamical models of the inner Solar System.
4.1 Forced LaplaceLagrange dynamics
For low eccentricities and inclinations, the secular Hamiltonian (14) is a perturbation of the integrable Laplace–Lagrange (LL) problem. Its truncation at degree 2 in Poincaré’s variables gives the quadratic form (26)
where we have defined the column vectors and , T is the transposition operator, and the dagger stands for the Hermitian transposition, that is, . The matrices M and N are real and symmetric. Because 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 substituting the time dependence of the outer planet variables given in Eq. (15), we obtain (27)
where we have defined the inner and outer planet column vectors, and , respectively, and the real matrix M is written as a block matrix, with M_{ii} and M_{oo} symmetric 4 × 4 matrices. When terms that only depend on time are discarded, the Hamiltonian reads (28)
where Re stands for the real part of a complex quantity. The Hamiltonian (28) corresponds to a forced Laplace–Lagrange system. Because it is quadratic in the inner vector x_{i}, we can analytically solve the corresponding equations of motion and introduce appropriate actionangle 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}, (29)
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 eigenvalues^{7}, g_{LL} being a column vector. The orthogonal matrix O_{M} induces the canonical change of variables^{8} defined as (30)
The transformed Hamiltonian reads (31)
with . The corresponding Hamilton equations are given by (32)
and constitute a firstorder inhomogeneous matrix ordinary differential equation. The general solution can be written as (33)
The free solution is the general integral to the associated homogeneous equation , representingthe autonomous perihelia precession of the inner orbits, while the forced solution is a particular integral of the complete Eq. (32), which arises from the gravitational forcing of the outer planets. We define, once and for all, the forced solution as (34)
When the decomposition given in Eq. (15) is employed, the components of the forced solution are (35)
The (constant) denominators appearing in Eq. (35) are far from zero because the inner planets are not involved in the corresponding secular resonances. The forced solution is thus well defined.
4.2 Proper modes
The following derivation shows that a canonical transformation exists that depends on time, , such that the transformed LL Hamiltonian reads . As in the case of an autonomous Laplace–Lagrange system, the new canonical variables u will physically correspond to the free part of the solution in Eq. (33) (e.g. Morbidelli 2002). We thus begin by defining a new set of variables u such that (36)
By using the fact that is a solutionto Eq. (32), and discarding terms only depending on time, the forced LL Hamiltonian in Eq. (31) can be written as (37)
We now ask that the change of variables derive from a timedepending generating function that satisfies (38)
so that the last two terms in Eq. (37) cancel out in the transformed Hamiltonian. By integrating the above equation with respect to time, we obtain (39)
where is an unknown function that has to be determined. Because the generating function F depends on the old momenta and the new coordinates − ju^{†}, it must verify the relations (e.g. Landau & Lifshitz 1969) (40)
With the choice , these relations are both equivalent to Eq. (36).
The previous derivations imply that the variable transformation is canonical, with
We call the complex variables u the proper modes of thePoincaré’s variables x_{i}. In the forced LaplaceLagrange dynamics, they simply rotate in the complex plane, with constant angular frequencies given by g_{LL}. According to Eq. (42), the corresponding Poincaré variables thus result in a superposition of these independent harmonic oscillations and those arising from the forcing of the outer planets. When higherdegree terms of the Hamiltonian (25) are taken into account, the dynamics of the proper modes becomes coupled, and their frequency spectrum is no longer monochromatic. Physically, the proper modes u correspond to the variables defined in Laskar (1990c). Nevertheless, it is important to emphasize that their mathematical definitions do not coincide. The transformation matrix O_{M} does not include contributions at the second order in planet masses, as that considered in Laskar (1990c) does. Moreover, transformation (42) is nonlinear in the forcing of the outer planets, in the sense that harmonics of order^{9} m greater thanone are present in the decomposition (15) (see Table D.5) and thus in the forced solution (35).
Similarderivations allow us to define the proper modes v of the Poincaré variables y_{i} by means of a timedependent canonical change of variables ,
the columns of O_{N} being the eigenvectors of N_{ii}, while the diagonal entries of D_{N} = −diag(s_{LL}) are the corresponding real eigenvalues. Using the quasiperiodic decomposition in Eq. (15), the components of the forced solution are (47)
where .
When the proper modes u, v have been defined, it is very useful to establish a corresponding truncated Hamiltonian of the form (48)
where is the LL Hamiltonian. To define higherdegree truncations, we first considered the nonquadratic part of the secular Hamiltonian (14), truncated at a given total degree in Poincaré’s variables , (49)
Then, the proper modes u, v were injected in (49) by means of Eqs. (42) and (44), and the Poincaré variables of the outer planets were replaced by their quasiperiodic decompositions given in Eq. (15). However, these substitutions do not conserve the degree of the terms in the expansion because nonlinear harmonics are present in Eqs. (15), (35), and (47), that is, harmonics of order m or n greater than one. These harmonics are in principle much smaller than the linear ones and generate higherdegree terms in the substitution process. To define the terms in a consistent way, we therefore introduced a fictitious real variable ϵ to redefine the quasiperiodic decompositions (15) as (50)
meaning that each harmonic is treated to be of degree m_{kℓ} or n_{kℓ} for the purpose of the expansion^{10}. The forced solutions (35) and (47) were modified accordingly. The series resulting from the substitution of variables in (49) can thus be truncated at total degree 2n with respect to . Finally, addition of the quadratic terms by definition gives the truncated Hamiltonian (48). Truncation is effectively performed in TRIP in parallel with variable substitutions thanks to a dedicated monomial truncated product (Gastineau & Laskar 2011), which allows us to minimize the computational cost of the expansion. It is worthwhile to note that at the end of the expansion process, the Hamiltonian terms that only depend on time can be removed because 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 in the Hamiltonian according the degree of truncation. We also report the typical computation time of the series in TRIP, related to the transformation . We note the dramatic increase in the number of terms caused by this change of variables. We point out that the precise number of monomials in depends on the particular quasiperiodic decomposition (15) we used in this work (see Appendix D).
4.3 Actionangle variables
Actionangle variables corresponding to the proper modes u, v are introduced by standard canonical transformations (u, −ju^{†}) → (X, χ) and (v, −jv^{†}) → (Ψ, ψ) defined as (51)
for k ∈{1, 2, 3, 4}. The Laplace–Lagrange Hamiltonian thus reads and it is trivially integrable. These transformations allow us to expand the Hamiltonian as a Fourier series of the angle variables χ, ψ and ϕ(t) = ω_{o}t, (52)
where are complex amplitudes, and we employed a compact notation for the actionangle variables, (53)
Only a finite number of harmonics have nonzero amplitude in Eq. (52) because we deal with a truncated Hamiltonian. We also recall that because the Hamiltonian is a real function, we have for all k, ℓ. In Table 2, we show the number of harmonics in Eq. (52) according to the truncation degree of the Hamiltonian, that is, the number of different wave vectors (k, ℓ) up to a global minus sign.
The explicit time dependence of the Hamiltonian can be easily absorbed in a phasespace extension through the definition of actionangle variables for the trivial degrees of freedom of the outer orbits. We introduce (Φ, ϕ) such that the new Hamiltonian reads (54)
with . The dynamics of the additional angles is thus consistently given by , while that of the actions Φ is irrelevant. This autonomous formulation is useful in the context of canonical perturbation theory.
Frequency in arcsec yr^{−1} of the largestamplitude Fourier harmonic for each proper mode of the Poincaré variables over the timeinterval [0, 20] Myr.
5 Comparison with reference dynamical models
When averaged models are constructed, it is essential to compare the resulting trajectories to those of the nominal Hamiltonian or to the solutions of a more comprehensive dynamical model. This allows validating 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 was used to predetermine the secular motion of the outer planets (Sect. 2.3 and Appendix D). We also consider the secular solution of Laskar (1990c), denoted as La90 from now on, because it constitutes the most precise secular model to this day.
The comparison of the different models begins at short secular timescales in Fig. 1, where we show the inner planet eccentricities and inclinations over 250 000 yr in the future. The plots indicate that our 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 shorttime (orbital) oscillations of the direct integration, as the La90 solution does as well. In particular, the behaviour of the solutions at the origin of time shows that the initial conditions for the secular variables of this work are correctly determined. Some small periodic differences between thetwo secular models appear in some of the plots, with the La90 solution generally being more precise in reproducing the longterm average of the direct integration. However, these 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 our solution remains close to the direct integration. The solution La90 nearly coincides with LaX13b in these plots (see Laskar et al. 2004)and is therefore 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 noticeable. First of all, the orbital oscillations sometimes show 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 we consider 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 because it translates into 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 these deviations. To this end, we performed a frequency analysis of the proper modes of the Poincaré variables for the different orbital solutions over the [0, 20] Myr time interval. For solutions LaX13b and La90, we employed the proper modes defined in Laskar (1990c), and for our solution we used the proper modes (using gives the same results at the numerical precision of 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 ours is a few hundredths of an arcsecond per year, with a maximum of 0.06′′ yr^{−1} in the case of s_{3} and s_{4}. This agrees with the behaviour shown in Fig. 2 and the general expectation from the discussion of the model precision in Sect. 2.1. When combinations of the fundamental frequencies (which determine the resonant structure) are considered, the differences between the two models can be even smaller. 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 are only 0.01′′ yr^{−1} and 0.004′′ yr^{−1}, respectively.
Table 3 also reports the fundamental frequencies of the truncated model (Eq. 25) for different total degrees of truncation. When combinations of frequencies (g_{1} − g_{5}) − (s_{1} − s_{2}) and 2(g_{3} − g_{4}) − (s_{3} − s_{4}) are considered, models and almost show the same deviations from solution LaX13b as the nontruncated forced model . 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 models and nearly coincide with that of the nontruncated forced model over the first 20 Myr.
Fig. 1 Eccentricities (left side) and inclinations (right side) over the invariant plane J2000 (see Appendix D.2) of the inner planets over 250 000 yr in the future. The solid black line stands for the full (i.e. nonfiltered) direct integration LaX13b, while the red line is the solution of our forced model in Eq. (16). The secular integration La90 is represented by the dashed blue line. The maximum relative root mean square of our model residuals is 2% for the eccentricities and 0.7% for theinclinations. 
Fig. 2 Eccentricities (left side) and inclinations (right side) over the invariant plane J2000 (see Appendix D.2) of the inner planets over10 Myr in the future. The solid black line stands for the full (i.e. nonfiltered) direct integration LaX13b, and the red line is the solution of our forced model in Eq. (16). The secular solution La90 nearly coincides with the direct integration in these plots and it is therefore not shown. 
6 Finitetime maximum Lyapunov exponent
When the adequacy of our dynamical model has been established over a few million years, we have to assess whether it is able to correctly reproduce the resonant structure of the inner Solar System over longer times. Even though the fundamental frequencies in Table 3 suggest that this is the case, the computation of the maximum Lyapunov exponent (MLE) remains an essential test because 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 nonnull probability of an unstable evolution of the inner planets (Laskar 1994, 2008; Batygin & Laughlin 2008; Laskar & Gastineau 2009) prevents the existence of the MLE as an infinitetime limit, as in its usual mathematical definition (Oseledec 1968); the consideration of a finitetime MLE (FTMLE) is pertinent to our case. Because we do not have an efficient numerical algorithm for evaluating the variational equations of Gauss’s dynamics atour disposal, we considered the FTMLE given by the twoparticle algorithm proposed in Benettin et al. (1976), which was 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) again brings the shadow particle at a distance d_{0} from the reference particle. The resulting FTMLE depends in an 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 computing it acquires full physical significance only for an ensemble of trajectories.
In Fig. 3, we show the computation of the FTMLE of the forced inner Solar System over 5 Gyr in the future for an ensemble of 1148 stable orbital solutions with initial conditions very close to the nominal ones (see Sect. 7 for the definition of an unstable solution in the framework of our model). On the left vertical axis, we report the FTMLE expressed as an angular frequency in arcsec yr^{−1}, while on the right axis, we show the corresponding Lyapunov time, defined as FTMLE^{−1}, 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 FTMLE, which by definition encloses 50% of the solutions around the median (shown by the red line). In the same manner, the light grey region corresponds to the [1th, 99th] percentile range of the PDF, enclosing 98% of the probability. The two FTMLEs 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 conditions. They were chosen by taking the relative variation of each coordinate of the nominal phasespace 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 phasespace position of the system, with a diagonal covariance matrix. The width of this distribution is of the same order of magnitude as that in Laskar (2008, Table 1). The initial position of each shadow particle around the corresponding reference particle was then chosen according to the same kind of Gaussian distribution, with a relative standard deviation of 10^{−8}. Such a value shouldminimize the accumulation over long timescales of numerical errors, due to the floatingpoint implementation of the algorithm ofBenettin et al. (1976), when working in double precision (Tancredi et al. 2001; Mei & Huang 2018). In Fig. 3 the computation of the FTMLE 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 our FTMLEs are asymptotically independent of d_{0}, the asymptotic regime being practically reached in a few hundred million years. In the same manner, our computation was tested to be asymptotically independent of the renormalization time τ, set to 5 Myr in our computation, and of the norm chosen for the phasespace vectors, taken here to be Euclidean as usual.
For times shorter than 100 Myr, Fig. 3 shows that the distribution of the FTMLE only reflects the choice of different initial tangent vectors d_{0} around essentially the same reference trajectory, that is, the nominal trajectory. Therefore this part of the plot can be discarded because of its nonasymptotic character. We also note that for t < 100 Myr, the PDF of the FTMLE tends to shrinks with increasing time, as expected by the asymptotic behaviour of the method (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 FTMLE, which chaotically depends on its initial conditions. The timeasymptotic regions shown in Fig. 3 thus represent the probabilistic knowledge of the FTMLE 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 estimate (Rosenblatt 1956; Parzen 1962) of the PDF of the Lyapunov time FTMLE^{−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 to select the optimal bandwidth (Silverman 1986). We also show through the light blue region the pointwise confidence interval at the 98% level of the estimated PDF, obtained through nonparametric bootstrap, that is, by resampling with replacement of the original data (Efron 1979). The value of about 5 Myrtypically reported in literature (Laskar 1989; Sussman & Wisdom 1992) is right in the bulk of our 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 thatthe asymptotic PDFs of the two secular models should largely overlap. Finally, we note that the value of 1.1 Myr numerically found in Batygin et al. (2015) as the Lyapunov time of a simplified Mercury dynamics does not agree with our findings. Figures 3 and 4 show that even considering the diffusion of the planet orbits over 5 billion years, a high FTMLE like this is practically never reached by our dynamical model, even though the latter includes all the Hamiltonian harmonics taken into account in Batygin et al. (2015).
Fig. 3 Finitetime maximum Lyapunov exponent (FTMLE) and corresponding Lyapunov time FTMLE^{−1} of the forced inner Solar System over 5 Gyr from an ensemble of 1148 stable orbital solutions with very close initial conditions. 
Fig. 4 Probability density function of the Lyapunov time FTMLE^{−1} at 5 Gyr. The kernel density estimate is shown as a dark blue line and the corresponding CDF is shown on the upper horizontal axis. The pointwise confidence interval from bootstrap is at the 98% level and is shown as a light blue region. 
7 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 with the aim to assess if our model is able to reproduce the rate of the high Mercury eccentricities observed in Laskar (2008) and 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, that is, the intersections of each orbit with the orbital plane of the other planet. When two orbits cross, a pair of mutual nodes exchanges their positions along the line of nodes (i.e. the intersection of the two orbital planes). This computation is inexpensive in the context of Gauss’s dynamics and can be performed at each time step of the integration scheme. When a secular collision occurs, our averaged model is no longer a good longterm approximation of the original Nbody 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 solution. The statistics of the firstcollision time τ_{c} should 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 are not involved in extrapolating our secular dynamics out of its validity range, and the use of a fixedtimestep integration scheme (see Sect. 3) turns out to be sufficient and even suitable. We note that a secular collision of 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 performed an ensemble of 10 560 numerical integrations of the forced inner system over 5 Gyr. The initial conditions were chosen according to the same multivariate Gaussian distribution centred at the nominal values as in Sect. 6. The relative standard deviation of each phasespace coordinate is 10^{−9}, which is the same order of magnitude as that in Laskar (2008, Table 1). Table 4 shows the observed probability of having a maximum Mercury eccentricity greater than a given value over the entire time. The statistics of our dynamical model is compared to the statistics resulting from Laskar (2008, Table 3) and Laskar & Gastineau (2009, see Table C.1 of this paper). The statistical bounds represent the Wilson (1927) score interval, corresponding to a 98% confidence level. At our 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 Nbody dynamics, even if they are somewhat lower. This agrees with the general expectations from an averaged model, which has fewer degrees of freedom (Laskar 1994, 2008). When the midpoints of the statistical intervals are compared, Table 4 might also suggest that the true probabilities of our model are slightly lower than those in Laskar (2008). If this is true, it would agree with the expectation of a higher regularity of the forced model because the degrees of freedom dominated by the outer planets (their fundamental frequencies in particular) are frozen into a quasiperiodic time dependence. Moreover, the Hamiltonian terms at second order in planet masses could also contribute to somewhat larger chaotic excursions of the inner orbits.
Confidence interval at the 98% level of the probability P(sup_{t≤5 Gyr}e_{1}(t) > e_{max}) in percent, where e_{1} is the eccentricity of Mercury.
Fig. 5 Evolution of the fundamental frequency g_{1}, as defined by Eq. (56), in the 42 collisional orbital solutions (top panel). Six solutions are isolated in the bottom panel. The horizontal line stands for the constant fundamental frequency g_{5}. 
7.1 Secular collisions
In 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), six collisions between Mercury and the Sun, nine collisions between Mercury and Venus, and one collision between Mars and the Sun are found (see Appendix C). These 16 events correspond to a 98% confidence interval (0.36%, 1.13%) for the probability of a physical collision. Our definition of a secular collision thus provides a very good estimate of the probability of a physical collision in the nonaveraged 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 pair, with a maximum Mercury eccentricity ranging from 0.797 to 0.966. In particular, in 41 out of 42 collisional solutions, it is higher 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 might not correspond to a physical collision in the nonaveraged dynamics. We emphasize that as a consequence, the statistics of the secular collisions other than that of the Mercury–Venus pair cannot be correctly reproduced by our numerical computation because events like this, as observed in Laskar & Gastineau (2009), always require a very high 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 first been established numerically (Laskar 2008; Batygin & Laughlin 2008) and was then confirmed by an analytical study (Boué et al. 2012). We conclude this section by showing that the g_{1} − g_{5} resonance isindeed reached by the collisional orbital solutions of our ensemble. While g_{5} ≈ 4.257′′ yr^{−1} is constant in our dynamical model, the frequency g_{1}, which dominates the spectrum of the Poincaré 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é variables x with a time step Δt = 1000 yr, that is, we consider the time series x(t_{n}) with t_{n} = nΔt and . 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, that is, as a continuous function of time. The instantaneous (angular) frequency ω_{1} of the proper mode u_{1} is given by (Cohen 1995) (55)
We could analytically compute the derivative of Eq. (42) to obtain the time series ω_{1} (t_{n}). In practice, we employed numerical differentiation via the Lagrange fivepoint formula (Olver et al. 2020) because very high numerical precision is unnecessary here. As it is, the frequency ω_{1} (t_{n}) presents shorttime oscillations that hide its longterm evolution. To suppress these highfrequency fluctuations, we therefore defined the time series g_{1}(t_{n}) as the output of the lowpass Kolmogorov–Zurbenko (KZ) filter (Yang & Zurbenko 2010, and references therein) applied to ω_{1} (t_{n}), (56)
The KZ filter is defined as an iteration of the common moving average and is comprehensively characterized in Appendix B. For our application, we used 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} in 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} in 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 halfwidth of more than 1′′ yr^{−1}, which has been reported in Laskar (2008, Fig. 6(a2))^{12}. In 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.
8 Numerical experiment over 100 Gyr
In this last section, we present the results of a numerical experiment in which we integrated an ensemble of orbital solutions of the forced inner system over 100 Gyr in the future. Through these simulations, we aim to characterize the type of stochastic process that arises from the Hamiltonian (16) and drives the destabilization of the inner system in a regime where the highly excited orbits no longer represent rare events, as they do in Sect. 7. This type of simulation clearly does not correspond in any way to the real future evolution of the orbits of the inner planets over much more than 5 Gyr. 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 semimajor 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 longterm stochastic behaviour of the dynamical system (16), without inferring any conclusion on the real evolution of the inner system over these timescales.
We numerically integrated an ensemble of 1042 orbital solutions of the forced inner system over 100 Gyr. The initial conditions were 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 was 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. 6represents the PDF estimated through the kernel density estimation (KDE) method (Rosenblatt 1956; Parzen 1962). We employed a standard Gaussian kernel and Silverman’s rule of thumb to select the optimal bandwidth (Silverman 1986). The PDF was normalized in such a way that its integral over 100 Gyr is equal to the overall percentage of secular collisions obtained in the experiment^{13}, that is, 890∕1042 ≈ 85.4%. The blue region in the lower panel of Fig. 6 represents the pointwise confidence interval for the PDF at the 98% level, obtained via nonparametric bootstrap, that is, by resampling with replacement of the original data (Efron 1979). This region thus corresponds pointwisely to the [1th, 99th] percentile range of the bootstrapped estimated PDF. We emphasize that for simplicity, we did not correct for the bias arising in the kernel density estimation, so that strictly speaking, we show a confidence region for the expectation of the KDE and not for the true subjacent PDF (Chen 2017). When referred to the true PDF, this confidence region typically represents an undercoverage at the points where the function bends, that is, 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 that we deal with over the first few billion years. These corrections of the confidence intervals are sufficiently small, however, to not affect at all the conclusions we infer in this section. To give a further idea of the statistical variance, in Fig. 6 we also show in red 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.
Fig. 6 Kernel density estimation of the PDF of the time τ_{c} of the first secular collision from 1042 integrations of the forced inner system over 100 Gyr (blue line in the lower panel). The estimate of the CDF is shown as a blue line in the upper panel. The blue regions represent the pointwise confidence intervals at the 98% level from bootstrap. In red we show the PDF estimate over 5 Gyr from the 10 560 orbital solutions of Sect. 7. The dashed black lines stand for the analytical PDF (57) with T_{1} = 27.6 Gyr and α = 0.9 and the corresponding CDF. 
8.1 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 metastability that characterizes the inner Solar System according to the definition of Laskar (1996). Over a timescale comparable to its age, that is, 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 this probability increases very fast with time, so that the system rapidly becomes statistically unstable over longer timescales, meaning that unstable evolutions no longer represent rare events. In particular, the probability of an instability is about 20% already at 20 Gyr and is 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 ofmarginal stability, that is, it should be practically stable over a timescale comparable to its age, while stronginstabilities arise very fast over longer times, resulting in more stable configurations of the surviving planets.
8.2 Destabilization of the inner Solar System as a stochastic process
In the framework of the simplified model introduced in Batygin et al. (2015), Woillez & Bouchet (2020) recently applied the theory of the whitenoise limit for slowfast dynamical systems (e.g. Gardiner 1985) to describe the longterm dynamics ofMercury through an effective stochastic process. They showed that the instability timescale of the Mercury orbit is reproduced by the hitting time of a onedimensional Wiener process with a reflecting barrier. The corresponding PDF peaks at 1.3 Gyr (Woillez & Bouchet 2020, Fig. 4), in agreement with the analysis in Batygin et al. (2015). This timescale is an order of magnitude shorter than that in Fig. 6, and the simplified dynamics of Mercury predicts a high 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 the longterm wanderings of the fundamental frequency g_{1} towards the destabilizing secular resonance g_{1} − g_{5} (see Sect. 7.2) in an approximate way. Numerically, the PDF of g_{1} is asymmetric, with a longer tail that diffuses towards low values, while its upper end is nearly fixed at about 6′′ yr^{−1} (Hoang et al. 2021, see also Fig. 5). Therefore we assumed that g_{1} performs a random walk described by a Wiener process starting at g_{1,0} = g_{1}(t = 0), and considered an upper reflecting barrier at g_{1,max} > g_{1,0}. We assumed that a secular collision rapidly occurs whenever the secular resonance g_{1} − g_{5} is reached, as shown in Fig. 5; this was taken into account by considering an absorbing barrier at g_{1} = g_{5}. The PDF ρ(τ) of the hitting time is thus given by (e.g. Schwarz 1992; Woillez & Bouchet 2020) (57)
where α = (g_{1,0} − g_{5})∕(g_{1,max} − g_{5}) and , with D being the diffusion coefficient ofthe Wiener process. After a leastsquares search, we plot in Fig. 6 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 our statistical precision and apart from a certain excess of events in the tail of the distribution^{14}. On the one hand, this shows that at short times, the PDF behaves as (58)
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 decay of the density probability for τ → 0 in Eq. (58) characterizes the rare destabilizations of the Mercury orbit over the first few billion years. In particular, the instanton phenomenology described in Woillez & Bouchet (2020) could indeed apply to the real Solar System (this scenario still needs to be tested through integration of realistic models such as ours.). On the other hand, Eq. (57) shows that for τ → +∞ 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 because 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 that is more than ten times lower than that in Fig. 6, that is, 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^{15}. As discussed in Sect. 7, averaged models generally tend 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 that are ten times faster than in our dynamics, even though our model includes all the dynamical interactions considered in the simplified model. This 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 the Mercury orbit can be estimated from the maximum Lyapunov exponent (Fig. 3) if we assume that no dynamical constraints prevent g_{1} from diffusing 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 resonances affecting the frequency g_{1} (Laskar 1990c). We may estimate an upper bound on the diffusion coefficient D as (59)
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 thusgiven by (60)
with g_{1,0} = 5.577′′ yr^{−1} (Table 3) and g_{5} = 4.257′′ yr^{−1}. Equation (60) gives essentially the same diffusion timescale as in Batygin et al. (2015) and Woillez & Bouchet (2020). This is already two orders of magnitude greater than the Lyapunov time of the system, but it is still insufficient to explain the findings of Fig. 6. This means that some dynamical constraints effectively cause the diffusion towards the g_{1} − g_{5} resonance to be slower than in the estimates (59) and (60). Certain secondary resonant harmonics, in particular, must play a determining role in the destabilization of the inner Solar System.
9 Conclusions
We introduced 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 low masses of the inner planets, and the absence of strong meanmotion 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, that is, 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 such as 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 our 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 mechanism 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 after it is truncated at degree 4 in planet eccentricities and inclinations by using the corresponding expression of the twobody secular Hamiltonian reported in Laskar & Robutel (1995, Appendix) and the quasiperiodic secular solution for the outer planets given in Appendix D.
Acknowledgements
FM acknowledges the invaluable support of M. Gastineau in the implementation of the algebraic computations in TRIP and the fruitful discussions with N. H. Hoàng. FM has been supported by a PSL postdoctoral fellowship and by a grant of the French Agence Nationale de la Recherche (AstroMeso ANR19CE31000201). This project has been supported from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo885250). This work was granted access to the HPC resources of MesoPSL financed by the Région ÎledeFrance and the project Equip@Meso (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
Appendix A Indirect part of the twobody Hamiltonian perturbing function
We outline the derivation of the coefficients of the indirect part of the twobody perturbation in the Fourier expansion (5) by adapting the presentations of (Laskar 1991; Laskar & Robutel 1995; Laskar & Boué 2010). From Sect. 2, we have (A.1)
where denotes the velocity vector tangent to the instantaneous Keplerian ellipse of a given planet. In the reference frame of the Keplerian orbit, with the origin at the Sun, the xaxis directed towards the pericenter, and the zaxis normal to the orbital plane, the components of the vector v read (A.2)
where E is the eccentric anomaly, is the mean motion, and r = r is the heliocentric orbital distance. After transformation to the fixed reference frame, the velocity components are given by (A.3)
where the rotation matrices and are given in (Laskar 1991; Laskar & Robutel 1995).
We introduce the eccentric longitude F = E + ϖ, which obeysthe modified Kepler equation (A.4)
where z = eE^{jϖ} and . By using the fact that , and after some algebra, we find (A.5)
By employing the complex formalism of (Laskar & Robutel 1995), we can thus write (A.6)
The application of matrix in Eq. (A.3) then gives (A.7)
with and . The scalar product in Eq. (A.1) thus reads (A.8)
The coefficients are readily derived afterthe quantities a∕r and (a∕r)E^{jF}, appearing in Eq. (A.6), are expanded in Fourier series of the mean longitude λ. We have (A.9)
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, (A.10)
can be straightforwardly expanded in series of and , and is well defined at zero eccentricity. In a similar manner, we obtain (A.11)
Appendix B KolmogorovZurbenko filter
Lowpass filters are employed in timeseries analysis to extract the longterm (lowfrequency) components of a given signal. A dedicated lowpass filter has been constructed in (Carpino et al. 1987) to recover the secular changes of the giant planet orbits from the numerical output of an Nbody integration. Here we propose the use of the KolmogorovZurbenko filter (Yang & Zurbenko 2010) as an outofthebox, computationally advantageous, and still effective choice.
We consider the realvalued time series (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 discretetime Fourier transform (B.2)
which is a periodic function with period ω_{s} and can thusbe restricted to the interval ω≤ ω_{s}∕2. The action of a finiteimpulseresponse filter F on the time series is defined as the discrete convolution of the signal ξ_{n} with a given finite sequence d_{n}, (B.3)
M being a nonnegative 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 discretetime Fourier transform of its impulse response, which is called the frequency response of the filter, and gives the spectrum of the output signal as (B.4)
Ideally, a lowpass filter should be characterized by a null response above a given cutoff frequency, that is, for ω > ω_{c}, while the response should be unitary at lower frequencies, that is, for ω ≤ ω_{c}. Moreover, the filter must not change the phase of the Fourier components of the signal, that is, . This is achieved by requiring (Carpino et al. 1987). Nevertheless, a real filter is rather characterized by a passband, for ω≤ ω_{p}, ρ being a given maximum loss; a stopband, for ω≥ ω_{c}, α being a given maximum gain; and a transition band between ω_{p} and ω_{c}, where the response function smoothly decreases from ≈1 to ≈0.
Fig. B.1 Impulse response of the KolmogorovZurbenko filter (the moving average corresponds to the case k = 1). The axes have been normalized to show the asymptotic Gaussian behaviour of the function. The standard normal distribution N(0, 1) is represented by the solid black curve. 
One of the simplest lowpass filters is the moving average (MA), (B.5)
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 this problem by applying the moving average iteratively (Yang & Zurbenko 2010, and reference therein). The KolmogorovZurbenko filter is thus defined as (B.6)
The output signal of this filter is given by (B.7)
where is the number of ways of choosing k integers in the interval [−M, M] such that their sum is equal to s. The numbers can be expressed as the coefficients of the finite Laurent series (Yang & Zurbenko 2010), that is, (B.8)
The impulse response (IR) of the KZ filter given in Eq. (B.7), , can be interpreted as a discrete probability distribution function (PDF) over the interval [−kM, kM], resulting from the convolution of k uniform distribution over [−M, M]. The central limit theorem implies that for k ≫ 1, this PDF isasymptotically Gaussian, with zero mean and variance equal to , being the variance of the discrete uniform distribution over [−M, M]. This is shown in Fig. B.1. The width of the impulse response therefore scales as for M, k ≫ 1, even if the filter length L_{kM} is linear in k.
Fig. B.2 Absolute value of the frequency response for the KolmogorovZurbenko filter (the moving average corresponds to the case k = 1). The vertical dashed black line stands for the cutoff frequency f_{c} given by Eq. (B.14), while the horizontal coloured lines represent the maximum gains in the stopband α_{k} = η^{k}. The vertical dotted coloured lines show the passband frequencies f_{p} given by Eq. (B.17) for a loss level of 5%, which is represented by the horizontal dotted black line. 
B.1 Frequency response
The following derivations describe the behaviour of the KZ filter in the M ≫ 1 regime. Its frequency response (FR) is simply the k^{th} power of that of the moving average, (B.9)
where f = ωΔt∕2π is the angular frequency in units of the sampling rate ω_{s}. The zeros ofthe frequency response are given by (B.10)
as shown in Fig. B.2. Following (Carpino et al. 1987), we define the dimensionless cutoff frequency f_{c} as (B.11)
where we have defined the maximum gain in the stopband , being the frequency of the first local maximum of FR_{M,k}(f) for f > 0. In the regime M ≫ 1, this frequency is given by the first positive solution x^{⋆} of the equation (B.12)
where x = πfL_{M}. Numerically, we find x^{⋆} ≈ 4.49. Therefore we obtain , independently of k, and the gain α_{k} = η^{k}, with η = sin(x^{⋆})∕x^{⋆} ≈ 0.217, which is thus asymptotically independent of M. These results are shown in Fig. B.2. The cutoff frequency f_{c} can be computed from Eq. (B.11) through the first positive solution x_{c} of the equation (B.13)
which is found numerically to be x_{c} ≈ 2.55. The cutoff frequency of the KZ filter is thus given by (B.14)
independently of M and k in the asymptotic regime M ≫ 1. Figure B.2 and Table B.1 show that with a small number k of iterations, the KZ filter is already able to provide a strong attenuation in the stopband. Given its straightforward numerical implementation through Eq. (B.6), the KZ filter is an effective general choice for a lowpass filter. However, choosing a very large number k of iterations would reduce the passband in the frequency response of the filter. We define the dimensionless passband f_{p} as (B.15)
where 0 < ρ ≪ 1 is a given maximum loss. In the regime M ≫ 1, by developing the frequency response in a Taylor series about x = 0, we have (B.16)
By neglecting terms of order , the pertinent solution is given by (B.17)
At first order in , we obtain (B.18)
which shows that the passband shrinks as k^{−1∕2}, as shown in Fig. B.2. This dependence on k is fortunately sublinear. Generally speaking, a reasonable choice of the number of iterations should be restricted in our opinion to 3 ≤ k ≲ 5, as suggested by Table B.1. We note that the resulting transition band can easily span more than one decade of frequency, and it is thus not narrow, as in more complex filters (e.g. Carpino et al. 1987). This is the main drawback of the simplicity ofthe KZ filter.
Frequency response of the KZ filter for M ≫ 1. Maximum gain in the stopband α_{k} = η^{k}. Ratio of the cutoff frequency f_{c} to the passband frequency f_{p} for two different maximum loss ρ.
Appendix C Additional statistics for (Laskar & Gastineau 2009)
In (Laskar & Gastineau 2009), 2501 numerical integrations of the full Solar System were performed over 5 Gyr using the SABA4 highorder symplectic integrator (Laskar & Robutel 2001). In most cases, the step size was constant, but for the integrations that included close planetary encounters or effective collisions, the step size was very much reduced in order to conserve the accuracy of the integration. This led to very long integration times for some of the runs, and not all the solutions were completed at the time of the publication of the paper. All 2501 runs were started on August 7, 2008. Of these, 2472 were completed in 2008, but 22 ended in 2009 and another 7 in 2010, the last was completed on June 15, 2010. We present here the full statistics resulting from these numerical integrations,which have since been presented in many conferences, but were never published. In Table C.1 we give the total number of solutions for which the eccentricity of Mercury reaches a given value over 500 Myr to 5000 Myr. This table is similar to Table S1 of (Laskar & Gastineau 2009), but provides the total number of solutions instead of the percentage. In the last column, we also provide the values reached after full completion of the numerical integrations in June 2010. All the other numbers were unchanged.
Number of solutions that reach a given eccentricity of Mercury (e_{m0}) over a given time (500, 1000, 1500, 2000, 3000, 4000, and 5000 Myr). The statistics are made over 2501 orbital solutions as publishedin (Laskar & Gastineau 2009). In the three 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. Of these 21 solutions, 6 present a collision of Mercury with Venus (i.e. their centre 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 (e.g. 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 less than 794 km), followed by a collision of Mars with the Sun.
Appendix D Quasiperiodic secular solution for the outer planets
We use a quasiperiodic 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 was deduced from a numerically integrated full solution of the Solar System. The solution of reference (LaX13b) used an identical model as in (Laskar et al. 2011), but was initially fitted to the improved INPOP13b highprecision planetary ephemeris (Verma et al. 2014; Fienga et al. 2014), extended over 1 Myr. The model was then integrated over 30 Myr to derive a quasiperiodic approximation through frequency analysis (Laskar 1988, 2005). The variables considered here are the Poincaré complex canonical variables in Eq. (2), scaled to suppress the semimajor axis (or Λ) dependence,that is, (D.1)
It should be noted that and where z = eE^{jϖ} and ζ = sin(i∕2)E^{jΩ} are the classical, noncanonical, complex elliptic elements.
The derivation of the QPSO model needs some care. We wish 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 meanmotion resonance among Jupiter and Saturn and the proximity of the 2:1 mean motion resonance in the UranusNeptune system, several shortperiod 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. Only the terms for which the angular argument is recognized in an unambiguous way as a combination of the fundamental frequencies (Table D.3) were selected because we wished to derive an analytic model for the outer planets variables. These are expressed in the form (D.2)
where stands for or , alternatively, ω_{o} = (g_{5}, g_{6}, g_{7}, g_{8}, s_{6}, s_{7}, s_{8}) is the vectorof the fundamental secular frequencies of the outer planet system, given in Table D.3, and is a septuple of integers.
In order to test this model, we compared it to the full solution LaX13b (Figs. D.1 and D.2). Only the real parts of and are represented because the imaginary parts lead to very similar plots. The full solution from LaX13b is plotted in purple, while the residuals after the analytical model QPSO was removed are plotted in green. For the inclination variables (), these residuals arevery small and appear as a straight line in the plots (Fig. D.2). This is not the case for the eccentricity variables (), where a significant band of residuals appears in green (Fig. D.1). This is mostly due to the shortperiod terms 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).
D.1 Secular initial conditions
To complete the QPSO solution or to integrate the secular equations by Gauss’s averaging method (Sec. 3), the values of the secular semimajor axes â_{k} are needed, or more precisely, those of the secular canonical variables , which are constant in the secular system. We defined the secular semimajor axes as the square of the secular averages in the LaX13b solution (Table D.1), i.e. . The mean mean motions N_{k} (secular averages of the mean motions) and the values of the mean longitudes at the origin λ_{0k} were also derived from the LaX13b solution and are provided in Table D.1. The remaining initial conditions of the secular system, corresponding to the complex variables , were obtained by a leastsquares fit of a polynomial of degree 5 in time to the full LaX13b solution over the first few thousand years after the shortperiod component was removed by Fourier filtering (Table D.2).
D.2 La2004 invariant reference frame
All the solutions given here are established in the invariant reference frame, which would ideally be the (x, y, z) reference frame, whose zaxis is aligned with the total angular momentum of the system, and with the xaxis pointing towards the equinox J2000. Nevertheless, this convention is not very practical because this invariant reference frame would change for any small variation of planetary masses or number of objects in the system. In the numerical longterm integrations of (Laskar & Gastineau 2009; Laskar et al. 2011), we therefore adopted a fixed reference frame, which is the conventional invariant reference frame of the orbital solution La2004 of (Laskar et al. 2004), which has been widely used in the paleoclimate community. We call this reference frame the La2004 invariant reference frame. It is thus a fixed reference frame that is derived from the ICRF equatorial reference frame (Ma et al. 1998) by the fixed transformation (D.3)
where R_{x} and R_{z} are rotation matrices defined as (D.4)
Secular average of the semimajor axes (⟨a_{k}⟩), mean mean motions (N_{k}), and mean longitudes at the origin (λ_{0k}) in the full solution LaX13b in the La2004 invariant reference frame (see Sec. D.2) derived by frequency analysis. The secular average of the actionlike variable is also provided because it cannot be directly obtained from ⟨a_{k}⟩: due to the contribution of the shortperiod terms, the average ⟨f(x)⟩ of f(x) is not f(⟨x⟩).
Initial conditions for the secular eccentricity and inclination variables at time J2000 in the La2004 invariant reference frame (see Sec. D.2), derived from the full LaX13b solution. In the first column,k is the index of the planet, and in columns 2 and 3, and are the real and imaginary part of the eccentricity variables . In columns 4and 5, and are the same quantities for the inclination variables.
Values of the fundamental secular frequencies of the outer Solar System derived by frequency analysis of the LaX13b solution over 30 Myr. These values were used in the QPSO model. They may slightly differ from the corresponding values given in (Laskar et al. 2004, 2011).
Quasiperiodic decomposition (Eq. D.2) of the secular solution for the outer planets for the eccentricity variable . 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 . 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}).
Quasiperiodic decomposition (Eq. D.2) of the secular solution for the outer planets for the eccentricity variable . 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 . 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}).
Fig. D.1 Real part of the eccentricity variable (Eq. D.1). The purple curve is the full LaX13b solution. The green curve is the residual after the contribution of the secular model QPSO is removed from LaX13b. The black curve is the same after the main shortperiod terms are removed from the residuals. 
Fig. D.2 Real part of the inclination variable (Eq. D.1). The purple curve is the full LaX13b solution. The green curve is the residual after the contribution of the secular model QPSO is removed from LaX13b. 
References
 Batygin, K., & Laughlin, G. 2008, ApJ, 683, 1207 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120 [CrossRef] [Google Scholar]
 Benettin, G., Galgani, L., & Strelcyn, J.M. 1976, Phys. Rev. A, 14, 2338 [Google Scholar]
 Boué, G., Laskar, J., & Farago, F. 2012, A&A, 548, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bour, E. 1855, PhD thesis, Faculté des Sciences de Paris [Google Scholar]
 Carpino, M., Milani, A., & Nobili, A. M. 1987, A&A, 181, 182 [NASA ADS] [Google Scholar]
 Chen, Y.C. 2017, ArXiv eprints [arXiv:1704.03924] [Google Scholar]
 Chirikov, B. V. 1979, Phys. Rep., 52, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, L. 1995, TimeFrequency Analysis (Upper Saddle River: NJ Prentice Hall PTR) [Google Scholar]
 Deprit, A. 1969, Celest. Mech., 1, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Efron, B. 1979, Ann. Stat., 7, 1 [Google Scholar]
 Fienga, A., Manche, H., Laskar, J., Gastineau, M., & Verma, A. 2014, arXiv eprints, [arXiv:1405.0484] [Google Scholar]
 Fouvry, J.B., Dehnen, W., Tremaine, S., & BarOr, B. 2020, AAS J., submitted, [arXiv:2011.01673] [Google Scholar]
 Fukushima, T. 2015, J. Comput. Appl. Math., 282, 71 [CrossRef] [Google Scholar]
 Gardiner, C. W. 1985, Stochastic Methods (Berlin, SpringerVerlag) [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra [Google Scholar]
 Gastineau, M., & Laskar, J. 2020, TRIP 1.4.94, TRIP Reference manual, IMCCE, Paris Observatory, https://www.imcce.fr/trip/ [Google Scholar]
 Gauss, C. 1818, Werke, 3, 331 [Google Scholar]
 Hill, G. W. 1882, Astron. Pap. Am. Ephemeris, 1, 315 [NASA ADS] [Google Scholar]
 Hoang, N. H., Mogavero, F., & Laskar, J. 2021, A&A, 654, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hori, G. 1966, PASJ, 18, 287 [NASA ADS] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1969, Mechanics (Pergamon Press) [Google Scholar]
 Laplace, P. S. 1785, Mém. Acad. Roy. Sci. Paris, 11, 95 [Google Scholar]
 Laskar, J. 1984, PhD thesis, Observatoire de Paris [Google Scholar]
 Laskar, J. 1985, A&A, 144, 133 [NASA ADS] [Google Scholar]
 Laskar, J. 1988, A&A, 198, 341 [NASA ADS] [Google Scholar]
 Laskar, J. 1989, Nature, 338, 237 [Google Scholar]
 Laskar, J. 1990a, in Modern Methods in Celestial Mechanics, eds. D. Benest, & C. Froeschle, 89 [Google Scholar]
 Laskar, J. 1990b, in Modern Methods in Celestial Mechanics, eds. D. Benest, & C. Froeschle (Editions Frontières, Gif SurYvette), 63 [Google Scholar]
 Laskar, J. 1990c, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1991, in NATO Advanced Study Institute (ASI) Series B, 272, Predictability, Stability, and Chaos in NBody Dynamical Systems, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1992, in Chaos, Resonance, and Collective Dynamical Phenomena in the Solar System, ed. S. FerrazMello, IAU Symposium, 152, 1 [Google Scholar]
 Laskar, J. 1993, Physica D, 67, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1994, A&A, 287, L9 [NASA ADS] [Google Scholar]
 Laskar, J. 1996, Celest. Mech. Dyn. Astron., 64, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1999, Phil. Trans. Roy. Soc. Lond. A, 357, 1735 [CrossRef] [Google Scholar]
 Laskar, J. 2005, in Hamiltonian Systems and Fourier Analysis: New Prospects for Gravitational Dynamics, eds. D. Benest, C. Froeschlé, & E. Lega (Cambridge Scientific Publishers Ltd), 93 [Google Scholar]
 Laskar, J. 2008, Icarus, 196, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [Google Scholar]
 Laskar, J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Simon, J. L. 1988, Celest. Mech., 43, 37 [NASA ADS] [Google Scholar]
 Laskar, J., Froeschlé, C., & Celletti, A. 1992, Physica D, 56, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011, A&A, 532, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lithwick, Y., & Wu, Y. 2011, ApJ, 739, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, C., Arias, E. F., Eubanks, T. M., et al. 1998, AJ, 116, 516 [Google Scholar]
 Mei, L., & Huang, L. 2018, Comput. Phys. Commun., 224, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Mogavero, F. 2017, A&A, 606, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (Taylor & Francis) [Google Scholar]
 Musen, P. 1970, Celest. Mech., 2, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Nobili, A. M., Milani, A., & Carpino, M. 1989, A&A, 210, 313 [NASA ADS] [Google Scholar]
 Olver, F. W. J., Olde Daalhuis, A. B., Lozier, D. W., et al. 2020, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.28 of 20200915 [Google Scholar]
 Oseledec, V. I. 1968, Trans. Moscow Math. Soc., 19, 197 [Google Scholar]
 Parzen, E. 1962, Ann. Math. Statist., 33, 1065 [CrossRef] [Google Scholar]
 Poincaré, H. 1896, Compt. Rend. Hebdomad. Acad. Sci. Paris, 123, 1031 [Google Scholar]
 Rein, H., & Tamayo, D. 2018, MNRAS, 473, 3351 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenblatt, M. 1956, Ann. Math. Stat., 27, 832 [Google Scholar]
 Rybicki, K. R., & Denis, C. 2001, Icarus, 151, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Sackmann, I. J., Boothroyd, A. I., & Kraemer, K. E. 1993, ApJ, 418, 457 [CrossRef] [Google Scholar]
 Saha, P., & Tremaine, S. 1992, AJ, 104, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Schröder, K. P., & Smith, R. C. 2008, MNRAS, 386, 155 [CrossRef] [Google Scholar]
 Schwarz, W. 1992, J. Appl. Prob., 29, 597 [Google Scholar]
 Silverman, B. W. 1986, Density Estimation for Statistics and Data Analysis (London: Chapman & Hall/CRC) [Google Scholar]
 Sussman, G. J., & Wisdom, J. 1992, Science, 257, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Tancredi, G., Sánchez, A., & Roig, F. 2001, AJ, 121, 1171 [Google Scholar]
 Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 [NASA ADS] [CrossRef] [Google Scholar]
 Verma, A. K., Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2014, A&A, 561, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wilson, E. B. 1927, J. Am. Stat. Assoc., 22, 209 [CrossRef] [Google Scholar]
 Woillez, E., & Bouchet, F. 2020, Phys. Rev. Lett., 125, 021101 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, W., & Zurbenko, I. 2010, WIREs Comput. Stat., 2, 340 [CrossRef] [Google Scholar]
Even the secondorder secular dynamics in Laskar (1985) needed an ad hoc correction of 0.27′′ yr^{−1} for the Saturndominated eccentricity mode g_{6} (Laskar 1988).
The symmetries of the planetary Hamiltonian (1) imply constraints on the exponents known as D’Alembert rules (e.g. Laskar & Robutel 1995; Morbidelli 2002). The rotational invariance requires , while it follows from planar symmetry that the terms of the series are even with respect to the variables .
The terms in only involving outer planet variables can be discarded in Eq. (16) because they depend on time alone and no longer affect the inner planet dynamics.
All the first secular collisions, with the exception of a single EarthMars event, involved the MercuryVenus pair (see Sect. 7.1)
All Tables
Frequency in arcsec yr^{−1} of the largestamplitude Fourier harmonic for each proper mode of the Poincaré variables over the timeinterval [0, 20] Myr.
Confidence interval at the 98% level of the probability P(sup_{t≤5 Gyr}e_{1}(t) > e_{max}) in percent, where e_{1} is the eccentricity of Mercury.
Frequency response of the KZ filter for M ≫ 1. Maximum gain in the stopband α_{k} = η^{k}. Ratio of the cutoff frequency f_{c} to the passband frequency f_{p} for two different maximum loss ρ.
Number of solutions that reach a given eccentricity of Mercury (e_{m0}) over a given time (500, 1000, 1500, 2000, 3000, 4000, and 5000 Myr). The statistics are made over 2501 orbital solutions as publishedin (Laskar & Gastineau 2009). In the three 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.
Secular average of the semimajor axes (⟨a_{k}⟩), mean mean motions (N_{k}), and mean longitudes at the origin (λ_{0k}) in the full solution LaX13b in the La2004 invariant reference frame (see Sec. D.2) derived by frequency analysis. The secular average of the actionlike variable is also provided because it cannot be directly obtained from ⟨a_{k}⟩: due to the contribution of the shortperiod terms, the average ⟨f(x)⟩ of f(x) is not f(⟨x⟩).
Initial conditions for the secular eccentricity and inclination variables at time J2000 in the La2004 invariant reference frame (see Sec. D.2), derived from the full LaX13b solution. In the first column,k is the index of the planet, and in columns 2 and 3, and are the real and imaginary part of the eccentricity variables . In columns 4and 5, and are the same quantities for the inclination variables.
Values of the fundamental secular frequencies of the outer Solar System derived by frequency analysis of the LaX13b solution over 30 Myr. These values were used in the QPSO model. They may slightly differ from the corresponding values given in (Laskar et al. 2004, 2011).
Quasiperiodic decomposition (Eq. D.2) of the secular solution for the outer planets for the eccentricity variable . 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 . 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}).
Quasiperiodic decomposition (Eq. D.2) of the secular solution for the outer planets for the eccentricity variable . 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 . 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}).
All Figures
Fig. 1 Eccentricities (left side) and inclinations (right side) over the invariant plane J2000 (see Appendix D.2) of the inner planets over 250 000 yr in the future. The solid black line stands for the full (i.e. nonfiltered) direct integration LaX13b, while the red line is the solution of our forced model in Eq. (16). The secular integration La90 is represented by the dashed blue line. The maximum relative root mean square of our model residuals is 2% for the eccentricities and 0.7% for theinclinations. 

In the text 
Fig. 2 Eccentricities (left side) and inclinations (right side) over the invariant plane J2000 (see Appendix D.2) of the inner planets over10 Myr in the future. The solid black line stands for the full (i.e. nonfiltered) direct integration LaX13b, and the red line is the solution of our forced model in Eq. (16). The secular solution La90 nearly coincides with the direct integration in these plots and it is therefore not shown. 

In the text 
Fig. 3 Finitetime maximum Lyapunov exponent (FTMLE) and corresponding Lyapunov time FTMLE^{−1} of the forced inner Solar System over 5 Gyr from an ensemble of 1148 stable orbital solutions with very close initial conditions. 

In the text 
Fig. 4 Probability density function of the Lyapunov time FTMLE^{−1} at 5 Gyr. The kernel density estimate is shown as a dark blue line and the corresponding CDF is shown on the upper horizontal axis. The pointwise confidence interval from bootstrap is at the 98% level and is shown as a light blue region. 

In the text 
Fig. 5 Evolution of the fundamental frequency g_{1}, as defined by Eq. (56), in the 42 collisional orbital solutions (top panel). Six solutions are isolated in the bottom panel. The horizontal line stands for the constant fundamental frequency g_{5}. 

In the text 
Fig. 6 Kernel density estimation of the PDF of the time τ_{c} of the first secular collision from 1042 integrations of the forced inner system over 100 Gyr (blue line in the lower panel). The estimate of the CDF is shown as a blue line in the upper panel. The blue regions represent the pointwise confidence intervals at the 98% level from bootstrap. In red we show the PDF estimate over 5 Gyr from the 10 560 orbital solutions of Sect. 7. The dashed black lines stand for the analytical PDF (57) with T_{1} = 27.6 Gyr and α = 0.9 and the corresponding CDF. 

In the text 
Fig. B.1 Impulse response of the KolmogorovZurbenko filter (the moving average corresponds to the case k = 1). The axes have been normalized to show the asymptotic Gaussian behaviour of the function. The standard normal distribution N(0, 1) is represented by the solid black curve. 

In the text 
Fig. B.2 Absolute value of the frequency response for the KolmogorovZurbenko filter (the moving average corresponds to the case k = 1). The vertical dashed black line stands for the cutoff frequency f_{c} given by Eq. (B.14), while the horizontal coloured lines represent the maximum gains in the stopband α_{k} = η^{k}. The vertical dotted coloured lines show the passband frequencies f_{p} given by Eq. (B.17) for a loss level of 5%, which is represented by the horizontal dotted black line. 

In the text 
Fig. D.1 Real part of the eccentricity variable (Eq. D.1). The purple curve is the full LaX13b solution. The green curve is the residual after the contribution of the secular model QPSO is removed from LaX13b. The black curve is the same after the main shortperiod terms are removed from the residuals. 

In the text 
Fig. D.2 Real part of the inclination variable (Eq. D.1). The purple curve is the full LaX13b solution. The green curve is the residual after the contribution of the secular model QPSO is removed from LaX13b. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.