翌日。
老貓和唐躍花了一上午的時間，在一平方公里的區域內再次取了二十個樣本。
崑崙站坐落在火星北半球的伊希地平原中，這是一個直徑一千五百公里的圓形盆地，在三十九億年前被一顆隕石撞了出來，崑崙站的地址之所以設置在此處，是因為火星登陸任務的專家們懷疑伊希地平原中存在大量水冰，崑崙站的後期任務就是尋找並發掘這些潛在的水資源，為進一步的登陸計劃做準備。
下一次火星登陸任務中，地球方面會帶來有能力支持長途行駛的改良型火星車，接下來的科考隊員們就能駕著車子奔馳在廣袤的平原上。
不過很不幸，地球消失了，整個計劃也就戛然而止，到此中斷了。
唐躍在實驗艙內弄了一下午，終於在二十個土壤樣本中找到了合格品。
「這個是到目前為止條件最好的土質，唐躍你看，這個樣本的主要成分是斜長石和高鈣單斜輝石，高矽玻璃的含量較低，含鐵以及鎂矽酸鹽的黏土非常豐富，這種土質和地球上的土壤成分比較接近，PH值比較居中。」麥冬很欣喜。
「能用來種番茄嗎？」唐躍問。
「先做一個Vis-NIR光譜檢測，看看結果。」麥冬指示。
唐躍點點頭，他把土壤樣本浸入溶液中，過濾提取出浸出液，放進光度計裡。
結果很快就出來了。
「嗯……看這個近紅外光譜的數據結果，指標很理想！」麥冬盯著圖譜看了許久，慢慢地點頭，「唐躍，這個樣本是在哪裡採到的？」
唐躍端起玻璃皿看標籤，「E2，是一條古河道的河床上，距離我們不遠，老貓挖了很深才取的樣本。」
「對，伊希地平原是形成於諾亞紀（Noachian）晚期的隕石坑，這個時期有大量的硫酸鹽和水合化合物通過火山噴發沉積在了地表上，你們真厲害，就這個！」麥冬拍板，「唐躍，這個可以用作種植番茄的基質！」
唐躍大喜，轉身朝著車頭大喊:
「老貓！掉頭回去掉頭回去！咱們找到合適的土壤啦！」
話音剛落，車身一震，唐躍的身子也跟著一歪。
老貓扭過頭來，悶悶地說:「叫毛叫啊，車輪陷坑裡啦。」
·
·
·
老貓和唐躍用塑膠密封袋取了十幾公斤的泥土回來，取土的位置距離崑崙站大概一公里左右，是一條早已乾涸的古河道，老貓掃除了表面積累的沙子和礫石，用鏟子往下深刨了一公尺，才找到了這些合格的泥土。
他們用火星流浪狗把泥土運回崑崙站。
唐躍氣喘吁吁地把袋子堆在車庫門前，拍拍巴掌。
「這些泥土起碼有十五六公斤，夠不夠？如果不夠我們可以再去挖。」
「用於初步實驗性的種植應該足夠了。」麥冬說，「不過這些泥土不能直接用來種番茄，它們只是基質，缺乏植物生長所必須的營養成分，所以接下來我們還得給它加入……」
「金坷垃！」唐躍眼睛一亮。
·
·
·
肥料是個很神奇的東西，它是人類自學會刀耕火種之後最重要的發明，透過人工手段改善土壤的性質，提高土壤的肥力，以此極大地提升農作物的產量——在人類漫長的農耕歷史上，學會使用肥料是一個革命性的進步，它是農耕社會得以發展壯大的基礎。
肥料是農業科學中永恆不變的重要課題，社會與科技的進步，特別是化學工業的發展，讓人們可以定量地分析農作物所需的營養成分，並人為地為它們提供這些基本元素，這些技術最終催生出了現代農業中最重要的組成部分，化肥。
現代社會人口得以如此迅速地膨脹，究其根本原因之一是化肥的使用，透過科學地為農作物供給其所必需的礦質元素，讓農業生產迅速工業化，促進糧食產量的逐年攀升，糧食產量的提高是人口增長的基礎。
當然，金坷垃也是一個很神奇的東西。
俗話說肥料摻了金坷垃，一袋能抵兩袋……
崑崙站內很顯然是沒什麼東西可以用來當肥料的，但好在唐躍就是一台有機肥製造機，這大概就是唐躍存在的意義——早在數千年前，中國人的老祖宗們就學會了如何使用有機肥料來提高土壤的肥力，促進莊稼的生長，這大概是根植在中國人骨子裡的先天技能。
千年的農耕文明把對土地的熱愛深深地鐫刻在了每個人的基因中。
要不然怎麼那麼多人喜歡看種田文呢。
老貓和唐躍把廢棄物收集箱內的乾燥糞便全部都翻了出來，這些大便全部經過馬桶的乾燥和抽吸，干硬得像是石頭，用真空包裝包好了——按照規定，崑崙站內產生的一切垃圾都不允許隨意扔在火星表面，必須全部帶回地球。
不過現在這條規定已經作廢。
唐躍想去哪拉屎，就去哪拉屎。
兩人把所有的糞便都帶進了車庫，這種搗屎的活肯定不能在主站內干，否則崑崙站還住不住人了。
唐躍把乾燥糞便倒在車庫地板上，隨意清點了一下，發現他這三個月以來，排便還算均勻順暢，這裡所有的大便都是他自己的，再往前其他人的糞便和垃圾都已經被獵戶座一號帶走了。
老貓蹲下來，手中捏著一根不知道哪裡找來的棍子，饒有趣味地戳了戳地板上包裝好的大便，「唐躍，我覺得你可能嚴重便秘且大便乾燥，你看你拉的翔硬得跟大理石一樣。」
唐躍戴上口罩，並不想搭理老貓這個多話。
老貓還在戳地上的大便，翻過來覆過去地戳。
「唐躍你看，這坨屎像不像一顆真空包裝的茶葉蛋？你是怎麼拉出這麼圓的屎蛋蛋來的？能不能示範一下？」
「還有這個，這坨屎大，我估計一下，起碼得有五兩重吧……」
「這坨很有藝術氣息，看起來像是梵谷的星空。」
「哎唐躍！你來看這個，這坨屎長得很像你欸！簡直就是一個模子裡刻出來的，你們真是一對父子……」
唐躍惱怒地抄起一塊乾燥的大便砸了過去。
------------
對火星軌道變化問題的最後解釋
作者君在作品相關中其實已經解釋過這個問題。
不過仍然有人質疑——「你說得太含糊了」，「火星軌道的變化比你想像要大得多！」
那好吧，既然作者君的簡單解釋不夠有力，那咱們就看看嚴肅的東西，反正這本書寫到現在，嚷嚷著本書BUG一大堆，用國高中物理在書中挑刺的人也不少。
以下是文章內容：
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.
2.4 Error estimation
2.4.1 Relative errors in total energy and angular momentum
According to one of the basic properties of symplectic integrators， which conserve the physically conservative quantities well (total orbital energy and angular momentum)， our long-term numerical integrations seem to have been performed with very small errors. The averaged relative errors of total energy (∼10−9) and of total angular momentum (∼10−11) have remained nearly constant throughout the integration period (Fig. 1). The special startup procedure， warm start， would have reduced the averaged relative error in total energy by about one order of magnitude or more.
Relative numerical error of the total angular momentum δA/A0 and the total energy δE/E0 in our numerical integrationsN± 1，2，3， where δE and δA are the absolute change of the total energy and total angular momentum， respectively， andE0andA0are their initial values. The horizontal unit is Gyr.
Note that different operating systems， different mathematical libraries， and different hardware architectures result in different numerical errors， through the variations in round-off error handling and numerical algorithms. In the upper panel of Fig. 1， we can recognize this situation in the secular numerical error in the total angular momentum， which should be rigorously preserved up to machine-ε precision.
2.4.2 Error in planetary longitudes
Since the symplectic maps preserve total energy and total angular momentum of N-body dynamical systems inherently well， the degree of their preservation may not be a good measure of the accuracy of numerical integrations， especially as a measure of the positional error of planets， i.e. the error in planetary longitudes. To estimate the numerical error in the planetary longitudes， we performed the following procedures. We compared the result of our main long-term integrations with some test integrations， which span much shorter periods but with much higher accuracy than the main integrations. For this purpose， we performed a much more accurate integration with a stepsize of 0.125 d (1/64 of the main integrations) spanning 3 × 105 yr， starting with the same initial conditions as in the N−1 integration. We consider that this test integration provides us with a ‘pseudo-true’ solution of planetary orbital evolution. Next， we compare the test integration with the main integration， N−1. For the period of 3 × 105 yr， we see a difference in mean anomalies of the Earth between the two integrations of ∼0.52°(in the case of the N−1 integration). This difference can be extrapolated to the value ∼8700°， about 25 rotations of Earth after 5 Gyr， since the error of longitudes increases linearly with time in the symplectic map. Similarly， the longitude error of Pluto can be estimated as ∼12°. This value for Pluto is much better than the result in Kinoshita & Nakai (1996) where the difference is estimated as ∼60°.
3 Numerical results – I. Glance at the raw data
In this section we briefly review the long-term stability of planetary orbital motion through some snapshots of raw numerical data. The orbital motion of planets indicates long-term stability in all of our numerical integrations: no orbital crossings nor close encounters between any pair of planets took place.
3.1 General description of the stability of planetary orbits
First， we briefly look at the general character of the long-term stability of planetary orbits. Our interest here focuses particularly on the inner four terrestrial planets for which the orbital time-scales are much shorter than those of the outer five planets. As we can see clearly from the planar orbital configurations shown in Figs 2 and 3， orbital positions of the terrestrial planets differ little between the initial and final part of each numerical integration， which spans several Gyr. The solid lines denoting the present orbits of the planets lie almost within the swarm of dots even in the final part of integrations (b) and (d). This indicates that throughout the entire integration period the almost regular variations of planetary orbital motion remain nearly the same as they are at present.
Vertical view of the four inner planetary orbits (from the z -axis direction) at the initial and final parts of the integrationsN±1. The axes units are au. The xy -plane is set to the invariant plane of Solar system total angular momentum.(a) The initial part ofN+1 ( t = 0 to 0.0547 × 10 9 yr).(b) The final part ofN+1 ( t = 4.9339 × 10 8 to 4.9886 × 10 9 yr).(c) The initial part of N−1 (t= 0 to −0.0547 × 109 yr).(d) The final part ofN−1 ( t =−3.9180 × 10 9 to −3.9727 × 10 9 yr). In each panel， a total of 23 684 points are plotted with an interval of about 2190 yr over 5.47 × 107 yr . Solid lines in each panel denote the present orbits of the four terrestrial planets (taken from DE245).
The variation of eccentricities and orbital inclinations for the inner four planets in the initial and final part of the integration N+1 is shown in Fig. 4. As expected， the character of the variation of planetary orbital elements does not differ significantly between the initial and final part of each integration， at least for Venus， Earth and Mars. The elements of Mercury， especially its eccentricity， seem to change to a significant extent. This is partly because the orbital time-scale of the planet is the shortest of all the planets， which leads to a more rapid orbital evolution than other planets; the innermost planet may be nearest to instability. This result appears to be in some agreement with Laskar's (1994， 1996) expectations that large and irregular variations appear in the eccentricities and inclinations of Mercury on a time-scale of several 109 yr. However， the effect of the possible instability of the orbit of Mercury may not fatally affect the global stability of the whole planetary system owing to the small mass of Mercury. We will mention briefly the long-term orbital evolution of Mercury later in Section 4 using low-pass filtered orbital elements.
The orbital motion of the outer five planets seems rigorously stable and quite regular over this time-span (see also Section 5).
3.2 Time–frequency maps
Although the planetary motion exhibits very long-term stability defined as the non-existence of close encounter events， the chaotic nature of planetary dynamics can change the oscillatory period and amplitude of planetary orbital motion gradually over such long time-spans. Even such slight fluctuations of orbital variation in the frequency domain， particularly in the case of Earth， can potentially have a significant effect on its surface climate system through solar insolation variation (cf. Berger 1988).
To give an overview of the long-term change in periodicity in planetary orbital motion， we performed many fast Fourier transformations (FFTs) along the time axis， and superposed the resulting periodgrams to draw two-dimensional time–frequency maps. The specific approach to drawing these time–frequency maps in this paper is very simple – much simpler than the wavelet analysis or Laskar's (1990， 1993) frequency analysis.
Divide the low-pass filtered orbital data into many fragments of the same length. The length of each data segment should be a multiple of 2 in order to apply the FFT.
Each fragment of the data has a large overlapping part: for example， when the ith data begins from t=ti and ends at t=ti+T， the next data segment ranges from ti+δT≤ti+δT+T， where δT?T. We continue this division until we reach a certain number N by which tn+T reaches the total integration length.
We apply an FFT to each of the data fragments， and obtain n frequency diagrams.
In each frequency diagram obtained above， the strength of periodicity can be replaced by a grey-scale (or colour) chart.
We perform the replacement， and connect all the grey-scale (or colour) charts into one graph for each integration. The horizontal axis of these new graphs should be the time， i.e. the starting times of each fragment of data (ti， where i= 1，…， n). The vertical axis represents the period (or frequency) of the oscillation of orbital elements.
We have adopted an FFT because of its overwhelming speed， since the amount of numerical data to be decomposed into frequency components is terribly huge (several tens of Gbytes).
A typical example of the time–frequency map created by the above procedures is shown in a grey-scale diagram as Fig. 5， which shows the variation of periodicity in the eccentricity and inclination of Earth in N+2 integration. In Fig. 5， the dark area shows that at the time indicated by the value on the abscissa， the periodicity indicated by the ordinate is stronger than in the lighter area around it. We can recognize from this map that the periodicity of the eccentricity and inclination of Earth only changes slightly over the entire period covered by the N+2 integration. This nearly regular trend is qualitatively the same in other integrations and for other planets， although typical frequencies differ planet by planet and element by element.
4.2 Long-term exchange of orbital energy and angular momentum
We calculate very long-periodic variation and exchange of planetary orbital energy and angular momentum using filtered Delaunay elements L， G， H. G and H are equivalent to the planetary orbital angular momentum and its vertical component per unit mass. L is related to the planetary orbital energy E per unit mass as E=−μ2/2L2. If the system is completely linear， the orbital energy and the angular momentum in each frequency bin must be constant. Non-linearity in the planetary system can cause an exchange of energy and angular momentum in the frequency domain. The amplitude of the lowest-frequency oscillation should increase if the system is unstable and breaks down gradually. However， such a symptom of instability is not prominent in our long-term integrations.
In Fig. 7， the total orbital energy and angular momentum of the four inner planets and all nine planets are shown for integration N+2. The upper three panels show the long-periodic variation of total energy (denoted asE- E0)， total angular momentum ( G- G0)， and the vertical component ( H- H0) of the inner four planets calculated from the low-pass filtered Delaunay elements.E0， G0， H0 denote the initial values of each quantity. The absolute difference from the initial values is plotted in the panels. The lower three panels in each figure showE-E0，G-G0 andH-H0 of the total of nine planets. The fluctuation shown in the lower panels is virtually entirely a result of the massive jovian planets.
Comparing the variations of energy and angular momentum of the inner four planets and all nine planets， it is apparent that the amplitudes of those of the inner planets are much smaller than those of all nine planets: the amplitudes of the outer five planets are much larger than those of the inner planets. This does not mean that the inner terrestrial planetary subsystem is more stable than the outer one: this is simply a result of the relative smallness of the masses of the four terrestrial planets compared with those of the outer jovian planets. Another thing we notice is that the inner planetary subsystem may become unstable more rapidly than the outer one because of its shorter orbital time-scales. This can be seen in the panels denoted asinner 4 in Fig. 7 where the longer-periodic and irregular oscillations are more apparent than in the panels denoted astotal 9. Actually， the fluctuations in theinner 4 panels are to a large extent as a result of the orbital variation of the Mercury. However， we cannot neglect the contribution from other terrestrial planets， as we will see in subsequent sections.
4.4 Long-term coupling of several neighbouring planet pairs
Let us see some individual variations of planetary orbital energy and angular momentum expressed by the low-pass filtered Delaunay elements. Figs 10 and 11 show long-term evolution of the orbital energy of each planet and the angular momentum in N+1 and N−2 integrations. We notice that some planets form apparent pairs in terms of orbital energy and angular momentum exchange. In particular， Venus and Earth make a typical pair. In the figures， they show negative correlations in exchange of energy and positive correlations in exchange of angular momentum. The negative correlation in exchange of orbital energy means that the two planets form a closed dynamical system in terms of the orbital energy. The positive correlation in exchange of angular momentum means that the two planets are simultaneously under certain long-term perturbations. Candidates for perturbers are Jupiter and Saturn. Also in Fig. 11， we can see that Mars shows a positive correlation in the angular momentum variation to the Venus–Earth system. Mercury exhibits certain negative correlations in the angular momentum versus the Venus–Earth system， which seems to be a reaction caused by the conservation of angular momentum in the terrestrial planetary subsystem.
It is not clear at the moment why the Venus–Earth pair exhibits a negative correlation in energy exchange and a positive correlation in angular momentum exchange. We may possibly explain this through observing the general fact that there are no secular terms in planetary semimajor axes up to second-order perturbation theories (cf. Brouwer & Clemence 1961; Boccaletti & Pucacco 1998). This means that the planetary orbital energy (which is directly related to the semimajor axis a) might be much less affected by perturbing planets than is the angular momentum exchange (which relates to e). Hence， the eccentricities of Venus and Earth can be disturbed easily by Jupiter and Saturn， which results in a positive correlation in the angular momentum exchange. On the other hand， the semimajor axes of Venus and Earth are less likely to be disturbed by the jovian planets. Thus the energy exchange may be limited only within the Venus–Earth pair， which results in a negative correlation in the exchange of orbital energy in the pair.
As for the outer jovian planetary subsystem， Jupiter–Saturn and Uranus–Neptune seem to make dynamical pairs. However， the strength of their coupling is not as strong compared with that of the Venus–Earth pair.
5 ± 5 × 1010-yr integrations of outer planetary orbits
Since the jovian planetary masses are much larger than the terrestrial planetary masses， we treat the jovian planetary system as an independent planetary system in terms of the study of its dynamical stability. Hence， we added a couple of trial integrations that span ± 5 × 1010 yr， including only the outer five planets (the four jovian planets plus Pluto). The results exhibit the rigorous stability of the outer planetary system over this long time-span. Orbital configurations (Fig. 12)， and variation of eccentricities and inclinations (Fig. 13) show this very long-term stability of the outer five planets in both the time and the frequency domains. Although we do not show maps here， the typical frequency of the orbital oscillation of Pluto and the other outer planets is almost constant during these very long-term integration periods， which is demonstrated in the time–frequency maps on our webpage.
In these two integrations， the relative numerical error in the total energy was ∼10−6 and that of the total angular momentum was ∼10−10.
5.1 Resonances in the Neptune–Pluto system
Kinoshita & Nakai (1996) integrated the outer five planetary orbits over ± 5.5 × 109 yr . They found that four major resonances between Neptune and Pluto are maintained during the whole integration period， and that the resonances may be the main causes of the stability of the orbit of Pluto. The major four resonances found in previous research are as follows. In the following description，λ denotes the mean longitude，Ω is the longitude of the ascending node and ϖ is the longitude of perihelion. Subscripts P and N denote Pluto and Neptune.
Mean motion resonance between Neptune and Pluto (3:2). The critical argument θ1= 3 λP− 2 λN−ϖP librates around 180° with an amplitude of about 80° and a libration period of about 2 × 104 yr.
The argument of perihelion of Pluto ωP=θ2=ϖP−ΩP librates around 90° with a period of about 3.8 × 106 yr. The dominant periodic variations of the eccentricity and inclination of Pluto are synchronized with the libration of its argument of perihelion. This is anticipated in the secular perturbation theory constructed by Kozai (1962).
The longitude of the node of Pluto referred to the longitude of the node of Neptune，θ3=ΩP−ΩN， circulates and the period of this circulation is equal to the period of θ2 libration. When θ3 becomes zero， i.e. the longitudes of ascending nodes of Neptune and Pluto overlap， the inclination of Pluto becomes maximum， the eccentricity becomes minimum and the argument of perihelion becomes 90°. When θ3 becomes 180°， the inclination of Pluto becomes minimum， the eccentricity becomes maximum and the argument of perihelion becomes 90° again. Williams & Benson (1971) anticipated this type of resonance， later confirmed by Milani， Nobili & Carpino (1989).
An argument θ4=ϖP−ϖN+ 3 (ΩP−ΩN) librates around 180° with a long period，∼ 5.7 × 108 yr.
In our numerical integrations， the resonances (i)–(iii) are well maintained， and variation of the critical arguments θ1，θ2，θ3 remain similar during the whole integration period (Figs 14–16 ). However， the fourth resonance (iv) appears to be different: the critical argument θ4 alternates libration and circulation over a 1010-yr time-scale (Fig. 17). This is an interesting fact that Kinoshita & Nakai's (1995， 1996) shorter integrations were not able to disclose.
6 Discussion
What kind of dynamical mechanism maintains this long-term stability of the planetary system? We can immediately think of two major features that may be responsible for the long-term stability. First， there seem to be no significant lower-order resonances (mean motion and secular) between any pair among the nine planets. Jupiter and Saturn are close to a 5:2 mean motion resonance (the famous ‘great inequality’)， but not just in the resonance zone. Higher-order resonances may cause the chaotic nature of the planetary dynamical motion， but they are not so strong as to destroy the stable planetary motion within the lifetime of the real Solar system. The second feature， which we think is more important for the long-term stability of our planetary system， is the difference in dynamical distance between terrestrial and jovian planetary subsystems (Ito & Tanikawa 1999， 2001). When we measure planetary separations by the mutual Hill radii (R_)， separations among terrestrial planets are greater than 26RH， whereas those among jovian planets are less than 14RH. This difference is directly related to the difference between dynamical features of terrestrial and jovian planets. Terrestrial planets have smaller masses， shorter orbital periods and wider dynamical separation. They are strongly perturbed by jovian planets that have larger masses， longer orbital periods and narrower dynamical separation. Jovian planets are not perturbed by any other massive bodies.
The present terrestrial planetary system is still being disturbed by the massive jovian planets. However， the wide separation and mutual interaction among the terrestrial planets renders the disturbance ineffective; the degree of disturbance by jovian planets is O(eJ)(order of magnitude of the eccentricity of Jupiter)， since the disturbance caused by jovian planets is a forced oscillation having an amplitude of O(eJ). Heightening of eccentricity， for example O(eJ)∼0.05， is far from sufficient to provoke instability in the terrestrial planets having such a wide separation as 26RH. Thus we assume that the present wide dynamical separation among terrestrial planets (> 26RH) is probably one of the most significant conditions for maintaining the stability of the planetary system over a 109-yr time-span. Our detailed analysis of the relationship between dynamical distance between planets and the instability time-scale of Solar system planetary motion is now on-going.
Although our numerical integrations span the lifetime of the Solar system， the number of integrations is far from sufficient to fill the initial phase space. It is necessary to perform more and more numerical integrations to confirm and examine in detail the long-term stability of our planetary dynamics.
——以上文段引自 Ito， T.& Tanikawa， K. Long-term integrations and stability of planetary orbits in our Solar System. Mon. Not. R. Astron. Soc. 336， 483–500 (2002)
這只是作者君參考的一篇文章，關於太陽系的穩定性。
還有其他論文，不過也都是英文的，相關課題的中文文獻很少，那些論文下載一篇要九美元（《Nature》真是暴利），作者君寫這篇文章的時候已經回家，不在檢測中心，所以沒有資料庫的使用權，下不起，就不貼上來了。
------------