刘文芳，胡诗杨,2，刘福窑

（1.上海工程技术大学 数理与统计学院，上海 201620；2.广西大学 物理科学与工程技术学院，南宁 530004）

近几年，激光干涉引力波天文台（LIGO）陆续发现数十个引力波事件[1]，证实了一个世纪前爱因斯坦的预言.目前，捕获的引力波均由两个大质量天体并和所辐射，即波源为双星系统，这类系统广泛存在于宇宙之中.事实上，宇宙中还可能存在由3 个天体组成的孤立系统，这类系统中的天体在相互绕转的过程中同样会辐射引力波，在引力波探测技术不断发展的今天，利用引力波探测三体系统成为一种可能.

要利用引力波这一媒介来研究三体问题，熟知三体模型及其轨道特征是先决条件.三体系统[2−3]的运动状态极不稳定，系统内天体容易发生相互碰撞或逃逸现象.目前只能通过数学手段推断三体系统存在“8 字”型、“蝴蝶”型轨道[4−6]，但在宇宙中并没有发现相应的光学、电磁对应体.当3 个天体中有1 个天体的质量远小于另外两个天体的质量，即两大质量天体不受第三天体的引力影响，这类三体系统称为限制性三体问题[7−8].限制性三体问题在宇宙中有迹可循，如地球、月球、人造卫星组成的三体系统中，人造卫星对地球和月球的引力可以忽略不计，此系统可看作限制性三体问题[9]；太阳、木星及木星的任何一颗卫星组成的三体也可作为限制性三体问题来处理.

限制性三体问题按照两主天体相互绕转的轨迹又可分为圆型限制性三体问题和椭圆型限制性三体问题[10].关于圆型限制性三体问题的研究较为热门.张汉清等[11]借助相空间中的一类流函数，提出一种计算圆型限制性三体问题的周期性轨道的方法；陈云龙等[12]利用力梯度辛算法求解圆型限制性三体问题的哈密顿正则方程，并分析该模型的混沌轨道；石绍伍等[13]研究受扁率摄动和辐射摄动影响的圆型限制性三体问题的轨道稳定性，并指出两种摄动会对混沌轨道造成影响.圆型限制性三体问题存在拉格朗日平动点，如两个三角平动点和3 个共线平动点.当航天器或天体运动至平动点时，达到引力平衡状态，此时平动点上的物体和两主天体保持相对静止.李翔宇等[14]分析圆型限制性三体问题的零速度曲面，研究平动点附近的不同周期性轨道.圆型限制性三体问题在现实生活中的应用极为广泛，航空航天的发展与该模型的轨道特征、拉格朗日平动点密切相关，所以，研究此类问题的轨道动力学对发展近地卫星、绕月卫星、深空探测器有重要意义.

本研究基于圆型限制性三体问题的哈密顿函数，利用数值积分算法模拟不同参数、初始条件的轨道，分析系统参数对零速度曲线的影响.同时，本研究寻找了大量模型的混沌轨道，研究混沌轨道与系统参数、初始条件之间的联系，并阐述混沌轨道产生的原因.

1 圆型限制性三体问题的动力学模型

三体问题由3 个天体构成，各天体的运动规律及整个体系的动力学状态由天体之间的万有引力决定.对于每一个天体，均需要3 维坐标和3 维动量来描述运动状态.所以，三体问题是一个包含18 个变量的动力学系统.在三体问题的基础上，规定其中两天体几乎占据体系所有质量，并绕质心做匀速圆周运动；第三天体在两主天体的轨道平面内运动，其质量及对两主天体的引力作用可以忽略不计.那么，三体问题转变为圆型限制性三体问题.

1.1 旋转坐标系下的三体运动方程

三体系统在旋转坐标系中的运动示意图如图1所示.图中，3 个天体分别标记为Ma、Mb和Mc，对应质量为ma、mb和mc(规定ma>mb≫mc)，那么体系的质量参量μ=mb/(ma+mb).取两主天体Ma和Mb的质心为旋转坐标系原点.两主天体分别位于(−μ,0)和(1 −μ,0)，以角速度ω绕坐标系原点做逆时针匀速圆周运动(蓝色虚线轨迹).同时，整个坐标系以角速度ω做逆时针旋转(灰色虚线轨迹).两主天体在运动过程中始终处于旋转坐标系x轴上，并与初始坐标完全一致.采用旋转坐标系来描述圆型限制性三体问题，可以忽略体系中两主天体的运动，只需关注第三天体Mc，将原本具有18 个变量的三体问题简化为只含有4 个变量的动力学模型，为进一步研究提供便利.

图1 三体系统在旋转坐标系中的运动示意图Fig.1 Sketch map of motion of three-body system in rotating frame

设第三天体Mc的坐标为(x,y)，可知Mc与两主天体Ma、Mb的距离分别为r1、r2，公式为

根据式(1)，写出第三天体在两主天体作用下的引力势能U(x,y)和有效势V(x,y)，公式为

式中：vx和vy分别为第三天体在x和y方向的速度，可得第三天体的拉格朗日量为

用广义动量px、py代替式(4)中的速度，并通过勒让德变换，可得第三天体的无量纲化哈密顿量为

上述哈密顿量的正则方程[15]为

式(6)也为Mc天体的演化方程.给定Mc的初始坐标及动量(x0,y0,px0,py0)，式(6)两边积分可得Mc任意时刻的运动状态.此外，观察式(5)不难发现，哈密顿量中含有ypx和xpy，这是旋转坐标系附加给第三天体的坐标，即动量耦合项(当采用惯性坐标系时，这一耦合部分将消失).由于存在此项，使得哈密顿量式(5)中的坐标、动量不可分离，最终导致式(6)不存在解析解.尽管如此，依然可以利用数值积分算法，如龙格−库塔法、辛算法求解正则方程式(6)，得到第三天体Mc的运动轨迹，进而分析相应的动力学特征.

1.2 雅可比积分和零速度曲线

在旋转坐标系下，第三天体的哈密顿量不显含时间t，所以此哈密顿量为守恒量，并与第三天体的能量等价.根据哈密顿量的全微分，可知哈密顿量与雅可比积分存在代数关系为

式中：Cj为雅可比常数，在一定程度上表征第三天体的能量大小.将雅可比积分中(px+y)和(py−x)固定为0(即速度为0)，可得到第三天体的零速度曲线方程为

由式(8)可得，当给定雅可比常数Cj时，无论第三天体如何演化，其坐标(x,y)必须满足零速度曲线方程，即雅可比常数Cj的取值会影响第三天体的运动范围.取雅可比常数Cj=3.200，质量参量μ=0.012 30，根据式(8)可得到第三天体的零速度曲线如图2(a)所示.

图2 不同参数下第三天体的运动区域和禁行域示意图Fig.2 Sketch maps of motion region and forbidden zone of the third body for different parameters

图中，黑色圆点为两主天体；红色区域的边界线为零速度曲线.当第三天体朝向红色区域运动并到达零速度曲线时，速度变为0.此时，第三天体将在两主天体的引力作用下远离红色区域.显然，第三天体不可能越过零速度曲线到达红色区域内部，所以此红色封闭区域也称为第三天体的禁行域.图2(a)中，禁行域将整个轨道平面大致分为3 个区域.其中，A 区域和B 区域分别为两主天体周围的轨道平面；C 区域为除A、B 两区域及禁行域以外的所有轨道平面.不难发现，A 区域和B 区域由(0.6,0)坐标点附近的“通道”相连，但与C 区域完全隔绝.这说明第三天体的初始位置决定后续的运动范围：当初始位置位于A、B 两区域内，第三天体的运动范围被限制在这两区域之中，无法产生逃逸行为；当初始位置位于C 区域时，第三天体的后续运动只能在C 区域内，此情况下可能发生逃逸现象.

不同参数对应不同的零速度曲线，相应的禁行域也会发生变化.在图2(a)的基础上，略微减小雅可比常数Cj可得此情况下的第三天体禁行域，如图2(b)所示.从图中可以发现，A、B、C 区域均被连通.这意味着即使第三天体的初始位置在两主天体附近，其在后续运动过程中也可以到达C 区域，甚至产生逃逸行为.同理，在图2(a)的基础上，略微减小质量参量μ，可得此参数下的禁行域，如图2(c)所示.这时禁行域将整个轨道平面分割为互不相通的3 个区域，说明第三天体的运动受到更多约束.综上所述，不同系统参数对应不同的零速度曲线和禁行域；较小的质量参量或较大的雅可比常数将缩小第三天体的运动区域；改变系统参数将影响第三天体的轨道特征.

2 圆型限制性三体问题的轨道动力学

雅可比常数Cj和质量参量μ只是在一定程度上限制了第三体的运动范围.对第三体的实际运动轨迹起决定作用的是该天体的初始条件.本节选取数条不同初始条件、不同参数的轨道，通过数值积分算法，模拟各轨道的演化，并分析其中部分轨道的动力学特征.

2.1 轨道构型

在模型(5)中，取12 条不同轨道，各轨道的雅可比常数Cj和初始横坐标x0见表1.每条轨道的质量参量μ固定为0.001 00，初始纵坐标y和动量px均为0，根据式(5)求出初始py值，采用八九阶变步长龙格−库塔法[RKF8(9)]①八九阶变步长龙格−库塔法[RKF8(9)]是一种高精度、可自适应步长的数值积分算法.此算法可以在不同的积分步数中，自动选择最优步长，尽可能减少舍入误差，并在一定积分时间内将系统守恒量（如本文的雅可比常数Cj）误差维持在机器精度.积分式(6)，可以得上述12 条轨道的演化轨迹，如图3 至图5所示.采用此算法可以得到高精度的数值解，避免在后续分析中对系统动力学特征产生误判.

图3 雅可比常数Cj=3.003 时不同轨道示意图Fig.3 Sketch maps of different orbits when Jacobi constant Cj=3.003

图4 雅可比常数Cj=3.010 时不同轨道示意图Fig.4 Sketch maps of different orbits when Jacobi constant Cj=3.010

图5 雅可比常数Cj=3.050 时不同轨道示意图Fig.5 Sketch maps of different orbits when the Jacobi constant Cj=3.050

表1 不同测试轨道的雅克比常数及初始横坐标参数Table 1 Parameters of Jacobian constant and initial x-coordinate values for different test orbits

图3 中，黑色圆点表示两主天体.其中，较大质量主天体MA位于(−0.001,0)；另一主天体MB(质量占整体的0.001)位于(0.999,0)；红色封闭区域和蓝色实线分别为第三天体的禁行域和运动轨迹.当第三天体的起始点靠近较大质量主天体MA时，第三天体会在该主天体附近运动，且能运动至距离该主天体较近的位置.由于远处的主天体MB质量较小，对第三天体的引力也较小，无法将第三天体吸引至自身周围运动.所以，尽管此时禁行域没有隔绝两主天体周围的轨道平面，第三天体也无法在x0较小的情况下在两主天体附近做往返运动，如图3(a)所示.将x0增大至0.35，此时第三天体不仅不能绕较小质量的主天体MB运动，也无法在MA的较低轨道上运行，如图3(b)所示.当第三天体的起始点位于(−1.1,0)时，其后续会沿着整个禁行域的边界运动.尽管禁行域在x轴负半平面留有诸多可供第三天体通过的“缺口”，第三天体仍然无法任意通过这些“缺口”，如图3(c)所示.当第三天体的起始点在MB附近时，第三天体的运动不同于前3 种情况，不再局限于某一区域内，而是存在于更广阔的轨道平面；且此时第三天体不会运动至MA附近，如图3(d)所示.综观图3 不难发现，此参数下第三天体的轨道构型大致分为3 类：环绕MA的拟周期轨道(包含高轨道和低轨道)，环绕禁行域的拟周期轨道，在禁行域外部(远离两主天体一侧)的有界轨道.

由图4 可见，当雅可比常数Cj由3.003 增大至3.010 时，第三天体的轨道构型的类别并没有明显改变.由图4(a)和(b)可见，第三天体在大质量主天体MA附近释放后，将被限制在该主天体周围，围绕该主天体做拟周期运动.当第三天体在靠近较小质量主天体MB附近释放时，其后续将沿着禁行域的边界运动，规律与图3(c)类似.当x0增大至1.05，即第三天体在MB右边释放.此时第三天体受到来自两主天体的引力均较小，可在自身具有一定初速度的情况下运动至较远区域，如图4(d)所示.继续增大雅可比常数Cj至3.050，此时第三天体的轨道构型与之前相比变化较大.由图5 可见，禁行域将整个轨道平面分割为3 个互相独立的区域，第三天体根据初始位置被限制在相应的区域内运动.当x0位于禁行域内侧(靠近MA一侧)时，第三天体的轨迹也只存在于此区域，如图5(a)和图5(b)所示.此参数下，第三天体存在只围绕较小质量主天体运动的轨迹，如图5(c)所示.

综上，圆型限制性三体问题中，第三天体的轨道构型与雅可比常数及初始坐标均有关联，分析这种关系可以提高航天器的入轨成功率及效率，在航空航天中具有重大意义[16].地球、月球、航天器组成的三体系统可近似作为圆型限制性三体问题来处理.地球、月球分别对应两主天体，航天器可看作第三天体，其对两主天体的引力作用可以忽略.雅可比常数和第三天体初始坐标分别对应航天器的能量和变轨位置，可以根据需要的轨道类型(绕月轨道、地球近地轨道等)，通过数值模拟来精确计算出航天器的关机点和变轨点，以及成功入轨所需要的能量.

2.2 混沌轨道

拉普拉斯决定论认为，一个系统在某个时刻的全部状态参数可以决定系统的演化过程.在哈密顿动力学系统中，如果轨道的演化规律由轨道的初始状态决定，那么就说此条轨道符合拉普拉斯决定论.这类轨道具有严格的周期性或拟周期性，称为有序轨道.还有一种轨道，在给定初始状态后，其后续的演化表现出随机性和不可预测性，这类轨道对初始条件极为敏感且不具备周期性，称为混沌轨道.混沌轨道(混沌状态)是目标物体的重要动力学特征之一，寻找这类轨道是分析目标物体运动规律的关键.

可以借助诸多混沌指标来判断轨道状态.常见的混沌指标有最大李雅普诺夫指数、快速李雅普诺夫指标、较小排列指标、余弦指标和功率谱等.其中，最大李雅普诺夫指数判断结果较为准确，被广泛应用于分析物体的非线性现象.常用于计算最大李雅普诺夫指数的方法有两种：一是利用两邻近轨道的相空间距离之差来计算的两粒子法[17−18]；二是通过求解运动方程的变分方程，利用变分来得到最大李雅普诺夫指数的变分法.本研究采用变分法计算最大李雅普诺夫指数λ，相应的表达式[19]为

式中：δ(t)和δ(0)分别为轨道在t时刻和0 时刻的切向量，且δ(0)=1.式(9)说明，最大李雅普诺夫指数λ是一种衡量轨道偏移程度的指标.若轨道处于混沌状态，那么轨道的切向量变化剧烈，这种剧烈程度会体现在变分δ(t)中，并导致λ趋于一个正数；有序轨道的δ(t)和δ(0)时刻保持近似相等，λ趋于0.因此，可以通过λ值判断目标轨道的状态.

根据刘维尔可积定理，哈密顿量式(5)是一个含有两个自由度的保守系统，但只存在一个首次积分式(7)，说明圆型限制性三体问题是一个不可积系统，可能存在混沌轨道①系统不可积是存在混沌轨道的必要非充分条件..固定雅可比常数Cj=3.050，质量参量μ=0.001 00，分别取4 条初始值不同的轨道，计算每条轨道的最大李雅普诺夫指数，如图6 所示.当log10λ随积分时间t的增加而不断减小时，λ趋于0，对应有序轨道，图中黑色轨道和红色轨道；当log10λ趋于稳定时，λ为正数，对应两混沌轨道，图中绿色轨道和蓝色轨道.图6 说明，在圆型限制性三体问题中，第三天体存在混沌轨道，且混沌轨道的出现与第三天体的初始条件有关.

图6 不同轨道的最大李雅普诺夫指数Fig.6 The maximum Lyapunov exponents for different orbits

由图6 可见，当t达到105时，log10λ=−3.5 可作为判断混沌轨道的临界值(图中粉色三角记号).在接下来的分析中，若log10λ≥ −3.5，可以判断轨道处于混沌状态；反之，轨道则为有序.固定雅可比常数为3.050，分别取质量参量μ=0.001 00、0.000 50、0.000 10 和0.000 05，并在区间[0.05,0.80]内扫描x0，可以得到这些轨道在演化时间达到105时的最大李雅普诺夫指数随x0的分布，如图7 所示.图中，红色虚线表示混沌轨道和有序轨道的分界线，位于虚线上方的轨道是混沌轨道，虚线以下的轨道处于有序状态.当μ=0.001 00 时，存在大量混沌轨道，混沌轨道约占轨道总数38%；混沌轨道的分布并不集中，而是较为均匀地散布在整个x0的取值区间，如图7(a)所示.将μ减小至0.000 50，混沌轨道的数目也随之下降，约占轨道总数的25%，如图7(b)所示.此情况下混沌轨道主要分布在区间[0.2,0.6]内.继续减小μ值，此时混沌轨道几乎均匀分布在两主天体连线的中点处，数目只占全体轨道数的6%，如图7(c)所示.图7(d)说明，存在一类较小μ值，只要x0在区间[0.05,0.8]内取值，第三天体均做有序运动.综观图7 不难发现，混沌轨道出现的概率随着μ的减小而减小，这是因为质量参量μ反映了两主天体的质量比：当μ较大时，两主天体质量差距不悬殊，可以在特定位置给第三天体提供方向相反、大小近似的万有引力；两力抗衡的情况下，第三天体容易产生混沌现象；当μ足够小时，两主天体的质量差距过大，较小质量主天体对第三天体的引力作用可以忽略，导致第三天体只在较大质量主天体引力的约束下做严格的有序运动.

图7 不同质量参量下最大李雅普诺夫指数λ 与初始坐标x0 的关系Fig.7 Relations between the maximum Lyapunov exponents λ and initial coordinates x0 for different mass parameters

综上所述，圆型限制性三体问题中，第三天体既存在有序轨道，也存在混沌轨道.混沌轨道产生的原因与两主天体对第三天体的引力有关.混沌轨道的分布范围与系统参数及轨道的初始条件有关，且存在特定的质量参量μ使得第三天体在一定的区域内不存在混沌运动.

3 结语

本研究利用数值积分算法求解不同参数、不同初始条件下圆型限制性三体问题的轨道，并讨论相应的禁行域及动力学特征.研究发现，此模型的轨道构型与雅可比常数和第三天体的初始位置密切相关.当固定质量参量时，较小的雅可比常数对应小范围存在的禁行域，第三天体有机会在较大轨道平面内运动；增大雅可比常数将扩大禁行域的影响范围，进而把第三天体限制在更狭小的区域内运动.同时，禁行域的形状只是限制第三天体运动范围的一个指标，第三天体的实际运动轨迹还与初始条件有关.当第三天体的初始位置位于较大质量的主天体附近时，其后续会环绕该主天体做有界运动.同时，存在一些特定的参数和初始条件，使得第三天体只环绕较小质量主天体运动.

混沌轨道是圆型限制性三体问题中的一类特殊轨道，其产生的原因及分布与系统参数和初始条件有关.在质量参量较大的情况下，两主天体都能给第三天体带来量级相当的引力，这时第三天体容易产生混沌运动.当质量参量足够小时，可以忽略小质量主天体对第三天体的引力作用，相当于第三天体只处于大质量主天体的强引力场中，这种情况一般不容易产生混沌现象.本研究发现的混沌随第三天体初始位置的分布关系为后续在其他参数下寻找混沌轨道的研究提供了指导方向.同时，混沌是宇宙中大部分多天体系统在形成及运行过程中需要经历的一个阶段，研究天体系统的混沌现象可以更深刻地理解天体系统的演化及引力定律.

研究中采用的圆型限制性三体问题只含有牛顿项，事实上，该类模型还含有后牛顿项.后牛顿限制性三体问题是广义相对论下对限制性三体系统更加精确的描述.部分学者研究了轨道半长径和轨道偏心率对后牛顿圆型限制性三体问题的动力学影响[20−21]，但质量参量和雅可比常数对系统动力学特征的影响尚不明确.所以，在后续的研究中，将采用此类模型，分析雅可比常数和质量参量对第三天体轨道的影响.