11institutetext: Astronomie et Systèmes Dynamiques, Institut de Mécanique
Céleste et de Calcul des Éphémérides
CNRS UMR 8028, Observatoire de Paris, Université PSL, Sorbonne Université, 77
Avenue Denfert-Rochereau, 75014 Paris, France
11email<EMAIL_ADDRESS>
# Long-term dynamics of the solar system inner planets

Federico Mogavero  Jacques Laskar

(Received ; accepted )

Although the discovery of the chaotic motion of the inner planets in the solar
system (Mercury to Mars) dates back to more than thirty years ago, the secular
chaos of their orbits still dares more analytical analyses. Apart from the
high-dimensional structure of the motion, this is probably related to the lack
of an adequately simple dynamical model. In this paper, we consider a new
secular dynamics for the inner planets, with the aim of retaining a
fundamental set of interactions responsible for their chaotic behaviour, while
being consistent with the predictions of the most precise orbital solutions
currently available. We exploit the regularity in the secular motion of the
outer planets (Jupiter to Neptune), to predetermine a quasi-periodic solution
for their orbits. This reduces the secular phase space to the degrees of
freedom dominated by the inner planets. On top of that, the smallness of the
inner planet masses and the absence of strong mean-motion resonances permits
to restrict ourselves to first-order secular averaging. The resulting dynamics
can be integrated numerically in a very efficient way through Gauss’s method,
while computer algebra allows for analytical inspection of planet
interactions, once the Hamiltonian is truncated at a given total degree in
eccentricities and inclinations. The new model matches very satisfactorily
reference orbital solutions of the solar system over timescales shorter than
or comparable to the Lyapunov time. It correctly reproduces the maximum
Lyapunov exponent of the inner system and the statistics of the high
eccentricities of Mercury over the next five billion years. The destabilizing
role of the $g_{1}-g_{5}$ secular resonance also arises. A numerical
experiment, consisting of a thousand orbital solutions over one hundred
billion years, reveals the essential properties of the stochastic process
driving the destabilization of the inner solar system and clarifies its
current metastable state.

###### Key Words.:

Solar system – Chaos – Celestial mechanics – Planets and satellites: dynamical
evolution and stability – Methods: analytical – Methods: numerical

##  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 if its chaotic behaviour
has been demonstrated numerically three decades ago, still no analytical study
corroborates or rules out the role in chaos generation of the secular
resonances $2(g_{3}-g_{4})-(s_{3}-s_{4})$ and $(g_{3}-g_{4})-(s_{3}-s_{4})$,
between the Earth and Mars fundamental precession frequencies, as proposed in
(Laskar, 1990c, 1992) and supported by (Laskar et al., 2004). Such a fact
probably relates to the complex network of dynamical interactions among the
inner planets: the fundamental precession frequencies of their orbits vary in
an intricate way over a 0.1$\arcsec\,{\textrm{yr}}^{-1}$ scale in a few tens
of million years (Myr), with the exception of the Venus-dominated eccentricity
mode $g_{2}$ which has somewhat smaller variations (Laskar, 1990c; Laskar et
al., 2004). This behaviour reveals the essential high-dimensional structure of
the inner solar system, which limits the possibility of faithfully modelling
its dynamics with few effective degrees of freedom (e.g., Lithwick & Wu, 2011;
Batygin et al., 2015).

Getting analytical insight into the motion of the inner planets requires an
appropriate dynamical modelling of the long-term evolution of their orbits. On
the one hand, such a model must be consistent with the predictions of the
reference numerical integrations available in literature (Laskar, 1990c;
Laskar et al., 2004; Laskar, 2008; Laskar & Gastineau, 2009), to ensure it
reproduces the dynamical features of the inner system with sufficient
precision. On the other hand, the corresponding Hamiltonian should be put in a
form suitable for the systematic application of canonical perturbation
techniques (Hori, 1966; Deprit, 1969), which is essential to an unbiased
analysis of such a high-dimensional dynamics. Moreover, the possibility to
numerically integrate the equations of motion in an efficient way is
fundamental to study the chaotic evolution of the orbital solutions in a
statistical way. Unfortunately, the construction of such a model turns out to
be a delicate task. In principle, one could just consider the full N-body
Hamiltonian of the Newtonian gravitational interactions among the solar system
planets, with the addition of the leading corrections coming from general
relativity and the Earth-Moon interaction. This already reproduces the
precession frequencies of the inner orbits with a precision better than
0.01$\arcsec\,{\textrm{yr}}^{-1}$ (Laskar, 1999). However, such a Hamiltonian
is unnecessarily complicated, as it includes short-time harmonics, with
periods of less than 5 000 yr (e.g., Carpino et al., 1987), which are known to
generate small quasi-periodic oscillations in the inner orbits, without being
implied in chaos generation. Indeed, the inner planets are not involved in any
relevant 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 to
consider a secular Hamiltonian to directly describe the slow movement of the
planet perihelia and nodes, after proper averaging over the short-time orbital
motion (Laskar, 1984, 1985). Secular dynamics includes the essential planet
interactions responsible for chaos in the inner solar system and allows to
perform the fastest long-term numerical integrations (Laskar, 1988, 1989,
1994, 2008). Indeed, this made the discovery of chaos in the inner system
possible before the use of symplectic integration schemes (Sussman & Wisdom,
1992). Unfortunately, an effective secular model for the entire solar system
has to be of high order in planet masses, principally because of the 5:2 near
mean-motion resonance between Jupiter and Saturn, the so-called Great
Inequality (Laplace, 1785; Laskar, 1996). A simple averaging of the N-body
Hamiltonian over the planet mean longitudes, resulting in a first-order
secular dynamics in planet masses, would reproduce very poorly the fundamental
frequencies $g_{5}$ and $g_{6}$, which dominate the perihelion precession of
Jupiter and Saturn, respectively (Laskar, 1988)111Even the second-order
secular dynamics in (Laskar, 1985) needed an ad-hoc correction of
0.27$\arcsec\,{\textrm{yr}}^{-1}$ for the Saturn-dominated eccentricity mode
$g_{6}$ (Laskar, 1988).. The construction of such higher-order models requires
the manipulation of large Poisson series and the use of sophisticated computer
algebra systems. This is probably the reason why they are still not widely
used, at least as a basis of extensive research.

This paper introduces a new secular dynamics for the orbits of the inner
planets. It is based on the practical long-term regularity of the outer planet
trajectories (Laskar, 1990c; Laskar et al., 2004; Hoang et al., 2021), the
smallness of the inner planet masses and the absence of relevant mean-motion
resonances in the inner system (Sect. 2). We show that the present model can
be numerically integrated via the so-called Gauss (1818)’s method (Sect. 3),
while its Hamiltonian is suitable for a systematic expansion in planet
eccentricities and inclinations, once one employs the algorithm of (Laskar &
Robutel, 1995) and a computer algebra system like TRIP (Gastineau & Laskar,
2011, 2020) (Sect. 4). We compare its orbital solution to a reference N-body
integration over short (secular) timescales (Sect. 5) and compute the
corresponding maximal Lyapunov exponent in a statistical fashion (Sect. 6). We
then determine the percentages of the high Mercury eccentricities over the
next 5 billion years and highlight the destabilizing role of the $g_{1}-g_{5}$
secular resonance (Sect. 7). Finally, we perform a new numerical simulation
involving a thousand orbital solutions over one hundred billion years, to
characterize the effective stochastic process which drives the destabilization
of the inner solar system, and discuss its implications on the conjecture of
marginal stability of a secularly evolving planetary system formulated in
(Laskar, 1996) (Sect. 8).

##  2 Dynamical model

We model the dynamics of the largest bodies in the solar system, by
considering the Sun and the $N=8$ planets as point masses $m_{0}$,
$(m_{k})_{k=1,N}$, indexed in order of increasing semi-major axis. The
barycentric coordinates of the bodies and the corresponding momenta are
denoted by $(\bm{R}_{0},\tilde{\bm{R}}_{0}=m_{0}\dot{\bm{R}}_{0})$,
$(\bm{R}_{k},\tilde{\bm{R}}_{k}=m_{k}\dot{\bm{R}}_{k})_{k=1,N}$. By employing
the canonical heliocentric variables of (Poincaré, 1896), the Hamiltonian of
the Newtonian gravitational interactions among the bodies reads

$H=\sum_{k=1}^{N}\left(\frac{\|\tilde{\bm{r}}_{k}\|^{2}}{2\mu_{k}}-G\frac{m_{0}m_{k}}{\|\bm{r}_{k}\|}\right)+\sum_{1\leq
k<\ell\leq
N}\left(\frac{\tilde{\bm{r}}_{k}\cdot\tilde{\bm{r}}_{\ell}}{m_{0}}-G\frac{m_{k}m_{\ell}}{\|\bm{r}_{k}-\bm{r}_{\ell}\|}\right),$
(1)

where $\bm{r}_{k}=\bm{R}_{k}-\bm{R}_{0}$ are the planet heliocentric
coordinates and $\tilde{\bm{r}}_{k}=\tilde{\bm{R}}_{k}$ their conjugated
momenta, $\mu_{k}=m_{0}m_{k}/(m_{0}+m_{k})$ are the reduced masses of the
planets and $G$ is the gravitational constant (Laskar, 1991; Laskar & Robutel,
1995). The Hamiltonian $H$ is a perturbation to the union of disjoint Kepler
problems, so that it is useful to introduce a set of canonical variables which
trivially integrates the unperturbed problems (e.g., Morbidelli, 2002). By
adopting from now on the momentum-coordinate ordering of conjugate pairs,
appropriate variables are $(\Lambda_{k},\lambda_{k};$ $x_{k},-j\bar{x}_{k};$
$y_{k},-j\bar{y}_{k})_{k=1,N}$, defined as

$\displaystyle\Lambda_{k}=\mu_{k}\sqrt{G(m_{0}+m_{k})a_{k}},$ (2)
$\displaystyle
x_{k}=\sqrt{\Lambda_{k}}\sqrt{1-\sqrt{1-e_{k}^{2}}}\,\mathrm{E}^{j\varpi_{k}},$
$\displaystyle
y_{k}=\sqrt{2\Lambda_{k}}\left(1-e_{k}^{2}\right)^{\frac{1}{4}}\sin(i_{k}/2)\,\mathrm{E}^{j\Omega_{k}},$

where $a_{k}$ are the planet semi-major axes, $\lambda_{k}$ the mean
longitudes, $e_{k}$ the eccentricities, $i_{k}$ the inclinations, $\varpi_{k}$
the longitudes of the perihelia and $\Omega_{k}$ the longitudes of the nodes
(Laskar, 1991; Laskar & Robutel, 1995). Throughout the paper, $j=\sqrt{-1}$
stands for the imaginary unit, $\mathrm{E}$ represents the exponential
operator and the overbar denotes the conjugate of a complex variable. The
variables $x_{k}$ and $y_{k}$ are the Poincaré’s rectangular coordinates in
complex form; we shall refer to them throughout the paper as Poincaré’s
complex variables, or simply Poincaré’s variables. With such a choice of
canonical variables, the integrable part of the Hamiltonian (1) reads

$H_{0}=\sum_{k=1}^{N}\left(\frac{\|\tilde{\bm{r}}_{k}\|^{2}}{2\mu_{k}}-G\frac{m_{0}m_{k}}{\|\bm{r}_{k}\|}\right)=-\sum_{k=1}^{N}\frac{G^{2}(m_{0}+m_{k})^{2}\mu_{k}^{3}}{2\Lambda_{k}^{2}},$
(3)

so that the Poincaré’s complex variables are constants of motion for the
Kepler problem.

In the regime of small orbital eccentricities and inclinations, which
characterizes the solar system, Eqs. (2) give
$x_{k}=\sqrt{\Lambda_{k}/2}e_{k}\mathrm{E}^{j\varpi_{k}}+\mathcal{O}(e_{k}^{3})$
and
$y_{k}=\sqrt{\Lambda_{k}/2}i_{k}\mathrm{E}^{j\Omega_{k}}+\mathcal{O}(e_{k}^{2}i_{k},i_{k}^{3})$,
and the Poincaré’s variables are also small. The principal part of the two-
body perturbation in Eq. (1) can be thus expanded as a Fourier series in the
planet mean longitudes, with polynomial coefficients depending on the
Poincaré’s variables (e.g., Laskar, 1990b, 1991; Laskar & Robutel, 1995),

$\displaystyle\mathcal{U}_{1}=-G\frac{mm^{\prime}}{\|\bm{r}-\bm{r}^{\prime}\|}=-G\frac{mm^{\prime}}{a^{\prime}}\sum_{\ell,\ell^{\prime}\in\mathbb{Z}}\widetilde{\mathcal{U}}_{\ell,\ell^{\prime}}\mathrm{E}^{j(\ell\lambda+\ell^{\prime}\lambda^{\prime})},$
(4)
$\displaystyle\widetilde{\mathcal{U}}_{\ell,\ell^{\prime}}=\sum\Gamma^{\ell,\ell^{\prime}}_{\mathcal{N}}\\!(\alpha)\mathcal{X}^{n}{\mathcal{X}^{\prime}}^{n^{\prime}}\bar{\mathcal{X}}^{\bar{n}}{\bar{\mathcal{X}}}^{\prime\bar{n}^{\prime}}\mathcal{Y}^{m}{\mathcal{Y}^{\prime}}^{m^{\prime}}\bar{\mathcal{Y}}^{\bar{m}}{\bar{\mathcal{Y}}}^{\prime\bar{m}^{\prime}},$

where
$\mathcal{N}=(n,n^{\prime},\bar{n},\bar{n}^{\prime},m,m^{\prime},\bar{m},\bar{m}^{\prime})$
is a tuple of non-negative integers222The symmetries of the planetary
Hamiltonian (1) imply constraints on the exponents
$(\ell,\ell^{\prime},\mathcal{N})$ known as D’Alembert rules (e.g., Laskar &
Robutel, 1995; Morbidelli, 2002). The rotational invariance requires
$n+n^{\prime}+m+m^{\prime}-\bar{n}-\bar{n}^{\prime}-\bar{m}-\bar{m}^{\prime}+\ell+\ell^{\prime}=0$,
while from the planar symmetry it follows that the terms of the series
$\widetilde{\mathcal{U}}_{\ell,\ell^{\prime}}$ are even with respect to the
variables
$(\mathcal{Y},\mathcal{Y}^{\prime},\bar{\mathcal{Y}},\bar{\mathcal{Y}}^{\prime})$..
Following (Laskar & Robutel, 1995), we have defined the dimensionless
Poincaré’s variables $\mathcal{X}=x\sqrt{2/\Lambda}$,
$\mathcal{Y}=y/\sqrt{2\Lambda}$ and the semi-major axis ratio
$\alpha=a/a^{\prime}$, with $a<a^{\prime}$. The analytical expression of the
coefficients $\Gamma^{\ell,\ell^{\prime}}_{\mathcal{N}}\\!(\alpha)$, 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,

$\mathcal{T}_{1}=\frac{\tilde{\bm{r}}\cdot\tilde{\bm{r}}^{\prime}}{m_{0}}=\sum_{\ell,\ell^{\prime}\in\mathbb{Z}}\widetilde{\mathcal{T}}_{\ell,\ell^{\prime}}(\Lambda,\Lambda^{\prime},\mathcal{X},\mathcal{X}^{\prime},\mathcal{Y},\mathcal{Y}^{\prime})\,\mathrm{E}^{j(\ell\lambda+\ell^{\prime}\lambda^{\prime})},$
(5)

The computation of the Fourier coefficients
$\widetilde{\mathcal{T}}_{\ell,\ell^{\prime}}$ is outlined in Appendix A.
Equations (4) and (5) allow to explicitly compute the Fourier expansion of the
Hamiltonian perturbing function,

$H_{1}=\sum_{1\leq k<\ell\leq
N}\left(\frac{\tilde{\bm{r}}_{k}\cdot\tilde{\bm{r}}_{\ell}}{m_{0}}-G\frac{m_{k}m_{\ell}}{\|\bm{r}_{k}-\bm{r}_{\ell}\|}\right)=\sum_{\bm{\ell}\in\mathbb{Z}^{N}}\widetilde{H}_{\bm{\ell}}\,\mathrm{E}^{j\bm{\ell}\cdot\bm{\lambda}},$
(6)

where $\bm{\lambda}$ stands for the vector of the planet mean longitudes,
$\bm{\lambda}=(\lambda_{1},\dots,\lambda_{8})$, and the coefficients
$\widetilde{H}_{\bm{\ell}}$ 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 ones, essentially consists of the slow precession of their perihelia and
nodes, driven by secular, i.e. orbit-averaged, gravitational interactions
(Laskar, 1990c; Laskar et al., 2004). A secular Hamiltonian, describing such
long-term motion, can be introduced in its simplest form by searching for a
change of variables that eliminate, at first order in the planet masses, all
the harmonics with non-null wavevectors $\bm{\ell}$ 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

$H_{1}+\\{S,H_{0}\\}=\langle H_{1}\rangle,$ (7)

where the braces represent the Poisson bracket. The angle-bracket operator
stands for averaging over the mean longitudes,

$\langle\cdot\rangle=\frac{1}{(2\pi)^{N}}\int_{\mathbb{T}^{N}}d\bm{\lambda}\,\,\cdot\,\,,$
(8)

with the integration defined over the hypertorus $\mathbb{T}^{N}$, at fixed
values of all the remaining canonical variables. This means that $\langle
H_{1}\rangle$ is the Fourier coefficient $\widetilde{H}_{\bm{0}}$
corresponding to the null harmonic in the expansion (6). The homologic
equation (7) gives the generating Hamiltonian $S$ as a formal Fourier series,

$S=-j\sum_{\bm{\ell}\in\mathbb{Z}^{N}\backslash\\{\bm{0}\\}}\frac{\widetilde{H}_{\bm{\ell}}}{\bm{\ell}\cdot\bm{n}}\mathrm{E}^{j\bm{\ell}\cdot\bm{\lambda}},$
(9)

where $\bm{n}=\partial H_{0}/\partial\bm{\Lambda}$ is the vector of the planet
mean motions. The secular Hamiltonian $\widehat{H}$ is formally given by the
Lie transform generated by the function $S$ and applied to the Hamiltonian
$H$,

$\widehat{H}=\mathrm{E}^{L_{S}}H\Big{|}_{\hat{\bm{\Lambda}},\hat{\bm{\lambda}},\hat{\bm{x}},\hat{\bm{y}}},\quad\mathrm{E}^{\pm
L_{S}}\,\cdot=\sum_{n=0}^{+\infty}\frac{(\pm 1)^{n}}{n!}L_{S}^{n}\,\cdot$ (10)

where $L_{S}\cdot=\\{S,\cdot\\}$ is the Lie derivative associated to the
generating Hamiltonian, $L_{S}^{0}$ is defined as the identity operator and
$L^{n}_{S}\cdot=L_{S}L_{S}^{n-1}\cdot$ for $n\geq 1$. The Hamiltonian
$\widehat{H}$ in Eq. (10) is expressed in the new canonical variables
$(\hat{\Lambda}_{k},\hat{\lambda}_{k};$ $\hat{x}_{k},-j\hat{\bar{x}}_{k};$
$\hat{y}_{k},-j\hat{\bar{y}}_{k})_{k=1,N}$, which we shall call the secular
variables. They are related to the original ones via the Lie transforms

$\Lambda_{k}=\mathrm{E}^{L_{S}}\hat{\Lambda}_{k},\quad\lambda_{k}=\mathrm{E}^{L_{S}}\hat{\lambda}_{k},\quad
x_{k}=\mathrm{E}^{L_{S}}\hat{x}_{k},\quad
y_{k}=\mathrm{E}^{L_{S}}\hat{y}_{k}.$ (11)

The original variables are therefore the superposition of the secular ones and
short-time oscillations generated by the $n\geq 1$ terms of the Lie
transforms. Since in the present study we are interested in the long-term
dynamics of the planets, we shall only focus on the secular variables.
Therefore, to keep simpler notation, we shall omit the hat on the secular
variables from now on.

Differently from the outer planets, the inner ones are not currently involved
in strong mean-motion resonances. More precisely, a maximal contribution of
only 0.07$\arcsec\,{\textrm{yr}}^{-1}$ to the fundamental precession
frequencies of the inner orbits in the Laplace-Lagrange solution arises from
the second order in the planet masses (Laskar, 1985, Table 8). This same
contribution is 0.9$\arcsec\,{\textrm{yr}}^{-1}$ for the outer planets, mainly
due to the Great Inequality. Building on this fact, in the present work we
choose to truncate the series (10) at first order in the planet masses, i.e.
we neglect quadratic and higher-order terms with respect to the Fourier
coefficients $\widetilde{H}_{\bm{\ell}}$. Indeed, we expect the main
contribution to the precession frequencies of the inner orbits to come from
the linear terms, given the small masses of the inner planets. In absence of
strong mean-motion resonances, in Eq. (9) the denominators
$\bm{\ell}\cdot\bm{n}$ involving at least one inner planet are sufficiently
far from zero. Under these assumptions, the higher-order terms only produce
small corrections to the dynamics generated by the leading ones. Using the
homologic equation (7), one thus obtains

$\widehat{H}=H_{0}+\langle H_{1}\rangle.$ (12)

The resulting secular Hamiltonian $\widehat{H}$ is simply the average of the
N-body Hamiltonian $\eqref{eq:NbodyHam}$ over the planet mean longitudes. The
averaging process is mathematically equivalent to replacing each planet by its
instantaneous Keplerian orbit, with the corresponding mass distributed along
it in a way that is inversely proportional to the local orbital speed of the
planet. The secular dynamics is thus the slow gravitational interaction of
such Keplerian rings. This equivalence was pointed out by Gauss (1818) and it
is thus referred to as Gauss’ dynamics. As a result of the averaging over the
mean longitudes, the $\Lambda_{k}$ variables are constant of motion in the
secular dynamics, and so are the semi-major axes of the planets. The resulting
Hamiltonian system consists of two degrees of freedom for each planet333The
entire planetary system is therefore described by 16 degrees of freedom, even
reduced to 15 when taking into account the conservation of the total angular
momentum., corresponding to the two Poincaré’s complex variables $x_{k}$ and
$y_{k}$.

##### General relativity and minor effects.

The dynamical interactions described by the Hamiltonian (1) are not sufficient
to finely reproduce the precession frequencies of the inner orbits; to this
end, additional physical effects must be taken into account (Laskar, 1999).
Indeed, it is well-known that general relativity contributes with
0.430$\arcsec\,{\textrm{yr}}^{-1}$ to the secular precession of Mercury
perihelion. Such a correction is critical for the statistics of the long-term
destabilization of the inner orbits (Laskar, 2008; Laskar & Gastineau, 2009),
since it moves the system away from the $g_{1}-g_{5}$ secular resonance, which
is responsible for the very high eccentricities of Mercury (Laskar, 2008;
Batygin & Laughlin, 2008; Boué et al., 2012). We shall therefore include in
the Hamiltonian (12) the leading secular contribution of general relativity,
which reads

$\widehat{H}_{\textrm{GR}}=\sum_{k=1}^{N}\frac{G^{2}m_{0}^{2}m_{k}}{c^{2}a_{k}^{2}}\left(\frac{15}{8}-\frac{3}{\sqrt{1-e_{k}^{2}}}\right),$
(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$\arcsec\,{\textrm{yr}}^{-1}$
contribution to the secular perihelion precession of the Earth (Laskar, 1999).
However, this is the order of magnitude of the contribution coming from the
Hamiltonian terms at second order in the planet masses, which are neglected in
the present model. Indeed, at degree 2 in planet eccentricities and
inclinations, (Laskar, 1985, Table 8) reported contributions of
0.073$\arcsec\,{\textrm{yr}}^{-1}$ and 0.062$\arcsec\,{\textrm{yr}}^{-1}$ to
the Earth-dominated and Mars-dominated eccentricity modes $g_{3}$ and $g_{4}$,
respectively, the other frequency corrections being at least ten times
smaller. Since the denominators $\bm{\ell}\cdot\bm{n}$ appearing in Eq. (9)
and involving the inner planets are not close to zero, we expect the
contributions from higher degrees to be generally smaller than those from
degree 2. We thus choose to not include the Earth-Moon interaction in the
present model, preferring to deal with a Hamiltonian whose dependence on
eccentricities and inclinations is exact, as compared to a non-decisive
increase of precision on the orbit precession frequencies. Moreover, we know
that the N-body Hamiltonian (1), corrected for general relativity, reproduces
the maximum Lyapunov exponent of the inner solar system (Rein & Tamayo, 2018).

Therefore, apart from irrelevant constant terms only depending on the semi-
major axes in Eqs. (3) and (13), the secular Hamiltonian of the $8$ solar
system planets considered in the present study reads

$\widehat{H}=-\sum_{k=1}^{8}G\frac{m_{0}m_{k}}{a_{k}}\left(\sum_{\ell=1}^{k-1}\frac{m_{\ell}}{m_{0}}\left<\frac{a_{k}}{\|\bm{r}_{k}-\bm{r}_{\ell}\|}\right>+3\frac{Gm_{0}}{c^{2}a_{k}}\frac{1}{\sqrt{1-e_{k}^{2}}}\right)$
(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_{5}$ and $g_{6}$, respectively. Fortunately enough, the
very small variations over billions of years of the precession frequencies of
the outer orbits, compared to those of the inner one (Laskar, 1990c; Laskar et
al., 2004; Hoang et al., 2021), naturally suggest a way to make out of the
Hamiltonian (14) a very effective model for the inner system. This is achieved
by choosing, once for all, an explicit quasi-periodic time dependence for the
orbits of the giant planets, i.e. by expressing the corresponding Poincaré’s
complex variables as finite Fourier series,

$\left.\begin{aligned}
&x_{k}(t)=\sum_{\ell=1}^{M_{k}}\tilde{x}_{k\ell}\,\mathrm{E}^{j\bm{m}_{k\ell}\cdot\bm{\phi}(t)}\\\
&y_{k}(t)=\sum_{\ell=1}^{N_{k}}\tilde{y}_{k\ell}\,\mathrm{E}^{j\bm{n}_{k\ell}\cdot\bm{\phi}(t)}\end{aligned}\;\right\\}\;k\in\\{5,6,7,8\\}\,,$
(15)

where $t$ denotes the time, $\tilde{x}_{k\ell}$ and $\tilde{y}_{k\ell}$ are
complex amplitudes, $\bm{m}_{k\ell}$ and $\bm{n}_{k\ell}$ integer vectors, and
$\bm{\phi}(t)=\bm{\omega}_{\textnormal{\scriptsize o}}t$, with
$\bm{\omega}_{\textnormal{\scriptsize o}}$ the vector of the (constant)
fundamental precession frequencies of the outer orbits, denoted as
$\bm{\omega}_{\textnormal{\scriptsize
o}}=(g_{5},g_{6},g_{7},g_{8},s_{6},s_{7},s_{8})$444The frequency of the
Jupiter-dominated inclination mode $s_{5}$ is null because the total angular
momentum of the system is conserved. (Laskar, 1990c). The number of harmonics
$M_{k},N_{k}$ appearing in the decompositions depends on the planet. When such
predetermined time dependence is injected in Eq. (14), one obtains the
Hamiltonian $\mathcal{H}$ of a _forced secular inner solar system_ ,

$\displaystyle\mathcal{H}[(x_{k},y_{k})_{k=1,4},t]\\!=\\!\widehat{H}[(x_{k},y_{k})_{k=1,4},(x_{k}\\!\\!=\\!\\!x_{k}(t),y_{k}\\!\\!=\\!\\!y_{k}(t))_{k=5,8}].$
(16)

The explicit time dependence in the Hamiltonian $\mathcal{H}$ physically means
that the inner planets interact with each other while moving in an external
gravitational potential generated by the giant planets555The terms in
$\widehat{H}$ only involving outer planet variables can be discarded in Eq.
(16), as they just depend on time and do not affect the inner planet dynamics
anymore.. 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
$\mathcal{H}$ possesses 8 degrees of freedom. As already stated, compared to
the second-order secular system of (Laskar, 1984, 1985), the present model
neglects terms of order higher than one in the planet masses. Nevertheless,
the precession frequencies of the outer orbits can be set to very precise
values in Eq. (15). Moreover, the Hamiltonian (16) is not truncated neither in
eccentricities nor in inclinations, so that the range of validity of its
dynamics extends to very excited states, provided that the Keplerian orbits of
the planets do not cross each other. This makes the present model perfectly
suited to explore the very long-term evolution of the inner system, even when
high eccentric and inclined orbits become statistically recurrent.

Table 1: Average wall-clock time on an Intel(R) Core(TM) i7-7700T CPU @ 2.90GHz for the numerical integration of the forced inner system according to the total degree of truncation of the Hamiltonian in eccentricities and inclinations ($\mathcal{H}_{2n}$ is given in Eq. (25), while $\mathcal{H}$ in Eq. (16)). $\mathcal{H}_{2n}$ | 4 | 6 | 8 | 10 | $\mathcal{H}$
---|---|---|---|---|---
100-Myr orbital solution | 1.3 s | 7.5 s | 35 s | 4 m 01 s | 2 m 48 s

###  2.3 Construction of the outer planet solution. Initial conditions

The quasi-periodic form of the outer orbits in Eq. (15) is established
numerically and explicitly reported in Appendix D. The full equations of
motion of the main bodies of the solar system are numerically integrated over
30 Myr in the future, following the comprehensive model of (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) extending over 1 Myr. Throughout the paper, we
shall refer to this direct numerical integration as LaX13b. Through frequency
analysis (Laskar, 1988; Laskar et al., 1992; Laskar, 1993, 2005) of the
orbital solution, the leading secular harmonics of the dimensionless
Poincaré’s variables of the outer planets
$(\mathcal{X}_{k},\mathcal{Y}_{k})_{k=5,8}$ are extracted to construct the
Fourier series (15), as in (Laskar, 1988, 1990c). This is performed up to a
numerical precision that does not allow anymore to recognize in an unambiguous
way new harmonics as linear combination of the fundamental frequencies
$\bm{\omega}_{\textnormal{\scriptsize o}}$. The precision of the outer planet
solution is shown in Table 10 of Appendix D. We report there the root mean
square of the dimensionless Poincaré’s variables of the outer planets in the
LaX13b solution, and that of the corresponding residuals after subtraction of
their quasi-periodic decomposition from Eq. (15) and filtering out of the non-
secular Fourier components. The effectiveness of the quasi-periodic
approximation is illustrated in Figs. 9 and 10.

The choice of the secular semi-major axes and that of the initial conditions
for the secular Poincaré’s variables of the inner planets deserve a
discussion. In principle, the initial secular variables have to be computed by
inverting the Lie transforms in Eq. (11). They should not be set to the
initial value of the corresponding original variable, as this could cause an
offset in the secular frequencies of the motion, due to the short-time (high-
frequency) oscillations generated by the Lie transform. This is a general rule
when constructing averaged dynamical systems (e.g., Laskar & Simon, 1988). To
avoid the explicit computation of the generating function $S$ in Eq. (9) and
the related Lie transforms, appropriate initial conditions for the secular
solar system can be effectively computed by filtering out the short-time
components of the numerical solution of the non-averaged system, as it is done
in the present study. The constant term in the frequency analysis of the
variables $\Lambda_{k}$, performed on the LaX13b solution, provides the value
of the secular semi-major axes of the planets (Table 7). 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é’s complex variables of the inner planets (Table 8).

##  3 Numerical integration. Gauss’s method in Hamiltonian formalism

The dimensionless Hamilton’s equations for the forced inner system (16) read

$\dot{\mathcal{X}}_{k}=-\frac{2j}{\Lambda_{k}}\frac{\partial\mathcal{H}}{\partial\bar{\mathcal{X}}_{k}},\quad\dot{\mathcal{Y}}_{k}=-\frac{j}{2\Lambda_{k}}\frac{\partial\mathcal{H}}{\partial\bar{\mathcal{Y}}_{k}},$
(17)

for $k\in\\{1,2,3,4\\}$. Because of the averaging over the mean longitudes in
Eq. (14), the above equations include double integrals of the form

$\displaystyle\frac{\partial}{\partial\bar{\mathcal{Z}}}\left<\frac{1}{\|\bm{r}-\bm{r}^{\prime}\|}\right>=\frac{1}{2\pi}\int_{0}^{2\pi}\bm{f}\cdot\frac{\partial\bm{r}}{\partial\bar{\mathcal{Z}}}d\lambda,$
(18)
$\displaystyle\bm{f}=\frac{1}{2\pi}\int_{0}^{2\pi}-\frac{\bm{r}-\bm{r}^{\prime}}{\|\bm{r}-\bm{r}^{\prime}\|^{3}}d\lambda^{\prime},$
(19)

where $\mathcal{Z}$ stands for $\mathcal{X}$ or $\mathcal{Y}$, alternatively.
The equations of motion (17) are thus non-algebraic, and the numerical
computation of the derivatives is in principle much more complex than in
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’s method
(Gauss, 1818; Bour, 1855; Hill, 1882). Indeed, the vector $\bm{f}$ in Eq. (19)
is proportional to the gravitational force exerted on a test particle by a
Keplerian ring, and can be analytically expressed in terms of complete
Legendre’s elliptic integrals (e.g., Olver et al., 2020, Chapter 19). This was
first shown by Gauss (1818), who introduced at same time the arithmetic-
geometric mean to numerically evaluate such special functions in a few
elementary iterations and with high precision. Building on its modern
derivation in (Musen, 1970), we implemented Gauss’s method into the present
complex Hamiltonian formalism (17). This allows, in particular, to eliminate
the degeneracy typically appearing at $e=0$ (circular orbits) and $i=0$
(equatorial orbits) (e.g., Fouvry et al., 2020), which is fundamental to
guarantee high numerical precision. We employ a state-of-the-art algorithm,
based on piecewise minimax rational function approximation, to numerically
compute the complete elliptic integrals in double floating-point precision at
the cost of elementary functions (Fukushima, 2015). The remaining simple
integral in Eq. (18) is then effectively evaluated via the trapezoidal rule,
which converges exponentially fast with the number of function evaluations,
because of the periodicity of the integrand (e.g., Touma et al., 2009). We
apply the trapezoidal rule in an adaptive way, by doubling the number of
function evaluations until the estimated relative error on the integral is
smaller than $10^{-12}$. It should be noted that, as the convergence of the
numerical integral is exponential, the resulting error is often orders of
magnitude smaller than this tolerance.

We integrate the Eqs. (17) using an Adams PECE method of order 12 and a
conservative timestep of 250 years, as in (Laskar, 1994). In absence of any
fundamental integral of motion (see Sect. 2.2), we estimate the integration
error after a time $T$ following (Laskar, 1994). Starting from the nominal
initial conditions, we integrate the dynamics over the time interval [0,
$T$/2] and then backwards to the initial time. By denoting the deviations from
the initial coordinates of the system in the phase space as
$(\delta\mathcal{X}_{k}(T),\delta\mathcal{Y}_{k}(T))_{k=1,4}$, we define the
relative integration error as

$\displaystyle\delta(T)=\max_{k\in\\{1,\dots,4\\}}\Bigg{\\{}$
$\displaystyle\left|\frac{\operatorname{Re}[\delta\mathcal{X}_{k}(T)]}{\operatorname{Re}[\mathcal{X}_{k}(0)]}\right|,\left|\frac{\operatorname{Im}[\delta\mathcal{X}_{k}(T)]}{\operatorname{Im}[\mathcal{X}_{k}(0)]}\right|,$
(20)
$\displaystyle\left|\frac{\operatorname{Re}[\delta\mathcal{Y}_{k}(T)]}{\operatorname{Re}[\mathcal{Y}_{k}(0)]}\right|,\left|\frac{\operatorname{Im}[\delta\mathcal{Y}_{k}(T)]}{\operatorname{Im}[\mathcal{Y}_{k}(0)]}\right|\Bigg{\\}},$

where $(\mathcal{X}_{k}(0),\mathcal{Y}_{k}(0))_{k=1,4}$ are the initial
conditions of the system. Over the time interval $T\in$ [0, 20] Myr, the
average integration error is well described by the power law
$\delta(T)=\delta_{0}[1+(T/T_{0})^{\alpha}]$, with $\delta_{0}=9\cdot
10^{-15}$, $T_{0}=0.6$ Myr and $\alpha=1.2$. Therefore, one has
$\delta(10\,\mathrm{Myr})\approx 3\cdot 10^{-13}$, and
$\delta(10\,\mathrm{Gyr})\approx 10^{-9}$ if one extrapolates by ignoring the
chaotic behaviour of the solution. The integration error thus turns out to be
similar to that of the secular system in (Laskar, 1994). We emphasize that the
present fixed-timestep integration scheme allows to reach very highly excited
orbital states. Indeed, violations of the precision goal on the integral
$\eqref{eq:gauss_double_integral_1}$ typically occurs only when the system is
already close to the first intersection of the instantaneous Keplerian orbits
of two planets (see Sect. 7), and are related to an intrinsic discontinuity in
the equations of motion (17) at orbit crossing (e.g., Touma et al., 2009).

From a computational perspective, the choice of a multistep method is well
adapted to the present dynamical system, as the derivative evaluation
constitutes by far the most expensive step of the integration scheme. Once
implemented in the C programming language, Gauss’s method allows to compute
100-Myr orbital solutions of the forced inner planets in a few minutes on a
PC, as shown in Table 1. Therefore, the computational cost turns out to be
comparable to that of the secular system of Laskar (1990c), in spite of the
different complexity of the corresponding Hamiltonians.

The rather involved analytical derivations adapting Gauss’s method to the
present complex Hamiltonian formalism shall be the subject of a forthcoming
paper (Mogavero, 2021), along with the release of the C program implementing
it.

##  4 Analytical expansion. Computer algebra

In addition to its fast numerical integration, the present dynamical model
allows for a systematic analytical development of its Hamiltonian by means of
computer algebra. In this study we employ TRIP, a computer algebra system
dedicated to perturbation series, specially those of celestial mechanics
(Laskar, 1990a; Gastineau & Laskar, 2011, 2020). The main objects of its
symbolic kernel are the Poisson series, i.e. multivariate Fourier series whose
coefficients are multivariate Laurent series,

$S(z_{1},\dots,z_{n},\varphi_{1},\dots,\varphi_{m})=\sum
C_{k,\ell}z_{1}^{k_{1}}\cdots
z_{n}^{k_{n}}\mathrm{E}^{j\left(\ell_{1}\varphi_{1}+\dots+\ell_{m}\varphi_{m}\right)},$
(21)

where $(z_{p})_{p=1,n}$ and $(\varphi_{p})_{p=1,m}$ are complex and real
variables, respectively, $k=(k_{p})_{p=1,n}\in\mathbb{Z}^{n}$,
$\ell=(\ell_{p})_{p=1,m}\in\mathbb{Z}^{m}$ and $C_{k,\ell}$ are complex
coefficients.

Following (Laskar & Robutel, 1995), the smallness of eccentricities and
inclinations in the solar system can be exploited to develop the averages
appearing in Eq. (14) as formal series in Poincaré’s variables
$(x_{k},\bar{x}_{k},y_{k},\bar{y}_{k})_{k=1,8}$. Indeed, Eq. (4) gives

$\left<\frac{a^{\prime}}{\|\bm{r}^{\prime}-\bm{r}\|}\right>=\sum\Gamma^{0,0}_{\mathcal{N}}(\alpha)\mathcal{X}^{n}{\mathcal{X}^{\prime}}^{n^{\prime}}\bar{\mathcal{X}}^{\bar{n}}{\bar{\mathcal{X}}}^{\prime\bar{n}^{\prime}}\mathcal{Y}^{m}{\mathcal{Y}^{\prime}}^{m^{\prime}}\bar{\mathcal{Y}}^{\bar{m}}{\bar{\mathcal{Y}}}^{\prime\bar{m}^{\prime}}.$
(22)

The rotational invariance of the secular Hamiltonian (14) requires
$n+n^{\prime}+m+m^{\prime}-\bar{n}-\bar{n}^{\prime}-\bar{m}-\bar{m}^{\prime}=0$,
while from the planar symmetry it follows that
$m+m^{\prime}+\bar{m}+\bar{m}^{\prime}$ is an even integer. These relations
imply that the monomials in Eq. (22) are even with respect to both
$(x,x^{\prime},\bar{x},\bar{x}^{\prime})$ and
$(y,y^{\prime},\bar{y},\bar{y}^{\prime})$ variables. The relativistic terms in
Eq. (14) are readily expanded,

$\frac{1}{\sqrt{1-e^{2}}}=\frac{1}{1-x\bar{x}/\Lambda}=\sum_{p=0}^{\infty}\left(\frac{x\bar{x}}{\Lambda}\right)^{p}.$
(23)

By truncating the series (22) and (23) at a given total degree $2n$
($n\in\mathbb{N}_{0}$) in Poincaré’s variables
$(x_{k},\bar{x}_{k},y_{k},\bar{y}_{k})_{k=1,8}$, and after substitution in Eq.
(14), one obtains a polynomial secular Hamiltonian for the ensemble of the
solar system planets,

$\widehat{H}_{2n}=\sum_{p=0}^{n}\widehat{H}_{(2p)}\,,$ (24)

where $\widehat{H}_{(2p)}$ groups all the monomials of same total degree $2p$.
Practically, the expansion (24) is readily obtained in TRIP, which implements
the algorithm of (Laskar & Robutel, 1995) to compute a truncation of the
series (22). The truncated Hamiltonian (24) is in the form of a Poisson series
(21) in Poincaré’s complex variables, with no angular dependencies. The
crucial point here is that its coefficients only depend on the planet masses
and secular semi-major axes, which are constant parameters in the present
dynamics. They can be thus numerically evaluated once for all, to obtain very
compact series. This is shown in Table 2, which gives the total number of
monomials in (24) with respect to the degree of truncation, along with the
disk usage of the corresponding series when stored in a plain-text file, and
their typical computation time in TRIP.

The expansion (24) straightforwardly provides a truncated Hamiltonian for the
forced inner system (16),

$\displaystyle\mathcal{H}_{2n}=\sum_{p=1}^{n}\mathcal{H}_{(2p)}\,,$ (25)
$\displaystyle\mathcal{H}_{(2p)}[(x_{k},y_{k})_{k=1,4},t]\\!=\\!\widehat{H}_{(2p)}[(x_{k},y_{k})_{k=1,4},(x_{k}\\!\\!=\\!\\!x_{k}(t),y_{k}\\!\\!=\\!\\!y_{k}(t))_{k=5,8}].$

The computation of the equations of motion for the truncated Hamiltonian (25)
only requires to take the derivatives of a multivariate polynomial, and can be
thus systematically performed in TRIP. We put the resulting polynomials in
Horner form, to ensure speed and stability of their numerical evaluation, and
wrap them into a C code to achieve the best computational
performance666Numerical integration of ODEs is also available directly in
TRIP.. The equations of motion are then integrated through the Adams PECE
method of order 12 employed for Gauss’s dynamics, with the same timestep of
250 years. The typical computational cost of a 100-Myr orbital solution is
shown in Table 1, according to the degree of truncation. In light of the
pertinence of the forced inner system that will be shown in Sect. 5, the
truncated Hamiltonian (25), at degrees 4 and 6 in particular, constitutes the
state of the art of very fast, still realistic, dynamical models of the inner
solar system.

Table 2: Number of monomials of variables $(x_{k},\bar{x}_{k},y_{k},\bar{y}_{k})_{k=1,8}$ in $\widehat{H}_{2n}$ (Eqs. 2, 24) and of variables $\\{(u_{k},\bar{u}_{k},v_{k},\bar{v}_{k})_{k=1,4},(\mathrm{E}^{j\phi_{k}(t)})_{k=1,7}\\}$ in $\mathbb{H}_{2n}$ (Eqs. 42, 44, 48) according to the truncation degree. It also shows the disk space in Bytes occupied by the series, their computation time on an Intel(R) Xeon(R) CPU E5-2698 v3 @ 2.30GHz (CPU time in parentheses) and the number of Fourier harmonics in $\mathbb{H}_{2n}$ according Eq. (52). | 2 | 4 | 6 | 8 | 10
---|---|---|---|---|---
$\widehat{H}_{2n}$ |  |  |  |  |
Monomials | 129 | 1 345 | 6 561 | 23 213 | 66 253
Disk usage | 21 kB | 215 kB | 1.0 MB | 3.7 MB | 10.6 MB
Wall-clock time | 2.4 s (137%) | 2.5 s (129%) | 3.0 s (154%) | 5.7 s (187%) | 17.1 s (272%)
$\mathbb{H}_{2n}$ |  |  |  |  |
Monomials | 8 | 6 304 | 188 024 | 3 394 892 | 42 817 100
Disk usage | – | 0.9 MB | 29.5 MB | 547.6 MB | 7.1 GB
Wall-clock time | – | 0.15 s (447%) | 3.4 s (773%) | 1 m 16 s (1038%) | 20 m 37 s (1343%)
Harmonics | 1 | 2 748 | 69 339 | 1 029 137 | 10 279 581

###  4.1 Forced Laplace-Lagrange dynamics

For small eccentricities and inclinations, the secular Hamiltonian (14) is a
perturbation to the integrable Laplace-Lagrange (LL) problem. Indeed, its
truncation at degree 2 in Poincaré’s variables gives the quadratic form

$\widehat{H}_{\textnormal{\scriptsize
LL}}=\widehat{H}_{(2)}=\bm{x}^{\dagger}\bm{M}\bm{x}+\bm{y}^{\dagger}\bm{N}\bm{y},$
(26)

where we have defined the column vectors $\bm{x}=(x_{1},\dots,x_{8})^{T}$ and
$\bm{y}=(y_{1},\dots,y_{8})^{T}$, $T$ is the transposition operator and the
dagger stands for Hermitian transposition, i.e.
$\bm{x}^{\dagger}=\bar{\bm{x}}^{T}$. The matrices $\bm{M}$ and $\bm{N}$ are
real and symmetric. Since the variables $\bm{x}$ and $\bm{y}$ are not coupled
in the LL Hamiltonian, the following derivations focus on the degrees of
freedom related to $\bm{x}$. A similar treatment holds for the $\bm{y}$
variables.

After substitution of the time dependence of the outer planet variables given
in Eq. (15), one obtains

$\mathcal{H}_{\textnormal{LL}(\bm{x})}=\bm{x}^{\dagger}\bm{M}\bm{x}\Big{|}_{(x_{k}=x_{k}(t))_{k=5,8}}=\left(\bm{x}_{\textnormal{\scriptsize
i}}^{\dagger},\bm{x}_{\textnormal{\scriptsize
o}}^{\dagger}(t)\right)\begin{pmatrix}\bm{M}_{\textnormal{\scriptsize
ii}}&\bm{M}_{\textnormal{\scriptsize io}}\\\ \bm{M}_{\textnormal{\scriptsize
io}}^{T}&\bm{M}_{\textnormal{\scriptsize
oo}}\end{pmatrix}\begin{pmatrix}\bm{x}_{\textnormal{\scriptsize i}}\\\
\bm{x}_{\textnormal{\scriptsize o}}(t)\end{pmatrix},$ (27)

where we have defined the inner and outer planet column vectors,
$\bm{x}_{\textnormal{\scriptsize i}}=(x_{1},x_{2},x_{3},x_{4})^{T}$ and
$\bm{x}_{\textnormal{\scriptsize
o}}(t)=(x_{5}(t),x_{6}(t),x_{7}(t),x_{8}(t))^{T}$, respectively, and the real
matrix $\bm{M}$ has been written as a block matrix, with
$\bm{M}_{\textnormal{\scriptsize ii}}$ and $\bm{M}_{\textnormal{\scriptsize
oo}}$ symmetric $4\times 4$ matrices. By discarding terms only depending on
time, the Hamiltonian reads

$\mathcal{H}_{\textnormal{LL}(\bm{x})}=\bm{x}_{\textnormal{\scriptsize
i}}^{\dagger}\bm{M}_{\textnormal{\scriptsize
ii}}\bm{x}_{\textnormal{\scriptsize
i}}+2\operatorname{Re}\left(\bm{x}_{\textnormal{\scriptsize
i}}^{\dagger}\bm{M}_{\textnormal{\scriptsize
io}}\bm{x}_{\textnormal{\scriptsize o}}(t)\right),$ (28)

where $\operatorname{Re}$ stands for the real part of a complex quantity. The
Hamiltonian (28) corresponds to a forced Laplace-Lagrange system. Since it is
quadratic in the inner vector $\bm{x}_{\textnormal{\scriptsize i}}$, one can
analytically solve the corresponding equations of motion and introduce
appropriate action-angle variables.

####  4.1.1 Solution of the equations of motion

The matrix $\bm{M}_{\textnormal{\scriptsize ii}}$ is real and symmetric and
can thus be diagonalized through an orthogonal matrix $\bm{O}_{\text{M}}$,

$\bm{M}_{\textnormal{\scriptsize
ii}}=\bm{O}_{\text{M}}\bm{D}_{\text{M}}\bm{O}_{\text{M}}^{T}.$ (29)

The columns of $\bm{O}_{\text{M}}$ are the eigenvectors of
$\bm{M}_{\textnormal{\scriptsize ii}}$, while the diagonal entries of
$\bm{D}_{\text{M}}=-\mathrm{diag}(\bm{g_{\textnormal{\scriptsize LL}}})$ are
the corresponding real eigenvalues777The minus sign in the definition of the
matrix $\bm{D}_{\text{M}}$ is such that the precession frequencies of the
perihelia $\dot{\bm{\varpi}}$ have the same sign as
$\bm{g_{\textnormal{\scriptsize LL}}}$., $\bm{g_{\textnormal{\scriptsize
LL}}}$ being a column vector. The orthogonal matrix $\bm{O}_{\text{M}}$
induces the canonical change of variables888In the present notation, the
canonical momenta constitute column vectors, while the coordinates are wrapped
in row vectors. $(\bm{x}_{\textnormal{\scriptsize
i}},-j{\bm{x}^{\dagger}_{\textnormal{\scriptsize
i}}})\rightarrow(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j{\bm{x}_{\textnormal{\scriptsize i}}^{\prime}}^{\dagger})$
defined as

$\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}=\bm{O}_{\text{M}}^{T}\bm{x}_{\textnormal{\scriptsize i}}.$ (30)

The transformed Hamiltonian reads

$\mathcal{H}_{\textnormal{LL}(\bm{x})}={\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger}\bm{D}_{\text{M}}\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}+2\operatorname{Re}\left({\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger}\bm{M}_{\textnormal{\scriptsize
io}}^{\prime}\bm{x}_{\textnormal{\scriptsize o}}(t)\right),$ (31)

with $\bm{M}_{\textnormal{\scriptsize
io}}^{\prime}=\bm{O}_{\text{M}}^{T}\bm{M}_{\textnormal{\scriptsize io}}$. The
corresponding Hamilton’s equations are given by

$\bm{\dot{x}}_{\textnormal{\scriptsize
i}}^{\prime}=-j\frac{\partial\mathcal{H}_{\textnormal{LL}(\bm{x})}}{\partial{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger}}=-j\left(\bm{D}_{\text{M}}\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}+\bm{M}_{\textnormal{\scriptsize
io}}^{\prime}\bm{x}_{\textnormal{\scriptsize o}}(t)\right),$ (32)

and constitute a first-order inhomogeneous matrix ordinary differential
equation. The general solution can be written as

$\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}(t)=\bm{x}_{\textnormal{\scriptsize
i,F}}^{\prime}(t)+\bm{x}_{\textnormal{\scriptsize i,f}}^{\prime}(t).$ (33)

The free solution $\bm{x}_{\textnormal{\scriptsize i,F}}^{\prime}(t)$ is the
general integral to the associated homogeneous equation
$\bm{\dot{x}}_{\textnormal{\scriptsize
i}}^{\prime}=-j\bm{D}_{\text{M}}\bm{x}_{\textnormal{\scriptsize i}}^{\prime}$,
representing the autonomous perihelia precession of the inner orbits, while
the forced solution $\bm{x}_{\textnormal{\scriptsize i,f}}^{\prime}(t)$ is a
particular integral of the complete Eq. (32), which arises from the
gravitational forcing of the outer planets. We define, once for all, the
forced solution to be:

$\bm{x}_{\textnormal{\scriptsize
i,f}}^{\prime}(t)=-j\,\mathrm{E}^{-jt\bm{D}_{\text{M}}}\int^{t}d\tau\,\mathrm{E}^{j\tau\bm{D}_{\text{M}}}\bm{M}_{\textnormal{\scriptsize
io}}^{\prime}\bm{x}_{\textnormal{\scriptsize o}}(\tau).$ (34)

By employing the decomposition given in Eq. (15), the components of the forced
solution are

$(\bm{x}_{\textnormal{\scriptsize
i,f}}^{\prime})_{k}=\sum_{\ell=5}^{8}\sum_{p=1}^{M_{\ell}}\frac{(\bm{M}_{\textnormal{\scriptsize
io}}^{\prime})_{k\ell}\tilde{x}_{\ell p}}{(\bm{g_{\textnormal{\scriptsize
LL}}})_{k}-\bm{m}_{\ell p}\cdot\bm{\omega}_{\textnormal{\scriptsize
o}}}\mathrm{E}^{j\bm{m}_{\ell p}\cdot\bm{\phi}}\quad k\in\\{1,2,3,4\\}.$ (35)

The (constant) denominators appearing in Eq. (35) are far from zero, as the
inner planets are not involved in the corresponding secular resonances. The
forced solution is thus well defined.

####  4.1.2 Proper modes

The following derivation shows that there exists a canonical transformation
depending on time, $(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger})\rightarrow(\bm{u},-j\bm{u}^{\dagger})$, such that the
transformed LL Hamiltonian reads
$\mathbb{H}_{\textnormal{LL}(\bm{u})}=\bm{u}^{\dagger}\bm{D}_{\text{M}}\bm{u}$.
As in the case of an autonomous Laplace-Lagrange system, the new canonical
variables $\bm{u}$ will physically correspond to the free part
$\bm{x}_{\textnormal{\scriptsize i,F}}^{\prime}(t)$ of the solution in Eq.
(33) (e.g., Morbidelli, 2002). We thus begin by defining a new set of
variables $\bm{u}$ such that

$\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}=\bm{u}+\bm{x}_{\textnormal{\scriptsize i,f}}^{\prime}(t).$ (36)

By using the fact that $\bm{x}_{\textnormal{\scriptsize i,f}}^{\prime}(t)$ is
a solution to Eq. (32), and discarding terms only depending on time, the
forced LL Hamiltonian in Eq. (31) can be written as

$\mathcal{H}_{\textnormal{LL}(\bm{x})}=\bm{u}^{\dagger}\bm{D}_{\text{M}}\bm{u}+j\bm{u}^{\dagger}\dot{\bm{x}}_{\textnormal{\scriptsize
i,f}}^{\prime}(t)-j{\dot{\bm{x}}_{\textnormal{\scriptsize
i,f}}}^{\prime\dagger}(t)\bm{x}_{\textnormal{\scriptsize i}}^{\prime}.$ (37)

We now ask that the change of variables $(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger})\rightarrow(\bm{u},-j\bm{u}^{\dagger})$ derive from a
time-depending generating function $F(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j\bm{u}^{\dagger},t)$, satisfying

$\frac{\partial F}{\partial
t}=-j\bm{u}^{\dagger}\dot{\bm{x}}_{\textnormal{\scriptsize
i,f}}^{\prime}(t)+j{\dot{\bm{x}}_{\textnormal{\scriptsize
i,f}}}^{\prime\dagger}(t)\bm{x}_{\textnormal{\scriptsize i}}^{\prime},$ (38)

so that the last two terms in Eq. (37) shall be killed in the transformed
Hamiltonian. By integrating the above equation with respect to time, one
obtains

$F(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j\bm{u}^{\dagger},t)=-j\bm{u}^{\dagger}\bm{x}_{\textnormal{\scriptsize
i,f}}^{\prime}(t)+j\bm{x}^{\prime\dagger}_{\textnormal{\scriptsize
i,f}}(t)\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}+f(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},\bm{u}^{\dagger}),$ (39)

where $f(\bm{x}_{\textnormal{\scriptsize i}}^{\prime},\bm{u}^{\dagger})$ is an
unknown function, which has to be determined. Since the generating function
$F$ depends on the old momenta $\bm{x}_{\textnormal{\scriptsize i}}^{\prime}$
and the new coordinates $-j\bm{u}^{\dagger}$, it must verify the relations
(e.g., Landau & Lifshitz, 1969)

$\frac{\partial F}{\partial\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}=j{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger},\quad\frac{\partial
F}{\partial\bm{u}^{\dagger}}=j\bm{u},$ (40)

which imply

$j\bm{x}^{\prime\dagger}_{\textnormal{\scriptsize i,f}}(t)+\frac{\partial
f}{\partial\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}=j{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger},\quad-j\bm{x}_{\textnormal{\scriptsize
i,f}}^{\prime}(t)+\frac{\partial f}{\partial\bm{u}^{\dagger}}=j\bm{u}.$ (41)

With the choice $f(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},\bm{u}^{\dagger})=j\bm{u}^{\dagger}\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}$, these relations are both equivalent to Eq. (36).

The previous derivations imply that the variable transformation
$(\bm{x}_{\textnormal{\scriptsize i}},-j{\bm{x}_{\textnormal{\scriptsize
i}}}^{\dagger})\rightarrow(\bm{x}_{\textnormal{\scriptsize
i}}^{\prime},-j{\bm{x}_{\textnormal{\scriptsize
i}}^{\prime}}^{\dagger})\rightarrow(\bm{u},-j\bm{u}^{\dagger})$ is canonical,
with

$\displaystyle\bm{u}=\bm{O}_{\text{M}}^{T}\bm{x}_{\textnormal{\scriptsize
i}}-\bm{x}_{\textnormal{\scriptsize i,f}}^{\prime}(t),$ (42)
$\displaystyle\mathbb{H}_{\textnormal{LL}(\bm{u})}=\mathcal{H}_{\textnormal{LL}(\bm{x})}+\frac{\partial
F}{\partial t}=\bm{u}^{\dagger}\bm{D}_{\text{M}}\bm{u}.$ (43)

We call the complex variables $\bm{u}$ the _proper modes_ of the Poincaré’s
variables $\bm{x}_{\textnormal{\scriptsize i}}$. In the forced Laplace-
Lagrange dynamics, they simply rotate in the complex plane, with constant
angular frequencies given by $\bm{g_{\textnormal{\scriptsize LL}}}$. According
to Eq. (42), the corresponding Poincaré’s variables thus result in a
superposition of such independent harmonic oscillations and those arising from
the forcing of the outer planets. When higher-degree terms of the Hamiltonian
(25) are taken into account, the dynamics of the proper modes becomes coupled,
and their frequency spectrum is no more monochromatic. Physically, the proper
modes $\bm{u}$ correspond to the variables $(z^{*}_{k})_{k=1,4}$ defined in
(Laskar, 1990c). Nevertheless, it is important to emphasize that their
mathematical definitions do not coincide. Indeed, the transformation matrix
$\bm{O}_{\text{M}}$ do not include contributions at the second order in planet
masses, as that considered in (Laskar, 1990c) does. Moreover, the
transformation (42) is nonlinear in the forcing of the outer planets, in the
sense that harmonics of order999We define the order of a harmonic as the
1-norm of its wavevector, i.e. $|\bm{m}|=\sum_{\ell=1}^{7}|m_{\ell}|$.
$|\bm{m}|$ greater than one are present in the decomposition (15) (see Table
11) and thus in the forced solution (35).

Similar derivations allow to define the proper modes $\bm{v}$ of the
Poincaré’s variables $\bm{y}_{\textnormal{\scriptsize i}}$, by means of a
time-dependent canonical change of variables $(\bm{y}_{\textnormal{\scriptsize
i}},-j{\bm{y}_{\textnormal{\scriptsize
i}}}^{\dagger})\rightarrow(\bm{v},-j\bm{v}^{\dagger})$,

$\displaystyle\bm{v}=\bm{O}_{\text{N}}^{T}\bm{y}_{\textnormal{\scriptsize
i}}-\bm{y}_{\textnormal{\scriptsize i,f}}^{\prime}(t),$ (44)
$\displaystyle\mathbb{H}_{\textnormal{LL}(\bm{v})}=\bm{v}^{\dagger}\bm{D}_{\text{N}}\bm{v}.$
(45)

We have defined

$\bm{N}=\begin{pmatrix}\bm{N}_{\textnormal{\scriptsize
ii}}&\bm{N}_{\textnormal{\scriptsize io}}\\\ \bm{N}_{\textnormal{\scriptsize
io}}^{T}&\bm{N}_{\textnormal{\scriptsize
oo}}\end{pmatrix},\quad\bm{N}_{\textnormal{\scriptsize
ii}}=\bm{O}_{\text{N}}\bm{D}_{\text{N}}\bm{O}_{\text{N}}^{T},$ (46)

the columns of $\bm{O}_{\text{N}}$ being the eigenvectors of
$\bm{N}_{\textnormal{\scriptsize ii}}$, while the diagonal entries of
$\bm{D}_{\text{N}}=-\mathrm{diag}(\bm{s_{\textnormal{\scriptsize LL}}})$ are
the corresponding real eigenvalues. Using the quasi-periodic decomposition in
Eq. (15), the components of the forced solution
$\bm{y}_{\textnormal{\scriptsize i,f}}^{\prime}(t)$ are

$(\bm{y}_{\textnormal{\scriptsize
i,f}}^{\prime})_{k}=\sum_{\ell=5}^{8}\sum_{p=1}^{N_{\ell}}\frac{(\bm{N}_{\textnormal{\scriptsize
io}}^{\prime})_{k\ell}\tilde{y}_{\ell p}}{(\bm{s_{\textnormal{\scriptsize
LL}}})_{k}-\bm{n}_{\ell p}\cdot\bm{\omega}_{\textnormal{\scriptsize
o}}}\mathrm{E}^{j\bm{n}_{\ell p}\cdot\bm{\phi}}\quad k\in\\{1,2,3,4\\}.$ (47)

where $\bm{N}_{\textnormal{\scriptsize
io}}^{\prime}=\bm{O}_{\text{N}}^{T}\bm{N}_{\textnormal{\scriptsize io}}$.

Table 3: Frequency in arcsec yr-1 of the largest-amplitude Fourier harmonic for each proper mode of the Poincaré’s variables $(x_{k},y_{k})_{k=1,4}$ over the time interval [0, 20] Myr. For the solutions LaX13b and La90 the proper modes $(\bm{z}_{k}^{\star},\bm{\zeta}_{k}^{\star})_{k=1,4}$ defined in (Laskar, 1990c) are employed, while for the present solutions we use the proper modes $(\bm{u}_{k},\bm{v}_{k})_{k=1,4}$. | $g_{1}$ | $g_{2}$ | $g_{3}$ | $g_{4}$ | $s_{1}$ | $s_{2}$ | $s_{3}$ | $s_{4}$
---|---|---|---|---|---|---|---|---
LaX13b | 5.546 | 7.457 | 17.388 | 17.928 | -5.617 | -6.983 | -18.844 | -17.759
La90 | 5.531 | 7.460 | 17.379 | 17.927 | -5.620 | -7.011 | -18.846 | -17.758
$\mathcal{H}$ | 5.577 | 7.446 | 17.390 | 17.932 | -5.588 | -6.993 | -18.901 | -17.817
$\mathcal{H}_{4}$ | 5.631 | 7.457 | 17.371 | 17.923 | -5.602 | -7.071 | -18.922 | -17.817
$\mathcal{H}_{6}$ | 5.560 | 7.453 | 17.396 | 17.935 | -5.628 | -7.002 | -18.894 | -17.816
$\mathcal{H}_{8}$ | 5.576 | 7.446 | 17.390 | 17.932 | -5.589 | -6.993 | -18.900 | -17.816
$\mathcal{H}_{10}$ | 5.577 | 7.446 | 17.390 | 17.932 | -5.587 | -6.993 | -18.901 | -17.817

Once the proper modes $\bm{u},\bm{v}$ have been defined, it is very useful to
establish a corresponding truncated Hamiltonian in the form

$\mathbb{H}_{2n}=\sum_{p=1}^{n}\mathbb{H}_{(2p)},$ (48)

where
$\mathbb{H}_{2}=\mathbb{H}_{\textnormal{LL}(\bm{u})}+\mathbb{H}_{\textnormal{LL}(\bm{v})}=\bm{u}^{\dagger}\bm{D}_{\text{M}}\bm{u}+\bm{v}^{\dagger}\bm{D}_{\text{N}}\bm{v}$
is the LL Hamiltonian. To define higher-degree truncations, we first consider
the non-quadratic part $\delta\widehat{H}$ of the secular Hamiltonian (14),
truncated at a given total degree in Poincaré’s variables
$(x_{k},\bar{x}_{k},y_{k},\bar{y}_{k})_{k=1,8}$,

$\delta\widehat{H}_{2n}=\sum_{p=2}^{n}\widehat{H}_{(2p)}.$ (49)

Then, the proper modes $\bm{u},\bm{v}$ are injected in (49) by means of Eqs.
(42) and (44), and the Poincaré’s variables of the outer planets replaced by
their quasi-periodic decompositions given in Eq. (15). However, such
substitutions do not conserve the degree of the terms in the expansion, since
non-linear harmonics are present in Eqs. (15), (35) and (47), i.e. harmonics
of order $|\bm{m}|$ or $|\bm{n}|$ greater than one. Such harmonics are in
principle much smaller than the linear ones, and generate higher-degree terms
in the substitution process. To define the terms $\mathbb{H}_{(2p)},p\geq 2$
in a consistent way, we thus introduce a fictitious real variable $\epsilon$
to redefine the quasi-periodic decompositions (15) as

$\left.\begin{aligned}
&x_{k}^{(\epsilon)}(t)=\sum_{\ell=1}^{M_{k}}\epsilon^{|\bm{m}_{k\ell}|}\tilde{x}_{k\ell}\,\mathrm{E}^{j\bm{m}_{k\ell}\cdot\bm{\phi}(t)}\\\
&y_{k}^{(\epsilon)}(t)=\sum_{\ell=1}^{N_{k}}\epsilon^{|\bm{n}_{k\ell}|}\tilde{y}_{k\ell}\,\mathrm{E}^{j\bm{n}_{k\ell}\cdot\bm{\phi}(t)}\end{aligned}\;\right\\}\;k\in\\{5,6,7,8\\},$
(50)

meaning that each harmonic is treated as of degree $|\bm{m}_{k\ell}|$ or
$|\bm{n}_{k\ell}|$ for the purpose of the expansion101010Indeed, this would be
the case if such terms came from some analytical non-linear secular theory..
The forced solutions (35) and (47) are modified accordingly. The series
resulting from the substitution of variables in (49) can be thus truncated at
total degree $2n$ with respect to
$\\{(u_{k},\bar{u}_{k},v_{k},\bar{v}_{k})_{k=1,4},\epsilon\\}$. Finally,
addition of the quadratic terms $\mathbb{H}_{2}$ gives, by definition, the
truncated Hamiltonian (48). Truncation is effectively performed in TRIP in
parallel with variable substitutions, thanks to a dedicated monomial truncated
product (Gastineau & Laskar, 2011), which allows to minimize the computational
cost of the expansion. It is worthwhile to note that, at the end of the
expansion process, the Hamiltonian terms only depending on time can be
removed, since they are dynamically irrelevant. When numerical evaluations
involving the terms of the Hamiltonian have to be performed, the fictitious
variable $\epsilon$ is simply set to $1$. In Table 2 we show the total number
of monomials of variables
$\\{(u_{k},\bar{u}_{k},v_{k},\bar{v}_{k})_{k=1,4},(\mathrm{E}^{j\phi_{k}(t)})_{k=1,7}\\}$
in the Hamiltonian $\mathbb{H}_{2n}$ according the degree of truncation. We
also report the typical computation time of the series in TRIP, related to the
transformation
$(x_{k},\bar{x}_{k},y_{k},\bar{y}_{k})_{k=1,8}\rightarrow\\{(u_{k},\bar{u}_{k},v_{k},\bar{v}_{k})_{k=1,4},(\mathrm{E}^{j\phi_{k}(t)})_{k=1,7}\\}$.
One may note the dramatic increase of the number of terms due to such change
of variables. We point out that the precise number of monomials in
$\mathbb{H}_{2n}$ depends on the particular quasi-periodic decomposition (15)
we used in this work (see Appendix D).

(a) Mercury

(b) Venus

(c) Earth

(d) Mars

Figure 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 black solid line stands for the full (i.e. non-filtered)
direct integration LaX13b, while the red one is the solution of the present
forced model in Eq. (16). The secular integration La90 is represented by the
dashed blue line. The maximal relative root mean square of the present model
residuals is 2% for the eccentricities and 0.7% for the inclinations.

(a) Mercury

(b) Venus

(c) Earth

(d) Mars

Figure 2: Eccentricities (left side) and inclinations (right side) over the
invariant plane J2000 (see Appendix D.2) of the inner planets over 10 Myr in
the future. The black solid line stands for the full (i.e. non-filtered)
direct integration LaX13b, while the red one is the solution of the present
forced model in Eq. (16). The secular solution La90 practically coincides with
the direct integration on these plots and it is thus not shown.

####  4.1.3 Action-angle variables

Action-angle variables corresponding to the proper modes $\bm{u},\bm{v}$ are
introduced by standard canonical transformations
$(\bm{u},-j\bm{u}^{\dagger})\rightarrow(\bm{X},\bm{\chi})$ and
$(\bm{v},-j\bm{v}^{\dagger})\rightarrow(\bm{\Psi},\bm{\psi})$ defined as

$\displaystyle u_{k}=\sqrt{X_{k}}\,\mathrm{E}^{-j\chi_{k}},$ (51)
$\displaystyle v_{k}=\sqrt{\Psi_{k}}\,\mathrm{E}^{-j\psi_{k}},$

for $k\in\\{1,2,3,4\\}$. The Laplace-Lagrange Hamiltonian thus reads
$\mathbb{H}_{2}=-\bm{g}_{\text{LL}}\cdot\bm{X}-\bm{s}_{\text{LL}}\cdot\bm{\Psi}$
and it is trivially integrable. Such transformations allow to expand the
Hamiltonian $\mathbb{H}_{2n}$ as a Fourier series of the angle variables
$\bm{\chi}$, $\bm{\psi}$ and
$\bm{\phi}(t)=\bm{\omega}_{\textnormal{\scriptsize o}}t$,

$\mathbb{H}_{2n}=\sum_{\bm{k}\in\mathbb{Z}^{8}}\sum_{\bm{\ell}\in\mathbb{Z}^{7}}\widetilde{\mathbb{H}}_{2n}^{\bm{k},\bm{\ell}}(\bm{I})\,\mathrm{E}^{j\left(\bm{k}\cdot\bm{\theta}+\bm{\ell}\cdot\bm{\omega}_{\textnormal{\scriptsize
o}}t\right)},$ (52)

where $\widetilde{\mathbb{H}}_{2n}^{\bm{k},\bm{\ell}}$ are complex amplitudes
and we have employed a compact notation for the action-angle variables,

$\bm{I}=(\bm{X},\bm{\Psi}),\quad\bm{\theta}=(\bm{\chi},\bm{\psi}).$ (53)

Only a finite number of harmonics have non-zero amplitude in Eq. (52), since
we deal with a truncated Hamiltonian. We also recall that, the Hamiltonian
being a real function, one has
$\widetilde{\mathbb{H}}_{2n}^{-\bm{k},-\bm{\ell}}=\overline{\widetilde{\mathbb{H}}_{2n}^{\bm{k},\bm{\ell}}}$
for all $\bm{k},\bm{\ell}$. In Table 2 we show the number of harmonics in Eq.
(52) according the truncation degree of the Hamiltonian, i.e. the number of
different wavevectors $(\bm{k},\bm{\ell})$ up to a global minus sign.

It is worthwhile to mention that the explicit time dependency of the
Hamiltonian $\mathbb{H}_{2n}$ can be easily absorbed in a phase-space
extension, through the definition of action-angle variables for the trivial
degrees of freedom of the outer orbits. Indeed, one introduces
$(\bm{\Phi},\bm{\phi})$ such that the new Hamiltonian reads

$\mathbb{H}^{\star}_{2n}=\bm{\omega}_{\textnormal{\scriptsize
o}}\cdot\bm{\Phi}+\sum_{\bm{k}\in\mathbb{Z}^{8}}\sum_{\bm{\ell}\in\mathbb{Z}^{7}}\widetilde{\mathbb{H}}_{2n}^{\bm{k},\bm{\ell}}(\bm{I})\,\mathrm{E}^{j\left(\bm{k}\cdot\bm{\theta}+\bm{\ell}\cdot\bm{\phi}\right)},$
(54)

with $\bm{\Phi}=(\Phi_{k})_{k=1,7}$. The dynamics of the additional angles is
thus consistently given by
$\dot{\bm{\phi}}=\bm{\omega}_{\textnormal{\scriptsize o}}$, while that of the
actions $\bm{\Phi}$ is irrelevant. Such an autonomous formulation is useful in
the context of canonical perturbation theory.

##  5 Comparison with reference dynamical models

It is essential, when constructing averaged models, to compare the resulting
trajectories to those of the nominal Hamiltonian, or of a more comprehensive
dynamical model. This allows to validate both the underlying averaging
approximations and the choice of the initial conditions for the secular
variables (see Sect. 2.3 and Appendix D). In this Section, we compare the
orbital solution of the forced inner system (Eq. 16) with the direct numerical
integration of the full equations of motion LaX13b, which has been used to
predetermine the secular motion of the outer planets (Sect. 2.3 and Appendix
D). We shall also employ in the comparison the orbital solution of (Laskar,
1990c, La90 from now on), as it includes all the dynamical interactions taken
into account in the present work and constitutes the most precise secular
system to this day.

The comparison between the different models begins at short secular timescales
in Fig. 1, where we show the inner planet eccentricities and inclinations over
250 000 years in the future. The plots indicate that the present model
correctly reproduces the secular behaviour of the LaX13b solution. This is
particularly manifest in the eccentricity plots, where the forced model
constitutes the average of the short-time (orbital) oscillations of the direct
integration, as the La90 solution also does. In particular, the behaviour of
the solutions at the origin of time shows that the initial conditions for the
secular variables of this work have been correctly determined. Some small
periodic differences between the two secular models appear on some of the
plots, with the La90 solution generally being more precise in reproducing the
long-term average of the direct integration. However, such deviations are
practically irrelevant at these scales.

In Fig. 2 the comparison is extended to a longer time interval of 10 Myr, to
evaluate how long the present solution remains close to the direct
integration. The solution La90 practically coincides with LaX13b on these
plots (see Laskar et al., 2004) and it is thus not shown. The general
agreement shown over the first few million years by the curves in Fig. 2 is
still very satisfactory. For longer times, the divergence of the two models
becomes important. First of all, the orbital oscillations show sometime
different amplitudes. Nevertheless, as we are not interested in the
construction of a secular ephemeris (the solution La90 is clearly more adapted
to this end), this is not the most relevant point. When considering that the
inner solar system has a Lyapunov time of about 5 Myr (see Sect. 6), the
crucial aspect is rather the slow dephasing of the solutions appearing from
this comparison, as it translates a difference in the fundamental frequencies
of the motion. To ensure that the forced model faithfully reproduces the
resonant structure of the inner system, it is essential to quantify such
deviations. To this end, we perform a frequency analysis of the proper modes
of the Poincaré’s variables $(x_{k},y_{k})_{k=1,4}$ for the different orbital
solutions over the [0, 20] Myr time interval. For the solutions LaX13b and
La90 we employ the proper modes
$(\bm{z}_{k}^{\star},\bm{\zeta}_{k}^{\star})_{k=1,4}$ defined in (Laskar,
1990c), while for the present solution the proper modes
$(\bm{u}_{k},\bm{v}_{k})_{k=1,4}$ are used (the employ of
$(\bm{z}_{k}^{\star},\bm{\zeta}_{k}^{\star})_{k=1,4}$ gives the same results
at the numerical precision used in the comparison). Table 3 shows the dominant
frequency in the Fourier spectrum for each proper mode, according the notation
of (Laskar, 1990c). The average absolute difference between the LaX13b
solution and the present one amounts to a few hundredths of arc second per
year, with a maximum of 0.06$\arcsec\,{\textrm{yr}}^{-1}$ in the case of
$s_{3}$ and $s_{4}$. This is in agreement with the behaviour shown in Fig. 2
and the general expectation from the discussion on the model precision in
Sect. 2.1. It is important to note that, when considering combinations of the
fundamental frequencies (which determine the resonant structure), the
differences between the two models can be even smaller. Indeed, in the case of
the two leading resonances $(g_{1}-g_{5})-(s_{1}-s_{2})$ and
$2(g_{3}-g_{4})-(s_{3}-s_{4})$ described in (Laskar, 1990c), the deviations
amount to just 0.01$\arcsec\,{\textrm{yr}}^{-1}$ and
0.004$\arcsec\,{\textrm{yr}}^{-1}$, respectively.

Table 3 also reports the fundamental frequencies of the truncated model
$\mathcal{H}_{2n}$ (Eq. 25) for different total degrees of truncation. When
considering the combinations of frequencies $(g_{1}-g_{5})-(s_{1}-s_{2})$ and
$2(g_{3}-g_{4})-(s_{3}-s_{4})$, the models $\mathcal{H}_{4}$ and
$\mathcal{H}_{6}$ practically show the same deviations as the full forced
model $\mathcal{H}$ with respect to the solution LaX13b. This suggests that
the truncated model should already realistically reproduce the resonant
structure of the inner system at the lowest degrees. We also note that the
orbital solutions of the models $\mathcal{H}_{8}$ and $\mathcal{H}_{10}$
practically coincide with that of the full forced model $\mathcal{H}$ over the
first 20 Myr.

Figure 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 1 148 stable orbital solutions with very close initial conditions.

##  6 Finite-time maximum Lyapunov exponent

Once the adequacy of the present dynamical model has been established over a
few million years, one has to assess if it is able to correctly reproduce the
resonant structure of the inner solar system over longer times. Even if the
fundamental frequencies in Table 3 suggest this is actually the case, the
computation of the maximum Lyapunov exponent (MLE) remains an essential test,
as its value is related to the width of the leading resonant harmonics of the
Hamiltonian, which are the dynamical sources of stochasticity (Chirikov,
1979). At this point, it is important to realize that the non-null probability
of an unstable evolution of the inner planets (Laskar, 1994, 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); pertinent to the present case is the consideration of a finite-time MLE
(FT-MLE). As we do not have at our disposal an efficient numerical algorithm
to evaluate the variational equations of Gauss’s dynamics, we consider the FT-
MLE given by the two-particle algorithm proposed in (Benettin et al., 1976),
which has been used in (Laskar, 1989). This method computes the divergence of
close trajectories of a dynamical system, by considering the motion of
reference and shadow particles, initially separated by a tiny vector
$\bm{d}_{0}$ in the phase space. At time intervals $\tau$, a renormalization
procedure applied to the trajectory separation $\bm{d}(t)$ brings again the
shadow particle at a distance $||\bm{d}_{0}||$ from the reference one. The
resulting FT-MLE depends in a intricate way on the initial position of the
reference particle in the phase space, so that its asymptotic evolution is
chaotic (and independent of the choice of the vector $\bm{d}_{0}$ (Benettin et
al., 1976)). Therefore, its computation acquires full physical significance
only for an ensemble of trajectories.

Figure 4: Probability density function (PDF) of the Lyapunov time FT-MLE-1 at
5 Gyr. Kernel density estimation in dark-blue line and corresponding
cumulative distribution function (CDF) on the upper horizontal axis. Pointwise
confidence interval at 98% level from bootstrap as light-blue region.

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 1 148 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 the present
model). On the left vertical axis, we report the FT-MLE expressed as an
angular frequency in arcsec yr-1, while on the right one the corresponding
Lyapunov time, defined as FT-MLE-1, is given in million years. The dark-grey
region represents, at each renormalization time, the interquantile range of
the observed probability density function (PDF) of the FT-MLE, which encloses
by definition 50% of the solutions around the median (shown by the red line).
In the same manner, the light-grey region corresponds to the
$[1^{\textrm{th}},99^{\textrm{th}}]$-percentile range of the PDF, enclosing
the 98% of the probability. The two FT-MLEs reported in (Laskar, 1989) are
also shown with black lines.

As stated above, the initial conditions of each solution in the ensemble are
very close to the nominal ones. They are chosen by taking the relative
variation of each coordinate of the nominal phase-space vector as a normal
random variable with zero mean and a standard deviation of $10^{-9}$. The
initial conditions thus follow a multivariate Gaussian distribution centred at
the nominal phase-space position of the system, with a diagonal covariance
matrix. The width of such distribution is of the same order of magnitude of
that in (Laskar, 2008, Table 1). The initial position of each shadow particle,
around the corresponding reference one, is then chosen according the same kind
of Gaussian distribution, with a relative standard deviation of $10^{-8}$.
Indeed, such a value should be close to minimize the accumulation over long
timescales of numerical errors, due to the floating-point implementation of
the algorithm of (Benettin et al., 1976), when working in double precision
(Tancredi et al., 2001; Mei & Huang, 2018). In Fig. 3 the computation of the
FT-MLE for the nominal initial conditions is shown by dotted blue lines for a
set of 16 different initial tangent vectors $\bm{d}_{0}$. This manifestly
shows that the present FT-MLEs are asymptotically independent of $\bm{d}_{0}$,
the asymptotic regime being practically reached in a few hundred million
years. In the same manner, our computation has been tested to be
asymptotically independent of the renormalization time $\tau$, 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
$\bm{d}_{0}$, around essentially the same reference trajectory, i.e. the
nominal one. Therefore, this part of the plot can be discarded because of its
non-asymptotic character. We also note that, for $t<$ 100 Myr, the PDF of the
FT-MLE tends to shrinks with increasing time, as expected by the asymptotic
behaviour of the method employed (see the blue lines and two curves from
(Laskar, 1989)). For longer times, the PDF begins to broaden, because each
orbital solution acquires a macroscopically different FT-MLE, which
chaotically depends on its initial conditions. The time-asymptotic regions
shown in Fig. 3 thus represent the probabilistic knowledge of the FT-MLE of
the inner solar system, which arises from its chaotic behaviour and our
determination of the current planet positions and velocities.

The dark-blue line in Fig. 4 shows the kernel density estimation (Rosenblatt,
1956; Parzen, 1962) of the PDF of the Lyapunov time FT-MLE-1 at 5 Gyr, along
with the corresponding cumulative distribution function (CDF) on the upper
horizontal axis. We employed the standard Gaussian kernel and Silverman’s rule
of thumb for the selection of the optimal bandwidth (Silverman, 1986). We also
show, through the light-blue region, the pointwise confidence interval at 98%
level on the estimated PDF, obtained via nonparametric bootstrap, i.e. by
resampling the original data with replacement (Efron, 1979). One can
appreciate that the value of about 5 Myr typically reported in literature
(Laskar, 1989; Sussman & Wisdom, 1992) is right in the bulk of the present
PDF, even though the distribution rather peaks at 3.7 Myr and the average
value is 4.3 Myr. The value of 6.5 Myr found in (Rein & Tamayo, 2018) can also
be accounted for by our computation. The curves of (Laskar, 1989) suggest that
the asymptotic PDFs of the two secular models should largely overlap. Finally,
we note that numerical value of 1.1 Myr reported in (Batygin et al., 2015), as
the Lyapunov time of a simplified Mercury dynamics, does not agree with the
present findings. Figures 3 and 4 show that, even considering the diffusion of
the planet orbits over 5 billion years, such a high FT-MLE is practically
never reached by the present dynamical model, even though the latter includes
all the Hamiltonian harmonics taken into account in (Batygin et al., 2015).

##  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, the aim being to assess if the
present model is able to reproduce the rate of the high Mercury eccentricities
observed in (Laskar, 2008; Laskar & Gastineau, 2009). Following (Laskar,
1994), we characterize the unstable evolutions of the system by defining a
_secular collision_ as the intersection of the instantaneous Keplerian orbits
of a pair of planets. Numerically, such an event is easily detected by
tracking the relative positions of their mutual nodes, i.e. the intersections
of each orbit with the orbital plane of the other planet. When two orbits
crosses, a pair of mutual nodes exchanges their positions along the line of
nodes (i.e. the intersection of the two orbital planes). Such a computation is
inexpensive in the context of Gauss’s dynamics and can be performed at each
timestep of the integration scheme. When a secular collision occurs, the
present averaged model is no longer a good long-term approximation of the
original N-body dynamics, in which the involved planets can experience a close
encounter, possibly leading to a physical collision or to the escape of one of
two planets from the solar system (Laskar, 1994). We thus denote by $\tau_{c}$
the time of the _first_ secular collision for a given orbital solution. The
statistics of the first-collision time $\tau_{c}$ shall provide an estimate of
the rate of catastrophic events in the real inner system. Moreover, by
stopping the numerical integration at the first secular collision, we will not
be involved in extrapolating the present secular dynamics out of its validity
range, and the use of a fixed-timestep integration scheme (see Sect. 3) turns
out to be sufficient end even suitable. We note that a secular collision
between a planet and the Sun can be simply defined to occur whenever its
pericenter distance $a(1-e)$ is smaller than the Sun radius.

$e_{\textrm{max}}$ | LG09 | La08 | This work
---|---|---|---
0.35 | (46.9, 51.5) | (37.5, 48.0) | (48.5, 50.8)
0.4 | (23.6, 27.6) | (20.0, 29.1) | (20.3, 22.1)
0.5 | (3.15, 4.97) | (1.90, 5.84) | (1.74, 2.38)
0.6 | (0.60, 1.53) | (0.39, 2.80) | (0.38, 0.71)
0.7 | (0.57, 1.48) | (0.39, 2.80) | (0.32, 0.62)
0.8 | (0.54, 1.43) | (0.28, 2.50) | (0.28, 0.57)
0.9 | (0.51, 1.38) | (0.03, 1.50) | (0.08, 0.26)
Table 4: Confidence interval at 98% level of the probability $P(\sup_{t\leq
5\,\textrm{Gyr}}e_{1}(t)>e_{\textrm{max}})$ in percentage, where $e_{1}$ is
the eccentricity of Mercury. LG09 represents the 2 501 orbital solutions in
(Laskar & Gastineau, 2009), while La08 stands for the 478 integrations in
(Laskar, 2008). The statistics of this work derives from an ensemble of 10 560
orbital solutions.

We perform an ensemble of 10 560 numerical integrations of the forced inner
system over 5 Gyr. The initial conditions are chosen according the same
multivariate Gaussian distribution centred at the nominal values as in Sect.
6. The relative standard deviation of each phase-space coordinate is
$10^{-9}$, of the same order of magnitude of that in (Laskar, 2008, Table 1).
Table 4 shows the observed probability of having, over the entire time span, a
maximum Mercury eccentricity greater than a given value. The statistics of the
present dynamical model is compared to those resulting from (Laskar, 2008,
Table 3) and (Laskar & Gastineau, 2009, see Table 6 in Appendix C of the
present paper). The statistical bounds represent the Wilson (1927) score
interval corresponding to a 98% confidence level. At the present statistical
precision111111As there is no reason for the true probabilities of the three
different dynamical models to be exactly the same, for a very large number of
orbital solutions their statistical intervals are not likely to overlap., the
probabilities of the high Mercury eccentricities in the forced inner system
are fully compatible with those arising in the secular model of (Laskar,
1990c). They also represent a very good estimate of the probabilities
resulting from N-body dynamics, even if somewhat lower. This is in agreement
with the general expectations from an averaged model, which has less degrees
of freedom (Laskar, 1994, 2008). When comparing the midpoints of the
statistical intervals, Table 4 could also suggest that the true probabilities
of the present model are slightly lower than those in (Laskar, 2008). If true,
such a fact would agree with the expectation of higher regularity of the
forced model, as the degrees of freedom dominated by the outer planets (their
fundamental frequencies in particular) are frozen into a quasi-periodic time
dependency. Moreover, the Hamiltonian terms at second order in planet masses
could also contribute to somewhat larger chaotic excursions of the inner
orbits.

###  7.1 Secular collisions

Among the 10 560 orbital solutions of our ensemble, we found 42 secular
collisions, corresponding to a 98% confidence interval $(0.28\%,0.57\%)$ for
the probability $P(\tau_{c}\leq 5\,\textrm{Gyr})$. Out of the 2501 integrated
solutions of (Laskar & Gastineau, 2009), there result 6 collisions between
Mercury and the Sun, 9 collisions between Mercury and Venus, and 1 collision
between Mars and the Sun (see Appendix C). These 16 events corresponds to a
98% confidence interval $(0.36\%,1.13\%)$ of probability. Our definition of
secular collision thus provides a very good estimate of the probability of a
physical collision in the non-averaged inner solar system over the next 5 Gyr,
probably somewhat lower in agreement with the previous discussion.

All the secular collisions in our ensemble of orbital solutions involve the
Mercury-Venus couple, with a maximum Mercury eccentricity ranging from 0.797
to 0.966. In particular, in 41 out of 42 collisional solutions it is larger
than 0.8, while in 15 solutions it exceeds 0.9. This implies that the
statistics corresponding to the cases $e_{\textrm{max}}=0.8$ and $0.9$ in
Table 4 is affected by the choice of stopping the numerical integration when
the first secular collisions is detected, independently of the fact that this
could not physically correspond to any unstable further evolution. We
emphasize that, as a consequence, the statistics of the secular collisions
other than the Mercury-Venus one cannot be correctly reproduced by the present
numerical computation, since such events, as observed in (Laskar & Gastineau,
2009), always require a very large Mercury eccentricity at some previous time.

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

###  7.2 Secular resonance $g_{1}-g_{5}$

The destabilizing role of the secular resonance $g_{1}-g_{5}$, able to drive
Mercury orbit to very high eccentricities, has been first established
numerically (Laskar, 2008; Batygin & Laughlin, 2008) and then confirmed by an
analytical study (Boué et al., 2012). We conclude this Section by showing that
the $g_{1}-g_{5}$ resonance is indeed reached by the collisional orbital
solutions of our ensemble. While $g_{5}\approx$
4.257$\arcsec\,{\textrm{yr}}^{-1}$ is constant in the present dynamical model,
the frequency $g_{1}$, which dominates the spectrum of the Poincaré’s variable
$x_{1}$ and thus the precession of Mercury perihelion, is numerically computed
in the following way. Given an orbital solution, we sample the corresponding
Poincaré’s variables $\bm{x}$ with a timestep $\Delta t$ = 1 000 yr, i.e. we
consider the time series $\bm{x}(t_{n})$ with $t_{n}=n\Delta t$ and
$n\in\mathbb{N}_{0}$. The proper mode $u_{1}(t_{n})$ is thus computed through
Eq. (42) and the corresponding angle variable $\chi_{1}(t_{n})$, following the
definition in Eq. (51), retrieved as an unwrapped phase, i.e. as a continuous
function of time. The instantaneous (angular) frequency $\omega_{1}$ of the
proper mode $u_{1}$ is given by (Cohen, 1995)

$\omega_{1}=-\frac{d\chi_{1}}{dt}.$ (55)

We could analytically compute the derivative of Eq. (42) to obtain the time
series $\omega_{1}(t_{n})$. In practice, we employ numerical differentiation
via the Lagrange five-point formula (Olver et al., 2020, Chapter 3), as very
high numerical precision is unnecessary here. As it is, the frequency
$\omega_{1}(t_{n})$ presents short-time oscillations which hide its long-term
evolution. To suppress such high-frequency fluctuations, we thus define the
time series $g_{1}(t_{n})$ as the output of the low-pass Kolmogorov-Zurbenko
(KZ) filter (Yang & Zurbenko, 2010, and references therein) applied to
$\omega_{1}(t_{n})$,

$g_{1}(t_{n})=\textrm{KZ}[\omega_{1}(t_{n})].$ (56)

The KZ filter is defined as an iteration of the common moving average and has
been comprehensively characterized in Appendix B. For the present application,
we use three iterations of the moving average and a cutoff frequency of 1
$\textrm{Myr}^{-1}$ (see Appendix B). This particular choice is motivated by
the typical duration of the libration episodes of the $g_{1}-g_{5}$ resonance,
which only last a few million years (Laskar, 2008; Batygin & Laughlin, 2008).

In the top panel of Fig. 5 we present the time evolution of the fundamental
frequency $g_{1}$ along the 42 collisional integrations. We isolate in the
bottom panel six solutions to show the typical evolution of a single curve
more clearly. As previously discussed, the curves stop at the first secular
collision. The horizontal line represents the constant fundamental frequency
$g_{5}$ and thus the location of the $g_{1}-g_{5}$ resonance. Within the first
2 Gyr, the chaotic diffusion of $g_{1}$ along the solutions of Fig. 5 is
limited to no more than 0.25$\arcsec\,{\textrm{yr}}^{-1}$ with respect to its
initial value. Nevertheless, when lower values of about
5$\arcsec\,{\textrm{yr}}^{-1}$ are eventually reached, the dynamics of $g_{1}$
starts to be characterized by large random variations. This seems to reproduce
the large chaotic zone related to the resonance $g_{1}-g_{5}$, with a half-
width of more than 1$\arcsec\,{\textrm{yr}}^{-1}$, which has been reported in
(Laskar, 2008, Fig. 6(a2))121212Although no analytical computations of the
$g_{1}-g_{5}$ half-width have been reported in literature, such a large value
would not be completely unexpected for a linear secular resonance.. Along all
the curves, the fluctuations eventually lead to a crossing of the resonance
(corresponding locally to a libration of the related argument), shortly before
the first secular collision.

##  8 A numerical experiment over 100 Gyr

In this last Section, we present the results of a numerical experiment
consisting in integrating an ensemble of orbital solutions of the forced inner
system over 100 Gyr in the future. Through such simulations, we aim to
characterize the kind of stochastic process which arises from the Hamiltonian
(16) and drives the destabilization of the inner system, in a regime where the
highly excited orbits do not represent rare events anymore, as they do in
Sect. 7. Obviously, this kind of simulation does not correspond in anyway to
the real future evolution of the orbits of the inner planets over much more
than 5 Gyr. Indeed, it is well known that over the next 7 Gyr, the Sun will
experience a significant mass loss, down to 0.54 $M_{\odot}$, as it leaves the
main sequence (Sackmann et al., 1993). As a consequence, the semi-major axes
of the planet will adiabatically inflate by a factor of 1.85. Even more
critically for the fate of the inner system, it is likely that the innermost
planets Mercury, Venus and the Earth, will be engulfed by the Sun as its
radius expands along the red giant branch, marking the end of their existence
(Rybicki & Denis, 2001; Schröder & Smith, 2008). Therefore, this Section
intends to investigate the very long-term stochastic behaviour of the
dynamical system (16), without inferring any conclusion on the real evolution
of the inner system over such timescales.

We numerically integrate an ensemble of 1 042 orbital solutions of the forced
inner system over 100 Gyr. The initial conditions are chosen according the
same Gaussian distribution as in Sect. 7, but with different realizations, to
obtain different evolutions over the first 5 Gyr. As in the previous Section,
the numerical integration is stopped at the first secular collision. In Fig. 6
we show the probability density function of the time $\tau_{c}$ of the first
secular collision. In the context of stochastic processes, this represents the
hitting time (or first hit time) corresponding to the subspace of the phase
space realizing a secular collisions between planets. The blue curve in the
lower panel of Fig. 6 represents the PDF estimated via the kernel density
estimation (KDE) method (Rosenblatt, 1956; Parzen, 1962). We employed a
standard Gaussian kernel and Silverman’s rule of thumb for the selection of
the optimal bandwidth (Silverman, 1986). The PDF is normalized in such a way
that its integral over 100 Gyr is equal to the overall percentage of secular
collisions obtained in the experiment, i.e. $890/1042\approx 85.4\%$. The blue
region in the lower panel of Fig. 6 represents the pointwise confidence
interval for the PDF at 98% level, obtained via nonparametric bootstrap, i.e.
by resampling the original data with replacement (Efron, 1979). This region
thus corresponds pointwisely to the
[$1^{\mathrm{th}}$,$99^{\mathrm{th}}$]-percentile range of the bootstrapped
estimated PDF. We emphasize that, for simplicity, we do not correct for the
bias arising in the kernel density estimation, so that it constitutes,
strictly speaking, a confidence region for the expectation of the KDE and not
for the true subjacent PDF (Chen, 2017). When referred to the true PDF, such a
confidence region typically represents an undercoverage at the points where
the function bends, i.e. where its second derivative is large. Even if we
could consider simple corrections for the KDE bias (Chen, 2017), this would
not be completely satisfactory given the very small number of events which we
deal with over the first few Gyr. Anyway, such corrections on the confidence
intervals turn out to be sufficiently small to not affect at all the kind of
conclusions we infer in this Section. To give a further idea of the
statistical variance, in Fig. 6 we also show in red colour the PDF of
$\tau_{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.

Figure 6: Kernel density estimation of the PDF of the time $\tau_{c}$ of the
first secular collision, from 1 042 integrations of the forced inner system
over 100 Gyr (blue line in the lower panel). The estimation of the cumulative
distribution function (CDF) is shown as a blue line in the upper panel. The
blue regions represent the pointwise confidence intervals at 98% level from
bootstrap. In red colour the PDF estimation over 5 Gyr from the 10 560 orbital
solutions of Sect. 7. The dashed black lines stand for the analytical PDF (57)
with $T_{1}$ = 27.6 Gyr and $\alpha$ = 0.9 and the corresponding CDF.

###  8.1 Marginal stability of the inner solar system

The PDF of the time $\tau_{c}$ of the first secular collision peaks at 17.6
Gyr, with a 98% confidence interval (15.4, 20.1) Gyr, while its median is 40.8
Gyr, with a 98% confidence interval (38.2, 43.5) Gyr. The percentage of
collisional evolutions over 100 Gyr has a 98% confidence interval (82.8%,
87.9%). These findings clarify the state of matastability which characterizes
the inner solar system, according to the definition of (Laskar, 1996). Over a
timescale comparable to its age, i.e. over the next 5 Gyr, the inner system is
statistically very stable, with a probability of only about 0.5% of an
unstable evolution (see Sect. 7.1). Nevertheless, here we find that such
probability increases very fast with time, so that the system rapidly becomes
statistically unstable over longer timescales, meaning that unstable
evolutions do not represent rare events anymore. In particular, the
probability of an instability is about 20% at 20 Gyr already, and 50% at 40
Gyr. In the conjecture of (Laskar, 1996), this picture should characterize any
secularly evolving planetary system: at each moment of its formation history,
a planetary system should be in a state of marginal stability, that is
practically stable over a timescale comparable to its age, while strong
instabilities arise very fast over longer times, resulting in more stable
configurations of the surviving planets.

###  8.2 Destabilization of the inner solar system as a stochastic process

Recently, in the framework of the simplified model introduced in (Batygin et
al., 2015), Woillez & Bouchet (2020) applied the theory of the white noise
limit for slow-fast dynamical systems (e.g., Gardiner, 1985), to describe the
long-term dynamics of Mercury through an effective stochastic process. They
show that the instability timescale of Mercury orbit is reproduced by the
hitting time of a one-dimensional Wiener process with a reflecting barrier.
The corresponding PDF peaks at 1.3 Gyr (Woillez & Bouchet, 2020, Figure 4), in
agreement with the analysis in (Batygin et al., 2015). Such a timescale is an
order of magnitude smaller than that in Fig. 6, and the simplified dynamics of
Mercury predicts a large probability of instability over the solar system age,
in contrast with the findings of realistic models.

In this Section, we reconsider the stochastic process discussed in (Woillez &
Bouchet, 2020) to characterize the overall structure of the PDF in Fig. 6. To
fix ideas, it could describe in a crude way the long-term wanderings of the
fundamental frequency $g_{1}$ towards the destabilizing secular resonance
$g_{1}-g_{5}$ (see Sect. 7.2). Indeed, numerically, the PDF of $g_{1}$ turns
out to be asymmetric, with a longer tail diffusing towards low values, while
its upper end is practically fixed at about
6$\arcsec\,{\textrm{yr}}^{-1}$(Hoang et al., 2021, see also Fig. 5).
Therefore, let suppose that $g_{1}$ performs a random walk described by a
Wiener process starting at $g_{1,0}=g_{1}(t=0)$, and consider an upper
reflecting barrier at $g_{1,\textrm{max}}>g_{1,0}$. We assume that a secular
collisions rapidly occurs whenever the secular resonance $g_{1}-g_{5}$ is
reached, as shown in Fig. 5; this is taken into account by considering an
absorbing barrier at $g_{1}=g_{5}$. The PDF $\rho(\tau)$ of the hitting time
$\tau=\inf_{t}\\{g_{1}(t)\leq g_{5}\\}$ is thus given by (e.g., Schwarz, 1992;
Woillez & Bouchet, 2020)

$\rho(\tau)=\frac{\pi}{2T_{1}}\sum_{n=0}^{\infty}\left(n+\frac{1}{2}\right)\sin\left[\pi\left(n+\frac{1}{2}\right)\alpha\right]\mathrm{E}^{-\pi^{2}\left(n+\frac{1}{2}\right)^{2}\tau/4T_{1}},$
(57)

where $\alpha=(g_{1,0}-g_{5})/(g_{1,\textrm{max}}-g_{5})$ and
$T_{1}=(g_{1,\textrm{max}}-g_{5})^{2}/4D$, with $D$ being the diffusion
coefficient of the Wiener process. After a least-squares search, we plot in
Fig. 6 as a dashed black line the curve (57) with parameters $T_{1}$ = 27.6
Gyr and $\alpha$ = 0.9. Even if the real random walk performed by the
frequency $g_{1}$, and the destabilization mechanism in general, are likely to
be much more complex, a Wiener process with a reflecting barrier is able to
reproduce the overall behaviour of the observed PDF, at the level of the
present statistical precision and apart from a certain excess of events in the
tail of the distribution131313We stress that such a result is independent of
the hypothesis of which precise dynamical quantity (if any) actually undergoes
the stochastic process considered here.. On the one hand, this shows that at
short times the PDF behaves as

$\rho(\tau)=\left(\frac{T_{0}}{\pi\tau^{3}}\right)^{\frac{1}{2}}\mathrm{E}^{-\frac{T_{0}}{\tau}},\quad\textrm{for
}\tau\ll T_{1},$ (58)

where $T_{0}=\alpha^{2}T_{1}$. In this regime, the PDF coincides with that of
the hitting time of a standard Wiener process (i.e., without reflecting
barrier). The very fast decaying of the density probability for
$\tau\rightarrow 0$ in Eq. (58) characterizes the rare destabilizations of
Mercury orbit over the first few billions of years. In particular, the
instanton phenomenology described in (Woillez & Bouchet, 2020) could indeed
apply to the real solar system (such a hint still needs to be tested through
integration of realistic models, as the present one.). On the other hand, Eq.
(57) shows that at large times the decay of the observed PDF is at least
exponentially fast, with a characteristic time equal to $4T_{1}$. This
behaviour results from the reflecting barrier taken into account in the
stochastic process; it would not be reproduced by a standard Wiener process,
as the corresponding PDF has a polynomial heavy tail proportional to
$\tau^{-3/2}$, as shown by Eq. (58). In particular, the average of the hitting
time as given by Eq. (57) is finite, differently from Eq. (58).

We note that (Woillez & Bouchet, 2020) reported $T_{0}$ = 1.56 Gyr, a value
which is more than ten times smaller than that in Fig. 6, i.e. $T_{0}$ = 22.4
Gyr, in agreement with the previous discussion about the validity of the
simplified dynamics of Mercury in (Batygin et al., 2015). Here we emphasize
that the incompatibility with realistic models seems to be related to the
limitations of some of the simplifying assumptions141414Freezing the
fundamental frequencies of the inner planet orbits other than $g_{1}$ and
$s_{1}$ is probably one of the sources of these discrepancies.. Indeed, as
discussed in Sect. 7, averaged models generally tends to produce slower
instabilities with respect to the full dynamics, because the number of degrees
of freedom is smaller. Nevertheless, the model of (Batygin et al., 2015)
produces instabilities which are ten times faster than in the present
dynamics, even though the latter includes all the dynamical interactions
considered in the simplified one. Such a lack of continuity in the
destabilization time with respect to the complexity of the dynamical model
seems rather severe.

We conclude this Section by noting that a lower bound for the typical
instability time of Mercury orbit can be estimated from the maximum Lyapunov
exponent (Fig. 3), if one assumes that no dynamical constraints prevent
$g_{1}$ to diffuse at a rate determined by the leading secular resonances. In
particular, the resonance $(g_{1}-g_{5})-(s_{1}-s_{2})$, with a libration
frequency of 0.12$\arcsec$ yr-1, is among the leading ones affecting the
frequency $g_{1}$ (Laskar, 1990c). We may estimate an upper bound on the
diffusion coefficient $D$ as

$D_{\textrm{max}}=\frac{\left(2\pi\,\textrm{MLE}\right)^{2}}{\textrm{MLE}^{-1}}\approx
0.77\left(\textrm{arcsec yr}^{-1}\right)^{2}\textrm{Gyr}^{-1},$ (59)

with $2\pi\,\textrm{MLE}$ = 0.1$\arcsec\,{\textrm{yr}}^{-1}$ (this is also the
lower limit for the MLE in Fig. 3). A lower bound on the typical
destabilization time of Mercury orbit is thus given by

$T_{0,\textrm{min}}=\frac{(g_{1,0}-g_{5})^{2}}{4D_{\textrm{max}}}\approx
0.56\,\textrm{Gyr},$ (60)

with $g_{1,0}$ = 5.577$\arcsec\,{\textrm{yr}}^{-1}$ (Table 3) and $g_{5}$ =
4.257$\arcsec\,{\textrm{yr}}^{-1}$ . Eq. (60) gives essentially the same
diffusion timescale as in (Batygin et al., 2015; Woillez & Bouchet, 2020).
This is already two order of magnitude greater than the Lyapunov time of the
system, but still insufficient to explain the findings of Fig. 6. This means
that some dynamical constraints effectively makes the diffusion towards the
$g_{1}-g_{5}$ resonance slower than in the estimates (59), (60). Certain
secondary resonant harmonics, in particular, must play a determinant role in
the destabilization of the inner solar system.

##  9 Conclusions

This work introduces a new secular model for the dynamics of the inner planets
of the solar system. It exploits the practical constancy of the fundamental
precession frequencies of the outer planet orbits over timescales of billion
years, the smallness of the inner planet masses and the absence of strong
mean-motion resonances in the inner system. The range of validity of the
resulting dynamics extends to very high eccentricities and inclinations, up to
a secular planetary collision, i.e. the geometric intersection of the
instantaneous Keplerian orbits of two planets. The new model can be studied
analytically by truncating the Hamiltonian at a given total degree in planet
eccentricities and inclinations, with the aid of a computer algebra system
like TRIP. It can also be integrated numerically in a very efficient way
through Gauss’s method. The orbital solution matches the predictions of a
comprehensive model of the solar system at a very satisfactory level over
timescales shorter than or comparable to the Lyapunov time. The new model also
correctly reproduces the maximal Lyapunov exponent of the inner system and the
statistics of the highly eccentric Mercury orbits over the next 5 billion
years. Moreover, the destabilizing role of the secular resonance $g_{1}-g_{5}$
clearly stands out. We performed a numerical experiment consisting of a
thousand orbital solutions of the inner solar system over one hundred billion
years, to explore a regime in which unstable orbits are statistically common.
We first pointed out the fast growth of the rate of orbit instabilities in the
framework of planetary system formation through successive metastable states.
We then showed that the PDF of the time of the first secular collision is
reasonably well reproduced by a Wiener process with a reflecting barrier,
which could be performed, for example, by the fundamental precession frequency
$g_{1}$. Given the robustness of the statistical predictions of the present
dynamical model over 5 Gyr, the main properties of this PDF, that is the
characteristic time of its peak and the behaviours at short and long times,
are likely to represent what would arise from the full dynamics of the
planets. We finally argued that a dynamical mechanics is needed to explain the
rarity of the large excursions of $g_{1}$ up to the $g_{1}-g_{5}$ secular
resonance within the next 5 Gyr.

We emphasize that the new dynamical model can be straightforwardly implemented
once truncated at degree 4 in planet eccentricities and inclinations, by using
the corresponding expression of the two-body secular Hamiltonian reported in
(Laskar & Robutel, 1995, Appendix) and the quasi-periodic secular solution for
the outer planets given in Appendix D.

###### 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. This project
has received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Advanced Grant
AstroGeo-885250) and from the French Agence Nationale de la Recherche (ANR)
(Grant AstroMeso ANR-19-CE31-0002-01). This work was granted access to the HPC
resources of MesoPSL financed by the Region Île-de-France and the project
Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements
d’Avenir supervised by the Agence Nationale pour la Recherche.

## References

  * Batygin & Laughlin (2008) Batygin, K. & Laughlin, G. 2008, ApJ, 683, 1207
  * Batygin et al. (2015) Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120
  * Benettin et al. (1976) Benettin, G., Galgani, L., & Strelcyn, J.-M. 1976, Phys. Rev. A, 14, 2338
  * Boué et al. (2012) Boué, G., Laskar, J., & Farago, F. 2012, A&A, 548, A43
  * Bour (1855) Bour, E. 1855, PhD thesis, Faculté des Sciences de Paris
  * Carpino et al. (1987) Carpino, M., Milani, A., & Nobili, A. M. 1987, A&A, 181, 182
  * Chen (2017) Chen, Y.-C. 2017, arXiv e-prints, arXiv:1704.03924
  * Chirikov (1979) Chirikov, B. V. 1979, Phys. Rep, 52, 263
  * Cohen (1995) Cohen, L. 1995, Time-Frequency Analysis (Prentice Hall PTR, Upper Saddle River, New Jersey 07458)
  * Deprit (1969) Deprit, A. 1969, Celestial Mechanics, 1, 12
  * Efron (1979) Efron, B. 1979, The Annals of Statistics, 7, 1
  * Fienga et al. (2014) Fienga, A., Manche, H., Laskar, J., Gastineau, M., & Verma, A. 2014, arXiv e-prints, arXiv:1405.0484
  * Fouvry et al. (2020) Fouvry, J.-B., Dehnen, W., Tremaine, S., & Bar-Or, B. 2020, arXiv e-prints, arXiv:2011.01673
  * Fukushima (2015) Fukushima, T. 2015, Journal of Computational and Applied Mathematics, 282, 71
  * Gardiner (1985) Gardiner, C. W. 1985, Stochastic Methods (Springer-Verlag, Berlin, Heidelberg, New York, Tokyo)
  * Gastineau & Laskar (2011) Gastineau, M. & Laskar, J. 2011, ACM Communications in Computer Algebra
  * Gastineau & Laskar (2020) Gastineau, M. & Laskar, J. 2020, TRIP 1.4.94, TRIP Reference manual, IMCCE, Paris Observatory, https://www.imcce.fr/trip/
  * Gauss (1818) Gauss, C. 1818, Werke, 3, 331
  * Hill (1882) Hill, G. W. 1882, Astronomical Papers of the American Ephemeris, 1, 315
  * Hoang et al. (2021) Hoang, N., Mogavero, F., & Laskar, J. 2021, ”Accepted for publication in A&A”
  * Hori (1966) Hori, G. 1966, PASJ, 18, 287
  * Landau & Lifshitz (1969) Landau, L. D. & Lifshitz, E. M. 1969, Mechanics (Pergamon Press)
  * Laplace (1785) Laplace, P. S. 1785, Mémoires de l’Académie royale des sciences de Paris, 11, 95
  * Laskar (1984) Laskar, J. 1984, PhD thesis, Observatoire de Paris
  * Laskar (1985) Laskar, J. 1985, A&A, 144, 133
  * Laskar (1988) Laskar, J. 1988, A&A, 198, 341
  * Laskar (1989) Laskar, J. 1989, Nature, 338, 237
  * Laskar (1990a) Laskar, J. 1990a, in Modern Methods in Celestial Mechanics, ed. D. Benest & C. Froeschle, 89
  * Laskar (1990b) Laskar, J. 1990b, in Modern Methods in Celestial Mechanics, ed. D. Benest & C. Froeschle (Editions Frontières, Gif -Sur-Yvette), 63–87
  * Laskar (1990c) Laskar, J. 1990c, Icarus, 88, 266
  * Laskar (1991) Laskar, J. 1991, in NATO Advanced Study Institute (ASI) Series B, Vol. 272, Predictability, Stability, and Chaos in N-Body Dynamical Systems, 93–114
  * Laskar (1992) Laskar, J. 1992, in IAU Symposium, Vol. 152, Chaos, Resonance, and Collective Dynamical Phenomena in the Solar System, ed. S. Ferraz-Mello, 1
  * Laskar (1993) Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257
  * Laskar (1994) Laskar, J. 1994, A&A, 287, L9
  * Laskar (1996) Laskar, J. 1996, Celestial Mechanics and Dynamical Astronomy, 64, 115
  * Laskar (1999) Laskar, J. 1999, Philosophical Transactions of the Royal Society of London Series A, 357, 1735
  * Laskar (2005) Laskar, J. 2005, in Hamiltonian Systems and Fourier Analysis: New Prospects For Gravitational Dynamics, ed. D. Benest, C. Froeschlé, & E. Lega (Cambridge Scientific Publishers Ltd), 93–114, arXiv: math/0305364
  * Laskar (2008) Laskar, J. 2008, Icarus, 196, 1
  * Laskar & Boué (2010) Laskar, J. & Boué, G. 2010, A&A, 522, A60
  * Laskar et al. (2011) Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011, A&A, 532, A89
  * Laskar et al. (1992) Laskar, J., Froeschlé, C., & Celletti, A. 1992, Physica D Nonlinear Phenomena, 56, 253
  * Laskar & Gastineau (2009) Laskar, J. & Gastineau, M. 2009, Nature, 459, 817
  * Laskar & Robutel (1995) Laskar, J. & Robutel, P. 1995, Celestial Mechanics and Dynamical Astronomy, 62, 193
  * Laskar & Robutel (2001) Laskar, J. & Robutel, P. 2001, Celestial Mechanics and Dynamical Astronomy, 80, 39
  * Laskar et al. (2004) Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261
  * Laskar & Simon (1988) Laskar, J. & Simon, J. L. 1988, Celestial Mechanics, 43, 37
  * Lithwick & Wu (2011) Lithwick, Y. & Wu, Y. 2011, ApJ, 739, 31
  * Ma et al. (1998) Ma, C., Arias, E. F., Eubanks, T. M., et al. 1998, AJ, 116, 516
  * Mei & Huang (2018) Mei, L. & Huang, L. 2018, Computer Physics Communications, 224, 108
  * Mogavero (2017) Mogavero, F. 2017, A&A, 606, A79
  * Mogavero (2021) Mogavero, F. 2021, ”In preparation”
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics: aspects of solar system dynamics (Taylor & Francis)
  * Musen (1970) Musen, P. 1970, Celestial Mechanics, 2, 41
  * Nobili et al. (1989) Nobili, A. M., Milani, A., & Carpino, M. 1989, A&A, 210, 313
  * Olver et al. (2020) Olver, F. W. J., Olde Daalhuis, A. B., Lozier, D. W., et al. 2020, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.28 of 2020-09-15
  * Oseledec (1968) Oseledec, V. I. 1968, Trans.Moscow Math.Soc., 19, 197
  * Parzen (1962) Parzen, E. 1962, Ann. Math. Statist., 33, 1065
  * Poincaré (1896) Poincaré, H. 1896, Comptes rendus hebdomadaires de l’Académie des sciences de Paris, 123, 1031
  * Rein & Tamayo (2018) Rein, H. & Tamayo, D. 2018, MNRAS, 473, 3351
  * Rosenblatt (1956) Rosenblatt, M. 1956, Ann. Math. Statist., 27, 832
  * Rybicki & Denis (2001) Rybicki, K. R. & Denis, C. 2001, Icarus, 151, 130
  * Sackmann et al. (1993) Sackmann, I. J., Boothroyd, A. I., & Kraemer, K. E. 1993, ApJ, 418, 457
  * Saha & Tremaine (1992) Saha, P. & Tremaine, S. 1992, AJ, 104, 1633
  * Schröder & Smith (2008) Schröder, K. P. & Smith, R. C. 2008, MNRAS, 386, 155
  * Schwarz (1992) Schwarz, W. 1992, Journal of Applied Probability, 29, 597
  * Silverman (1986) Silverman, B. W. 1986, Density estimation for statistics and data analysis (London: Chapman & Hall/CRC)
  * Sussman & Wisdom (1992) Sussman, G. J. & Wisdom, J. 1992, Science, 257, 56
  * Tancredi et al. (2001) Tancredi, G., Sánchez, A., & Roig, F. 2001, AJ, 121, 1171
  * Touma et al. (2009) Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085
  * Verma et al. (2014) Verma, A. K., Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2014, A&A, 561, A115
  * Wilson (1927) Wilson, E. B. 1927, Journal of the American Statistical Association, 22, 209
  * Woillez & Bouchet (2020) Woillez, E. & Bouchet, F. 2020, Phys. Rev. Lett., 125, 021101
  * Yang & Zurbenko (2010) Yang, W. & Zurbenko, I. 2010, WIREs Computational Statistics, 2, 340

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

We outline the derivation of the coefficients
$\widetilde{\mathcal{T}}_{\ell,\ell^{\prime}}$ 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, one has

$\mathcal{T}_{1}=\frac{\tilde{\bm{r}}\cdot\tilde{\bm{r}}^{\prime}}{m_{0}}=\frac{\mu\mu^{\prime}}{m_{0}}\bm{\varv}\cdot\bm{\varv}^{\prime}.$
(61)

where $\bm{\varv}=\tilde{\bm{r}}/\mu$ denotes the velocity vector tangent to
the instantaneous Keplerian ellipse of a given planet. In the reference frame
of the Keplerian orbit, with origin at the Sun, $x$-axis directed towards the
pericenter and $z$-axis normal to the orbital plane, the components of the
vector $\bm{\varv}$ read

$\bm{\varv}=\frac{na^{2}}{r}\begin{pmatrix}-\sin E\\\ \sqrt{1-e^{2}}\cos E\\\
0\end{pmatrix},$ (62)

where $E$ is the eccentric anomaly, $n=\sqrt{G(m_{0}+m)/a^{3}}$ the mean
motion and $r=||\bm{r}||$ the heliocentric orbital distance. After
transformation to the fixed reference frame, the velocity components are given
by

$\bm{\varv}=\frac{na^{2}}{r}\mathcal{R}(i,\Omega)\mathcal{R}_{3}(\varpi)\begin{pmatrix}-\sin
E\\\ \sqrt{1-e^{2}}\cos E\\\ 0\end{pmatrix},$ (63)

where the rotation matrices
$\mathcal{R}(i,\Omega)=\mathcal{R}_{3}(\Omega)\mathcal{R}_{1}(i)\mathcal{R}_{3}(-\Omega)$
and $\mathcal{R}_{3}(\varpi)$ are given in (Laskar, 1991; Laskar & Robutel,
1995).

We introduce the eccentric longitude $F=E+\varpi$, which obeys the modified
Kepler’s equation

$\lambda=F+\operatorname{Im}\left(z\,\mathrm{E}^{-jF}\right)=F+\gamma\operatorname{Im}\left(\mathcal{X}\,\mathrm{E}^{-jF}\right),$
(64)

where $z=e\mathrm{E}^{j\varpi}$ and
$\gamma=\sqrt{1-\mathcal{X}\bar{\mathcal{X}}/4}$. By using the fact that
$1-r/a=\gamma\operatorname{Re}\left(\mathcal{X}\,\mathrm{E}^{-jF}\right)$, and
after some algebra, one finds

$\mathcal{R}_{3}(\varpi)\begin{pmatrix}-\sin E\\\ \sqrt{1-e^{2}}\cos E\\\
0\end{pmatrix}=\begin{pmatrix}-\sin
F+(2\gamma)^{-1}\operatorname{Im}\left(\mathcal{X}\right)\left(1-r/a\right)\\\
\phantom{-}\cos
F-(2\gamma)^{-1}\operatorname{Re}\left(\mathcal{X}\right)\left(1-r/a\right)\\\
0\end{pmatrix}.$ (65)

By employing the complex formalism of (Laskar & Robutel, 1995), one can thus
write

$\displaystyle\frac{na^{2}}{r}\mathcal{R}_{3}(\varpi)\begin{pmatrix}-\sin E\\\
\sqrt{1-e^{2}}\cos E\\\ 0\end{pmatrix}=\begin{pmatrix}\mathcal{Z}\\\
0\end{pmatrix},$ (66)
$\displaystyle\mathcal{Z}=jna\left[\frac{a}{r}\mathrm{E}^{jF}-\frac{\mathcal{X}}{2\gamma}\left(\frac{a}{r}-1\right)\right].$

The application of matrix $\mathcal{R}(i,\Omega)$ in Eq. (63) then gives

$\mathcal{R}(i,\Omega)\begin{pmatrix}\mathcal{Z}\\\
0\end{pmatrix}=\operatorname{Re}\left(\mathcal{Z}\bm{\mu}\right),\quad\bm{\mu}=\delta^{-1}\begin{pmatrix}\eta^{2}+\bar{\mathcal{Y}}^{2}\\\
-j\left(\eta^{2}-\bar{\mathcal{Y}}^{2}\right)\\\
-2j\eta\bar{\mathcal{Y}}\end{pmatrix},$ (67)

with $\delta=1-\mathcal{X}\bar{\mathcal{X}}/2$ and
$\eta=\sqrt{1-\mathcal{X}\bar{\mathcal{X}}/2-\mathcal{Y}\bar{\mathcal{Y}}}$.
The scalar product in Eq. (61) thus reads

$\bm{\varv}\cdot\bm{\varv}^{\prime}=\frac{1}{2}\operatorname{Re}\left(\mathcal{Z}\mathcal{Z}^{\prime}\bm{\mu}\cdot\bm{\mu}^{\prime}+\mathcal{Z}\bar{\mathcal{Z}}^{\prime}\bm{\mu}\cdot\bar{\bm{\mu}}^{\prime}\right).$
(68)

The coefficients $\widetilde{\mathcal{T}}_{\ell,\ell^{\prime}}$ are readily
derived once the quantities $a/r$ and $(a/r)\mathrm{E}^{jF}$, appearing in Eq.
(66), are expanded in Fourier series of the mean longitude $\lambda$. One has

$\frac{a}{r}=\sum_{\ell=-\infty}^{+\infty}c_{\ell}\mathrm{E}^{j\ell\lambda},\quad\textrm{with}\quad
c_{\ell}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{a}{r}\mathrm{E}^{-j\ell\lambda}d\lambda.$
(69)

By using Eq. (64) and the fact that $d\lambda/dF=r/a$, a classical calculation
gives $c_{\ell}=J_{\ell}(\ell e)\,\mathrm{E}^{-j\ell\varpi}$, where the
$J_{\ell}$ are the Bessel functions of the first kind. The following
alternative expression of the Fourier coefficients,

$c_{\ell}=\sum_{m=-\infty}^{+\infty}(-j)^{m}J_{m}(\ell\gamma\operatorname{Im}\mathcal{X})\,J_{\ell-m}(\ell\gamma\operatorname{Re}\mathcal{X}),$
(70)

can be straightforwardly expanded in series of $\mathcal{X}$ and
$\bar{\mathcal{X}}$, and is well-defined at zero eccentricity. In a similar
manner, one obtains

$\frac{a}{r}\mathrm{E}^{jF}=\sum_{\ell=-\infty}^{+\infty}\left[\sum_{m=-\infty}^{+\infty}(-j)^{m}J_{m}(\ell\gamma\operatorname{Im}\mathcal{X})\,J_{\ell-1-m}(\ell\gamma\operatorname{Re}\mathcal{X})\right]\mathrm{E}^{j\ell\lambda}.$
(71)

##  Appendix B Kolmogorov-Zurbenko filter

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

Let consider the real-valued time series

$\xi_{n}=\xi(n\Delta t),\quad n\in\mathbb{Z},$ (72)

resulting from the discrete sampling of a continuous signal $\xi(t)$, with
constant sampling rate $\omega_{s}=2\pi/\Delta t$. In the frequency domain,
the time series is characterized by its discrete-time Fourier transform

$\tilde{\xi}(\omega)=\sum_{n=-\infty}^{+\infty}\xi_{n}\mathrm{E}^{-j\omega
n\Delta t},$ (73)

which is a periodic function with period $\omega_{s}$ and can be thus
restricted to the interval $|\omega|\leq\omega_{s}/2$. The action of a finite-
impulse-response filter F on the time series is defined as the discrete
convolution of the signal $\xi_{n}$ with a given finite sequence $d_{n}$,

$\textrm{F}[\xi]_{n}=(d\ast\xi)_{n}=\sum_{m=-M}^{M}d_{m}\xi_{n-m},$ (74)

$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 $\tilde{d}(\omega)$ of
its impulse response, which is called the frequency response of the filter and
gives the spectrum of the output signal as

$\widetilde{\textrm{F}}[\xi](\omega)=\tilde{d}(\omega)\tilde{\xi}(\omega).$
(75)

Ideally, a low-pass filter should be characterized by a null response above a
given cutoff frequency, i.e. $\tilde{d}(\omega)=0$ for $\omega>\omega_{c}$,
while the response should be unitary at lower frequencies, i.e.
$\tilde{d}(\omega)=1$ for $\omega\leq\omega_{c}$. Moreover, the filter must
not change the phase of the Fourier components of the signal, i.e.
$\operatorname{Im}(\tilde{d}(\omega))=0$. This is achieved by requiring
$d_{n}=d_{-n}\in\mathbb{R}$ (Carpino et al., 1987). Nevertheless, a real
filter is rather characterized by: a passband, $|\hat{d}(\omega)-1|\leq\rho$
for $|\omega|\leq\omega_{p}$, $\rho$ being a given maximum loss; a stopband,
$|\hat{d}(\omega)|\leq\alpha$ for $|\omega|\geq\omega_{c}$, $\alpha$ being a
given maximum gain; a transition band between $\omega_{p}$ and $\omega_{c}$,
where the response function smoothly decreases from $\approx$1 to $\approx$0.

Figure 7: 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)$ being represented by the black solid curve.

One of the simplest low-pass filters is the moving average (MA),

$\textrm{MA}_{M}[\xi]_{n}=\frac{1}{2M+1}\sum_{m=-M}^{M}\xi_{n-m},$ (76)

which simply constitutes the local unweighted average of the signal over a
time window $L_{M}\Delta 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. 8). Kolmogorov
proposed to bypass such a problem by applying the moving average iteratively
(Yang & Zurbenko, 2010, and reference therein). The Kolmogorov-Zurbenko filter
is thus defined as

$\displaystyle\textrm{KZ}_{M,k=1}[\xi]=\textrm{MA}_{M}[\xi],$ (77)
$\displaystyle\textrm{KZ}_{M,k}[\xi]=\textrm{MA}_{M}[\textrm{KZ}_{M,k-1}[\xi]],\quad
k\geq 2.$

The output signal of such a filter is given by

$\displaystyle\textrm{KZ}_{M,k}[\xi]_{n}=\frac{1}{(2M+1)^{k}}\sum_{m_{1}=-M}^{M}\cdots\sum_{m_{k}=-M}^{M}\xi_{n-\sum_{\ell=1}^{k}m_{\ell}}=$
(78) $\displaystyle=\sum_{s=-kM}^{kM}\frac{C_{s}^{M,k}}{(2M+1)^{k}}\xi_{n-s},$

where $C_{s}^{M,k}$ is the number of ways of choosing $k$ integers in the
interval $[-M,M]$ such that their sum is equal to $s$. The numbers
$C_{s}^{M,k}$ can be expressed as the coefficients of the finite Laurent
series $\left(\sum_{\ell=-M}^{M}z^{\ell}\right)^{k}$ (Yang & Zurbenko, 2010),
i.e.

$\sum_{s=-kM}^{kM}C_{s}^{M,k}z^{s}=\left(z^{-M}+\ldots+1+\ldots+z^{M}\right)^{k}.$
(79)

The impulse response (IR) of the KZ filter given in Eq. (78),
$\text{IR}^{M,k}_{s}=C_{s}^{M,k}/L_{M}^{k}$, can be interpret 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\gg 1$, such PDF is asymptotically
Gaussian, with zero mean and variance equal to $k\sigma^{2}_{M}$,
$\sigma^{2}_{M}=\left(L_{M}^{2}-1\right)/12$ being the variance of the
discrete uniform distribution over $[-M,M]$. This is shown in Fig. 7.
Therefore, the width of the impulse response scales as $\sqrt{k}M$ for $M,k\gg
1$, even if the filter length $L_{kM}$ is linear in $k$.

Figure 8: Absolute value of the frequency response for the Kolmogorov-Zurbenko
filter (the moving average corresponds to the case $k=1$). The vertical black
dashed line stands for the cutoff frequency $f_{c}$ given by Eq. (85), while
the horizontal coloured ones represent the maximum gains in the stopband
$\alpha_{k}=\eta^{k}$. The vertical coloured dotted lines show the passband
frequencies $f_{p}$ given by Eq. (88) for a loss level of 5%, which is
represented by the horizontal black dotted line.

###  B.1 Frequency response

Table 5: Frequency response of the KZ filter for $M\gg 1$. Maximum gain in the stopband $\alpha_{k}=\eta^{k}$. Ratio of the cutoff frequency $f_{c}$ to the passband frequency $f_{p}$ for two different maximum loss $\rho$. $k$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$\alpha_{k}$ | 2.17 $\cdot$ 10-1 | 4.72 $\cdot$ 10-2 | 1.03 $\cdot$ 10-2 | 2.23 $\cdot$ 10-3 | 4.84 $\cdot$ 10-4
$f_{c}/f_{p}$, $\rho=5\%$ | 4.63 | 6.53 | 7.99 | 9.22 | 10.3
$f_{c}/f_{p}$, $\rho=1\%$ | 10.4 | 14.7 | 18.0 | 20.8 | 23.3

The following derivations describe the behaviour of the KZ filter in the $M\gg
1$ regime. Its frequency response (FR) is simply the $k^{\textrm{th}}$ power
of that of the moving average,

$\textrm{FR}_{M,k}(f)=\left[\frac{\sin(\pi fL_{M})}{L_{M}\sin(\pi
f)}\right]^{k},\quad|f|\leq\frac{1}{2},$ (80)

where $f=\omega\Delta t/2\pi$ is the angular frequency in units of the
sampling rate $\omega_{s}$. The zeros of the frequency response are given by

$\textrm{FR}_{M,k}(f)=0\iff fL_{M}\in\mathbb{Z}\setminus\\{0\\},$ (81)

as shown in Fig. 8. Following (Carpino et al., 1987), we define the
dimensionless cutoff frequency $f_{c}$ as

$|\textrm{FR}_{M,k}(f)|\leq\alpha_{k}\quad\textrm{for }|f|\geq f_{c},$ (82)

where we have defined the maximum gain in the stopband
$\alpha_{k}=|\textrm{FR}_{M,k}(f^{\star}_{M})|$, $f^{\star}_{M}$ being the
frequency of the first local maximum of $|\textrm{FR}_{M,k}(f)|$ for $f>0$. In
the regime $M\gg 1$, such frequency is given by the first positive solution
$x^{\star}$ of the equation

$x=\tan(x),$ (83)

where $x=\pi fL_{M}$. Numerically, one finds $x^{\star}\approx 4.49$.
Therefore, we obtain $f^{\star}_{M}L_{M}=x^{\star}/\pi\approx 1.43$,
independently of $k$, and the gain $\alpha_{k}=\eta^{k}$, with
$\eta=\lvert\sin(x^{\star})\rvert/x^{\star}\approx 0.217$, which is thus
asymptotically independent of $M$. These results are shown in Fig. 8. The
cutoff frequency $f_{c}$ can be computed from Eq. (82) through the first
positive solution $x_{c}$ of the equation

$\eta x=\lvert\sin(x)\rvert,$ (84)

which is found numerically to be $x_{c}\approx 2.55$. The cutoff frequency of
the KZ filter is thus given by

$f_{c}L_{M}=x_{c}/\pi\approx 0.813,$ (85)

independently of $M$ and $k$ in the asymptotic regime $M\gg 1$. Figure 8 and
Table 5 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. (77), the KZ filter turns
out to be 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. Indeed, we define the dimensionless passband
$f_{p}$ as

$\textrm{FR}_{M,k}(f)\geq 1-\rho,\quad\textrm{for }|f|\leq f_{p},$ (86)

where $0<\rho\ll 1$ is a given maximum loss. In the regime $M\gg 1$, by
developing the frequency response in a Taylor series about $x=0$, one has

$1-\frac{x_{p}^{2}}{6}+\frac{x_{p}^{4}}{120}+\mathcal{O}(x_{p}^{6})=(1-\rho)^{1/k},\quad
x_{p}\ll 1.$ (87)

By neglecting terms of order $\mathcal{O}(x_{p}^{6})$, the pertinent solution
is given by

$f_{p}L_{M}=\frac{\sqrt{10}}{\pi}\sqrt{1-\sqrt{1-\frac{6}{5}\bigg{[}1-(1-\rho)^{1/k}\bigg{]}}}.$
(88)

At first order in $\sqrt{\rho/k}$, one obtains

$f_{p}L_{M}=\frac{\sqrt{6}}{\pi}\sqrt{\frac{\rho}{k}}+\mathcal{O}\left(\frac{\rho}{k}\right),$
(89)

which shows that the passband shrinks as $k^{-1/2}$, as shown in Fig. 8. Such
a dependence on $k$ is fortunately sublinear. Generally speaking, a reasonable
choice of the number of iterations should be restricted in our opinion to
$3\leq k\lesssim 5$, as suggested by Table 5. 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 major drawback of the simplicity of the 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 stepsize was constant,
but for the integrations that suffered close planetary encounters or effective
collisions, the stepsize 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 all the solutions were not completed at the time of the
publication of the paper. Indeed, all 2501 runs were started on August 7,
2008. Among them, 2472 were finished in 2008, but 22 ended in 2009 and still 7
in 2010, the last one ending on June, 15, 2010. We present here the full
statistics resulting from these numerical integrations, which have since been
presented in many conferences, but never published. In Table 6 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 the 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.