# Stability of the solar system

This article has not yet been published; it may contain inaccuracies, unapproved changes, or be unfinished.

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

Figure 1: Elliptical elements. At any given time, a planet ($$J$$) can be considered to move on an elliptical orbit, with semimajor axis $$a$$ and eccentricity $$e\ ,$$ with the sun at one focus ($$O$$). The orientation of this ellipse with respect to a fixed plane $$\Pi\ ,$$ and a direction of reference $$OX\ ,$$ is given by three angles:the inclination $$i\ ,$$ the longitude of the node $$\Omega\ ,$$ the longitude of perihelion $$\varpi=\Omega+\omega\ ,$$ where $$\omega$$ is the argument of perihelion ($$P$$). The position of the planet on this ellipse is given by the mean longitude $$\lambda=M+\varpi\ ,$$ where $$M$$ (mean anomaly) is an angle which is proportional to the area OPJ (third Kepler's law).

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.

Table 1.Fundamental frequencies of the precession motion of the solar system. These values are taken as the mean values over 20 Myr for the inner planets and 50 Myr for the outer planets (laskar et al., 2004). For the inner planets, due to chaotic diffusion, the frequencies can change significantly with time (Laskar, 1990, Laskar et al., 2004). This is why less significant digits are given for the inner planets secular frequencies.
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
Figure 2: The solutions of Laplace-Lagrange for the motion of the planets are combinations of circular and uniform motions with frequencies the precession frequencies $$g_i$$ and $$s_i$$ of the solar system (Table 1). The eccentricity $$e_3$$ of the Earth is given by $$OP\ ,$$ while the inclination of the Earth with respect to the invariant plane of the solar system ($$i_3$$) is $$OQ$$ (Laskar, 1992).

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

Figure 3: The eccentricity of the Earth (a) and Mars (b) during a 6 Ma timespan centered at the present. The solid line is the numerical solution from Quinn etal. (1991), and the dotted line the integration of the secular equations (Laskar, 1990). For clarity, the difference between the two solutions is also plotted (from Laskar, etal., 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$.

Figure 4: Numerical integration of the averaged equations of motion of the solar system 10 Ga backward and 15 Ga forward. For each planet, the maximum value obtained over intervals of 10 Ma for the eccentricity is plotted versus time. For clarity of the figures, Mercury, Venus and the Earth are plotted separately from Mars, Jupiter, Saturn, Uranus and Neptune. The large planets behavior is so regular that all the curves of maximum eccentricity appear as straight lines. On the contrary the corresponding curves of the inner planets show very large and irregular variations, which attest to their diffusion in the chaotic zone.(Laskar, 1994)

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.

Figure 5: Estimates of the zones possibly occupied by the inner planets of the solar system over 5 Ga. The circular orbits correspond to the bold lines, and the zones visited by each planet resulting from the possible increase of eccentricity are the shaded zones. In the case of Mercury and Venus, these shaded zones overlap. Mars can go as far as 1.9 AU, which roughly corresponds to the inner limit of the asteroid belt (Laskar, 1995).

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

Laskar, J., 2013, Is the Solar System stable, Progress in Mathematical Physics Volume 66, 2013, pp 239-270

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