HTTP_Request2_Exception Unable to connect to tcp://think-ws.ca.edps.org:85. Error: php_network_getaddresses: getaddrinfo failed: Name or service not known A simple model of the chaotic eccentricity of Mercury | Astronomy & Astrophysics (A&A)
Free Access
Issue
A&A
Volume 548, December 2012
Article Number A43
Number of page(s) 11
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201219991
Published online 20 November 2012

© 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 (ϖ1) and Jupiter (ϖ5) (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 g1g5 and (g1g5)−(s1s2) secular resonances between the precession frequency of the perihelion (gi) and the regression rate of the ascending node (si). 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 g1g5 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 (a1/ap)20 in the relativistic case and (a1/ap)50 in the Newtonian case, where (ap)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 (a1/ap) in order to reach the asymptotic evolution. We also show that the maximum eccentricity attained with an initial eccentricity of e1 = 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 u0, of Mercury u1, of Venus u2, etc. The heliocentric positions of the planets are similarly noted (rp)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 up 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 as1N=r12˜2Gm0r1Gp=28mp|r1rp(t)|+p=28mpm0˜r1·p(t),\begin{equation} \hat{H}_N = \frac{\tilde{\vec r}_1^2}{2} - G \frac{m_0}{r_1} - G\sum_{p=2}^8 \frac{m_p}{\abs{\vec r_1 - \vec r_p(t)}} + \sum_{p=2}^8 \frac{m_p}{m_0}{\dop{\tilde{\vec r}_1}{\dot{\vec u}_p(t)}} , \label{eq.Htime} \end{equation}(1)where ˜r1=1\hbox{$\tilde{\vec r}_1 = \dot{\vec u}_1$} is the conjugate momentum of r1, G the gravitational constant, m0 the mass of the Sun, and (mp)p = 2,8 the masses of the other planets. The first two terms of this Hamiltonian represent the Keplerian motion and are equal to −Gm0/(2a1), where a1 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 (r1/rp)p = 2,8N=Gm02a1Gp=28mprpn=0(r1rp)nPn(r1·rpr1rp)+p=28mpm0˜r1·p,\begin{equation} \hat{H}_N = -\frac{Gm_0}{2a_1} -G\sum_{p=2}^8 \frac{m_p}{r_p} \sum_{n=0}^{\infty} \left(\frac{r_1}{r_p}\right)^n P_n \left(\frac{\dop{\vec r_1}{\vec r_p}}{r_1 r_p}\right) + \sum_{p=2}^8 \frac{m_p}{m_0}{\dop{\tilde{\vec r}_1}{\dot{\vec u}_p}} , \label{eq.HamPl} \end{equation}(2)where Pn is the Legendre polynomial of order n.

A first approximating model could stop the expansion at the second Legendre polynomial P2. 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 (Mp)p = 1,8 of all planets. In the Hamiltonian (2), the terms −Gm0/rp obtained for n = 0 do not depend on Mercury’s elements; the terms (˜r1·up)\hbox{$({\tilde{\vec r}_1} \cdot {{\vec u}_p})$} vanish when averaged over M1; the term −Gm0/(2a1) becomes constant after averaging over M1, since a1 becomes constant. As such, all these terms are dropped in the following expressions. Finally, the averaged expression of the perturbing functions ap/|r1rp| are given in Laskar & Boué (2010). The resulting Hamiltonian of the coplanar problem is then N,plan=Gp=28mpapn=2(a1ap)nn(0,0)(e1,ep,ϖ1ϖp),\begin{equation} \moy{M_p}{\hat {H}_{N,{\rm plan}}} = -G\sum_{p=2}^8 \frac{m_p}{a_p} \sum_{n=2}^\infty \left(\frac{a_1}{a_p}\right)^n {\cal F}_n^{(0,0)} (e_1, e_p, \varpi_1-\varpi_p) , \label{eq.Hsec} \end{equation}(3)where n(0,0)(e1,ep)=ϵnfn,n2X0n,0(e1)X0(n+1),0(ep)\begin{eqnarray} {\cal F}_n^{(0,0)}\!\!\!\!\!\! &&(e_1, e_p, \varpi) = \epsilon_n f_{n,\frac{n}{2}} X_0^{n,0}(e_1) X_0^{-(n+1),0} (e_p) \nonumber\\ &&\!+\! \sum_{q=0}^{[(n-1)/2]}\! 2f_{n,q} X_0^{n,n-2q}(e_1) X_0^{-(n+1),n-2q}(e_p) \cos((n-2q)\varpi). \label{eq.Fn00} \end{eqnarray}(4)In these expressions, ap, ep, and ϖp are the semimajor axis, the eccentricity, and the longitude of the pericenter of the planet p, respectively. The Xkn,m(e)\hbox{$X_k^{n,m}(e)$} are the Hansen coefficients defined by (ra)neimv=k=Xkn,m(e)eikM,\begin{equation} \left(\frac{r}{a}\right)^n \e^{{\rm i}mv} = \sum_{k=-\infty}^{\infty} X_{k}^{n,m}(e) \e^{{\rm i} k M} , \end{equation}(5)where ϵn = 0 if n is odd, and ϵn = 1 if n is even, and fn,q=(2q)!(2n2q)!22n(q!)2((nq)!)2·\begin{equation} f_{n,q} = \frac{(2q)!(2n-2q)!}{2^{2n}(q!)^2((n-q)!)^2} \cdot \end{equation}(6)To simplify the expansion of the Hamiltonian (2), we take advantage of the eccentricities (ep)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 e1. With this approximation, up to the octupole order, the averaged Hamiltonian (2) reads as N,planG8p=28mp(a12ap3(2+3e12)158a13ap4e1ep(4+3e12)cos(ϖ1ϖp)).\begin{equation} \EQM{ \moy{M_p}{\hat {H}_{N, {\rm plan}}} \approx& -\frac{G}{8}\sum_{p=2}^8 m_p \left(\frac{a_1^2}{a_p^3} (2+3e_1^2)\right. \cr & \left. -\frac{15}{8}\frac{a_1^3}{a_p^4} e_1 e_p (4+3e_1^2) \cos(\varpi_1-\varpi_p) \right) . } \end{equation}(7)More generally, as long as the Hamiltonian is truncated at the first order in ep, its expression can be written as (see Appendix A) N,plan=Gp=28mpap(Pα(a1ap,e12)e1epQα(a1ap,e12)cos(ϖ1ϖp)),\begin{equation} \EQM{ \moy{M_p}{\hat {H}_{N, {\rm plan}}} =& -G \sum_{p=2}^8 \frac{m_p}{a_p} \Bigg( \Pa\left(\frac{a_1}{a_p},e_1^2\right) \cr & - e_1 e_p \Qa\left(\frac{a_1}{a_p}, e_1^2\right) \cos (\varpi_1-\varpi_p)\Bigg) , } \label{eq.HamPQ} \end{equation}(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 x0 in Pα(α,x) is the Taylor expansion of C1(α)1=(1/2)b1/2(0)(α)\hbox{$C_1(\alpha) -1 = (1/2) b_{1/2}^{(0)}(\alpha)$} −1, and the coefficient of x1 in Pα(α,x) is the Taylor series of C3(α)/2=(1/8)αb3/2(1)(α)\hbox{$C_3(\alpha)/2 = \linebreak (1/8)\alpha b_{3/2}^{(1)}(\alpha)$}, etc., where bs(k)(α)\hbox{$b_{s}^{(k)}(\alpha)$} are Laplace coefficients (Laskar & Robutel 1995).

In the Hamiltonian (8), ep and ϖp are functions of time. Neglecting the slow diffusion of the eccentricities of the inner planets, the evolution of each zp = epexpiϖp is described well by a quasiperiodic expansion (Laskar 1990). In this study, we focus on the secular resonance g1g5. In the quasiperiodic decomposition of the variables (zp)p = 2,8, we thus keep only the terms associated to the eigenmode z5=expig5t\hbox{$z_5^\star=\exp {\rm i} g_5 t$} with frequency g5 ≈ 4.25′′/yr. The others average out and disappear from the Hamiltonian.

Table 1

Amplitude of the quadrupole εpquad\hbox{$\varepsilon_p^{\rm quad}$} (14) and octupole εpoctu\hbox{$\varepsilon_p^{\rm octu}$} (13) terms, due to each planet p, amplitude Ap, and phase ϕp of the eigenmode with frequency g5 = 4.2488163′′/yr in the quasiperiodic decomposition of zp = epexpiϖp (taken from Laskar 1990).

The amplitude and the phase of the mode z5\hbox{$z_5^\star$} in the decomposition of each variable zp 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, ϕ7 = 30.6 + 180deg, it is equivalent to take ϕ7=30.6deg\hbox{$\varphi'_7=30.6\deg$} and A7=A7\hbox{$A'_7=-A_7$}. This is the convention that is followed hereafter, but the prime is omitted for clarity. We note ϖ5=g5t+ϕ5\hbox{$\varpi_5^\star=g_5 t + \varphi_5^\star$}, where ϕ5=30.6deg\hbox{$\varphi_5^\star=30.6\deg$}. We also use the fact that Λ1=Gm0a1\hbox{$\Lambda_1=\sqrt{G m_0 a_1}$} is constant in the secular problem to rescale the Hamiltonian by this quantity as it simplifies the following expressions. We note ˇHN,plan=N,plan/Λ1.\hbox{$ \check{H}_{N,{\rm plan}} = \moy{M_p}{\hat{H}_{N,{\rm plan}}}/\Lambda_1 . $}Given that Gmp/(Λ1ap) = n1a1mp/(apm0), the resonant Hamiltonian now reads as ˇHN,plan=n1p=28mpm0a1ap(Pα(a1ap,e12)e1ApQα(a1ap,e12)cos(ϖ1ϖ5)).\begin{equation} \EQM{ \check{H}_{N,{\rm plan}} =& -n_1 \sum_{p=2}^8 \frac{m_p}{m_0}\frac{a_1}{a_p} \left( \Pa\left(\frac{a_1}{a_p},e_1^2\right) \right. \cr & \left. - e_1 A_p \Qa\left(\frac{a_1}{a_p}, e_1^2\right) \cos (\varpi_1-\varpi_5^\star)\right) . } \label{eq.Hamres0} \end{equation}(9)With this convention, the time should also be rescaled by the same factor Λ1 to keep the canonical form of the equations of motion, unless the canonical variables are modified as follows. We let =Λ1(11e12),θ̂=ϖ1\begin{equation} \EQM{ \hat I = \Lambda_1 \left(1-\sqrt{1-e_1^2}\right) , \crm \hat \theta = -\varpi_1 } \end{equation}(10)be the canonical conjugated variables of the initial Hamiltonian ĤN,plan of (Eq. (8)). It can be easily shown that (ˇI,ˇθ)=(Λ1,θ̂)\begin{equation} (\check{I},\check{\theta}) = \left(\frac{\hat I}{\Lambda_1}, \hat \theta\right) \end{equation}(11)are canonical conjugated variables of the Hamiltonian ˇHN,plan\hbox{$\check{H}_{N,{\rm plan}}$} (Eq. (9)) without rescaling the time t. The equations of motion are thus dˇIdt=ˇHN,planˇθ(ˇI,ˇθ,t),dˇθdt=ˇHN,planˇI(ˇI,ˇθ,t).\begin{equation} \EQM{ \frac{{\rm d}\check{I}}{{\rm d}t} &= -\Dron{\check{H}_{N,{\rm plan}}}{\check \theta} (\check I,\check \theta,t) , \qquad & \frac{{\rm d}\check \theta}{{\rm d}t} &= \Dron{\check{H}_{N,{\rm plan}}}{\check I} (\check I,\check \theta,t) . } \end{equation}(12)At this point, it is interesting to evaluate the contribution of each planet to the octupole interaction with Mercury. Although the resonance g1g5 involves only the precession frequencies associated to Mercury and Jupiter, the eigenmode z5\hbox{$z_5^\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 (εpoctu)p=2,8\hbox{$(\varepsilon^{\rm octu}_p)_{p=2,8}$} given by εpoctu=1564(n1g5)(mpm0)(a1ap)4Ap.\begin{equation} \varepsilon_p^{\rm octu} = \frac{15}{64}\left(\frac{n_1}{g_5}\right)\left(\frac{m_p}{m_0}\right) \left(\frac{a_1}{a_p}\right)^4 A_p. \label{eq.eps_octu} \end{equation}(13)The values taken by these parameters are gathered in Table 1. The maximal amplitude is due to Venus with ε2octu0.0012\hbox{$\varepsilon_2^{\rm octu}\approx 0.0012$} followed by the Earth-Moon barycenter and Jupiter with ε3octuε5octu4×10-4\hbox{$\varepsilon_3^{\rm octu}\approx \varepsilon_5^{\rm octu}\approx 4\times10^{-4}$}. Thus, the strongest perturbation on Mercury’s orbit comes from the precession of Venus excited by Jupiter. a1/a2 ≈ 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 εpquad=18(n1g5)(mpm0)(a1ap)3·\begin{equation} \varepsilon_p^{\rm quad} = \frac{1}{8} \left(\frac{n_1}{g_5}\right)\left(\frac{m_p}{m_0}\right) \left(\frac{a_1}{a_p}\right)^3 \cdot \label{eq.eps_quad} \end{equation}(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 e1 alone, given by (Eq. (9)) P(x)=n1g5p=28mpm0a1apPα(a1ap,x),Q(x)=n1g5p=28mpm0a1apQα(a1ap,x)×Ap.\begin{equation} \EQM{ \Pb(x) &= \frac{n_1}{g_5} \sum_{p=2}^8 \frac{m_p}{m_0} \frac{a_1}{a_p} \Pa\left(\frac{a_1}{a_p}, x\right) , \crm \Qb(x) &= \frac{n_1}{g_5} \sum_{p=2}^8 \frac{m_p}{m_0} \frac{a_1}{a_p} \Qa\left(\frac{a_1}{a_p}, x\right) \times A_p . } \end{equation}(15)Then, the resonant Hamiltonian (9) reads as ˇHN,plan=g5(P(e12)e1Q(e12)cos(ϖ1ϖ5)).\begin{equation} \check{H}_{N,{\rm plan}} = -g_5 \left(\Pb(e_1^2) - e_1 \Qb(e_1^2) \cos(\varpi_1-\varpi_5^\star) \right) . \label{eq.Hamres} \end{equation}(16)The numerical values of the coefficients of the polynomials P and Q are given in Appendix B.

The Hamiltonian ˇHN,plan\hbox{$\check{H}_{N,{\rm plan}}$} (16) has one and a half degrees of freedom, but it can be reduced to a one degree of freedom Hamiltonian ˜HN,plan\hbox{$\tilde{H}_{N,{\rm plan}}$} after some modifications.For that, we first make the Hamiltonian autonomous by a adding a momentum ˇT\hbox{$\check T$} conjugated to time t. Then, we perform a canonical transformation where the old variables are (ˇI=11e12,ˇθ=ϖ1);(ˇT,t),\begin{equation} (\check I = 1-\sqrt{1-e_1^2}\ , \check \theta = -\varpi_1) ; \quad (\check T ,t) , \end{equation}(17)and the new ones are (I,Δϖ) and (˜T,˜t)\hbox{$(\tilde T, \tilde t)$}, defined by Δϖ=ˇθg5tϕ5ϖ1ϖ5,˜t=t.\begin{equation} \EQM{ \Delta\varpi = -\check \theta - g_5 t - \varphi_5^\star \equiv \varpi_1 - \varpi_5^\star , \\ \tilde t = t . } \end{equation}(18)For this transformation to be canonical, the momentums should verify ˇI=I11e12,ˇT=g5I+˜T.\begin{equation} \EQM{ \check I = - I \equiv 1-\sqrt{1-e_1^2} ,\\ \check T = -g_5 I + \tilde T . } \end{equation}(19)The new Hamiltonian ˜HN,plan\hbox{$\tilde{H}_{N,{\rm plan}}$} expressed in the new variables does not depend on the cyclic coordinate ˜t\hbox{$\tilde t$}. Thus, its conjugated momentum ˜T\hbox{$\tilde T$} is an integral of the motion and can be dropped. Up to a constant, the new Hamiltonian is then˜HN,plan=g51e12g5P(e12)+g5e1Q(e12)cosΔϖ.\begin{equation} \tilde{H}_{N,{\rm plan}} = -g_5\sqrt{1-e_1^2} -g_5 \Pb\left(e_1^2\right) + g_5 e_1 \Qb\left(e_1^2\right) \cos \Delta\varpi . \label{eq.Ham1} \end{equation}(20)This one degree of freedom Hamiltonian is integrable. The orbits are given by ˜HN,plan=Cte\hbox{$\tilde{H}_{N,{\rm plan}}={\rm Cte}$}. The temporal evolution of the eccentricity e1 and of the resonant angle Δϖ are deduced from the canonical equations of motion. This leads to de1dt=1e12e1˜HN,planΔϖ,ϖdt=1e12e1˜HN,plane1·\begin{equation} \EQM{ \frac{{\rm d}e_1}{{\rm d}t} &= \frac{\sqrt{1-e_1^2}}{e_1}\Dron{\tilde{H}_{N,{\rm plan}}}{\Delta\varpi}\ ,\crm \frac{{\rm d}\Delta\varpi}{{\rm d}t} &=-\frac{\sqrt{1-e_1^2}}{e_1}\Dron{\tilde{H}_{N,{\rm plan}}}{e_1} \cdot } \label{eq.dt} \end{equation}(21)

2.2. General relativistic precession

The secular effect of relativity is described by (e.g., Touma et al. 2009), HR=gr11e12,\begin{equation} H_R = -g_r \frac{1}{\sqrt{1-e_1^2}} , \end{equation}(22)where gr=3(Gm0)2/(Λ1a12c2)\hbox{$g_r = 3 (G m_0)^2/(\Lambda_1 a_1^2 c^2)$}, and c is the speed of light. In the case of Mercury, we have gr ≈ 0.41′′/yr. The total Hamiltonian ˜HR,plan=˜HN,plan+HR\hbox{$\tilde{H}_{R,{\rm plan}}=\tilde{H}_{N,{\rm plan}}+H_R$}, including the Newtonian interaction and general relativity, becomes ˜HR,plan=g5(1e12+P(e12)e1Q(e12)cosΔϖ)gr11e12·\begin{equation} \EQM{ \tilde{H}_{R,{\rm plan}} = & -g_5\left(\sqrt{1-e_1^2} + \Pb\left(e_1^2\right) - e_1 \Qb\left(e_1^2\right) \cos \Delta\varpi \right) \crm & -g_r \frac{1}{\sqrt{1-e_1^2}}\cdot } \label{eq.Hamtot} \end{equation}(23)

thumbnail 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.

thumbnail 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 g1g5 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 g1 = 5.54′′/yr in the Newtonian case and g1 = 5.96′′/yr in the relativistic case, whereas the resonance occurs in the vicinity of g1 = g5 = 4.25′′/yr. This result is slightly different, but representative of the present behavior of Mercury’s eccentricity. Indeed, the present value of g1 for the solar system with GR is g1 = 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 g1 by an increment δg. In addition, owing to the slow chaotic diffusion in the inner solar system, Mercury’s precession frequency g1 can change quasi-randomly and come close to the value of g5, 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 g1 around the resonant frequency g5, the above-mentioned correction is added to both ˜HN,plan\hbox{$\tilde{H}_{N,{\rm plan}}$} and ˜HR,plan\hbox{$\tilde{H}_{R,{\rm plan}}$}; i.e., HN,plan=˜HN,plan+12δge12\begin{equation} H_{N,{\rm plan}} = \tilde{H}_{N,{\rm plan}} + \frac{1}{2}\delta_g e_1^2 \label{eq.HamNplan} \end{equation}(24)and HR,plan=˜HR,plan+12δge12.\begin{equation} H_{R,{\rm plan}} = \tilde{H}_{R,{\rm plan}} + \frac{1}{2}\delta_g e_1^2 . \label{eq.HamRplan} \end{equation}(25)In Eqs. (24) and (25), the factor δg controls the precession frequency g1. For δg = 0, we recover the values obtained with ˜HN,plan\hbox{$\tilde{H}_{N,{\rm plan}}$} and ˜HR,plan\hbox{$\tilde{H}_{R,{\rm plan}}$}, respectively. Otherwise the value of g1 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 HN,plan (24) and HR,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 e1, but to show that e1 can actually reach values as high as observed by Laskar (2008) (see Fig. 3a) when the system is in the vicinity of the resonance g1 = g5. 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 (a1/ap)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., HN,plan(n)\hbox{$H_{N,{\rm plan}}^{(n)}$} or HR,plan(n)\hbox{$H_{R,{\rm plan}}^{(n)}$}. In a first step, we focus only on the relativistic Hamiltonian which is more complete.

Within the quadrupole approximation, the Hamiltonian HR,plan(2)\hbox{$H_{R,{\rm plan}}^{(2)}$} reduces to a simple expression HR,plan(2)=g5(1e12+εquad(2+3e12))gr11e12+12δge12,\begin{equation} \EQM{ H_{R,{\rm plan}}^{(2)} = & -g_5 \left(\sqrt{1-e_1^2} +\varepsilon^{\rm quad} \left(2+3e_1^2\right)\right) \crm & - g_r \frac{1}{\sqrt{1-e_1^2}} + \frac{1}{2}\delta_g e_1^2 , } \label{eq.Ham2} \end{equation}(26)where εquad=pεpquad\hbox{$\varepsilon^{\rm quad} = \sum_p \varepsilon^{\rm quad}_p$}. This Hamiltonian is independent of Δϖ, thus e1 remains constant. All the trajectories in the plane (x = e1cosΔϖ, y = e1sinΔϖ) are circles centered on the origin (0,0) with an eventually infinite period when \hbox{$\Delta\dot\varpi=0$}. Figure 1a illustrates this result for g1 = 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 ϖ1 and ϖ5\hbox{$\varpi_5^\star$} and thus, no possible increase in the eccentricity.

The next level of approximation is the octupole. The corresponding Hamiltonian reads as HR,plan(3)=g5(1e12+εquad(2+3e12)εoctu(4+3e12)cosΔϖ)gr11e12+δge12,\begin{equation} \EQM{ H_{R,{\rm plan}}^{(3)} =& -g_5 \left( \sqrt{1-e_1^2} + \varepsilon^{\rm quad} (2+3e_1^2)\right. \cr & - \varepsilon^{\rm octu} (4+3e_1^2) \cos\Delta\varpi \bigg) -g_r \frac{1}{\sqrt{1-e_1^2}} + \delta_g e_1^2 , } \end{equation}(27)with εoctu=pεpoctu\hbox{$\varepsilon^{\rm octu} = \sum_p \varepsilon^{\rm octu}_p$}. 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 g1 = 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 g1. This is depicted in the subpanel labeled α3 of Fig. 2. The curves represent the position on the x axis of the fixed points of the phase space as a function of g1. 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 HR,plan(n)\hbox{$H_{R,{\rm plan}}^{(n)}$} increases, the degrees of the polynomials Pα(α,e2) and Qα(α,e2) increase both for α and e2. Then, the topology of the phase space of the Hamiltonian HR,plan(n)\hbox{$H_{R,{\rm plan}}^{(n)}$} evolves as shown in Fig. 2. Up to n = 5, five fix points – three stable points and two unstable points – coexist within 4.5 ≲ g1 ≲ 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 e1cosΔϖ > 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 g1 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 e1 = 0.8. The separatrix disappears at g1 ≥ 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 g1 ≲ 4.5′′/yr. For example, with g1 = 4.0′′/yr (Fig. 1e), a trajectory starting within e1 ≤ 0.2 can reach e1 ≈ 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 e1 ≈ 0.2 up to e1 ≈ 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 HR,plan(20)\hbox{$H_{R,{\rm plan}}^{(20)}$} 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 dxdt=1e12HR,plan(20)∂y;dydt=1e12HR,plan(20)∂x,\begin{equation} \EQM{ \frac{{\rm d}x}{{\rm d}t} = &\sqrt{1-e_1^2}\, \Dron{H_{R,{\rm plan}}^{(20)}}{y} ;\crm \frac{{\rm d}y}{{\rm d}t} =-&\sqrt{1-e_1^2}\, \Dron{H_{R,{\rm plan}}^{(20)}}{x} , } \label{eq.integ} \end{equation}(28)where x = e1cosΔϖ and y = e1sinΔϖ. They are equivalent to the equations of motion (21), but without singularities at the origin e1 = 0.

thumbnail Fig. 3

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

thumbnail 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 e1 = 0.2, Δϖ = 110deg, and g1 = 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.

thumbnail 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 HN,plan(n)\hbox{$H_{N,{\rm plan}}^{(n)}$} (24). Figure 4 shows the level curves of HN,plan(20)\hbox{$H_{N,{\rm plan}}^{(20)}$} and HN,plan(50)\hbox{$H_{N,{\rm plan}}^{(50)}$}. The two figures are similar to the subpanel α20 of Fig. 1 within the region of low eccentricities (e1 ≲ 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 g1, it should be increased by gr = 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 HN,plan(50)\hbox{$H_{N,{\rm plan}}^{(50)}$}. 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 g1g5 resonance.

To do so, we add inclination in our system, and focus on the term related to (g1g5)−(s1s2). 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) ˇHR,inc=gr11e12n1p=28(mpm0(a1ap)n+1×n(0,0)(e1,Ap,i1,Ip,ϖ1ϖ5,Ω1Ω2)),\begin{equation} \EQM{ \check{H}_{R,{\rm inc}} =& - g_r\frac{1}{\sqrt{1-e_1^2}} -n_1 \sum_{p=2}^8 \left(\frac{m_p}{m_0}\left(\frac{a_1}{a_p}\right)^{n+1}\right. \crm & \left. \hspace*{-1.1cm} \phantom{\left(\frac{a_1}{a_p}\right)^{n+1}} \times\, {\cal F}_{n}^{(0,0)} \left(e_1,A_p,i_1,I_p,\varpi_1-\varpi_5^\star,\Omega_1-\Omega_2^\star\right) \right), } \label{eq.Haminc0} \end{equation}(29)where ip and Ωp are the inclination and the longitude of the ascending node of the planet p, and Ip and Ω2\hbox{$\Omega_2^\star$} are the amplitude and the argument of the term precessing at the frequency s2 in the quasiperiodic decomposition of sin(ip/2)exp(iΩp), respectively. In (29), both ϖ5=g5t+ϕ5\hbox{$\varpi_5^\star=g_5 t+\varphi_5^\star$} and Ω2=s2t+φ2\hbox{$\Omega_2^\star=s_2 t+\phi_2^\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 ˇT\hbox{$\check T$} conjugated to the time t. The old variables of the subsequent canonical transformation are ˇI=11e12,ˇθ=ϖ1,ˇJ=1e12(1cosi1),ˇψ=Ω1,ˇT,t,\begin{equation} \EQM{ \check I = 1-\sqrt{1-e_1^2} ,\quad & \check \theta = -\varpi_1 , \\ \check J = \sqrt{1-e_1^2}(1-\cos i_1) ,\quad & \check \psi = -\Omega_1 , \\ \check T\ , & t , } \label{eq.oldinc} \end{equation}(30)and the new variables are denoted (I,Δϖ), (J,ΔΩ), (˜T,˜t)\hbox{$(\tilde T, \tilde t)$}. We set Δϖ=ˇθg5tϕ5ϖ1ϖ5,ΔΩ=ˇψs2tφ2Ω1Ω2,˜t=t.\begin{equation} \EQM{ \Delta\varpi = -\check \theta - g_5 t - \varphi_5^\star \equiv \varpi_1-\varpi_5^\star ,\crm \Delta\Omega = -\check \psi - s_2 t - \phi_2^\star \equiv \Omega_1 - \Omega_2^\star ,\crm \tilde t = t . } \end{equation}(31)The associated conjugated momenta are given by ˇI=I,ˇJ=J,ˇT=g5Is2J+˜T.\begin{equation} \EQM{ \check I = -I ,\crm \check J = -J ,\crm \check T = -g_5 I -s_2 J + \tilde T . } \end{equation}(32)Since the new Hamiltonian ˜HR,inc\hbox{$\tilde H_{R, {\rm inc}}$} is independent of t, its conjugate momentum ˜T\hbox{$\tilde T$} is an integral of the motion. From the expressions of ˇI\hbox{$\check I$} and ˇJ\hbox{$\check J$} (30), we thus get, up to a constant, ˜HR,inc=(g5+2s2sin2i12)1e12gr11e12n1p=28mpm0(a1ap)n+1n(0,0)(e1,Ap,i1,Ip,Δϖ,ΔΩ).\begin{equation} \EQM{ \tilde{H}_{R,{\rm inc}} =& \left(-g_5+2s_2\sin^2\frac{i_1}{2}\right)\sqrt{1-e_1^2} -g_r \frac{1}{\sqrt{1-e_1^2}} \crm & -n_1\sum_{p=2}^8 \frac{m_p}{m_0} \left(\frac{a_1}{a_p}\right)^{n+1} {\cal F}_n^{(0,0)}\left(e_1,A_p,i_1,I_p,\Delta\varpi,\Delta\Omega\right) . } \end{equation}(33)To obtain the final Hamiltonian HR,inc in the spatial case, we add the control term and get HR,inc=˜HR,inc+12δge12.\begin{equation} H_{R,{\rm inc}} = \tilde{H}_{R,{\rm inc}} + \frac{1}{2}\delta_g e_1^2 . \label{eq.Haminc} \end{equation}(34)The Hamiltonian (34) now has two degrees of freedom. The variables are (e1ϖ) and (i1,ΔΩ). To perform numerical integrations, we use nonsingular rectangular coordinates k=e1cosΔϖ,h=e1sinΔϖq=sini12cosΔΩ,p=sini12sinΔΩ.\begin{equation} \EQM{ &k=e_1\cos\Delta\varpi , \quad &h=e_1\sin\Delta\varpi \crm &q=\sin\frac{i_1}{2}\cos\Delta\Omega , \quad &p=\sin\frac{i_1}{2}\sin\Delta\Omega . } \end{equation}(35)With χ=1e12\hbox{$\chi=\sqrt{1-e_1^2}$}, the equations of motion read as (e.g. Bretagnon 1974) dkdt=+χHinc∂h+h2χ(qHinc∂q+pHinc∂p),dhdt=χHinc∂kk2χ(qHinc∂q+pHinc∂p),dqdt=+14χHinc∂pq2χ(hHinc∂kkHinc∂h),dpdt=14χHinc∂qp2χ(hHinc∂kkHinc∂h)·\begin{equation} \EQM{ \frac{{\rm d}k}{{\rm d}t} &= +\chi\Dron{H_{\rm inc}}{h} + \frac{h}{2\chi} \left(q\Dron{H_{\rm inc}}{q}+p\Dron{H_{\rm inc}}{p}\right) , \crm \frac{{\rm d}h}{{\rm d}t} &= -\chi\Dron{H_{\rm inc}}{k} - \frac{k}{2\chi} \left(q\Dron{H_{\rm inc}}{q}+p\Dron{H_{\rm inc}}{p}\right) , \crm \frac{{\rm d}q}{{\rm d}t} &= +\frac{1}{4\chi}\Dron{H_{\rm inc}}{p} -\frac{q}{2\chi} \left(h\Dron{H_{\rm inc}}{k}-k\Dron{H_{\rm inc}}{h}\right) , \crm \frac{{\rm d}p}{{\rm d}t} &= -\frac{1}{4\chi}\Dron{H_{\rm inc}}{q} -\frac{p}{2\chi} \left(h\Dron{H_{\rm inc}}{k}-k\Dron{H_{\rm inc}}{h}\right) \cdot } \end{equation}(36)Since Venus is the only planet with a significative amplitude associated to the frequency s2, 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 (g1g5)−(s1s2)), in a state of slow chaos, with slow chaotic diffusion. This diffusion will quasi-randomly change the value of g1, which can then approach resonance with g5. 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 g1g5 resonance.

thumbnail 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 g1 = g5. 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 (g1g5)−(s1s2) 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).


1

For clarity, we drop the explicit time dependence in all equations after (1).

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

  1. Batygin, K., & Laughlin, G. 2008, ApJ, 683, 1207 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bretagnon, P. 1974, A&A, 30, 141 [NASA ADS] [Google Scholar]
  3. Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189 [NASA ADS] [CrossRef] [Google Scholar]
  4. Kozai, Y. 1962, AJ, 67, 579 [NASA ADS] [CrossRef] [Google Scholar]
  5. Laskar, J. 1984, Thèse, Observatoire de Paris http://tel.archives-ouvertes.fr/tel-00702723 [Google Scholar]
  6. Laskar, J. 1989, Nature, 338, 237 [NASA ADS] [CrossRef] [Google Scholar]
  7. Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
  8. Laskar, J. 1994, A&A, 287, L9 [NASA ADS] [Google Scholar]
  9. Laskar, J. 2008, Icarus, 196, 1 [NASA ADS] [CrossRef] [Google Scholar]
  10. Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  12. Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  13. Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Lidov, M. 1962, Planet. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
  15. Lidov, M. L., & Ziglin, S. L. 1976, Celest. Mech., 13, 471 [NASA ADS] [CrossRef] [Google Scholar]
  16. Lithwick, Y., & Wu, Y. 2011, ApJ, 739, 31 [NASA ADS] [CrossRef] [Google Scholar]
  17. Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 [NASA ADS] [CrossRef] [Google Scholar]

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 plan=Gmpapn=2αnn(0,0)(e1,ep,ϖ1ϖp).\appendix \setcounter{section}{1} \begin{equation} \hat{H}_{\rm plan} = -\frac{G m_p}{a_p} \sum_{n=2}^\infty \alpha^n {\cal F}_n^{(0,0)} (e_1, e_p, \varpi_1-\varpi_p) . \label{eq.AHam1} \end{equation}(A.1)The indices 1 and p refer to the inner massless body and to the outer massive planet, respectively; mp is the mass of the outer planet; α = a1/ap is the semimajor axis ratio; (ek)k = 1,p and (ϖk)k = 1,p are the eccentricities and the longitudes of periastron of the two planets.

In (A.1), the n(0,0)(e1,ep)\hbox{${\cal F}_n^{(0,0)} (e_1, e_p, \varpi)$} are given by Laskar & Boué (2010)n(0,0)(e1,ep)=ϵnfn,n2X0n,0(e1)X0(n+1),0(ep)+q=0[(n1)/2]2fn,qX0n,n2q(e1)X0(n+1),n2q(ep)cos((n2q)ϖ),\appendix \setcounter{section}{1} \begin{equation} \EQM{ {\cal F}_n^{(0,0)} (e_1, e_p, \varpi) = \epsilon_n f_{n,\frac{n}{2}} X_0^{n,0}(e_1) X_0^{-(n+1),0}(e_p) \crm +\sum_{q=0}^{[(n-1)/2]} 2 f_{n,q} X_0^{n,n-2q}(e_1) X_0^{-(n+1),n-2q}(e_p) \cos((n-2q)\varpi) , } \label{eq.AF00} \end{equation}(A.2)where ϵn = 1 if n is even and 0, otherwise fn,q=(2q)!(2n2q)!22n(q!)2((nq)!)2,\appendix \setcounter{section}{1} \begin{equation} f_{n,q} = \frac{(2q)!(2n-2q)!}{2^{2n}(q!)^2((n-q)!)^2} , \label{eq.fnq} \end{equation}(A.3)and X0n,m(e)\hbox{$X_0^{n,m}(e)$} are Hansen coefficients.

Now, we assume that the eccentricity of the outer planet remains low, and we only keep the linear terms in ep. Since the lowest power in eccentricity of X0(n+1),m(e)\hbox{$X_0^{-(n+1),m}(e)$} 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), X0n,m=1(1e2)n3/2×=0[(n2m)/2](n2)!!(m+)!(n2(m+2))!(e2)m+2\appendix \setcounter{section}{1} \begin{equation} \EQM{ X_0^{-n,m} = & \frac{1}{(1-e^2)^{n-3/2}} \crm & \times \sum_{\ell=0}^{[(n-2-m)/2]}\!\! \!\!\frac{(n-2)!}{\ell!(m+\ell)!(n-2-(m+2\ell))!} \left(\frac{e}{2}\right)^{m+2\ell} } \end{equation}(A.4)for n ≥ 2, one gets X0n,0=1+O(e2),X0n,1=n22e+O(e3).\appendix \setcounter{section}{1} \begin{equation} \EQM{ X_0^{-n,0} = 1 + O(e^2)\ , \crm X_0^{-n,1} = \frac{n-2}{2} e + O(e^3) . } \label{eq.ha} \end{equation}(A.5)Then, the expression of n(0,0)\hbox{${\cal F}_n^{(0,0)}$} simplifies as n(0,0)(e1,ep)=ϵnfn,n2X0n,0(e1)+(n1)(1ϵn)fn,n12X0n,1(e1)epcosϖ+O(ep2).\appendix \setcounter{section}{1} \begin{equation} \EQM{ {\cal F}_n^{(0,0)} (e_1, e_p, \varpi) = & \epsilon_n f_{n,\frac{n}{2}} X_0^{n,0}(e_1) \crm & +\,(n-1)(1-\epsilon_n) f_{n,\frac{n-1}{2}} X_0^{n,1}(e_1) e_p \cos\varpi \crm & +\, O(e_p^2) . } \end{equation}(A.6)Substituting this expression into (A.1), one obtains plan=Gmpap(Pα(α,e12)e1epQα(α,e12)cos(ϖ1ϖp))+O(ep2),\appendix \setcounter{section}{1} \begin{equation} \EQM{ \hat{H}_{\rm plan} = & -\frac{G m_p}{a_p} \Bigg( \Pa\left(\alpha, e_1^2\right) - e_1 e_p \Qa\left(\alpha, e_1^2\right) \cos (\varpi_1-\varpi_p) \Bigg) \crm & +\, O(e_p^2) , } \label{eq.AHam3} \end{equation}(A.7)where Pα(α,e2)=q=1α2qf2q,qX02q,0(e),\appendix \setcounter{section}{1} \begin{equation} \Pa\left(\alpha, e^2 \right) = \sum_{q=1}^\infty \alpha^{2q} f_{2q,q} X_0^{2q,0} (e) , \end{equation}(A.8)and Qα(α,e2)=1eq=12qα2q+1f2q+1,qX02q+1,1(e).\appendix \setcounter{section}{1} \begin{equation} \Qa\left(\alpha, e^2 \right) = -\frac{1}{e}\sum_{q=1}^{\infty} 2q \alpha^{2q+1} f_{2q+1,q} X_0^{2q+1,1}(e) . \end{equation}(A.9)To derive explicit expressions of these two quantities, we use the analytical formulae of the coefficients X0n,m(e)\hbox{$X_0^{n,m}(e)$} for n ≥ 2 (Laskar & Boué 2010), X0n,m(e)=(1)m(n+1m)!(n+1)!×=0[(n+1m)/2](n+1m)!!(m+)!(n+1m2)!(e2)m+2·\appendix \setcounter{section}{1} \begin{equation} \EQM{ X_0^{n,m} (e) = & (-1)^m \frac{(n+1-m)!}{(n+1)!} \crm & \times \sum_{\ell=0}^{[(n+1-m)/2]} \frac{(n+1-m)!}{\ell!(m\! +\! \ell)!(n\! +\! 1\! -\! m\! -\! 2\ell)!} \left(\frac{e}{2}\right)^{m+2\ell} \cdot } \end{equation}(A.10)Moreover, we invert the sums and separate the cases  = 0 and  ≥ 1 in Pα(α,x). This gives Pα(α,x)=q=1f2q,qα2q+=1(q=f2q,q(2q+1)!!!(2q+12)!α2q)(x4),\appendix \setcounter{section}{1} \begin{equation} \EQM{ \Pa(\alpha, x) & = \sum_{q=1}^{\infty} f_{2q,q} \alpha^{2q} \crm & + \sum_{\ell=1}^{\infty} \left( \sum_{q=\ell}^\infty f_{2q,q} \frac{(2q+1)!}{\ell!\ell!(2q+1-2\ell)!} \alpha^{2q} \right) \left(\frac{x}{4}\right)^\ell , } \label{eq.Pa} \end{equation}(A.11)and Qα(α,x)==0(q=f2q+1,q×q(2q+3)(2q+1)!!(+1)!(2q+12)!α2q+1)(x4)·\appendix \setcounter{section}{1} \begin{equation} \EQM{ \Qa(\alpha, x) = \sum_{\ell=0}^{\infty} \Bigg( & \sum_{q=\ell}^\infty f_{2q+1,q} \crm & \times\frac{q(2q+3) (2q+1)!}{\ell!(\ell+1)!(2q+1-2\ell)!} \alpha^{2q+1} \Bigg) \left(\frac{x}{4}\right)^\ell \cdot } \label{eq.Qa} \end{equation}(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 αn, 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

Table B.1

Coefficients of the polynomials P(x) =  ∑ px (B.4) and Q(x) =  ∑ qx (B.5) computed up to α50.

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 N,plan=p=28Gmpap(Pα(a1ap,e12)e1epQα(a1ap,e12)cos(ϖ1ϖp)).\appendix \setcounter{section}{2} \begin{equation} \EQM{ \moy{M_p}{\hat{H}_{N,{\rm plan}}} =& -\sum_{p=2}^8 \frac{G m_p}{a_p} \Bigg( \Pa\left(\frac{a_1}{a_p}, e_1^2\right) \crm & - e_1 e_p \Qa\left(\frac{a_1}{a_p}, e_1^2\right) \cos (\varpi_1-\varpi_p) \Bigg) . } \label{eq.BHam1} \end{equation}(B.1)Then, using the quasiperiodic expansion of the variables (epeiϖp)p = 2,8 and keeping only the terms at the fundamental frequency g5(Apeiϖ5)p=2,8\hbox{$(A_p\e^{{\rm i}\varpi^\star_5})_{p=2,8}$}, one gets N,plan=p=28Gmpap(Pα(a1ap,e12)e1ApQα(a1ap,e12)cos(ϖ1ϖ5)).\appendix \setcounter{section}{2} \begin{equation} \EQM{ \moy{M_p}{\hat{H}_{N,{\rm plan}}} = &-\sum_{p=2}^8 \frac{G m_p}{a_p} \Bigg( \Pa\left(\frac{a_1}{a_p}, e_1^2\right) \crm & - e_1 A_p \Qa\left(\frac{a_1}{a_p}, e_1^2\right) \cos (\varpi_1-\varpi^\star_5) \Bigg) . } \label{eq.BHam2} \end{equation}(B.2)We now define the polynomials P and Q such that ˇHN,planN,plan/Λ1=g5P(e12)+g5Q(e12)e1cos(ϖ1ϖ5),\appendix \setcounter{section}{2} \begin{equation} \EQM{ \check{H}_{N,{\rm plan}} & \equiv \moy{M_p}{\hat{H}_{N,{\rm plan}}}/\Lambda_1 \crm & = -g_5 \Pb\left(e_1^2\right) + g_5 \Qb\left(e_1^2\right) e_1 \cos(\varpi_1-\varpi^\star_5) , } \label{eq.BHam3} \end{equation}(B.3)where Λ1=Gm0a1\hbox{$\Lambda_1 = \sqrt{G m_0 a_1}$}. The new polynomials P and Q are derived from Pα  and Qα  through P(x)=n1g5p=28mpm0a1apPα(a1ap,x),\appendix \setcounter{section}{2} \begin{equation} \Pb(x) = \frac{n_1}{g_5} \sum_{p=2}^8 \frac{m_p}{m_0} \frac{a_1}{a_p} \Pa\left(\frac{a_1}{a_p}, x\right) , \label{eq.Pbnum} \end{equation}(B.4)and Q(x)=n1g5p=28mpm0a1apQα(a1ap,x)×Ap.\appendix \setcounter{section}{2} \begin{equation} \Qb(x) = \frac{n_1}{g_5} \sum_{p=2}^8 \frac{m_p}{m_0} \frac{a_1}{a_p} \Qa\left(\frac{a_1}{a_p}, x\right) \times A_p . \label{eq.Qbnum} \end{equation}(B.5)Their numerical values, summarized in Table B.1, have been computed up to the order α50 from Laskar (1990).

Appendix C: Frequencies at zero eccentricity and corrections

The precession frequency g1 of Mercury’s perihelia computed either from the Newtonian Hamiltonian (20), ˜HN,plan=g51e12g5P(e12)+g5e1Q(e12)cosΔϖ,\appendix \setcounter{section}{3} \begin{equation} \tilde{H}_{N,{\rm plan}} = -g_5\sqrt{1-e_1^2} -g_5 \Pb\left(e_1^2\right) + g_5 e_1 \Qb\left(e_1^2\right) \cos \Delta\varpi , \end{equation}(C.1)or from the relativistic Hamiltonian (23) ˜HR,plan=g5(1e12+P(e12)e1Q(e12)cosΔϖ)gr11e12\appendix \setcounter{section}{3} \begin{equation} \EQM{ \tilde{H}_{R,{\rm plan}} = & -g_5\left(\sqrt{1-e_1^2} + \Pb\left(e_1^2\right) - e_1 \Qb\left(e_1^2\right) \cos \Delta\varpi \right) \crm & -g_r \frac{1}{\sqrt{1-e_1^2}} } \end{equation}(C.2)is far enough from g5 ≈ 4.249′′/yr for the system not to be in secular resonance (see Table C.1). Nevertheless, due to the slow diffusion of the inner planets of the solar system, this frequency is subject to small variations, and Mercury can eventually reach the resonance. To model this change in frequency, we include an additional term in the Hamiltonians: δge12/2\hbox{$\delta_g e_1^2/2$}.

Table C.1

Frequencies at zero eccentricity and corrections for different orders of expansion of the perturbing function.

The values of δg putting the system in exact resonance at zero eccentricity are shown in Table C.1. One can observe that the correction is higher (in absolute value) when the relativistic precession is taken into account, which explains why relativity stabilizes the solar system.

Appendix D: Explicit expression of the Hamiltonian development in the spatial case

Here we develop the secular Hamiltonian of an inclined restricted two-planet system following Laskar & Boué (2010). We use the same kind of approximation as for the planar case by considering only the linear dependency on the eccentricities and inclinations of the outer planets. We note i1 and ip the inclinations of the massless planet and of the perturber respectively, and Ω1 and Ωp are their longitudes of ascending node. We also note c1 = cos(i1/2), s1 = sin(i1/2), cp = cos(ip/2), and sp = sin(ip/2), from which we define μ=(c1cpeiΩ1Ωp2+s1speiΩ1Ωp2)2,ν=(c1speiΩ1Ωp2s1cpeiΩ1Ωp2)2.\appendix \setcounter{section}{4} \begin{equation} \EQM{ \mu_* &= \left(c_1c_p\e^{\ii\frac{\Omega_1-\Omega_p}{2}} +s_1s_p\e^{-\ii\frac{\Omega_1-\Omega_p}{2}}\right)^2 , \crm \nu_* &= \left(c_1s_p\e^{\ii\frac{\Omega_1-\Omega_p}{2}} -s_1c_p\e^{-\ii\frac{\Omega_1-\Omega_p}{2}}\right)^2 . } \label{eq.munu} \end{equation}(D.1)According to Laskar & Boué (2010, Eq. (B.30)), the term in factor of αn in the perturbing function reads as n(0,0)=s=0nq=0n(Qs,q(n)˜(μ,ν)X0n,n2s(e1)X0(n+1),n2q(ep)×ei(n2s)ω1ei(n2q)ωp),\appendix \setcounter{section}{4} \begin{equation} \EQM{ {\cal F}_n^{(0,0)} = \sum_{s=0}^n\sum_{q=0}^n & \Big( \tilde Q_{s,q}^{(n)}(\mu_*,\nu_*) X_0^{n,n-2s}(e_1) X_0^{-(n+1),n-2q}(e_p) \crm & \times\, \e^{\ii(n-2s)\omega_1} \e^{\ii(n-2q)\omega_p} \Big) , } \label{eq.F00i} \end{equation}(D.2)Since we only consider the harmonics (ϖ1ϖp) and ((ϖ1ϖp)−(Ω1−Ωp)), we keep the terms such that q = ns, and drop the others. Consequently, one sum disappears from (D.2), and it remains n(0,0)=s=0n(Qs,ns(n)˜(μ,ν)X0n,n2s(e1)X0(n+1),n2s(ep)×ei(n2s)Δω).\appendix \setcounter{section}{4} \begin{equation} \EQM{ {\cal F}_n^{(0,0)} = \sum_{s=0}^n & \Big( \tilde Q_{s,n-s}^{(n)}(\mu_*,\nu_*) X_0^{n,n-2s}(e_1) X_0^{-(n+1),n-2s}(e_p) \crm & \times\, \e^{\ii(n-2s)\Delta\omega} \Big) . } \end{equation}(D.3)Then, we extract the linear terms in ep. Since X0(n+1),m(ep)e|m|p\hbox{$X_0^{-(n+1),m}(e_p)\propto e_p^{\abs{m}}$}, we consider only the values of s such that n−2s = 0, ± 1. Using the asymptotic expressions of the Hansen coefficients (A.5), we get 2n(0,0)=Qn,n(2n)˜(μ,ν)X02n,0(e1),\appendix \setcounter{section}{4} \begin{equation} {\cal F}_{2n}^{(0,0)} = \tilde Q_{n,n}^{(2n)}(\mu_*,\nu_*) X_0^{2n,0}(e_1) , \label{eq.F2n}\vspace*{-2mm} \end{equation}(D.4)and 2n+1(0,0)=Qn,n+1(2n+1)˜(μ,ν)X02n+1,1(e1)nepeω+cc,\appendix \setcounter{section}{4} \begin{equation} {\cal F}_{2n+1}^{(0,0)} = \tilde Q_{n,n+1}^{(2n+1)}(\mu_*,\nu_*) \,X_0^{2n+1,1}(e_1) \,n e_p \,\e^{\ii\Delta\omega} + cc , \label{eq.F2n1}\vspace*{-2mm} \end{equation}(D.5)where cc means complex conjugate. Then, we use the definition of the ˜Q\hbox{$\tilde Q$} functions (Laskar & Boué 2010, Eq. (B.31)). For s + q ≤ n, one has Qs,q(n)˜(μ,ν)=μqsνnqsAqs,nqs(n)(|ν|),\appendix \setcounter{section}{4} \begin{equation} \tilde Q_{s,q}^{(n)} (\mu_*,\nu_*) = \mu_*^{q-s} \nu_*^{n-q-s} A_{q-s,n-q-s}^{(n)}(\abs{\nu_*}) ,\vspace*{-2mm} \end{equation}(D.6)with Laskar & Boué (2010), Eq. (B.32) Aqs,nqs(n)(x)=122n(2s)!(2n2q)!s!(ns)!q!(nq)!×k=02s(1)k(2n2s+k)!(2sk)!(2n2q2s+k)!xkk!·\appendix \setcounter{section}{4} \begin{equation} \EQM{ A_{q-s,n-q-s}^{(n)} (x) = & \frac{1}{2^{2n}} \frac{(2s)!(2n-2q)!}{s!(n-s)!q!(n-q)!} \crm & \times \sum_{k=0}^{2s} (-1)^k \frac{(2n-2s+k)!}{(2s\! -\! k)!(2n\! -\! 2q\! -\! 2s\! +\! k)!} \frac{x^k}{k!} \cdot }\vspace*{-2mm} \end{equation}(D.7)We thus have Qn,n(2n)˜(μ,ν)=A0,02n(|ν|)=124n(2n)!(2n)!(n!)4k=0(2n)(1)k(2n+k)!(2nk)!(k!)2|ν|k,\appendix \setcounter{section}{4} \begin{equation} \EQM{ \tilde Q_{n,n}^{(2n)}(\mu_*,\nu_*) &= A_{0,0}^{2n}(\abs{\nu_*}) \crm &= \frac{1}{2^{4n}}\frac{(2n)!(2n)!}{(n!)^4} \sum_{k=0}^{(2n)} (-1)^k \frac{(2n+k)!}{(2n-k)!(k!)^2} \abs{\nu_*}^k , }\vspace*{-2mm} \end{equation}(D.8)or in a more condensed form, Qn,n(2n)˜(μ,ν)=f2n,nF(2n+1,2n,1;|ν|),\appendix \setcounter{section}{4} \begin{equation} \tilde Q_{n,n}^{(2n)}(\mu_*,\nu_*) = f_{2n,n} F(2n+1,-2n,1;\abs{\nu_*} ) , \label{eq.Qnn} \end{equation}(D.9)

where F is the hypergeometric function and fp,q is given in (A.3). In the same way, Qn,n+1(2n+1)˜=μf2n+1,nF(2n+3,2n,1;|ν|).\appendix \setcounter{section}{4} \begin{equation} \tilde Q_{n,n+1}^{(2n+1)} = \mu_* f_{2n+1,n} F(2n+3,-2n,1;\abs{\nu_*} ) . \label{eq.Qnn1} \end{equation}(D.10)Now, we substitute the expressions of μ ∗  and ν ∗  (D.1) into those of ˜Qn,n2n\hbox{$\tilde Q_{n,n}^{2n}$} (D.9) and of ˜Qn,n+12n+1\hbox{$\tilde Q_{n,n+1}^{2n+1}$} (D.10). Keeping only the linear terms in inclinations, μc12eiΔΩ+2c1s1sp,\appendix \setcounter{section}{4} \begin{equation} \mu_* \approx c_1^2 \e^{\ii\Delta\Omega} + 2 c_1 s_1 s_p , \label{eq.mu*} \end{equation}(D.11)and |ν|ks12k2kc1s12k1spcosΔΩ,\appendix \setcounter{section}{4} \begin{equation} \abs{\nu_*}^k \approx s_1^{2k} - 2 k c_1 s_1^{2k-1} s_p \cos\Delta\Omega , \label{eq.nu*} \end{equation}(D.12)with ΔΩ = Ω1−Ωp, and retaining only the terms involved in the secular resonances, we get 2n(0,0)=f2n,nF(2n+1,2n,1;s12)X02n,0(e1);2n+1(0,0)=2f2n+1,n(c12F(2n+3,2n,1;s12)cosΔϖ+c1s1spG2n+1,n(s12)cos(ΔϖΔΩ))X02n+1,1(e1)nep,\appendix \setcounter{section}{4} \begin{equation} \EQM{ {\cal F}_{2n}^{(0,0)} &=& f_{2n,n} F(2n+1,-2n,1;s_1^2) X_0^{2n,0}(e_1) ; \crm {\cal F}_{2n+1}^{(0,0)} &=& 2 f_{2n+1,n} \Big( c_1^2 F(2n+3,-2n,1;s_1^2) \cos\Delta\varpi \crm && \! \! \!\!\! + c_1 s_1 s_p G_{2n+1,n}(s_1^2) \cos(\Delta\varpi\! -\! \Delta\Omega) \Big)\, X_0^{2n+1,1}(e_1)\,ne_p , } \label{eq.Finc} \end{equation}(D.13)where G2n+1,n(x)=2F(2n+3,2n,1;x)(1x)F(2n+3,2n,1;x),\appendix \setcounter{section}{4} \begin{equation} \EQM{ G_{2n+1,n}(x) =& 2 F(2n+3,-2n,1;x) \crm & -(1-x)\,F'(2n+3,-2n,1;x) , } \end{equation}(D.14)and F(a,b,c;x)=abcF(a+1,b+1,c+1;x).\appendix \setcounter{section}{4} \begin{equation} F'(a,b,c;x) = \frac{ab}{c}F(a+1,b+1,c+1;x) . \end{equation}(D.15)The Eqs. (D.13) generalize the expression (4) obtained in the coplanar problem. The Hamiltonian of the inclined system is then inc=Gmpapn=2αnn(0,0)(e1,ep,ii,ip,ϖ1ϖp),\appendix \setcounter{section}{4} \begin{equation} \hat{H}_{\rm inc} = -\frac{G m_p}{a_p} \sum_{n=2}^\infty \alpha^n {\cal F}_n^{(0,0)} (e_1, e_p, i_i, i_p, \varpi_1-\varpi_p) , \label{eq.CHam} \end{equation}(D.16)with n(0,0)\hbox{${\cal F}_n^{(0,0)}$} given in (D.13).

All Tables

Table 1

Amplitude of the quadrupole εpquad\hbox{$\varepsilon_p^{\rm quad}$} (14) and octupole εpoctu\hbox{$\varepsilon_p^{\rm octu}$} (13) terms, due to each planet p, amplitude Ap, and phase ϕp of the eigenmode with frequency g5 = 4.2488163′′/yr in the quasiperiodic decomposition of zp = epexpiϖp (taken from Laskar 1990).

Table B.1

Coefficients of the polynomials P(x) =  ∑ px (B.4) and Q(x) =  ∑ qx (B.5) computed up to α50.

Table C.1

Frequencies at zero eccentricity and corrections for different orders of expansion of the perturbing function.

All Figures

thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 3

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

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.