这是用户在 2025-4-14 21:15 为 https://app.immersivetranslate.com/pdf-pro/396260b7-65cf-4f9d-8adc-51d36a69cfe8/ 保存的双语快照页面，由 沉浸式翻译 提供双语支持。了解如何保存？

A long-term numerical solution for the insolation quantities of the Earth
地球日照量的长期数值解

J. Laskar 1 1 ^(1){ }^{1}, P. Robutel 1 1 ^(1){ }^{1}, F. Joutel 1 1 ^(1){ }^{1}, M. Gastineau 1 1 ^(1){ }^{1}, A. C. M. Correia 1 , 2 1 , 2 ^(1,2){ }^{1,2}, and B. Levrard 1 1 ^(1){ }^{1}
J. Laskar 1 1 ^(1){ }^{1} 、P. Robutel 1 1 ^(1){ }^{1} 、F. Joutel 1 1 ^(1){ }^{1} 、M. Gastineau 1 1 ^(1){ }^{1} 、ACM Correia 1 , 2 1 , 2 ^(1,2){ }^{1,2} 和 B. Levrard 1 1 ^(1){ }^{1} 1 Astronomie et Systèmes Dynamiques, IMCCE - CNRS UMR8028, 77 Av. Denfert-Rochereau, 75014 Paris, France e-mail: laskar@imcce.fr
1 天文学和动力系统，IMCCE - CNRS UMR8028，77 Av.丹费尔罗什洛，75014 巴黎，法国电子邮件：laskar@imcce.fr 2 2 ^(2){ }^{2} Departamento de Física da Universidade de Aveiro, Campus Universitário de Santiago, 3810-193 Aveiro, Portugal
2 2 ^(2){ }^{2} 阿威罗大学物理系，圣地亚哥大学校园，3810-193 阿威罗，葡萄牙

Received 23 May 2004 / Accepted 11 August 2004
收到日期：2004 年 5 月 23 日 / 接受日期：2004 年 8 月 11 日

Abstract  抽象的

We present here a new solution for the astronomical computation of the insolation quantities on Earth spanning from - 250 Myr to 250 Myr . This solution has been improved with respect to La93 (Laskar et al. 1993) by using a direct integration of the gravitational equations for the orbital motion, and by improving the dissipative contributions, in particular in the evolution of the Earth-Moon System. The orbital solution has been used for the calibration of the Neogene period (Lourens et al. 2004), and is expected to be used for age calibrations of paleoclimatic data over 40 to 50 Myr, eventually over the full Palaeogene period ( 65 Myr ) with caution. Beyond this time span, the chaotic evolution of the orbits prevents a precise determination of the Earth’s motion. However, the most regular components of the orbital solution could still be used over a much longer time span, which is why we provide here the solution over 250 Myr . Over this time interval, the most striking feature of the obliquity solution, apart from a secular global increase due to tidal dissipation, is a strong decrease of about 0.38 degree in the next few millions of years, due to the crossing of the s 6 + g 5 − g 6 s 6 + g 5 − g 6 s_(6)+g_(5)-g_(6)s_{6}+g_{5}-g_{6} resonance (Laskar et al. 1993). For the calibration of the Mesozoic time scale (about 65 to 250 Myr ), we propose to use the term of largest amplitude in the eccentricity, related to g 2 − g 5 g 2 − g 5 g_(2)-g_(5)g_{2}-g_{5}, with a fixed frequency of 3.200 ′ ′ / yr 3.200 ′ ′ / yr 3.200^('')//yr3.200^{\prime \prime} / \mathrm{yr}, corresponding to a period of 405000 yr . The uncertainty of this time scale over 100 Myr should be about 0.1 % 0.1 % 0.1%0.1 \%, and 0.2 % 0.2 % 0.2%0.2 \% over the full Mesozoic era.
我们在此提出一种新的解决方案，用于天文计算地球日照量，计算时间跨度从 -250 Myr 到 250 Myr。此解决方案相对于 La93（Laskar 等人，1993 年）有所改进，通过使用轨道运动引力方程的直接积分，并改进耗散贡献，特别是在地球-月球系统演化中。轨道解决方案已经用于新近纪的校准（Lourens 等人，2004 年），预计用于 40 到 50 Myr 的古气候数据的年龄校准，最终将谨慎用于整个古近纪（65 Myr）。超过这个时间跨度，轨道的混乱演化会妨碍对地球运动的精确测定。但是，轨道解中最规则的分量仍可用于更长的时间跨度，这就是我们在此提供 2.5 亿年解的原因。在此时间间隔内，除了由于潮汐耗散导致的全球长期增加之外，倾角解最显著的特征是，由于跨越 s 6 + g 5 − g 6 s 6 + g 5 − g 6 s_(6)+g_(5)-g_(6)s_{6}+g_{5}-g_{6} 共振 (Laskar et al. 1993)，在接下来的几百万年中，倾角强烈下降了约 0.38 度。对于中生代时间尺度 (约 65 至 2.5 亿年) 的校准，我们建议使用偏心率中的最大振幅项，与 g 2 − g 5 g 2 − g 5 g_(2)-g_(5)g_{2}-g_{5} 相关，固定频率为 3.200 ′ ′ / yr 3.200 ′ ′ / yr 3.200^('')//yr3.200^{\prime \prime} / \mathrm{yr} ，对应于 405000 年的周期。此时间尺度在 1 亿年内的不确定性应约为 0.1 % 0.1 % 0.1%0.1 \% ，在整个中生代时期为 0.2 % 0.2 % 0.2%0.2 \% 。

Key words. chaos - celestial mechanics - ephemerides - Earth
关键词：混沌 - 天体力学 - 星历表 - 地球

1. Introduction  1. 简介

Due to gravitational planetary perturbations, the elliptical elements of the orbit of the Earth are slowly changing in time, as is the orientation of the planet’s spin axis. These changes induce variations of the insolation received on the Earth’s surface. The first computation of the secular variations of the Earth’s orbital elements were made by Lagrange (1781, 1782), and then Pontécoulant (1834), but it was the work of Agassiz (1840), showing geological evidence of past ice ages, that triggered the search for a correlation between the geological evidence of large climatic changes, and the variations of the Earth’s astronomical parameters. Shortly after, Adhémar (1842) proposed that these climatic variations originated from the precession of the Earth’s rotation axis.
由于行星引力摄动，地球轨道的椭圆根数随时间缓慢变化，行星自转轴的方向亦然。这些变化导致地球表面日照量发生变化。拉格朗（1781 年、1782 年）和庞特库朗（1834 年）首次计算了地球轨道根数的长期变化，但阿加西（1840 年）的工作揭示了过去冰河时代的地质证据，引发了人们对巨大气候变化的地质证据与地球天文参数变化之间关联性的探索。不久之后，阿德马尔（1842 年）提出，这些气候变化源于地球自转轴的进动。
After the publication of a more precise solution of the Earth by Le Verrier (1856), that took into account the secular perturbations of all main planets except for Neptune, Croll (1890) proposed that the variation of the Earth’s eccentricity was also an important parameter for the understanding of the past climates of the Earth.
勒威耶（1856 年）发表了更为精确的地球解决方案，其中考虑了除海王星之外的所有主要行星的长期摄动，克罗尔（1890 年）提出地球偏心率的变化也是了解地球过去气候的重要参数。
The first computations of the variations of the obliquity (angle between the equator and orbital plane) due to the secular variations of the orbital plane of the Earth are due to
由于地球轨道平面的长期变化而导致的倾角（赤道和轨道平面之间的角度）变化的首次计算是由于
Pilgrim (1904), and were later used by Milankovitch (1941) to establish his theory of the Earth’s insolation parameters. Since then, the understanding of the climate response to the orbital forcing has evolved, but all the necessary ingredients for the insolation computations were present in Milankovitch’s work.
Pilgrim (1904) 提出了一些理论，后来被 Milankovitch (1941) 用来建立他的地球日照参数理论。此后，人们对气候对轨道强迫响应的理解不断发展，但计算日照所需的所有要素在 Milankovitch 的著作中都已存在。
The revival of the Milankovitch theory of paleoclimate can be related to the landmark work of Hays et al. (1976), that established a correlation between astronomical forcing and the δ 18 O δ 18 O delta^(18)O\delta^{18} \mathrm{O} records over the past 500 kyr . The Milankovitch theory has since been confirmed overall with variations in the climate response to the insolation forcing (see Imbrie & Imbrie 1979; Imbrie 1982, for more historical details; and Zachos et al. 2001; Grastein et al. 2004, for a recent review on the astronomical calibration of geological data).
米兰科维奇古气候理论的复兴可以与 Hays 等人（1976 年）的里程碑式研究联系起来，他们建立了天文强迫与过去 50 万年 δ 18 O δ 18 O delta^(18)O\delta^{18} \mathrm{O} 记录之间的相关性。米兰科维奇理论此后得到了全面证实，气候对日照强迫的响应存在差异（更多历史细节请参见 Imbrie & Imbrie 1979 年；Imbrie 1982 年；以及 Zachos 等人 2001 年；Grastein 等人 2004 年关于地质数据天文校准的最新综述）。
Since the work of Pilgrim (1904) and Milankovitch (1941), the orbital and precession quantities of the Earth have undergone several improvements. Le Verrier’s solution (1856) consisted in the linearized equations for the secular evolution of the planetary orbits. Stockwell (1873), and Harzer (1895) added the planet Neptune to Le Verrier’s computations. A significant improvement is due to Hill (1897) who discovered that the proximity of a resonance in Jupiter and Saturn’s motion induces some important complements at the second order with respect to the masses. The solution of
自皮尔格里姆（Pilgrim，1904）和米兰科维奇（Milankovitch，1941）的工作以来，地球的轨道和进动量经历了多次改进。勒威耶（L​​e Verrier）的解（1856）包含行星轨道长期演化的线性化方程。斯托克韦尔（Stockwell，1873）和哈泽尔（Harzer，1895）在勒威耶的计算中加入了海王星。希尔（Hill，1897）的发现带来了显著的改进，他发现木星和土星运动中共振的接近性会在质量方面引发一些重要的二阶补项。解
Brouwer & Van Woerkom (1950) is essentially the solution of Le Verrier and Stockwell, complemented by the higher order contributions computed by Hill. This solution was used for insolation computations by Sharav & Boudnikova (1967a,b) with updated values of the parameters, and later on by Vernekar (1972). The computations of Vernekar were actually used by Hays et al. (1976).
Brouwer & Van Woerkom (1950) 的解法本质上是勒威耶和斯托克韦尔的解法，并辅以希尔计算的高阶贡献。Sharav & Boudnikova (1967a,b) 使用该解法并更新了参数值来计算日照量，后来 Vernekar (1972) 也采用了该解法。Vernekar 的计算结果实际上被 Hays 等人 (1976) 所采用。
The next improvement in the computation of the orbital data is by Bretagnon (1974), who computed the terms of second order and degree 3 in eccentricity and inclination in the secular equations, but omitted the terms of degree 5 in the JupiterSaturn system from Hill (1897). This solution was then used by Berger (1978) for the computation of the precession and insolation quantities of the Earth, following Sharav & Boudnikova (1967a,b). All these works assumed implicitly that the motion of the Solar system was regular and that the solutions could be obtained as quasiperiodic series, using perturbation theory.
轨道数据计算的下一个改进来自 Bretagnon (1974)，他计算了长期方程中偏心率和倾角的二阶和三阶项，但省略了 Hill (1897)提出的木星-土星系统中五阶项。随后，Berger (1978) 沿用 Sharav & Boudnikova (1967a,b)的解法，计算了地球的岁差和日照量。所有这些工作都隐含地假设太阳系的运动是规则的，并且可以利用摄动理论以准周期序列的形式获得解。
When Laskar ( 1984 , 1985 , 1986 ) ( 1984 , 1985 , 1986 ) (1984,1985,1986)(1984,1985,1986) computed in an extensive way the secular equations for the Solar system, including all terms up to order 2 in the masses and 5 in eccentricity and inclination, he realized that the traditional perturbative methods could not be used for the integration of the secular equations, due to strong divergences that became apparent in the system of the inner planets (Laskar 1984). This difficulty was overcome by switching to a numerical integration of the secular equations, which could be done in a very effective way, with a very large stepsize of 500 years. These computations provided a much more accurate solution for the orbital motion of the Solar system (Laskar 1986, 1988), which also included a full solution for the precession and obliquity of the Earth for insolation computations over 10 Myr (million of years). Extending his integration to 200 Myr , Laskar ( 1989 , 1990 ) ( 1989 , 1990 ) (1989,1990)(1989,1990) demonstrated that the orbital motion of the planets, and especially of the terrestrial planets, is chaotic, with an exponential divergence corresponding to an increase of the error by a factor of 10 every 10 Myr , thus destroying the hope to obtain a precise astronomical solution for paleoclimate studies over more than a few tens of million of years (Laskar 1999).
当拉斯卡 ( 1984 , 1985 , 1986 ) ( 1984 , 1985 , 1986 ) (1984,1985,1986)(1984,1985,1986) 广泛计算太阳系的长期方程时，包括所有高达 2 阶的质量项以及 5 阶的偏心率和倾角项，他意识到传统的微扰法无法用于长期方程的积分，因为在内行星系统中出现了明显的强烈发散（Laskar 1984）。通过对长期方程进行数值积分，克服了这一困难，这种方法可以非常有效地完成，步长可达 500 年。这些计算为太阳系的轨道运动提供了更精确的解（Laskar 1986, 1988），其中还包括地球进动和倾角的完整解，可用于计算超过 10 百万年的日照。 Laskar ( 1989 , 1990 ) ( 1989 , 1990 ) (1989,1990)(1989,1990) 将他的积分延伸至 2 亿年，证明行星（尤其是类地行星）的轨道运动是混乱的，其指数发散度对应于每 10 亿年误差增加 10 倍，从而摧毁了在几千万年以上的时间里获得古气候研究的精确天文解的希望（Laskar 1999）。
The first long term direct numerical integration (without averaging) of a realistic model of the Solar system, together with the precession and obliquity equations, was made by Quinn et al. (1991) over 3 Myr. Over its range, this solution presented very small differences with the updated secular solution of Laskar et al. (1993) that was computed over 20 Myr , and has since been extensively used for paleoclimate computations under the acronym La93. The orbital motion of the full Solar system has also been computed over 100 Myr by Sussman & Wisdom (1992), using a symplectic integrator with mixed variables (Wisdom & Holman 1991), confirming the chaotic behaviour found by Laskar ( 1989 , 1990 ) ( 1989 , 1990 ) (1989,1990)(1989,1990). Following the improvement of computer technology, long term integrations of realistic models of the Solar system become more easy to perform, but are still challenging when the time interval is the age of the Solar system or if the accuracy of the model is comparable to the precision of short term planetary ephemeris. The longest integration was made over several billions of years by Ito & Tanikawa (2002) with a Newtonian model that did not contain general relativity or lunar contributions, while a recent
Quinn 等人（1991 年）首次对太阳系真实模型以及岁差和轨道倾角方程进行了超过 300 万年的长期直接数值积分（未取平均值）。在其范围内，该解与 Laskar 等人（1993 年）更新的长期解（计算时间超过 2000 万年）差别很小，该长期解后来被广泛用于古气候计算，缩写为 La93。Sussman 和 Wisdom（1992 年）使用混合变量的辛积分器（Wisdom 和 Holman 1991 年）也计算了超过 1 亿年的整个太阳系的轨道运动，证实了 Laskar 发现的混沌行为 ( 1989 , 1990 ) ( 1989 , 1990 ) (1989,1990)(1989,1990) 。随着计算机技术的进步，对太阳系真实模型进行长期积分变得更加容易，但当时间间隔为太阳系的年龄，或模型精度与短期行星历表的精度相当时，仍然具有挑战性。最长的积分是由 Ito & Tanikawa (2002) 使用不包含广义相对论或月球贡献的牛顿模型进行的，跨越数十亿年，而最近的
long term integration of the orbital motion of the Solar system including general relativity, and where the Moon is resolved was made by Varadi et al. (2003) over about 50 Myr.
Varadi 等人（2003）对约 5000 万年的太阳系轨道运动（包括广义相对论）进行了长期积分，并确定了月球的位置。
Here we present the result of the efforts that we have conducted in the past years in our group in order to obtain a new solution for Earth paleoclimate studies over the Neogene period ( ≈ 23 Myr ≈ 23 Myr ~~23Myr\approx 23 \mathrm{Myr} ) and beyond. After a detailed analysis of the main limiting factors for long term integrations (Laskar 1999), and the design of new symplectic integrators (Laskar & Robutel 2001), we have obtained a new numerical solution for both orbital and rotational motion of the Earth that can be used over about 50 Myr for paleoclimate studies, and even over longer periods of time if only the most stable features of the solution are used. Apart from the use of a very complete dynamical model for the orbital motion of the planets, and of the new symplectic integrator of Laskar & Robutel (2001) (see next section), a major change with respect to La93 (Laskar et al. 1993) is the consideration, in the precession solution, of a more complete model for the Earth-Moon tidal interactions (Sect. 4.1).
这里我们展示了我们小组过去几年的努力成果，旨在为新近纪（ ≈ 23 Myr ≈ 23 Myr ~~23Myr\approx 23 \mathrm{Myr} ）及以后的地球古气候研究获得新的解决方案。在详细分析了长期积分的主要限制因素（Laskar 1999）并设计了新的辛积分器（Laskar & Robutel 2001）之后，我们获得了地球轨道和自转运动的新数值解，该解可用于约 50 百万年的古气候研究，如果仅使用解中最稳定的特征，甚至可以用于更长的时间周期。除了使用非常完整的行星轨道运动动力学模型和 Laskar & Robutel (2001) 的新辛积分器（参见下一节）之外，相对于 La93（Laskar 等人，1993 年）的一个主要变化是在岁差解决方案中考虑了更完整的地球-月球潮汐相互作用模型（第 4.1 节）。
The present solution (La2004) or some of its previous variants have already been distributed (Laskar 2001) and used for calibration of sedimentary records over extended time interval (Pälike 2002), and in particular for the construction of a new astronomically calibrated geological time scale for the Neogene period (Lourens et al. 2004). The full solution is available on the Web site www.imcce.fr/ Equipes/ASD/insola/earth/earth.html, together with a set of routines for the computation of the insolation quantities following (Laskar et al. 1993).
本方案（La2004）或其先前的一些变体已发布（Laskar 2001），并用于校准长时间间隔的沉积记录（Pälike 2002），特别是用于构建新近纪新的天文校准地质年代表（Lourens 等人 2004）。完整方案可在网站 www.imcce.fr/Equipes/ASD/insola/earth/earth.html 获取，同时还提供了一套用于计算后续日照量的程序（Laskar 等人 1993）。
In addition to the numerical output of the orbital and rotational parameter of the Earth that is available on the Web site, we have made a special effort to provide in Sect. 8 very compact analytical approximations for the orbital and rotational quantities of the Earth, that can be used in many cases for a better analytical understanding of the insolation variations. Finally, in the last sections, the stability of the solution La2004 is discussed, as well as the possible chaotic transitions of the arguments present in the main secular resonances.
除了网站上提供的地球轨道和自转参数的数值输出外，我们还特别在第 8 节中提供了地球轨道和自转量的非常简洁的解析近似值，这些近似值在许多情况下可用于更好地分析理解日照变化。最后，在最后几节中，我们讨论了 La2004 解的稳定性，以及主要长期共振中可能出现的参数混沌跃迁。

2. Numerical model  2. 数值模型

The orbital solutions La90-93 (Laskar 1990; Laskar et al. 1993) were obtained by a numerical integration of the averaged equations of the Solar system, including the main general relativity and Lunar perturbations. The averaging process was performed using dedicated computer algebra routines. The resulting equations were huge, with about 150000 polynomial terms, but as the short period terms were no longer present, these equations could be integrated with a step size of 200 to 500 years, allowing very extensive long term orbital computations for the Solar system.
La90-93 轨道解（Laskar 1990；Laskar 等人 1993）是通过对太阳系平均方程（包括主要的广义相对论和月球摄动）进行数值积分得到的。平均过程使用专用的计算机代数程序执行。得到的方程非常庞大，包含约 150,000 个多项式项，但由于不再存在短周期项，这些方程可以以 200 到 500 年的步长进行积分，从而可以对太阳系进行非常广泛的长期轨道计算。
Although these averaged equations solutions could be improved by some new adjustment of the initial conditions and parameters (Laskar et al. 2004), it appears that because of the improvements of computer technology, it becomes now possible to obtain more precise results over a few tens of million of years using a more direct numerical integration of the gravitational equations, and this is how the present new
虽然这些平均方程的解可以通过对初始条件和参数进行一些新的调整而得到改进（Laskar et al. 2004），但由于计算机技术的进步，现在可以通过对引力方程进行更直接的数值积分，在几千万年的时间里获得更精确的结果，这就是目前新的
Table 1. Main constants used in La2004. IAU76 refers to the resolutions of the International Astronmical Union of 1976, IERS1992, and IERS2000, refers to the IERS conventions (McCarthy 1992; McCarthy & Petit 2004).
表 1. La2004 中使用的主要常数。IAU76 指 1976 年国际天文学联合会的决议，IERS1992 和 IERS2000 指 IERS 约定（McCarthy 1992；McCarthy & Petit 2004）。
Symbol  象征 Value  价值 Name  姓名 Ref.  参考。
ϵ 0 ϵ 0 epsilon_(0)\epsilon_{0} 84381.448 Obliquity of the ecliptic at J2000.0
黄道倾角 J2000.0		 IAU76
ω 0 ω 0 omega_(0)\omega_{0} 7.292115 × 10 − 5 rad s − 1 7.292115 × 10 − 5 rad s − 1 7.292115 xx10^(-5)rads^(-1)7.292115 \times 10^{-5} \mathrm{rad} \mathrm{s}^{-1} Mean angular velocity of the Earth at J2000.0
J2000.0 时的地球平均角速度		 IERS2000
ψ 0 ψ 0 psi_(0)\psi_{0} 5029.0966 ′ ′ / cy 5029.0966 ′ ′ / cy 5029.0966^('')//cy5029.0966^{\prime \prime} / \mathrm{cy} Precession constant  岁差常数 IERS2000
δ ψ 0 δ ψ 0 deltapsi_(0)\delta \psi_{0} − 0.29965 ′ ′ / cy − 0.29965 ′ ′ / cy -0.29965^('')//cy-0.29965^{\prime \prime} / \mathrm{cy} Correction to the precession
校正岁差		 IERS2000
k 2 k 2 k2k 2 0.305 k 2 k 2 k_(2)k_{2} of the Earth
地球的 k 2 k 2 k_(2)k_{2} 		(Lambeck 1988)  （Lambeck 1988）
k 2 M k 2 M k2_(M)k 2_{\mathrm{M}} 0.0302 k 2 k 2 k_(2)k_{2} of the Moon
月亮的 k 2 k 2 k_(2)k_{2} 		(Yoder 1995)  （Yoder 1995）
Δ t Δ t Delta t\Delta t 639 s  639秒 Time lag of the Earth
地球的时间滞后
Δ t M Δ t M Deltat_(M)\Delta t_{\mathrm{M}} 7055 s  7055秒 Time lag of the Moon
月球的时间滞后
J 2 E J 2 E J_(2E)J_{2 \mathrm{E}} 0.0010826362 Dynamical form-factor for the Earth
地球的动态形状因子		 IERS1992
J 2 M J 2 M J_(2M)J_{2 \mathrm{M}} 0.000202151 Dynamical form-factor for the Moon
月球的动力学形状因子		 IERS1992
J 2 S J 2 S J_(2S)J_{2 \mathrm{~S}} 2 . × 10 − 7 2 . × 10 − 7 2.xx10^(-7)2 . \times 10^{-7} Dynamical form-factor for the Sun
太阳的动力学形状因子		 IERS2000
R E R E R_(E)R_{\mathrm{E}} 6378.1366 km  6378.1366 公里 Equatorial radius of the Earth
地球赤道半径		 IERS2000
R M R M R_(M)R_{\mathrm{M}} 1738 km  1738公里 Equatorial radius of the Moon
月球赤道半径		 IAU76
R sun R sun  R_("sun ")R_{\text {sun }} 696000000 m  696000000 米 Equatorial radius of the Sun
太阳赤道半径		 IAU76
G S G S GSG S 1.32712442076 × 10 20 m 3 s − 2 1.32712442076 × 10 20 m 3 s − 2 1.32712442076 xx10^(20)m^(3)s^(-2)1.32712442076 \times 10^{20} \mathrm{~m}^{3} \mathrm{~s}^{-2} Heliocentric gravitational constant
日心引力常数		 IERS2000
M mer M mer  M_("mer ")M_{\text {mer }} 6023600 Sun - Mercury mass ratio
太阳-水星质量比		 IERS1992
M ven M ven  M_("ven ")M_{\text {ven }} 408523.71 Sun - Venus mass ratio
太阳-金星质量比		 IERS1992
M ear M ear  M_("ear ")M_{\text {ear }} 328900.56 Sun - Earth and Moon mass ratio
太阳-地球和月球的质量比		 IERS1992
M mar M mar  M_("mar ")M_{\text {mar }} 3098708 Sun - Mars mass ratio
太阳-火星质量比		 IERS1992
M jup M jup  M_("jup ")M_{\text {jup }} 1047.3486 Sun - Jupiter and satellites mass ratio
太阳-木星及其卫星的质量比		 IERS1992
M sat M sat  M_("sat ")M_{\text {sat }} 3497.9 Sun - Saturn and satellites mass ratio
太阳-土星及其卫星的质量比		 IERS1992
M ura M ura  M_("ura ")M_{\text {ura }} 22902.94 Sun - Uranus and satellites mass ratio
太阳-天王星及其卫星的质量比		 IERS1992
M nep M nep  M_("nep ")M_{\text {nep }} 19412.24 Sun - Neptune and satellites mass ratio
太阳-海王星及其卫星的质量比		 IERS1992
M plu M plu  M_("plu ")M_{\text {plu }} 135000000 Sun - Pluto and Charon mass ratio
太阳-冥王星和冥卫一的质量比		 IERS1992
μ μ mu\mu 0.0123000383 Moon-Earth mass ratio  月球与地球的质量比 IERS2000
Symbol Value Name Ref. epsilon_(0) 84381.448 Obliquity of the ecliptic at J2000.0 IAU76 omega_(0) 7.292115 xx10^(-5)rads^(-1) Mean angular velocity of the Earth at J2000.0 IERS2000 psi_(0) 5029.0966^('')//cy Precession constant IERS2000 deltapsi_(0) -0.29965^('')//cy Correction to the precession IERS2000 k2 0.305 k_(2) of the Earth (Lambeck 1988) k2_(M) 0.0302 k_(2) of the Moon (Yoder 1995) Delta t 639 s Time lag of the Earth Deltat_(M) 7055 s Time lag of the Moon J_(2E) 0.0010826362 Dynamical form-factor for the Earth IERS1992 J_(2M) 0.000202151 Dynamical form-factor for the Moon IERS1992 J_(2S) 2.xx10^(-7) Dynamical form-factor for the Sun IERS2000 R_(E) 6378.1366 km Equatorial radius of the Earth IERS2000 R_(M) 1738 km Equatorial radius of the Moon IAU76 R_("sun ") 696000000 m Equatorial radius of the Sun IAU76 GS 1.32712442076 xx10^(20)m^(3)s^(-2) Heliocentric gravitational constant IERS2000 M_("mer ") 6023600 Sun - Mercury mass ratio IERS1992 M_("ven ") 408523.71 Sun - Venus mass ratio IERS1992 M_("ear ") 328900.56 Sun - Earth and Moon mass ratio IERS1992 M_("mar ") 3098708 Sun - Mars mass ratio IERS1992 M_("jup ") 1047.3486 Sun - Jupiter and satellites mass ratio IERS1992 M_("sat ") 3497.9 Sun - Saturn and satellites mass ratio IERS1992 M_("ura ") 22902.94 Sun - Uranus and satellites mass ratio IERS1992 M_("nep ") 19412.24 Sun - Neptune and satellites mass ratio IERS1992 M_("plu ") 135000000 Sun - Pluto and Charon mass ratio IERS1992 mu 0.0123000383 Moon-Earth mass ratio IERS2000 | Symbol | Value | Name | Ref. | | :---: | :---: | :--- | ---: | | $\epsilon_{0}$ | 84381.448 | Obliquity of the ecliptic at J2000.0 | IAU76 | | $\omega_{0}$ | $7.292115 \times 10^{-5} \mathrm{rad} \mathrm{s}^{-1}$ | Mean angular velocity of the Earth at J2000.0 | IERS2000 | | $\psi_{0}$ | $5029.0966^{\prime \prime} / \mathrm{cy}$ | Precession constant | IERS2000 | | $\delta \psi_{0}$ | $-0.29965^{\prime \prime} / \mathrm{cy}$ | Correction to the precession | IERS2000 | | $k 2$ | 0.305 | $k_{2}$ of the Earth | (Lambeck 1988) | | $k 2_{\mathrm{M}}$ | 0.0302 | $k_{2}$ of the Moon | (Yoder 1995) | | $\Delta t$ | 639 s | Time lag of the Earth | | | $\Delta t_{\mathrm{M}}$ | 7055 s | Time lag of the Moon | | | $J_{2 \mathrm{E}}$ | 0.0010826362 | Dynamical form-factor for the Earth | IERS1992 | | $J_{2 \mathrm{M}}$ | 0.000202151 | Dynamical form-factor for the Moon | IERS1992 | | $J_{2 \mathrm{~S}}$ | $2 . \times 10^{-7}$ | Dynamical form-factor for the Sun | IERS2000 | | $R_{\mathrm{E}}$ | 6378.1366 km | Equatorial radius of the Earth | IERS2000 | | $R_{\mathrm{M}}$ | 1738 km | Equatorial radius of the Moon | IAU76 | | $R_{\text {sun }}$ | 696000000 m | Equatorial radius of the Sun | IAU76 | | $G S$ | $1.32712442076 \times 10^{20} \mathrm{~m}^{3} \mathrm{~s}^{-2}$ | Heliocentric gravitational constant | IERS2000 | | $M_{\text {mer }}$ | 6023600 | Sun - Mercury mass ratio | IERS1992 | | $M_{\text {ven }}$ | 408523.71 | Sun - Venus mass ratio | IERS1992 | | $M_{\text {ear }}$ | 328900.56 | Sun - Earth and Moon mass ratio | IERS1992 | | $M_{\text {mar }}$ | 3098708 | Sun - Mars mass ratio | IERS1992 | | $M_{\text {jup }}$ | 1047.3486 | Sun - Jupiter and satellites mass ratio | IERS1992 | | $M_{\text {sat }}$ | 3497.9 | Sun - Saturn and satellites mass ratio | IERS1992 | | $M_{\text {ura }}$ | 22902.94 | Sun - Uranus and satellites mass ratio | IERS1992 | | $M_{\text {nep }}$ | 19412.24 | Sun - Neptune and satellites mass ratio | IERS1992 | | $M_{\text {plu }}$ | 135000000 | Sun - Pluto and Charon mass ratio | IERS1992 | | $\mu$ | 0.0123000383 | Moon-Earth mass ratio | IERS2000 | | | | | |
solution La2004 is obtained. A detailed account of the dynamical equations and numerical integrator is not in the scope of the present paper, but will be made in a forthcoming publication. Nevertheless, we will report here on the main issues concerning the numerical integrator, and will more extensively discuss on the dissipative terms in the integration (see Sect. 4.1).
得到了 La2004 的解。本文不详细阐述动力学方程和数值积分器，将在即将出版的出版物中介绍。尽管如此，我们将在此报告有关数值积分器的主要问题，并将更广泛地讨论积分中的耗散项（参见第 4.1 节）。

2.1. Dynamical model  2.1 动力学模型

The orbital model differs from La93, as it comprises now all 9 planets of the Solar System, including Pluto. The postNewtonian general relativity corrections of order 1 / c 2 1 / c 2 1//c^(2)1 / c^{2} due to the Sun are included following Saha & Tremaine (1994).
该轨道模型与 La93 不同，因为它包含了太阳系全部 9 颗行星，包括冥王星。后牛顿广义相对论中，由于太阳的影响， 1 / c 2 1 / c 2 1//c^(2)1 / c^{2} 阶的修正遵循了 Saha & Tremaine (1994) 的计算结果。
The Moon is treated as a separate object. In order to obtain a realistic evolution of the Earth-Moon system, we also take into account the most important coefficient ( J 2 S ) J 2 S (J_(2)^(S))\left(J_{2}^{\mathrm{S}}\right) in the gravitational potential of the Earth and of the Moon (Table 1), and the tidal dissipation in the Earth-Moon System (see Sect. 4.1). We also integrate at the same time the precession and obliquity equations for the Earth and the evolution of its rotation
月球被视为一个独立的天体。为了获得地月系统的真实演化，我们还考虑了地球和月球引力势中最重要的系数 ( J 2 S ) J 2 S (J_(2)^(S))\left(J_{2}^{\mathrm{S}}\right) （表1），以及地月系统中的潮汐耗散（参见第4.1节）。同时，我们还积分了地球的岁差和地轴倾角方程以及地球自转的演化。
period in a comprehensive and coherent way, following the lines of Néron de Surgy & Laskar (1997), Correia et al. (2003) (see Sect. 3).
按照 Néron de Surgy & Laskar (1997)、Correia 等人 (2003) 的研究思路（见第 3 节），以全面、连贯的方式对这一时期进行了阐述。

2.2. Numerical integrator
2.2. 数值积分器

In order to minimize the accumulation of roundoff errors, the numerical integration was performed with the new symplectic integrator scheme S A B A C 4 S A B A C 4 SABAC_(4)S A B A C_{4} of (Laskar & Robutel 2001), with a correction step for the integration of the Moon. This integrator is particularly adapted to perturbed systems where the Hamiltonian governing the equations of motion can be written in the form H = A + ε B H = A + ε B H=A+epsi BH=A+\varepsilon B, as the sum of an integrable part A A AA (the Keplerian equations of the planets orbiting the Sun and of the Moon around the Earth), and a small perturbation potential ε B ε B epsi B\varepsilon B (here the small parameter ϵ ϵ epsilon\epsilon is of the order of the planetary masses). Using this integrator with step size τ τ tau\tau is then equivalent to integrate exactly a close by Hamiltonian H ~ H ~ tilde(H)\tilde{H}, where the error of method H − H ~ H − H ~ H- tilde(H)H-\tilde{H} is of the order of O ( τ 8 ϵ ) + O ( τ 2 ϵ 2 ) O τ 8 ϵ + O τ 2 ϵ 2 O(tau^(8)epsilon)+O(tau^(2)epsilon^(2))O\left(\tau^{8} \epsilon\right)+O\left(\tau^{2} \epsilon^{2}\right), and even O ( τ 8 ϵ ) + O ( τ 4 ϵ 2 ) O τ 8 ϵ + O τ 4 ϵ 2 O(tau^(8)epsilon)+O(tau^(4)epsilon^(2))O\left(\tau^{8} \epsilon\right)+O\left(\tau^{4} \epsilon^{2}\right) when the correction step is added,
为了最小化舍入误差的累积，数值积分采用 (Laskar & Robutel 2001) 提出的新型辛积分方案 S A B A C 4 S A B A C 4 SABAC_(4)S A B A C_{4} 进行，并对月球的积分进行了修正步骤。该积分器特别适用于扰动系统，其中控制运动方程的哈密顿量可以写成 H = A + ε B H = A + ε B H=A+epsi BH=A+\varepsilon B 的形式，即可积部分 A A AA （行星绕太阳运行的开普勒方程和月球绕地球运行的开普勒方程）与小扰动势 ε B ε B epsi B\varepsilon B （此处小参数 ϵ ϵ epsilon\epsilon 为行星质量量级）之和。使用步长为 τ τ tau\tau 的积分器，相当于对一个接近的哈密顿量 H ~ H ~ tilde(H)\tilde{H} 进行积分，其中方法 H − H ~ H − H ~ H- tilde(H)H-\tilde{H} 的误差为 O ( τ 8 ϵ ) + O ( τ 2 ϵ 2 ) O τ 8 ϵ + O τ 2 ϵ 2 O(tau^(8)epsilon)+O(tau^(2)epsilon^(2))O\left(\tau^{8} \epsilon\right)+O\left(\tau^{2} \epsilon^{2}\right) 量级，如果加上修正步骤，甚至会达到 O ( τ 8 ϵ ) + O ( τ 4 ϵ 2 ) O τ 8 ϵ + O τ 4 ϵ 2 O(tau^(8)epsilon)+O(tau^(4)epsilon^(2))O\left(\tau^{8} \epsilon\right)+O\left(\tau^{4} \epsilon^{2}\right) 量级。
while the same quantity is of the order O ( τ 2 ϵ ) O τ 2 ϵ O(tau^(2)epsilon)O\left(\tau^{2} \epsilon\right), in the widely used symplectic integrator of (Wisdom & Holman 1991). More precisely, the integration step S 0 ( τ ) S 0 ( τ ) S_(0)(tau)S_{0}(\tau) with S A B A 4 S A B A 4 SABA_(4)S A B A_{4} is
而在广泛使用的辛积分器中，相同量的阶为 O ( τ 2 ϵ ) O τ 2 ϵ O(tau^(2)epsilon)O\left(\tau^{2} \epsilon\right) (Wisdom & Holman 1991)。更准确地说，积分步长 S 0 ( τ ) S 0 ( τ ) S_(0)(tau)S_{0}(\tau) 与 S A B A 4 S A B A 4 SABA_(4)S A B A_{4} 为
S 0 ( τ ) = e τ c 1 L A e τ d 1 L B e τ c 2 L A e τ d 2 L B e τ c 3 L A × e τ d 2 L B e τ c 2 L A e τ d 1 L B e τ c 1 L A S 0 ( τ ) = e τ c 1 L A e τ d 1 L B e τ c 2 L A e τ d 2 L B e τ c 3 L A × e τ d 2 L B e τ c 2 L A e τ d 1 L B e τ c 1 L A {:[S_(0)(tau)=e^(tauc_(1)L_(A))e^(taud_(1)L_(B))e^(tauc_(2)L_(A))e^(taud_(2)L_(B))e^(tauc_(3)L_(A))],[ xxe^(taud_(2)L_(B))e^(tauc_(2)L_(A))e^(taud_(1)L_(B))e^(tauc_(1)L_(A))]:}\begin{aligned} S_{0}(\tau)= & \mathrm{e}^{\tau c_{1} L_{A}} \mathrm{e}^{\tau d_{1} L_{B}} \mathrm{e}^{\tau c_{2} L_{A}} \mathrm{e}^{\tau d_{2} L_{B}} \mathrm{e}^{\tau c_{3} L_{A}} \\ & \times \mathrm{e}^{\tau d_{2} L_{B}} \mathrm{e}^{\tau c_{2} L_{A}} \mathrm{e}^{\tau d_{1} L_{B}} \mathrm{e}^{\tau c_{1} L_{A}} \end{aligned}
where τ τ tau\tau is the integration step, L f L f L_(f)L_{f} is the differential operator defined by L f g = { f , g } L f g = { f , g } L_(f)g={f,g}L_{f} g=\{f, g\}, where { f , g } { f , g } {f,g}\{f, g\} is the usual Poisson bracket, and the coefficients c i c i c_(i)c_{i} and d i d i d_(i)d_{i} are
其中 τ τ tau\tau 是积分步长， L f L f L_(f)L_{f} 是由 L f g = { f , g } L f g = { f , g } L_(f)g={f,g}L_{f} g=\{f, g\} 定义的微分算子，其中 { f , g } { f , g } {f,g}\{f, g\} 是通常的泊松括号，系数 c i c i c_(i)c_{i} 和 d i d i d_(i)d_{i} 是
c 1 = 1 / 2 − 525 + 70 30 / 70 c 1 = 1 / 2 − 525 + 70 30 / 70 c_(1)=1//2-sqrt(525+70sqrt30)//70c_{1}=1 / 2-\sqrt{525+70 \sqrt{30}} / 70
c 2 = ( 525 + 70 30 − 525 − 70 30 ) / 70 c 2 = ( 525 + 70 30 − 525 − 70 30 ) / 70 c_(2)=(sqrt(525+70sqrt30)-sqrt(525-70sqrt30))//70c_{2}=(\sqrt{525+70 \sqrt{30}}-\sqrt{525-70 \sqrt{30}}) / 70
c 3 = ( 525 − 70 30 ) / 35 c 3 = ( 525 − 70 30 ) / 35 c_(3)=(sqrt(525-70sqrt30))//35c_{3}=(\sqrt{525-70 \sqrt{30}}) / 35
d 1 = 1 / 4 − 30 / 72 d 1 = 1 / 4 − 30 / 72 d_(1)=1//4-sqrt30//72d_{1}=1 / 4-\sqrt{30} / 72
d 2 = 1 / 4 + 30 / 72 d 2 = 1 / 4 + 30 / 72 d_(2)=1//4+sqrt30//72d_{2}=1 / 4+\sqrt{30} / 72.
For the integration of the Earth-Moon system, we have added a correction step (see Laskar & Robutel 2001), and the full integrator steps S 1 ( τ ) S 1 ( τ ) S_(1)(tau)S_{1}(\tau) of S A B A C 4 S A B A C 4 SABAC_(4)S A B A C_{4} then becomes
对于地月系统的积分，我们添加了一个校正步骤（参见 Laskar & Robutel 2001），然后 S A B A C 4 S A B A C 4 SABAC_(4)S A B A C_{4} 的完整积分器步骤 S 1 ( τ ) S 1 ( τ ) S_(1)(tau)S_{1}(\tau) 变为
S 1 ( τ ) = e − τ 3 ε 2 b / 2 L C S 0 ( τ ) e − τ 3 ε 2 b / 2 L C S 1 ( τ ) = e − τ 3 ε 2 b / 2 L C S 0 ( τ ) e − τ 3 ε 2 b / 2 L C S_(1)(tau)=e^(-tau^(3)epsi^(2)b//2L_(C))S_(0)(tau)e^(-tau^(3)epsi^(2)b//2L_(C))S_{1}(\tau)=\mathrm{e}^{-\tau^{3} \varepsilon^{2} b / 2 L_{C}} S_{0}(\tau) \mathrm{e}^{-\tau^{3} \varepsilon^{2} b / 2 L_{C}}
where b = 0.00339677504820860133 b = 0.00339677504820860133 b=0.00339677504820860133b=0.00339677504820860133 and C = { { A , B } , B } C = { { A , B } , B } C={{A,B},B}C=\{\{A, B\}, B\}. The step size used in the integration was τ = 5 × 10 − 3 yr = τ = 5 × 10 − 3 yr = tau=5xx10^(-3)yr=\tau=5 \times 10^{-3} \mathrm{yr}= 1.82625 days. The initial conditions of the integration were least square adjusted to the JPL ephemeris DE406, in order to compensate for small differences in the model. In particular, we do not take into account the effect of the minor planets, and the modelization of the body interactions in the Earth-Moon system is more complete in DE406 (see Williams et al. 2001). The integration time for our complete model, with τ = 5 × 10 − 3 yr τ = 5 × 10 − 3 yr tau=5xx10^(-3)yr\tau=5 \times 10^{-3} \mathrm{yr} is about one day per 5 Myr on a Compaq Alpha (ev68, 833 Mhz ) workstation.
其中 b = 0.00339677504820860133 b = 0.00339677504820860133 b=0.00339677504820860133b=0.00339677504820860133 和 C = { { A , B } , B } C = { { A , B } , B } C={{A,B},B}C=\{\{A, B\}, B\} 。积分中使用的步长为 τ = 5 × 10 − 3 yr = τ = 5 × 10 − 3 yr = tau=5xx10^(-3)yr=\tau=5 \times 10^{-3} \mathrm{yr}= 1.82625 天。为了补偿模型中的细微差异，积分的初始条件已根据 JPL 星历表 DE406 进行最小二乘调整。特别地，我们没有考虑小行星的影响，并且地月系统中天体相互作用的模型在 DE406 中更为完整（参见 Williams 等人，2001 年）。在 Compaq Alpha（ev68，833 Mhz）工作站上，我们完整模型（包含 τ = 5 × 10 − 3 yr τ = 5 × 10 − 3 yr tau=5xx10^(-3)yr\tau=5 \times 10^{-3} \mathrm{yr} ）的积分时间约为每 5 百万年一天。

2.3. Numerical roundoff error
2.3. 数值舍入误差

In Fig. 1a is plotted the evolution of the total energy of the system from -250 Myr to +250 Myr . The energy presents a secular trend that corresponds to the dissipation in the Earth-Moon system. Indeed, after removing the computed energy change due to the secular evolution of the Earth-Moon semi-major axis, the secular trend in the energy evolution disappears and we are left with a residual that is smaller than 2.5 × 10 − 10 2.5 × 10 − 10 2.5 xx10^(-10)2.5 \times 10^{-10} after 250 Myr , and seems to behave as a random walk (Fig. 1b). The normal component of the angular momentum is conserved over the same time with a relative error of less than 1.3 × 10 − 10 1.3 × 10 − 10 1.3 xx10^(-10)1.3 \times 10^{-10} (Fig. 1c).
图 1a 绘制了系统总能量从 -2.5 亿年到 +2.5 亿年的变化。能量呈现出一种长期趋势，与地月系统的耗散相对应。事实上，在去除由于地月半长轴长期演化而计算出的能量变化后，能量演化的长期趋势消失了，我们在 2.5 亿年之后得到的残差小于 2.5 × 10 − 10 2.5 × 10 − 10 2.5 xx10^(-10)2.5 \times 10^{-10} ，似乎表现为随机游动（图 1b）。角动量的法向分量在相同时间内守恒，相对误差小于 1.3 × 10 − 10 1.3 × 10 − 10 1.3 xx10^(-10)1.3 \times 10^{-10} （图 1c）。
In order to test the randomness of the numerical error in the integration, we have plotted the distribution of the difference of total energy from one 1000 year step to the next. After normalization to 1 of the area of the distribution curve (Fig. 2), we obtain a curve that is almost identical to the normalized Gaussian function
为了检验积分中数值误差的随机性，我们绘制了每 1000 年一阶的总能量差值分布图。将分布曲线面积归一化为 1 后（图 2），我们得到了一条与归一化高斯函数几乎相同的曲线。
f ( t ) = 1 σ 2 π exp ⁡ ( − x 2 2 σ 2 ) f ( t ) = 1 σ 2 π exp ⁡ − x 2 2 σ 2 f(t)=(1)/(sigmasqrt(2pi))exp(-(x^(2))/(2sigma^(2)))f(t)=\frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{x^{2}}{2 \sigma^{2}}\right)
Fig. 1. Conservation of integrals. a) Relative variation of the total energy of the system versus time (in Myr) from -250 Myr to +250 Myr . b) The same after correction of the secular trend due to the tidal dissipation in the Earth-Moon system. c) Relative variation of the normal component of the total angular momentum.
图 1。积分守恒。a）系统总能量随时间（以 Myr 为单位）的相对变化，从 -250 Myr 到 +250 Myr。b）在校正由于地月系统潮汐耗散导致的长期趋势后的结果相同。c）总角动量法向分量的相对变化。
Fig. 2. Normalized repartition of the energy numerical error for 1000 yr steps, over the whole integration, from -250 Myr to +250 Myr . Practically superposed to this curve is the computed normal distribution given by Eqs. (4) and (5).
图2. 在整个积分过程中，以1000年为步长，从-2.5亿年到+2.5亿年，对能量数值误差进行归一化重分配。实际上与该曲线叠加的是通过公式(4)和公式(5)计算得出的正态分布。
where σ = 2.6756 × 10 − 13 σ = 2.6756 × 10 − 13 sigma=2.6756 xx10^(-13)\sigma=2.6756 \times 10^{-13} is the computed standard deviation
其中 σ = 2.6756 × 10 − 13 σ = 2.6756 × 10 − 13 sigma=2.6756 xx10^(-13)\sigma=2.6756 \times 10^{-13} 是计算出的标准差
σ = ∑ i ( x i − m ) 2 / ( N − 1 ) σ = ∑ i   x i − m 2 / ( N − 1 ) sigma=sqrt(sum_(i)(x_(i)-m)^(2))//(N-1)\sigma=\sqrt{\sum_{i}\left(x_{i}-m\right)^{2}} /(N-1)
while the mean m = − 2.6413 × 10 − 18 m = − 2.6413 × 10 − 18 m=-2.6413 xx10^(-18)m=-2.6413 \times 10^{-18} is neglected. The integration stepsize is 5 × 10 − 3 yr 5 × 10 − 3 yr 5xx10^(-3)yr5 \times 10^{-3} \mathrm{yr}, the standard deviation per step is thus σ 1 = 5.9828 × 10 − 16 σ 1 = 5.9828 × 10 − 16 sigma_(1)=5.9828 xx10^(-16)\sigma_{1}=5.9828 \times 10^{-16}, that is about 2.7 ε M 2.7 ε M 2.7epsi_(M)2.7 \varepsilon_{\mathrm{M}},
而平均值 m = − 2.6413 × 10 − 18 m = − 2.6413 × 10 − 18 m=-2.6413 xx10^(-18)m=-2.6413 \times 10^{-18} 被忽略。积分步长为 5 × 10 − 3 yr 5 × 10 − 3 yr 5xx10^(-3)yr5 \times 10^{-3} \mathrm{yr} ，因此每步的标准差为 σ 1 = 5.9828 × 10 − 16 σ 1 = 5.9828 × 10 − 16 sigma_(1)=5.9828 xx10^(-16)\sigma_{1}=5.9828 \times 10^{-16} ，即约为 2.7 ε M 2.7 ε M 2.7epsi_(M)2.7 \varepsilon_{\mathrm{M}} ，
Fig. 3. Fundamental planes for the definition of precession and obliquity. E q t E q t E_(q_(t))E_{q_{t}} and E c t E c t E_(c_(t))E_{c_{t}} are the mean equator and ecliptic at date t . E c 0 t . E c 0 t.E_(c_(0))t . E_{c_{0}} is the fixed ecliptic at Julian date J2000, with equinox γ 0 γ 0 gamma_(0)\gamma_{0}. The general precession in longitude ψ ψ psi\psi is defined by ψ = Λ − Ω . ω ψ = Λ − Ω . ω psi=Lambda-Omega.omega\psi=\Lambda-\Omega . \omega is the longitude of the node, and i i ii the inclination. The angle ε ε epsi\varepsilon between E q t E q t E_(q_(t))E_{q_{t}} and E c t E c t E_(c_(t))E_{c_{t}} is the obliquity.
图 3. 定义岁差和黄赤交角的基本平面。 E q t E q t E_(q_(t))E_{q_{t}} 和 E c t E c t E_(c_(t))E_{c_{t}} 分别为平赤道， t . E c 0 t . E c 0 t.E_(c_(0))t . E_{c_{0}} 时刻的黄道面为儒略日 J2000 的固定黄道面，春分点为 γ 0 γ 0 gamma_(0)\gamma_{0} 。经度 ψ ψ psi\psi 上的广义岁差由以下公式定义： ψ = Λ − Ω . ω ψ = Λ − Ω . ω psi=Lambda-Omega.omega\psi=\Lambda-\Omega . \omega 为交点经度， i i ii 为黄赤交角。 E q t E q t E_(q_(t))E_{q_{t}} 和 E c t E c t E_(c_(t))E_{c_{t}} 之间的夹角 ε ε epsi\varepsilon 为黄赤交角。
where ε M ≈ 2.22 × 10 − 16 ε M ≈ 2.22 × 10 − 16 epsi_(M)~~2.22 xx10^(-16)\varepsilon_{M} \approx 2.22 \times 10^{-16} is the machine Epsilon in double precision. Assuming that the machine error is uniformly distributed in the interval [ − ε M / 2 , + ε M / 2 ] − ε M / 2 , + ε M / 2 [-epsi_(M)//2,+epsi_(M)//2]\left[-\varepsilon_{\mathrm{M}} / 2,+\varepsilon_{\mathrm{M}} / 2\right], the standard error of an elementary operation is σ 0 = ε M / 12 σ 0 = ε M / 12 sigma_(0)=epsi_(M)//sqrt12\sigma_{0}=\varepsilon_{\mathrm{M}} / \sqrt{12}, and the error over one step thus corresponds to ( σ 1 / σ 0 ) 2 ≈ 87 σ 1 / σ 0 2 ≈ 87 (sigma_(1)//sigma_(0))^(2)~~87\left(\sigma_{1} / \sigma_{0}\right)^{2} \approx 87 elementary operations. This value is thus extremely small, and we believe that we will not be able to improve σ 1 σ 1 sigma_(1)\sigma_{1} unless we decrease σ 0 σ 0 sigma_(0)\sigma_{0} by switching to higher machine precision. For example, in quadruple precision ( 128 bits for the representation of a real number), ε M ≈ 1.93 × 10 − 34 ε M ≈ 1.93 × 10 − 34 epsi_(M)~~1.93 xx10^(-34)\varepsilon_{\mathrm{M}} \approx 1.93 \times 10^{-34}. Moreover, the Gaussian distribution of the increments of the energy (Fig. 2) ensures that we do not have systematic trends in our computations, and that the numerical errors behave as a random variable with zero mean. The total energy error reflects mostly the behavior of the outer planets, but we can assume that the numerical error will behave in roughly the same for all the planets, although for Mercury, due to its large eccentricity, some slight increase of the error may occur. This should not be the case for the error of method, resulting from the truncation in power of τ τ tau\tau in the integrator, where the error should decrease as the period of the planet increases.
其中 ε M ≈ 2.22 × 10 − 16 ε M ≈ 2.22 × 10 − 16 epsi_(M)~~2.22 xx10^(-16)\varepsilon_{M} \approx 2.22 \times 10^{-16} 是双精度机器的 Epsilon。假设机器误差在区间 [ − ε M / 2 , + ε M / 2 ] − ε M / 2 , + ε M / 2 [-epsi_(M)//2,+epsi_(M)//2]\left[-\varepsilon_{\mathrm{M}} / 2,+\varepsilon_{\mathrm{M}} / 2\right] 内均匀分布，则基本运算的标准误差为 σ 0 = ε M / 12 σ 0 = ε M / 12 sigma_(0)=epsi_(M)//sqrt12\sigma_{0}=\varepsilon_{\mathrm{M}} / \sqrt{12} ，因此一步的误差对应于 ( σ 1 / σ 0 ) 2 ≈ 87 σ 1 / σ 0 2 ≈ 87 (sigma_(1)//sigma_(0))^(2)~~87\left(\sigma_{1} / \sigma_{0}\right)^{2} \approx 87 次基本运算。因此，该值非常小，我们认为除非通过切换到更高的机器精度来降低 σ 0 σ 0 sigma_(0)\sigma_{0} ，否则无法改善 σ 1 σ 1 sigma_(1)\sigma_{1} 。例如，在四倍精度（128 位用于表示实数）下， ε M ≈ 1.93 × 10 − 34 ε M ≈ 1.93 × 10 − 34 epsi_(M)~~1.93 xx10^(-34)\varepsilon_{\mathrm{M}} \approx 1.93 \times 10^{-34} 。此外，能量增量的高斯分布（图 2）确保我们的计算中不存在系统性趋势，并且数值误差表现为均值为零的随机变量。总能量误差主要反映了外行星的行为，但我们可以假设所有行星的数值误差表现大致相同，尽管对于水星来说，由于其偏心率较大，误差可能会略有增加。对于方法误差来说，情况不应该如此，它是由积分器中 τ τ tau\tau 幂的截断引起的，其中误差应该随着行星周期的增加而减少。

3. Precession equations  3.进动方程

We suppose here that the Earth is a homogeneous rigid body with moments of inertia A < B < C A < B < C A < B < CA<B<C and we assume that its spin axis is also the principal axis of inertia. The precession ψ ψ psi\psi and obliquity ε ε epsi\varepsilon (Fig. 3) equations for the rigid Earth in the presence of planetary perturbations are given by Kinoshita (1977), Laskar (1986), Laskar et al. (1993), Néron de Surgy & Laskar (1997).
我们在此假设地球是一个均匀的刚体，其惯性矩为 A < B < C A < B < C A < B < CA<B<C ，并假设其自转轴也是惯性主轴。在行星摄动作用下，刚性地球的进动 ψ ψ psi\psi 和倾角 ε ε epsi\varepsilon （图 3）方程由 Kinoshita（1977）、Laskar（1986）、Laskar 等人（1993）以及 Néron de Surgy 和 Laskar（1997）给出。
{ d X d t = L 1 − X 2 L 2 ( B ( t ) sin ⁡ ψ − A ( t ) cos ⁡ ψ ) d ψ d t = α X L − X L 1 − X 2 L 2 ( A ( t ) sin ⁡ ψ + B ( t ) cos ⁡ ψ ) − 2 C ( t ) d X d t = L 1 − X 2 L 2 ( B ( t ) sin ⁡ ψ − A ( t ) cos ⁡ ψ ) d ψ d t = α X L − X L 1 − X 2 L 2 ( A ( t ) sin ⁡ ψ + B ( t ) cos ⁡ ψ ) − 2 C ( t ) {[(dX)/((d)t)=Lsqrt(1-(X^(2))/(L^(2)))(B(t)sin psi-A(t)cos psi)],[(dpsi)/((d)t)=(alpha X)/(L)-(X)/(Lsqrt(1-(X^(2))/(L^(2))))(A(t)sin psi+B(t)cos psi)-2C(t)]:}\left\{\begin{array}{l} \frac{\mathrm{d} X}{\mathrm{~d} t}=L \sqrt{1-\frac{X^{2}}{L^{2}}}(\mathcal{B}(t) \sin \psi-\mathcal{A}(t) \cos \psi) \\ \frac{\mathrm{d} \psi}{\mathrm{~d} t}=\frac{\alpha X}{L}-\frac{X}{L \sqrt{1-\frac{X^{2}}{L^{2}}}}(\mathcal{A}(t) \sin \psi+\mathcal{B}(t) \cos \psi)-2 C(t) \end{array}\right.
With X = cos ⁡ ε , L = C ω X = cos ⁡ ε , L = C ω X=cos epsi,L=C omegaX=\cos \varepsilon, L=C \omega, where ω ω omega\omega is the spin rate of the Earth, and
其中 X = cos ⁡ ε , L = C ω X = cos ⁡ ε , L = C ω X=cos epsi,L=C omegaX=\cos \varepsilon, L=C \omega 为地球自转速度， ω ω omega\omega 为
{ A ( t ) = 2 1 − p 2 − q 2 [ q ˙ + p ( q p ˙ − p q ˙ ) ] B ( t ) = 2 1 − p 2 − q 2 [ p ˙ − q ( q p ˙ − p q ˙ ) ] C ( t ) = q p ˙ − p q ˙ A ( t ) = 2 1 − p 2 − q 2 [ q ˙ + p ( q p ˙ − p q ˙ ) ] B ( t ) = 2 1 − p 2 − q 2 [ p ˙ − q ( q p ˙ − p q ˙ ) ] C ( t ) = q p ˙ − p q ˙ {[A(t)=(2)/(sqrt(1-p^(2)-q^(2)))[q^(˙)+p(qp^(˙)-pq^(˙))]],[B(t)=(2)/(sqrt(1-p^(2)-q^(2)))[p^(˙)-q(qp^(˙)-pq^(˙))]],[C(t)=qp^(˙)-pq^(˙)]:}\left\{\begin{array}{l}\mathcal{A}(t)=\frac{2}{\sqrt{1-p^{2}-q^{2}}}[\dot{q}+p(q \dot{p}-p \dot{q})] \\ \mathcal{B}(t)=\frac{2}{\sqrt{1-p^{2}-q^{2}}}[\dot{p}-q(q \dot{p}-p \dot{q})] \\ C(t)=q \dot{p}-p \dot{q}\end{array}\right.
where q = sin ⁡ ( i / 2 ) cos ⁡ Ω q = sin ⁡ ( i / 2 ) cos ⁡ Ω q=sin(i//2)cos Omegaq=\sin (i / 2) \cos \Omega and p = sin ⁡ ( i / 2 ) sin ⁡ Ω p = sin ⁡ ( i / 2 ) sin ⁡ Ω p=sin(i//2)sin Omegap=\sin (i / 2) \sin \Omega, and where α α alpha\alpha is the “precession constant”:
其中 q = sin ⁡ ( i / 2 ) cos ⁡ Ω q = sin ⁡ ( i / 2 ) cos ⁡ Ω q=sin(i//2)cos Omegaq=\sin (i / 2) \cos \Omega 和 p = sin ⁡ ( i / 2 ) sin ⁡ Ω p = sin ⁡ ( i / 2 ) sin ⁡ Ω p=sin(i//2)sin Omegap=\sin (i / 2) \sin \Omega ，且 α α alpha\alpha 是“进动常数”：
α = 3 G 2 ω [ m ⊙ ( a ⊙ 1 − e ⊙ 2 ) 3 + m M ( a M 1 − e M 2 ) 3 ( 1 − 3 2 sin 2 ⁡ i M ) ] E d . α = 3 G 2 ω m ⊙ a ⊙ 1 − e ⊙ 2 3 + m M a M 1 − e M 2 3 1 − 3 2 sin 2 ⁡ i M E d . {:[alpha=(3G)/(2omega)[(m_(o.))/((a_(o.)sqrt(1-e_(o.)^(2)))^(3)):}],[{:+(m_(M))/((a_(M)sqrt(1-e_(M)^(2)))^(3))(1-(3)/(2)sin^(2)i_(M))]E_(d).]:}\begin{aligned} \alpha= & \frac{3 G}{2 \omega}\left[\frac{m_{\odot}}{\left(a_{\odot} \sqrt{1-e_{\odot}{ }^{2}}\right)^{3}}\right. \\ & \left.+\frac{m_{\mathrm{M}}}{\left(a_{\mathrm{M}} \sqrt{1-e_{\mathrm{M}}^{2}}\right)^{3}}\left(1-\frac{3}{2} \sin ^{2} i_{\mathrm{M}}\right)\right] E_{\mathrm{d}} . \end{aligned}
For a fast rotating planet like the Earth, the dynamical ellipticity E d = ( 2 C − A − B ) / C E d = ( 2 C − A − B ) / C E_(d)=(2C-A-B)//CE_{\mathrm{d}}=(2 C-A-B) / C can be considered as proportional to ω 2 ω 2 omega^(2)\omega^{2}; this corresponds to the hydrostatic equilibrium (see for example Lambeck 1980). In this approximation, α α alpha\alpha is thus proportional to ω ω omega\omega. The quantities A , B A , B A,B\mathcal{A}, \mathcal{B} and C C CC are related to the secular evolution of the orbital plane of the Earth and are given by the integration of the planetary motions.
对于像地球这样快速旋转的行星，动力学椭圆率 E d = ( 2 C − A − B ) / C E d = ( 2 C − A − B ) / C E_(d)=(2C-A-B)//CE_{\mathrm{d}}=(2 C-A-B) / C 可以被认为与 ω 2 ω 2 omega^(2)\omega^{2} 成正比；这对应于流体静力平衡（例如，参见 Lambeck 1980）。因此，在这个近似中， α α alpha\alpha 与 ω ω omega\omega 成正比。物理量 A , B A , B A,B\mathcal{A}, \mathcal{B} 和 C C CC 与地球轨道平面的长期演化有关，由行星运动的积分给出。

4. Contributions of dissipative effects
4.耗散效应的贡献

4.1. The body tides  4.1. 体潮

4.1.1. Rotational evolution
4.1.1 旋转演化

Following (Darwin 1880; Mignard 1979), we assume that the torque resulting from tidal friction is proportional to the time lag Δ t lag Δ t lag Delta t\operatorname{lag} \Delta t needed for the deformation to reach the equilibrium. This time lag is supposed to be constant, and the angle between the direction of the tide-raising body and the direction of the high tide (which is carried out of the former by the rotation of the Earth) is proportional to the speed of rotation. Such a model is called “viscous”, and corresponds to the case for which 1 / Q 1 / Q 1//Q1 / Q is proportional to the tidal frequency.
遵循（Darwin 1880；Mignard 1979）的论述，我们假设潮汐摩擦产生的扭矩与变形达到平衡所需的时间 lag Δ t lag Δ t lag Delta t\operatorname{lag} \Delta t 成正比。该时间差应为常数，而潮汐体方向与高潮方向（由地球自转将潮汐体引向高潮）之间的夹角与地球自转速度成正比。这种模型被称为“粘性”模型，对应于 1 / Q 1 / Q 1//Q1 / Q 与潮汐频率成正比的情况。
After averaging over the mean anomaly of the Moon and Earth (indices M M MM and ⊕ ⊕ o+\oplus ), and over the longitude of node and perigee of the Moon, we have at second order in eccentricity for the contributions of the solar tides
在对月球和地球的平均异常值（指标 M M MM 和 ⊕ ⊕ o+\oplus ）以及月球节点和近地点经度进行平均后，我们得到了太阳潮汐贡献的二阶偏心率
{ d L d t = − 3 G m ⊙ 2 R 5 k 2 Δ t 2 a ⊙ 6 × [ ( 1 + 15 2 e ⊙ 2 ) ( 1 + X 2 L 2 ) L C − 2 ( 1 + 27 2 e ⊙ 2 ) X L n ⊙ ] d X d t = − 3 G m ⊙ 2 R 5 k 2 Δ t 2 a ⊙ 6 × [ 2 ( 1 + 15 2 e ⊙ 2 ) X C − 2 ( 1 + 27 2 e ⊙ 2 ) n ⊙ ] d L d t = − 3 G m ⊙ 2 R 5 k 2 Δ t 2 a ⊙ 6 × 1 + 15 2 e ⊙ 2 1 + X 2 L 2 L C − 2 1 + 27 2 e ⊙ 2 X L n ⊙ d X d t = − 3 G m ⊙ 2 R 5 k 2 Δ t 2 a ⊙ 6 × 2 1 + 15 2 e ⊙ 2 X C − 2 1 + 27 2 e ⊙ 2 n ⊙ {[(dL)/((d)t)=-(3Gm_(o.)^(2)R^(5)k_(2)Delta t)/(2a_(o.)^(6))],[ xx[(1+(15)/(2)e_(o.)^(2))(1+(X^(2))/(L^(2)))(L)/(C)-2(1+(27)/(2)e_(o.)^(2))(X)/(L)n_(o.)]],[(dX)/((d)t)=-(3Gm_(o.)^(2)R^(5)k_(2)Delta t)/(2a_(o.)^(6))],[ xx[2(1+(15)/(2)e_(o.)^(2))(X)/(C)-2(1+(27)/(2)e_(o.)^(2))n_(o.)]]:}\left\{\begin{aligned} \frac{\mathrm{d} L}{\mathrm{~d} t}= & -\frac{3 G m_{\odot}^{2} R^{5} k_{2} \Delta t}{2 a_{\odot}{ }^{6}} \\ & \times\left[\left(1+\frac{15}{2} e_{\odot}{ }^{2}\right)\left(1+\frac{X^{2}}{L^{2}}\right) \frac{L}{C}-2\left(1+\frac{27}{2} e_{\odot}^{2}\right) \frac{X}{L} n_{\odot}\right] \\ \frac{\mathrm{d} X}{\mathrm{~d} t}= & -\frac{3 G m_{\odot}^{2} R^{5} k_{2} \Delta t}{2 a_{\odot}{ }^{6}} \\ & \times\left[2\left(1+\frac{15}{2} e_{\odot}^{2}\right) \frac{X}{C}-2\left(1+\frac{27}{2} e_{\odot}{ }^{2}\right) n_{\odot}\right] \end{aligned}\right.
where n ⊙ n ⊙ n_(o.)n_{\odot} is the mean motion of the Sun (index ⊙ ⊙ o.\odot ) around the Earth, while the contribution of the lunar tides is
其中 n ⊙ n ⊙ n_(o.)n_{\odot} 是太阳绕地球的平均运动（指数 ⊙ ⊙ o.\odot ），而月球潮汐的贡献是
{ d L d t = − 3 G m M 2 R 5 k 2 Δ t 2 a M 6 × [ 1 2 ( 1 + 15 2 e M 2 ) [ 3 − cos 2 ⁡ i M + ( 3 cos 2 ⁡ i M − 1 ) X 2 L 2 ] L C − 2 ( 1 + 27 2 e M 2 ) X L n M cos ⁡ i M ] d X d t = − 3 G m M 2 R 5 k 2 Δ t 2 a M 6 × [ ( 1 + 15 2 e M 2 ) ( 1 + cos 2 ⁡ i M ) X C − 2 ( 1 + 27 2 e M 2 ) n M cos ⁡ i M ] d L d t = − 3 G m M 2 R 5 k 2 Δ t 2 a M 6 × 1 2 1 + 15 2 e M 2 3 − cos 2 ⁡ i M + 3 cos 2 ⁡ i M − 1 X 2 L 2 L C − 2 1 + 27 2 e M 2 X L n M cos ⁡ i M d X d t = − 3 G m M 2 R 5 k 2 Δ t 2 a M 6 × 1 + 15 2 e M 2 1 + cos 2 ⁡ i M X C − 2 1 + 27 2 e M 2 n M cos ⁡ i M {[(dL)/((d)t)=-(3Gm_(M)^(2)R^(5)k_(2)Delta t)/(2a_(M)^(6))],[xx[(1)/(2)(1+(15)/(2)e_(M)^(2))[3-cos^(2)i_(M)+(3cos^(2)i_(M)-1)(X^(2))/(L^(2))](L)/(C):}],[{:-2(1+(27)/(2)e_(M)^(2))(X)/(L)n_(M)cos i_(M)]],[(dX)/((d)t)=-(3Gm_(M)^(2)R^(5)k_(2)Delta t)/(2a_(M)^(6))],[xx[(1+(15)/(2)e_(M)^(2))(1+cos^(2)i_(M))(X)/(C):}],[{:-2(1+(27)/(2)e_(M)^(2))n_(M)cos i_(M)]]:}\left\{\begin{array}{l} \frac{\mathrm{d} L}{\mathrm{~d} t}=-\frac{3 G m_{\mathrm{M}}^{2} R^{5} k_{2} \Delta t}{2 a_{\mathrm{M}}^{6}} \\ \times\left[\frac{1}{2}\left(1+\frac{15}{2} e_{\mathrm{M}}^{2}\right)\left[3-\cos ^{2} i_{\mathrm{M}}+\left(3 \cos ^{2} i_{\mathrm{M}}-1\right) \frac{X^{2}}{L^{2}}\right] \frac{L}{C}\right. \\ \left.-2\left(1+\frac{27}{2} e_{\mathrm{M}}^{2}\right) \frac{X}{L} n_{\mathrm{M}} \cos i_{\mathrm{M}}\right] \\ \frac{\mathrm{d} X}{\mathrm{~d} t}=-\frac{3 G m_{\mathrm{M}}^{2} R^{5} k_{2} \Delta t}{2 a_{\mathrm{M}}^{6}} \\ \times\left[\left(1+\frac{15}{2} e_{\mathrm{M}}^{2}\right)\left(1+\cos ^{2} i_{\mathrm{M}}\right) \frac{X}{C}\right. \\ \left.-2\left(1+\frac{27}{2} e_{\mathrm{M}}^{2}\right) n_{\mathrm{M}} \cos i_{\mathrm{M}}\right] \end{array}\right.
The “cross tides”, are obtained by the perturbation of the tidal bulge raised by the Sun or the Moon on the other body. Their importance was noticed by Touma & Wisdom (1994). These tidal effects tend to drive the equator towards the ecliptic, and have both the same magnitude
“交叉潮汐”是由太阳或月亮对另一天体产生的潮汐隆起的扰动而产生的。Touma & Wisdom (1994) 注意到了它们的重要性。这些潮汐效应倾向于将赤道推向黄道，并且两者的震级相同
{ d L d t = − 3 G m ⊙ m M R 5 k 2 Δ t 4 a ⊙ 3 a M 3 × ( 1 + 3 2 e ⊙ 2 ) ( 1 + 3 2 e M 2 ) ( 3 cos 2 ⁡ i M − 1 ) ( 1 − X 2 L 2 ) L C d X d t = 0 . d L d t = − 3 G m ⊙ m M R 5 k 2 Δ t 4 a ⊙ 3 a M 3 × 1 + 3 2 e ⊙ 2 1 + 3 2 e M 2 3 cos 2 ⁡ i M − 1 1 − X 2 L 2 L C d X d t = 0 . {[(dL)/((d)t)=-(3Gm_(o.)m_(M)R^(5)k_(2)Delta t)/(4a_(o.)^(3)a_(M)^(3))],[xx(1+(3)/(2)e_(o.)^(2))(1+(3)/(2)e_(M)^(2))(3cos^(2)i_(M)-1)(1-(X^(2))/(L^(2)))(L)/(C)],[((d)X)/((d)t)=0.]:}\left\{\begin{array}{l} \frac{\mathrm{d} L}{\mathrm{~d} t}=-\frac{3 G m_{\odot} m_{\mathrm{M}} R^{5} k_{2} \Delta t}{4 a_{\odot}{ }^{3} a_{\mathrm{M}}{ }^{3}} \\ \times\left(1+\frac{3}{2} e_{\odot}{ }^{2}\right)\left(1+\frac{3}{2} e_{\mathrm{M}}{ }^{2}\right)\left(3 \cos ^{2} i_{\mathrm{M}}-1\right)\left(1-\frac{X^{2}}{L^{2}}\right) \frac{L}{C} \\ \frac{\mathrm{~d} X}{\mathrm{~d} t}=0 . \end{array}\right.

4.1.2. Orbital evolution  4.1.2 轨道演化

At second order in eccentricity, we have for the orbital evolution of the Moon (Mignard 1979, 1980, 1981; Néron de Surgy & Laskar 1997)
在偏心率的二阶下，我们得到月球轨道的演化（Mignard 1979、1980、1981；Néron de Surgy & Laskar 1997）
{ d a M d t = 6 G m M 2 R 5 k 2 Δ t μ a M 7 × [ ( 1 + 27 2 e M 2 ) X C n M cos ⁡ i M − ( 1 + 23 e M 2 ) ] d e M d t = 3 G m M 2 R 5 k 2 Δ t e M μ a M 8 11 2 X C n M cos ⁡ i M − 9 ] dcos i M d t = 3 G m M 2 R 5 k 2 Δ t 2 μ a M 8 ( 1 + 8 e M 2 ) X C n M sin 2 ⁡ i M d a M d t = 6 G m M 2 R 5 k 2 Δ t μ a M 7 × 1 + 27 2 e M 2 X C n M cos ⁡ i M − 1 + 23 e M 2 d e M d t = 3 G m M 2 R 5 k 2 Δ t e M μ a M 8 11 2 X C n M cos ⁡ i M − 9 dcos i M d t = 3 G m M 2 R 5 k 2 Δ t 2 μ a M 8 1 + 8 e M 2 X C n M sin 2 ⁡ i M {[(da_(M))/((d)t)=(6Gm_(M)^(2)R^(5)k_(2)Delta t)/(mua_(M)^(7))],[ xx[(1+(27)/(2)e_(M)^(2))(X)/(Cn_(M))cos i_(M)-(1+23e_(M)^(2))]],[(de_(M))/((d)t)={:(3Gm_(M)^(2)R^(5)k_(2)Delta te_(M))/(mua_(M)^(8))(11)/(2)(X)/(Cn_(M))cos i_(M)-9]],[(dcosi_(M))/((d)t)=(3Gm_(M)^(2)R^(5)k_(2)Delta t)/(2mua_(M)^(8))(1+8e_(M)^(2))(X)/(Cn_(M))sin^(2)i_(M)]:}\left\{\begin{aligned} \frac{\mathrm{d} a_{\mathrm{M}}}{\mathrm{~d} t}= & \frac{6 G m_{\mathrm{M}}^{2} R^{5} k_{2} \Delta t}{\mu a_{\mathrm{M}}^{7}} \\ & \times\left[\left(1+\frac{27}{2} e_{\mathrm{M}}^{2}\right) \frac{X}{C n_{\mathrm{M}}} \cos i_{\mathrm{M}}-\left(1+23 e_{\mathrm{M}}^{2}\right)\right] \\ \frac{\mathrm{d} e_{\mathrm{M}}}{\mathrm{~d} t}= & \left.\frac{3 G m_{\mathrm{M}}^{2} R^{5} k_{2} \Delta t e_{\mathrm{M}}}{\mu a_{\mathrm{M}}^{8}} \frac{11}{2} \frac{X}{C n_{\mathrm{M}}} \cos i_{\mathrm{M}}-9\right] \\ \frac{\mathrm{dcos} i_{\mathrm{M}}}{\mathrm{~d} t}= & \frac{3 G m_{\mathrm{M}}^{2} R^{5} k_{2} \Delta t}{2 \mu a_{\mathrm{M}}^{8}}\left(1+8 e_{\mathrm{M}}^{2}\right) \frac{X}{C n_{\mathrm{M}}} \sin ^{2} i_{\mathrm{M}} \end{aligned}\right.
where μ = m ⊕ m M / ( m ⊕ + m M ) μ = m ⊕ m M / m ⊕ + m M mu=m_(o+)m_(M)//(m_(o+)+m_(M))\mu=m_{\oplus} m_{\mathrm{M}} /\left(m_{\oplus}+m_{\mathrm{M}}\right), and for the Sun
其中 μ = m ⊕ m M / ( m ⊕ + m M ) μ = m ⊕ m M / m ⊕ + m M mu=m_(o+)m_(M)//(m_(o+)+m_(M))\mu=m_{\oplus} m_{\mathrm{M}} /\left(m_{\oplus}+m_{\mathrm{M}}\right) ，对于太阳
{ d a ⊙ d t = 6 G m ⊙ 2 R 5 k 2 Δ t μ ′ a ⊙ 7 × [ ( 1 + 27 2 e ⊙ 2 ) X C n ⊙ − ( 1 + 23 e ⊙ 2 ) ] d e ⊙ d t = 3 G m ⊙ 2 R 5 k 2 Δ t e ⊙ μ ′ a ⊙ 8 [ 11 2 X C n ⊙ − 9 ] d a ⊙ d t = 6 G m ⊙ 2 R 5 k 2 Δ t μ ′ a ⊙ 7 × 1 + 27 2 e ⊙ 2 X C n ⊙ − 1 + 23 e ⊙ 2 d e ⊙ d t = 3 G m ⊙ 2 R 5 k 2 Δ t e ⊙ μ ′ a ⊙ 8 11 2 X C n ⊙ − 9 {[(da_(o.))/(dt)=(6Gm_(o.)^(2)R^(5)k_(2)Delta t)/(mu^(')a_(o.)^(7))],[ xx[(1+(27)/(2)e_(o.)^(2))(X)/(Cn_(o.))-(1+23e_(o.)^(2))]],[(de_(o.))/(dt)=(3Gm_(o.)^(2)R^(5)k_(2)Delta te_(o.))/(mu^(')a_(o.)^(8))[(11)/(2)(X)/(Cn_(o.))-9]]:}\left\{\begin{aligned} \frac{\mathrm{d} a_{\odot}}{\mathrm{d} t}= & \frac{6 G m_{\odot}{ }^{2} R^{5} k_{2} \Delta t}{\mu^{\prime} a_{\odot}{ }^{7}} \\ & \times\left[\left(1+\frac{27}{2} e_{\odot}{ }^{2}\right) \frac{X}{C n_{\odot}}-\left(1+23 e_{\odot}{ }^{2}\right)\right] \\ \frac{\mathrm{d} e_{\odot}}{\mathrm{d} t}= & \frac{3 G m_{\odot}{ }^{2} R^{5} k_{2} \Delta t e_{\odot}}{\mu^{\prime} a_{\odot}{ }^{8}}\left[\frac{11}{2} \frac{X}{C n_{\odot}}-9\right] \end{aligned}\right.
where μ ′ = m ⊕ m ⊙ / ( m ⊕ + m ⊙ ) μ ′ = m ⊕ m ⊙ / m ⊕ + m ⊙ mu^(')=m_(o+)m_(o.)//(m_(o+)+m_(o.))\mu^{\prime}=m_{\oplus} m_{\odot} /\left(m_{\oplus}+m_{\odot}\right). However, both last variations are negligible: about 3 meters per Myr for d a ⊙ / d t d a ⊙ / d t da_(o.)//dt\mathrm{d} a_{\odot} / \mathrm{d} t and 10 − 12 10 − 12 10^(-12)10^{-12} per Myr for d e ⊙ / d t d e ⊙ / d t de_(o.)//dt\mathrm{d} e_{\odot} / \mathrm{d} t.
其中 μ ′ = m ⊕ m ⊙ / ( m ⊕ + m ⊙ ) μ ′ = m ⊕ m ⊙ / m ⊕ + m ⊙ mu^(')=m_(o+)m_(o.)//(m_(o+)+m_(o.))\mu^{\prime}=m_{\oplus} m_{\odot} /\left(m_{\oplus}+m_{\odot}\right) 。然而，最后两种变化都可以忽略不计：对于 d a ⊙ / d t d a ⊙ / d t da_(o.)//dt\mathrm{d} a_{\odot} / \mathrm{d} t ，每百万年约为 3 米；对于 d e ⊙ / d t d e ⊙ / d t de_(o.)//dt\mathrm{d} e_{\odot} / \mathrm{d} t ，每百万年约为 10 − 12 10 − 12 10^(-12)10^{-12} 。

4.1.3. Tides on the Moon
4.1.3 月球上的潮汐

We have also taken into account the tidal effect raised by the Earth on the Moon, following Néron de Surgy & Laskar (1997), with the following simplifying assumptions that the Moon is locked in synchronous spin-orbit resonance with the Earth (i.e. ω M = n M ω M = n M omega_(M)=n_(M)\omega_{\mathrm{M}}=n_{\mathrm{M}} ), and expanding only at first order in the Moon inclination ( 6.41 ∘ 6.41 ∘ 6.41^(@)6.41^{\circ} at t = J 2000 t = J 2000 t=J2000t=\mathrm{J} 2000 ). We thus obtain the additional contributions
我们还考虑了地球对月球产生的潮汐效应，遵循 Néron de Surgy & Laskar (1997) 的研究，并做了以下简化假设：月球与地球锁定在同步自旋轨道共振状态（即 ω M = n M ω M = n M omega_(M)=n_(M)\omega_{\mathrm{M}}=n_{\mathrm{M}} ），且月球倾角仅在一级膨胀（在 t = J 2000 t = J 2000 t=J2000t=\mathrm{J} 2000 处为 6.41 ∘ 6.41 ∘ 6.41^(@)6.41^{\circ} ）。因此，我们得到了额外的贡献
{ d a M d t = − 57 G m ⊕ 2 R M 5 k 2 M Δ t M e M 2 μ a M 7 d e M d t = − 21 G m ⊕ 2 R M 5 k 2 M Δ t M e M 2 μ a M 8 . d a M d t = − 57 G m ⊕ 2 R M 5 k 2 M Δ t M e M 2 μ a M 7 d e M d t = − 21 G m ⊕ 2 R M 5 k 2 M Δ t M e M 2 μ a M 8 . {[(da_(M))/(dt)=-(57 Gm_(o+)^(2)R_(M)^(5)k_(2M)Deltat_(M)e_(M)^(2))/(mua_(M)^(7))],[((d)e_(M))/(dt)=-(21 Gm_(o+)^(2)R_(M)^(5)k_(2M)Deltat_(M)e_(M))/(2mua_(M)^(8)).]:}\left\{\begin{array}{l}\frac{\mathrm{d} a_{\mathrm{M}}}{\mathrm{d} t}=-\frac{57 G m_{\oplus}^{2} R_{\mathrm{M}}{ }^{5} k_{2 \mathrm{M}} \Delta t_{\mathrm{M}} e_{\mathrm{M}}{ }^{2}}{\mu a_{\mathrm{M}}{ }^{7}} \\ \frac{\mathrm{~d} e_{\mathrm{M}}}{\mathrm{d} t}=-\frac{21 G m_{\oplus}{ }^{2} R_{\mathrm{M}}{ }^{5} k_{2 \mathrm{M}} \Delta t_{\mathrm{M}} e_{\mathrm{M}}}{2 \mu a_{\mathrm{M}}{ }^{8}} .\end{array}\right.
With the numerical values from Table 1, these contributions represent 1.2 % 1.2 % 1.2%1.2 \% of the total d a M / d t d a M / d t da_(M)//dt\mathrm{d} a_{\mathrm{M}} / \mathrm{d} t and 30 % 30 % 30%30 \% of the total d e M / d t d e M / d t de_(M)//dt\mathrm{d} e_{\mathrm{M}} / \mathrm{d} t for present conditions. Finally, the tides raised on the Moon by the Sun can be neglected because the ratio of the magnitudes of solar and terrestrial tides on the Moon is
根据表 1 中的数值，在当前条件下，这些贡献占总潮汐力的 1.2 % 1.2 % 1.2%1.2 \% ，占总潮汐力的 30 % 30 % 30%30 \% 。最后，太阳在月球上引起的潮汐可以忽略不计，因为月球上太阳潮汐力和地球潮汐力的比值为
( m ⊙ m ⊕ ) 2 ( a ⊙ a M ) 6 ≃ 3.2 × 10 − 5 m ⊙ m ⊕ 2 a ⊙ a M 6 ≃ 3.2 × 10 − 5 ((m_(o.))/(m_(o+)))^(2)((a_(o.))/(a_(M)))^(6)≃3.2 xx10^(-5)\left(\frac{m_{\odot}}{m_{\oplus}}\right)^{2}\left(\frac{a_{\odot}}{a_{\mathrm{M}}}\right)^{6} \simeq 3.2 \times 10^{-5}.

4.2. Other dissipative effects
4.2 其他耗散效应

4.2.1. Core-mantle friction
4.2.1 核幔摩擦

The core and the mantle have different dynamical ellipticities, so they tend to have different precession rates. This trend produces a viscous friction at the core-mantle boundary (Poincaré 1910; Rochester 1976; Lumb & Aldridge 1991). We will consider here that an effective viscosity v v vv can account for a weak laminar friction, as well as a strong turbulent one which thickens the boundary layer (e.g. Correia et al. 2003).
地核和地幔具有不同的动力学椭圆率，因此它们往往具有不同的进动速率。这种趋势在地核-地幔边界产生了粘性摩擦（Poincaré 1910；Rochester 1976；Lumb & Aldridge 1991）。我们在此考虑有效粘性 v v vv 既可以解释弱层流摩擦，也可以解释强湍流摩擦（后者会使边界层增厚）（例如 Correia 等人 2003）。
The core-mantle friction tends to slow down the rotation and to bring the obliquity down to 0 ∘ 0 ∘ 0^(@)0^{\circ} if ε < 90 ∘ ε < 90 ∘ epsi < 90^(@)\varepsilon<90^{\circ} and up to 180 ∘ 180 ∘ 180^(@)180^{\circ} otherwise, which contrasts with the effect of the tides. Furthermore, one can see that, despite the strong coupling by pressure forces, there can be a substantial contribution to the variation of the spin for high viscosities and moderate speeds of rotation.
核幔摩擦往往会减缓自转，并使倾角在 ε < 90 ∘ ε < 90 ∘ epsi < 90^(@)\varepsilon<90^{\circ} 时降至 0 ∘ 0 ∘ 0^(@)0^{\circ} ，否则升高至 180 ∘ 180 ∘ 180^(@)180^{\circ} ，这与潮汐的影响形成对比。此外，我们可以看出，尽管压力作用有很强的耦合作用，但在高粘度和中等旋转速度下，自旋的变化仍会受到相当大的影响。

4.2.2. Atmospheric tides  4.2.2 大气潮汐

The Earth’s atmosphere also undergoes some torques which can be transmitted to the surface by friction: a torque caused by the gravitational tides raised by the Moon and the Sun, a magnetic one generated by interactions between the magnetosphere and the solar wind. Both effects are negligible; see respectively Chapman & Lindzen (1970) and Volland (1978).
地球大气层也会受到一些扭矩的影响，这些扭矩可以通过摩擦传递到地表：由月球和太阳引起的引力潮汐力引起的扭矩，以及由磁层和太阳风相互作用产生的磁力扭矩。这两种影响都可以忽略不计；分别参见 Chapman & Lindzen (1970) 和 Volland (1978)。
Finally, a torque is produced by the daily solar heating which induces a redistribution of the air pressure, mainly driven by a semidiurnal wave, hence the so-called thermal atmospheric tides (Chapman & Lindzen 1970). The axis of symmetry of the resulting bulge of mass is permanently shifted out of the direction of the Sun by the Earth’s rotation. As for the body tides, this loss of symmetry is responsible for the torque which, in the present conditions, tends to accelerate the spin (Correia & Laskar 2003).
最后，每日太阳加热会产生扭矩，导致气压重新分布，主要由半日波驱动，因此形成了所谓的热大气潮汐（Chapman & Lindzen 1970）。由于地球自转，由此产生的质量凸起的对称轴永久地偏离了太阳方向。与体潮汐一样，这种对称性的丧失是扭矩产生的根源，在目前的条件下，扭矩往往会加速自转（Correia & Laskar 2003）。
Volland (1978) showed that this effect is not negligible, reducing the Earth’s despinning by about 7.5%. But, although the estimate of their long term contributions deserves careful attention, we have not taken the atmospheric tides into account in the computations presented in the next sections, assuming that, as the uncertainty on some of the other factors is still large, the global results obtained here will not differ much when taking this additional effect into consideration.
Volland (1978) 指出，这种效应不可忽略，它使地球自转减慢了约 7.5%。尽管对其长期贡献的估算值得认真考虑，但我们在下一节的计算中并未考虑大气潮汐，因为其他一些因素的不确定性仍然很大，即使考虑到这一额外效应，本文得出的全球结果也不会有太大差异。

4.2.3. Mantle convection  4.2.3 地幔对流

Several internal geophysical processes are able to affect the moments of inertia of the Earth and thereby its precession constant value. Redistribution of mass within the Earth can occur as a consequence of plate subduction, upwelling plumes, mantle avalanches or density anomalies driven by mantle convection. Typical time-scales of these phenomena range from ∼ 10 ∼ 10 ∼10\sim 10 to 100 Myr . However, their impact is still poorly known. Forte & Mitrovica (1997) computed convection-induced perturbations in the dynamical ellipticity of the Earth over the last 20 Myr , using different sets of 3D seismic heterogeneity models and mantle viscosity profiles. Most of the models indicate that mantle convection led to a mean relative increase of the precession constant close to ∼ 0.5 % ∼ 0.5 % ∼0.5%\sim 0.5 \% since the beginning of the Neogene period ( ∼ 23 Myr ∼ 23 Myr ∼23Myr\sim 23 \mathrm{Myr} ago). Although this value is similar to our estimate of the tidal dissipation effect (but in the opposite direction, see Sect. 9), we did not take account of its impact in our computations. Because of its still large uncertainty, we feel that it is more appropriate to include its possible contribution as a perturbation of the tidal dissipation parameter.
有多种地球物理过程能够影响地球的转动惯量，从而影响其岁差常数值。板块俯冲、上涌柱、地幔雪崩或地幔对流引起的密度异常都可能导致地球内部质量的重新分布。这些现象的典型时间尺度从 ∼ 10 ∼ 10 ∼10\sim 10 到 1 亿年。然而，人们对其影响仍然知之甚少。Forte 和 Mitrovica (1997) 使用不同的 3D 地震非均匀性模型和地幔粘度剖面，计算了过去 20 亿年中对流引起的地球动态椭圆度扰动。大多数模型表明，自新近纪开始（ ∼ 23 Myr ∼ 23 Myr ∼23Myr\sim 23 \mathrm{Myr} 以前），地幔对流导致岁差常数平均相对增加接近 ∼ 0.5 % ∼ 0.5 % ∼0.5%\sim 0.5 \% 。虽然该值与我们对潮汐耗散效应的估计值相似（但方向相反，参见第 9 节），但我们在计算中并未考虑其影响。由于其不确定性仍然很大，我们认为将其可能的贡献作为潮汐耗散参数的扰动更为合适。

4.2.4. Climate friction  4.2.4. 气候摩擦

Climate friction is a positive and dissipative feedback between obliquity variations and climate which may cause a secular drift of the spin axis. In response to quasi-periodic variations in the obliquity, glacial and interglacial conditions drive transport of water into and out of the polar regions, affecting the dynamical ellipticity of the Earth. A significant fraction of the surface loading is compensated by viscous flow within the Earth, but delayed responses in both climatic and viscous relaxation processes may introduce a secular term in the obliquity evolution (Rubincam 1990, 1995). Because both the ice load history and the visco-elastic structure of the Earth are not strongly constrained, it is difficult to produce accurate simulations and predictions of the coupled response of the entire system. Nevertheless, simplifying assumptions have been used to estimate the magnitude and the direction of the secular drift. Several analyses have examined this phenomenon, suggesting that climate friction could have changed the Earth’s obliquity by more that 20 ∘ 20 ∘ 20^(@)20^{\circ} over its whole geological history (Bills 1994; Ito et al. 1995; Williams et al. 1998). A more detailed theoretical and numerical treatment of climate friction has been developped by Levrard & Laskar (2003). Using available constraints on the ice volume response to obliquity forcing based on δ 18 O δ 18 O delta^(18)O\delta^{18} \mathrm{O} oxygen-isotope records, they showed that climate friction impact is likely negligible. Over the last 3 Ma ,
气候摩擦是地轴倾角变化与气候之间的正耗散反馈，可能导致自转轴的长期漂移。由于地轴倾角的准周期性变化，冰川期和间冰期条件会驱动水流进出两极地区，从而影响地球的动力学椭圆度。地球表面负荷的很大一部分由地球内部的粘性流动补偿，但气候和粘性松弛过程的延迟响应可能会在地轴倾角的演化中引入一个长期项（Rubincam，1990，1995）。由于冰荷载历史和地球的粘弹性结构都未受到严格约束，因此很难对整个系统的耦合响应进行精确的模拟和预测。尽管如此，人们还是会使用简化假设来估计长期漂移的幅度和方向。多项分析已对这一现象进行了检验，表明气候摩擦可能在整个地质历史时期内使地球的自转轴倾角改变超过 20 ∘ 20 ∘ 20^(@)20^{\circ} （Bills 1994；Ito 等人 1995；Williams 等人 1998）。Levrard 和 Laskar（2003）对气候摩擦进行了更详细的理论和数值处理。他们利用基于 δ 18 O δ 18 O delta^(18)O\delta^{18} \mathrm{O} 氧同位素记录的冰体积对自转轴倾角强迫响应的现有约束，表明气候摩擦的影响可能可以忽略不计。在过去的 3 百万年里，
Fig. 4. Differences La2003-DE406 over the full range of DE406 ( -5000 yr to +1000 from J2000). The units for semi-major axis (a) are AU , arcsec for mean longitude ( l l ll ) and longitude of perihelion (perihelion). The eccentricity ( e ) ( e ) (e)(e) and the inclinations variables ( q = q = q=q= sin ⁡ ( i / 2 ) cos ⁡ ( Ω ) , p = sin ⁡ ( i / 2 ) sin ⁡ ( Ω ) sin ⁡ ( i / 2 ) cos ⁡ ( Ω ) , p = sin ⁡ ( i / 2 ) sin ⁡ ( Ω ) sin(i//2)cos(Omega),p=sin(i//2)sin(Omega)\sin (i / 2) \cos (\Omega), p=\sin (i / 2) \sin (\Omega) ), where i i ii and Ω Ω Omega\Omega are the inclination and node from the ecliptic and equinox J2000.
图 4. La2003 至 DE406 在整个 DE406 范围内的差异（从 J2000 开始，-5000 年至 +1000 年）。半长轴 (a) 的单位为天文单位 (AU)，平均经度 ( l l ll ) 和近日点经度 ( 近日点 ) 的单位为弧秒。偏心率 ( e ) ( e ) (e)(e) 和倾角变量 ( q = q = q=q= sin ⁡ ( i / 2 ) cos ⁡ ( Ω ) , p = sin ⁡ ( i / 2 ) sin ⁡ ( Ω ) sin ⁡ ( i / 2 ) cos ⁡ ( Ω ) , p = sin ⁡ ( i / 2 ) sin ⁡ ( Ω ) sin(i//2)cos(Omega),p=sin(i//2)sin(Omega)\sin (i / 2) \cos (\Omega), p=\sin (i / 2) \sin (\Omega) )，其中 i i ii 和 Ω Ω Omega\Omega 分别是相对于黄道和春分点 J2000 的倾角和交点。
corresponding to the intensification of the Northern Hemisphere glaciation, a maximal absolute drift of only ∼ 0.04 ∘ / Myr ∼ 0.04 ∘ / Myr ∼0.04^(@)//Myr\sim 0.04^{\circ} / \mathrm{Myr} has been estimated for a realistic lower mantle viscosity of 10 22 Pa 10 22 Pa 10^(22)Pa10^{22} \mathrm{~Pa}. s. Before 3 Ma , climate friction impact is expected to be much more negligible because of the lack of massive ice caps. In this context, this effect was not taken into account in our long-term computations. Indeed, the uncertainty on the tidal dissipation is still important and it is expected that most of the possible effect of climate friction would be absorbed by a small change of the main tidal dissipation parameters.
与北半球冰川作用的加剧相对应，对于实际的下地幔粘度 10 22 Pa 10 22 Pa 10^(22)Pa10^{22} \mathrm{~Pa} . s-1，最大绝对漂移估计仅为 ∼ 0.04 ∘ / Myr ∼ 0.04 ∘ / Myr ∼0.04^(@)//Myr\sim 0.04^{\circ} / \mathrm{Myr} 。在 300 万年前，由于缺乏大规模冰盖，预计气候摩擦的影响可以忽略不计。因此，我们在长期计算中并未考虑这种影响。事实上，潮汐耗散的不确定性仍然很重要，预计气候摩擦的大部分可能影响将被主要潮汐耗散参数的微小变化所吸收。
Nevertheless, glaciation-induced perturbations in the dynamical ellipticity of the Earth change the precession frequency and thereby the obliquity frequencies leading to a temporal shift between a perturbated and a nominal solution (see Sect. 9). Estimates of this offset have been performed by Mitrovica & Forte (1995) and Levrard & Laskar (2003) for a large set of mantle viscosity profiles. It may reach between ∼ 0.7 ∼ 0.7 ∼0.7\sim 0.7 and 8 kyr after 3 Myr , corresponding respectively to extreme lower mantle viscosities of 3 × 10 21 3 × 10 21 3xx10^(21)3 \times 10^{21} and 10 23 Pa 10 23 Pa 10^(23)Pa10^{23} \mathrm{~Pa}. s. In the last case, the Earth is nearly rigid on orbital timescales and the corresponding offset can be thus reasonably considered as a maximal value.
尽管如此，冰川作用引起的地球动态椭圆率扰动会改变岁差频率，进而改变地轴倾角频率，导致扰动解和标准解之间的时间偏移（见第 9 节）。Mitrovica & Forte（1995）和 Levrard & Laskar（2003）利用大量地幔粘度剖面对这一偏移进行了估算。它可能在 3 百万年后达到 ∼ 0.7 ∼ 0.7 ∼0.7\sim 0.7 到 8 千年之间，分别对应于极端的下地幔粘度 3 × 10 21 3 × 10 21 3xx10^(21)3 \times 10^{21} 和 10 23 Pa 10 23 Pa 10^(23)Pa10^{23} \mathrm{~Pa} .s。在最后一种情况下，地球在轨道时间尺度上几乎是刚性的，因此相应的偏移量可以合理地视为最大值。

5. Comparison with DE406  5. 与 DE406 的比较

Using a direct numerical integrator, our goal is to provide a long term solution for the orbital and precessional elements of the Earth with a precision that is comparable with the usual accuracy of a short time ephemeris. We have thus compared our solution with the most advanced present numerical integration, DE406, that was itself adjusted to the observations (Standish 1998). Over the full range of DE406, that is from -5000 yr to +1000 yr from the present date, the maximum difference in the position of the Earth-Moon barycenter is less than 0.09 arcsec in longitude, and the difference in eccentricity less than 10 − 8 10 − 8 10^(-8)10^{-8}. The difference in the longitude of the Moon are more important, as the dissipative models that are used in the two integrations are slightly different. They amount to 240 arcsec after 5000 years, and the eccentricity difference is of the order of 2 × 10 − 5 2 × 10 − 5 2xx10^(-5)2 \times 10^{-5}, but it should be noted that only the averaged motion of the Moon will have some significant influence on the precession and obliquity of the Earth. Table 2 summarizes the maximum differences between La2004 and DE406 for the main orbital elements, over the whole interval of DE406, for different time intervals.
我们的目标是使用直接数值积分器，为地球的轨道和岁差元素提供长期解，其精度可与短时星历表的通常精度相媲美。因此，我们将我们的解与目前最先进的数值积分 DE406 进行了比较，后者本身已根据观测结果进行了调整（Standish 1998）。在 DE406 的整个范围内，即从现在起 -5000 年至 +1000 年，地月质心位置的最大经度差异小于 0.09 角秒，偏心率差异小于 10 − 8 10 − 8 10^(-8)10^{-8} 。月球经度的差异更为重要，因为这两个积分中使用的耗散模型略有不同。 5000 年后，它们的差异达到 240 角秒，偏心率差异约为 2 × 10 − 5 2 × 10 − 5 2xx10^(-5)2 \times 10^{-5} ，但需要注意的是，只有月球的平均运动才会对地球的岁差和倾角产生显著的影响。表 2 总结了 La2004 和 DE406 在 DE406 整个周期内不同时间间隔的主要轨道根数的最大差异。

6. Comparison with La93  6. 与 La93 的比较

We have compared in Figs. 5 and 6 the eccentricity and inclination of the secular solution La93 and the new present solution La2004. Over 10 Myr , the two orbital solutions La93 and La2004 are nearly identical, and do not differ significantly over 20 Myr . The differences in eccentricity increase regularly with time, amounting to about 0.02 in the eccentricity after 20 Myr , about 1 / 3 1 / 3 1//31 / 3 of the total amplitude ≈ 0.063 ≈ 0.063 ~~0.063\approx 0.063, while the difference in the orbital inclination reaches 1 degree, to compare to a maximum variation of 4.3 degrees. These differences result mostly from a small difference in the main secular frequency g 6 g 6 g_(6)g_{6} from Jupiter and Saturn that is now g 6 = 28.2450 arcsec / year g 6 = 28.2450 arcsec / year g_(6)=28.2450arcsec//yearg_{6}=28.2450 \mathrm{arcsec} / \mathrm{year} (Table 3) instead of 28.2207 arcsec / year 28.2207 arcsec / year 28.2207arcsec//year28.2207 \mathrm{arcsec} / \mathrm{year} in the previous solution. Beyond 20 million years, the differences between the two solutions become more noticeable. The solutions in obliquity and climatic precession are plotted in Figs. 7 and 8. The difference in obliquity over 20 Myr amounts to about 2 degrees, which means that the obliquity cycles of the two solutions become nearly out of phase after this date. The differences La2004-La93 are plotted in (a), while in (b) only the precession model has been changed and in © only the orbital model has been changed. Comparing 7 b and 7 c shows that most of the differences between the La93 and La2004 obliquity solutions result from the change in the dissipative model of the Earth-Moon system, while the only change of the orbital solution would lead to a much smaller change of only 0.6 degree in obliquity after 20 Myr . Similar conclusions can be made for the climatic precession (Figs. 8a-c).
在图 5 和图 6 中，我们比较了长期解 La93 和新解 La2004 的偏心率和轨道倾角。在 10 Myr 内，两个轨道解 La93 和 La2004 几乎相同，在 20 Myr 内没有显著差异。偏心率的差异随着时间的推移而有规律地增加，在 20 Myr 之后偏心率约为 0.02，约为总振幅 ≈ 0.063 ≈ 0.063 ~~0.063\approx 0.063 的 1 / 3 1 / 3 1//31 / 3 ，而轨道倾角的差异达到 1 度，而最大变化为 4.3 度。这些差异主要是由于木星和土星的主要长期频率 g 6 g 6 g_(6)g_{6} 的微小差异造成的，现在为 g 6 = 28.2450 arcsec / year g 6 = 28.2450 arcsec / year g_(6)=28.2450arcsec//yearg_{6}=28.2450 \mathrm{arcsec} / \mathrm{year} （表 3），而不是以前解中的 28.2207 arcsec / year 28.2207 arcsec / year 28.2207arcsec//year28.2207 \mathrm{arcsec} / \mathrm{year} 。在 2000 万年之后，两个解之间的差异变得更加明显。倾角和气候岁差的解绘制在图中。 7 和 8。20 Myr 后倾角的差异约为 2 度，这意味着在此日期之后，两个解的倾角周期几乎不同步。 (a) 绘制了 La2004-La93 的差异，而 (b) 中仅改变了岁差模型，(c) 中仅改变了轨道模型。比较 7 b 和 7 c 表明，La93 和 La2004 倾角解之间的大部分差异来自地月系统耗散模型的变化，而仅仅改变轨道解会导致 20 Myr 后倾角的变化小得多，仅为 0.6 度。对于气候岁差也可以得出类似的结论（图 8a-c）。

7. Secular frequencies  7. 世俗频率

Over long time scales, the determination of longitude of the planet on its orbit is not very important, and we are more concerned by the slow evolution of the Earth orbit under secular planetary perturbations. The semi-major axis of the Earth
在长期尺度上，行星在其轨道上的经度的确定并不十分重要，我们更关心的是地球轨道在长期行星摄动下的缓慢演变。地球的半长轴
Table 2. Maximum difference between La2004 and DE406 over the whole time interval of DE406 ( -5000 yr to +1000 yr with origin at J2000); Col. 1: -100 to +100 yr ; Col. 2 : − 1000 2 : − 1000 2:-10002:-1000 to +1000 yr ; Col. 3: -5000 to +1000 yr. EMB is the Earth-Moon barycenter.
表 2. DE406 整个时间间隔内 La2004 与 DE406 之间的最大差异（-5000 年至 +1000 年，起源于 J2000）；第 1 列：-100 至 +100 年；第 2 : − 1000 2 : − 1000 2:-10002:-1000 列至 +1000 年；第 3 列：-5000 至 +1000 年。EMB 是地月质心。
λ ( ′ ′ × 1000 ) λ ′ ′ × 1000 lambda(^('')xx1000)\lambda\left({ }^{\prime \prime} \times 1000\right)
Mercury  汞 8 60 90
Venus  金星 12 107 143
EMB 3 28 79
Mars  行进 19 171 277
Jupiter  木星 5 69 69
Saturn  土星 1 4 7
Uranus  天王星 1 3 4
Neptune  海王星 2 7 16
Pluto  冥王星 3 9 17
Moon  月亮 13470 85873 221724
Earth  地球 10 88 238
a ( UA × 10 10 ) a UA × 10 10 a(UAxx10^(10))a\left(\mathrm{UA} \times 10^{10}\right)
Mercury  汞 4 43 58
Venus  金星 13 23 27
EMB 29 62 62
Mars  行进 27 76 76
Jupiter  木星 268 470 565
Saturn  土星 743 1073 1168
Uranus  天王星 1608 3379 3552
Neptune  海王星 3315 5786 6458
Pluto  冥王星 4251 10603 10603
Moon  月亮 23 204 520
Earth  地球 487 3918 10238
e ( × 10 10 ) e × 10 10 e(xx10^(10))e\left(\times 10^{10}\right)
Mercury  汞 33 52 144
Venus  金星 12 33 103
EMB 30 80 98
Mars  行进 62 130 468
Jupiter  木星 30 95 158
Saturn  土星 69 107 138
Uranus  天王星 99 135 144
Neptune  海王星 65 177 201
Pluto  冥王星 92 215 215
Moon  月亮 11152 252832 252832
Earth  地球 460 3669 8694
sin ⁡ ( i / 2 ) ( × 10 10 ) sin ⁡ ( i / 2 ) × 10 10 sin(i//2)(xx10^(10))\sin (i / 2)\left(\times 10^{10}\right)
Mercury  汞 2 18 88
Venus  金星 5 55 144
EMB 17 190 414
Mars  行进 48 492 2053
Jupiter  木星 1 11 122
Saturn  土星 1 6 100
Uranus  天王星 2 5 5
Neptune  海王星 2 6 6
Pluto  冥王星 5 13 17
Moon  月亮 5765 97275 600789
Earth  地球 19 202 841
lambda(^('')xx1000) Mercury 8 60 90 Venus 12 107 143 EMB 3 28 79 Mars 19 171 277 Jupiter 5 69 69 Saturn 1 4 7 Uranus 1 3 4 Neptune 2 7 16 Pluto 3 9 17 Moon 13470 85873 221724 Earth 10 88 238 a(UAxx10^(10)) Mercury 4 43 58 Venus 13 23 27 EMB 29 62 62 Mars 27 76 76 Jupiter 268 470 565 Saturn 743 1073 1168 Uranus 1608 3379 3552 Neptune 3315 5786 6458 Pluto 4251 10603 10603 Moon 23 204 520 Earth 487 3918 10238 e(xx10^(10)) Mercury 33 52 144 Venus 12 33 103 EMB 30 80 98 Mars 62 130 468 Jupiter 30 95 158 Saturn 69 107 138 Uranus 99 135 144 Neptune 65 177 201 Pluto 92 215 215 Moon 11152 252832 252832 Earth 460 3669 8694 sin(i//2)(xx10^(10)) Mercury 2 18 88 Venus 5 55 144 EMB 17 190 414 Mars 48 492 2053 Jupiter 1 11 122 Saturn 1 6 100 Uranus 2 5 5 Neptune 2 6 6 Pluto 5 13 17 Moon 5765 97275 600789 Earth 19 202 841| | $\lambda\left({ }^{\prime \prime} \times 1000\right)$ | | | | :---: | :---: | :---: | :---: | | Mercury | 8 | 60 | 90 | | Venus | 12 | 107 | 143 | | EMB | 3 | 28 | 79 | | Mars | 19 | 171 | 277 | | Jupiter | 5 | 69 | 69 | | Saturn | 1 | 4 | 7 | | Uranus | 1 | 3 | 4 | | Neptune | 2 | 7 | 16 | | Pluto | 3 | 9 | 17 | | Moon | 13470 | 85873 | 221724 | | Earth | 10 | 88 | 238 | | | $a\left(\mathrm{UA} \times 10^{10}\right)$ | | | | Mercury | 4 | 43 | 58 | | Venus | 13 | 23 | 27 | | EMB | 29 | 62 | 62 | | Mars | 27 | 76 | 76 | | Jupiter | 268 | 470 | 565 | | Saturn | 743 | 1073 | 1168 | | Uranus | 1608 | 3379 | 3552 | | Neptune | 3315 | 5786 | 6458 | | Pluto | 4251 | 10603 | 10603 | | Moon | 23 | 204 | 520 | | Earth | 487 | 3918 | 10238 | | | $e\left(\times 10^{10}\right)$ | | | | Mercury | 33 | 52 | 144 | | Venus | 12 | 33 | 103 | | EMB | 30 | 80 | 98 | | Mars | 62 | 130 | 468 | | Jupiter | 30 | 95 | 158 | | Saturn | 69 | 107 | 138 | | Uranus | 99 | 135 | 144 | | Neptune | 65 | 177 | 201 | | Pluto | 92 | 215 | 215 | | Moon | 11152 | 252832 | 252832 | | Earth | 460 | 3669 | 8694 | | | $\sin (i / 2)\left(\times 10^{10}\right)$ | | | | Mercury | 2 | 18 | 88 | | Venus | 5 | 55 | 144 | | EMB | 17 | 190 | 414 | | Mars | 48 | 492 | 2053 | | Jupiter | 1 | 11 | 122 | | Saturn | 1 | 6 | 100 | | Uranus | 2 | 5 | 5 | | Neptune | 2 | 6 | 6 | | Pluto | 5 | 13 | 17 | | Moon | 5765 | 97275 | 600789 | | Earth | 19 | 202 | 841 |
Fig. 5. Eccentricity of the Earth over 25 Myr in negative time from J2000. The solid line stands for the present solution La2004, while the dotted line is the eccentricity in the La93 solution (Laskar et al. 1993). The differences of the two solutions becomes noticeable after 10 Myr , and significantly different after 15-20 Myr.
图 5. 自 2000 年 J 以来，地球偏心率在 25 百万年内的变化（负值）。实线代表当前 La2004 解，虚线代表 La93 解（Laskar et al. 1993）中的偏心率。两种解的差异在 10 百万年后变得明显，在 15-20 百万年后差异更为显著。
Fig. 6. Inclination (in degrees) of the Earth with respect to the fixed ecliptic J2000 over 25 Myr in negative time from J2000. The solid line stands for the present solution La2004, while the dotted line is the eccentricity in the La93 solution (Laskar et al. 1993). The differences of the two solutions become noticeable after 15 Myr , and significantly different after 20 Myr .
图 6. 地球相对于固定黄道 J2000 的倾角（以度为单位），自 J2000 算起，在 25 百万年内呈负值。实线代表当前 La2004 解，虚线代表 La93 解（Laskar 等人，1993）中的偏心率。两种解的差异在 15 百万年后变得明显，在 20 百万年后则显著不同。
has only very small variations (Fig. 11). Indeed, as shown by Laplace (1773) and Lagrange (1776), there are no secular variations of the semi-major axis of the planets at first order with respect to the masses, while some terms of higher order can be present (Haretu 1885). These terms are very small for the
只有非常小的变化（图 11）。事实上，正如拉普拉斯（1773）和拉格朗日（1776）所指出的，行星半长轴相对于质量不存在一阶长期变化，而可以存在一些高阶项（Haretu 1885）。这些项对于
inner planets, but more visible in the solutions of the outer planets where the proximity of some mean motion resonances increases their contribution at second order with respect to the masses (Milani et al. 1987; Bretagnon & Simon 1990). For the solution of the Earth, they are so small that they will not induce
内行星的解中，这一点更为明显，但在外行星的解中更为明显，因为一些平均运动共振的接近增加了它们相对于质量的二阶贡献 (Milani et al. 1987; Bretagnon & Simon 1990)。对于地球的解来说，它们太小了，以至于不会诱导
Table 3. Main secular frequencies g i g i g_(i)g_{i} and s i s i s_(i)s_{i} of La2004 determined over 20 Ma for the four inner planets, and over 50 Ma for the 5 outer planets (in arcsec yr − 1 yr − 1 yr^(-1)\mathrm{yr}^{-1} ). Δ 100 Δ 100 Delta_(100)\Delta_{100} and Δ 250 Δ 250 Delta_(250)\Delta_{250} are the observed variations of the frequencies over respectively 100 and 250 Myr . In the last column, the period of the secular term are given.
表 3. La2004 的主要长期频率 g i g i g_(i)g_{i} 和 s i s i s_(i)s_{i} ，针对四颗内行星，在 20 Ma 内确定，针对五颗外行星，在 50 Ma 内确定（以角秒 yr − 1 yr − 1 yr^(-1)\mathrm{yr}^{-1} 为单位）。 Δ 100 Δ 100 Delta_(100)\Delta_{100} 和 Δ 250 Δ 250 Delta_(250)\Delta_{250} 分别是 1 亿年和 2.5 亿年内观测到的频率变化。最后一列给出了长期周期。
( ′ ′ / yr ) ′ ′ / yr (^('')//yr)\left({ }^{\prime \prime} / \mathrm{yr}\right) Δ 100 Δ 100 Delta_(100)\Delta_{100} Δ 250 Δ 250 Delta_(250)\Delta_{250} Period (yr)  期间（年）
g 1 g 1 g_(1)g_{1} 5.59 0.13 0.20 231843
g 2 g 2 g_(2)g_{2} 7.452 0.019 0.023 173913
g 3 g 3 g_(3)g_{3} 17.368 0.20 0.20 74620
g 4 g 4 g_(4)g_{4} 17.916 0.20 0.20 72338
g 5 g 5 g_(5)g_{5} 4.257452 0.000030 0.000030 304407
g 6 g 6 g_(6)g_{6} 28.2450 0.0010 0.0010 45884
g 7 g 7 g_(7)g_{7} 3.087951 0.000034 0.000048 419696
g 8 g 8 g_(8)g_{8} 0.673021 0.000015 0.000021 1925646
g 9 g 9 g_(9)g_{9} -0.34994 0.00063 0.0012 3703492
s 1 s 1 s_(1)s_{1} -5.59 0.15 0.16 231843
s 2 s 2 s_(2)s_{2} -7.05 0.19 0.25 183830
s 3 s 3 s_(3)s_{3} -18.850 0.066 0.11 68753
s 4 s 4 s_(4)s_{4} -17.755 0.064 0.14 72994
s 5 s 5 s_(5)s_{5} 0.00000013 0.0000001 0.0000001
s 6 s 6 s_(6)s_{6} -26.347855 0.000076 0.000087 49188
s 7 s 7 s_(7)s_{7} -2.9925259 0.000025 0.000025 433079
s 8 s 8 s_(8)s_{8} -0.691736 0.000010 0.000011 1873547
s 9 s 9 s_(9)s_{9} -0.34998 0.00051 0.0011 3703069
(^('')//yr) Delta_(100) Delta_(250) Period (yr) g_(1) 5.59 0.13 0.20 231843 g_(2) 7.452 0.019 0.023 173913 g_(3) 17.368 0.20 0.20 74620 g_(4) 17.916 0.20 0.20 72338 g_(5) 4.257452 0.000030 0.000030 304407 g_(6) 28.2450 0.0010 0.0010 45884 g_(7) 3.087951 0.000034 0.000048 419696 g_(8) 0.673021 0.000015 0.000021 1925646 g_(9) -0.34994 0.00063 0.0012 3703492 s_(1) -5.59 0.15 0.16 231843 s_(2) -7.05 0.19 0.25 183830 s_(3) -18.850 0.066 0.11 68753 s_(4) -17.755 0.064 0.14 72994 s_(5) 0.00000013 0.0000001 0.0000001 s_(6) -26.347855 0.000076 0.000087 49188 s_(7) -2.9925259 0.000025 0.000025 433079 s_(8) -0.691736 0.000010 0.000011 1873547 s_(9) -0.34998 0.00051 0.0011 3703069| | $\left({ }^{\prime \prime} / \mathrm{yr}\right)$ | $\Delta_{100}$ | $\Delta_{250}$ | Period (yr) | | :--- | :--- | :--- | :--- | ---: | | $g_{1}$ | 5.59 | 0.13 | 0.20 | 231843 | | $g_{2}$ | 7.452 | 0.019 | 0.023 | 173913 | | $g_{3}$ | 17.368 | 0.20 | 0.20 | 74620 | | $g_{4}$ | 17.916 | 0.20 | 0.20 | 72338 | | $g_{5}$ | 4.257452 | 0.000030 | 0.000030 | 304407 | | $g_{6}$ | 28.2450 | 0.0010 | 0.0010 | 45884 | | $g_{7}$ | 3.087951 | 0.000034 | 0.000048 | 419696 | | $g_{8}$ | 0.673021 | 0.000015 | 0.000021 | 1925646 | | $g_{9}$ | -0.34994 | 0.00063 | 0.0012 | 3703492 | | $s_{1}$ | -5.59 | 0.15 | 0.16 | 231843 | | $s_{2}$ | -7.05 | 0.19 | 0.25 | 183830 | | $s_{3}$ | -18.850 | 0.066 | 0.11 | 68753 | | $s_{4}$ | -17.755 | 0.064 | 0.14 | 72994 | | $s_{5}$ | 0.00000013 | 0.0000001 | 0.0000001 | | | $s_{6}$ | -26.347855 | 0.000076 | 0.000087 | 49188 | | $s_{7}$ | -2.9925259 | 0.000025 | 0.000025 | 433079 | | $s_{8}$ | -0.691736 | 0.000010 | 0.000011 | 1873547 | | $s_{9}$ | -0.34998 | 0.00051 | 0.0011 | 3703069 |
Fig. 7. Differences in obliquity (in degrees) La2004-La93(1, 1) versus time (in Myr) over 20 Myr from present a). In La2004 both orbital and precession models have been changed. In b b b\mathbf{b} ) only the precession model has been changed, while in c) only the orbital motion has been changed. It is thus clear that the main change in La2004 arise from the change of precession model.
图 7. 自现在起 20 百万年以来，La2004-La93(1, 1) 的倾角（度）随时间（百万年）的变化 a)。在 La2004 中，轨道和进动模型都发生了变化。在 b b b\mathbf{b} ) 中，只有进动模型发生了变化；而在 c) 中，只有轨道运动发生了变化。由此可见，La2004 的主要变化源于进动模型的变化。
any noticeable change of the mean Earth-Sun distance, or of the duration of the year in the geological past, at least over the past 250 millions of years.
地球与太阳的平均距离或地质历史上一年的长度的任何明显变化，至少在过去的 2.5 亿年里。
It was first shown by Lagrange ( 1774 , 1777 ) ( 1774 , 1777 ) (1774,1777)(1774,1777) that the inclination and nodes of the planets suffer long term quasiperiodic variations. Shortly after, Laplace (1776) demonstrated that this was the same for the eccentricity and longitudes. Both computations were made using the linearization of the first order
拉格朗日 ( 1774 , 1777 ) ( 1774 , 1777 ) (1774,1777)(1774,1777) 首次证明了行星的倾角和节点会经历长期的准周期变化。不久之后，拉普拉斯（1776）证明了偏心率和经度也存在同样的变化。这两次计算均采用一阶线性化方法
Fig. 8. Differences in climatic precession ( e sin ⁡ ω e sin ⁡ ω e sin omegae \sin \omega, where ω = ϖ + ψ ω = ϖ + ψ omega=ϖ+psi\omega=\varpi+\psi is the perihelion from the moving equinox) La2004-La93(1, 1) versus time (in Myr) over 20 Myr from present a). In La2004 both orbital and precession models have been changed. In b) only the precession model has been changed, while in c) only the orbital motion has been changed. It is thus clear that the main change in La2004 arise from the change of precession model.
图 8. 气候岁差（ e sin ⁡ ω e sin ⁡ ω e sin omegae \sin \omega ，其中 ω = ϖ + ψ ω = ϖ + ψ omega=ϖ+psi\omega=\varpi+\psi 为距移动春分点的近日点）La2004-La93(1, 1) 随时间（以百万年为单位）的变化，自现在起已超过 20 百万年。a) 在 La2004 中，轨道和岁差模型都发生了变化。b) 中仅岁差模型发生了变化，而 c) 中仅轨道运动发生了变化。由此可见，La2004 的主要变化源于岁差模型的变化。
averaged equations of motion. In this approximation, if we use the complex notations
平均运动方程。在这个近似中，如果我们使用复数符号
z k = e k exp ⁡ ( i ϖ k ) ; ζ k = sin ⁡ ( i k / 2 ) exp ⁡ ( i Ω k ) z k = e k exp ⁡ i ϖ k ; ζ k = sin ⁡ i k / 2 exp ⁡ i Ω k z_(k)=e_(k)exp(iϖ_(k));quadzeta_(k)=sin(i_(k)//2)exp(iOmega_(k))z_{k}=e_{k} \exp \left(i \varpi_{k}\right) ; \quad \zeta_{k}=\sin \left(i_{k} / 2\right) \exp \left(i \Omega_{k}\right),
the equations of motion are given by a linear equation
运动方程由线性方程给出
d [ x ] d t = i A [ x ] d [ x ] d t = i A [ x ] (d[x])/(dt)=iA[x]\frac{\mathrm{d}[x]}{\mathrm{d} t}=i A[x]
where [ x ] [ x ] [x][x] is the column vector ( z 1 , … , z 9 , ζ 1 , … , ζ 9 ) z 1 , … , z 9 , ζ 1 , … , ζ 9 (z_(1),dots,z_(9),zeta_(1),dots,zeta_(9))\left(z_{1}, \ldots, z_{9}, \zeta_{1}, \ldots, \zeta_{9}\right), and
其中 [ x ] [ x ] [x][x] 是列向量 ( z 1 , … , z 9 , ζ 1 , … , ζ 9 ) z 1 , … , z 9 , ζ 1 , … , ζ 9 (z_(1),dots,z_(9),zeta_(1),dots,zeta_(9))\left(z_{1}, \ldots, z_{9}, \zeta_{1}, \ldots, \zeta_{9}\right) ，并且
A = ( A 1 0 0 A 2 ) A = A 1 0 0 A 2 A=([A_(1),0],[0,A_(2)])A=\left(\begin{array}{cc}A_{1} & 0 \\ 0 & A_{2}\end{array}\right)
where the 9 × 9 9 × 9 9xx99 \times 9 matrices (for 9 planets) A 1 , A 2 A 1 , A 2 A_(1),A_(2)A_{1}, A_{2} depend on the planetary masses and semi-major axis of the planets. The resolution of these equations is now classical, and is obtained through the diagonalization of the matrix A A AA by the change of coordinates into proper modes u u uu
其中， 9 × 9 9 × 9 9xx99 \times 9 矩阵（针对 9 颗行星） A 1 , A 2 A 1 , A 2 A_(1),A_(2)A_{1}, A_{2} 取决于行星质量和行星的半长轴。这些方程的解析式现已成为经典方程，可通过将矩阵 A A AA 对角化（将坐标变换为本征模态 u u uu ）获得。
[ x ] = S [ u ] [ x ] = S [ u ] [x]=S[u][x]=S[u].
The solution in the proper modes [ u ] = ( z 1 ∙ , … , z 9 ∙ , ζ 1 ∙ , … , ζ 9 ∙ ) [ u ] = z 1 ∙ , … , z 9 ∙ , ζ 1 ∙ , … , ζ 9 ∙ [u]=(z_(1)^(∙),dots,z_(9)^(∙),zeta_(1)^(∙),dots,zeta_(9)^(∙))[u]=\left(z_{1}^{\bullet}, \ldots, z_{9}^{\bullet}, \zeta_{1}^{\bullet}, \ldots, \zeta_{9}^{\bullet}\right) is then given by the diagonal system
然后，通过对角系统给出适当模式 [ u ] = ( z 1 ∙ , … , z 9 ∙ , ζ 1 ∙ , … , ζ 9 ∙ ) [ u ] = z 1 ∙ , … , z 9 ∙ , ζ 1 ∙ , … , ζ 9 ∙ [u]=(z_(1)^(∙),dots,z_(9)^(∙),zeta_(1)^(∙),dots,zeta_(9)^(∙))[u]=\left(z_{1}^{\bullet}, \ldots, z_{9}^{\bullet}, \zeta_{1}^{\bullet}, \ldots, \zeta_{9}^{\bullet}\right) 的解
d [ u ] d t = i D [ u ] d [ u ] d t = i D [ u ] (d[u])/(dt)=iD[u]\frac{\mathrm{d}[u]}{\mathrm{d} t}=i D[u]
where D = S − 1 A S D = S − 1 A S D=S^(-1)ASD=S^{-1} A S is a diagonal matrix D = D = D=D= Diag ( g 1 , g 2 , … , g 9 , s 1 , s 2 , … , s 9 ) Diag g 1 , g 2 , … , g 9 , s 1 , s 2 , … , s 9 Diag(g_(1),g_(2),dots,g_(9),s_(1),s_(2),dots,s_(9))\operatorname{Diag}\left(g_{1}, g_{2}, \ldots, g_{9}, s_{1}, s_{2}, \ldots, s_{9}\right), where all the eigenvalues g k , s k g k , s k g_(k),s_(k)g_{k}, s_{k} are real. The solutions of the proper modes z k ∗ , ζ k ∙ z k ∗ , ζ k ∙ z_(k)^(**),zeta_(k)^(∙)z_{k}^{*}, \zeta_{k}^{\bullet} are
其中 D = S − 1 A S D = S − 1 A S D=S^(-1)ASD=S^{-1} A S 是对角矩阵 D = D = D=D= Diag ( g 1 , g 2 , … , g 9 , s 1 , s 2 , … , s 9 ) Diag g 1 , g 2 , … , g 9 , s 1 , s 2 , … , s 9 Diag(g_(1),g_(2),dots,g_(9),s_(1),s_(2),dots,s_(9))\operatorname{Diag}\left(g_{1}, g_{2}, \ldots, g_{9}, s_{1}, s_{2}, \ldots, s_{9}\right) ，所有特征值 g k , s k g k , s k g_(k),s_(k)g_{k}, s_{k} 均为实数。本征模态 z k ∗ , ζ k ∙ z k ∗ , ζ k ∙ z_(k)^(**),zeta_(k)^(∙)z_{k}^{*}, \zeta_{k}^{\bullet} 的解为
z k ∗ ( t ) = z k ∗ ( 0 ) exp ⁡ ( i g k t ) ; ζ k ∗ ( t ) = ζ k ∗ ( 0 ) exp ⁡ ( i s k t ) z k ∗ ( t ) = z k ∗ ( 0 ) exp ⁡ i g k t ; ζ k ∗ ( t ) = ζ k ∗ ( 0 ) exp ⁡ i s k t z_(k)^(**)(t)=z_(k)^(**)(0)exp(ig_(k)t);quadzeta_(k)^(**)(t)=zeta_(k)^(**)(0)exp(is_(k)t)z_{k}^{*}(t)=z_{k}^{*}(0) \exp \left(i g_{k} t\right) ; \quad \zeta_{k}^{*}(t)=\zeta_{k}^{*}(0) \exp \left(i s_{k} t\right).
The solutions in the elliptical variables z k , ζ k z k , ζ k z_(k),zeta_(k)z_{k}, \zeta_{k} are then given as linear combinations of the proper modes (19), and becomes quasiperiodic functions of the time t t tt, that is sums of pure periodic terms with independent frequencies.
然后，椭圆变量 z k , ζ k z k , ζ k z_(k),zeta_(k)z_{k}, \zeta_{k} 中的解作为适当模式 (19) 的线性组合给出，并成为时间 t t tt 的准周期函数，即具有独立频率的纯周期项之和。
z k = ∑ k = 1 9 α k j e i g k t ; ζ k = ∑ k = 1 9 β k j e i s k t z k = ∑ k = 1 9   α k j e i g k t ; ζ k = ∑ k = 1 9   β k j e i s k t z_(k)=sum_(k=1)^(9)alpha_(kj)e^(ig_(k)t);quadzeta_(k)=sum_(k=1)^(9)beta_(kj)e^(is_(k)t)z_{k}=\sum_{k=1}^{9} \alpha_{k j} \mathrm{e}^{i g_{k} t} ; \quad \zeta_{k}=\sum_{k=1}^{9} \beta_{k j} \mathrm{e}^{i s_{k} t}.
Lagrange ( 1781 , 1782 ) ( 1781 , 1782 ) (1781,1782)(1781,1782) actually did this computation for a planetary system including Mercury, Venus, Earth, Mars, Jupiter and Saturn. Pontécoulant (1834) extended it to the 7 known planets of the time, but with several errors in the numerical constants. The first precise solution for the Earth was given by Le Verrier (1856) for 7 planets, and was later on extended by Stockwell (1873) with the addition of Neptune. The validity of the linear approximation (17) was also questioned and subsequent improvements introduced some additional contribution of higher degree in eccentricity and inclination in Eq. (17), as well as terms of higher order with respect to the planetary masses in the averaged equations (Hill 1897; Brouwer & Van Woerkom 1950; Bretagnon 1974; Duriez 1977; Laskar 1985). With the addition of additional non linear terms B ( x , x ¯ ) B ( x , x ¯ ) B(x, bar(x))B(x, \bar{x}) (of degree ≥ 3 ≥ 3 >= 3\geq 3 in eccentricity and inclination), the secular equations
拉格朗 ( 1781 , 1782 ) ( 1781 , 1782 ) (1781,1782)(1781,1782) 实际上对包括水星、金星、地球、火星、木星和土星在内的行星系统进行了这一计算。庞特库朗 (1834) 将其扩展到当时已知的 7 颗行星，但在数值常数方面存在一些错误。勒威耶 (1856) 给出了地球 7 颗行星的第一个精确解，后来斯托克韦尔 (1873) 扩展了这一解，加入了海王星。线性近似 (17) 的有效性也受到质疑，随后的改进在方程 (17) 中引入了一些更高阶的偏心率和倾角贡献，以及在平均方程中引入了关于行星质量的高阶项（Hill 1897；Brouwer & Van Woerkom 1950；Bretagnon 1974；Duriez 1977；Laskar 1985）。通过添加额外的非线性项 B ( x , x ¯ ) B ( x , x ¯ ) B(x, bar(x))B(x, \bar{x}) （偏心率和倾角为 ≥ 3 ≥ 3 >= 3\geq 3 度），长期方程
d [ x ] d t = i A [ x ] + B ( x , x ¯ ) d [ x ] d t = i A [ x ] + B ( x , x ¯ ) (d[x])/(dt)=iA[x]+B(x, bar(x))\frac{\mathrm{d}[x]}{\mathrm{d} t}=i A[x]+B(x, \bar{x})
are no longer integrable. Nevertheless, one can formally perform a Birkhoff normalization of Eq. (23) and after truncation, obtain some quasiperiodic expressions
不再可积。然而，我们可以对方程 (23) 进行 Birkhoff 正则化，并在截断后得到一些准周期表达式
z j ∗ ( t ) = ∑ k = 1 9 a j k e i g k t + ∑ ( k ) a j ( k ) e i < ( k ) , v > t z j ∗ ( t ) = ∑ k = 1 9   a j k e i g k t + ∑ ( k )   a j ( k ) e i < ( k ) , v > t z_(j)^(**)(t)=sum_(k=1)^(9)a_(jk)e^(ig_(k)t)+sum_((k))a_(j(k))e^(i < (k),v > t)z_{j}^{*}(t)=\sum_{k=1}^{9} a_{j k} \mathrm{e}^{i g_{k} t}+\sum_{(k)} a_{j(k)} \mathrm{e}^{i<(k), v>t}
ζ j ∗ ( t ) = ∑ k = 1 9 b j k e i s k t + ∑ ( k ) b j ( k ) e i < ( k ) , v > t ζ j ∗ ( t ) = ∑ k = 1 9   b j k e i s k t + ∑ ( k )   b j ( k ) e i < ( k ) , v > t zeta_(j)^(**)(t)=sum_(k=1)^(9)b_(jk)e^(is_(k)t)+sum_((k))b_(j(k))e^(i < (k),v > t)\zeta_{j}^{*}(t)=\sum_{k=1}^{9} b_{j k} \mathrm{e}^{i s_{k} t}+\sum_{(k)} b_{j(k)} \mathrm{e}^{i<(k), v>t}
where ( k ) = ( k 1 , … , k 18 ) ∈ Z 18 , v = ( g 1 , … , g 9 , s 1 , … , s 9 ) ( k ) = k 1 , … , k 18 ∈ Z 18 , v = g 1 , … , g 9 , s 1 , … , s 9 (k)=(k_(1),dots,k_(18))inZ^(18),v=(g_(1),dots,g_(9),s_(1),dots,s_(9))(k)=\left(k_{1}, \ldots, k_{18}\right) \in \mathbb{Z}^{18}, v=\left(g_{1}, \ldots, g_{9}, s_{1}, \ldots, s_{9}\right), and ⟨ ( k ) , v ⟩ = k 1 g 1 + … + k 9 g 9 + k 10 s 1 + … + k 18 s 9 ⟨ ( k ) , v ⟩ = k 1 g 1 + … + k 9 g 9 + k 10 s 1 + … + k 18 s 9 (:(k),v:)=k_(1)g_(1)+dots+k_(9)g_(9)+k_(10)s_(1)+dots+k_(18)s_(9)\langle(k), v\rangle=k_{1} g_{1}+\ldots+k_{9} g_{9}+k_{10} s_{1}+\ldots+k_{18} s_{9}. The quasiperiodic expressions z j ∗ ( t ) , ζ j ∗ ( t ) z j ∗ ( t ) , ζ j ∗ ( t ) z_(j)^(**)(t),zeta_(j)^(**)(t)z_{j}^{*}(t), \zeta_{j}^{*}(t), have been computed at various orders and degrees in eccentricity and inclination of expansion (Hill 1897; Bretagnon 1974; Duriez 1977; Laskar 1984, 1985), but as forecasted by Poincaré (1893), it was found that when all the main planets are taken into account, these formal expansions do not converge (Laskar 1984), for initial conditions close to the initial conditions of the Solar System, even at the lowest orders, and thus will not provide good approximation of the Solar System motion over very long time. In other terms, it appears that the actual solution of the Solar System is not close to a KAM tori (see Arnold et al. 1988) of quasiperiodic solutions.
其中 ( k ) = ( k 1 , … , k 18 ) ∈ Z 18 , v = ( g 1 , … , g 9 , s 1 , … , s 9 ) ( k ) = k 1 , … , k 18 ∈ Z 18 , v = g 1 , … , g 9 , s 1 , … , s 9 (k)=(k_(1),dots,k_(18))inZ^(18),v=(g_(1),dots,g_(9),s_(1),dots,s_(9))(k)=\left(k_{1}, \ldots, k_{18}\right) \in \mathbb{Z}^{18}, v=\left(g_{1}, \ldots, g_{9}, s_{1}, \ldots, s_{9}\right) 和 ⟨ ( k ) , v ⟩ = k 1 g 1 + … + k 9 g 9 + k 10 s 1 + … + k 18 s 9 ⟨ ( k ) , v ⟩ = k 1 g 1 + … + k 9 g 9 + k 10 s 1 + … + k 18 s 9 (:(k),v:)=k_(1)g_(1)+dots+k_(9)g_(9)+k_(10)s_(1)+dots+k_(18)s_(9)\langle(k), v\rangle=k_{1} g_{1}+\ldots+k_{9} g_{9}+k_{10} s_{1}+\ldots+k_{18} s_{9} 。准周期表达式 z j ∗ ( t ) , ζ j ∗ ( t ) z j ∗ ( t ) , ζ j ∗ ( t ) z_(j)^(**)(t),zeta_(j)^(**)(t)z_{j}^{*}(t), \zeta_{j}^{*}(t) 已在不同的偏心率和膨胀倾角阶数和次数下计算过（Hill 1897；Bretagnon 1974；Duriez 1977；Laskar 1984，1985），但正如庞加莱（1893）所预测的那样，当考虑到所有主要行星时，这些正式的展开式在接近太阳系初始条件的情况下，即使在最低阶的情况下也不会收敛（Laskar 1984），因此无法很好地近似太阳系在很长一段时间内的运动。换句话说，太阳系的实际解似乎并不接近准周期解的 KAM 环面（参见 Arnold 等人 1988）。
Over a finite time of a few millions of years, one can still numerically integrate the equation of motion, and then search through frequency analysis (Laskar 1990, 1999, 2003) for a quasiperiodic approximation of the solution of the form (23). The non regular behavior of the solution will induce a drift of the frequencies in time, related to the chaotic diffusion of the trajectories (Laskar 1990, 1993, 1999).
在几百万年的有限时间内，人们仍然可以对运动方程进行数值积分，然后通过频率分析（Laskar 1990, 1999, 2003）寻找形式（23）解的准周期近似。解的非规则行为将导致频率随时间漂移，这与轨迹的混沌扩散有关（Laskar 1990, 1993, 1999）。
In Table 3, the fundamental frequencies of the solution La2004 have been computed by frequency analysis over 20 Myr of the proper modes ( z 1 ∙ , … , z 4 ∙ , ζ 1 ∙ , … , ζ 4 ∙ ) z 1 ∙ , … , z 4 ∙ , ζ 1 ∙ , … , ζ 4 ∙ (z_(1)^(∙),dots,z_(4)^(∙),zeta_(1)^(∙),dots,zeta_(4)^(∙))\left(z_{1}^{\bullet}, \ldots, z_{4}^{\bullet}, \zeta_{1}^{\bullet}, \ldots, \zeta_{4}^{\bullet}\right) related to the inner planets (Mercury, Venus, Earth, Mars), and over 50 Myr for the proper modes ( z 5 ∙ , … , z 9 ∙ , ζ 5 ∙ , … , ζ 9 ∙ ) z 5 ∙ , … , z 9 ∙ , ζ 5 ∙ , … , ζ 9 ∙ (z_(5)^(∙),dots,z_(9)^(∙),zeta_(5)^(∙),dots,zeta_(9)^(∙))\left(z_{5}^{\bullet}, \ldots, z_{9}^{\bullet}, \zeta_{5}^{\bullet}, \ldots, \zeta_{9}^{\bullet}\right) associated to the outer planets (Jupiter, Saturn, Uranus, Neptune, Pluto) (Table 3). This difference reflects the fact that the main
在表 3 中，La2004 解的基频是通过对内行星（水星、金星、地球、火星）相关固有模式 ( z 1 ∙ , … , z 4 ∙ , ζ 1 ∙ , … , ζ 4 ∙ ) z 1 ∙ , … , z 4 ∙ , ζ 1 ∙ , … , ζ 4 ∙ (z_(1)^(∙),dots,z_(4)^(∙),zeta_(1)^(∙),dots,zeta_(4)^(∙))\left(z_{1}^{\bullet}, \ldots, z_{4}^{\bullet}, \zeta_{1}^{\bullet}, \ldots, \zeta_{4}^{\bullet}\right) 进行 20 Myr 以上频率分析，以及对外行星（木星、土星、天王星、海王星、冥王星）相关固有模式 ( z 5 ∙ , … , z 9 ∙ , ζ 5 ∙ , … , ζ 9 ∙ ) z 5 ∙ , … , z 9 ∙ , ζ 5 ∙ , … , ζ 9 ∙ (z_(5)^(∙),dots,z_(9)^(∙),zeta_(5)^(∙),dots,zeta_(9)^(∙))\left(z_{5}^{\bullet}, \ldots, z_{9}^{\bullet}, \zeta_{5}^{\bullet}, \ldots, \zeta_{9}^{\bullet}\right) 进行 50 Myr 以上频率分析计算得出的（表 3）。这种差异反映了以下事实：
source of chaotic behavior comes from overlap of secular resonances in the inner planet system, and that the outer planets secular systems is more regular than the inner one. Actually, in Cols. 2 and 3 of Table 3, we give the maximum observed in the variation of the secular frequencies over 100 and 250 Myr , in both positive and negative time, for our nominal solution La2004. These figures are in very good agreement with the equivalent quantities given in Table X of (Laskar 1990) for the integration of the secular equations. The number of digits of the fundamental frequencies in Table 3 is different for each frequency and reflects the stability of the considered frequency with time. Additionally, we have given in Figs. 9 and 10 the actual evolution of the secular frequencies obtained with a sliding window with a step size of 1 Myr (compare to Laskar 1990, Figs. 8 and 9).
混沌行为的根源来自内行星系统中长期共振的重叠，而外行星的长期系统比内行星的长期系统更有规律。实际上，在表 3 的第 2 列和第 3 列中，我们给出了在 100 和 250 Myr 期间观察到的长期频率变化的最大值（正时间和负时间），这是我们的标称解 La2004。这些数字与 (Laskar 1990) 表 X 中给出的长期方程积分的等效量非常吻合。表 3 中基频的位数对于每个频率都是不同的，这反映了所考虑频率随时间的稳定性。此外，我们在图 9 和图 10 中给出了使用步长为 1 Myr 的滑动窗口获得的长期频率的实际演变（与 Laskar 1990 的图 8 和图 9 比较）。

8. Analytical approximations
8. 解析近似

8.1. Orbital motion  8.1 轨道运动

It is interesting for practical use to have an analytical expression for the main orbital quantities of the Earth. From the numerical values of La2004, we have performed a frequency analysis (Laskar 1990, 1999, 2003) in order to obtain a quasiperiodic approximation of the solutions over a few Myr. In order to be consistent with the remaining part of the paper, we chose a time interval covering 20 Myr . As we are mostly interested over negative time, we made the analysis from -15 Myr to +5 Myr , as usually, the precision of the approximation decreases at the edges of the time interval. The analysis of the eccentricity variables z = e exp ⁡ i ϖ z = e exp ⁡ i ϖ z=e exp iϖz=e \exp i \varpi is given in Table 4. z z zz is obtained as
对地球主要轨道量的解析表达式在实际应用中很有意思。我们根据 La2004 的数值进行了频率分析（Laskar 1990、1999、2003），以获得几百万年的准周期近似解。为了与本文的其余部分保持一致，我们选择了 20 百万年的时间间隔。由于我们主要对负时间感兴趣，因此我们从 -15 百万年到 +5 百万年进行了分析，因为通常近似的精度在时间间隔的边缘会降低。偏心率变量 z = e exp ⁡ i ϖ z = e exp ⁡ i ϖ z=e exp iϖz=e \exp i \varpi 的分析在表 4 中给出。 z z zz 的计算公式为
z = ∑ k = 1 26 b k e i ( μ k t + φ k ) z = ∑ k = 1 26   b k e i μ k t + φ k z=sum_(k=1)^(26)b_(k)e^(i(mu_(k)t+varphi_(k)))z=\sum_{k=1}^{26} b_{k} \mathrm{e}^{i\left(\mu_{k} t+\varphi_{k}\right)}.
It should be noted that in the present complex form, the approximation is much more precise than what would give an equivalent quasiperiodic approximation of the eccentricity. Even with twice the number of periodic terms, the direct approximation of the eccentricity is not as good as the approximation obtained in complex variables. Moreover, here we obtain at the same time the approximation of the longitude of perihelion of the Earth from the fixed J2000 equinox ( ϖ ϖ ϖ\varpi ). The comparison of this approximation, limited to only 26 terms, to the actual solution for the eccentricity of the Earth, is plotted in Fig. 12 (top).
需要注意的是，在目前的复数形式下，该近似值比等效的准周期偏心率近似值精确得多。即使周期项的数量增加一倍，直接近似偏心率也不如用复变量获得的近似值好。此外，我们同时获得了从固定的 J2000 春分点（ ϖ ϖ ϖ\varpi ）获得的地球近日点经度的近似值。图 12（顶部）绘制了该近似值（仅限于 26 项）与地球偏心率的实际解的比较。
Nevertheless, for information, we will give also here the leading terms in the expansion of the eccentricity as they may be useful in paleoclimate studies (Table 6). The three leading terms in eccentricity are well known in paleoclimate studies g 2 − g 5 g 2 − g 5 g_(2)-g_(5)g_{2}-g_{5} ( 405 kyr period), g 4 − g 5 g 4 − g 5 g_(4)-g_(5)g_{4}-g_{5} ( 95 kyr period), and g 4 − g 2 g 4 − g 2 g_(4)-g_(2)g_{4}-g_{2} (124 kyr period).
尽管如此，为了提供信息，我们还将在此列出偏心率扩展中的主项，因为它们可能在古气候研究中有用（表6）。偏心率中的三个主项在古气候研究中是众所周知的： g 2 − g 5 g 2 − g 5 g_(2)-g_(5)g_{2}-g_{5} （405 千年周期）、 g 4 − g 5 g 4 − g 5 g_(4)-g_(5)g_{4}-g_{5} （95 千年周期）和 g 4 − g 2 g 4 − g 2 g_(4)-g_(2)g_{4}-g_{2} （124 千年周期）。
The solution for the inclination variables ζ = sin ⁡ i / 2 exp ⁡ i Ω ζ = sin ⁡ i / 2 exp ⁡ i Ω zeta=sin i//2exp i Omega\zeta=\sin i / 2 \exp i \Omega where i , Ω i , Ω i,Omegai, \Omega are the Earth inclination and longitude of node with respect to the fixed ecliptic and equinox J2000, limited to 24 quasiperiodic terms
倾角变量 ζ = sin ⁡ i / 2 exp ⁡ i Ω ζ = sin ⁡ i / 2 exp ⁡ i Ω zeta=sin i//2exp i Omega\zeta=\sin i / 2 \exp i \Omega 的解，其中 i , Ω i , Ω i,Omegai, \Omega 是相对于固定黄道和春分点 J2000 的地球倾角和节点经度，限制为 24 个准周期项
ζ = ∑ k = 1 24 a k e i ( v k t + ϕ k ) ζ = ∑ k = 1 24   a k e i v k t + ϕ k zeta=sum_(k=1)^(24)a_(k)e^(i(v_(k)t+phi_(k)))\zeta=\sum_{k=1}^{24} a_{k} \mathrm{e}^{\mathrm{i}\left(v_{k} t+\phi_{k}\right)},
Fig. 9. Variation of the secular frequencies g 1 − 9 g 1 − 9 g_(1-9)g_{1-9} from -250 to +250 Myr. The frequencies are computed over 20 Myr for g 1 − 4 g 1 − 4 g_(1-4)g_{1-4} and over 50 Myr for g 5 − 9 g 5 − 9 g_(5-9)g_{5-9}, after transformation of elliptical elements to proper modes (Laskar 1990).
图 9. 长期频率 g 1 − 9 g 1 − 9 g_(1-9)g_{1-9} 在-250 百万年到+250 百万年之间的变化。在将椭圆元素变换为本征模态后，计算出 g 1 − 4 g 1 − 4 g_(1-4)g_{1-4} 在 20 百万年内的频率，以及 g 5 − 9 g 5 − 9 g_(5-9)g_{5-9} 在 50 百万年内的频率（Laskar 1990）。
is given in Table 5. The comparison of this quasiperiodic approximation with the complete solution is given in Fig. 12 (bottom). It should be noted that in Table 5, the first 22 terms are the terms of largest amplitude in the frequency decomposition of the inclination variable ζ ζ zeta\zeta, but the two last terms, of frequency v 23 v 23 v_(23)v_{23} and v 24 v 24 v_(24)v_{24} are of much smaller amplitude. They have been kept in the solution, as they are close to the resonance with the main precession frequency.
列于表5。图12（底部）给出了该准周期近似与完整解的比较。需要注意的是，表5中，前22项是倾角变量 ζ ζ zeta\zeta 频率分解中振幅最大的项，但最后两个频率为 v 23 v 23 v_(23)v_{23} 和 v 24 v 24 v_(24)v_{24} 的项振幅要小得多。由于它们接近主进动频率的共振点，因此保留在解中。

Fig. 10. Variation of the secular frequencies s 1 − 9 s 1 − 9 s_(1-9)s_{1-9} from -250 to +250 Myr . The frequencies are computed over 20 Myr for s 1 − 4 s 1 − 4 s_(1-4)s_{1-4} and over 50 Myr for s 6 − 9 s 6 − 9 s_(6-9)s_{6-9}, after transformation of elliptical elements to proper modes (Laskar 1990).
图 10。长期频率 s 1 − 9 s 1 − 9 s_(1-9)s_{1-9} 从-250 百万年到+250 百万年的变化。在将椭圆元素变换为本征模式后，计算出 s 1 − 4 s 1 − 4 s_(1-4)s_{1-4} 在 20 百万年内的频率，以及 s 6 − 9 s 6 − 9 s_(6-9)s_{6-9} 在 50 百万年内的频率（Laskar 1990）。

8.2. Obliquity and precession
8.2 倾角和进动

The solutions for precession and obliquity are obtained through the precession Eq. (6). In absence of planetary perturbations, the obliquity is constant ( ε = ε 0 ) ε = ε 0 (epsi=epsi_(0))\left(\varepsilon=\varepsilon_{0}\right), and
通过进动方程(6)可得到进动和倾角的解。当没有行星摄动时，倾角为常数 ( ε = ε 0 ) ε = ε 0 (epsi=epsi_(0))\left(\varepsilon=\varepsilon_{0}\right) ，
d ψ d t = α cos ⁡ ε 0 d ψ d t = α cos ⁡ ε 0 (dpsi)/(dt)=alpha cos epsi_(0)\frac{\mathrm{d} \psi}{\mathrm{d} t}=\alpha \cos \varepsilon_{0}.
We have thus ψ = ψ 0 + p t ψ = ψ 0 + p t psi=psi_(0)+pt\psi=\psi_{0}+p t where p = α cos ⁡ ε 0 p = α cos ⁡ ε 0 p=alpha cos epsi_(0)p=\alpha \cos \varepsilon_{0} is the precession frequency. This can be considered as a solution of order zero.
因此，我们得到 ψ = ψ 0 + p t ψ = ψ 0 + p t psi=psi_(0)+pt\psi=\psi_{0}+p t ，其中 p = α cos ⁡ ε 0 p = α cos ⁡ ε 0 p=alpha cos epsi_(0)p=\alpha \cos \varepsilon_{0} 是进动频率。这可以被视为零阶解。
Table 4. Frequency decomposition of z = e exp ⁡ i ϖ z = e exp ⁡ i ϖ z=e exp iϖz=e \exp i \varpi for the Earth on the time interval [ − 15 , + 5 ] Myr [ − 15 , + 5 ] Myr [-15,+5]Myr[-15,+5] \mathrm{Myr} (Eq. (25)).
表 4. 时间间隔 [ − 15 , + 5 ] Myr [ − 15 , + 5 ] Myr [-15,+5]Myr[-15,+5] \mathrm{Myr} 内地球的 z = e exp ⁡ i ϖ z = e exp ⁡ i ϖ z=e exp iϖz=e \exp i \varpi 频率分解（公式 (25)）。
n n nn μ k ( ′ ′ / yr ) μ k ′ ′ / yr mu_(k)(^('')//yr)\mu_{k}\left({ }^{\prime \prime} / \mathrm{yr}\right) b k b k b_(k)b_{k} φ k ( φ k ( varphi_(k)(\varphi_{k}( degree ) ) ))   φ k ( φ k ( varphi_(k)(\varphi_{k}( 度 ) ) ))
1 g 5 g 5 g_(5)g_{5} 4.257564 0.018986 30.739
2 g 2 g 2 g_(2)g_{2} 7.456665 0.016354 -157.801
3 g 4 g 4 g_(4)g_{4} 17.910194 0.013055 140.577
4 g 3 g 3 g_(3)g_{3} 17.366595 0.008849 -55.885
5 g 1 g 1 g_(1)g_{1} 5.579378 0.004248 77.107
6 17.112064 0.002742 47.965
7 17.654560 0.002386 61.025
8 6.954694 0.001796 -159.137
9 2 g 3 − g 4 2 g 3 − g 4 2g_(3)-g_(4)2 g_{3}-g_{4} 16.822731 0.001908 105.698
10 g 6 g 6 g_(6)g_{6} 28.244959 0.001496 127.609
11 18.452681 0.001325 156.113
12 5.446730 0.001266 -124.155
13 18.203535 0.001165 -92.112
14 7.325741 0.001304 -175.144
15 7.060100 0.001071 2.423
16 7.587926 0.000970 34.169
17 5.723726 0.001002 112.564
18 16.564805 0.000936 24.637
19 16.278532 0.000781 -90.999
20 18.999613 0.000687 0.777
21 3.087424 0.000575 120.376
22 17.223297 0.000577 -113.456
23 6.778196 0.000651 117.184
24 6.032715 0.000416 174.987
25 6.878659 0.000497 95.424
26 5.315656 0.000392 33.255
n mu_(k)(^('')//yr) b_(k) varphi_(k)( degree ) 1 g_(5) 4.257564 0.018986 30.739 2 g_(2) 7.456665 0.016354 -157.801 3 g_(4) 17.910194 0.013055 140.577 4 g_(3) 17.366595 0.008849 -55.885 5 g_(1) 5.579378 0.004248 77.107 6 17.112064 0.002742 47.965 7 17.654560 0.002386 61.025 8 6.954694 0.001796 -159.137 9 2g_(3)-g_(4) 16.822731 0.001908 105.698 10 g_(6) 28.244959 0.001496 127.609 11 18.452681 0.001325 156.113 12 5.446730 0.001266 -124.155 13 18.203535 0.001165 -92.112 14 7.325741 0.001304 -175.144 15 7.060100 0.001071 2.423 16 7.587926 0.000970 34.169 17 5.723726 0.001002 112.564 18 16.564805 0.000936 24.637 19 16.278532 0.000781 -90.999 20 18.999613 0.000687 0.777 21 3.087424 0.000575 120.376 22 17.223297 0.000577 -113.456 23 6.778196 0.000651 117.184 24 6.032715 0.000416 174.987 25 6.878659 0.000497 95.424 26 5.315656 0.000392 33.255 | $n$ | | $\mu_{k}\left({ }^{\prime \prime} / \mathrm{yr}\right)$ | $b_{k}$ | $\varphi_{k}($ degree $)$ | | ---: | :---: | ---: | ---: | ---: | | 1 | $g_{5}$ | 4.257564 | 0.018986 | 30.739 | | 2 | $g_{2}$ | 7.456665 | 0.016354 | -157.801 | | 3 | $g_{4}$ | 17.910194 | 0.013055 | 140.577 | | 4 | $g_{3}$ | 17.366595 | 0.008849 | -55.885 | | 5 | $g_{1}$ | 5.579378 | 0.004248 | 77.107 | | 6 | | 17.112064 | 0.002742 | 47.965 | | 7 | | 17.654560 | 0.002386 | 61.025 | | 8 | | 6.954694 | 0.001796 | -159.137 | | 9 | $2 g_{3}-g_{4}$ | 16.822731 | 0.001908 | 105.698 | | 10 | $g_{6}$ | 28.244959 | 0.001496 | 127.609 | | 11 | | 18.452681 | 0.001325 | 156.113 | | 12 | | 5.446730 | 0.001266 | -124.155 | | 13 | | 18.203535 | 0.001165 | -92.112 | | 14 | | 7.325741 | 0.001304 | -175.144 | | 15 | | 7.060100 | 0.001071 | 2.423 | | 16 | | 7.587926 | 0.000970 | 34.169 | | 17 | | 5.723726 | 0.001002 | 112.564 | | 18 | | 16.564805 | 0.000936 | 24.637 | | 19 | | 16.278532 | 0.000781 | -90.999 | | 20 | | 18.999613 | 0.000687 | 0.777 | | 21 | | 3.087424 | 0.000575 | 120.376 | | 22 | | 17.223297 | 0.000577 | -113.456 | | 23 | | 6.778196 | 0.000651 | 117.184 | | 24 | | 6.032715 | 0.000416 | 174.987 | | 25 | | 6.878659 | 0.000497 | 95.424 | | 26 | | 5.315656 | 0.000392 | 33.255 | | | | | | |
If we keep only the terms of degree one in inclination in Eqs. (6), we obtain the solution of order one
如果我们只保留方程（6）中倾斜度为一阶的项，我们得到一阶解
d ε d t = 2 ( p ˙ sin ⁡ ψ − q ˙ cos ⁡ ψ ) = 2 R e ( ζ ˙ i ψ ψ ) d ε d t = 2 ( p ˙ sin ⁡ ψ − q ˙ cos ⁡ ψ ) = 2 R e ζ ˙ i ψ ψ (depsi)/(dt)=2(p^(˙)sin psi-q^(˙)cos psi)=2Re(zeta^(˙)^(i psi psi))\frac{\mathrm{d} \varepsilon}{\mathrm{d} t}=2(\dot{p} \sin \psi-\dot{q} \cos \psi)=2 \mathcal{R} e\left(\dot{\zeta}^{i \psi \psi}\right).
With the quasiperiodic approximation
利用准周期近似
ζ = ∑ k = 1 N a k e i ( v k t + ϕ k ) ζ = ∑ k = 1 N   a k e i v k t + ϕ k zeta=sum_(k=1)^(N)a_(k)e^(i(v_(k)t+phi_(k)))\zeta=\sum_{k=1}^{N} a_{k} \mathrm{e}^{i\left(v_{k} t+\phi_{k}\right)},
of Table 5, the first order solution in obliquity will be a similar quasiperiodic function
表 5 中，倾角的一阶解将是类似的准周期函数
ε = ε 0 + 2 ∑ k = 1 N a k v k v k + p cos ⁡ ( ( v k + p ) t + ϕ k + ψ 0 ) ε = ε 0 + 2 ∑ k = 1 N   a k v k v k + p cos ⁡ v k + p t + ϕ k + ψ 0 epsi=epsi_(0)+2sum_(k=1)^(N)(a_(k)v_(k))/(v_(k)+p)cos((v_(k)+p)t+phi_(k)+psi_(0))\varepsilon=\varepsilon_{0}+2 \sum_{k=1}^{N} \frac{a_{k} v_{k}}{v_{k}+p} \cos \left(\left(v_{k}+p\right) t+\phi_{k}+\psi_{0}\right),
with a similar form for the solution of precession ψ ψ psi\psi. This first order approximation will only be valid when we are far from resonance, that is | v k + p | ≫ 0 v k + p ≫ 0 |v_(k)+p|≫0\left|v_{k}+p\right| \gg 0. This is not the case for the last term v 23 = s 6 + g 5 − g 6 = − 50.336259 arcsec yr − 1 v 23 = s 6 + g 5 − g 6 = − 50.336259 arcsec yr − 1 v_(23)=s_(6)+g_(5)-g_(6)=-50.336259 arcsecyr^(-1)v_{23}=s_{6}+g_{5}-g_{6}=-50.336259 \operatorname{arcsec} \mathrm{yr}^{-1} in Table 5, as p ≈ 50.484 arcsec yr − 1 p ≈ 50.484 arcsec yr − 1 p~~50.484 arcsecyr^(-1)p \approx 50.484 \operatorname{arcsec} \mathrm{yr}^{-1}. This resonant term can lead to a
对于进动 ψ ψ psi\psi 的解，其形式类似。这种一阶近似只有在远离共振时才有效，即 | v k + p | ≫ 0 v k + p ≫ 0 |v_(k)+p|≫0\left|v_{k}+p\right| \gg 0 。表5中的最后一项 v 23 = s 6 + g 5 − g 6 = − 50.336259 arcsec yr − 1 v 23 = s 6 + g 5 − g 6 = − 50.336259 arcsec yr − 1 v_(23)=s_(6)+g_(5)-g_(6)=-50.336259 arcsecyr^(-1)v_{23}=s_{6}+g_{5}-g_{6}=-50.336259 \operatorname{arcsec} \mathrm{yr}^{-1} 并非如此，因为 p ≈ 50.484 arcsec yr − 1 p ≈ 50.484 arcsec yr − 1 p~~50.484 arcsecyr^(-1)p \approx 50.484 \operatorname{arcsec} \mathrm{yr}^{-1} 。这个共振项可以导致
Fig.11. Variation of the semi-major axis of the Earth-Moon barycenter (in AU) from -250 to +250 Myr .
图11. 地月质心半长轴的变化（以天文单位计），从-250百万年到+250百万年。
Table 5. Frequency decomposition of ζ = sin ⁡ i / 2 exp ⁡ i Ω ζ = sin ⁡ i / 2 exp ⁡ i Ω zeta=sin i//2exp i Omega\zeta=\sin i / 2 \exp i \Omega for the Earth on the time interval [ − 15 , + 5 ] Myr [ − 15 , + 5 ] Myr [-15,+5]Myr[-15,+5] \operatorname{Myr} (Eq. (26)).
表 5. 时间间隔 [ − 15 , + 5 ] Myr [ − 15 , + 5 ] Myr [-15,+5]Myr[-15,+5] \operatorname{Myr} 内地球的 ζ = sin ⁡ i / 2 exp ⁡ i Ω ζ = sin ⁡ i / 2 exp ⁡ i Ω zeta=sin i//2exp i Omega\zeta=\sin i / 2 \exp i \Omega 频率分解（公式 (26)）。
n n nn v k ( ′ ′ / yr ) v k ′ ′ / yr v_(k)(^('')//yr)v_{k}\left({ }^{\prime \prime} / \mathrm{yr}\right) a k a k a_(k)a_{k} ϕ k ϕ k phi_(k)\phi_{k} (degree)   ϕ k ϕ k phi_(k)\phi_{k} （度）
1 s 5 s 5 s_(5)s_{5} -0.000001 0.01377449 107.581
2 s 3 s 3 s_(3)s_{3} -18.845166 0.00870353 -111.310
3 s 1 s 1 s_(1)s_{1} -5.605919 0.00479813 4.427
4 s 2 s 2 s_(2)s_{2} -7.050665 0.00350477 130.826
5 s 4 s 4 s_(4)s_{4} -17.758310 0.00401601 -77.666
6 -18.300842 0.00262820 -93.287
7 -7.178497 0.00189297 -65.292
8 -6.940312 0.00164641 -66.106
9 -6.817771 0.00155131 -46.436
10 -19.389947 0.00144918 50.158
11 -5.481324 0.00154505 23.486
12 s 6 s 6 s_(6)s_{6} -26.347880 0.00133215 127.306
13 -19.111118 0.00097602 159.658
14 -2.992612 0.00088810 140.098
15 -6.679492 0.00079348 -43.800
16 -5.835373 0.00073058 -160.927
17 -0.691879 0.00064133 23.649
18 -6.125755 0.00049598 143.759
19 -5.684257 0.00054320 138.236
20 -18.976417 0.00040739 -105.209
21 -7.290771 0.00041850 106.192
22 -5.189414 0.00033868 50.827
23 s 6 + g 5 − g 6 s 6 + g 5 − g 6 s_(6)+g_(5)-g_(6)s_{6}+g_{5}-g_{6} -50.336259 0.00000206 -150.693
24 -47.144010 0.00000023 -167.522
n v_(k)(^('')//yr) a_(k) phi_(k) (degree) 1 s_(5) -0.000001 0.01377449 107.581 2 s_(3) -18.845166 0.00870353 -111.310 3 s_(1) -5.605919 0.00479813 4.427 4 s_(2) -7.050665 0.00350477 130.826 5 s_(4) -17.758310 0.00401601 -77.666 6 -18.300842 0.00262820 -93.287 7 -7.178497 0.00189297 -65.292 8 -6.940312 0.00164641 -66.106 9 -6.817771 0.00155131 -46.436 10 -19.389947 0.00144918 50.158 11 -5.481324 0.00154505 23.486 12 s_(6) -26.347880 0.00133215 127.306 13 -19.111118 0.00097602 159.658 14 -2.992612 0.00088810 140.098 15 -6.679492 0.00079348 -43.800 16 -5.835373 0.00073058 -160.927 17 -0.691879 0.00064133 23.649 18 -6.125755 0.00049598 143.759 19 -5.684257 0.00054320 138.236 20 -18.976417 0.00040739 -105.209 21 -7.290771 0.00041850 106.192 22 -5.189414 0.00033868 50.827 23 s_(6)+g_(5)-g_(6) -50.336259 0.00000206 -150.693 24 -47.144010 0.00000023 -167.522| $n$ | | $v_{k}\left({ }^{\prime \prime} / \mathrm{yr}\right)$ | $a_{k}$ | $\phi_{k}$ (degree) | | ---: | ---: | ---: | ---: | ---: | | 1 | $s_{5}$ | -0.000001 | 0.01377449 | 107.581 | | 2 | $s_{3}$ | -18.845166 | 0.00870353 | -111.310 | | 3 | $s_{1}$ | -5.605919 | 0.00479813 | 4.427 | | 4 | $s_{2}$ | -7.050665 | 0.00350477 | 130.826 | | 5 | $s_{4}$ | -17.758310 | 0.00401601 | -77.666 | | 6 | | -18.300842 | 0.00262820 | -93.287 | | 7 | | -7.178497 | 0.00189297 | -65.292 | | 8 | | -6.940312 | 0.00164641 | -66.106 | | 9 | | -6.817771 | 0.00155131 | -46.436 | | 10 | | -19.389947 | 0.00144918 | 50.158 | | 11 | | -5.481324 | 0.00154505 | 23.486 | | 12 | $s_{6}$ | -26.347880 | 0.00133215 | 127.306 | | 13 | | -19.111118 | 0.00097602 | 159.658 | | 14 | | -2.992612 | 0.00088810 | 140.098 | | 15 | | -6.679492 | 0.00079348 | -43.800 | | 16 | | -5.835373 | 0.00073058 | -160.927 | | 17 | | -0.691879 | 0.00064133 | 23.649 | | 18 | | -6.125755 | 0.00049598 | 143.759 | | 19 | | -5.684257 | 0.00054320 | 138.236 | | 20 | | -18.976417 | 0.00040739 | -105.209 | | 21 | | -7.290771 | 0.00041850 | 106.192 | | 22 | | -5.189414 | 0.00033868 | 50.827 | | 23 | $s_{6}+g_{5}-g_{6}$ | -50.336259 | 0.00000206 | -150.693 | | 24 | | -47.144010 | 0.00000023 | -167.522 |
significative change in the obliquity variations, and have been studied in detail in Laskar et al. (1993).
倾角变化的显著变化，已由 Laskar 等人 (1993) 进行了详细研究。
Another difficulty, arising with the obliquity solution, is due to the dissipation in the Earth-Moon system, which induces a significant variation of the precession frequency with time, as the Earth rotation ω ω omega\omega slows down, and the Earth-Moon distance a M a M a_(M)a_{\mathrm{M}} increases.
倾角解决方案产生的另一个困难是由于地月系统的耗散，随着地球自转 ω ω omega\omega 减慢和地月距离 a M a M a_(M)a_{\mathrm{M}} 增加，这会导致进动频率随时间发生显著变化。
We have plotted in Fig. 14 the evolution of the obliquity of the Earth from -250 to +250 Myr . The effect of the tidal dissipation is clearly visible on this timescale, and there is a general increase of the obliquity from -250 Myr to present time. In positive time, there is an obvious singularity in the obliquity, with a decrease of about 0.4 degree. This results from the
图14绘制了地球自转轴倾角从-2.5亿年到+2.5亿年的变化。潮汐消散的影响在这个时间尺度上清晰可见，自-2.5亿年到现在，自转轴倾角总体呈增加趋势。在正时间尺度上，自转轴倾角出现了一个明显的奇点，减小了约0.4度。这是由于
Fig. 12. (Top) Eccentricity of the Earth from -11 to +1 Myr. The solid curve is the solution La2004. Almost completely hidden behind the solid curve is a dashed curve representing the quasiperiodic approximation of Table 4 with 26 periodic terms. (Bottom) Inclination of the Earth from -11 to +1 Myr. The solid curve is the solution La2004. The dashed line (almost identical to the solid curve) is the quasiperiodic approximation of Table 5 with 24 periodic terms. On both plots, as the two curves are nearly identical, the difference of the two solutions is also plotted in dotted curve.
图 12. （上图）地球偏心率从 -11 百万年到 +1 百万年。实线为 La2004 解。几乎完全隐藏在实线后面的是一条虚线，表示表 4 的准周期近似，包含 26 个周期项。（下图）地球倾角从 -11 百万年到 +1 百万年。实线为 La2004 解。虚线（与实线几乎相同）表示表 5 的准周期近似，包含 24 个周期项。由于两幅图上的两条曲线几乎相同，因此两个解的差异也以虚线绘制。
crossing of the resonance with the term s 6 + g 5 − g 6 s 6 + g 5 − g 6 s_(6)+g_(5)-g_(6)s_{6}+g_{5}-g_{6} that was already described in a different context in Laskar et al. (1993). After this sudden decrease, the obliquity increases again, with roughly the same trend. An additional (but smaller singularity is observable around +150 Myr , that corresponds to the crossing of the resonance with the term v 24 v 24 v_(24)v_{24} (Table 5). The averaged obliquity can be well approximated in negative time with the expression
共振与项 s 6 + g 5 − g 6 s 6 + g 5 − g 6 s_(6)+g_(5)-g_(6)s_{6}+g_{5}-g_{6} 的交叉，这在 Laskar 等人（1993）的不同论述中已有描述。在这次突然下降之后，倾角再次增加，趋势大致相同。在 +150 Myr 左右可以观察到一个额外的（但较小的）奇点，它对应于共振与项 v 24 v 24 v_(24)v_{24} 的交叉（表 5）。平均倾角可以用以下表达式很好地近似于负时间
ε M = 23.270773 + 2.011295 T ε M = 23.270773 + 2.011295 T epsi_(M)=23.270773+2.011295 T\varepsilon_{\mathrm{M}}=23.270773+2.011295 T
where ε M ε M epsi_(M)\varepsilon_{\mathrm{M}} is in degree, and T T TT in billion of years (Gyr). In the interval [ + 20 Myr , + 130 Myr ] [ + 20 Myr , + 130 Myr ] [+20Myr,+130Myr][+20 \mathrm{Myr},+130 \mathrm{Myr}] this approximation becomes
其中 ε M ε M epsi_(M)\varepsilon_{\mathrm{M}} 表示度， T T TT 表示十亿年 (Gyr)。在区间 [ + 20 Myr , + 130 Myr ] [ + 20 Myr , + 130 Myr ] [+20Myr,+130Myr][+20 \mathrm{Myr},+130 \mathrm{Myr}] 中，这个近似值变为
ε M = 22.893171 + 1.952715 T ε M = 22.893171 + 1.952715 T epsi_(M)=22.893171+1.952715 T\varepsilon_{\mathrm{M}}=22.893171+1.952715 T
so the main effect of the crossing of the resonance is a decrease of 0.377602 degree of the obliquity. This effect can be easily understood with the analytical integration of a simple model, in terms of Fresnel integrals (see Appendix A).
因此，共振相交的主要效应是使倾角减小 0.377602 度。这一效应可以通过菲涅尔积分这个简单模型的解析积分轻松理解（参见附录 A）。

8.3. Approximation for the obliquity
8.3 倾角近似

As for the solutions of the orbital elements, we have obtained an approximation formula for the evolution of the obliquity of the Earth over the period from -15 Myr to +5 Myr . Obtaining this solution presents some additional difficulties because of the dissipative effects in the Earth-Moon system that induces a significative change of the precession frequency. The obliquity is obtained on the form
至于轨道要素的解，我们得到了地球自转轴倾角在-15百万年到+5百万年期间演变的近似公式。由于地月系统的耗散效应会导致进动频率发生显著变化，因此获得该解会面临一些额外的困难。自转轴倾角的形式如下
ε = ε 0 + ∑ k = 1 N a k ′ cos ⁡ ( ( v k ′ + p 1 t ) t + ϕ k ′ ) ε = ε 0 + ∑ k = 1 N   a k ′ cos ⁡ v k ′ + p 1 t t + ϕ k ′ epsi=epsi_(0)+sum_(k=1)^(N)a_(k)^(')cos((v_(k)^(')+p_(1)t)t+phi_(k)^('))\varepsilon=\varepsilon_{0}+\sum_{k=1}^{N} a_{k}^{\prime} \cos \left(\left(v_{k}^{\prime}+p_{1} t\right) t+\phi_{k}^{\prime}\right),
where ε 0 = 23.254500 ε 0 = 23.254500 epsi_(0)=23.254500\varepsilon_{0}=23.254500 degrees, and the coefficients a k ′ , v k ′ , ϕ k ′ a k ′ , v k ′ , ϕ k ′ a_(k)^('),v_(k)^('),phi_(k)^(')a_{k}^{\prime}, v_{k}^{\prime}, \phi_{k}^{\prime} are given in Table 7. The coefficent p 1 p 1 p_(1)p_{1} (Eq. (35)) results from the secular change of the precession frequency (34) due to tidal dissipation (see next).
其中 ε 0 = 23.254500 ε 0 = 23.254500 epsi_(0)=23.254500\varepsilon_{0}=23.254500 度，系数 a k ′ , v k ′ , ϕ k ′ a k ′ , v k ′ , ϕ k ′ a_(k)^('),v_(k)^('),phi_(k)^(')a_{k}^{\prime}, v_{k}^{\prime}, \phi_{k}^{\prime} 在表 7 中给出。系数 p 1 p 1 p_(1)p_{1} （公式 (35)）是由于潮汐耗散导致的进动频率 (34) 的长期变化（见下文）。
The comparison with the full solution is given in Fig. 15, over selected periods of time, from -11 Myr to the present.
图 15 给出了与完整解决方案的比较，比较范围涵盖了从 -11 Myr 到现在的选定时间段。
Fig. 13. a) Evolution of length of the day for the Earth, in hours, from -250 to +250 Myr . b) Residuals with the fit of the averaged solution with the polynomial expression (41).
图 13. a) 地球一天的长度演变，以小时为单位，从 -250 百万年到 +250 百万年。b) 与多项式表达式 (41) 拟合的平均解的残差。
In Table 7, the terms of frequency v k ′ v k ′ v_(k)^(')v_{k}^{\prime} are recognized as combinations of the form p + v k p + v k p+v_(k)p+v_{k}, where p p pp is the precession frequency, and v k v k v_(k)v_{k} a secular frequency from the secular motion of the orbital plane (from Table 5). The approximation is made over the interval [ − 15 Myr , + 2 Myr ] [ − 15 Myr , + 2 Myr ] [-15Myr,+2Myr][-15 \mathrm{Myr},+2 \mathrm{Myr}] with p = p 0 + p 1 T p = p 0 + p 1 T p=p_(0)+p_(1)Tp=p_{0}+p_{1} T, where p 0 p 0 p_(0)p_{0} and p 1 p 1 p_(1)p_{1} are taken from Eq. (35).
在表 7 中，频率项 v k ′ v k ′ v_(k)^(')v_{k}^{\prime} 被识别为形式 p + v k p + v k p+v_(k)p+v_{k} 的组合，其中 p p pp 是进动频率， v k v k v_(k)v_{k} 是轨道平面长期运动的长期频率（来自表 5）。近似值在区间 [ − 15 Myr , + 2 Myr ] [ − 15 Myr , + 2 Myr ] [-15Myr,+2Myr][-15 \mathrm{Myr},+2 \mathrm{Myr}] 和 p = p 0 + p 1 T p = p 0 + p 1 T p=p_(0)+p_(1)Tp=p_{0}+p_{1} T 上进行，其中 p 0 p 0 p_(0)p_{0} 和 p 1 p 1 p_(1)p_{1} 取自方程 (35)。

8.4. Approximation for the precession
8.4. 进动的近似

The full approximate solution for the precession is more complicated because of the presence of the secular resonance s 6 + s 6 + s_(6)+s_{6}+ g 5 − g 6 g 5 − g 6 g_(5)-g_(6)g_{5}-g_{6}. In fact, the proximity of this resonance is responsible for the largest periodic term in the expression of the precession
由于存在长期共振 s 6 + s 6 + s_(6)+s_{6}+ g 5 − g 6 g 5 − g 6 g_(5)-g_(6)g_{5}-g_{6} ，进动的完整近似解更加复杂。事实上，这种共振的接近性导致了进动表达式中最大的周期项
Table 6. Frequency decomposition of the eccentricity of the Earth on the time interval [ − 15 , + 5 ] [ − 15 , + 5 ] [-15,+5][-15,+5] Myr. For all computations, the data of Table 4 should be preferred. Here, e = e 0 + ∑ k = 1 20 b k ′ cos ⁡ ( μ k ′ t + φ k ) e = e 0 + ∑ k = 1 20   b k ′ cos ⁡ μ k ′ t + φ k e=e_(0)+sum_(k=1)^(20)b_(k)^(')cos(mu_(k)^(')t+varphi_(k))e=e_{0}+\sum_{k=1}^{20} b_{k}^{\prime} \cos \left(\mu_{k}^{\prime} t+\varphi_{k}\right) with e 0 = 0.0275579 e 0 = 0.0275579 e_(0)=0.0275579e_{0}=0.0275579. For each term, in the first column is the corresponding combination of frequencies where g i g i g_(i)g_{i} are the fundamental frequencies (Table 3), and μ j μ j mu_(j)\mu_{j} the frequencies of the terms in z z zz from Table 4 .
表 6. 地球偏心率在时间间隔 [ − 15 , + 5 ] [ − 15 , + 5 ] [-15,+5][-15,+5] Myr 内的频率分解。所有计算均应优先使用表 4 中的数据。此处， e = e 0 + ∑ k = 1 20 b k ′ cos ⁡ ( μ k ′ t + φ k ) e = e 0 + ∑ k = 1 20   b k ′ cos ⁡ μ k ′ t + φ k e=e_(0)+sum_(k=1)^(20)b_(k)^(')cos(mu_(k)^(')t+varphi_(k))e=e_{0}+\sum_{k=1}^{20} b_{k}^{\prime} \cos \left(\mu_{k}^{\prime} t+\varphi_{k}\right) 与 e 0 = 0.0275579 e 0 = 0.0275579 e_(0)=0.0275579e_{0}=0.0275579 。对于每一项，第一列是相应的频率组合，其中 g i g i g_(i)g_{i} 为基频（表 3）， μ j μ j mu_(j)\mu_{j} 为表 4 中 z z zz 项的频率。
k k kk μ k ′ ( ′ ′ / yr ) μ k ′ ′ ′ / yr mu_(k)^(')(^('')//yr)\mu_{k}^{\prime}\left({ }^{\prime \prime} / \mathrm{yr}\right) P ( yr ) P ( yr ) P(yr)P(\mathrm{yr}) b k ′ b k ′ b_(k)^(')b_{k}^{\prime} φ k ′ ( φ k ′ ( varphi_(k)^(')(\varphi_{k}^{\prime}( degree ) ) ))   φ k ′ ( φ k ′ ( varphi_(k)^(')(\varphi_{k}^{\prime}( 度 ) ) ))
1 g 2 − g 5 g 2 − g 5 g_(2)-g_(5)g_{2}-g_{5} 3.199279 405091 0.010739 170.739
2 g 4 − g 5 g 4 − g 5 g_(4)-g_(5)g_{4}-g_{5} 13.651920 94932 0.008147 109.891
3 g 4 − g 2 g 4 − g 2 g_(4)-g_(2)g_{4}-g_{2} 10.456224 123945 0.006222 -60.044
4 g 3 − g 5 g 3 − g 5 g_(3)-g_(5)g_{3}-g_{5} 13.109803 98857 0.005287 -86.140
5 g 3 − g 2 g 3 − g 2 g_(3)-g_(2)g_{3}-g_{2} 9.909679 130781 0.004492 100.224
6 μ 7 − μ 6 μ 7 − μ 6 mu_(7)-mu_(6)\mu_{7}-\mu_{6} 0.546076 2373298 0.002967 -168.784
7 g 1 − g 5 g 1 − g 5 g_(1)-g_(5)g_{1}-g_{5} 1.325696 977600 0.002818 57.718
8 g 4 − g 1 g 4 − g 1 g_(4)-g_(1)g_{4}-g_{1} 12.325286 105150 0.002050 49.546
9 g 2 − g 5 + μ 6 − μ 7 g 2 − g 5 + μ 6 − μ 7 g_(2)-g_(5)+mu_(6)-mu_(7)g_{2}-g_{5}+\mu_{6}-\mu_{7} 2.665308 486248 0.001971 148.774
10 g 2 − g 1 g 2 − g 1 g_(2)-g_(1)g_{2}-g_{1} 1.883616 688038 0.001797 137.155
11 μ 6 − g 5 μ 6 − g 5 mu_(6)-g_(5)\mu_{6}-g_{5} 12.856520 100805 0.002074 24.487
12 g 2 + g 4 − 2 g 5 g 2 + g 4 − 2 g 5 g_(2)+g_(4)-2g_(5)g_{2}+g_{4}-2 g_{5} 16.851127 76909 0.001525 102.380
13 μ 6 − g 2 μ 6 − g 2 mu_(6)-g_(2)\mu_{6}-g_{2} 9.650041 134300 0.001491 -167.676
14 2 g 3 − g 4 − g 5 2 g 3 − g 4 − g 5 2g_(3)-g_(4)-g_(5)2 g_{3}-g_{4}-g_{5} 12.563233 103158 0.001316 69.234
15 2 g 4 − g 2 − g 3 2 g 4 − g 2 − g 3 2g_(4)-g_(2)-g_(3)2 g_{4}-g_{2}-g_{3} 10.975914 118077 0.001309 -51.163
16 g 1 − g 3 g 1 − g 3 g_(1)-g_(3)g_{1}-g_{3} 11.788709 109936 0.001300 -146.081
17 μ 7 − g 2 μ 7 − g 2 mu_(7)-g_(2)\mu_{7}-g_{2} 10.194846 127123 0.001306 -139.827
18 g 3 + g 4 − g 2 − g 5 g 3 + g 4 − g 2 − g 5 g_(3)+g_(4)-g_(2)-g_(5)g_{3}+g_{4}-g_{2}-g_{5} 23.562689 55002 0.001261 30.098
19 μ 7 − g 5 μ 7 − g 5 mu_(7)-g_(5)\mu_{7}-g_{5} 13.398519 96727 0.001344 27.612
20 g 2 + g 4 − g 3 − g 5 g 2 + g 4 − g 3 − g 5 g_(2)+g_(4)-g_(3)-g_(5)g_{2}+g_{4}-g_{3}-g_{5} 3.742221 346318 0.001058 178.662
k mu_(k)^(')(^('')//yr) P(yr) b_(k)^(') varphi_(k)^(')( degree ) 1 g_(2)-g_(5) 3.199279 405091 0.010739 170.739 2 g_(4)-g_(5) 13.651920 94932 0.008147 109.891 3 g_(4)-g_(2) 10.456224 123945 0.006222 -60.044 4 g_(3)-g_(5) 13.109803 98857 0.005287 -86.140 5 g_(3)-g_(2) 9.909679 130781 0.004492 100.224 6 mu_(7)-mu_(6) 0.546076 2373298 0.002967 -168.784 7 g_(1)-g_(5) 1.325696 977600 0.002818 57.718 8 g_(4)-g_(1) 12.325286 105150 0.002050 49.546 9 g_(2)-g_(5)+mu_(6)-mu_(7) 2.665308 486248 0.001971 148.774 10 g_(2)-g_(1) 1.883616 688038 0.001797 137.155 11 mu_(6)-g_(5) 12.856520 100805 0.002074 24.487 12 g_(2)+g_(4)-2g_(5) 16.851127 76909 0.001525 102.380 13 mu_(6)-g_(2) 9.650041 134300 0.001491 -167.676 14 2g_(3)-g_(4)-g_(5) 12.563233 103158 0.001316 69.234 15 2g_(4)-g_(2)-g_(3) 10.975914 118077 0.001309 -51.163 16 g_(1)-g_(3) 11.788709 109936 0.001300 -146.081 17 mu_(7)-g_(2) 10.194846 127123 0.001306 -139.827 18 g_(3)+g_(4)-g_(2)-g_(5) 23.562689 55002 0.001261 30.098 19 mu_(7)-g_(5) 13.398519 96727 0.001344 27.612 20 g_(2)+g_(4)-g_(3)-g_(5) 3.742221 346318 0.001058 178.662| $k$ | | $\mu_{k}^{\prime}\left({ }^{\prime \prime} / \mathrm{yr}\right)$ | $P(\mathrm{yr})$ | $b_{k}^{\prime}$ | $\varphi_{k}^{\prime}($ degree $)$ | | ---: | :---: | ---: | ---: | ---: | ---: | | 1 | $g_{2}-g_{5}$ | 3.199279 | 405091 | 0.010739 | 170.739 | | 2 | $g_{4}-g_{5}$ | 13.651920 | 94932 | 0.008147 | 109.891 | | 3 | $g_{4}-g_{2}$ | 10.456224 | 123945 | 0.006222 | -60.044 | | 4 | $g_{3}-g_{5}$ | 13.109803 | 98857 | 0.005287 | -86.140 | | 5 | $g_{3}-g_{2}$ | 9.909679 | 130781 | 0.004492 | 100.224 | | 6 | $\mu_{7}-\mu_{6}$ | 0.546076 | 2373298 | 0.002967 | -168.784 | | 7 | $g_{1}-g_{5}$ | 1.325696 | 977600 | 0.002818 | 57.718 | | 8 | $g_{4}-g_{1}$ | 12.325286 | 105150 | 0.002050 | 49.546 | | 9 | $g_{2}-g_{5}+\mu_{6}-\mu_{7}$ | 2.665308 | 486248 | 0.001971 | 148.774 | | 10 | $g_{2}-g_{1}$ | 1.883616 | 688038 | 0.001797 | 137.155 | | 11 | $\mu_{6}-g_{5}$ | 12.856520 | 100805 | 0.002074 | 24.487 | | 12 | $g_{2}+g_{4}-2 g_{5}$ | 16.851127 | 76909 | 0.001525 | 102.380 | | 13 | $\mu_{6}-g_{2}$ | 9.650041 | 134300 | 0.001491 | -167.676 | | 14 | $2 g_{3}-g_{4}-g_{5}$ | 12.563233 | 103158 | 0.001316 | 69.234 | | 15 | $2 g_{4}-g_{2}-g_{3}$ | 10.975914 | 118077 | 0.001309 | -51.163 | | 16 | $g_{1}-g_{3}$ | 11.788709 | 109936 | 0.001300 | -146.081 | | 17 | $\mu_{7}-g_{2}$ | 10.194846 | 127123 | 0.001306 | -139.827 | | 18 | $g_{3}+g_{4}-g_{2}-g_{5}$ | 23.562689 | 55002 | 0.001261 | 30.098 | | 19 | $\mu_{7}-g_{5}$ | 13.398519 | 96727 | 0.001344 | 27.612 | | 20 | $g_{2}+g_{4}-g_{3}-g_{5}$ | 3.742221 | 346318 | 0.001058 | 178.662 |
angle (Fig. 16). The precession angle can thus be approximate by ( ψ ψ psi\psi in arcsec, t t tt in years)
角度（图 16）。因此，岁差角可以近似为（ ψ ψ psi\psi 角秒， t t tt 年）
ψ = + 49086 + p 0 t + p 1 t 2 + 42246 v 0 2 ( v 0 + 2 p 1 t ) 2 cos ⁡ ( v 0 t + p 1 t 2 + ϕ 0 ) ψ = + 49086 + p 0 t + p 1 t 2 + 42246 v 0 2 v 0 + 2 p 1 t 2 cos ⁡ v 0 t + p 1 t 2 + ϕ 0 {:[psi=+49086+p_(0)t+p_(1)t^(2)],[+42246(v_(0)^(2))/((v_(0)+2p_(1)t)^(2))cos(v_(0)t+p_(1)t^(2)+phi_(0))]:}\begin{aligned} \psi= & +49086+p_{0} t+p_{1} t^{2} \\ & +42246 \frac{v_{0}^{2}}{\left(v_{0}+2 p_{1} t\right)^{2}} \cos \left(v_{0} t+p_{1} t^{2}+\phi_{0}\right) \end{aligned}
where  在哪里
v 0 = 0.150019 ′ ′ / yr v 0 = 0.150019 ′ ′ / yr v_(0)=0.150019^('')//yrv_{0}=0.150019^{\prime \prime} / \mathrm{yr}
p 0 = 50.467718 ′ ′ / yr p 0 = 50.467718 ′ ′ / yr p_(0)=50.467718^('')//yrp_{0}=50.467718^{\prime \prime} / \mathrm{yr}
p 1 = − 13.526564 − 9 ′ ′ / yr 2 p 1 = − 13.526564 − 9 ′ ′ / yr 2 p_(1)=-13.526564^(-9'')//yr^(2)p_{1}=-13.526564^{-9 \prime \prime} / \mathrm{yr}^{2}
ϕ 0 = 171.424 ϕ 0 = 171.424 phi_(0)=171.424\phi_{0}=171.424 degree.   ϕ 0 = 171.424 ϕ 0 = 171.424 phi_(0)=171.424\phi_{0}=171.424 度。
It should be noted that in this formula, the cosine term is an approximation for a more complicated expression that arises from a double integration of the term cos ⁡ ( v 0 t + p 1 t 2 + ϕ 0 ) cos ⁡ v 0 t + p 1 t 2 + ϕ 0 cos(v_(0)t+p_(1)t^(2)+phi_(0))\cos \left(v_{0} t+p_{1} t^{2}+\phi_{0}\right) in the expression of the derivative of the obliquity (6). In order to obtain an expression as simple as possible, without expressions involving Fresnel integrals (Appendix A), we had to restrain its interval of validity to [ − 15 , + 2 ] [ − 15 , + 2 ] [-15,+2][-15,+2] Myr.
需要注意的是，该公式中的余弦项是对更复杂表达式的近似，该表达式源于倾角导数表达式 (6) 中项 cos ⁡ ( v 0 t + p 1 t 2 + ϕ 0 ) cos ⁡ v 0 t + p 1 t 2 + ϕ 0 cos(v_(0)t+p_(1)t^(2)+phi_(0))\cos \left(v_{0} t+p_{1} t^{2}+\phi_{0}\right) 的双重积分。为了获得一个尽可能简单的表达式，不包含涉及菲涅尔积分的表达式（附录 A），我们必须将其有效区间限制为 [ − 15 , + 2 ] [ − 15 , + 2 ] [-15,+2][-15,+2] Myr。
For paleoclimate studies, the usual quantity that relates more directly to insolation is the climatic precession e sin ⁡ ω ¯ e sin ⁡ ω ¯ e sin bar(omega)e \sin \bar{\omega}, where ω ¯ = ψ + ϖ ω ¯ = ψ + ϖ bar(omega)=psi+ϖ\bar{\omega}=\psi+\varpi is the longitude of the perihelion from the moving equinox. We have thus
对于古气候研究，与日照更直接相关的通常量是气候岁差 e sin ⁡ ω ¯ e sin ⁡ ω ¯ e sin bar(omega)e \sin \bar{\omega} ，其中 ω ¯ = ψ + ϖ ω ¯ = ψ + ϖ bar(omega)=psi+ϖ\bar{\omega}=\psi+\varpi 是从移动春分点到近日点的经度。因此，我们有
e sin ⁡ ω ¯ = J ( z e i ψ ( t ) ) e sin ⁡ ω ¯ = J z e i ψ ( t ) e sin bar(omega)=J(ze^(i psi(t)))e \sin \bar{\omega}=\mathfrak{J}\left(z \mathrm{e}^{i \psi(t)}\right),
where the decomposition of z = e exp ⁡ ( i ϖ ) z = e exp ⁡ ( i ϖ ) z=e exp(iϖ)z=e \exp (i \varpi) is provided in Table 4, and the precession angle ψ ψ psi\psi is given by Eq. (34).
其中 z = e exp ⁡ ( i ϖ ) z = e exp ⁡ ( i ϖ ) z=e exp(iϖ)z=e \exp (i \varpi) 的分解见表4，进动角 ψ ψ psi\psi 由公式（34）给出。
Fig. 14. Evolution of the obliquity of the Earth in degrees, from -250 to +250 Myr. The grey zone is the actual obliquity, while the black curve is the averaged value of the obliquity over 0.5 Myr time intervals. The dotted line is a straight line fitted to the average obliquity in the past.
图 14. 地球自转轴倾角的演变（以度为单位），从 -250 百万年到 +250 百万年。灰色区域为实际倾角，黑色曲线为 0.5 百万年时间间隔内的倾角平均值。虚线为拟合过去平均倾角的直线。
The comparison of this approximation with the complete solution from -11 to +1 Myr is given in Fig. 17. The solution for climatic precession is thus on the form
图 17 给出了该近似值与 -11 至 +1 Myr 的完整解的比较。因此，气候岁差的解形式为
e sin ⁡ ω ¯ = ∑ k = 1 20 b k sin ⁡ ( μ k t + φ k + ψ ( t ) ) e sin ⁡ ω ¯ = ∑ k = 1 20   b k sin ⁡ μ k t + φ k + ψ ( t ) e sin bar(omega)=sum_(k=1)^(20)b_(k)sin(mu_(k)t+varphi_(k)+psi(t))e \sin \bar{\omega}=\sum_{k=1}^{20} b_{k} \sin \left(\mu_{k} t+\varphi_{k}+\psi(t)\right),
where b k , μ k , φ k b k , μ k , φ k b_(k),mu_(k),varphi_(k)b_{k}, \mu_{k}, \varphi_{k} are from Table 4 , and where ψ ( t ) ψ ( t ) psi(t)\psi(t) is given by Eq. (34). The frequencies of the terms of the climatic precession are thus of the form
其中 b k , μ k , φ k b k , μ k , φ k b_(k),mu_(k),varphi_(k)b_{k}, \mu_{k}, \varphi_{k} 来自表 4， ψ ( t ) ψ ( t ) psi(t)\psi(t) 由公式 (34) 给出。因此，气候岁差项的频率形式为
μ k ′ ′ = μ k + p + p 1 t μ k ′ ′ = μ k + p + p 1 t mu_(k)^('')=mu_(k)+p+p_(1)t\mu_{k}^{\prime \prime}=\mu_{k}+p+p_{1} t,
where p , p 1 p , p 1 p,p_(1)p, p_{1} are given in Eqs. (34) and (35).
其中 p , p 1 p , p 1 p,p_(1)p, p_{1} 在公式 (34) 和 (35) 中给出。
Table 7. Approximation for the obliquity of the Earth, following Eq. (33). This expression is not strictly quasiperiodic, because of the presence of the dissipative term p 1 p 1 p_(1)p_{1} in the evolution of the precesion frequency (33).
表7. 地球倾角近似值，遵循方程(33)。该表达式并非严格准周期的，因为在进动频率(33)的演化过程中存在耗散项 p 1 p 1 p_(1)p_{1} 。
k k kk v k ′ ( ′ ′ / yr ) v k ′ ′ ′ / yr v_(k)^(')(^('')//yr)v_{k}^{\prime}\left({ }^{\prime \prime} / \mathrm{yr}\right) P ( yr ) P ( yr ) P(yr)P(\mathrm{yr}) a k ′ a k ′ a_(k)^(')a_{k}^{\prime} ϕ k ′ ( d ) ϕ k ′ ( d ) phi_(k)^(')(d)\phi_{k}^{\prime}(\mathrm{d})
1 p + s 3 p + s 3 p+s_(3)p+s_{3} 31.626665 40978 0.582412 86.645
2 p + s 4 p + s 4 p+s_(4)p+s_{4} 32.713667 39616 0.242559 120.859
3 p + s 6 p + s 6 p+s_(6)p+s_{6} 24.124241 53722 0.163685 -35.947
4 p + v 6 p + v 6 p+v_(6)p+v_{6} 32.170778 40285 0.164787 104.689
5 p + v 10 p + v 10 p+v_(10)p+v_{10} 31.081475 41697 0.095382 -112.872
6 p + v 20 p + v 20 p+v_(20)p+v_{20} 31.493347 41152 0.094379 60.778
7 p + s 6 + g 5 − g 6 p + s 6 + g 5 − g 6 p+s_(6)+g_(5)-g_(6)p+s_{6}+g_{5}-g_{6} 0.135393 9572151 0.087136 39.928
8 p + s 2 p + s 2 p+s_(2)p+s_{2} 43.428193 29842 0.064348 -15.130
9 p + s 1 p + s 1 p+s_(1)p+s_{1} 44.865444 28886 0.072451 -155.175
10 31.756641 40810 0.080146 -70.983
11 p + v 13 p + v 13 p+v_(13)p+v_{13} 31.365950 41319 0.072919 10.533
12 32.839446 39465 0.033666 -31.614
13 32.576100 39784 0.033722 77.554
14 32.035200 40455 0.030677 71.757
15 p + v 8 p + v 8 p+v_(8)p+v_{8} 43.537092 29768 0.039351 145.835
16 p + v 7 p + v 7 p+v_(7)p+v_{7} 43.307432 29926 0.030375 160.109
17 p + v 9 p + v 9 p+v_(9)p+v_{9} 43.650496 29690 0.024733 144.926
18 31.903983 40622 0.025201 -173.656
19 30.945195 41880 0.021615 -144.933
20 23.986877 54030 0.021565 -79.670
21 24.257837 53426 0.021270 -178.441
22 32.312463 40108 0.021851 -24.566
23 44.693687 28997 0.014725 124.744
k v_(k)^(')(^('')//yr) P(yr) a_(k)^(') phi_(k)^(')(d) 1 p+s_(3) 31.626665 40978 0.582412 86.645 2 p+s_(4) 32.713667 39616 0.242559 120.859 3 p+s_(6) 24.124241 53722 0.163685 -35.947 4 p+v_(6) 32.170778 40285 0.164787 104.689 5 p+v_(10) 31.081475 41697 0.095382 -112.872 6 p+v_(20) 31.493347 41152 0.094379 60.778 7 p+s_(6)+g_(5)-g_(6) 0.135393 9572151 0.087136 39.928 8 p+s_(2) 43.428193 29842 0.064348 -15.130 9 p+s_(1) 44.865444 28886 0.072451 -155.175 10 31.756641 40810 0.080146 -70.983 11 p+v_(13) 31.365950 41319 0.072919 10.533 12 32.839446 39465 0.033666 -31.614 13 32.576100 39784 0.033722 77.554 14 32.035200 40455 0.030677 71.757 15 p+v_(8) 43.537092 29768 0.039351 145.835 16 p+v_(7) 43.307432 29926 0.030375 160.109 17 p+v_(9) 43.650496 29690 0.024733 144.926 18 31.903983 40622 0.025201 -173.656 19 30.945195 41880 0.021615 -144.933 20 23.986877 54030 0.021565 -79.670 21 24.257837 53426 0.021270 -178.441 22 32.312463 40108 0.021851 -24.566 23 44.693687 28997 0.014725 124.744| $k$ | | $v_{k}^{\prime}\left({ }^{\prime \prime} / \mathrm{yr}\right)$ | $P(\mathrm{yr})$ | $a_{k}^{\prime}$ | $\phi_{k}^{\prime}(\mathrm{d})$ | | :---: | :---: | :---: | :---: | :---: | :---: | | 1 | $p+s_{3}$ | 31.626665 | 40978 | 0.582412 | 86.645 | | 2 | $p+s_{4}$ | 32.713667 | 39616 | 0.242559 | 120.859 | | 3 | $p+s_{6}$ | 24.124241 | 53722 | 0.163685 | -35.947 | | 4 | $p+v_{6}$ | 32.170778 | 40285 | 0.164787 | 104.689 | | 5 | $p+v_{10}$ | 31.081475 | 41697 | 0.095382 | -112.872 | | 6 | $p+v_{20}$ | 31.493347 | 41152 | 0.094379 | 60.778 | | 7 | $p+s_{6}+g_{5}-g_{6}$ | 0.135393 | 9572151 | 0.087136 | 39.928 | | 8 | $p+s_{2}$ | 43.428193 | 29842 | 0.064348 | -15.130 | | 9 | $p+s_{1}$ | 44.865444 | 28886 | 0.072451 | -155.175 | | 10 | | 31.756641 | 40810 | 0.080146 | -70.983 | | 11 | $p+v_{13}$ | 31.365950 | 41319 | 0.072919 | 10.533 | | 12 | | 32.839446 | 39465 | 0.033666 | -31.614 | | 13 | | 32.576100 | 39784 | 0.033722 | 77.554 | | 14 | | 32.035200 | 40455 | 0.030677 | 71.757 | | 15 | $p+v_{8}$ | 43.537092 | 29768 | 0.039351 | 145.835 | | 16 | $p+v_{7}$ | 43.307432 | 29926 | 0.030375 | 160.109 | | 17 | $p+v_{9}$ | 43.650496 | 29690 | 0.024733 | 144.926 | | 18 | | 31.903983 | 40622 | 0.025201 | -173.656 | | 19 | | 30.945195 | 41880 | 0.021615 | -144.933 | | 20 | | 23.986877 | 54030 | 0.021565 | -79.670 | | 21 | | 24.257837 | 53426 | 0.021270 | -178.441 | | 22 | | 32.312463 | 40108 | 0.021851 | -24.566 | | 23 | | 44.693687 | 28997 | 0.014725 | 124.744 |

8.5. Earth-Moon system  8.5. 地月系统

In Fig. 18, the evolution of the Earth-Moon semi-major a M a M a_(M)a_{\mathrm{M}} axis is plotted over 250 Myr in positive and negative time. The grey area corresponds to the short period variations of the EarthMoon distance, obtained in the full integration of the system, while the black curve is the integration of the averaged equations that is used for the precession computations. The agreement between the two solutions is very satisfying.
图18绘制了地月半长轴 a M a M a_(M)a_{\mathrm{M}} 在2.5亿年内的正负时间演化曲线。灰色区域对应于系统完全积分得到的地月距离的短周期变化，而黑色曲线是用于进动计算的平均方程的积分。两个解之间的一致性非常令人满意。
We have searched for a polynomial approximation for the secular evolution of a M a M a_(M)a_{\mathrm{M}}, in order to provide a useful formula for simple analytical computations. In fact, the Earth-Moon evolution is driven by the tidal dissipation equations (Sect. 4.1), that involve the obliquity of the planet. The important singularity of the obliquity in the future (see previous section) prevents from using a single approximation over the whole time interval, and we have used a different polynomial in positive and negative time.
我们寻找了 a M a M a_(M)a_{\mathrm{M}} 长期演化的多项式近似值，以便为简单的解析计算提供有用的公式。事实上，地月演化是由潮汐耗散方程（第4.1节）驱动的，该方程涉及行星的倾角。由于倾角在未来存在重要的奇点（参见上一节），我们无法在整个时间间隔内使用单一近似值，因此我们在正负时间使用了不同的多项式。
a M − = 60.142611 a M + = + 60.142611 − 6.100887 T + 6.120902 T + 1.369407 T 2 − 2.727887 T 2 − 1.484062 T 4 + 1.614481 T 3 − 0.926406 T 4 a M − = 60.142611 a M + = + 60.142611 − 6.100887 T + 6.120902 T + 1.369407 T 2 − 2.727887 T 2 − 1.484062 T 4 + 1.614481 T 3 − 0.926406 T 4 {:[a_(M)^(-)=,60.142611,a_(M)^(+)=],[,+60.142611],[,-6.100887 T,+6.120902 T],[,+1.369407T^(2),-2.727887T^(2)],[,-1.484062T^(4),+1.614481T^(3)],[,-0.926406T^(4)]:}\begin{array}{rlr} a_{\mathrm{M}}^{-}= & 60.142611 & a_{\mathrm{M}}^{+}= \\ & +60.142611 \\ & -6.100887 T & +6.120902 T \\ & +1.369407 T^{2} & -2.727887 T^{2} \\ & -1.484062 T^{4} & +1.614481 T^{3} \\ & -0.926406 T^{4} \end{array}
where T T TT is in billions of years ( Gyr ), and a M + − a M + − a_(M)^(+-)a_{\mathrm{M}}^{+-}in Earth radius (see Table 1). In the residuals plot (Fig. 18b), we can check that this polynomial approximation is very precise for negative time, but in positive time, we have an additional singularity around 150 Myr , resulting from the secular resonance of the precession frequency with v 24 v 24 v_(24)v_{24} (Table 5).
其中 T T TT 以十亿年（Gyr）为单位， a M + − a M + − a_(M)^(+-)a_{\mathrm{M}}^{+-} 以地球半径为单位（见表 1）。在残差图（图 18b）中，我们可以验证，这个多项式近似对于负时间非常精确，但在正时间中，我们在 1.5 亿年左右发现了一个额外的奇点，这是由于进动频率与 v 24 v 24 v_(24)v_{24} 的长期共振造成的（表 5）。
In the same way, it is possible to approximate the evolution of the precession frequency of the Earth p p pp as the Moon goes away and the rotation of the Earth slows down with time. The polynomial approximation of the precession frequency is then in arcsec yr − 1 arcsec yr − 1 arcsecyr^(-1)\operatorname{arcsec} \mathrm{yr}^{-1}, with T T TT in Gyr
同理，随着月球逐渐远离地球，地球自转速度逐渐减慢，我们可以近似计算地球进动频率 p p pp 的演变。进动频率的多项式近似值为 arcsec yr − 1 arcsec yr − 1 arcsecyr^(-1)\operatorname{arcsec} \mathrm{yr}^{-1} ，其中 T T TT 为 Gyr
p − = 50.475838 p + = 50.475838 − 26.368583 T + 21.890862 T 2 − 27.000654 T + 15.603265 T 2 p − = 50.475838 p + = 50.475838 − 26.368583 T + 21.890862 T 2 − 27.000654 T + 15.603265 T 2 {:[p^(-)=,50.475838,p^(+)=,50.475838],[,-26.368583 T],[,+21.890862T^(2)quad,-27.000654 T],[,+15.603265T^(2)]:}\begin{array}{rlrl} p^{-}= & 50.475838 & p^{+}= & 50.475838 \\ & -26.368583 T \\ & +21.890862 T^{2} \quad & -27.000654 T \\ & +15.603265 T^{2} \end{array}
It should be said that the initial value is kept equal for both positive and negative time to avoid discontinuities. It is clear from Fig. 19 that for the period between +30 and +130 Myr , a better approximation could be obtained by adding to the constant part the offset +0.135052 arcsec yr − 1 yr − 1 yr^(-1)\mathrm{yr}^{-1}. It is interesting to note again the dramatic influence of the resonance v 23 v 23 v_(23)v_{23} on the
应该指出的是，为了避免不连续性，初始值在正负时间上都保持相等。从图19可以清楚地看出，在+30至+130百万年之间的时间段内，通过在常数部分添加偏移量+0.135052角秒 yr − 1 yr − 1 yr^(-1)\mathrm{yr}^{-1} 可以获得更好的近似值。值得再次注意的是，共振 v 23 v 23 v_(23)v_{23} 对
Fig. 15. Comparison of the solution of the obliquity La2004 (solid line) with its approximation using Table 7 (dotted line) with 26 periodic terms. The difference of two solutions (+22 degrees) is also plotted.
图 15. 倾角 La2004 的解（实线）与使用表 7 的近似值（虚线）（包含 26 个周期项）的比较。图中还绘制了两个解之间的差异（+22 度）。
evolution of the precession frequency. The equivalent formulas for the length of the day are (in hours)
岁差频率的演变。一天的长度的等效公式是（以小时为单位）
L O D − = 23.934468 L O D + = 23.934468 + 7.432167 T + 7.444649 T − 0.727046 T 2 − 0.715049 T 2 . + 0.409572 T 3 + 0.458097 T 3 − 0.589692 T 4 L O D − = 23.934468 L O D + = 23.934468 + 7.432167 T + 7.444649 T − 0.727046 T 2 − 0.715049 T 2 . + 0.409572 T 3 + 0.458097 T 3 − 0.589692 T 4 {:[LOD^(-)=,23.934468,LOD^(+)=,23.934468],[,+7.432167 T,,+7.444649 T],[,-0.727046T^(2),,-0.715049T^(2).],[,+0.409572T^(3),,+0.458097T^(3)],[,-0.589692T^(4),,]:}\begin{array}{rlrl} L O D^{-}= & 23.934468 & L O D^{+}= & 23.934468 \\ & +7.432167 T & & +7.444649 T \\ & -0.727046 T^{2} & & -0.715049 T^{2} . \\ & +0.409572 T^{3} & & +0.458097 T^{3} \\ & -0.589692 T^{4} & & \end{array}
One should note that Eqs. (39) and (41) correspond at the origin J2000 to a change in the semi major axis of the Moon of about 3.89 cm / yr 3.89 cm / yr 3.89cm//yr3.89 \mathrm{~cm} / \mathrm{yr}, and a change of the LOD of 2.68 ms / 2.68 ms / 2.68ms//2.68 \mathrm{~ms} / century. These values slightly differ from the current ones of Dickey et al. (1994) ( 3.82 ± 0.07 cm / yr 3.82 ± 0.07 cm / yr 3.82+-0.07cm//yr3.82 \pm 0.07 \mathrm{~cm} / \mathrm{yr} for the receding of the Moon). This is understandable, as the Earth-Moon model in our integration that was adjusted to DE406 is simplified with respect to the one of DE406.
需要注意的是，方程 (39) 和 (41) 在原点 J2000 处对应于月球半长轴约 3.89 cm / yr 3.89 cm / yr 3.89cm//yr3.89 \mathrm{~cm} / \mathrm{yr} 的变化，以及 LOD 2.68 ms / 2.68 ms / 2.68ms//2.68 \mathrm{~ms} / 个世纪的变化。这些值与 Dickey 等人 (1994) 的当前值（月球退行距离为 3.82 ± 0.07 cm / yr 3.82 ± 0.07 cm / yr 3.82+-0.07cm//yr3.82 \pm 0.07 \mathrm{~cm} / \mathrm{yr} ）略有不同。这是可以理解的，因为我们积分中针对 DE406 进行调整的地月模型，相对于 DE406 的模型进行了简化。

9. Stability of the solution
9. 溶液的稳定性

In such long term numerical computations, the sources of errors are the error of method due to the limitation of the numerical integration algorithm, the numerical roundoff error resulting from the limitation of the representation of the real numbers
在这种长期数值计算中，误差来源是由于数值积分算法的限制而产生的方法误差，由于实数表示的限制而产生的数值舍入误差
Fig. 16. Comparison of the solution of the precession La2004 with its approximation from -15 to +2 Myr . The grey curve is obtained after removing uniquely the secular trend ψ 0 + p 0 t + p 1 t 2 ψ 0 + p 0 t + p 1 t 2 psi_(0)+p_(0)t+p_(1)t^(2)\psi_{0}+p_{0} t+p_{1} t^{2}. The dark curve are the residual (in radians) after removing the resonant term (Eq. (34)), which is also displayed in solid line.
图 16. La2004 进动的解与其在-15 百万年至+2 百万年期间的近似值的比较。灰色曲线是在唯一去除长期趋势 ψ 0 + p 0 t + p 1 t 2 ψ 0 + p 0 t + p 1 t 2 psi_(0)+p_(0)t+p_(1)t^(2)\psi_{0}+p_{0} t+p_{1} t^{2} 后获得的。深色曲线是去除共振项（方程(34)）后的残差（以弧度为单位），该残差也以实线显示。
Fig. 17. Comparison of the solution of the climatic precession of La 2004 with its approximation from -11 to +1 Myr . The grey curve is the full climatic precession e sin ⁡ ω ¯ e sin ⁡ ω ¯ e sin bar(omega)e \sin \bar{\omega}. The dark curve are the residual after removing the approximation given by Eq. (37).
图 17. La 2004 气候岁差解与其-11 至+1 百万年近似值的比较。灰色曲线表示完整的气候岁差 e sin ⁡ ω ¯ e sin ⁡ ω ¯ e sin bar(omega)e \sin \bar{\omega} 。深色曲线表示去除方程(37)给出的近似值后的残差。
in the computer, and the errors of the initial conditions and of the model which are due to our imperfect knowledge of all the physical parameters in the Solar system, and to the necessary limitations of the model. Some of the main sources of uncertainty in the model for long term integrations were reviewed in (Laskar 1999). The effect of all these errors is amplified by the chaotic behaviour of the system, with an exponential increase of the difference between two solutions with different settings, until saturation due to the limited range of the considered variables. A summary of these limitations, extracted from (Laskar 1999) is given in Table 8, but it should be noted that these times of validity are probably pessimistic, as they correspond to analytical estimates based on the “worst case” situation. Our requirement for the precision of a long term solution is also not the same as for precise short term ephemeris. For the long term solutions, aimed at paleoclimate or qualitative studies, two solutions can be considered as similar, as long as the pattern of the two solutions (eccentricity for example), resulting (at first order) from the combination of various proper modes remain similar. This will last until some of the main proper modes get completely out of phase (see Sect. 7). It will be the same for the solutions in obliquity and precession.
在计算机中，初始条件和模型的误差是由于我们对太阳系所有物理参数的了解不够完善，以及模型的必要限制。长期积分模型中一些主要的不确定性来源已在（Laskar 1999）中进行了综述。所有这些误差的影响都被系统的混沌行为放大，不同设置的两个解之间的差异呈指数增长，直至因所考虑变量的范围有限而达到饱和。表 8 给出了摘自（Laskar 1999）的这些限制的摘要，但应该注意，这些有效时间可能是悲观的，因为它们对应于基于“最坏情况”的分析估计。我们对长期解决方案的精度要求与精确短期星历表不同。对于针对古气候或定性研究的长期解，只要由各种固有模态组合（在一阶上）得到的两个解（例如偏心率）的模式保持相似，就可以认为这两个解相似。这种情况将持续到某些主要的固有模态完全失相（参见第 7 节）。对于行星倾角和岁差的解，情况也是如此。
Fig. 18. a) Evolution of the Earth-Moon semi major axis (in Earth radii) from -250 to +250 Myr . The grey zone is the result of the integration of the full equations, while the black curve is the integration of the averaged equations, as used in the precession computations; b) residuals with the fit of the averaged solution with the polynomial expression (39).
图 18. a) 地月半长轴（以地球半径为单位）从-250 百万年到+250 百万年的演变。灰色区域是完整方程积分的结果，黑色曲线是平均方程积分的结果，用于进动计算；b) 平均解与多项式表达式（39）拟合的残差。
Fig. 19. a) Evolution of the precession frequency of the Earth p p pp (in arcsec / yr arcsec / yr arcsec//yr\operatorname{arcsec} / \mathrm{yr} ) -250 to +250 Myr ; b) residuals with the fit of the averaged solution with the polynomial expression (40).
图 19. a) 地球进动频率的演变 p p pp （在 arcsec / yr arcsec / yr arcsec//yr\operatorname{arcsec} / \mathrm{yr} 中）-250 至 +250 Myr；b) 与多项式表达式 (40) 的平均解拟合的残差。

9.1. Orbital motion  9.1 轨道运动

Practically, we have tested the stability of our solution by comparison of the eccentricity (Fig. 20) and inclination (Fig. 21) with different solutions. In Figs. 20a and 21a, the nominal solution La2004 (with stepsize τ = 5 × 10 − 3 τ = 5 × 10 − 3 tau=5xx10^(-3)\tau=5 \times 10^{-3} years) is compared to an alternate solution, La2004*, with the same dynamical
实际上，我们通过比较不同解的偏心率（图 20）和倾角（图 21）来测试我们解的稳定性。在图 20a 和 21a 中，我们将标称解 La2004（步长为 τ = 5 × 10 − 3 τ = 5 × 10 − 3 tau=5xx10^(-3)\tau=5 \times 10^{-3} 年）与具有相同动力学参数的替代解 La2004* 进行了比较。
Table 8. Main sources of uncertainty in the orbital solution (from Laskar 1999). For each limiting factor, an analytical estimate of the time of validity of the solution T V T V T_(V)T_{\mathrm{V}} is given (in Myr), taking into account the exponential growth of the error.
表 8. 轨道解的主要不确定性来源（源自 Laskar 1999）。对于每个限制因素，给出了解的有效时间 T V T V T_(V)T_{\mathrm{V}} 的分析估计（以百万年为单位），同时考虑到误差的指数增长。
Limiting factor  限制因素 T V T V T_(V)T_{\mathrm{V}}
Uncertainty on the masses and intial conditions
质量和初始条件的不确定性		 38 Myr  38 马来西亚林吉特
Contribution of the main Galilean satellites
主要伽利略卫星的贡献		 35 Myr  35 马来西亚林吉特
Uncertainty in the Earth-Moon system evolution
地月系统演化的不确定性		 40 Myr  40 马来西亚林吉特
Effect of the main asteroids
主要小行星的影响		 32 Myr  32 马来西亚林吉特
Mass loss of the Sun
太阳质量损失		 50 Myr  50 马来西亚林吉特
Uncertainty of 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} on the J 2 J 2 J_(2)J_{2} of the Sun
太阳 J 2 J 2 J_(2)J_{2} 的不确定性为 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} 		26 Myr  26 马来西亚林吉特
Limiting factor T_(V) Uncertainty on the masses and intial conditions 38 Myr Contribution of the main Galilean satellites 35 Myr Uncertainty in the Earth-Moon system evolution 40 Myr Effect of the main asteroids 32 Myr Mass loss of the Sun 50 Myr Uncertainty of 2xx10^(-7) on the J_(2) of the Sun 26 Myr| Limiting factor | $T_{\mathrm{V}}$ | | :--- | ---: | | Uncertainty on the masses and intial conditions | 38 Myr | | Contribution of the main Galilean satellites | 35 Myr | | Uncertainty in the Earth-Moon system evolution | 40 Myr | | Effect of the main asteroids | 32 Myr | | Mass loss of the Sun | 50 Myr | | Uncertainty of $2 \times 10^{-7}$ on the $J_{2}$ of the Sun | 26 Myr |

Fig. 20. Stability of the solution for eccentricity of the Earth. a) Difference of the nominal solution La2004 with stepsize τ = τ = tau=\tau= 5. × 10 − 3 × 10 − 3 xx10^(-3)\times 10^{-3} years, and La2004*, obtained with τ ∗ = 4.8828125 × 10 − 3 τ ∗ = 4.8828125 × 10 − 3 tau^(**)=4.8828125 xx10^(-3)\tau^{*}=4.8828125 \times 10^{-3} years; b) difference of the nominal solution with the solution obtained while setting J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 for the S S SS un (instead of 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} in the nominal solution).
图 20. 地球偏心率解的稳定性。a) 步长为 τ = τ = tau=\tau= 5. × 10 − 3 × 10 − 3 xx10^(-3)\times 10^{-3} 年的标称解 La2004 与步长为 τ ∗ = 4.8828125 × 10 − 3 τ ∗ = 4.8828125 × 10 − 3 tau^(**)=4.8828125 xx10^(-3)\tau^{*}=4.8828125 \times 10^{-3} 年的 La2004* 之间的差异；b) 标称解与将 S S SS un 设置 J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 （而不是标称解中的 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} ）时获得的解之间的差异。
model, and a very close stepsize τ ∗ = 4.8828125 × 10 − 3 τ ∗ = 4.8828125 × 10 − 3 tau^(**)=4.8828125 xx10^(-3)\tau^{*}=4.8828125 \times 10^{-3} years (Table 9). This special value was chosen in order that our output time span h = 1000 h = 1000 h=1000h=1000 years corresponds to an integer number (204 800) of steps, in order to avoid any interpolation problems in the check of the numerical accuracy. This is thus a test of the time of validity for obtaining a precise numerical solution, resulting from method and roundoff errors of the integrator. It is thus limited here to about 60 Myr .
模型，步长非常接近 τ ∗ = 4.8828125 × 10 − 3 τ ∗ = 4.8828125 × 10 − 3 tau^(**)=4.8828125 xx10^(-3)\tau^{*}=4.8828125 \times 10^{-3} 年（表 9）。选择这个特殊值是为了使我们的输出时间跨度 h = 1000 h = 1000 h=1000h=1000 年对应于整数步长（204 800），以避免在检查数值精度时出现任何插值问题。因此，这是对获得精确数值解的有效性时间的测试，该时间是由积分器的方法和舍入误差引起的。因此，这里将其限制在约 60 Myr。
This limitation of 60 Myr is a limitation for the time of validity of an orbital solution, independently of the precision of the dynamical model. In order to go beyond this limit, the only way will be to increase the numerical accuracy of our computations, by improving the numerical algorithm, or with an extended precision for the number representation in the computer.
6000万年的这个限制是轨道解有效时间的限制，与动力学模型的精度无关。为了突破这个限制，唯一的方法是提高计算的数值精度，例如改进数值算法，或者提高计算机中数值表示的精度。
A second test is made on the uncertainty of the model. Several sources of uncertainty were listed in (Laskar 1999), and it was found that one of the main sources of uncertainty was due to the imprecise knowledge of the Solar
对模型的不确定性进行了第二次检验。Laskar 1999 列出了几种不确定性来源，发现其中一个主要不确定性来源是由于对太阳活动规律的了解不够精确。

Fig. 21. Stability of the solution for inclination of the Earth. a) Difference of sin ⁡ i / 2 sin ⁡ i / 2 sin i//2\sin i / 2 for the nominal solution La2004 with stepsize τ = 5 . × 10 − 3 τ = 5 . × 10 − 3 tau=5.xx10^(-3)\tau=5 . \times 10^{-3} years, and for La2004*, obtained with τ ∗ = τ ∗ = tau^(**)=\tau^{*}= 4.8828125 × 10 − 3 4.8828125 × 10 − 3 4.8828125 xx10^(-3)4.8828125 \times 10^{-3} years; b) difference of sin ⁡ i / 2 sin ⁡ i / 2 sin i//2\sin i / 2 for the nominal solution with the solution obtained while setting J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 for the Sun (instead of 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} in the nominal solution).
图 21. 地球倾斜度解的稳定性。a) 步长为 τ = 5 . × 10 − 3 τ = 5 . × 10 − 3 tau=5.xx10^(-3)\tau=5 . \times 10^{-3} 年的标称解 La2004 的差异为 sin ⁡ i / 2 sin ⁡ i / 2 sin i//2\sin i / 2 ，而步长为 τ ∗ = τ ∗ = tau^(**)=\tau^{*}= 4.8828125 × 10 − 3 4.8828125 × 10 − 3 4.8828125 xx10^(-3)4.8828125 \times 10^{-3} 年的 La2004* 的差异为 sin ⁡ i / 2 sin ⁡ i / 2 sin i//2\sin i / 2 ；b) 标称解与将太阳设置为 J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 （而不是标称解中的 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} ）时获得的解的差异为 sin ⁡ i / 2 sin ⁡ i / 2 sin i//2\sin i / 2 。
Table 9. Variants of the La2004 solutions. The nominal solution La2004 is obtained for a stepsize τ = 5 × 10 − 3 τ = 5 × 10 − 3 tau=5xx10^(-3)\tau=5 \times 10^{-3}, and a Solar oblateness J 2 S = 2.0 × 10 − 7 J 2 S = 2.0 × 10 − 7 J_(2)^(S)=2.0 xx10^(-7)J_{2}^{\mathrm{S}}=2.0 \times 10^{-7}. The other solutions, with different setting have been computed to test the stability of the nominal solution. In the last column γ γ gamma\gamma is the tidal dissipation factor.
表 9. La2004 解的变体。标称解 La2004 是在步长为 τ = 5 × 10 − 3 τ = 5 × 10 − 3 tau=5xx10^(-3)\tau=5 \times 10^{-3} 、太阳扁率为 J 2 S = 2.0 × 10 − 7 J 2 S = 2.0 × 10 − 7 J_(2)^(S)=2.0 xx10^(-7)J_{2}^{\mathrm{S}}=2.0 \times 10^{-7} 的情况下获得的。为了测试标称解的稳定性，我们计算了其他不同设置的解。最后一列 γ γ gamma\gamma 是潮汐耗散因子。
Name  姓名 τ τ tau\tau (yr)   τ τ tau\tau （年份） J 2 S ― J 2 S ¯ bar(J_(2)^(S))\overline{J_{2}^{\mathrm{S}}} γ γ gamma\gamma
La2004 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 2.0 × 10 − 7 2.0 × 10 − 7 2.0 xx10^(-7)2.0 \times 10^{-7} 1
La2004*  2004年* 4.8828125 × 10 − 3 4.8828125 × 10 − 3 4.8828125 xx10^(-3)4.8828125 \times 10^{-3} 2.0 × 10 − 7 2.0 × 10 − 7 2.0 xx10^(-7)2.0 \times 10^{-7} 1
La2004**  2004** 5.00000005 × 10 − 3 5.00000005 × 10 − 3 5.00000005 xx10^(-3)5.00000005 \times 10^{-3} 2.0 × 10 − 7 2.0 × 10 − 7 2.0 xx10^(-7)2.0 \times 10^{-7} 1
La2004 0 0 _(0){ }_{0} 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 0 1
La2004 1.0 1.0 _(1.0){ }_{1.0} 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 1.0 × 10 − 7 1.0 × 10 − 7 1.0 xx10^(-7)1.0 \times 10^{-7} 1
La2004 1.5 1.5 _(1.5){ }_{1.5} 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 1.5 × 10 − 7 1.5 × 10 − 7 1.5 xx10^(-7)1.5 \times 10^{-7} 1
La2004 0.95 0.95 ^(0.95){ }^{0.95} 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 2.0 × 10 − 7 2.0 × 10 − 7 2.0 xx10^(-7)2.0 \times 10^{-7} 0.95
La2004 0.90 0.90 ^(0.90){ }^{0.90} 5 × 10 − 3 5 × 10 − 3 5xx10^(-3)5 \times 10^{-3} 2.0 × 10 − 7 2.0 × 10 − 7 2.0 xx10^(-7)2.0 \times 10^{-7} 0.90
Name tau (yr) bar(J_(2)^(S)) gamma La2004 5xx10^(-3) 2.0 xx10^(-7) 1 La2004* 4.8828125 xx10^(-3) 2.0 xx10^(-7) 1 La2004** 5.00000005 xx10^(-3) 2.0 xx10^(-7) 1 La2004 _(0) 5xx10^(-3) 0 1 La2004 _(1.0) 5xx10^(-3) 1.0 xx10^(-7) 1 La2004 _(1.5) 5xx10^(-3) 1.5 xx10^(-7) 1 La2004 ^(0.95) 5xx10^(-3) 2.0 xx10^(-7) 0.95 La2004 ^(0.90) 5xx10^(-3) 2.0 xx10^(-7) 0.90| Name | $\tau$ (yr) | $\overline{J_{2}^{\mathrm{S}}}$ | $\gamma$ | | :---: | :---: | :---: | :---: | | La2004 | $5 \times 10^{-3}$ | $2.0 \times 10^{-7}$ | 1 | | La2004* | $4.8828125 \times 10^{-3}$ | $2.0 \times 10^{-7}$ | 1 | | La2004** | $5.00000005 \times 10^{-3}$ | $2.0 \times 10^{-7}$ | 1 | | La2004 ${ }_{0}$ | $5 \times 10^{-3}$ | 0 | 1 | | La2004 ${ }_{1.0}$ | $5 \times 10^{-3}$ | $1.0 \times 10^{-7}$ | 1 | | La2004 ${ }_{1.5}$ | $5 \times 10^{-3}$ | $1.5 \times 10^{-7}$ | 1 | | La2004 ${ }^{0.95}$ | $5 \times 10^{-3}$ | $2.0 \times 10^{-7}$ | 0.95 | | La2004 ${ }^{0.90}$ | $5 \times 10^{-3}$ | $2.0 \times 10^{-7}$ | 0.90 |
oblateness ( J 2 S ) J 2 S (J_(2)^(S))\left(J_{2}^{\mathrm{S}}\right) (Table 8). Its most recent determination, obtained with the SOHO and GONG helioseismic data give J 2 S = J 2 S = J_(2)^(S)=J_{2}^{\mathrm{S}}= ( 2.18 ± 0.06 ) × 10 − 7 ( 2.18 ± 0.06 ) × 10 − 7 (2.18+-0.06)xx10^(-7)(2.18 \pm 0.06) \times 10^{-7} (Pijpers 1998), with a very similar value adopted in DE406 ( J 2 S = 2 × 10 − 7 J 2 S = 2 × 10 − 7 J_(2)^(S)=2xx10^(-7)J_{2}^{S}=2 \times 10^{-7}, Standish 1998), while in DE200, it was not taken into account (Newhall et al. 1983). We will thus consider that comparing the nominal solution La2004 (with J 2 S = 2 × 10 − 7 J 2 S = 2 × 10 − 7 J_(2)^(S)=2xx10^(-7)J_{2}^{\mathrm{S}}=2 \times 10^{-7} ) with an alternate solution (La2004 0 0 _(0){ }_{0} ) with J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 is representative of the uncertainty of the dynamical model for our long term integrations 1 1 ^(1){ }^{1} (Table 9).
扁率 ( J 2 S ) J 2 S (J_(2)^(S))\left(J_{2}^{\mathrm{S}}\right) （表 8）。其最新测定结果，利用 SOHO 和 GONG 日震数据得出 J 2 S = J 2 S = J_(2)^(S)=J_{2}^{\mathrm{S}}= ( 2.18 ± 0.06 ) × 10 − 7 ( 2.18 ± 0.06 ) × 10 − 7 (2.18+-0.06)xx10^(-7)(2.18 \pm 0.06) \times 10^{-7} （Pijpers 1998），与 DE406（ J 2 S = 2 × 10 − 7 J 2 S = 2 × 10 − 7 J_(2)^(S)=2xx10^(-7)J_{2}^{S}=2 \times 10^{-7} ，Standish 1998）采用的值非常相似，而在 DE200 中，并未将其考虑在内（Newhall 等人，1983）。因此，我们认为，将标准解 La2004（ J 2 S = 2 × 10 − 7 J 2 S = 2 × 10 − 7 J_(2)^(S)=2xx10^(-7)J_{2}^{\mathrm{S}}=2 \times 10^{-7} ）与备选解（La2004 0 0 _(0){ }_{0} ） J 2 S = 0 J 2 S = 0 J_(2)^(S)=0J_{2}^{\mathrm{S}}=0 进行比较，可以代表我们长期积分动力学模型的不确定性 1 1 ^(1){ }^{1} （表 9）。
The results La2004 - La2004 for the Earth’s eccentricity and inclination are displayed in (Figs. 20b and 21b)
地球偏心率和倾角的结果 La2004 - La2004 显示在（图 20b 和 21b）中
over 100 Myr. As for Mars (Laskar et al. 2004), the effect of J 2 S J 2 S J_(2)^(S)J_{2}^{\mathrm{S}} becomes noticeable after about 30 Myr ( 26 Myr were predicted with an analytical estimate in (Laskar 1999), and the solution remains very similar over 40 Myr , and totally out of phase after 45 Myr . We will thus consider here that 40 Myr is about the time of validity of our present orbital solution for the Earth. Other integrations made with J 2 S = 1.5 × 10 − 7 J 2 S = 1.5 × 10 − 7 J_(2)^(S)=1.5 xx10^(-7)J_{2}^{\mathrm{S}}=1.5 \times 10^{-7} and J 2 S = 1.0 × 10 − 7 J 2 S = 1.0 × 10 − 7 J_(2)^(S)=1.0 xx10^(-7)J_{2}^{\mathrm{S}}=1.0 \times 10^{-7} (Table 9) gave practically the same estimates. It should be noted that with our present algorithms, we are much more limited by the precision of the model (Figs. 20b and 21b) than by the numerical accuracy (Figs. 20a and 21a).
超过 1 亿年。对于火星 (Laskar 等人 2004)， J 2 S J 2 S J_(2)^(S)J_{2}^{\mathrm{S}} 的影响在大约 30 百万年后变得明显（Laskar 1999 中的分析估计预测了 26 百万年），并且该解决方案在 40 百万年后仍然非常相似，在 45 百万年之后完全不同步。因此，我们在这里认为 40 百万年大约是我们目前的地球轨道解决方案的有效时间。使用 J 2 S = 1.5 × 10 − 7 J 2 S = 1.5 × 10 − 7 J_(2)^(S)=1.5 xx10^(-7)J_{2}^{\mathrm{S}}=1.5 \times 10^{-7} 和 J 2 S = 1.0 × 10 − 7 J 2 S = 1.0 × 10 − 7 J_(2)^(S)=1.0 xx10^(-7)J_{2}^{\mathrm{S}}=1.0 \times 10^{-7} 进行的其他积分（表 9）给出了几乎相同的估计值。应该注意，使用我们目前的算法，我们受模型精度（图 20b 和 21b）的限制比受数值精度（图 20a 和 21a）的限制更大。

9.2. Obliquity and precession
9.2 倾角和进动

For the obliquity and precession, as was already mentioned in (Laskar et al. 1993), the main source of uncertainty is the evolution of the tidal dissipation and other related factors (see Sect. 4.1). Among these factors, one can mention the variations of the dynamical ellipticity of the Earth resulting from changes of the ice caps, or mantle convections. The maximal contribution of these effects can be estimated (see Levrard & Laskar 2003), but it is very difficult to obtain at present a precise value of their contributions. Moreover, the largest dissipative term, resulting from body tides can evolve substantially as a result of the variation of the distribution of the continents, as most of the tidal dissipation occurs in shallow seas (Jeffreys 1920, 1976; Egbert & Ray 2000).
对于地轴倾角和岁差，正如 Laskar et al. 1993 中所述，主要的不确定性来源是潮汐耗散的演变和其他相关因素（见 4.1 节）。在这些因素中，值得一提的是冰盖变化或地幔对流变化导致的地球动力学椭圆度的变化。这些影响的最大贡献可以估算（见 Levrard & Laskar 2003），但目前很难获得其贡献的精确值。此外，由于大部分潮汐耗散发生在浅海中，由体潮汐引起的最大耗散项可能会因大陆分布的变化而发生显著变化（Jeffreys 1920, 1976; Egbert & Ray 2000）。
In Fig. 22, we have explicited the result of a difference of 5% (Figs. 22a and c) and 10 % 10 % 10%10 \% (Figs. 22b and d) in the tidal dissipation term of the Earth and the Moon ( Δ t Δ t Delta t\Delta t and Δ t M Δ t M Deltat_(M)\Delta t_{\mathrm{M}} in Table 1). With 5 % 5 % 5%5 \% error ( La 2004 ( 0.95 ) ) La 2004 ( 0.95 ) (La2004^((0.95)))\left(\operatorname{La} 2004^{(0.95)}\right), the solution of obliquity is valid over about 20 Myr , but with a 10 % 10 % 10%10 \% error ( La 2004 ( 0.90 ) ) La 2004 ( 0.90 ) (La2004^((0.90)))\left(\mathrm{La} 2004^{(0.90)}\right), the solution is out of phase after 20 Myr. Nevertheless, the uncertainty in the tidal dissipation appears mostly as a small change of the precession frequency p p pp, which reveals in the obliquity solution as a time offset that does not change much the obliquity pattern. Indeed, in Figs. 22e and f, we have made a readjustment of the time scale for the solutions La2004 ( 0.95 ) ( 0.95 ) ^((0.95)){ }^{(0.95)} and La2004 ( 0.90 ) ( 0.90 ) ^((0.90)){ }^{(0.90)} of respectively
在图 22 中，我们明确了地球和月球潮汐耗散项（表 1 中的 Δ t Δ t Delta t\Delta t 和 Δ t M Δ t M Deltat_(M)\Delta t_{\mathrm{M}} ）存在 5%（图 22a 和 c）和 10 % 10 % 10%10 \% （图 22b 和 d）差异的结果。当存在 5 % 5 % 5%5 \% 误差 ( La 2004 ( 0.95 ) ) La 2004 ( 0.95 ) (La2004^((0.95)))\left(\operatorname{La} 2004^{(0.95)}\right) 时，倾角解在约 20 Myr 内有效，但当存在 10 % 10 % 10%10 \% 误差 ( La 2004 ( 0.90 ) ) La 2004 ( 0.90 ) (La2004^((0.90)))\left(\mathrm{La} 2004^{(0.90)}\right) 时，倾角解在 20 Myr 之后出现相位差。尽管如此，潮汐耗散的不确定性主要表现为进动频率 p p pp 的微小变化，这在倾角解中表现为时间偏移，但不会显著改变倾角模式。实际上，在图 22e 和 f 中，我们分别重新调整了 La2004 ( 0.95 ) ( 0.95 ) ^((0.95)){ }^{(0.95)} 和 La2004 ( 0.90 ) ( 0.90 ) ^((0.90)){ }^{(0.90)} 解的时间尺度
δ t ( 0.95 ) = 1.3346 × 10 − 5 t − 1.9362 × 10 − 5 t 2 δ t ( 0.95 ) = 1.3346 × 10 − 5 t − 1.9362 × 10 − 5 t 2 deltat_((0.95))=1.3346 xx10^(-5)t-1.9362 xx10^(-5)t^(2)\delta t_{(0.95)}=1.3346 \times 10^{-5} t-1.9362 \times 10^{-5} t^{2}
δ t ( 0.90 ) = 1.4358 × 10 − 5 t − 3.9666 × 10 − 5 t 2 δ t ( 0.90 ) = 1.4358 × 10 − 5 t − 3.9666 × 10 − 5 t 2 deltat_((0.90))=1.4358 xx10^(-5)t-3.9666 xx10^(-5)t^(2)\delta t_{(0.90)}=1.4358 \times 10^{-5} t-3.9666 \times 10^{-5} t^{2}
where t t tt is in Myr. With this new time scales, the different solutions become very close over 40 Myr , that is over the full range of validity of the orbital solution (Figs. 22e and f). An uncertainty of 10 % 10 % 10%10 \% in the tidal dissipation term thus corresponds mostly to an uncertainty in the timescale of 16 kyr after 20 Myr , and 63 kyr after 40 Myr .
其中 t t tt 以百万年为单位。采用新的时间尺度，不同的解在 40 百万年内变得非常接近，也就是在轨道解的整个有效范围内（图 22e 和 f）。因此，潮汐耗散项中 10 % 10 % 10%10 \% 的不确定性主要对应于 20 百万年之后 16 千年以及 40 百万年之后 63 千年的时间尺度上的不确定性。

10. Resonant angles  10. 共振角

10.1. Secular resonances  10.1. 长期共鸣

Using the secular equations, Laskar (1990, 1992), demonstrated that the chaotic behavior of the Solar system arises from multiple secular resonances in the inner solar system. In particular, the critical argument associated to
Laskar (1990, 1992) 利用长期方程证明了太阳系的混沌行为源于太阳系内部的多个长期共振。特别是，与
θ = ( s 4 − s 3 ) − 2 ( g 4 − g 3 ) θ = s 4 − s 3 − 2 g 4 − g 3 theta=(s_(4)-s_(3))-2(g_(4)-g_(3))\theta=\left(s_{4}-s_{3}\right)-2\left(g_{4}-g_{3}\right)
Fig. 22. Difference in obliquity for different tidal dissipation factor (Table 9): a) La2004-La2004 ( 0.95 ) ( 0.95 ) ^((0.95)){ }^{(0.95)}; b) La2004-La2004 ( 0.90 ) ( 0.90 ) ^((0.90)){ }^{(0.90)}. In c) and d d d\mathbf{d} ), the nominal solution of the obliquity is plotted in solid line, while La2004 ( 0.95 ) c