Stability of the solar system
| Jacques Laskar (2024), Scholarpedia, 19(4):5216. | doi:10.4249/scholarpedia.5216 | revision #150752 [link to/cite this article] |
Dr. Jacques Laskar accepted the invitation on 25 September 2007 (self-imposed deadline: 31 October 2014).
This article will briefly cover: Historical and current aspects of the problem of the stability of the solar system
Introduction
The problem of the stability of the solar system has fascinated astronomers and mathematicians since antiquity, when it was observed that among the fixed stars, there were also wandering stars, the planets. Efforts were first focused on finding a regularity in the motion of these wanderers, so their movement among the fixed stars could be predicted. For Hipparchus and Ptolemy, the ideal model was a combination of uniform circular motions, the epicycles, which were continually adjusted over the centuries to conform to the observed course of the planets.
From 1609 to 1618, Kepler fixed the planets' trajectories: having assimilated the lessons of Copernicus, he placed the Sun at the center of the universe and, based on the observations of Tycho Brahe, showed that the planets described ellipses around the Sun. At the end of a revolution, each planet found itself back where it started and so retraced the same ellipse. This vision of a perfectly stable solar system in which all orbits were periodic would not remain unchallenged for long.
In 1687 Newton announced the law of universal gravitation. By restricting this law to the interactions of planets with the Sun alone, one obtains Kepler's phenomenology. But Newton's law applies to all interactions: Jupiter is attracted by the Sun, as is Saturn, but Jupiter and Saturn also attract each other. There is no reason to assume that the planets' orbits are fixed invariant ellipses, and Kepler's beautiful regularity is destroyed.
In Newton's view, the perturbations among the planets were strong enough to destroy the stability of the solar system, and divine intervention was required from time to time to restore planets' orbits to their place. Moreover, Newton's law did not yet enjoy its present status, and astronomers wondered if it was truly enough to account for the observed movements of bodies in the solar system.
The problem of solar system stability was a real one, since after Kepler, Halley was able to show, by analyzing the Chaldean observations transmitted by Ptolemy, that Saturn was moving away from the Sun while Jupiter was moving closer. By crudely extrapolating these observations,one finds that six million years ago Jupiter and Saturn were at the same distance from the Sun. In the 18th century, Laplace took up one of these observations, which he dated March 1st, 228 BC: At 4:23 am, mean Paris time, Saturn was observed two fingers under Gamma in Virgo.
The variations of planetary orbits were such that, in order to predict the planets' positions in the sky, de LaLande was required to introduce artificial "secular" terms in his ephemeris tables. Could these terms be accounted for by Newton's law?
Laplace-Lagrange stability of the Solar System
The problem of these discrepancies between computations and observations remained open until the end of the 18th century, when Lagrange and Laplace correctly formulated the equations of motion (for a detailed historical account of this problem and more complete references, see (Laskar, 2013)). Lagrange started from the fact that the motion of a planet remains close, over a short duration, to a Keplerian ellipse, and so could use this ellipse as the basis for a coordinate system (Fig.1). Lagrange then wrote the differential equations that govern the variations in this elliptic motion under the effect of perturbations from other planets, thus inaugurating the methods of classical celestial mechanics. Laplace and Lagrange, whose work converged on this point, calculated secular variations, in other words long-term variations in the planets' semi-major axes under the effects of perturbations by the other planets. Their calculations showed that, up to first order in the masses of the planets, these variations vanish (Poisson, Haretu and Poincaré later showed that this result remains true through second order in the masses of the planets, but not through third order).
This result seemed to contradict Ptolemy's observations from antiquity, but by examining the periodic perturbations between Jupiter and Saturn, Laplace discovered a quasi-resonant term (\(2\lambda_{Jupiter} - 5\lambda_{Saturn}\)) in their longitudes. This term has an amplitude of \(46'50''\) in Saturn's longitude, and a period of about 900 years. This exlpains why observations taken in 228 BC and then in 1590 and 1650 could give the impression of a secular term.
Laplace then calculated many other periodic terms, and established a theory of motion for Jupiter and Saturn in very good agreement with 18th century observations. Above all, using the same theory, he was able to account for Ptolemy's observations to within one minute of arc, without additional terms in his calculations. He thus showed that Newton's law was in itself sufficient to explain the movement of the planets throughout known history, and this exploit no doubt partly accounted for Laplace's determinism.
This result, where Laplace and Lagrange demonstrated that the planets'
semi-major axes undergo only small
oscillations, and do not have secular terms was the first major result
of Stability for the Solar System. At the same time,
Laplace firmly established Newton's law as the universal explanation
for the motion of the celestial bodies.
The problem of the eccentricities
The stability of the semi-major axis of the planets is not sufficient to insure the stability of the Solar System. Indeed, if the eccentricity of the Earth becomes larger than 0.1, and the eccentricity of Mars becomes larger than 0.3, then collisions between these two planets can occur. The problem of the stability of the eccentricities and inclination of the planets was addressed by Laplace and Lagrange in an additional set of papers.
Taking into account only terms of first order in the perturbation series, they showed that the system of equations describing the mean motions of eccentricities and inclinations may be reduced to a system of linear differential equations with constant coefficients depending on the planetary masses and semi-major axes.
\( {d \over dt} \left[\begin{matrix} z_1\\ \vdots\\ z_k \end{matrix}\right] = \sqrt{-1} \left[ \begin{matrix} A_k & 0_k\\ 0_k&B_k \end{matrix} \right] \left[\begin{matrix} z_1\\ \vdots\\ z_k \end{matrix}\right] \)
where for each planet \(j\ ,\) \(z_j=e_j\,\exp\sqrt{-1}\varpi_j , \zeta_j=\sin(i_j/2)\,\exp\sqrt{-1}\Omega_j\ ,\) \(A_k\) and \(B_k\) are \((k,k)\) matrices with real coefficients depending on the values of the planetary masses and semi-major axes. \(0_k\) is the \((k,k)\) zero matrix, and \((a_j,e_j,i_j,\lambda_j, \varpi_j,\Omega_j)\) classical elliptical elements (Fig.1).
Using the conservation of angular momentum, Laplace demonstrated that, provided all the planets rotate around the Sun in the same direction, polynomial or exponential solutions cannot exist for this system. He concluded that the eigenvalues \(g_i, s_i\) of \(A_k\) and \(B_k\) are real, and that the solution for this differential system are quasiperiodic expressions of the form \( z_j = \sum_{i=1}^k \alpha_{ij}\ e^{ig_it} \) \( \zeta_j = \sum_{i=1}^k \beta_{ij}\ e^{is_it} \)
where \(\alpha_{ij}\) and \(\beta_{ij}\) are complex quantities. The frequencies \(g_i, s_i\) are called the secular frequencies of the Solar System, and their values, as computed with a more complete model (Laskar, 1990, Laskar et al, 2004), are given in Table 1.
| perihelion frequencies (arcsec/yr) | node frequencies (arcsec/yr) | ||
|---|---|---|---|
| g_1 | 5.59 | s_1 | -5.59 |
| g_2 | 7.452 | s_2 | -7.05 |
| g_3 | 17.368 | s_3 | -18.85 |
| g_4 | 17.916 | s_4 | -17.755 |
| g_5 | 4.257452 | s_5 | 0 |
| g_6 | 28.245 | s_6 | -26.347855 |
| g_7 | 3.087951 | s_7 | -2.9925259 |
| g_8 | 0.673021 | s_8 | -0.691736 |
| g_9 | -0.34994 | s_9 | -0.34998 |
The inclinations and eccentricities of the orbits
are therefore subject to only small variations about their mean values.
But it must be stressed that Laplace's solutions are very different from
Kepler's, because the orbits are no longer fixed. They are subject to a
double precessional motion with periods ranging from 45,000 to several
million years: precession of the perihelion, which is the slow rotation
of the orbit in its plane, and precession of the nodes, which is the
rotation of the plane of the orbit in space.
Later, Le Verrier, famed for his discovery in 1846 of the planet Neptune through calculations based on observations of irregularities in the movement of Uranus, took up Laplace and Lagrange's calculations and considered the effects of higher order terms in the series (Le Verrier, 1840, 1841). He showed that these terms produced significant corrections and that Laplace's and Lagrange's calculations "could not be used for an indefinite length of time." He then challenged future mathematicians to find exact solutions, without approximations. The difficulty posed by "small divisors" showed that the convergence of the series depended on initial conditions, and the proof of the stability of the solar system remained an open problem (see Laskar, 1992 for more details on this point).
Between 1892 and 1899 Poincaré formulated a negative response to Le Verrier's question. In so doing he rethought the methods of celestial mechanics along the lines of Jacobi's and Hamilton's work. In his memoir On the three body problem and the equations of dynamics, Poincaré showed that it is not possible to integrate the equations of motion of three bodies subject to mutual interaction, and not possible to find an analytic solution representing the movement of the planets valid over an infinite time interval, since the series used by astronomers to calculate the movement of the planets were not convergent.
In the 1950s and 60s, the mathematicians Kolmogorov and Arnold, took up Poincaré's work and showed that, for certain values of the initial conditions, it was nonetheless possible to obtain convergent series. If the masses, eccentricities, and inclinations of the planets are small enough, then many initial conditions lead to quasiperiodic planetary trajectories, similar to the Laplace-Lagrange solutions. But the actual masses of the planets are much too large for this result (known as the KAM theorem) to apply directly to the solar system and thereby prove its stability (In 1966, Michel Henon computed that the masses of the planets needed to be smaller than $10^{-320}$ of the solar mass in order to be able to apply Arnold's theorem on the stability of planetary systems (see Laskar, 2014)).
Although the constants required for the application of Arnold's theorem correspond to extremely small values of the planetary masses, this result reinforced one more time the idea that the Solar System was stable, by any reasonable acceptance of this term, on a time comparable to its age.
The results obtained through numerical integration in the past two decades will show the contrary.
Chaos in the Solar System
In the past decades, the problem of Solar System stability has advanced considerably, due largely to computers which allow extensive analytic calculations and numerical integrations over model time scales approaching the age of the solar system.
One part of these efforts consists of direct numerical integration of the equations of motion (Newton's equations, sometimes with additional relativistic corrections or perturbations due to the Moon). Initial studies were limited to the motion of the outer planets, from Jupiter to Pluto. In fact, the more rapid the orbital movement of a planet, the more difficult it is to numerically integrate its motion. To integrate the orbit of Jupiter with a conventional integrator, a step-size of 40 days will suffice, while a step-size of 0.5 days is required to integrate the motion of the whole solar system (Cohen etal., 1973, Kinoshita and Nakai, 1984, Carpino etal., 1987). These studies, reaching 100 million years essentially confirmed the stability of the system, and the validity of the old perturbative approach of Laplace and Lagrange. At about the same time, calculations of the same system were carried out at MIT over even longer periods, corresponding to times of 210 and 875 million years. These calculations were carried out on "Orrery," a vectorized computer specially designed for the task (Applegate etal., 1986, Sussman and Wisdom, 1988). This latter integration showed that the motion of Pluto is chaotic, showing exponential divergence with respect to the initial conditions, with a characteristic (Lyapunov) time of 20 Ma. But since the mass of Pluto is very small, (1/130,000,000 the mass of the Sun), this does not induce macroscopic instabilities in the rest of the solar system, which appeared relatively stable in these numerical studies.
The other possibility, in order to overcome some of the limitations
of numerical integrations, consists in a semi-analytical approach.
Using perturbation methods
developed by Lagrange, Laplace and Le Verrier, (Laskar, 1989) derived an extended averaged
system for the whole solar system except Pluto, including all
contributions up to second order with respect to the masses, and through
degree 5 in eccentricity and inclination. For the outer planets, some
estimated corrections of third order were also included. The system of
equations thus obtained comprises some 150,000 terms, and does not model
the motion of the planets, but rather the averaged motion of their orbits.
It thus can be integrated numerically on a computer using a very large
step-size, on the order of 500 years. An integration over 200 million
years showed that the solar system, and more particularly the system of
inner planets (Mercury, Venus, Earth, and Mars), is chaotic, with a
Lyapunov time of 5 million years (Laskar, 1989). An error of 15 m
in the Earth's initial position gives rise to an error of about 150 m
after 10 Ma; but this same error grows to 150 million km
after 100 Ma. It is thus possible to construct ephemerides
over a 10 million year period, but it becomes essentially impossible to
predict the motion of the planets with precision beyond 100 million
years.
This chaotic behavior essentially originates in the presence of two secular resonances among the planets: $\theta= 2(g_4-g_3)-(s_4-s_3)$, which is related to Mars and the Earth, and $\sigma = (g_1-g_5)-(s_1-s_2)$, related to Mercury, Venus, and Jupiter (the $g_i$ are the secular frequencies related to the perihelions of the planets, while the $s_i$ are the secular frequencies of the nodes) (Laskar, 1990). The two corresponding arguments change several times from libration to circulation over 200 million years, which is also a characteristic of chaotic behavior.
The improvement of computer speed and the development of new
methods for numerical integration allowed to confirm most of
these results by the direct integration of Newton's equations
(Fig. 3) (Quinn etal., 1991, Laskar
etal., 1992b, Sussman and Wisdom, 1992).
Evolution of planetary orbits.
Over less than one million years, although the motion of the Solar System is chaotic, a quasi periodic model as the one of Laplace-Lagrange (fig. 2) gives a fair representation of the evolution of the planetary orbits. This linear model, although not very precise, provides in particular a good account of the variations of the eccentricity and inclination of the Earth which will be at the origin of the variation of the orientation of its axis of rotation, and thus of the insolation on its surface. Indeed, a similar model derived by Le Verrier (1840, 1841) was used by M. Milankovitch for the establishment of his astronomical theory of paleoclimates.
Over a longer period, of a few millions of years, a quasiperiodic approximation
of the solution is still possible, but it should take into account
the effect of the resonances between the secular motion of inner planets.
It will thus not be possible to obtain it with the classical perturbative method of
Le Verrier and its successors.
On the other hand, such approximation can be obtained by some refined Fourier
techniques, after the numerical integration of the averaged equations (Laskar, 1988, 1990).
The question of the maximum possible variations of the planetary orbits
over the age of the Solar System becomes now even more difficult to answer,
as because of the exponential divergence of the orbits, we know that
it will not be possible to obtain precisely the orbital evolution of the
Solar System after more than 100 Ma.
The computation of the evolution of the Solar system over 5 billion years, may thus appear illusory, but one does not seek here to predict the precise evolution of the system, but to look only for its possible behavior. With this intention, the integration of the orbits was even pushed over durations going well beyond the age of the Solar system (Laskar, 1994, 1995). The results (fig. 4) provide a very clear vision of the stability of the planetary orbits. In this figure is represented the computed evolution of the eccentricity of the orbits of planets of the Solar System over a duration of 25 billion years (from $-10$ to $+15$ billion years).
In fact, for better clarity, the plotted curve represents only the variation of
the maximum eccentricity reached by the planetary orbits over
intervals of 10 million years.
Indeed, the oscillations of the eccentricity resulting from the linear coupling
of the solutions (fig. 3) are removed by this procedure. In doing so, the only variations
which appear in figure 4 are thus the variations due to the chaotic diffusion of the orbits.
For all external planets, the maximum eccentricity is almost constant. That reflects the fact that these trajectories are very close to regular and quasiperiodic trajectories; possible instabilities are insensitive with the scale of the drawing.
For Venus and the Earth, one observes moderated variations, but still significant. The maximum eccentricity of the Earth reached through chaotic diffusion reaches about $0.08$, whereas its current variations are approximately $0.06$. It is about the same for Venus.
The two curves of the maximum eccentricity of the Earth and Venus are very
similar, because of the linear coupling between these two planets. The evolutions
of the orbits of Mars and Mercury are very spectacular. The diffusion of the
eccentricity of Mars can bring this one to $0.2$ in a few billion years, whereas
the variations of the Mercury orbit can lead its eccentricity to values exceeding
$0.5$.
In fact, the system is still constrained by the angular momentum conservation, which constrains
strongly the most massive planets, and it is remarkable to note that in the system of interior
planets, the less one planet is massive, the larger is the possible diffusion of its orbit. The behaviors of
the inclinations are very similar to those of the eccentricities.
Because of the chaotic character of the orbits, a very small modification of the initial conditions
will lead to a solution different from the preceding one after a few hundreds of million years, but
the general aspect of the solutions will undoubtedly remain the same one. To evaluate which are the
possible variations maximum for the orbits of planets on 5 billion years, one can then seek, by very
small modifications of the initial conditions, the trajectory which leads to the strongest
variations of the orbits.
More systematically, Laskar (1994) calculated 5 trajectories of very close initial conditions over
500 Ma. The trajectory leading to the strongest Mercury eccentricity is then retained, and in the
vicinity of this maximum, 5 new trajectories are computed again for a new duration of 500 Ma. For
Mercury, this method then makes it possible to obtain in tens of such stages an orbit which comes to
cut Venus' orbit in less than 3.5 billion years. Let us note that, brought back to the initial
position, because of the exponential divergence of the orbits, the displacements of initial
conditions correspond to a displacement of the position of the Earth smaller than Planck's length
($\approx 10^{-33}$cm).
It should however be noted that to arrive at this possible collision between Mercury and Venus, the
model was used beyond its rigorous field of validity, which does not includes the vicinity of
collisions. In addition, the solution was carefully chosen, so in any case, it is surely not a very
probable one, and the majority of the solutions of close initial conditions will not lead to this
possible collision.
The same method was applied to all other planets, but the chaotic diffusion of their orbits did not
allowed for a collision in less than 5 billion years. The planet which, with Mercury, has the most
unstable orbit is the planet Mars, whose eccentricity can, by this same method reached approximately
$0.25$ in less than 5 billion of years, while the Earth's eccentricity barely reached 0.1.
Marginal stability of the Solar System.
If one summarizes the results obtained with these integrations of the secular equations on a plane graph representing the zone swept by the planetary orbits for the maximum values of their eccentricity, (fig. 5), one notes that the Solar system interns is "full": there is no place for an additional body. One needs at least 3.5 billion years to allow a collision between Mercury and Venus, but an additional body placed in this system will probably collide more rapidly with one of already existing planets.
This observation leads then to the concept of marginal stability for the Solar
system: the Solar system is unstable, but catastrophic phenomena leading to the
destruction of the System in its current form can take place only in a time
comparable with its age, that is to say approximately 5 billion years. The
observation of this present state then makes it possible to suppose that it
always was thus for the Solar system, since the end of its formation. At that
time, it could have remained some other bodies than the current planets, but in
this case, the System would have been much more unstable, and a collision or an
ejection could have taken place (an example could be the impactor of the Earth
which was at the origin of the formation of the Moon). After this event, the
remaining System becomes much more stable. We thus obtain a self-organization of
the System towards increasingly stable states which are always states of marginal
stability.
This vision is in agreement with the models of formation of planets by accretion
of planetesimals (Safronov, 1969), because it shows how the residual bodies could
disappear, in particular in the internal Solar system. It is remarkable that the
zone swept by the orbit of Mars to its maximum eccentricity reaches the limits of
the belt of asteroids.
Concerning the system of the outer planets, the things are appreciably different,
because the direct gravitational short period perturbations are more significant.
The recent numerical simulations show that particles placed among the outer
planets do not remain beyond a few hundreds of million years, apart for some
particular zones of stability or beyond Neptune, in the Kuiper belt, where
objects explicitly were found.
Finally, these observations also make it possible to have an idea of the general
aspect of a planetary system around a star. Indeed, if the process of formation
planetary from planetesimals is correct, it becomes possible that the
planetary systems will always be in a state of marginal stability, like our own
Solar system. At the end of the phase of formation of the system, a large number of bodies can remain, but in this case the system is strongly unstable,
which will
led to a collision or an ejection. After this event, the system becomes more
stable, with constantly, a time of stability comparable with its age.
Planetary collisions in the Solar System
The approach of (Laskar, 1994) allowed extremely fast computations but had some limitations, because the approximation obtained by the averaged equations decreases in accuracy as one approaches the collision. A study using the complete, non-averaged equations was therefore necessary to confirm these results. Actually, because of the chaotic nature of the solutions, A single trajectory will not provide the evolution of the Solar system over more than a few tens of Myr, and a statistical analysis, over an ensemble of solutions is required.
Indeed, before 2009, no direct integration of a single trajectory of the Solar System had yet been published using a realistic model, including the effect of the Moon and general relativity. To approach this problem, (Laskar, 2008) has carried out a statistical study, using the averaged equations (whose numerical integration is about 1000 times faster than for the full equations), for 1000 different solutions that were integrated over 5 Gyr. This study showed that the probability to reach very high eccentricities for Mercury (> 0.6) is on the order of 1%. When neglecting the contribution of general relativity, the same system of averaged equations, gave very surprising results, as in this case, more than half of the trajectories led to an increase of the eccentricity of more than 0.9 in less than 5 Gyr. These results were confirmed by a direct (non averaged) numerical integration, using a symplectic integrator (Laskar et al., 2004), of a pure Newtonian planetary model, for 10 trajectories with close initial conditions. The result was consistent with the results of the secular system since 4 trajectories out of 10 led to eccentricity values for Mercury larger than 0.9 (Laskar, 2008). This large excursion of the eccentricity of Mercury is explained by the presence of a resonance between the perihelion of Mercury and Jupiter, which is made easier in the absence of general relativity (GR). It is known that GR increases the precession speed of the perihelion of mercury by 0.43′′/𝑦𝑟. This moves it from 5.15′′/𝑦𝑟 to 5.58′′/𝑦𝑟, and thus send it further from the value of the perihelion speed of Jupiter (4.25′′/𝑦𝑟”). Independently, Batygin and Laughlin (2008) published similar results. They resumed the calculation of (Laskar, 1994) on a system of non-relativistic equations, and also demonstrated the possibility of collisions between Mercury and Venus.
These results were still incomplete. Indeed, as the relativistic system is much more stable than the non-relativistic system, it is much more difficult to exhibit an orbit of collision between Mercury and Venus in the realistic (relativistic) system than in the non-relativistic system taken into consideration in these two previous numerical studies. The real challenge was therefore in the estimation of the probability of collision of Mercury and Venus for a realistic, relativistic, non averaged model.
Collisions of Mercury, Mars and Venus with the Earth
In order to confirm the results obtained in 1994 and 2008 with the averaged equations, Laskar and Gastineau (2009) have then undertaken a massive computation of orbital solutions for the Solar System motion, using a non-averaged model consistent with the short-term highly accurate INPOP planetary ephemeris that we had developed in the past years (Fienga et al., 2008). Thanks to the installation of the JADE supercomputer at CINES, near Montpellier, with more than 12000 cores, Laskar and Gastineau could benefit of a large amount of computer time during the testing period of this machine. They started the computations as soon as the machine was switched on, in early August 2008, using 2501 cores, with one trajectory being computed on each core. The computations were finalised in about 6 months, totalling about 7 million hours of single node CPU time.
With the JADE machine, they were able to simulate 2501 different solutions of the movement of the planets of the whole Solar System on 5 billion years, corresponding to the life expectancy of the system, before the Sun becomes a red giant. The 2501 computed solutions are all compatible with our current knowledge of the Solar System. They should thus be considered as equiprobable outcomes of the future of the Solar System. In most of the solutions, the trajectories continue to evolve as in the current few millions of years: the planetary orbits are deformed and precess under the influence of the mutual perturbations of the planets but without the possibility of collisions or ejections of planets outside the Solar System. Nevertheless, as predicted by the secular equations, in 1% of the cases, the eccentricity of Mercury increases considerably. In many cases, this deformation of the orbit of Mercury then leads to a collision with Venus, or with the Sun in less than 5 Ga, while the orbit of the Earth remained little affected. However, for one of these orbits, the increase in the eccentricity of Mercury is followed by an increase in the eccentricity of Mars, and a complete internal destabilisation of the inner Solar System (Mercury, Venus, Earth, Mars) in about 3.4 Gyr. Out of 201 additional cases studied in the vicinity of this destabilisation at about 3.4 Gyr, 5 ended by an ejection of Mars out of the Solar System. Others lead to collisions between the planets, or between a planet and the Sun in less than 100 million years. One case resulted in a collision between Mercury and Earth, 29 cases in a collision between Mars and the Earth and 18 in a collision between Venus and the Earth (Laskar and Gastineau, 2009).
Beyond this spectacular aspect, these results also validated the methods of semi-analytical averaging developed for more than 20 years and which had allowed to show the possibility of collision between Mercury and Venus (Laskar, 1994). These results also answer to the question raised more than 300 years ago by Newton, by showing that collisions among planets or ejections are actually possible within the life expectancy of the Sun, that is, in less than 5 Gyr. The main surprise that comes from the numerical simulations of the recent years is that the probability for this catastrophic events to occur is relatively high, of the order of 1%, and thus not just a mathematical curiosity with extremely low probability values. At the same time, 99% of the trajectories will behave in a similar way as in the recent past millions of years, which is coherent with our common understanding that the Solar System has not much evolved in the past 4 Gyr. What is more surprising is that if we consider a pure Newtonian world, starting with the present initial conditions, the probability of collisions within 5 Gyr grows to 60%, which can thus be considered as an additional indirect confirmation of general relativity.
References
Applegate, J.H., Douglas, M.R., Gursel, Y., Sussman, G.J. and Wisdom, J.: 1986, The Solar System for 200 million years. Astron. J., 92, 176–194
Arnold, V.: 1963a, Proof of Kolmogorov’s theorem on the preservation of quasi-periodic motions under small perturbations of the Hamiltonian. Rus. Math. Surv., 18, N6, 9–36
Arnold, V.I.: 1963b, Small denominators and problems of stability of motion in classical celestial mechanics. Russian Math. Surveys, 18, 6, 85–193
Batygin, K., Laughlin, G.: 2008, On the Dynamical Stability of the Solar System. ApJ, 683, 1207–1216
Carpino, M., Milani, A. and Nobili, A.M.: 1987, Long-term numerical integrations and synthetic theories for the motion of the outer planets. Astron. Astrophys., 181, 182–194
Cohen, C.J., Hubbard, E.C., Oesterwinter, C.: 1973, Astron. Papers Am. Ephemeris XXII,1
Fienga, A., Manche, H., Laskar, J., Gastineau, M.: 2008, INPOP06. A new numerical planetary ephemeris, A&A, 477, 315-327
Kinoshita, H., Nakai, H.: 1984, Motions of the perihelion of Neptune and Pluto. Cel. Mech., 34, 203
Kolmogorov, A.N.: 1954, On the conservation of conditionally periodic motions under small perturbation of the Hamiltonian. Dokl. Akad. Nauk. SSSR, 98, 469
Laskar, J. : 1988, Secular evolution of the Solar System over 10 million years, Astron. Astrophys., 198, 341-362
Laskar, J.: 1989, A numerical experiment on the chaotic behavior of the Solar System Nature, 338, 237-238
Laskar, J.: 1990, The chaotic motion of the solar system. A numerical estimate of the size of the chaotic zones, Icarus, 88, 266-291
Laskar, J.: 1992a, La stabilité du Système Solaire, in Chaos et Déteminisme, A. Dahan etal, eds., Seuil, Paris, partially translated and reprinted as: Laskar, J.: 1995 The Stability of the Solar System from Laplace to the Present, in General History of Astronomy, R. Taton et Curtis Wilson eds.,vol. 2B, pp. 240-248
Laskar, J., Quinn, T., Tremaine, S.: 1992b, Confirmation of Resonant Structure in the Solar System. Icarus, 95, 148–152
Laskar, J.: 1994, Large scale chaos in the Solar System. Astron. Astrophys., 287, L9–L12
Laskar, J., 1995. Large scale chaos and Marginal stability of the solar system, the XIth ICMP meeting (Paris July 1994), International Press, pp. 75-120. also in Celestial Mechanics, 64, 115-162
Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levrard, B.: 2004, A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, pp. 261-285
Laskar, J. : 2008: Chaotic diffusion in the Solar System, Icarus, 196, 1-15
Laskar, J., Gastineau, M. : 2009: Existence of collisional trajectories of Mercury, Mars and Venus with the Earth, Nature, 459, 817-819
Laskar, J., 2013, Is the Solar System stable ?, Progress in Mathematical Physics, 66, 2013, pp 239-270
Laskar, J. 2014, Michel Henon and the Stability of the Solar System, http://arxiv.org/abs/1411.4930
LeVerrier U.J.J.: 1840, Mémoire sur les variations séculaires des éléments des orbites pour les sept planètes principales, Mercure, Vénus, la Terre, Mars, Jupiter, Saturne et Uranus. Presented at the Academy of Sciences on September 16, 1839, Additions à la Connaissance des temps pour l’an 1843 Paris, Bachelier, pp. 3–66
LeVerrier U.J.J.: 1841, Mémoire sur les inégalités séculaires des planètes. Presented at the Academy of Sciences on December 14, 1840, Additions à la Connaissance des temps pour l’an 1844 Paris, Bachelier, pp. 28–110
Quinn, T.R., Tremaine, S., Duncan, M.: 1991, 'A three million year integration of the Earth's orbit,' Astron. J. 101, 2287-2305
Safronov, V.S.: 1969, 1972, Evolution of the protoplanetary cloud and formation of the earth and planets., by Safronov, V. S.. Translated from Russian. Jerusalem (Israel) Keter Publishing House, 212 p.
Sussman, G.J., and Wisdom, J.: 1988, 'Numerical evidence that the motion of Pluto is chaotic.' Science 241, 433-437
Sussman, G.J., and Wisdom, J.: 1992, 'Chaotic evolution of the solar system', Science 257, 56-62
