© ESO, 2012

1. Introduction

After the discovery of the chaotic motion of the planets in the solar system (Laskar 1989, 1990), it has been demonstrated that the eccentricity of Mercury can rise to very high values (Laskar 1994), allowing for collision of the planet with Venus. This was confirmed later on by direct numerical integration in simplified models that neglect general relativity (GR, Batygin & Laughlin 2008; Laskar 2008), and even in the full model that includes the GR contribution (Laskar & Gastineau 2009). Quite surprisingly, the behavior of the system depends strongly on the GR contribution. Indeed, the probability that the eccentricity increases beyond 0.7 in less than 5 Gyr is about 1% in the full model, while it rises to more than 60% when GR is neglected (Laskar 2008; Laskar & Gastineau 2009). This behavior is due to the presence of a secular resonance between the perihelion motions of Mercury (ϖ₁) and Jupiter (ϖ₅) (Batygin & Laughlin 2008; Laskar 2008). Although the origin of the instability in Mercury’s eccentricity is identified, the precise mechanism of the behavior of Mercury’s orbit has not yet been explained in detail. A first attempt was provided by Lithwick & Wu (2011), who analyzed the overlap of resonances in the truncated secular Hamiltonian of degree four in eccentricity and inclination. This model reproduces the small diffusion of Mercury’s eccentricity well, and confirms the important role of the g₁−g₅ and (g₁−g₅)−(s₁−s₂) secular resonances between the precession frequency of the perihelion (g_i) and the regression rate of the ascending node (s_i). However, Lithwick & Wu (2011) does not explain the steady increase in Mercury’s eccentricity beyond 0.7 as observed in Laskar (2008).

In the present paper, our approach is different, since we do not want to be bound to the limitations given by expansions in eccentricity as in Lithwick & Wu (2011) or in the previous work of Laskar (1984, 1989, 1990, 2008). As we are looking for a large excursion of the eccentricity of Mercury, we prefer to derive expressions that are valid for all eccentricities, using the averaged expansions derived in Laskar & Boué (2010). We provide two simple analytical, secular models. The first one is coplanar with only one degree of freedom. It is thus integrable. This model shows that Mercury’s eccentricity can reach values as high as 0.8 if the system is in the vicinity of the g₁−g₅ resonance. The second model includes inclination and has two degrees of freedom. This spatial model exhibits the chaotic behavior of Mercury’s eccentricity in the neighborhood of the separatrix with possibilities of switching between a regime of moderate eccentricity to a regime of large oscillation where Mercury’s eccentricity reaches very high values. In both models, Mercury is treated as a massless particle, while the motion of the other planets in the solar system are given by a single term in their quasiperiodic decompositions taken from Laskar (1990). In simplified terms, the chaos in Mercury’s orbit is due to the proximity of the resonance with the precession motion of Venus excited at the frequency of Jupiter’s precession frequency.

The secular models that we use in this paper rely on a multipolar expansion of the perturbing function, up to order (a₁/a_p)²⁰ in the relativistic case and (a₁/a_p)⁵⁰ in the Newtonian case, where (a_p)_p = 2,8 are the semimajor axes of the planets of the solar system in increasing order. This expansion, subsequently averaged over the mean anomalies of all the planets, allows for arbitrary inclinations and eccentricities as long as no orbit crossing occurs.

In Sect. 2, we derive the equations of motion of the coplanar model. This model is very simple with only one degree of freedom associated to Mercury’s eccentricity and longitude of perihelion. Then, in Sect. 3, we derive the possible trajectories in the phase space using level curves of the Hamiltonian for different orders of expansion. We show that it is necessary to develop the perturbing function up to high orders in (a₁/a_p) in order to reach the asymptotic evolution. We also show that the maximum eccentricity attained with an initial eccentricity of e₁ = 0.2 is on the order of 0.8 within 4 Myr as reported in Laskar (2008). The spatial case is treated in Sect. 4. This nonintegrable model illustrates the generation of chaos in the vicinitiy of the hyperbolic fixed point of the planar model. We conclude in the last section.

2. Coplanar model

2.1. Newtonian interaction

Using the traditional notations where planets have increasing indices with respect to their semimajor axis, we note the barycentric position of the Sun u₀, of Mercury u₁, of Venus u₂, etc. The heliocentric positions of the planets are similarly noted (r_p)_p = 1,8. Since the mass of Mercury is much lower than the masses of the other planets in the solar system, we model Mercury as a massless particle.

Our goal is to study the behavior of Mercury under a known planetary perturbation. As such, all of the u_p and their derivatives, for p ≠ 1, are considered as given functions of the time t. Using the Poincaré heliocentric canonical variables (Laskar & Robutel 1995), the Hamiltonian describing the evolution of Mercury’s trajectory reads as¹$\mathit{Ĥ}\mathit{N}\mathrm{=}\frac{\stackrel{}{{{{r}}_{\mathrm{1}}^{\mathrm{2}}}_{\mathrm{˜}}}}{\mathrm{2}}\mathrm{-}\mathit{G}\frac{{\mathit{m}}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{1}}}\mathrm{-}\mathit{G}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{\left|{{r}}_{\mathrm{1}}\mathrm{-}{{r}}_{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right|}\mathrm{+}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\begin{array}{}\mathrm{˜}\\ {{r}}_{\mathrm{1}}\end{array}\mathrm{·}{u̇}\mathit{p}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(1)where $\begin{array}{}{}_{\mathrm{˜}}\\ {{r}}_{\mathrm{1}}\end{array}\mathrm{=}{u̇}{}_{\mathrm{1}}$ is the conjugate momentum of r₁, G the gravitational constant, m₀ the mass of the Sun, and (m_p)_p = 2,8 the masses of the other planets. The first two terms of this Hamiltonian represent the Keplerian motion and are equal to −Gm₀/(2a₁), where a₁ is the semimajor axis of Mercury.

Considering that Mercury is the closest planet to the Sun, we now expand Ĥ_N formally as a series in (r₁/r_p)_p = 2,8$\mathit{Ĥ}\mathit{N}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathrm{0}}}{\mathrm{2}{\mathit{a}}_{\mathrm{1}}}\mathrm{-}\mathit{G}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{r}}_{\mathit{p}}}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\left(\frac{{\mathit{r}}_{\mathrm{1}}}{{\mathit{r}}_{\mathit{p}}}\right)}^{\mathit{n}}{\mathit{P}}_{\mathit{n}}\left(\frac{{{r}}_{\mathrm{1}}\mathrm{·}{{r}}_{\mathit{p}}}{{\mathit{r}}_{\mathrm{1}}{\mathit{r}}_{\mathit{p}}}\right)\mathrm{+}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\begin{array}{}\mathrm{˜}\\ {{r}}_{\mathrm{1}}\end{array}\mathrm{·}{u̇}\mathit{p}\mathit{,}$(2)where P_n is the Legendre polynomial of order n.

A first approximating model could stop the expansion at the second Legendre polynomial P₂. However, it is well known that the resulting double-averaged secular Hamiltonian does not depend on the perihelia of the outer bodies at this order (Farago & Laskar 2010; Kozai 1962; Lidov 1962; Lidov & Ziglin 1976), so no secular resonance between the perihelia of Mercury and any other planet can be seen there. We thus push the expansion to higher orders. Furthermore, we see in the following that it is necessary to extend the sum at least up to n ≈ 10.

The next step consists in averaging the Hamiltonian (2) over the mean anomalies (M_p)_p = 1,8 of all planets. In the Hamiltonian (2), the terms −Gm₀/r_p obtained for n = 0 do not depend on Mercury’s elements; the terms $\mathrm{(}\begin{array}{}{}_{\mathrm{˜}}\\ {{r}}_{\mathrm{1}}\end{array}\mathrm{·}{{u}}_{\mathit{p}}\mathrm{)}$ vanish when averaged over M₁; the term −Gm₀/(2a₁) becomes constant after averaging over M₁, since a₁ becomes constant. As such, all these terms are dropped in the following expressions. Finally, the averaged expression of the perturbing functions a_p/|r₁−r_p| are given in Laskar & Boué (2010). The resulting Hamiltonian of the coplanar problem is then $\mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathrm{=}\mathrm{-}\mathit{G}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}\sum _{\mathit{n}\mathrm{=}\mathrm{2}}^{\mathrm{\infty }}{\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\right)}^{\mathit{n}}{\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)}\mathit{,}$(3)where $\begin{array}{ccc}{\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}& & \mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,\varpi }\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathit{n}}{\mathit{f}}_{\mathit{n,}\frac{\mathit{n}}{\mathrm{2}}}{\mathit{X}}_{\mathrm{0}}^{\mathit{n,}\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathrm{\right)}{\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathit{p}}\mathrm{\right)}\\ & & \end{array}$(4)In these expressions, a_p, e_p, and ϖ_p are the semimajor axis, the eccentricity, and the longitude of the pericenter of the planet p, respectively. The ${\mathit{X}}_{\mathit{k}}^{\mathit{n,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ are the Hansen coefficients defined by ${\left(\frac{\mathit{r}}{\mathit{a}}\right)}^{\mathit{n}}{\mathrm{e}}^{\mathrm{i}\mathit{mv}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathit{n,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{kM}}\mathit{,}$(5)where ϵ_n = 0 if n is odd, and ϵ_n = 1 if n is even, and ${\mathit{f}}_{\mathit{n,q}}\mathrm{=}\frac{\mathrm{\left(}\mathrm{2}\mathit{q}\mathrm{\right)}\mathrm{!}\mathrm{(}\mathrm{2}\mathit{n}\mathrm{-}\mathrm{2}\mathit{q}\mathrm{)}\mathrm{!}}{{\mathrm{2}}^{\mathrm{2}\mathit{n}}\mathrm{(}\mathit{q}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{(}\mathit{n}\mathrm{-}\mathit{q}\mathrm{)}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(6)To simplify the expansion of the Hamiltonian (2), we take advantage of the eccentricities (e_p)_p = 2,8 of all the planets beyond Mercury remaining low to only keep the linear terms in these eccentricities. In contrast, the expansion remains exact in the eccentricity e₁. With this approximation, up to the octupole order, the averaged Hamiltonian (2) reads as $\begin{array}{}\\ \mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathrm{\approx }\mathrm{-}\frac{\mathit{G}}{\mathrm{8}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}{\mathit{m}}_{\mathit{p}}(\frac{{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{a}}_{\mathit{p}}^{\mathrm{3}}}\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{3}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}\\ \mathrm{-}\frac{\mathrm{15}}{\mathrm{8}}\frac{{\mathit{a}}_{\mathrm{1}}^{\mathrm{3}}}{{\mathit{a}}_{\mathit{p}}^{\mathrm{4}}}{\mathit{e}}_{\mathrm{1}}{\mathit{e}}_{\mathit{p}}\mathrm{(}\mathrm{4}\mathrm{+}\mathrm{3}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)})\mathit{.}\\ \end{array}$(7)More generally, as long as the Hamiltonian is truncated at the first order in e_p, its expression can be written as (see Appendix A) $\begin{array}{}\\ \mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathrm{=}\mathrm{-}\mathit{G}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}({\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\\ \mathrm{-}{\mathit{e}}_{\mathrm{1}}{\mathit{e}}_{\mathit{p}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)})\mathit{,}\\ \end{array}$(8)where P_α(α,x) and Q_α(α,x) are polynomials of degree  [n/2]  and  [(n−1)/2]  in x and of degree 2 [n/2]  and 2 [(n + 1)/2] −1 in α, with n the degree of expansion of the perturbing functions as in Eq. (3). It can be noted that the coefficient of x⁰ in P_α(α,x) is the Taylor expansion of ${\mathit{C}}_{\mathrm{1}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{-}\mathrm{1}\mathrm{=}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}{\mathit{b}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$ −1, and the coefficient of x¹ in P_α(α,x) is the Taylor series of ${\mathit{C}}_{\mathrm{3}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathit{/}\mathrm{2}\mathrm{=}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{8}\mathrm{)}\mathit{\alpha }{\mathit{b}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$, etc., where ${\mathit{b}}_{\mathit{s}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}$ are Laplace coefficients (Laskar & Robutel 1995).

In the Hamiltonian (8), e_p and ϖ_p are functions of time. Neglecting the slow diffusion of the eccentricities of the inner planets, the evolution of each z_p = e_pexpiϖ_p is described well by a quasiperiodic expansion (Laskar 1990). In this study, we focus on the secular resonance g₁−g₅. In the quasiperiodic decomposition of the variables (z_p)_p = 2,8, we thus keep only the terms associated to the eigenmode ${\mathit{z}}_{\mathrm{5}}^{\mathit{\star }}\mathrm{=}\exp \mathrm{i}{\mathit{g}}_{\mathrm{5}}\mathit{t}$ with frequency g₅ ≈ 4.25′′/yr. The others average out and disappear from the Hamiltonian.

The amplitude and the phase of the mode ${\mathit{z}}_{\mathrm{5}}^{\mathit{\star }}$ in the decomposition of each variable z_p are provided in Table 1. We notice that, except for Uranus (p = 7), all the phases are the same and equal to  ≈ 30.6deg. Since, ϕ₇ = 30.6 + 180deg, it is equivalent to take ${\mathit{\varphi }}_{\mathrm{7}}^{\mathrm{\prime }}\mathrm{=}\mathrm{30.6}\mathrm{deg}$ and ${\mathit{A}}_{\mathrm{7}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}{\mathit{A}}_{\mathrm{7}}$. This is the convention that is followed hereafter, but the prime is omitted for clarity. We note ${\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{=}{\mathit{g}}_{\mathrm{5}}\mathit{t}\mathrm{+}{\mathit{\varphi }}_{\mathrm{5}}^{\mathit{\star }}$, where ${\mathit{\varphi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{=}\mathrm{30.6}\mathrm{deg}$. We also use the fact that ${\mathrm{\Lambda }}_{\mathrm{1}}\mathrm{=}\sqrt{\mathit{G}{\mathit{m}}_{\mathrm{0}}{\mathit{a}}_{\mathrm{1}}}$ is constant in the secular problem to rescale the Hamiltonian by this quantity as it simplifies the following expressions. We note $\begin{array}{}{}_{\mathrm{ˇ}}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{=}\mathit{Ĥ}{}_{\mathit{N,}\mathrm{plan}}\mathit{/}{\mathrm{\Lambda }}_{\mathrm{1}}\mathit{.}$Given that Gm_p/(Λ₁a_p) = n₁a₁m_p/(a_pm₀), the resonant Hamiltonian now reads as $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{=}\mathrm{-}{\mathit{n}}_{\mathrm{1}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}({\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\\ \mathrm{-}{\mathit{e}}_{\mathrm{1}}{\mathit{A}}_{\mathit{p}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{)})\mathit{.}\\ \end{array}$(9)With this convention, the time should also be rescaled by the same factor Λ₁ to keep the canonical form of the equations of motion, unless the canonical variables are modified as follows. We let $\begin{array}{}\\ \mathit{Î}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{1}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\right)\mathit{,}\\ \mathit{\theta ̂}\mathrm{=}\mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\\ \end{array}$(10)be the canonical conjugated variables of the initial Hamiltonian Ĥ_N,plan of (Eq. (8)). It can be easily shown that $\mathrm{\left(}\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathrm{\right)}\mathrm{=}\left(\frac{\mathit{Î}}{{\mathrm{\Lambda }}_{\mathrm{1}}}\mathit{,}\mathit{\theta ̂}\right)$(11)are canonical conjugated variables of the Hamiltonian $\begin{array}{}{\mathrm{ˇ}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ (Eq. (9)) without rescaling the time t. The equations of motion are thus $\begin{array}{}\\ \frac{\mathrm{d}\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }\begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}}{\mathit{\partial }\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}}\mathrm{\left(}\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathit{,t}\mathrm{\right)}\mathit{,} \frac{\mathrm{d}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial }\begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}}{\mathit{\partial }\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}}\mathrm{\left(}\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathit{,t}\mathrm{\right)}\mathit{.}\\ \end{array}$(12)At this point, it is interesting to evaluate the contribution of each planet to the octupole interaction with Mercury. Although the resonance g₁−g₅ involves only the precession frequencies associated to Mercury and Jupiter, the eigenmode ${\mathit{z}}_{\mathrm{5}}^{\mathit{\star }}$ is present in the quasiperiodic decomposition of the motion of all the planets. Furthermore, the amplitudes of this mode are very similar and vary only within a factor 2.5 between 0.019 and 0.044, except for Neptune whose amplitude is only 0.002 (Table 1). From the expression of the Hamiltonian (9), knowing that the lowest degree of Q(α,x) in α is 3, the contribution of the octupole terms can be estimated using the parameters $\mathrm{(}{\mathit{\epsilon }}_{\mathit{p}}^{\mathrm{octu}}{\mathrm{)}}_{\mathit{p}\mathrm{=}\mathrm{2}\mathit{,}\mathrm{8}}$ given by ${\mathit{\epsilon }}_{\mathit{p}}^{\mathrm{octu}}\mathrm{=}\frac{\mathrm{15}}{\mathrm{64}}\left(\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\right)\left(\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\right){\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\right)}^{\mathrm{4}}{\mathit{A}}_{\mathit{p}}\mathit{.}$(13)The values taken by these parameters are gathered in Table 1. The maximal amplitude is due to Venus with ${\mathit{\epsilon }}_{\mathrm{2}}^{\mathrm{octu}}\mathrm{\approx }\mathrm{0.0012}$ followed by the Earth-Moon barycenter and Jupiter with ${\mathit{\epsilon }}_{\mathrm{3}}^{\mathrm{octu}}\mathrm{\approx }{\mathit{\epsilon }}_{\mathrm{5}}^{\mathrm{octu}}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-4}$. Thus, the strongest perturbation on Mercury’s orbit comes from the precession of Venus excited by Jupiter. a₁/a₂ ≈ 0.53 is not very small explains why it is necessary to perform the expansion of the perturbing function up to a high order in the ratio of the semimajor axes. For completeness, Table 1 also provides the quadrupole contribution of each planet given by ${\mathit{\epsilon }}_{\mathit{p}}^{\mathrm{quad}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{8}}\left(\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\right)\left(\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\right){\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\right)}^{\mathrm{3}}\mathrm{·}$(14)Since we have shown that several planets play a role in the evolution of Mercury’s eccentricity, in the following we always consider all the perturbers from Venus to Neptune. To simplify the notations, we define two new polynomials P and Q of the eccentricity e₁ alone, given by (Eq. (9)) $\begin{array}{}\\ \mathrm{P}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}{\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,x}\right)\mathit{,}\\ \mathrm{Q}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,x}\right)\mathrm{\times }{\mathit{A}}_{\mathit{p}}\mathit{.}\\ \end{array}$(15)Then, the resonant Hamiltonian (9) reads as $\begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathrm{\left(}\mathrm{P}\mathrm{\right(}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{e}}_{\mathrm{1}}\mathrm{Q}\mathrm{(}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{\left)}\cos \mathrm{\right(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{.}$(16)The numerical values of the coefficients of the polynomials P and Q are given in Appendix B.

The Hamiltonian $\begin{array}{}{\mathrm{ˇ}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ (16) has one and a half degrees of freedom, but it can be reduced to a one degree of freedom Hamiltonian $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ after some modifications.For that, we first make the Hamiltonian autonomous by a adding a momentum $\begin{array}{}{\mathrm{ˇ}}_{}\\ \mathit{T}\end{array}$ conjugated to time t. Then, we perform a canonical transformation where the old variables are $\mathrm{(}\begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathrm{=}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{)}\mathrm{;} \mathrm{\left(}\begin{array}{}\mathrm{ˇ}\\ \mathit{T}\end{array}\mathit{,t}\mathrm{\right)}\mathit{,}$(17)and the new ones are (I,Δϖ) and $\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{T}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{t}\end{array}\mathrm{\right)}$, defined by $\begin{array}{}\\ \mathrm{\Delta }\mathit{\varpi }\mathrm{=}\mathrm{-}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathit{t}\mathrm{-}{\mathit{\varphi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{\equiv }{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{t}\end{array}\mathrm{=}\mathit{t}\mathit{.}\\ \end{array}$(18)For this transformation to be canonical, the momentums should verify $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathrm{=}\mathrm{-}\mathit{I}\mathrm{\equiv }\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{T}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathit{I}\mathrm{+}\begin{array}{}\mathrm{˜}\\ \mathit{T}\end{array}\mathit{.}\\ \end{array}$(19)The new Hamiltonian $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ expressed in the new variables does not depend on the cyclic coordinate $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{t}\end{array}$. Thus, its conjugated momentum $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{T}\end{array}$ is an integral of the motion and can be dropped. Up to a constant, the new Hamiltonian is then$\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathrm{P}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{g}}_{\mathrm{5}}{\mathit{e}}_{\mathrm{1}}\mathrm{Q}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\cos \mathrm{\Delta }\mathit{\varpi }\mathit{.}$(20)This one degree of freedom Hamiltonian is integrable. The orbits are given by $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{=}\mathrm{Cte}$. The temporal evolution of the eccentricity e₁ and of the resonant angle Δϖ are deduced from the canonical equations of motion. This leads to $\begin{array}{}\\ \frac{\mathrm{d}{\mathit{e}}_{\mathrm{1}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}{{\mathit{e}}_{\mathrm{1}}}\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}}{\mathit{\partial }\mathrm{\Delta }\mathit{\varpi }}\mathit{,}\\ \frac{\mathrm{d\Delta }\mathit{\varpi }}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}{{\mathit{e}}_{\mathrm{1}}}\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}}{\mathit{\partial }{\mathit{e}}_{\mathrm{1}}}\mathrm{·}\\ \end{array}$(21)

2.2. General relativistic precession

The secular effect of relativity is described by (e.g., Touma et al. 2009), ${\mathit{H}}_{\mathit{R}}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\mathit{,}$(22)where ${\mathit{g}}_{\mathit{r}}\mathrm{=}\mathrm{3}\mathrm{(}\mathit{G}{\mathit{m}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{(}{\mathrm{\Lambda }}_{\mathrm{1}}{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\mathrm{)}$, and c is the speed of light. In the case of Mercury, we have g_r ≈ 0.41′′/yr. The total Hamiltonian $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}\end{array}\mathrm{=}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{+}{\mathit{H}}_{\mathit{R}}$, including the Newtonian interaction and general relativity, becomes $\begin{array}{}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\left(\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}\mathrm{P}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\mathit{e}}_{\mathrm{1}}\mathrm{Q}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\cos \mathrm{\Delta }\mathit{\varpi }\right)\\ \mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{·}\\ \end{array}$(23)

Fig. 1 Level curves of the Hamiltonian in the plane (e1cosΔϖ,e1sinΔϖ). a) Quadrupole expansion with g1 = 5′′/yr. b) Octupole expansion with g1 = 5′′/yr. c), d), e), f) expansion up to the order n = 20 with g1 = 2.5,3.6,4.0,5.0′′/yr, respectively.

Fig. 2 Evolution of the fixed points as a function of the precession frequency of Mercury’s orbit at zero eccentricity and for different expansion orders. Dashed and dotted curves correspond to unstable points, while the solid ones show the positions of the stable points. The equilibrium points are labeled with roman letters as in Fig. 1. The thin vertical line indicates the resonance g1 = g5.

2.3. Additional control term

Using the Newtonian Hamiltonian (20), or the relativistic one (23), the system is far from the g₁−g₅ resonance, in a configuration where the amplitude of oscillation of the eccentricity is small (see next section). Indeed, taking an order of expansion n = 50, we get g₁ = 5.54′′/yr in the Newtonian case and g₁ = 5.96′′/yr in the relativistic case, whereas the resonance occurs in the vicinity of g₁ = g₅ = 4.25′′/yr. This result is slightly different, but representative of the present behavior of Mercury’s eccentricity. Indeed, the present value of g₁ for the solar system with GR is g₁ = 5.59′′/yr (Laskar et al. 2004), the differences with the present model being due to the simplifications that are made here. To recover a model that is dynamically close to the real solar system, we add a correction to the Hamiltonian which changes the value of g₁ by an increment δ_g. In addition, owing to the slow chaotic diffusion in the inner solar system, Mercury’s precession frequency g₁ can change quasi-randomly and come close to the value of g₅, then leading to high unstability (Batygin & Laughlin 2008; Laskar 1990, 2008; Laskar & Gastineau 2009). To explore the evolution of Mercury’s eccentricity for different values of g₁ around the resonant frequency g₅, the above-mentioned correction is added to both $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}\end{array}$; i.e., ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\delta }}_{\mathit{g}}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}$(24)and ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}\end{array}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\delta }}_{\mathit{g}}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathit{.}$(25)In Eqs. (24) and (25), the factor δ_g controls the precession frequency g₁. For δ_g = 0, we recover the values obtained with $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}\end{array}$, respectively. Otherwise the value of g₁ is incremented by δ_g. Appendix C provides the values of δ_g that put the system in exact resonance for all orders of expansion of the Hamiltonians.

3. Eccentricity behavior

Here, we analyze the possible trajectories and evolutions of Mercury’s eccentricity described by the Hamiltonians H_N,plan (24) and H_R,plan (25) obtained in the previous section. Since these Hamiltonians are integrable (they have only one degree of freedom), all the orbits are necessarily regular. Thus, the goal of this section is not to reproduce the chaotic behavior of Mercury’s eccentricity e₁, but to show that e₁ can actually reach values as high as observed by Laskar (2008) (see Fig. 3a) when the system is in the vicinity of the resonance g₁ = g₅. The chaotic evolution will be analyzed in the spatial case (see Sect. 4), where a new degree of freedom is added.

3.1. Dependency with the order of expansion

The Hamiltonians of the previous section have been computed as series in power of the semimajor axis ratios (a₁/a_p)_p = 2,8 for any order n. Here, we study the effect of the truncation of these series. Hereafter, we note with a superscript (n) any Hamiltonian expanded up to the order n, e.g., ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$ or ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$. In a first step, we focus only on the relativistic Hamiltonian which is more complete.

Within the quadrupole approximation, the Hamiltonian ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ reduces to a simple expression $\begin{array}{}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\left(\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}{\mathit{\epsilon }}^{\mathrm{quad}}\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{3}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\right)\\ \mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\delta }}_{\mathit{g}}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathit{,}\\ \end{array}$(26)where ${\mathit{\epsilon }}^{\mathrm{quad}}\mathrm{=}{\sum }_{\mathit{p}}{\mathit{\epsilon }}_{\mathit{p}}^{\mathrm{quad}}$. This Hamiltonian is independent of Δϖ, thus e₁ remains constant. All the trajectories in the plane (x = e₁cosΔϖ, y = e₁sinΔϖ) are circles centered on the origin (0,0) with an eventually infinite period when . Figure 1a illustrates this result for g₁ = 5′′/yr. The position of the fixed points is represented by red dashed curves. Within this approximation, as noted before, there is no possible resonance between ϖ₁ and ${\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}$ and thus, no possible increase in the eccentricity.

The next level of approximation is the octupole. The corresponding Hamiltonian reads as $\begin{array}{}\\ {\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}(\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}{\mathit{\epsilon }}^{\mathrm{quad}}\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{3}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}\\ \mathrm{-}{\mathit{\epsilon }}^{\mathrm{octu}}\mathrm{(}\mathrm{4}\mathrm{+}\mathrm{3}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}\cos \mathrm{\Delta }\mathit{\varpi })\mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{+}{\mathit{\delta }}_{\mathit{g}}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathit{,}\\ \end{array}$(27)with ${\mathit{\epsilon }}^{\mathrm{octu}}\mathrm{=}{\sum }_{\mathit{p}}{\mathit{\epsilon }}_{\mathit{p}}^{\mathrm{octu}}$. In that case the phase space can be much more complicated with five fixed points and two sets of separatrices (see Fig. 1b obtained with g₁ = 5′′/yr). The stable fixed points are represented by black dots, and the unstable ones are noted with crosses. The dashed curves correspond to the separatrices. The positions and existence of the fixed points depend on the value of the precession frequency g₁. This is depicted in the subpanel labeled α³ of Fig. 2. The curves represent the position on the x axis of the fixed points of the phase space as a function of g₁. The labels in roman numerals qualifying each fixed point are identical to those in Fig. 1.

It is interesting to have a look at all of Fig. 2. Indeed, as the order n of the expansion of the Hamiltonian ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$ increases, the degrees of the polynomials P_α(α,e²) and Q_α(α,e²) increase both for α and e². Then, the topology of the phase space of the Hamiltonian ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$ evolves as shown in Fig. 2. Up to n = 5, five fix points – three stable points and two unstable points – coexist within 4.5 ≲ g₁ ≲ 5.0′′/yr, while the system contains at most three fix points – two stable points and one unstable point – for n ≥ 6. Moreover, the hyperbolic point (labeled IV) located at e₁cosΔϖ > 0, i.e. Δϖ = 0, disappears when n ≥ 8. The asymptotic topology is reached at n ≈ 10 from which positions of the fixed points do not evolve significantly up to n = 20.

The Figs. 1c, 1d, 1e, and 1f display the allowed trajectories of Mercury’s eccentricity for n = 20 which is assumed to be representative of the asymptotic behavior. The values of g₁ are 2.5, 3.6, 4.0, and 5.0′′/yr, respectively. At this order n = 20, orbits are very similar to those of a simple pendulum with a large resonant island allowing for large increases in eccentricity beyond e₁ = 0.8. The separatrix disappears at g₁ ≥ 3.84′′/yr but high eccentricity excursions are still possible, to a smaller extent, since the elliptic fixed point (V) remains significantly offset with respect to the center of the phase space as long as g₁ ≲ 4.5′′/yr. For example, with g₁ = 4.0′′/yr (Fig. 1e), a trajectory starting within e₁ ≤ 0.2 can reach e₁ ≈ 0.7.

To conclude, the asymptotic phase space, which must reproduce Mercury’s eccentricity behavior as faithfully as possible, is obtained for n ≈ 20. Moreover, at this order of expansion, Mercury’s eccentricity is able to increase from e₁ ≈ 0.2 up to e₁ ≈ 0.8 as observed in the numerical simulations reported by Laskar (2008) (see Fig. 3).

3.2. Temporal evolution

The level curves of the Hamiltonian ${\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{20}\mathrm{\right)}}$ provide the trajectories of Mercury’s eccentricity in the phase space. In the previous section, we saw that the eccentricity is able to vary between 0.2 and 0.8 as observed in Laskar (2008) (Fig. 3a). We now check whether the timescale of the evolution matches the 4 Myr found by Laskar (2008). For this purpose, we integrate the following equations $\begin{array}{}\\ \frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{20}\mathrm{\right)}}}{\mathit{\partial y}}\mathrm{;}\\ \frac{\mathrm{d}\mathit{y}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{H}}_{\mathit{R,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{20}\mathrm{\right)}}}{\mathit{\partial x}}\mathit{,}\\ \end{array}$(28)where x = e₁cosΔϖ and y = e₁sinΔϖ. They are equivalent to the equations of motion (21), but without singularities at the origin e₁ = 0.

Fig. 3 Comparison between the eccentricity evolution observed in Laskar (2008)(a), and produced by the simple model of this work (b).

Fig. 4 Level curves of the Hamiltonian without relativity in the plane (e1cosΔϖ,e1sinΔϖ) for two orders of expansion: n = 20 (left), n = 50 (right). In both cases, g1 = 3.6′′/yr.

An example of temporal evolution given by the numerical integration of Eq. (28) is displayed in Fig. 3b. The initial conditions are e₁ = 0.2, Δϖ = 110deg, and g₁ = 3.68′′/yr. The comparison of Fig. 3a with Fig. 3b shows that the simple one degree of freedom model described in this study is able to reproduce the main evolution of Mercury’s eccentricity with the correct amplitude and timescale. Only once the eccentricity reaches a value close to 0.8, Mercury’s orbit becomes highly unstable due to close encounters with Venus, and the evolution observed in Laskar (2008) starts to differ significantly from the simple model. It should be noted that the small oscillations present in the integration of the full solar system (Fig. 3a) do not appear in the simple model since we keep only one single eigen frequency to describe the evolution of the outer planets eccentricities.

Fig. 5 Top: evolution of the fixed points of H(50) without relativity. Bottom: evolution of the eccentricity with time, using the simple model without relativity.

3.3. Without relativity

We now consider the evolution of Mercury’s eccentricity without general relativity, i.e., described by the Hamiltonian ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$ (24). Figure 4 shows the level curves of ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{20}\mathrm{\right)}}$ and ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{50}\mathrm{\right)}}$. The two figures are similar to the subpanel α²⁰ of Fig. 1 within the region of low eccentricities (e₁ ≲ 0.6). But for higher eccentricities, two new fixed points appear in the Newtonian case with n = 20. It is then necessary to extend the expansion up to the order n ≈ 50 to get the asymptotic topology of the phase space. Once the limit is closely reached, the result matches the relativistic ones well. The only difference is in the value of δ_g: for a given frequency g₁, it should be increased by g_r = 0.41′′/yr with respect to the relativistic case (see Table C.1).

Figure 5 shows the positions of the fixed points and the temporal evolution of the eccentricity driven by ${\mathit{H}}_{\mathit{N,}\mathrm{plan}}^{\mathrm{\left(}\mathrm{50}\mathrm{\right)}}$. There is no significant difference with respect to the relativistic case.

4. Spatial model

The evolutions studied in the previous sections are perfectly regular. Indeed, the Hamiltonian contains only one degree of freedom and is thus integrable. The goal was to show that it is possible to explain the large increase in Mercury’s eccentricity observed in numerical simulations (e.g., Laskar 2008). Our aim is now to introduce an additional degree of freedom to make the system non integrable and to reproduce a chaotic evolution expected in the vicinity of the separatrix of the g₁−g₅ resonance.

To do so, we add inclination in our system, and focus on the term related to (g₁−g₅)−(s₁−s₂). This resonant term is at present in libration, and was identified as one of the main source of chaos in the solar system (Laskar 1990). It was indeed also computed analytically in Laskar (1984), where it was identified as the major obstacle for the convergence of the secular perturbation series. The Hamiltonian of the inclined problem is (see Appendix D) $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{R,}\mathrm{inc}}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\mathrm{-}{\mathit{n}}_{\mathrm{1}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}(\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}{\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\right)}^{\mathit{n}\mathrm{+}\mathrm{1}}\\ \mathrm{\times }\hspace{0.17em}{\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{A}}_{\mathit{p}}\mathit{,}{\mathit{i}}_{\mathrm{1}}\mathit{,}{\mathit{I}}_{\mathit{p}}\mathit{,}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathit{,}{\mathrm{\Omega }}_{\mathrm{1}}\mathrm{-}{{\mathrm{\Omega }}_{\mathrm{2}}^{\mathit{\star }}}^{\mathrm{)}})\mathit{,}\\ \end{array}$(29)where i_p and Ω_p are the inclination and the longitude of the ascending node of the planet p, and I_p and ${\mathrm{\Omega }}_{\mathrm{2}}^{\mathit{\star }}$ are the amplitude and the argument of the term precessing at the frequency s₂ in the quasiperiodic decomposition of sin(i_p/2)exp(iΩ_p), respectively. In (29), both ${\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{=}{\mathit{g}}_{\mathrm{5}}\mathit{t}\mathrm{+}{\mathit{\varphi }}_{\mathrm{5}}^{\mathit{\star }}$ and ${\mathrm{\Omega }}_{\mathrm{2}}^{\mathit{\star }}\mathrm{=}{\mathit{s}}_{\mathrm{2}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathrm{2}}^{\mathit{\star }}$ are functions of time t. To reduce the number of degrees of freedom from 2.5 to 2, we apply transformations similar to those of the planar case (Sect. 2). We first add to the Hamiltonian the momentum $\begin{array}{}{\mathrm{ˇ}}_{}\\ \mathit{T}\end{array}$ conjugated to the time t. The old variables of the subsequent canonical transformation are $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathrm{=}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathit{,} \begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}{\mathit{\varpi }}_{\mathrm{1}}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{J}\end{array}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}\cos {\mathit{i}}_{\mathrm{1}}\mathrm{)}\mathit{,} \begin{array}{}\mathrm{ˇ}\\ \mathit{\psi }\end{array}\mathrm{=}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{1}}\mathit{,}\begin{array}{}\mathrm{ˇ}\\ \mathit{T}\end{array}\mathit{,}\mathit{t,}\\ \end{array}$(30)and the new variables are denoted (I,Δϖ), (J,ΔΩ), $\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{T}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{t}\end{array}\mathrm{\right)}$. We set $\begin{array}{}\\ \mathrm{\Delta }\mathit{\varpi }\mathrm{=}\mathrm{-}\begin{array}{}\mathrm{ˇ}\\ \mathit{\theta }\end{array}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathit{t}\mathrm{-}{\mathit{\varphi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{\equiv }{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathit{,}\\ \mathrm{\Delta \Omega }\mathrm{=}\mathrm{-}\begin{array}{}\mathrm{ˇ}\\ \mathit{\psi }\end{array}\mathrm{-}{\mathit{s}}_{\mathrm{2}}\mathit{t}\mathrm{-}{\mathit{\phi }}_{\mathrm{2}}^{\mathit{\star }}\mathrm{\equiv }{\mathrm{\Omega }}_{\mathrm{1}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{2}}^{\mathit{\star }}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ \mathit{t}\end{array}\mathrm{=}\mathit{t}\mathit{.}\\ \end{array}$(31)The associated conjugated momenta are given by $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ \mathit{I}\end{array}\mathrm{=}\mathrm{-}\mathit{I,}\\ \begin{array}{}\mathrm{ˇ}\\ \mathit{J}\end{array}\mathrm{=}\mathrm{-}\mathit{J,}\\ \begin{array}{}\mathrm{ˇ}\\ \mathit{T}\end{array}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathit{I}\mathrm{-}{\mathit{s}}_{\mathrm{2}}\mathit{J}\mathrm{+}\begin{array}{}\mathrm{˜}\\ \mathit{T}\end{array}\mathit{.}\\ \end{array}$(32)Since the new Hamiltonian $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{H}}_{\mathit{R,}\mathrm{inc}}\end{array}$ is independent of t, its conjugate momentum $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{T}\end{array}$ is an integral of the motion. From the expressions of $\begin{array}{}{\mathrm{ˇ}}_{}\\ \mathit{I}\end{array}$ and $\begin{array}{}{\mathrm{ˇ}}_{}\\ \mathit{J}\end{array}$ (30), we thus get, up to a constant, $\begin{array}{}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{R,}\mathrm{inc}}\end{array}\mathrm{=}\left(\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathrm{+}\mathrm{2}{\mathit{s}}_{\mathrm{2}}{\sin }^{\mathrm{2}}\frac{{\mathit{i}}_{\mathrm{1}}}{\mathrm{2}}\right)\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}{\mathit{g}}_{\mathit{r}}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}}\\ \mathrm{-}{\mathit{n}}_{\mathrm{1}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}{\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\right)}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{A}}_{\mathit{p}}\mathit{,}{\mathit{i}}_{\mathrm{1}}\mathit{,}{\mathit{I}}_{\mathit{p}}\mathit{,}\mathrm{\Delta }\mathit{\varpi ,}\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{)}}\mathit{.}\\ \end{array}$(33)To obtain the final Hamiltonian H_R,inc in the spatial case, we add the control term and get ${\mathit{H}}_{\mathit{R,}\mathrm{inc}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{H}}_{\mathit{R,}\mathrm{inc}}\end{array}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\delta }}_{\mathit{g}}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathit{.}$(34)The Hamiltonian (34) now has two degrees of freedom. The variables are (e₁,Δϖ) and (i₁,ΔΩ). To perform numerical integrations, we use nonsingular rectangular coordinates $\begin{array}{}\\ \mathit{k}\mathrm{=}{\mathit{e}}_{\mathrm{1}}\cos \mathrm{\Delta }\mathit{\varpi ,} \mathit{h}\mathrm{=}{\mathit{e}}_{\mathrm{1}}\sin \mathrm{\Delta }\mathit{\varpi }\\ \mathit{q}\mathrm{=}\sin \frac{{\mathit{i}}_{\mathrm{1}}}{\mathrm{2}}\cos \mathrm{\Delta \Omega }\mathit{,} \mathit{p}\mathrm{=}\sin \frac{{\mathit{i}}_{\mathrm{1}}}{\mathrm{2}}\sin \mathrm{\Delta \Omega }\mathit{.}\\ \end{array}$(35)With $\mathit{\chi }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}$, the equations of motion read as (e.g. Bretagnon 1974) $\begin{array}{}\\ \frac{\mathrm{d}\mathit{k}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{+}\mathit{\chi }\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial h}}\mathrm{+}\frac{\mathit{h}}{\mathrm{2}\mathit{\chi }}\left(\mathit{q}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial q}}\mathrm{+}\mathit{p}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial p}}\right)\mathit{,}\\ \frac{\mathrm{d}\mathit{h}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{\chi }\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial k}}\mathrm{-}\frac{\mathit{k}}{\mathrm{2}\mathit{\chi }}\left(\mathit{q}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial q}}\mathrm{+}\mathit{p}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial p}}\right)\mathit{,}\\ \frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}\mathit{\chi }}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial p}}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}\mathit{\chi }}\left(\mathit{h}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial k}}\mathrm{-}\mathit{k}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial h}}\right)\mathit{,}\\ \frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{4}\mathit{\chi }}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial q}}\mathrm{-}\frac{\mathit{p}}{\mathrm{2}\mathit{\chi }}\left(\mathit{h}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial k}}\mathrm{-}\mathit{k}\frac{\mathit{\partial }{\mathit{H}}_{\mathrm{inc}}}{\mathit{\partial h}}\right)\mathrm{·}\\ \end{array}$(36)Since Venus is the only planet with a significative amplitude associated to the frequency s₂, it is also the only planet for which the inclination is taken into account. The initial inclination of Mercury is taken from Laskar (1990) where all the terms in its frequency decomposition, except the constant one, have been added. Doing so, Mercury’s initial inclination is measured with respect to the invariant plane.

Figure 6 shows the results of a numerical integration of the inclined system. A large chaotic zone appears in the vicinity of the separatrix, as expected. Mercury’s eccentricity now switches randomly between low-amplitude circulations and high-amplitude excursions, reaching values that are close to 1.

It should be noted that the solar system is not presently in this chaotic zone. As described elsewhere (Batygin & Laughlin 2008; Laskar 1990, 2008; Laskar & Gastineau 2009; Lithwick & Wu 2011), the system is at present in a secular resonance (with resonant argument (g₁−g₅)−(s₁−s₂)), in a state of slow chaos, with slow chaotic diffusion. This diffusion will quasi-randomly change the value of g₁, which can then approach resonance with g₅. The phase diagram of Fig. 6 describes the behavior of Mercury’s eccentricity once this slow diffusion has brought the system into the vicinity of the g₁−g₅ resonance.

Fig. 6 Top: trajectory of Mercury’s eccentricity in the plane (e1cosΔϖ,e1sinΔϖ) for the inclined system. Bottom: evolution of Mercury’s eccentricity as a function of time in the same system.

5. Conclusion

In this study, we developed a simple model to account for the increase in Mercury’s eccentricity due to the resonance g₁ = g₅. This coplanar model is based on the expansion of the perturbing function with respect to the semimajor axis ratios, and exact in eccentricity. In the resulting Hamiltonian, we kept only one term in the quasiperiodic evolution of the outer planets’ eccentricity. We found that this approximation is sufficient to reproduce an increase in Mercury’s eccentricity up to 0.8. But we also noticed that it is necessary to extend the expansion up to the order n ≈ 20 and n ≈ 50 in the relativistic and in the Newtonian cases, respectively. The explicit form of this secular resonant Hamiltonian is provided in Appendix B.

The asymptotic topologies of the phase space are very similar no matter whether the relativity is taken into account or not, so these two cases just depend on the value of the resonant frequency (Table C.1). In both cases, the system is a one degree of freedom system that is integrable with a separatrix (Figs. 1d, 4). The eccentricity of the trajectories in the vicinity of this separatrix can rise to very high values, up to 0.8. Moreover, the timescale of Mercury’s evolution is in very good agreement with the one observed in the numerical integration of the full model Laskar (2008) (Fig. 3). With this integrable model, the behavior of the numerical solutions computed in Laskar (2008), Batygin & Laughlin (2008), Laskar & Gastineau (2009) are understood, but not the transition from a low-eccentricity regime to a high-eccentricity regime. To obtain such a transition, it is necessary to include an additional degree of freedom in the system, which transforms the separatrix into a chaotic zone.

Since we know that the resonant term associated with (g₁−g₅)−(s₁−s₂) has large amplitude (Laskar 1984, 1990; Lithwick & Wu 2011), we added this single term in the spatial problem which corresponds to our second model (Sect. 4). As expected, in this nonintegrable problem, the separatrix is replaced by a significant chaotic zone where transitions from low-eccentricity to high-eccentricity regimes occur quasi-randomly, as observed in the full system, with maximal eccentricity close to 1 (Fig. 6).

Acknowledgments

This work has been supported by PNP-CNRS, by the CS of the Paris Observatory, by PICS05998 France-Portugal program, by the European Research Council/European Community under the FP7 through a Starting Grant, as well as in the form of grant reference PTDC/CTE-AST/098528/2008, funded by the Fundação para a Ciência e a Tecnologia (FCT), Portugal.

References

Appendix A: Explicit expression of the Hamiltonian development in the planar case

The secular Hamiltonian Ĥ_plan of a restricted two-planet system where the massless body is on the inner orbit reads as $\mathit{Ĥ}\mathrm{plan}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}\sum _{\mathit{n}\mathrm{=}\mathrm{2}}^{\mathrm{\infty }}{\mathit{\alpha }}^{\mathit{n}}{\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)}\mathit{.}$(A.1)The indices 1 and p refer to the inner massless body and to the outer massive planet, respectively; m_p is the mass of the outer planet; α = a₁/a_p is the semimajor axis ratio; (e_k)_k = 1,p and (ϖ_k)_k = 1,p are the eccentricities and the longitudes of periastron of the two planets.

In (A.1), the ${\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,\varpi }\mathrm{\right)}$ are given by Laskar & Boué (2010)$\begin{array}{}\\ {\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,\varpi }\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathit{n}}{\mathit{f}}_{\mathit{n,}\frac{\mathit{n}}{\mathrm{2}}}{\mathit{X}}_{\mathrm{0}}^{\mathit{n,}\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathrm{\right)}{\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathit{p}}\mathrm{\right)}\\ \mathrm{+}\sum _{\mathit{q}\mathrm{=}\mathrm{0}}^{\mathrm{\left[}\mathrm{\right(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{]}}\mathrm{2}{\mathit{f}}_{\mathit{n,q}}{\mathit{X}}_{\mathrm{0}}^{\mathit{n,n}\mathrm{-}\mathrm{2}\mathit{q}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathrm{\right)}{\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,n}\mathrm{-}\mathrm{2}\mathit{q}}\mathrm{\left(}{\mathit{e}}_{\mathit{p}}\mathrm{\right)}\cos \mathrm{\left(}\mathrm{\right(}\mathit{n}\mathrm{-}\mathrm{2}\mathit{q}\mathrm{\left)}\mathit{\varpi }\mathrm{\right)}\mathit{,}\\ \end{array}$(A.2)where ϵ_n = 1 if n is even and 0, otherwise ${\mathit{f}}_{\mathit{n,q}}\mathrm{=}\frac{\mathrm{\left(}\mathrm{2}\mathit{q}\mathrm{\right)}\mathrm{!}\mathrm{(}\mathrm{2}\mathit{n}\mathrm{-}\mathrm{2}\mathit{q}\mathrm{)}\mathrm{!}}{{\mathrm{2}}^{\mathrm{2}\mathit{n}}\mathrm{(}\mathit{q}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{(}\mathit{n}\mathrm{-}\mathit{q}\mathrm{)}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(A.3)and ${\mathit{X}}_{\mathrm{0}}^{\mathit{n,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ are Hansen coefficients.

Now, we assume that the eccentricity of the outer planet remains low, and we only keep the linear terms in e_p. Since the lowest power in eccentricity of ${\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ is e|m|, all the terms in the sum (A.2) are dropped, except those for which n−2q ≤ 1. Furthermore, from the explicit expressions of the Hansen coefficients (Laskar & Boué 2010), $\begin{array}{}\\ {\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathit{n,m}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{n}\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}}\\ \mathrm{\times }\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{\left[}\mathrm{\right(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathit{m}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{]}}\frac{\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{!}}{\mathit{\ell }\mathrm{!}\mathrm{(}\mathit{m}\mathrm{+}\mathit{\ell }\mathrm{)}\mathrm{!}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{2}\mathrm{-}\mathrm{(}\mathit{m}\mathrm{+}\mathrm{2}\mathit{\ell }\mathrm{\left)}\mathrm{\right)}\mathrm{!}}{\left(\frac{\mathit{e}}{\mathrm{2}}\right)}^{\mathit{m}\mathrm{+}\mathrm{2}\mathit{\ell }}\\ \end{array}$(A.4)for n ≥ 2, one gets $\begin{array}{}\\ {\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathit{n,}\mathrm{0}}\mathrm{=}\mathrm{1}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}\mathit{,}\\ {\mathit{X}}_{\mathrm{0}}^{\mathrm{-}\mathit{n,}\mathrm{1}}\mathrm{=}\frac{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{2}}\mathit{e}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{e}}^{\mathrm{3}}\mathrm{\right)}\mathit{.}\\ \end{array}$(A.5)Then, the expression of ${\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}$ simplifies as $\begin{array}{}\\ {\mathrm{ℱ}}_{\mathit{n}}^{\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathit{,}{\mathit{e}}_{\mathit{p}}\mathit{,\varpi }\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathit{n}}{\mathit{f}}_{\mathit{n,}\frac{\mathit{n}}{\mathrm{2}}}{\mathit{X}}_{\mathrm{0}}^{\mathit{n,}\mathrm{0}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathrm{\right)}\\ \mathrm{+}\hspace{0.17em}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{ϵ}}_{\mathit{n}}\mathrm{)}{\mathit{f}}_{\mathit{n,}\frac{\mathit{n}\mathrm{-}\mathrm{1}}{\mathrm{2}}}{\mathit{X}}_{\mathrm{0}}^{\mathit{n,}\mathrm{1}}\mathrm{\left(}{\mathit{e}}_{\mathrm{1}}\mathrm{\right)}{\mathit{e}}_{\mathit{p}}\cos \mathit{\varpi }\\ \mathrm{+}\hspace{0.17em}\mathit{O}\mathrm{\left(}{\mathit{e}}_{\mathit{p}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}\\ \end{array}$(A.6)Substituting this expression into (A.1), one obtains $\begin{array}{}\\ \mathit{Ĥ}\mathrm{plan}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}({\mathrm{P}}_{\mathit{\alpha }}\mathrm{(}\mathit{\alpha ,}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\mathit{e}}_{\mathrm{1}}{\mathit{e}}_{\mathit{p}}{\mathrm{Q}}_{\mathit{\alpha }}\mathrm{(}\mathit{\alpha ,}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)})\\ \mathrm{+}\hspace{0.17em}\mathit{O}\mathrm{\left(}{\mathit{e}}_{\mathit{p}}^{\mathrm{2}}\mathrm{\right)}\mathit{,}\\ \end{array}$(A.7)where ${\mathrm{P}}_{\mathit{\alpha }}\mathrm{\left(}\mathit{\alpha ,}{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{q}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{\alpha }}^{\mathrm{2}\mathit{q}}{\mathit{f}}_{\mathrm{2}\mathit{q,q}}{\mathit{X}}_{\mathrm{0}}^{\mathrm{2}\mathit{q,}\mathrm{0}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{,}$(A.8)and ${\mathrm{Q}}_{\mathit{\alpha }}\mathrm{\left(}\mathit{\alpha ,}{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{e}}\sum _{\mathit{q}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\mathrm{2}\mathit{q}{\mathit{\alpha }}^{\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}}{\mathit{f}}_{\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathit{,q}}{\mathit{X}}_{\mathrm{0}}^{\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathit{,}\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{.}$(A.9)To derive explicit expressions of these two quantities, we use the analytical formulae of the coefficients ${\mathit{X}}_{\mathrm{0}}^{\mathit{n,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ for n ≥ 2 (Laskar & Boué 2010), $\begin{array}{}\\ {\mathit{X}}_{\mathrm{0}}^{\mathit{n,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{m}}\frac{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{-}\mathit{m}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\\ \mathrm{\times }\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{\left[}\mathrm{\right(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{-}\mathit{m}\mathrm{)}\mathit{/}\mathrm{2}\mathrm{]}}\frac{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{-}\mathit{m}\mathrm{)}\mathrm{!}}{\mathit{\ell }\mathrm{!}\mathrm{(}\mathit{m}\mathrm{+}\mathit{\ell }\mathrm{)}\mathrm{!}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{-}\mathit{m}\mathrm{-}\mathrm{2}\mathit{\ell }\mathrm{)}\mathrm{!}}{\left(\frac{\mathit{e}}{\mathrm{2}}\right)}^{\mathit{m}\mathrm{+}\mathrm{2}\mathit{\ell }}\mathrm{·}\\ \end{array}$(A.10)Moreover, we invert the sums and separate the cases ℓ = 0 and ℓ ≥ 1 in P_α(α,x). This gives $\begin{array}{}\\ {\mathrm{P}}_{\mathit{\alpha }}\mathrm{\left(}\mathit{\alpha ,x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{q}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{f}}_{\mathrm{2}\mathit{q,q}}{\mathit{\alpha }}^{\mathrm{2}\mathit{q}}\\ \mathrm{+}\sum _{\mathit{\ell }\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}\left(\sum _{\mathit{q}\mathrm{=}\mathit{\ell }}^{\mathrm{\infty }}{\mathit{f}}_{\mathrm{2}\mathit{q,q}}\frac{\mathrm{(}\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathit{\ell }\mathrm{!}\mathit{\ell }\mathrm{!}\mathrm{(}\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\ell }\mathrm{)}\mathrm{!}}{\mathit{\alpha }}^{\mathrm{2}\mathit{q}}\right){\left(\frac{\mathit{x}}{\mathrm{4}}\right)}^{\mathit{\ell }}\mathit{,}\\ \end{array}$(A.11)and $\begin{array}{}\\ {\mathrm{Q}}_{\mathit{\alpha }}\mathrm{\left(}\mathit{\alpha ,x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell }\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}(\sum _{\mathit{q}\mathrm{=}\mathit{\ell }}^{\mathrm{\infty }}{\mathit{f}}_{\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathit{,q}}\\ \mathrm{\times }\frac{\mathit{q}\mathrm{(}\mathrm{2}\mathit{q}\mathrm{+}\mathrm{3}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathit{\ell }\mathrm{!}\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}\mathrm{(}\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\ell }\mathrm{)}\mathrm{!}}{\mathit{\alpha }}^{\mathrm{2}\mathit{q}\mathrm{+}\mathrm{1}}){\left(\frac{\mathit{x}}{\mathrm{4}}\right)}^{\mathit{\ell }}\mathrm{·}\\ \end{array}$(A.12)It should be noted that in (A.11) and (A.12), the upper limits of ℓ and q are infinite because we consider here the infinite expansion of the perturbing function in series of α. Nevertheless, when the Hamiltonian is truncated at the order αⁿ, the sums become finite. In P_α , the maximum value taken by ℓ and q is  [n/2] , and in Q_α , it is  [(n−1)/2] .

Appendix B: Numerical coefficients of the polynomials P and Q

Here, we consider the system composed of Mercury perturbed by the seven outer planets from Venus to Neptune (p = 2,8). The secular Hamiltonian of this problem, truncated at the first order in the outer eccentricities reads as $\begin{array}{}\\ \mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathrm{=}\mathrm{-}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{\mathit{G}{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}({\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\\ \mathrm{-}{\mathit{e}}_{\mathrm{1}}{\mathit{e}}_{\mathit{p}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathit{p}}\mathrm{)})\mathit{.}\\ \end{array}$(B.1)Then, using the quasiperiodic expansion of the variables (e_pe^iϖ_p)_p = 2,8 and keeping only the terms at the fundamental frequency g₅$\mathrm{(}{\mathit{A}}_{\mathit{p}}{\mathrm{e}}^{\mathrm{i}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}}{\mathrm{)}}_{\mathit{p}\mathrm{=}\mathrm{2}\mathit{,}\mathrm{8}}$, one gets $\begin{array}{}\\ \mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathrm{=}\mathrm{-}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{\mathit{G}{\mathit{m}}_{\mathit{p}}}{{\mathit{a}}_{\mathit{p}}}({\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\\ \mathrm{-}{\mathit{e}}_{\mathrm{1}}{\mathit{A}}_{\mathit{p}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\right)\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{)})\mathit{.}\\ \end{array}$(B.2)We now define the polynomials P and Q such that $\begin{array}{}\\ \begin{array}{}\mathrm{ˇ}\\ {\mathit{H}}_{\mathit{N,}\mathrm{plan}}\end{array}\mathrm{\equiv }\mathit{Ĥ}\mathit{N,}\mathrm{plan}\mathit{/}{\mathrm{\Lambda }}_{\mathrm{1}}\\ \mathrm{=}\mathrm{-}{\mathit{g}}_{\mathrm{5}}\mathrm{P}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{g}}_{\mathrm{5}}\mathrm{Q}\mathrm{(}{{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{e}}_{\mathrm{1}}\cos \mathrm{(}{\mathit{\varpi }}_{\mathrm{1}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{5}}^{\mathit{\star }}\mathrm{)}\mathit{,}\\ \end{array}$(B.3)where ${\mathrm{\Lambda }}_{\mathrm{1}}\mathrm{=}\sqrt{\mathit{G}{\mathit{m}}_{\mathrm{0}}{\mathit{a}}_{\mathrm{1}}}$. The new polynomials P and Q are derived from P_α  and Q_α  through $\mathrm{P}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}{\mathrm{P}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,x}\right)\mathit{,}$(B.4)and $\mathrm{Q}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{1}}}{{\mathit{g}}_{\mathrm{5}}}\sum _{\mathit{p}\mathrm{=}\mathrm{2}}^{\mathrm{8}}\frac{{\mathit{m}}_{\mathit{p}}}{{\mathit{m}}_{\mathrm{0}}}\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}{\mathrm{Q}}_{\mathit{\alpha }}\left(\frac{{\mathit{a}}_{\mathrm{1}}}{{\mathit{a}}_{\mathit{p}}}\mathit{,x}\right)\mathrm{\times }{\mathit{A}}_{\mathit{p}}\mathit{.}$(B.5)Their numerical values, summarized in Table B.1, have been computed up to the order α⁵⁰ from Laskar (1990).