# Chaotic diffusion in the Solar System.

###### Abstract

The discovery of the chaotic behavior of the planetary orbits in the Solar System (Laskar, 1989, 1990) was obtained using numerical integration of averaged equations. In (Laskar, 1994), these same equations are integrated over several Gyr and show the evidence of very large possible increase of the eccentricity of Mercury through chaotic diffusion. On the other hand, in the direct numerical integration of (Ito and Tanikawa, 2002) performed without general relativity over Gyr, the eccentricity of Mercury presented some chaotic diffusion, but with a maximal excursion smaller than about 0.35. In the present work, a statistical analysis is performed over more than 1001 different integrations of the secular equations over 5 Gyr. This allows to obtain for each planet, the probability for the eccentricity to reach large values. In particular, we obtain that the probability of the eccentricity of Mercury to increase beyond in 5 Gyr is about 1 to 2 %, which is relatively large. In order to compare with (Ito and Tanikawa, 2002), we have performed the same analysis without general relativity, and obtained even more orbits of large eccentricity for Mercury. In order to clarify these differences, we have performed as well a direct integration of the planetary orbits, without averaging, for a dynamical model that do not include the Moon or general relativity with 10 very close initial conditions over 3 Gyr. The statistics obtained with this reduced set are comparable to the statistics of the secular equations, and in particular we obtain two trajectories for which the eccentricity of Mercury increases beyond in less than 1.3 Gyr and 2.8 Gyr respectively. These strong instabilities in the orbital motion of Mecury results from secular resonance beween the perihelion of Jupiter and Mercury that are facilitated by the absence of general relativity. The statistical analysis of the 1001 orbits of the secular equations also provides probability density functions (PDF) for the eccentricity and inclination of the terrestrial planets (Mercury, Venus, the Earth and Mars) that are very well approximated by Rice PDF. This provides a very simple representation of the planetary PDF over 5 Gyr. On this time-scale the evolution of the PDF of the terrestrial planets is found to be similar to the one of a diffusive process. As shown in (Laskar, 1994), the outer planets orbital elements do not present significant diffusion, and the PDFs of their eccentricities and inclinations are well represented by the PDF of quasiperiodic motions with a few periodic terms.

Keywords : Celestial Mechanics, Planetary Dynamics, Chaotic diffusion

^{1}

^{1}footnotetext: E-mail address: laskar@imcce.fr

## 1 Introduction

The discovery of the chaotic behavior of the planetary orbits in the Solar System (Laskar, 1989, 1990) was obtained using numerical integration of averaged equations. Since then, this chaotic behavior has been confirmed through direct numerical simulation, without averaging (Quinn et al., 1991, Sussman and Wisdom, 1992). More recently, integrations of more accurate planetary equations have been performed over 100 to 250 Myr (Varadi et al., 2003, Laskar et al., 2004a, b). Due to the chaotic evolution of the system, and to the uncertainty on the model and initial conditions, the interval of precise validity of these solutions is limited to about 40 Myr (Laskar et al., 2004a), even if some components of the solutions can be used over longer time for paleoclimate studies, as the 405 kyr oscillation of the Earth eccentricity (Laskar et al., 2004a).

Over periods of time longer than 100 Myr, it becomes hopeless to search for a precise solution for the orbital parameters of the inner planets (Mercury, Venus, Earth and Mars). On the other hand, it is important to understand the possible behavior of these solutions, and in particular of the possible variations of the action variables of the orbits (semi major axis, eccentricity and inclination).

A first study of the chaotic diffusion of the planetary orbits was made in (Laskar, 1994) is order to search for the maximum possible variations of the eccentricities and inclinations of the planets, using the secular equations.

In the present work, a statistical view of the chaotic evolution of the planetary orbits is described for all planets of the Solar System (Pluto is excluded). This study is based on the numerical integrations of 1001 different solutions of the averaged equations of the Solar System using very close initial conditions, compatible with our present knowledge. Some results concerning the orbital evolution of Mars were already presented in (Laskar et al., 2004b) and we will denote this paper PI in the following, since in many cases, we will refer to this previous paper to avoid duplication.

## 2 The secular equations

In order to investigate the diffusion of the orbits over 5 Gyr, we will use the secular equations of Laskar (1990), with some small modifications. The secular equations are obtained by averaging the equations of motion over the fast angles that are the mean longitudes of the planets. The averaging of the equation of motion is obtained by expanding the perturbations of the Keplerian orbits in Fourier series of the angles, where the coefficients themselves are expanded in series of the eccentricities and inclinations. This averaging process was conducted in a very extensive way, up to second order with respect to the masses, and through degree 5 in eccentricity and inclination, leading to truncated secular equations of the Solar System of the form

(1) |

where , with , ( is the longitude of the perihelion). The matrix is the linear Lagrange-Laplace system, while and gather the terms of degree 3 and 5.

The system of equations thus obtained contains some 150000 terms, but can be considered as a simplified system, as its main frequencies are now the precession frequencies of the orbits of the planets, and no longer comprise their orbital periods. The full system can thus be numerically integrated with a very large step-size of 200 to 500 years. Contributions due to the secular perturbation of the Moon and general relativity are also included (see Laskar 1990, 1996, and laskar et al., 2004b (PI) for more details and references).

This secular system is then simplified and reduced to about 50000 terms, after neglecting terms of very small value (Laskar 1994). Finally, a small correction of the terms of the matrix of (1), after diagonalization, is performed in order to adjust the linear frequencies, in a similar way as it was done in (Laskar 1990). Indeed, in the outer planetary system, terms of higher order are of some importance, but their main effect will be to slightly modify the values of the main frequencies of the system. The correction that is done here is a simple way to correct for this effect. With the present small adjustment, the secular solutions are very close to the direct numerical integration La2004 (Laskar et al., 2004a, b) over about 35 Myr (PI, Figs. 15 and 16). As noted in PI, this time is about the time over which the direct numerical solution itself is valid (Laskar et al., 2004a, Figs. 20, 21), because of the imperfections of the model. Moreover, as the step-size used in the secular equations is 200 years instead of days for the direct integration, over very long times the numerical noise will be smaller. It is thus legitimate to investigate the diffusion of the orbital motion over long times using the secular equations. The major advantage, besides reducing the roundoff errors, resides in the computation speed : the integration of the secular equations is 2000 times faster than the integration of the non-averaged equations, and we can compute a 5 Gyr solution for the Solar System in 12 hours on a Compaq alpha workstation (833 Mhz). We are thus able to make statistics over many solutions with close initial conditions. In these computations, our main limitation will be the huge amount of data generated by these numerical integrations.

variable | planet | offset |
---|---|---|

Mars | to | |

Earth | to | |

Venus | to | |

Jupiter | to | |

Earth | to |

500 | 1000 | 1500 | 2000 | 3000 | 4000 | 5000 | |
---|---|---|---|---|---|---|---|

0.35 | 14 | 44 | 86 | 127 | 228 | 328 | 426 |

0.40 | 2 | 8 | 17 | 37 | 81 | 153 | 219 |

0.50 | 0 | 0 | 0 | 1 | 8 | 28 | 48 |

0.60 | 0 | 0 | 0 | 0 | 2 | 10 | 21 |

0.70 | 0 | 0 | 0 | 0 | 1 | 8 | 14 |

0.80 | 0 | 0 | 0 | 0 | 1 | 8 | 12 |

0.90 | 0 | 0 | 0 | 0 | 0 | 6 | 9 |

500 | 1000 | 1500 | 2000 | 3000 | 4000 | 5000 | |
---|---|---|---|---|---|---|---|

0.35 | 25 | 75 | 128 | 165 | 280 | 366 | 427 |

0.40 | 4 | 21 | 38 | 52 | 113 | 180 | 243 |

0.50 | 0 | 0 | 0 | 0 | 6 | 19 | 33 |

0.60 | 0 | 0 | 0 | 0 | 0 | 6 | 10 |

0.70 | 0 | 0 | 0 | 0 | 0 | 6 | 10 |

0.80 | 0 | 0 | 0 | 0 | 0 | 2 | 8 |

0.90 | 0 | 0 | 0 | 0 | 0 | 0 | 2 |

## 3 Maximum excursion

For the present analysis, we have integrated 1001 orbital trajectories of the Solar System over 5 Gyr in negative time with very small variations of the initial conditions with respect to the nominal solution. The initial conditions of the secular system nominal solution are the same as Table 1 and 2 from (Laskar, 1986), and are derived from the initial conditions of the VSOP82 solution of (Bretagnon, 1982). The phase space of the secular system (1) is of real dimension 32. The various initial conditions for the 1001 cases are obtained with a small variation of the initial value of a single secular variable or , for a single planet, according to Table 1, leaving the other 31 initial variables equal to their nominal values.

As we have performed 1001 numerical simulations of the whole Solar System over 5 Gyr, it is impossible to display the detailed results of these integrations, and we have chosen here to describe the most significant features of these integrations. One important point, that was addressed in (Laskar, 1994) is the maximum value reached by the eccentricity or inclination, as a result of the chaotic diffusion of the trajectories. In (Laskar, 1994), only 5 solutions were followed, after some small change of initial conditions in every interval of 500 Myr. This was an economical way to obtain the maximal possible value for each parameter, but without any estimate of the probability to reach these values. In particular, I could demonstrate that it was possible for Mercury to reach very high values for its eccentricity, allowing eventually a close encounter with Venus, while similar crossings were not reached for the other planets.

Here we perform more conventional statistics, and all 1001 integrations have been followed over the whole 5 Gyr interval. As it was mentioned before, (Laskar, 1990, 1994), there is practically no diffusion for the outer planet system that behaves nearly as a quasiperiodic and regular system. On the opposite, there is a significant diffusion of the eccentricities and inclinations of the inner planets. The statistics on the maximum values reached by the eccentricity and inclination over different time intervals of 50, 100, 250, 500, 1000, 2000, 3000, 4000, and 5000 Myr are displayed in figure 1 for Mercury, Venus, The Earth and Mars.

Over 50 Myr, the effect of the diffusion is not yet noticeable, and all solutions reach the same maximum value. This is reflected by a vertical curve in figure 1. But beyond this time interval, significant differences appear. In order to make these statistics more readable, the lower part of these maximum graphs, corresponding to the 5 percents of the solutions with largest variations of the eccentricities and inclination has been enlarged

### 3.1 Discussion

An important aspect of this study is the estimate of the probability for the eccentricity of Mercury to reach high values, allowing a possible close encounter with Venus. Over 1 Gyr, all solutions remained with , and only 0.1% of the solutions went beyond over 2 Gyr (Table 2). On the other hand, over longer time interval, the chaotic diffusion allowed Mercury’s eccentricity to reach very high values. In our simulations, 0.9% of the solutions reached an eccentricity of 0.9 within 5 Gyr, and 0.6% within 4 Gyr. At this point, one should remind that a close encounter of Mercury and Venus is only possible if the eccentricity of Mercury reaches values of about 0.75 (assuming an eccentricity of Venus of 0.06). We still have in our simulations about 1 % of the solutions that would allow for a close encounter with Venus.

In order to test the stability of the results displayed in Table 2,
I have also plotted in Table 3 the results obtained with a similar
experiment, but performed in positive time.
The equations are the same, but the integration step is now positive.
The simulation is made over 478 different orbits^{1}^{1}1This odd number of cases is purely accidental.
Our dedicated parallel computer broke down, and as this positive time integrations are used only as a check,
we considered that 478 cases were sufficient for this purpose., but the results are scaled
to 1000 for easier comparison with Table 2.
In a same way as for a diffusive process, the results in
positive time are very similar to the results obtained in negative time,
which is what was expected.

In doing these estimates, we need to keep in mind that the equations that are integrated here are the averaged equations of motions, where the disturbing function of the mutual interactions of the planets is expanded in series of eccentricities and inclinations. This expansion is divergent for high eccentricities. Indeed, already in the two-body problem, the expansion of the eccentric anomaly in powers of excentricity has a radius of convergence (see Wintner, 1947). Additionnally, the inverse of the distance becomes singular at collision, that is for an eccentricity of Mercury of about , depending on the eccentricity of Venus. Moreover, in the vicinity of collisions, numerous mean motion resonances overlap, giving rise to mean motion chaotic behavior, with possible changes in the semi-major axis of the planets, while these are constant in the secular system (see Fig. 8).

Nevertheless, although the truncations of series expansions involved in the construction of the secular system are made without estimates of the remainders, as it is usually the case in astronomy, I conjecture here that up to an eccentricity of Mercury of about , the dynamics of the full system is well represented by the dynamics of the secular system. Indeed, it has been observed that in many cases, the range of validity of secular systems extends much beyond theoretical estimates (see for example Libert and Henrard, 2006).

Moreover, my assumption is that in general, the full non averaged system is less regular than the secular system. I can thus assume that for eccentricities of Mercury below the data displayed in Table 2 provide good estimates of the actual probabilities to reach these eccentricity values, while for eccentricities above , the data given in Table 2 are only lower bounds of these probabilities.

It is certainly desirable that the same statistical studies should be conducted with the full equations of motions, without averaging, although this will require a considerable amount of CPU time that is still difficult to obtain.

On the other hand, the data displayed in Table 2 provide substantially more information than in (Laskar, 1994) where only the possibility of reaching high values of Mercury’s eccentricity allowing for a collision with Venus was demonstrated. Here we show that the probability to reach these high values is relatively large (about 1 to 2%). It thus becomes conceivable to prepare in the near future a full scale numerical simulation of the same problem with the non averaged equations.

### 3.2 Comparison with direct integration

It would be interesting to compare with an integration of the full equations of motion over Gyr time scale, but the author is not aware of a single numerical integration of the kind with a comparable dynamical model, including the contribution of the Moon and of general relativity. The long term integrations of (Varadi et al.2003, Laskar et al., 2004a,b) use precise models but span only a few 100 Myr. On the other hand, the long term integration of Ito and Tanikawa (2002) is performed over a few Gyr but do not include the relativistic contribution.

#### 3.2.1 The integration of Ito and Tanikawa

The numerical integration of Ito and Tanikawa (2002) is an integration of the Newtonian planetary equations, without relativistic contributions. It does not comprise either the effect of the Moon as a separate body. In fact, the mass of the Moon has been added to the mass of the Earth for five integrations () that span from 3.9 to 5 Gyr in the past () or in the future (). The initial conditions are taken from the JPL ephemeris DE245 (Standish, 1990). The integrator is a second order symplectic integrator (Wisdom and Holman 1991) with a step size of 8 days. From the figure 1 of (Ito and Tanikawa 2002) we can estimate that the maximum relative error in angular momentum is about for all solutions except for one for which it is about , the difference being probably due to different hardware. The maximum relative error in energy is about .

The authors do not provide very precise details on the maximum values reached by the eccentricities, but mention that the general behavior of the eccentricity of Mercury is similar with the results of (Laskar, 1994, 1996), although they obtained a maximum eccentricity for Mercury of about 0.35 over Gyr.

If we consider the numerical experiment that we conducted here with the secular system in negative time (Table 2) over 1001 solutions, we have only 32.8 % of the solutions that lead to an eccentricity of Mercury larger than 0.35, while 71.4 % had a maximal eccentricity larger than 0.32, and only 26.4 % an eccentricity larger than 0.36. There is indeed a very rapid decrease of the probability to reach a given eccentricity in the vicinity of 0.35.

Due to the lack of precision on the maximal value reached by the eccentricity of Mercury in (Ito and Tanikawa 2002), we can consider that the event () reached in their simulations has, according to our experiment on the secular system, a probability of about 75% to occur. As they made 5 simulations, the resulting probability would be , that is about 24%. This is not very large, but not unrealistic.

Nevertheless, there is a difference in the two experiments, as our secular equations comprise the relativistic contributions. One should thus also test the behavior of the secular system in absence of the relativistic contribution.

#### 3.2.2 The secular system without relativity

I have thus repeated the previous simulations in absence of relativity, expecting to find a more stable system. But the result was the opposite, and most of the solutions lead to high values of the eccentricity in 5 Gyr (Table 4). One can see that in this case, the probability of the event () becomes less than 25 %, and the probability of having 5 solutions with this behavior (if they are not related) becomes totally unrealistic with a value smaller than , that is about .

500 | 1000 | 1500 | 2000 | 3000 | 4000 | 5000 | |
---|---|---|---|---|---|---|---|

0.35 | 130 | 341 | 478 | 558 | 692 | 763 | 812 |

0.40 | 75 | 249 | 373 | 449 | 589 | 684 | 747 |

0.50 | 24 | 118 | 226 | 306 | 442 | 552 | 640 |

0.60 | 16 | 76 | 169 | 238 | 364 | 476 | 564 |

0.70 | 14 | 67 | 150 | 218 | 343 | 454 | 541 |

0.80 | 12 | 63 | 141 | 209 | 331 | 442 | 531 |

0.90 | 12 | 61 | 138 | 202 | 325 | 441 | 530 |

### 3.3 A new direct integration without relativity

In order to clarify the situation, it was thus necessary to make a new direct numerical integration. Indeed, the probability of reaching high values of the eccentricity of Table 4 are so high that one should be able to obtain solutions with high eccentricity with a moderate number or trials. I thus decided to integrate 10 orbits, , with very close initial conditions over 3 Gyr.

500 | 1000 | 1500 | 2000 | 3000 | |
---|---|---|---|---|---|

0.35 | 1 | 2 | 4 | 4 | 6 |

0.40 | 0 | 1 | 3 | 3 | 5 |

0.50 | 0 | 0 | 1 | 1 | 2 |

0.60 | 0 | 0 | 1 | 1 | 2 |

0.70 | 0 | 0 | 1 | 1 | 2 |

0.80 | 0 | 0 | 1 | 1 | 2 |

The model comprises the 8 major planets and the dwarf planet Pluto with Newtonian interactions. The Earth-Moon barycenter with the sum of the masses of the Earth and Moon is used instead of the Earth. The integrator is the SABA4 symplectic integrator of (Laskar and Robutel, 2001) that was already used in (Laskar et al., 2004a,b) with a step size of years. For a perturbed Hamiltonian of the form , using this integrator with a step is equivalent to integrate exactly a close by Hamiltonian where the error of method is of the order . The integration is conducted in extended precision on Itanium II processors with 80 bit arithmetics.

All parameters and initial conditions of the nominal solution are the same as the ones used in the new high precision planetary ephemeris INPOP06 that has been developped in our group. The reader should refer to the associated publication (Fienga et al., 2008) and web site (www.imcce.fr/inpop) for more details. The 10 different solutions () have the same initial conditions as the nominal solution, except for a shift of cm in the initial position of the Earth.

The maximum values reached by the eccentricity of Mercury in these 10 solutions have been plotted in figure 3. One can see that the present results differ substantially from the results of Ito and Tanikawa (2002), as the eccentricity can reach much higher values than in their simulations. In particular, for 5 solutions, the eccentricity went beyond , and for two of them, and , the eccentricity increased to very high values, beyond , in respectively 2776 and 1286.8 Myr (Fig. 5).

In the symplectic integration, the total angular momentum is conserved. The relative variation of the total angular momentum is thus an estimate of the roundoff error in the integration. For all solutions, this error remains below over the total length of the integration (Fig.2). The relative error in energy is in general below , except when the eccentricity of Mercury reaches high values (Fig.4). For the solution that reaches very high values of the eccentricity, the step size was changed to yr on the interval Gyr, and then to yr over the interval Gyr in order to achieve a good conservation of the total energy. With these settings, even for an eccentricity of 0.8, the relative error of the energy remains below over the whole interval of integration (Fig. 5).

When we gather the results of these 10 numerical simulations in a table (Tab. 5) similar to the tables 2, 3, 4, one can see that in our numerical results, the number of high eccentricities are somewhat lower, but quite comparable to the results of our table 4, especially if we consider the small number of events (10) on which the present statistics are made.

The present results, that are conducted with slightly more accurate integrations than the ones of Ito and Tanikawa (2002), demonstrate clearly that the large excursions of the eccentricity of Mercury described in (Laskar, 1994) are actually present in the full integration of the Newtonian equations. Moreover, the statistics obtained with the secular equations in absence of general relativity (Tab. 4) are very close to the same statistics obtained with the direct integration (Tab. 5).

We can thus assume that the statistics obtained over 1001 integrations of the secular equations (Tab. 2) are representative of the full system. Actually, as was already stated in (Laskar, 1994), we expect even that the full equations of motion will be slightly more unstable than the secular equations, especially for high eccentricities, when the overlap of mean motion resonances will increase the chaoticity of the orbits.

### 3.4 The effect of relativity

The difference of behavior of the secular system with (Table 2) and without general relativity (GR) (Table 4) is impressive. The results of the direct integrations performed without GR are also in good agreement with the results of the secular equation in absence of GR. The contribution of GR is thus essential in order to ensure the relative stability of Mercury. It is therefore important to understand the contribution of GR in this problem.

The most obvious effect of GR is to increase the perihelion frequencies of the planets and especially of Mercury. Indeed, the value of the secular frequency related to Mercury is in the vicinity of the origin with a contribution of from general relativity. The contribution of GR to the perihelion velocity decreases to /yr for Venus, for the Earth, for Mars, and only for Jupiter (e. g. Laskar, 1999, Table 4). The main frequency of the longitude of perihelion of Jupiter, (Laskar et al., 2004a, Table 3) is thus not changed much by GR and the main effect of GR is to increase the difference from in absence of GR to in presence of GR.

In order to analyse the effect of this change in the dynamics of Mercury, I have performed several integrations of the secular system, with and without GR, for various values of the initial eccentricity of Mercury. The initial eccentricity of Mercury in the nominal solution is , while in the present experiment, this value is varied with a step of from 0 to about until the integration crashes rapidly. The integrations are computed over 40 Myr. A frequency analysis is performed in order to compute the secular frequencies of the system. For each orbit, This frequency analysis (Laskar, 1990, 1993) is made over slidings 20 Myr intervals with a 1 Myr offset. The values obtained for are reported in figure 6 with (a1) and without (a2) GR.

The differences between the two plots are striking. With GR, when the eccentricity increases, is modified, but mostly increases and although presents large variations resulting from the chaotic behavior of the system, the system never gets trapped into the resonance. On the opposite, in absence of GR, the values are smaller, and as the eccentricity increases beyond , there exists a large chaotic zone related to the secular resonance , the secular frequency being plotted as an horizontal line in Figs 6a1,a2. When the system is in this resonance, the eccentricity of Mercury is driven to very large excentricities, beyond 1 (the eccentricity is not bounded by 1 in the secular system). In figure 6b2, one trajectory starting at about also increases to high values. Indeed, most probably, the chaotic region extends in the region of eccenricity, but it will take some additional time to reach the region of very strong chaos beyond .

This study of the secular system is made here to explore rapidly the phase space of the system, and to guide our intuition. Once we see the importance of the resonance in the secular system, we can test its contribution on the full equations with minimal computations.

We can now verify that the large increase of eccentricity of the solution plotted in figure 5 is actually triggered by the secular resonance . In figure 7, is plotted the difference of the longitudes of perihelion of Mercury () and Jupiter () from Myr to Myr. It appears that this angle is circulating until the two perihelions get locked from about Myr to Myr, that is over the time interval that corresponds to a steady increase of Mercury’s eccentricity. In a similar way, the solution presents as well a large increase of Mercury’s eccentricity from to about in only 2.5 Myr that corresponds to a locking of the perihelions of Mercury and Jupiter from to Myr.

It should be noted that for Mercury and Jupiter, and contrarily to the Earth, the relation with the longitudes of perihelion and secular frequencies , respectively, is straightforward. Indeed, in both cases, the secular frequency is clearly the leading periodic term of the quasiperiodic expansion of (see Table II from Laskar, 1990).

In order to see the real effectiveness of this resonance, I have performed an additional single numerical experiment, with the full equations of motion, including GR, and starting with the nominal intitial conditions, that is the initial conditions of the planetary ephemeris INPOP06 (Fienga et al., 2008), without any change of the initial eccentricity of Mercury. But in order to set the system into the resonance from the starting time of integration, I have changed the value of the post Newtonian parameter to instead of in standard GR, while keeping (Will, 2006). The contribution factor of perihelion shift is then instead of . The effect of this modified relativity is thus to artificially decrease the secular frequency by , instead of adding as in standard GR. The system is thus from the begining in the resonance, and the effect is immediate as Mercury’s eccentricity increases beyond 0.7 in less than 4 Myr (Fig. 8). Once in this region of high eccentricity, strong chaotic behavior due to short period resonances induces significant changes in the semi-major axis of Mercury (Fig. 8).

## 4 Density functions

The maximal possible value of the planet eccentricity is important for the analysis of the system stability, but it concerns the most exceptional orbits of the system, and not necessarily the most probable behavior of the Solar System over its age. For the understanding of the general behavior of the Solar System, the density function for the eccentricity and inclination of the planets is a complementary information that can be very valuable for the analysis of many physical parameters during the evolution of the Solar System. It was for example useful for the understanding of the capture of Mercury into the 3/2 spin orbit resonance (Correia and Laskar, 2004), or for the analysis of the past climate evolution of Mars (Laskar et al., 2004b). Here we have systemized the approach elaborated in PI for all the planets of the Solar System, with statistics over 1001 different orbits in order to obtain a complete statistical view of the variations of the orbital parameters of the Solar System.

Because of the chaotic behavior of the system, we know that it will never be possible to retrieve precisely the past (or future) orbital evolution of the Solar System over more than a few tens of millions of years (Laskar, 1989, 1990). On the other hand, beyond about 250 Myr one can obtain very smooth density functions for the eccentricity of the planets that will tell us the main general behavior of the orbits. Moreover, these density functions are some of the only accurate informations one can obtain beyond a few 100 Myr.

As in PI, we have divided here the time interval in 250 Myr intervals, and statistics are done over each 250 Myr interval, using the set of 1001 orbital solutions in negative time with the initial conditions of Table 1, for which the output has been recorded with a 1000 yr step size. The first 250 Myr interval is discarded, as the randomization due to the chaotic evolution of the system has not yet taken place (see PI for a more complete analysis of Mars eccentricity solution), and the normalized probability distributions functions (PDF), for the eccentricities and inclinations of all the planets are displayed in figures 9 and 10.

### 4.1 Discussion

The difference in the density functions for the eccentricity (or inclination ) of the inner (Mercury, Venus, Earth, Mars) or outer planets (Jupiter, Saturn, Uranus, Neptune) is striking on these plots. First of all, the shape of the PDF is different. For the inner planets, the PDF is similar to a Gaussian curve, although slightly different. In particular, all the curves have zero values for , with a positive, linear slope, while they have some long tail for large values of the eccentricity.

On the opposite, for the outer planets the density curves are strictly confined between two non zero values, and the curves have two peaks close two the minimum and maximum value of the eccentricity (resp. inclination).

Another striking feature is the fact that for the inner planets, we see clearly that numerous curves are displayed in each plot, as the density curve is different for each time interval of 250 Myr. This is the result of a significant diffusion that occurs in the eccentricity and inclination (Laskar, 1994), while for the outer, as the diffusion is practically non existent, all the 19 density curves that are actually plotted in figures 9 and 10 are virtually identical and superpose nearly exactly.

The difference between the two kinds of PDF will even be more striking in the next section, as we will attempt to associate to each PDF the density function of a simple dynamical system. For the inner planets, we found that the densities were very well modelized by a Rice distribution, (Rice, 1945) that is characteristic of either the modulus of a random walk in two dimensions, or alternatively to a quasiperiodic signal with noise (Rice 1945). We will denotes these density functions Rice densities.

On the opposite, the outer planets density function is in fact representative of a density function of a quasiperiodic signal with a moderate number of periodic terms. We will call such a density function a quasiperiodic density.

In the next section, we will examine more closely these two different densities.

## 5 Probability density functions

### 5.1 Rice densities

The Rice distribution is a continuous probability distribution with density function

(2) |

where is the modified Bessel function of the first kind obtained with its series expansion

(3) |

This distribution is obtained for example for the modulus of where are two Gaussian independent variables with variance and mean , with .

These PDF could thus be well suited for the eccentricity, if the rectangular variables and become practically Gaussian random variables because of the chaotic diffusion.

A Rice PDF has been successfully fitted to the eccentricity and inclination of the inner planets for each of the 250 Myr time intervals. It is impossible to display here graphics demonstrating the quality of the adjustment for all of these intervals, so we will reproduce here only the most representative cases, for the eccentricity of Mercury (Fig.11), Venus (Fig.12), Earth (Fig.13), and Mars (Fig.14), at 3 different epochs, while all cases can be summarized in Table 6.

It is indeed remarkable that the PDF of the eccentricity and inclination of the planets are such smooth functions. It is as well remarkable that they are so well approximated by simple 2 parameters PDF () of a single curve family with relatively simple expression (2).

### 5.2 Diffusion over 5 Gyr

1 | 0.1875 | 2.070e-03 | 1.043e-03 |
---|---|---|---|

2 | 4.197e-04 | 5.110e-06 | |

3 | 3.181e-04 | 3.323e-06 | |

4 | 1.002e-03 | 9.127e-05 | |

1 | 4.9896 | 9.040e+00 | 1.272e+00 |

2 | 1.5864 | 1.492e+00 | 2.377e-02 |

3 | 1.5803 | 1.063e+00 | 2.308e-02 |

4 | 3.623e+00 | 2.398e-01 |

As the PDF of eccentricity and inclination are well approximated for various epochs by Rice PDF, we have derived by least square fit, the values of the parameters , of these PDF for all variables and all time intervals of 250 Myr. Except for Mars, the value (related to the mean) do not present large variations over time (Figs. 15, 16, 17, 18). On the opposite, the standard deviation parameter increases with time, and we have a nearly linear relation

(4) |

similar to a diffusion process. In Table 6 are gathered the fitted values of over 5 Gyr.

We have thus in equation 2 and Table 6 some very simple formulas that will allow one to represent the PDF of the eccentricity and inclination of the inner planets over very long times, of several Gyr. Quite remarkably, although it is impossible to predict the precise evolution of the individual trajectories, we are able to give very simple expressions that fit well with the observed eccentricity and inclination PDF of the inner planets.

These formulas can then be used for the analysis of the past or future behavior of the Solar System. They could be used for example to compute an estimate of the capture probability of Mercury in the 3/2 resonance without requiring heavy numerical computations (Correia and Laskar, 2004).

We have used here a different PDF than in (Laskar et al., 2004b), where a model of a random walk with an absorbing edge at zero was used for Mars eccentricity. In fact, the results obtained with the two PDF are very similar in this case, but we prefer the Rice formulation that can be more easily interpreted. Moreover, we see here that the evolution of the system is well describe by a diffusive process for the rectangular coordinates that behave like random Gaussian variables on long time scale. Actually, the variance of the eccentricity evolves linearly with time. This is somewhat different from the observation made in (Laskar, 2004) for Mars eccentricity, but this is because here we are searching for models that fit more generally with the behavior of all the inner planets, for both eccentricity and inclination.

## 6 Outer planets

(”/yr) | (deg) | |||
---|---|---|---|---|

0 | 0.045570 | |||

1 | 23.987561 | 0.015372 | -82.213 | |

2 | 47.975121 | 0.001849 | 15.573 | |

3 | 1.169512 | 0.001680 | -90.509 | |

4 | 25.157074 | 0.000503 | -172.658 | |

5 | 22.818063 | 0.000491 | -171.688 |

The PDF of the outer planets (Figs. 9, 10) are very different from the PDF of the inner planets. Moreover, as it was already stated, the diffusion is practically inexistent, and all PDF curves over the different time intervals are virtually identical. The shape of these outer planets PDF can be understood as the PDF of some quasiperiodic functions of time with a small number of harmonics. Indeed, for a simple periodic function, , the PDF is

(5) |

where is the characteristic function^{2}^{2}2
The characteristic function of the interval is defined
as for , and if .
of the interval (Fig.19).

One can now understand that when for , instead of a single sine term, we have several periodic terms, the PDF of will become slightly distorted with respect to a pure sine function, and we will recover the specific form of the PDF of the outer planets eccentricities and inclinations (Figs. 9, 10). In order to illustrate this, we will approximate the eccentricity of Jupiter with a few periodic terms, obtained through frequency map analysis (Laskar, 1990, 2005).

With a frequency analysis of the eccentricity of Jupiter over 50 Myr, we obtain a quasiperiodic approximation of the eccentricity with 5 periodic terms on the form

(6) |

where the values of are given in table 7. The eccentricity of Jupiter over 50 Myr is plotted in figure 20, as well as its difference with this quasiperiodic approximation. The density function obtained with this quasiperiodic approximation is plotted in figure 21 in dotted line together with the density function of the full numerical integration (in full line). The two PDF are very close. It will be the same for the other outer planets. For these planets with a motion that is very close to quasiperiodic, the quasiperiodic decomposition obtained over 50 Myr as in (Laskar, 1990) provides in fact a very good representation of the orbit over several Gyr. It is thus not necessary to reproduce these expressions here, and we would rather provide in a forthcoming publication some full and accurate quasiperiodic representations for the motion of the outer planets, in agreement with our latest adjustment of the planetary ephemeris to observations (Fienga et al., 2007).

## 7 Conclusions

In (Laskar 1994), I showed the possibility of a very large increase of the eccentricity of Mercury, allowing for a close encounter or a collision with Venus. But in this work, there was no estimate of the probability for such an event to take place within a few Gyr. Here, such an estimate is given, by the extensive study of more than 1001 orbits. Indeed, it is found that the probability for Mercury’s eccentricity to exceed 0.6 within 5 Gyr is about 1 to 2 % which can be considered as a large value for such an important event.

When the contributions of general relativity and the Moon are not taken into account, this probability increases in a large amount and in the present numerical simulations, nearly half of the orbits went beyond 0.6 in less than 4 Gyr. As this appear to be in contradiction with the numerical results of (Ito and Tanikawa 2002), I have also performed a direct numerical simulation of the Newtonian equations, with 10 nearby initial conditions, without the Moon or general relativity over 3 Gyr, and found results that are in good agreement with the results of the integration of the secular system. In particular, two orbit were found for which the eccentricity of Mercury rises to very large values, beyond 0.8, thus allowing for a close encounter or a collision with Venus.

This direct numerical simulation is performed with slightly better accuracy than the numerical integration of (Ito and Tanikawa 2002). One should wonder about the reason of such a different behavior in the present computations and in the simulation of (Ito and Tanikawa 2002). The most probable reason for this is that the two integrations do not have the same initial conditions or model. The solutions of (Ito and Tanikawa 2002) may evolve in a slightly more regular region of the phase space. On the other hand, the secular equations are in very good agreement with the present direct integrations. This comfort their reliability and usefulness as for the secular equations, it was possible to perform an extended statistical study over 1001 solutions. One could perform additional direct numerical integrations in order to verify the probability law for the excursion of Mercury’s eccentricity, but this is of limited interest if it concerns the model of pure Newtonian equations, and not the full model including general relativity.

Indeed, as we have demonstrated here in section 3.4, the contribution of general relativity changes in a considerable manner the behavior of the Solar System dynamics. Indeed, in absence of GR, the secular frequency of the perihelion of Mercury decrease by , and becomes closer to the secular frequency of the perihelion of Jupiter . As the eccentricity varies under the secular planetary perturbations, the system can enter into the secular resonance that can drive Mercury’s eccentricity to very high values, beyond , where additional short period chaotic behavior occurs, inducing changes of semi-major axis of the planet (Fig. 8). This is indeed the mechanism that is present in our two solutions of the Newtonian equations for which Mercury’s eccentricity increased beyond (Fig. 7).

An important result of the present paper is to show this possibility of very large increase of Mercury’s eccentricity, beyond 0.8, allowing for a possible collision with Venus. It was actually possible to see that with such a large value of the eccentricity, the planet entered a region of mean motion resonances and chaos inducing changes in its semi-major axis.

It is clearly desirable to conduct now a full scale numerical experiment with the complete equations including general relativity and the contribution of the Moon, in order to provide some precise results on the chaotic behavior and possible evolution of our Solar System. We are indeed planning such a numerical study for the near future.

Apart from demonstrating the possibility of very large values for the eccentricity of Mercury, as was forecasted in (Laskar, 1994), I provide here probability density functions (PDF) for the eccentricities and inclinations of the terrestrial planets. It is in fact quite remarkable that the PDF of all terrestrial planets, when computed over long time, beyond about 500 Myr, become very smooth functions that are well approximated by a very simple two parameters function, namely the Rice distribution. This provides very compact formulas for the evolution of these PDF over time for all the terrestrial planets. The evolution of the PDF with time is similar to the one of a diffusion process with a linear increase of the variance variable with time.

For the outer planets, we have seen here that the chaotic diffusion is extremely small, without practical change of the PDF over time. The PDF being well approximated by the PDF of a quasiperiodic approximation of the solution. In the full equations of motion of the outer Solar System, there is some chaotic behavior that has been described (Sussman and Wisdom, 1992, Murray and Holman, 1999, Guzzo, 2005), but although a more detailed analysis should be made, It seems that the diffusion induced by this intrinsic chaotic behavior resulting from mean motion interactions should be smaller than the diffusion induced by the chaotic behavior of the inner Solar System that is already included in the present work. It is thus doubtful that the PDF for the outer planets that have been computed here will be much changed by the consideration of the full dynamics of the outer Solar System.

More generally, the PDF that are given here are the PDF obtained with the secular equations. Nevertheless, they should be very stable with respect to small changes in the model or in the involved parameters, and we can conjecture as well that they are very close to the PDF of the full system.

### Acknowledgments

The assistance of Mickael Gastineau is greatly acknowledged. A large part of the computations were made at IDRIS-CNRS. This work also benefited from support from PNP-CNRS, and CS from Paris Observatory.

.

## References

- 1 Bretagnon, P.: 1982, Theory for the motion of all the planets - The VSOP82 solution A&A, 114, 278–288
- 2 Correia, A., Laskar, J.: 2004, Mercury’s capture into the 3/2 spin-orbit resonance as a result of its chaotic dynamics Nature, 429, 848-850
- 3 Fienga, A., Manche, H., Laskar, J., Gastineau, M.: 2008, INPOP06: a new numerical planetary ephemeris, A&A, 477, 315–327
- 4 Guzzo, M. 2005. The web of three-planet resonances in the outer Solar System. Icarus 174, 273–284.
- 5 Ito, T., Tanikawa, K. 2002. Long-term integrations and stability of planetary orbits in our Solar System. MNRAS 336, 483–500.
- 6 Laskar, J.: 1986, Secular terms of classical planetary theories using the results of general theory A&A, 157, 59–70
- 7 Laskar, J.: 1989, A numerical experiment on the chaotic behaviour of the Solar System, Nature, 338, 237–238
- 8 Laskar, J.: 1990, The chaotic motion of the Solar System. A numerical estimate of the size of the chaotic zones, Icarus, 88, 266–291
- 9 Laskar, J.: 1993, Frequency analysis for multi-dimensional systems. Global dynamics and diffusion, Physica D, 67, 257–281
- 10 Laskar, J.: 1994, ’Large scale chaos in the Solar System’, Astron. Astrophys., 287, L9–L12
- 11 Laskar, J.: 1999, The limits of Earth orbital calculations for geological time scale use, Phil. Trans. R. Soc. Lond. A., 357, 1735–1759
- 12 Laskar, J., Robutel, P.: 2001, High order symplectic integrators for perturbed Hamiltonian systems, Cel. Mech., 80, 39–62
- 13 Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., Levrard, B.: 2004a, A long term numerical solution for the insolation quantities of the Earth, A&A, 428, 261–285
- 14 Laskar, J., Correia, A., Gastineau, M., Joutel, F., Levrard, B., Robutel, P.: 2004b, Long term evolution and chaotic diffusion of the insolation quantities of Mars, Icarus, 170, 343–364
- 15 Laskar, J: 2005, Frequency Map analysis and quasi periodic decompositions, in ”Hamiltonian systems and Fourier analysis”, Benest et al., eds, Taylor and Francis
- 16 Libert, A-S., Henrard, J. 2006. Secular apsidal configuration of non-resonant exoplanetary systems. Icarus 183, 186–192.
- 17 Murray, N., Holman, M. 1999. The Origin of Chaos in the Outer Solar System. Science 283, 1877–1881.
- 18 Quinn, T. R., Tremaine, S., and Duncan, M. 1991. A three million year integration of the earth’s orbit. Astron. J. 101, 2287–2305.
- 19 Rice, S.O. 1945. Mathematical Analysis of Random Noise. Bell System Technical Journal 24, 46–156.
- 20 Standish E.M. 1990. The observational basis for JPL’s DE 200, the planetary ephemerides of the Astronomical Almanac. A&A 233, 252-271.
- 21 Sussman, G. J. and J. Wisdom 1992. Chaotic Evolution of the Solar System. Science 257, 56–62.
- 22 Varadi F., Bunnegar B., Ghil M. 2003. Successive refinements in long-term integrations of planetary orbits. Astrophysical J. 592, 620-630.
- 23 Will, C. M.: 2006, The Confrontation between General Relativity and Experiment, Living Rev. Relativity 9, (2006), http://www.livingreviews.org/lrr-2006-3
- 24 Wintner, A.: 1947, The analytical foundations of celestial mechanics, Princeton, Princeton Univ Press
- 25 Wisdom, J., Holman, M.: 1991, ‘Symplectic Maps for the N-Body Problem’, Astron. J., 102(4), 1528–1538.