© F. Mogavero and J. Laskar 2021

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₃ − g₄) − (s₃ − s₄) and (g₃ − g₄) − (s₃ − s₄), 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⁻¹ scale in a few dozen million years (Myr), with the exception of the Venus-dominated eccentricity mode g₂, which has somewhat smaller variations (Laskar 1990c; Laskar et al. 2004). This behaviour reveals the essential high-dimensional structure of the inner Solar System, which limits the possibility of faithfully modelling its dynamics with 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 long-term evolution of their orbits. On the one hand, such a model must be consistent with the predictions of the reference numerical integrations available in literature (Laskar 1990c, 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 high-dimensional 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 N-body Hamiltonian of the Newtonian gravitational interactions among the Solar System planets, with the addition of the leading corrections coming from general relativity and the Earth–Moon interaction. This already reproduces the precession frequencies of the inner orbits with a precision better than 0.01′′ yr⁻¹ (Laskar 1999). However, this Hamiltonian is unnecessarily complicated because it includes short-time harmonics, with periods shorter than 5000 yr (e.g. Carpino et al. 1987), which are known to generate small quasi-periodic oscillations in the inner orbits, without being implied in chaos generation. The inner planets are not, indeed, involved in any relevant mean-motion resonance. At the same time, the long-term numerical integration of the corresponding equations of motion is very time consuming, the resulting solutions needing to be filtered to extract the secular trend of the orbits (Carpino et al. 1987; Nobili et al. 1989). These facts suggest considering a secular Hamiltonian to directly describe the slow movement of the planet perihelia and nodes after proper averaging over the short-time orbital motion (Laskar 1984, 1985). Secular dynamics includes the essential planet interactions responsible for chaos in the inner Solar System and allows performing the fastest long-term 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 mean-motion resonance between Jupiter and Saturn, the so-called great inequality (Laplace 1785; Laskar 1996). A simple averaging of the N-body Hamiltonian over the planet mean longitudes, resulting in a first-order secular dynamics in planet masses, would reproduce the fundamental frequencies g₅ and g₆ (which dominate the perihelion precession of Jupiter and Saturn, respectively) very poorly¹ (Laskar 1988). The construction of higher-order models requires the manipulation of large Poisson series and the use of sophisticated computer algebra systems. This is probably the reason why they are still not widely used, at least as a basis of extensive research.

This paper introduces a new secular dynamics for the orbits of the inner planets. It is based on the practical long-term regularity 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 mean-motion resonances in the inner system (Sect. 2). We show that our model can be numerically integrated with the so-called 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 N-body 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₁ − g₅ 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₀, , indexed in order of increasing semi-major 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₀ are the planet heliocentric coordinates and their conjugated momenta, μ_k = m₀m_k∕(m₀ + 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 momentum-coordinate 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 semi-major axes, λ_k the mean longitudes, e_k the eccentricities, i_k the inclinations, ϖ_k the longitudes of the perihelia, and Ω_k the longitudes of the nodes (Laskar 1991; Laskar & Robutel 1995). Throughout the paper, 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 two-body 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 non-negative integers². Following Laskar & Robutel (1995), we defined the dimensionless Poincaré variables , and the semi-major axis ratio α = a∕a′, with a < a′. The analytical expression of the coefficients , only depending on the semi-major axis ratio, is given in Laskar & Robutel (1995) in terms of Laplace coefficients. The indirect part of the two-body 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, λ = (λ₁, …, λ₈), and the coefficients depend on all the remaining canonical variables.

2.1 Secular dynamics

The long-term 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, orbit-averaged, gravitational interactions (Laskar 1990c; Laskar et al. 2004). A secular Hamiltonian describing this long-term motion can be introduced in its simplest form by searching for a change of variables that eliminate at first order in the planet masses all the harmonics with non-null 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 time-1 flow of a generating Hamiltonian S satisfying the homologic equation (7)

where the braces represent the Poisson bracket. The angle-bracket 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₁ ⟩ 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₀∕∂Λ 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 short-time oscillations generated by the n ≥ 1 terms of the Lie transforms. Because we are interested in the long-term 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 mean-motion resonances. More precisely, a maximum contribution of only 0.07′′ yr⁻¹ to the fundamental precession frequencies of the inner orbits in the Laplace-Lagrange solution arises from the second order in the planet masses (Laskar 1985, Table 8). This same contribution is 0.9′′ yr⁻¹ 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 higher-order 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 mean-motion resonances, in Eq. (9) the denominators ℓ ⋅ n involving at least one inner planet are sufficiently far from zero. Under these assumptions, the higher-order terms only produce small corrections to the dynamics generated by the leading terms. Using the homologic Eq. (7), we thus obtain (12)

The resulting secular Hamiltonian is simply the average of the N-body Hamiltonian (1) over the planet mean longitudes. The averaging process is mathematically equivalent to replacing each planet by its instantaneous Keplerian orbit, with the corresponding mass distributed along it in a way that is inversely proportional to the local orbital speed of the planet. The secular dynamics is thus the slow gravitational interaction of 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 semi-major axes of the planets.The resulting Hamiltonian system consists of two degrees of freedom for each planet³, 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⁻¹ to the secular precession of the Mercury perihelion. This correction is critical for the statistics of the long-term destabilization of the inner orbits (Laskar 2008; Laskar & Gastineau 2009) because it moves the system away from the g₁ − g₅ 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 semi-major axis a_k and is thus a constant quantity in the secular dynamics. The next largest effect to be taken into account would be the Earth-Moon gravitational interaction, accounting for a 0.077′′ yr⁻¹ 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⁻¹ and 0.062′′ yr⁻¹ to the Earth-dominated and Mars-dominated eccentricity modes g₃ and g₄, 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 Earth-Moon interaction in our model, preferring to work with a Hamiltonian whose dependence on eccentricities and inclinations is exact, as compared to a non-decisive increase in precision on the orbit precession frequencies. Moreover, we know that the N-body Hamiltonian (1), corrected for general relativity, reproduces the maximum Lyapunov exponent of the inner Solar System (Rein & Tamayo 2018).

Therefore, apart from irrelevant constant terms only depending on the semi-major axes in Eqs. (3) and (13), the secular 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 Saturn-dominated eccentricity modes g₅ and g₆, 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 quasi-periodic 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) = ω_ot, with ω_o the vector of the (constant) fundamental precession frequencies of the outer orbits, denoted as ω_o = (g₅, g₆, g₇, g₈, s₆, s₇, s₈)⁴ (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⁵. 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 second-order 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 long-term 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 quasi-periodic 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 high-precision 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 quasi-periodic decomposition from Eq. (15) and filtering out of the non-secular Fourier components. The effectiveness of the quasi-periodic approximation is illustrated in Figs. D.1 and D.2.

The choice of the secular semi-major axes and that of the initial conditions for the secular Poincaré 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 short-time (high-frequency) 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 short-time components of the numerical solution of the non-averaged 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 semi-major 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 non-algebraic, and the numerical computation of the derivatives is in principle much more complex than in N-body 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 so-called 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 arithmetic-geometric 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 floating-point 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⁻¹². 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 δ₀ = 9 × 10⁻¹⁵, T₀ = 0.6 Myr and α = 1.2. Therefore we have δ(10 Myr) ≈ 3 × 10⁻¹³, and δ(10 Gyr) ≈ 10⁻⁹ 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 fixed-time-step 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 semi-major 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 plain-text 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⁶. 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 Laplace-Lagrange 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 action-angle variables.

Solution of the equations of motion

The matrix M_ii is real and symmetric and can thus be diagonalized through an orthogonal matrix O_M, (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⁷, g_LL being a column vector. The orthogonal matrix O_M induces the canonical change of variables⁸ defined as (30)

The transformed Hamiltonian reads (31)

with . The corresponding Hamilton equations are given by (32)

and constitute a first-order 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 time-depending 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)

which imply (41)

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 Laplace-Lagrange dynamics, they simply rotate in the complex plane, with constant angular frequencies given by g_LL. According to Eq. (42), the corresponding Poincaré variables thus result in a superposition of these independent harmonic oscillations and those arising from the forcing of the outer planets. When higher-degree terms of the Hamiltonian (25) are taken into account, the dynamics of the proper modes becomes coupled, and their frequency spectrum is no 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 non-linear in the forcing of the outer planets, in the sense that harmonics of order⁹ |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 time-dependent canonical change of variables ,

We have defined (46)

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 quasi-periodic 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 higher-degree truncations, we first considered the non-quadratic 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 quasi-periodic decompositions given in Eq. (15). However, these substitutions do not conserve the degree of the terms in the expansion because non-linear 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 higher-degree terms in the substitution process. To define the terms in a consistent way, we therefore introduced a fictitious real variable ϵ to redefine the quasi-periodic 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¹⁰. 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 quasi-periodic decomposition (15) we used in this work (see Appendix D).

4.3 Action-angle variables

Action-angle 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) = ω_ot, (52)

where are complex amplitudes, and we employed a compact notation for the action-angle variables, (53)

Only a finite number of harmonics have non-zero 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 phase-space extension through the definition of action-angle 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.

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 short-time (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 long-term 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⁻¹ in the case of s₃ and s₄. 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₁ − g₅) − (s₁ − s₂) and 2(g₃ − g₄) − (s₃ − s₄) described in Laskar (1990c), the deviations are only 0.01′′ yr⁻¹ and 0.004′′ yr⁻¹, 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₁ − g₅) − (s₁ − s₂) and 2(g₃ − g₄) − (s₃ − s₄) are considered, models and almost show the same deviations from solution LaX13b as the non-truncated 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 non-truncated 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. non-filtered) 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. non-filtered) 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 Finite-time 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 non-null 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 infinite-time limit, as in its usual mathematical definition (Oseledec 1968); the consideration of a finite-time MLE (FT-MLE) 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 FT-MLE given by the two-particle 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₀ 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₀|| from the reference particle. The resulting FT-MLE 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₀, 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 FT-MLE 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 FT-MLE expressed as an angular frequency in arcsec yr⁻¹, while on the right axis, we show the corresponding Lyapunov time, defined as FT-MLE⁻¹, given in million years. The dark grey region represents at each renormalization time the interquantile range of the observed probability density function (PDF) of the FT-MLE, which 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 FT-MLEs reported in Laskar (1989) are also shown with black lines.

As stated above, the initial conditions of each solution in the ensemble are very close to the nominal conditions. They were chosen by taking the relative variation of each coordinate of the nominal phase-space vector as a normal random variable with zero mean and a standard deviation of 10⁻⁹. The initial conditions thus follow a multivariate Gaussian distribution centred at the nominal phase-space position of the system, with a diagonal covariance matrix. The width of 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⁻⁸. Such a value shouldminimize the accumulation over long timescales of numerical errors, due to the floating-point 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 FT-MLE for the nominal initial conditions is shown by dotted blue lines for a set of 16 different initial tangent vectors d₀. This manifestly shows that our FT-MLEs are asymptotically independent of d₀, 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 phase-space vectors, taken here to be Euclidean as usual.

For times shorter than 100 Myr, Fig. 3 shows that the distribution of the FT-MLE only reflects the choice of different initial tangent vectors d₀ around essentially the same reference trajectory, that is, the nominal trajectory. Therefore this part of the plot can be discarded because of its non-asymptotic character. We also note that for t < 100 Myr, the PDF of the FT-MLE tends to shrinks with increasing time, as expected by the asymptotic behaviour of the method (see the blue lines and two curves from Laskar 1989). For longer times, the PDF begins to broaden because each orbital solution acquires a macroscopically different FT-MLE, which chaotically depends on its initial conditions. The time-asymptotic regions shown in Fig. 3 thus represent the probabilistic knowledge of the FT-MLE of the inner Solar System, which arises from its chaotic behaviour and our determination of the current planet positions and velocities.

The dark blue line in Fig. 4 shows the kernel density estimate (Rosenblatt 1956; Parzen 1962) of the PDF of the Lyapunov time FT-MLE⁻¹ 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 FT-MLE 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 Finite-time maximum Lyapunov exponent (FT-MLE) and corresponding Lyapunov time FT-MLE−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 FT-MLE−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 long-term approximation of the original N-body dynamics, in which the involved planets can experience a close encounter, possibly leading to a physical collision or to the escape of one of two planets from the Solar System (Laskar 1994). We thus denote by τ_c the time of the first secular collision for a given orbital solution. The statistics of the first-collision 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 fixed-time-step 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 phase-space coordinate is 10⁻⁹, 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¹¹, the probabilities of the high Mercury eccentricities in the forced inner system are fully compatible with those arising in the secular model of Laskar (1990c). They also represent a very good estimate of the probabilities resulting from N-body dynamics, even if 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 quasi-periodic time dependence. Moreover, the Hamiltonian terms at second order in planet masses could also contribute to somewhat larger chaotic excursions of the inner orbits.

Fig. 5 Evolution of the fundamental frequency g1, 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 g5.

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 non-averaged inner Solar System over the next 5 Gyr, probably somewhat lower, in agreement with the previous discussion.

All the secular collisions in our ensemble of orbital solutions involve the Mercury–Venus 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 non-averaged 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₁ – g₅

The destabilizing role of the secular resonance g₁ − g₅, 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₁ − g₅ resonance isindeed reached by the collisional orbital solutions of our ensemble. While g₅ ≈ 4.257′′ yr⁻¹ is constant in our dynamical model, the frequency g₁, which dominates the spectrum of the Poincaré variable x₁ 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₁(t_n) is thus computed through Eq. (42) and the corresponding angle variable χ₁ (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 ω₁ of the proper mode u₁ is given by (Cohen 1995) (55)

We could analytically compute the derivative of Eq. (42) to obtain the time series ω₁ (t_n). In practice, we employed numerical differentiation via the Lagrange five-point formula (Olver et al. 2020) because very high numerical precision is unnecessary here. As it is, the frequency ω₁ (t_n) presents short-time oscillations that hide its long-term evolution. To suppress these high-frequency fluctuations, we therefore defined the time series g₁(t_n) as the output of the low-pass Kolmogorov–Zurbenko (KZ) filter (Yang & Zurbenko 2010, and references therein) applied to ω₁ (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⁻¹ (see Appendix B). This particular choice is motivated by the typical duration of the libration episodes of the g₁ − g₅ 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₁ 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₅ and thus the location of the g₁ − g₅ resonance. Within the first 2 Gyr, the chaotic diffusion of g₁ in the solutions of Fig. 5 is limited to no more than 0.25′′ yr⁻¹ with respect to its initial value. Nevertheless, when lower values of about 5′′ yr⁻¹ are eventually reached, the dynamics of g₁ starts to be characterized by large random variations. This seems to reproduce the large chaotic zone related to the resonance g₁ − g₅, with a half-width of more than 1′′ yr⁻¹, which has been reported in Laskar (2008, Fig. 6(a2))¹². 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 semi-major axes of the planet will adiabatically inflate by a factor of 1.85. Even more critically for the fate of the inner system, it is likely that the innermost planets Mercury, Venus, and the Earth will be engulfed by the Sun as its radius expands along the red giant branch, marking the end of their existence (Rybicki & Denis 2001; Schröder & Smith 2008). Therefore this section intends to investigate the very long-term stochastic behaviour of the dynamical system (16), without inferring any conclusion on the real evolution of the inner system over 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¹³, 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 T1 = 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 white-noise limit for slow-fast dynamical systems (e.g. Gardiner 1985) to describe the long-term 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 one-dimensional 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 long-term wanderings of the fundamental frequency g₁ towards the destabilizing secular resonance g₁ − g₅ (see Sect. 7.2) in an approximate way. Numerically, the PDF of g₁ is asymmetric, with a longer tail that diffuses towards low values, while its upper end is nearly fixed at about 6′′ yr⁻¹ (Hoang et al. 2021, see also Fig. 5). Therefore we assumed that g₁ performs a random walk described by a Wiener process starting at g_1,0 = g₁(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₁ − g₅ is reached, as shown in Fig. 5; this was taken into account by considering an absorbing barrier at g₁ = g₅. The PDF ρ(τ) of the hitting time is thus given by (e.g. Schwarz 1992; Woillez & Bouchet 2020) (57)

where α = (g_1,0 − g₅)∕(g_1,max − g₅) and , with D being the diffusion coefficient ofthe Wiener process. After a least-squares search, we plot in Fig. 6 the curve (57) with parameters T₁ = 27.6 Gyr and α =0.9. Even if the real random walk performed by the frequency g₁ 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¹⁴. On the one hand, this shows that at short times, the PDF behaves as (58)

where T₀ = α²T₁. 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₁. 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₀ = 1.56 Gyr, a value that is more than ten times lower than that in Fig. 6, that is, T₀ = 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¹⁵. 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₁ from diffusing at a rate determined by the leading secular resonances. In particular, the resonance (g₁ − g₅) − (s₁ − s₂), with a libration frequency of 0.12′′ yr⁻¹, is among the leading resonances affecting the frequency g₁ (Laskar 1990c). We may estimate an upper bound on the diffusion coefficient D as (59)

with 2π MLE = 0.1′′ yr⁻¹ (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⁻¹ (Table 3) and g₅ = 4.257′′ yr⁻¹. 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₁ − g₅ 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 mean-motion resonances in the inner system. The range of validity of the resulting dynamics extends to very high eccentricities and inclinations, up to a secular planetary collision, 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₁ − g₅ 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₁. 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₁ up to the g₁ − g₅ 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 two-body secular Hamiltonian reported in Laskar & Robutel (1995, Appendix) and the quasi-periodic 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 post-doctoral fellowship and by a grant of the French Agence Nationale de la Recherche (AstroMeso ANR-19-CE31-0002-01). This project has been supported from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo-885250). This work was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

Appendix A Indirect part of the two-body Hamiltonian perturbing function

We outline the derivation of the coefficients of the indirect part of the two-body 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 x-axis directed towards the pericenter, and the z-axis 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 Kolmogorov-Zurbenko filter

Low-pass filters are employed in time-series analysis to extract the long-term (low-frequency) components of a given signal. A dedicated low-pass filter has been constructed in (Carpino et al. 1987) to recover the secular changes of the giant planet orbits from the numerical output of an N-body integration. Here we propose the use of the Kolmogorov-Zurbenko filter (Yang & Zurbenko 2010) as an out-of-the-box, computationally advantageous, and still effective choice.

We consider the real-valued 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 discrete-time 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 finite-impulse-response 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 non-negative integer. The sequence d_n is the impulse response function of the filter, and L_M = 2M + 1 is the filter length. A liner filter is characterized by the discrete-time Fourier 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 low-pass 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 Kolmogorov-Zurbenko 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 low-pass 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 Kolmogorov-Zurbenko 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 Kolmogorov-Zurbenko filter (the moving average corresponds to the case k = 1). The vertical dashed black line stands for the cutoff frequency fc 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 fp 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 low-pass 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.

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 high-order 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.

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 Quasi-periodic secular solution for the outer planets

We use a quasi-periodic solution for the secular evolution of the orbits of the outer planets (Jupiter to Uranus), denoted hereafter as QPSO. In order to achieve a realistic modelling, this solution 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 high-precision 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 quasi-periodic approximation through frequency analysis (Laskar 1988, 2005). The variables considered here are the Poincaré complex canonical variables in Eq. (2), scaled to suppress the semi-major axis (or Λ) dependence,that is, (D.1)

It should be noted that and where z = eE^jϖ and ζ = sin(i∕2)E^jΩ are the classical, non-canonical, complex elliptic elements.

The derivation of the QPSO model needs some care. 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 mean-motion resonance among Jupiter and Saturn and the proximity of the 2:1 mean motion resonance in the Uranus-Neptune system, several short-period terms (with arguments involving the planetary mean motions) have also to be discarded. The solution comprises a small number of terms and is fully given in Tables D.5 and D.6. 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₅, g₆, g₇, g₈, s₆, s₇, s₈) 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 short-period 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 semi-major axes â_k are needed, or more precisely, those of the secular canonical variables , which are constant in the secular system. We defined the secular semi-major 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 least-squares fit of a polynomial of degree 5 in time to the full LaX13b solution over the first few thousand years after the short-period 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 z-axis is aligned with the total angular momentum of the system, and with the x-axis 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 long-term 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)

and (D.5)

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 short-period 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

All Tables

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. non-filtered) 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. non-filtered) 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 Finite-time maximum Lyapunov exponent (FT-MLE) and corresponding Lyapunov time FT-MLE−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 FT-MLE−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 g1, 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 g5.
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 T1 = 27.6 Gyr and α = 0.9 and the corresponding CDF.

In the text

	Fig. B.1 Impulse response of the Kolmogorov-Zurbenko 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 Kolmogorov-Zurbenko filter (the moving average corresponds to the case k = 1). The vertical dashed black line stands for the cutoff frequency fc 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 fp 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 short-period 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