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/00046361/201219991  
Published online  20 November 2012 
A simple model of the chaotic eccentricity of Mercury
^{1}
Centro de Astrofísica, Universidade do Porto,
Rua das Estrelas,
4150762
Porto,
Portugal
^{2}
Astronomie et Systèmes Dynamiques, IMCCECNRS UMR8028,
Observatoire de Paris, UPMC, 77 Av.
DenfertRochereau, 75014
Paris,
France
^{3}
Department of Astronomy and Astrophysics, University of
Chicago, 5640 South Ellis
Avenue, Chicago,
IL
60637,
USA
email: boue@oddjob.uchicago.edu
Received:
11
July
2012
Accepted:
17
October
2012
Mercury’s eccentricity is chaotic and can increase so much that collisions with Venus or the Sun become possible. This chaotic behavior results from an intricate network of secular resonances, but in this paper, we show that a simple integrable model with only one degree of freedom is actually able to reproduce the large variations in Mercury’s eccentricity, with the correct amplitude and timescale. We show that this behavior occurs in the vicinity of the separatrices of the resonance g_{1}−g_{5} between the precession frequencies of Mercury and Jupiter. However, the main contribution does not come from the direct interaction between these two planets. It is due to the excitation of Venus’ orbit at Jupiter’s precession frequency g_{5}. We use a multipolar model that is not expanded with respect to Mercury’s eccentricity, but because of the proximity of Mercury and Venus, the Hamiltonian is expanded up to order 20 and more in the ratio of semimajor axis. When the effects of Venus’ inclination are added, the system becomes nonintegrable and a chaotic zone appears in the vicinity of the separatrices. In that case, Mercury’s eccentricity can chaotically switch between two regimes characterized by either lowamplitude circulations or highamplitude librations.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: individual: Mercury / chaos / methods: analytical / methods: numerical / celestial mechanics
© 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 g_{1}−g_{5} and (g_{1}−g_{5})−(s_{1}−s_{2}) 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_{1}−g_{5} 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_{1}/a_{p})^{20} in the relativistic case and (a_{1}/a_{p})^{50} 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_{1}/a_{p}) in order to reach the asymptotic evolution. We also show that the maximum eccentricity attained with an initial eccentricity of e_{1} = 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_{0}, of Mercury u_{1}, of Venus u_{2}, 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^{1}(1)where is the conjugate momentum of r_{1}, G the gravitational constant, m_{0} 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_{0}/(2a_{1}), where a_{1} 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_{1}/r_{p})_{p = 2,8}(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_{2}. However, it is well known that the resulting doubleaveraged 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_{0}/r_{p} obtained for n = 0 do not depend on Mercury’s elements; the terms vanish when averaged over M_{1}; the term −Gm_{0}/(2a_{1}) becomes constant after averaging over M_{1}, since a_{1} becomes constant. As such, all these terms are dropped in the following expressions. Finally, the averaged expression of the perturbing functions a_{p}/r_{1}−r_{p} are given in Laskar & Boué (2010). The resulting Hamiltonian of the coplanar problem is then (3)where (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 are the Hansen coefficients defined by (5)where ϵ_{n} = 0 if n is odd, and ϵ_{n} = 1 if n is even, and (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_{1}. With this approximation, up to the octupole order, the averaged Hamiltonian (2) reads as (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) (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^{0} in P_{α}(α,x) is the Taylor expansion of −1, and the coefficient of x^{1} in P_{α}(α,x) is the Taylor series of , etc., where 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_{p}expiϖ_{p} is described well by a quasiperiodic expansion (Laskar 1990). In this study, we focus on the secular resonance g_{1}−g_{5}. In the quasiperiodic decomposition of the variables (z_{p})_{p = 2,8}, we thus keep only the terms associated to the eigenmode with frequency g_{5} ≈ 4.25′′/yr. The others average out and disappear from the Hamiltonian.
Amplitude of the quadrupole (14) and octupole (13) terms, due to each planet p, amplitude A_{p}, and phase ϕ_{p} of the eigenmode with frequency g_{5} = 4.2488163′′/yr in the quasiperiodic decomposition of z_{p} = e_{p}expiϖ_{p} (taken from Laskar 1990).
The amplitude and the phase of the mode 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, ϕ_{7} = 30.6 + 180deg, it is equivalent to take and . This is the convention that is followed hereafter, but the prime is omitted for clarity. We note , where . We also use the fact that is constant in the secular problem to rescale the Hamiltonian by this quantity as it simplifies the following expressions. We note Given that Gm_{p}/(Λ_{1}a_{p}) = n_{1}a_{1}m_{p}/(a_{p}m_{0}), the resonant Hamiltonian now reads as (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 (10)be the canonical conjugated variables of the initial Hamiltonian Ĥ_{N,plan} of (Eq. (8)). It can be easily shown that (11)are canonical conjugated variables of the Hamiltonian (Eq. (9)) without rescaling the time t. The equations of motion are thus (12)At this point, it is interesting to evaluate the contribution of each planet to the octupole interaction with Mercury. Although the resonance g_{1}−g_{5} involves only the precession frequencies associated to Mercury and Jupiter, the eigenmode 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 given by (13)The values taken by these parameters are gathered in Table 1. The maximal amplitude is due to Venus with followed by the EarthMoon barycenter and Jupiter with . Thus, the strongest perturbation on Mercury’s orbit comes from the precession of Venus excited by Jupiter. a_{1}/a_{2} ≈ 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 (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_{1} alone, given by (Eq. (9)) (15)Then, the resonant Hamiltonian (9) reads as (16)The numerical values of the coefficients of the polynomials P and Q are given in Appendix B.
The Hamiltonian (16) has one and a half degrees of freedom, but it can be reduced to a one degree of freedom Hamiltonian after some modifications.For that, we first make the Hamiltonian autonomous by a adding a momentum conjugated to time t. Then, we perform a canonical transformation where the old variables are (17)and the new ones are (I,Δϖ) and , defined by (18)For this transformation to be canonical, the momentums should verify (19)The new Hamiltonian expressed in the new variables does not depend on the cyclic coordinate . Thus, its conjugated momentum is an integral of the motion and can be dropped. Up to a constant, the new Hamiltonian is then(20)This one degree of freedom Hamiltonian is integrable. The orbits are given by . The temporal evolution of the eccentricity e_{1} and of the resonant angle Δϖ are deduced from the canonical equations of motion. This leads to (21)
2.2. General relativistic precession
The secular effect of relativity is described by (e.g., Touma et al. 2009), (22)where , and c is the speed of light. In the case of Mercury, we have g_{r} ≈ 0.41′′/yr. The total Hamiltonian , including the Newtonian interaction and general relativity, becomes (23)
Fig. 1 Level curves of the Hamiltonian in the plane (e_{1}cosΔϖ,e_{1}sinΔϖ). a) Quadrupole expansion with g_{1} = 5′′/yr. b) Octupole expansion with g_{1} = 5′′/yr. c), d), e), f) expansion up to the order n = 20 with g_{1} = 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 g_{1} = g_{5}. 
2.3. Additional control term
Using the Newtonian Hamiltonian (20), or the relativistic one (23), the system is far from the g_{1}−g_{5} 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_{1} = 5.54′′/yr in the Newtonian case and g_{1} = 5.96′′/yr in the relativistic case, whereas the resonance occurs in the vicinity of g_{1} = g_{5} = 4.25′′/yr. This result is slightly different, but representative of the present behavior of Mercury’s eccentricity. Indeed, the present value of g_{1} for the solar system with GR is g_{1} = 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_{1} by an increment δ_{g}. In addition, owing to the slow chaotic diffusion in the inner solar system, Mercury’s precession frequency g_{1} can change quasirandomly and come close to the value of g_{5}, 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_{1} around the resonant frequency g_{5}, the abovementioned correction is added to both and ; i.e., (24)and (25)In Eqs. (24) and (25), the factor δ_{g} controls the precession frequency g_{1}. For δ_{g} = 0, we recover the values obtained with and , respectively. Otherwise the value of g_{1} 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_{1}, but to show that e_{1} 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_{1} = g_{5}. 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_{1}/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., or . In a first step, we focus only on the relativistic Hamiltonian which is more complete.
Within the quadrupole approximation, the Hamiltonian reduces to a simple expression (26)where . This Hamiltonian is independent of Δϖ, thus e_{1} remains constant. All the trajectories in the plane (x = e_{1}cosΔϖ, y = e_{1}sinΔϖ) are circles centered on the origin (0,0) with an eventually infinite period when . Figure 1a illustrates this result for g_{1} = 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 and thus, no possible increase in the eccentricity.
The next level of approximation is the octupole. The corresponding Hamiltonian reads as (27)with . 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_{1} = 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_{1}. 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 g_{1}. 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 increases, the degrees of the polynomials P_{α}(α,e^{2}) and Q_{α}(α,e^{2}) increase both for α and e^{2}. Then, the topology of the phase space of the Hamiltonian 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_{1} ≲ 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_{1}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_{1} 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_{1} = 0.8. The separatrix disappears at g_{1} ≥ 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_{1} ≲ 4.5′′/yr. For example, with g_{1} = 4.0′′/yr (Fig. 1e), a trajectory starting within e_{1} ≤ 0.2 can reach e_{1} ≈ 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_{1} ≈ 0.2 up to e_{1} ≈ 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 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 (28)where x = e_{1}cosΔϖ and y = e_{1}sinΔϖ. They are equivalent to the equations of motion (21), but without singularities at the origin e_{1} = 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 (e_{1}cosΔϖ,e_{1}sinΔϖ) for two orders of expansion: n = 20 (left), n = 50 (right). In both cases, g_{1} = 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_{1} = 0.2, Δϖ = 110deg, and g_{1} = 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 (24). Figure 4 shows the level curves of and . The two figures are similar to the subpanel α^{20} of Fig. 1 within the region of low eccentricities (e_{1} ≲ 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_{1}, 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 . 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_{1}−g_{5} resonance.
To do so, we add inclination in our system, and focus on the term related to (g_{1}−g_{5})−(s_{1}−s_{2}). 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) (29)where i_{p} and Ω_{p} are the inclination and the longitude of the ascending node of the planet p, and I_{p} and are the amplitude and the argument of the term precessing at the frequency s_{2} in the quasiperiodic decomposition of sin(i_{p}/2)exp(iΩ_{p}), respectively. In (29), both and 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 conjugated to the time t. The old variables of the subsequent canonical transformation are (30)and the new variables are denoted (I,Δϖ), (J,ΔΩ), . We set (31)The associated conjugated momenta are given by (32)Since the new Hamiltonian is independent of t, its conjugate momentum is an integral of the motion. From the expressions of and (30), we thus get, up to a constant, (33)To obtain the final Hamiltonian H_{R,inc} in the spatial case, we add the control term and get (34)The Hamiltonian (34) now has two degrees of freedom. The variables are (e_{1},Δϖ) and (i_{1},ΔΩ). To perform numerical integrations, we use nonsingular rectangular coordinates (35)With , the equations of motion read as (e.g. Bretagnon 1974) (36)Since Venus is the only planet with a significative amplitude associated to the frequency s_{2}, 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 lowamplitude circulations and highamplitude 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_{1}−g_{5})−(s_{1}−s_{2})), in a state of slow chaos, with slow chaotic diffusion. This diffusion will quasirandomly change the value of g_{1}, which can then approach resonance with g_{5}. 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_{1}−g_{5} resonance.
Fig. 6 Top: trajectory of Mercury’s eccentricity in the plane (e_{1}cosΔϖ,e_{1}sinΔϖ) 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_{1} = g_{5}. 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 loweccentricity regime to a higheccentricity 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_{1}−g_{5})−(s_{1}−s_{2}) 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 loweccentricity to higheccentricity regimes occur quasirandomly, as observed in the full system, with maximal eccentricity close to 1 (Fig. 6).
For clarity, we drop the explicit time dependence in all equations after (1).
Acknowledgments
This work has been supported by PNPCNRS, by the CS of the Paris Observatory, by PICS05998 FrancePortugal 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/CTEAST/098528/2008, funded by the Fundação para a Ciência e a Tecnologia (FCT), Portugal.
References
 Batygin, K., & Laughlin, G. 2008, ApJ, 683, 1207 [NASA ADS] [CrossRef] [Google Scholar]
 Bretagnon, P. 1974, A&A, 30, 141 [NASA ADS] [Google Scholar]
 Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 579 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1984, Thèse, Observatoire de Paris http://tel.archivesouvertes.fr/tel00702723 [Google Scholar]
 Laskar, J. 1989, Nature, 338, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J. 1994, A&A, 287, L9 [NASA ADS] [Google Scholar]
 Laskar, J. 2008, Icarus, 196, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Boué, G. 2010, A&A, 522, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lidov, M. 1962, Planet. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L., & Ziglin, S. L. 1976, Celest. Mech., 13, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., & Wu, Y. 2011, ApJ, 739, 31 [NASA ADS] [CrossRef] [Google Scholar]
 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 twoplanet system where the massless body is on the inner orbit reads as (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_{1}/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 are given by Laskar & Boué (2010)(A.2)where ϵ_{n} = 1 if n is even and 0, otherwise (A.3)and 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 is em, 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), (A.4)for n ≥ 2, one gets (A.5)Then, the expression of simplifies as (A.6)Substituting this expression into (A.1), one obtains (A.7)where (A.8)and (A.9)To derive explicit expressions of these two quantities, we use the analytical formulae of the coefficients for n ≥ 2 (Laskar & Boué 2010), (A.10)Moreover, we invert the sums and separate the cases ℓ = 0 and ℓ ≥ 1 in P_{α}(α,x). This gives (A.11)and (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
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 (B.1)Then, using the quasiperiodic expansion of the variables (e_{p}e^{iϖp})_{p = 2,8} and keeping only the terms at the fundamental frequency g_{5}, one gets (B.2)We now define the polynomials P and Q such that (B.3)where . The new polynomials P and Q are derived from P_{α} and Q_{α} through (B.4)and (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 g_{1} of Mercury’s perihelia computed either from the Newtonian Hamiltonian (20), (C.1)or from the relativistic Hamiltonian (23) (C.2)is far enough from g_{5} ≈ 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: .
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 twoplanet 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 i_{1} and i_{p} the inclinations of the massless planet and of the perturber respectively, and Ω_{1} and Ω_{p} are their longitudes of ascending node. We also note c_{1} = cos(i_{1}/2), s_{1} = sin(i_{1}/2), c_{p} = cos(i_{p}/2), and s_{p} = sin(i_{p}/2), from which we define (D.1)According to Laskar & Boué (2010, Eq. (B.30)), the term in factor of α^{n} in the perturbing function reads as (D.2)Since we only consider the harmonics (ϖ_{1}−ϖ_{p}) and ((ϖ_{1}−ϖ_{p})−(Ω_{1}−Ω_{p})), we keep the terms such that q = n−s, and drop the others. Consequently, one sum disappears from (D.2), and it remains (D.3)Then, we extract the linear terms in e_{p}. Since , 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 (D.4)and (D.5)where cc means complex conjugate. Then, we use the definition of the functions (Laskar & Boué 2010, Eq. (B.31)). For s + q ≤ n, one has (D.6)with Laskar & Boué (2010), Eq. (B.32) (D.7)We thus have (D.8)or in a more condensed form, (D.9)
where F is the hypergeometric function and f_{p,q} is given in (A.3). In the same way, (D.10)Now, we substitute the expressions of μ_{ ∗ } and ν_{ ∗ } (D.1) into those of (D.9) and of (D.10). Keeping only the linear terms in inclinations, (D.11)and (D.12)with ΔΩ = Ω_{1}−Ω_{p}, and retaining only the terms involved in the secular resonances, we get (D.13)where (D.14)and (D.15)The Eqs. (D.13) generalize the expression (4) obtained in the coplanar problem. The Hamiltonian of the inclined system is then (D.16)with given in (D.13).
All Tables
Amplitude of the quadrupole (14) and octupole (13) terms, due to each planet p, amplitude A_{p}, and phase ϕ_{p} of the eigenmode with frequency g_{5} = 4.2488163′′/yr in the quasiperiodic decomposition of z_{p} = e_{p}expiϖ_{p} (taken from Laskar 1990).
Frequencies at zero eccentricity and corrections for different orders of expansion of the perturbing function.
All Figures
Fig. 1 Level curves of the Hamiltonian in the plane (e_{1}cosΔϖ,e_{1}sinΔϖ). a) Quadrupole expansion with g_{1} = 5′′/yr. b) Octupole expansion with g_{1} = 5′′/yr. c), d), e), f) expansion up to the order n = 20 with g_{1} = 2.5,3.6,4.0,5.0′′/yr, respectively. 

In the text 
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 g_{1} = g_{5}. 

In the text 
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 
Fig. 4 Level curves of the Hamiltonian without relativity in the plane (e_{1}cosΔϖ,e_{1}sinΔϖ) for two orders of expansion: n = 20 (left), n = 50 (right). In both cases, g_{1} = 3.6′′/yr. 

In the text 
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 
Fig. 6 Top: trajectory of Mercury’s eccentricity in the plane (e_{1}cosΔϖ,e_{1}sinΔϖ) 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 (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.