主带小行星的动力学模拟

第29卷第1期天文学进展VoL 29, No. 12011年1月PROGRESS IN ASTRONOMY文章编号:10008349(2011)01010512主带小行星的动力学模拟李嘉中国科学院上海天文台,上海200030)摘要:采用接近真实太阳系的动力学模型,对主带小行星的动力学演化进行了数值模拟。计算的起始时间是儒略日JD=24540005×10°,计算的时间长度为100万年。力学模型采用n+m体模型,计算程序基于小行星轨道演化的软件 Orbit9。对演化结果进行分析可以发现测试粒子与木星的平运动共振对测试粒子稳定性的不同作用,以及在2:3、3:4共振处不同初始值对测试粒子演化结果的影响关键词:共振;主带小行星;轨道演化中图分类号:P181引言1766年 Johann titius提出 Titius-Bode法则,即以天文单位(AU)计算6颗当时已知的大行星到太阳的距离,存在下面的经验公式d≈0.4+0.3×22,i=-0,0,1,2,4,5分析上面的公式,Bode指出在火星和木星中间的i=3处应该存在一个尚未被发现的大行星。1801年意大利人 Piazzi在火星和木星之间发现第一颗小行星谷神星,在随后200多年时间里人类在这个区域陆续发现的小行星多达37万余颗,其中已经被编号的有23万颗左右。这些小行星间轨道半长径95%集中在217~364AU内,平均偏心率大约015,平均轨道倾角大约8.6°,被称为主带小行星。因为小行星主带可能保留了太阳系形成初期的一些信息,所以研究它的化学组成和动力学演化对于了解太阳系的形成和演化很有意义。而小行星带的起源至今没有统一的说法,主要有三种假设:大行星爆炸说,半成品说和彗星形成说。大行星爆炸学说认为小行星是颗爆炸了的行星所留下的残片;半成品说则认为小H中国煤化工星的前身,由CNMHG收稿日期:2010-03-19;修回日期:2010-04-15资助项目:国家自然科学基金(10633030)天文学进展卷于在太阳系形成的初期受到木星的引力摄动而停留在星子阶段,最后演化成小行星带;彗星形成说认为彗星经过长时间的消耗后,其彗发变短,成为小行星,已有观测证据表明某些彗星可能演变成小行星。1867年 Kirkwood对主带小行星的数目按轨道半长径进行统计时发现了 Kirkwood空隙( Kirkwood gap),在木星轨道周期和小行星轨道周期成4:1、3:1、5:2、7:3、2:1等整数比的平运动共振处明显没有小行星分布;在周期比为1:1、3:2、4:3等平运动共振处则有一定的聚集现象,3:2共振处聚集的小行星被称作Hlda群,1:1处的称作 Trojan群。自Kirkwood空隙被发现以后,解释小行星带的这种动力学结构一直是个难题,随后发展出了统计学假说、碰撞假说、宇宙学假说和引力假说等不同的假设。20世纪80年代以后,随着计算机技术的发展,人们可以进行较长时间的轨道计算,1982年 Widows利用映射方法研究1:3共振处的演化,将演化时间推近到几百万年量级后发现小行星轨道是混沌的,轨道偏心率可以达到03,从而受到火星摄动而离开小行星主带(最近的研究表明1:3共振处的小行星经过一百万年的演化后,偏心率e可以接近1,从而使小行星轨道穿过地球轨道,由此可认为与木星处于平运动共振的小行星带可能是近地小天体(NEA)的来源之-)。1984年 Hadjidemetriou等人在限制性三体问题的框架下研究一些特殊的共振,发现1:3和3:5共振处会产生空隙而在1:2、2:3共振处会有小行星存在Dvorak等人分别在199年和1995年,对主带小行星的动力学演化进行了模拟。在1993年的文章中,他们采用了太阳和木星的椭圆型限制性三体问题模型;在1995年的文章中,他们引入土星,考虑太阳、木星、土星间的相互作用。通过两次对不同偏心率和轨道半长径的小行星的动力学模拟,他们认为主带的很多动力学结构(1:3、1:2)共振可以通过限制性三体问题解决,主带内侧的结构由小行星和木星的长期共振决定,从1:2共振位置开始的小行星才会因过于靠近木星而受到木星的直接作用等本文将.用接近真实太阳系的动力学模型,利用Orb9程序对主带小行星进行动力学模拟。同时我们重点关注小行星在初始状态时的轨道取向对于它本身演化的影响。2演化模型演化的起始时间为儒略日JD=24540005×10°(公元2006年9月22日0点0分0秒),大行星的质量和此时的轨道根数从JPL星表中获得。太阳系的演化历史长达几十亿年,但是由于受制于计算机条件,本文中的演化时间长度没定为10°a太阳系大天体的动力学模型采用n=7体模型,即将水星、金星质量并入太阳质量以太阳2项形式出现,同时计算新的太闻质心;氢略亡中国煤化量、土星、天王星、海王星与太阳间的相互作用。本文将着重关注初始状态下测试粒子的不同轨追对化绍影啊,所以我们将选取若干组不同的近点角距u值;同时测试粒子的轨道半长径要覆盖整个主带的宽度且均1期李嘉:主带小行星的动力学模拟107匀分布。我们将测试粒子的初始轨道半长径选取在15~45AU内,以6×10-4AU为间隔均匀分布;同时在0°~360°范围内以等差间隔选取5个不同的初始近点角距u值,分别为通过对已知主带小行星的统计,得到主带小行星平均偏心率为0.15左右,平均轨道倾角为8.6°左右,所以初始条件选取的偏心率为0.1,轨道倾角为8°。升交点经度Ω和平近点角M都设为0。综上所述,模拟小行星的初始条件为1)半长径a均匀分布在主带15~4.5AU的范围内(因为木星轨道半长径在52AU左右),间隔为6×104AU,偏心率e=01,轨道倾角讠=8°,?=0°,a=0°、72°、144°、216°、288°,M=0。(2)共取m=25000颗测试粒子,每个值对应5000颗。(3)用的模拟程序是Obt93 Orbit9程序介绍Orbit9是一款精密计算太阳系天体轨道演化的数值软件,由 Andrea milani, Nobili和Carpino等人共同编写,可进行n+m体问题的计算。即考虑太阳和n-1个大行星的相互作用,m个质量较小的小行星,小行星只考虑其他n体对它们的引力,不考虑它们对n体的引力。对于不同区域的小行星 Orbit9可以应用不同的动力学模型。例如对于在主带之外的小行星的计算,可以只利用4个外行星作为摄动源,此时对于内行星的作用采用质心修正加入到太阳质量中,同时这些并入太阳质量的内行星对小行星的引力作用将以太阳J2项的形式出现(软件包中带有执行程序,用以计算某一时刻内太阳系质心修正和考虑内行星的J2项此时的太阳和内行星的质量、根数等条件都从JPL行星历表中获得)。同时,Orbt9还考虑相对论修正。软件作者的目的是尽可能地使Orb9的动力学模型接近实际的太阳系情况。Orbit9的积分器以单步法( (implicit Runge-kutt-gaus)作为开始,主要的演化计算由Adams-Cowel多步法完成。多步法有固定的步长,使用者可以采用程序自动计算所得的步长或者自己设定固定步长。Orb9涉及的坐标系有多种形式,使用者可以利用程序所提供的相应程序进行转换,般情况下,小行星输入时的初始条件是太阳系质心坐标系下的开普勒根数形式,最后计算结束的输出结果为日心坐标系下的第二类无奇点根数形式。4演化结果中国煤化工CNMHG经过计算,我们得到了测试粒子经过10a演化后的轨道根数。对这些根数进行统计,我108天文学进展29卷们得到了不同初始半长径对应的粒子演化后的半长径,以及偏心率e、轨道倾角i在初始半长径上的分布,并和目前观测得到实际情况进行比较。41统计演化后测试粒子的数量分布从图1可以看到:1401:41:32:5;1:22:33:4050藏351:4132:512P334522.533.544.51.5225335445a/AUJAU501:4132:51:22:33:41;41:32:51:22:33:425252005051.522.533.5L522.533.544.55a/AUa/AU351:41:32:51:241:32:5122:33:4嶽2520050505051051522.533.544.551.522.533.544.5JAUa /AU图1演化后的测试粒子数目在半长径上的分布z轴是轨道半长径,y轴是该初始半长径处对应的粒子数目.第一个图是所有粒子放在一起的统计结果,后面5个分别为初始=0°,72°,14°,218°,288°的统计结果(1)在1:3、1:4、2:5、1:2等平运动共振处,中国煤化工(2)对于2:3共振,初始值山=288时测试粒CNMHG有较多的测试粒子存在。1期李嘉:主带小行星的动力学模拟109(3)对于3:4共振,初始值=144°时还有较多测试粒子存在,其他w值则很少有测试粒子剩下4)在靠近木星的4~45AU附近,已没有测试粒子存在,主带内侧(22~3AU附近)的测试粒子的数目则保持稳定。与图2实际情况下的小行星数目在轨道半长径上的分布相比,模拟的结果中各个不同初始山值的1:4、1:3、2:5等共振处粒子数量剧烈减少,4~4.5AU附近靠近粒子全部逃逸ω=1449时3:4共振处除外,这些结果与实际情况比较接近,重现了平运动共振和靠近木星的长期摄动对小行星带的影响。140012001:4M1:38006005a /AU图2实际情况下小行星数目在轨道半长径上的分布x轴是轨道半长径,y轴是该初始半长径处对应的粒子数目。4.2统计演化后的测试粒子的轨道半长径通过图3可以发现:(1)在各个共振处粒子的半长径大多演化离开了初始值,例如1:4、1:3共振处的大量测试粒子都进入了火星和地球轨道附近;1:2共振处对轨道半长径的影响在百万年量级并不大。(2)非共振处,远离木星的22~35AU左右处的粒子演化后的半长径仍然稳定在初始位置附近。(3)在靠近木星的4~45AU附近,测试粒子由于木星的摄动已经全部逃逸;35~4AU附近的部分测试粒子轨道半长径演化后有较大变化,较多粒子亭开初始位置43统计演化后测试粒子的偏心率e中国煤化工半长径由图4和图5可以发现CNMHG(1)在大部分平运动共振处,部分测试粒子的偏心率很多有较大的增长,例如1:4、110天文学进展D半升554532.52.5!只游1.50.51:41:32:51:22:33:4051:41:32:51:22:33:4初始半长径a/AU初始半长径a/AU4.54.5:;升名吗3525151:41:32:51:22:33:41:41:32:51:22:33:4522.53.544.552.54.55初始半长径a/AU初始半长径a/AU553525255453521.50.51:41:32:51:22:33:40.51:41:32:5122:33:422.533.544.551.522.533.544.55初始半长径a/AU初始半长径a/AU图3演化后的粒子轨道半长径情况x轴是粒子的初始半长径,y轴是该初始半长径处对应的粒子演化后的轨道半长径。第一个图为所有测试粒子的结果,后面5个图为各个不同初始ω值的统计结果,对于演化后超过5AU的轨道半长径一律将其标注在图中5AU处,对于小于0的轨道半长径(非开普勒轨道)一律标注在0AU处中国煤化工CNMHG李嘉:主带小行星的动力学模拟11100.98765432080.70000002102:33:4:41:32:51:22:33;4522.533.544.55522.533.544.55初始半长径a/AU初始半长径a/AU0.840.7·0.7到:0.6060.50.50003211:41:32:51:22:33:4210111:41:32:51:22:33:41.522.53.544.5初始半长径a/AU初始半长径a/AU000987600000.8760.50.5b0求030.11:41:32:51:22:33:441:32:51:22:33:41.522.533.544.5522.533.544.55初始半长径a/AU初始半长径a/AU图4演化后的偏心率在初始轨道半长径上统计情况x轴是粒子的初始半长径,y轴是该初始半长径处对应的粒子演化后的偏心率。第1个图为所有测试粒子的结果,后面5图分别为初始w=0°,72°,144°,218°,288°的统计结果中国煤化工CNMHG112天文学进展29卷06(a)09}(b)0.80.70500.21.522.533.5445521.522.533.544.55a/AU演化后的半长径a/AU图5(a)实际观测到的主带小行星偏心率在轨道半长径上的分布;(b)模拟之后的所有测试粒子偏心率在演化后轨道半长径上的分布1:3、1:2共振处。(2)在2:3共振处,除了初始值山=288°的粒子的偏心率都被激发离开初始值e=01之外,其他山值还有较多的测试粒子的偏心率仍保持在e=0.3以下。(3)在3:4共振处,初始值山=144时仍然有较多粒子的偏心率维持在e=0.3以下,其他值则大部分都被激发到很高的值。(4)非平运动共振处,2.2~35AU附近的偏心率较稳定,而靠近木星的38~45AU附近的偏心率则被激发到很高的值,4~45AU附近的测试粒子全部逃逸(5)1.5~2AU区域内的测试粒子的偏心率演化后可达0.4,有可能进入地球轨道44统计演化后的测试粒子的轨道倾角由图6和图7可以发现:(1)在部分平运动共振处,测试粒子的轨道倾角很多有较大的增长,例如1:4、1:3、1:2共振处。(2)在2:3共振处,初始值ω=28的粒子的轨道倾角被激发很少维持在初始值8°附近,其他初始凵值还有较多的测试粒子的轨道倾角在8°附近。(3)在3:4共振处,初始值山=144时仍然有较多粒子的轨道倾角维持在8°附近,其他u值则大部分都被激发得离开这个区域(4)非平运动共振处,22~35AU处的轨道倾角较稳定,而靠近木星3~45AU的轨道倾角则有很大的变化(5)1.5~2AU区域内的测试粒子的轨道倾角出现不稳定情况,相比于初始值有较大变化。综合43与44与实际情况相比较,模拟结果(图5()图7(b)与实际情况(图5(a)图7(a))有一定符合。经过10°a量级的动力学演化1中国煤化工AU附近,测试粒子由于木星的摄动几乎全部逃逸(初始u=144CNMHG35~4AU附近的测试粒子尽管没有全部逃逸,但是在木星的摄动下已经有部分粒子的偏心率和轨道倾角1期李嘉:主带小行星的动力学模拟113188641001:41:32:51:22:31:41:32:51:22:33:4.522.533.544.522.533.544.55初始半长径a/AU初始半长径a/AU116:a14121:41:32:51:22:33:41:41:32:51:22:33:41.522.533.54.55初始半长径a/AU初始半长径a/AU18归1612堡本4208OF21:41:32:51:22:33:41:41:32:51:22:33:4-21.522.533.544.551.522.533.54初始半长径a/AU初始半长径a/AU图6演化后的轨道倾角i在初始轨道半长径上统计情况z轴是粒子的初始半长径,v轴是该初始半长径处对应的粒子演化后的轨道倾角第一个分图为所有测试粒子的结果,5个子图分别为初始w=0°,72°,144°,218°,288°的统计结果中国煤化工CNMHG114天文学进展29卷16留601:41:32:51:22:33:522.533.544.5半长径a/AU演化后的半长径a/AU图7(a)实际观测到的主带小行星轨道倾角在轨道半长径上的分布;(b)模拟之后的所有测试粒子轨道倾角在演化后轨道半长径上的分布被激发;主带内侧(22~3AU附近,不包括1:3、2:5等共振点)的测试粒子的偏心率和轨道倾角以及轨道半长径都保持稳定;主带1:3等平运动共振处,测试粒子发生逃逸。靠近木星的不稳定区域由于平运动共振2:3、3:4却仍然有测试粒子存在,同时不同的初始u值对于演化结果有一定影响:对于2:3共振,在=288°时测试粒子全部逃离共振区,其他4个值则有较多的粒子仍然稳定在2:3共振区域;对于3:4共振,只有w=144°时有较多粒子仍然留在共振区域,其他值则基本全部逃逸。0.250.2雲0.1593.913923.933943.953.963.973.983.99a /AU图8实际观测到的2:3平运动共振处附近的小行V凵中国煤化工分布CNMHG本文仅做了106a的演化模拟,测试粒子采用的初始模型也仅考虑了轨道半长径的均匀李嘉:主带小行星的动力学模拟分布以及特定的5个近点角距等特殊条件,尽管演化的结果在某些方面与实际统计结果有大致的符合,但是考虑到主带小行星的演化时间是几十亿年,而且主带小行星的形成原因目前尚无定论,所以结果只能作为一个参考。实际情况下,偏心率和轨道倾角在半长径上的分布往往呈现某些规律,这种规律的分布下所隐藏的动力学成因值得我们进行深一步的研究以2:3共振处的偏心率为例(如图8所示),实际观测到的在2:3平运动共振处附近的小行星的偏心率与其轨道半长径存在明显的函数关系。5结论本文回顾了主带小行星研究的一些历史,并粗略介绍了 Orbit9软件。利用 Orbit9软件,我们模拟了不同u值和轨道半长径的25000颗粒子在接近实际情况的太阳系模型下的演化模拟时间为10°a。通过对演化结果的各种统计分析我们重现了平运动共振和木星的长期摄动对小行星带的影响,例如1:4、1:3共振处在百万年量级即可激发测试粒子的偏心率e和轨道倾角i,使测试粒子进入火星轨道以内,4~45AU范围内的小行星在木星的摄动下全部逃逸;同时发现不同的初始u值对于2:3、3:4演化结果有一定影响,希望为解释目前这两个共振处仍然存在小行星提供一定的参考。致谢感谢中科院上海天文台研究员廖新浩老师的细致指导,也非常感谢两位审稿老师提出的宝贵意见。参考文献[2]http://hamilton.dm.unipi.it/astdys/index.php?pc=1.0.0,2010cambridgeUinversityPress,1999:5[1] Murray C D, Dermott S F Solar System Dynamics, Cambridge3]戴文赛,胡中为,阎林山等太阳系演化学(下册)上海:上海科学技术出版社,1986:18044ttp://www.nasa.gov/centers/ames/research/exploringtheuniverse/comets.html20055 SYLVIO FERRAZ-MELLO. IAUS, 1994: 175[6] Greenberg R, Scholl H. Asteroids. Tucson: University of Arizona Press, 1979: 3107] Wisdom j.AJ,1982,87:577[8] Morbidelli A, Moons M. Icarus, 1995, 115:60[9 Hadjidemetriou J D, Ichtiaroglou S. A&A, 1984, 131 20[10 Muller P, Dvorak R, Kallrath J. A& A, 1993, 274: 62711 Muller P, Dvorak R. A&A, 1995, 300: 289[12]http://hamilton.dm.unipi.it/astdys/index.phppc=5,2010[13 Nobili A M, Milani A, Carpino M, A&A, 1989, 210: 316中国煤化工CNMHG116天文学进展29卷Dynamical Simulation of Asteroids in the Main BeltLⅠJiaanghai Astronomical Observatory, Chinese Academy of Science, Shanghai 20030, China)Abstract: By using a dynamical model approach of real solar system, a numerical simulationof the dynamical evolution of the main belt asteroids in solar system is made. The evolutionof the test particles with a time span of 1 million years from the beginning of Julian dayJD=2.454 5 x 10 by using the Orbit9 integrator, which is designed specially to computethe orbit evolution of small bodies, is calculated. The dynamical model used in this paperis a system of n+m bodies, where n equals 7 which means the mass of Mercury and Venusare merged into the mass of Sun and the mass of Pluto is ignored, meanwhile the interactionamong the Sun and the rest 6 planets is considered; the number of test particles, m, equals25 000 in this model. The model for the initial distribution of the test particles is: there are5 groups of test particles with different argument of pericenter ( w) which are 0, 720, 1440216 and 288, the semi major axis(a) of each group distributes between 1.5 4.5 aU evenlywith an interval of 6X10-4 AU, the eccentricity(e) and orbital inclination(i)of all particlesare set to be the same values of e=0. 1 and i=8 as they distribute in the solar system, thelongitude of Ascending Node 2=0, and the mean anomaly M=0. The results show thedifferent effects of mean motion resonance between the Jupiter and the test particles on thestability of the test particles and the infiuence of different initial value of w at 2: 3 and 3: 4resonant zones on the evolution of the test particlesKey words: resonance; main belt asteroids; orbital evolution中国煤化工CNMHG