Long-term integrations and stability of planetary orbits in our Solar system
Abstract
We present the results of very long-term numerical integrations of planetary orbital motions over 109 -yr time-spans including all nine planets. A quick inspection of our numerical data shows that the planetary motion， at least in our simple dynamical model， seems to be quite stable even over this very long time-span. A closer look at the lowest-frequency oscillations using a low-pass filter shows us the potentially diffusive character of terrestrial planetary motion， especially that of Mercury. The behaviour of the eccentricity of Mercury in our integrations is qualitatively similar to the results from Jacques Laskar's secular perturbation theory (e.g. emax～ 0.35 over ～± 4 Gyr). However， there are no apparent secular increases of eccentricity or inclination in any orbital elements of the planets， which may be revealed by still longer-term numerical integrations. We have also performed a couple of trial integrations including motions of the outer five planets over the duration of ± 5 × 1010 yr. The result indicates that the three major resonances in the Neptune–Pluto system have been maintained over the 1011-yr time-span.
1 Introduction
1.1Definition of the problem
The question of the stability of our Solar system has been debated over several hundred years， since the era of Newton. The problem has attracted many famous mathematicians over the years and has played a central role in the development of non-linear dynamics and chaos theory. However， we do not yet have a definite answer to the question of whether our Solar system is stable or not. This is partly a result of the fact that the definition of the term ‘stability’ is vague when it is used in relation to the problem of planetary motion in the Solar system. Actually it is not easy to give a clear， rigorous and physically meaningful definition of the stability of our Solar system.
Among many definitions of stability， here we adopt the Hill definition (Gladman 1993): actually this is not a definition of stability， but of instability. We define a system as becoming unstable when a close encounter occurs somewhere in the system， starting from a certain initial configuration (Chambers， Wetherill & Boss 1996; Ito & Tanikawa 1999). A system is defined as experiencing a close encounter when two bodies approach one another within an area of the larger Hill radius. Otherwise the system is defined as being stable. Henceforward we state that our planetary system is dynamically stable if no close encounter happens during the age of our Solar system， about ±5 Gyr. Incidentally， this definition may be replaced by one in which an occurrence of any orbital crossing between either of a pair of planets takes place. This is because we know from experience that an orbital crossing is very likely to lead to a close encounter in planetary and protoplanetary systems (Yoshinaga， Kokubo & Makino 1999). Of course this statement cannot be simply applied to systems with stable orbital resonances such as the Neptune–Pluto system.
1.2Previous studies and aims of this research
In addition to the vagueness of the concept of stability， the planets in our Solar system show a character typical of dynamical chaos (Sussman & Wisdom 1988， 1992). The cause of this chaotic behaviour is now partly understood as being a result of resonance overlapping (Murray & Holman 1999; Lecar， Franklin & Holman 2001). However， it would require integrating over an ensemble of planetary systems including all nine planets for a period covering several 10 Gyr to thoroughly understand the long-term evolution of planetary orbits， since chaotic dynamical systems are characterized by their strong dependence on initial conditions.
From that point of view， many of the previous long-term numerical integrations included only the outer five planets (Sussman & Wisdom 1988; Kinoshita & Nakai 1996). This is because the orbital periods of the outer planets are so much longer than those of the inner four planets that it is much easier to follow the system for a given integration period. At present， the longest numerical integrations published in journals are those of Duncan & Lissauer (1998). Although their main target was the effect of post-main-sequence solar mass loss on the stability of planetary orbits， they performed many integrations covering up to ～1011 yr of the orbital motions of the four jovian planets. The initial orbital elements and masses of planets are the same as those of our Solar system in Duncan & Lissauer's paper， but they decrease the mass of the Sun gradually in their numerical experiments. This is because they consider the effect of post-main-sequence solar mass loss in the paper. Consequently， they found that the crossing time-scale of planetary orbits， which can be a typical indicator of the instability time-scale， is quite sensitive to the rate of mass decrease of the Sun. When the mass of the Sun is close to its present value， the jovian planets remain stable over 1010 yr， or perhaps longer. Duncan & Lissauer also performed four similar experiments on the orbital motion of seven planets (Venus to Neptune)， which cover a span of ～109 yr. Their experiments on the seven planets are not yet comprehensive， but it seems that the terrestrial planets also remain stable during the integration period， maintaining almost regular oscillations.
On the other hand， in his accurate semi-analytical secular perturbation theory (Laskar 1988)， Laskar finds that large and irregular variations can appear in the eccentricities and inclinations of the terrestrial planets， especially of Mercury and Mars on a time-scale of several 109 yr (Laskar 1996). The results of Laskar's secular perturbation theory should be confirmed and investigated by fully numerical integrations.
In this paper we present preliminary results of six long-term numerical integrations on all nine planetary orbits， covering a span of several 109 yr， and of two other integrations covering a span of ± 5 × 1010 yr. The total elapsed time for all integrations is more than 5 yr， using several dedicated PCs and workstations. One of the fundamental conclusions of our long-term integrations is that Solar system planetary motion seems to be stable in terms of the Hill stability mentioned above， at least over a time-span of ± 4 Gyr. Actually， in our numerical integrations the system was far more stable than what is defined by the Hill stability criterion: not only did no close encounter happen during the integration period， but also all the planetary orbital elements have been confined in a narrow region both in time and frequency domain， though planetary motions are stochastic. Since the purpose of this paper is to exhibit and overview the results of our long-term numerical integrations， we show typical example figures as evidence of the very long-term stability of Solar system planetary motion. For readers who have more specific and deeper interests in our numerical results， we have prepared a webpage (access )， where we show raw orbital elements， their low-pass filtered results， variation of Delaunay elements and angular momentum deficit， and results of our simple time–frequency analysis on all of our integrations.
In Section 2 we briefly explain our dynamical model， numerical method and initial conditions used in our integrations. Section 3 is devoted to a description of the quick results of the numerical integrations. Very long-term stability of Solar system planetary motion is apparent both in planetary positions and orbital elements. A rough estimation of numerical errors is also given. Section 4 goes on to a discussion of the longest-term variation of planetary orbits using a low-pass filter and includes a discussion of angular momentum deficit. In Section 5， we present a set of numerical integrations for the outer five planets that spans ± 5 × 1010 yr. In Section 6 we also discuss the long-term stability of the planetary motion and its possible cause.
2 Description of the numerical integrations
（本部分涉及比较复杂的积分计算，作者君就不贴上来了，贴上来了起点也不一定能成功显示。）
2.3 Numerical method
We utilize a second-order Wisdom–Holman symplectic map as our main integration method (Wisdom & Holman 1991; Kinoshita， Yoshida & Nakai 1991) with a special start-up procedure to reduce the truncation error of angle variables，‘warm start’(Saha & Tremaine 1992， 1994).
The stepsize for the numerical integrations is 8 d throughout all integrations of the nine planets (N±1，2，3)， which is about 1/11 of the orbital period of the innermost planet (Mercury). As for the determination of stepsize， we partly follow the previous numerical integration of all nine planets in Sussman & Wisdom (1988， 7.2 d) and Saha & Tremaine (1994， 225/32 d). We rounded the decimal part of the their stepsizes to 8 to make the stepsize a multiple of 2 in order to reduce the accumulation of round-off error in the computation processes. In relation to this， Wisdom & Holman (1991) performed numerical integrations of the outer five planetary orbits using the symplectic map with a stepsize of 400 d， 1/10.83 of the orbital period of Jupiter. Their result seems to be accurate enough， which partly justifies our method of determining the stepsize. However， since the eccentricity of Jupiter (～0.05) is much smaller than that of Mercury (～0.2)， we need some care when we compare these integrations simply in terms of stepsizes.
In the integration of the outer five planets (F±)， we fixed the stepsize at 400 d.
We adopt Gauss' f and g functions in the symplectic map together with the third-order Halley method (Danby 1992) as a solver for Kepler equations. The number of maximum iterations we set in Halley's method is 15， but they never reached the maximum in any of our integrations.
The interval of the data output is 200 000 d (～547 yr) for the calculations of all nine planets (N±1，2，3)， and about 8000 000 d (～21 903 yr) for the integration of the outer five planets (F±).
Although no output filtering was done when the numerical integrations were in process， we applied a low-pass filter to the raw orbital data after we had completed all the calculations. See Section 4.1 for more detail.
第(1/3)页