# Stability of the solar system

Dr. Jacques Laskar accepted the invitation on 25 September 2007 (self-imposed deadline: 25 March 2008).

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. 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 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), 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, 1856). 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.

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 decade will show the contrary.

## Recent studies

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., 1992a, Sussman and Wisdom, 1992).

## Evolution of planetary orbits.

Over less than one million years, 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 (1856) 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 reach approximately
$0.25$ in less than 5 billion of years, while the Earth's eccentricity barely reaches 0.1.

## Marginal stability of the Solar System.

If one summarizes these results 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 it can remain a
great number of bodies, but in this case the system is strongly unstable, and is
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.

In particular, a system with only one planet, or even with two planets like
Jupiter and Saturn, will not be able to exist, because this system would be too
stable so that gravitational instabilities could not evacuate the totality of the
other bodies initially present. More precisely, if such a System exists, it will
have to also remain there a multitude of small bodies, which will not have been
evacuated by these gravitational instabilities.

## References

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

Quinn, T.R., Tremaine, S., Duncan, M.: 1991, `A three million year
integration of the Earth's orbit,' *Astron. J.* **101**,
2287-2305

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