第29卷第1期天文學(xué)進展VoL 29, No. 12011年1月PROGRESS IN ASTRONOMY文章編號:10008349(2011)01010512主帶小行星的動力學(xué)模擬李嘉中國科學(xué)院上海天文臺,上海200030)摘要:采用接近真實太陽系的動力學(xué)模型,對主帶小行星的動力學(xué)演化進行了數(shù)值模擬。計算的起始時間是儒略日J(rèn)D=24540005×10°,計算的時間長度為100萬年。力學(xué)模型采用n+m體模型,計算程序基于小行星軌道演化的軟件 Orbit9。對演化結(jié)果進行分析可以發(fā)現(xiàn)測試粒子與木星的平運動共振對測試粒子穩(wěn)定性的不同作用,以及在2:3、3:4共振處不同初始值對測試粒子演化結(jié)果的影響關(guān)鍵詞:共振;主帶小行星;軌道演化中圖分類號:P181引言1766年 Johann titius提出 Titius-Bode法則,即以天文單位(AU)計算6顆當(dāng)時已知的大行星到太陽的距離,存在下面的經(jīng)驗公式d≈0.4+0.3×22,i=-0,0,1,2,4,5分析上面的公式,Bode指出在火星和木星中間的i=3處應(yīng)該存在一個尚未被發(fā)現(xiàn)的大行星。1801年意大利人 Piazzi在火星和木星之間發(fā)現(xiàn)第一顆小行星谷神星,在隨后200多年時間里人類在這個區(qū)域陸續(xù)發(fā)現(xiàn)的小行星多達37萬余顆,其中已經(jīng)被編號的有23萬顆左右。這些小行星間軌道半長徑95%集中在217~364AU內(nèi),平均偏心率大約015,平均軌道傾角大約8.6°,被稱為主帶小行星。因為小行星主帶可能保留了太陽系形成初期的一些信息,所以研究它的化學(xué)組成和動力學(xué)演化對于了解太陽系的形成和演化很有意義。而小行星帶的起源至今沒有統(tǒng)一的說法,主要有三種假設(shè):大行星爆炸說,半成品說和彗星形成說。大行星爆炸學(xué)說認(rèn)為小行星是顆爆炸了的行星所留下的殘片;半成品說則認(rèn)為小H中國煤化工星的前身,由CNMHG收稿日期:2010-03-19;修回日期:2010-04-15資助項目:國家自然科學(xué)基金(10633030)天文學(xué)進展卷于在太陽系形成的初期受到木星的引力攝動而停留在星子階段,最后演化成小行星帶;彗星形成說認(rèn)為彗星經(jīng)過長時間的消耗后,其彗發(fā)變短,成為小行星,已有觀測證據(jù)表明某些彗星可能演變成小行星。1867年 Kirkwood對主帶小行星的數(shù)目按軌道半長徑進行統(tǒng)計時發(fā)現(xiàn)了 Kirkwood空隙( Kirkwood gap),在木星軌道周期和小行星軌道周期成4:1、3:1、5:2、7:3、2:1等整數(shù)比的平運動共振處明顯沒有小行星分布;在周期比為1:1、3:2、4:3等平運動共振處則有一定的聚集現(xiàn)象,3:2共振處聚集的小行星被稱作Hlda群,1:1處的稱作 Trojan群。自Kirkwood空隙被發(fā)現(xiàn)以后,解釋小行星帶的這種動力學(xué)結(jié)構(gòu)一直是個難題,隨后發(fā)展出了統(tǒng)計學(xué)假說、碰撞假說、宇宙學(xué)假說和引力假說等不同的假設(shè)。20世紀(jì)80年代以后,隨著計算機技術(shù)的發(fā)展,人們可以進行較長時間的軌道計算,1982年 Widows利用映射方法研究1:3共振處的演化,將演化時間推近到幾百萬年量級后發(fā)現(xiàn)小行星軌道是混沌的,軌道偏心率可以達到03,從而受到火星攝動而離開小行星主帶(最近的研究表明1:3共振處的小行星經(jīng)過一百萬年的演化后,偏心率e可以接近1,從而使小行星軌道穿過地球軌道,由此可認(rèn)為與木星處于平運動共振的小行星帶可能是近地小天體(NEA)的來源之-)。1984年 Hadjidemetriou等人在限制性三體問題的框架下研究一些特殊的共振,發(fā)現(xiàn)1:3和3:5共振處會產(chǎn)生空隙而在1:2、2:3共振處會有小行星存在Dvorak等人分別在199年和1995年,對主帶小行星的動力學(xué)演化進行了模擬。在1993年的文章中,他們采用了太陽和木星的橢圓型限制性三體問題模型;在1995年的文章中,他們引入土星,考慮太陽、木星、土星間的相互作用。通過兩次對不同偏心率和軌道半長徑的小行星的動力學(xué)模擬,他們認(rèn)為主帶的很多動力學(xué)結(jié)構(gòu)(1:3、1:2)共振可以通過限制性三體問題解決,主帶內(nèi)側(cè)的結(jié)構(gòu)由小行星和木星的長期共振決定,從1:2共振位置開始的小行星才會因過于靠近木星而受到木星的直接作用等本文將.用接近真實太陽系的動力學(xué)模型,利用Orb9程序?qū)χ鲙⌒行沁M行動力學(xué)模擬。同時我們重點關(guān)注小行星在初始狀態(tài)時的軌道取向?qū)τ谒旧硌莼挠绊憽?演化模型演化的起始時間為儒略日J(rèn)D=24540005×10°(公元2006年9月22日0點0分0秒),大行星的質(zhì)量和此時的軌道根數(shù)從JPL星表中獲得。太陽系的演化歷史長達幾十億年,但是由于受制于計算機條件,本文中的演化時間長度沒定為10°a太陽系大天體的動力學(xué)模型采用n=7體模型,即將水星、金星質(zhì)量并入太陽質(zhì)量以太陽2項形式出現(xiàn),同時計算新的太聞質(zhì)心;氫略亡中國煤化量、土星、天王星、海王星與太陽間的相互作用。本文將著重關(guān)注初始狀態(tài)下測試粒子的不同軌追對化紹影啊,所以我們將選取若干組不同的近點角距u值;同時測試粒子的軌道半長徑要覆蓋整個主帶的寬度且均1期李嘉:主帶小行星的動力學(xué)模擬107勻分布。我們將測試粒子的初始軌道半長徑選取在15~45AU內(nèi),以6×10-4AU為間隔均勻分布;同時在0°~360°范圍內(nèi)以等差間隔選取5個不同的初始近點角距u值,分別為通過對已知主帶小行星的統(tǒng)計,得到主帶小行星平均偏心率為0.15左右,平均軌道傾角為8.6°左右,所以初始條件選取的偏心率為0.1,軌道傾角為8°。升交點經(jīng)度Ω和平近點角M都設(shè)為0。綜上所述,模擬小行星的初始條件為1)半長徑a均勻分布在主帶15~4.5AU的范圍內(nèi)(因為木星軌道半長徑在52AU左右),間隔為6×104AU,偏心率e=01,軌道傾角讠=8°,?=0°,a=0°、72°、144°、216°、288°,M=0。(2)共取m=25000顆測試粒子,每個值對應(yīng)5000顆。(3)用的模擬程序是Obt93 Orbit9程序介紹Orbit9是一款精密計算太陽系天體軌道演化的數(shù)值軟件,由 Andrea milani, Nobili和Carpino等人共同編寫,可進行n+m體問題的計算。即考慮太陽和n-1個大行星的相互作用,m個質(zhì)量較小的小行星,小行星只考慮其他n體對它們的引力,不考慮它們對n體的引力。對于不同區(qū)域的小行星 Orbit9可以應(yīng)用不同的動力學(xué)模型。例如對于在主帶之外的小行星的計算,可以只利用4個外行星作為攝動源,此時對于內(nèi)行星的作用采用質(zhì)心修正加入到太陽質(zhì)量中,同時這些并入太陽質(zhì)量的內(nèi)行星對小行星的引力作用將以太陽J2項的形式出現(xiàn)(軟件包中帶有執(zhí)行程序,用以計算某一時刻內(nèi)太陽系質(zhì)心修正和考慮內(nèi)行星的J2項此時的太陽和內(nèi)行星的質(zhì)量、根數(shù)等條件都從JPL行星歷表中獲得)。同時,Orbt9還考慮相對論修正。軟件作者的目的是盡可能地使Orb9的動力學(xué)模型接近實際的太陽系情況。Orbit9的積分器以單步法( (implicit Runge-kutt-gaus)作為開始,主要的演化計算由Adams-Cowel多步法完成。多步法有固定的步長,使用者可以采用程序自動計算所得的步長或者自己設(shè)定固定步長。Orb9涉及的坐標(biāo)系有多種形式,使用者可以利用程序所提供的相應(yīng)程序進行轉(zhuǎn)換,般情況下,小行星輸入時的初始條件是太陽系質(zhì)心坐標(biāo)系下的開普勒根數(shù)形式,最后計算結(jié)束的輸出結(jié)果為日心坐標(biāo)系下的第二類無奇點根數(shù)形式。4演化結(jié)果中國煤化工CNMHG經(jīng)過計算,我們得到了測試粒子經(jīng)過10a演化后的軌道根數(shù)。對這些根數(shù)進行統(tǒng)計,我108天文學(xué)進展29卷們得到了不同初始半長徑對應(yīng)的粒子演化后的半長徑,以及偏心率e、軌道傾角i在初始半長徑上的分布,并和目前觀測得到實際情況進行比較。41統(tǒng)計演化后測試粒子的數(shù)量分布從圖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演化后的測試粒子數(shù)目在半長徑上的分布z軸是軌道半長徑,y軸是該初始半長徑處對應(yīng)的粒子數(shù)目.第一個圖是所有粒子放在一起的統(tǒng)計結(jié)果,后面5個分別為初始=0°,72°,14°,218°,288°的統(tǒng)計結(jié)果(1)在1:3、1:4、2:5、1:2等平運動共振處,中國煤化工(2)對于2:3共振,初始值山=288時測試粒CNMHG有較多的測試粒子存在。1期李嘉:主帶小行星的動力學(xué)模擬109(3)對于3:4共振,初始值=144°時還有較多測試粒子存在,其他w值則很少有測試粒子剩下4)在靠近木星的4~45AU附近,已沒有測試粒子存在,主帶內(nèi)側(cè)(22~3AU附近)的測試粒子的數(shù)目則保持穩(wěn)定。與圖2實際情況下的小行星數(shù)目在軌道半長徑上的分布相比,模擬的結(jié)果中各個不同初始山值的1:4、1:3、2:5等共振處粒子數(shù)量劇烈減少,4~4.5AU附近靠近粒子全部逃逸ω=1449時3:4共振處除外,這些結(jié)果與實際情況比較接近,重現(xiàn)了平運動共振和靠近木星的長期攝動對小行星帶的影響。140012001:4M1:38006005a /AU圖2實際情況下小行星數(shù)目在軌道半長徑上的分布x軸是軌道半長徑,y軸是該初始半長徑處對應(yīng)的粒子數(shù)目。4.2統(tǒng)計演化后的測試粒子的軌道半長徑通過圖3可以發(fā)現(xiàn):(1)在各個共振處粒子的半長徑大多演化離開了初始值,例如1:4、1:3共振處的大量測試粒子都進入了火星和地球軌道附近;1:2共振處對軌道半長徑的影響在百萬年量級并不大。(2)非共振處,遠離木星的22~35AU左右處的粒子演化后的半長徑仍然穩(wěn)定在初始位置附近。(3)在靠近木星的4~45AU附近,測試粒子由于木星的攝動已經(jīng)全部逃逸;35~4AU附近的部分測試粒子軌道半長徑演化后有較大變化,較多粒子亭開初始位置43統(tǒng)計演化后測試粒子的偏心率e中國煤化工半長徑由圖4和圖5可以發(fā)現(xiàn)CNMHG(1)在大部分平運動共振處,部分測試粒子的偏心率很多有較大的增長,例如1:4、110天文學(xué)進展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軸是該初始半長徑處對應(yīng)的粒子演化后的軌道半長徑。第一個圖為所有測試粒子的結(jié)果,后面5個圖為各個不同初始ω值的統(tǒng)計結(jié)果,對于演化后超過5AU的軌道半長徑一律將其標(biāo)注在圖中5AU處,對于小于0的軌道半長徑(非開普勒軌道)一律標(biāo)注在0AU處中國煤化工CNMHG李嘉:主帶小行星的動力學(xué)模擬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演化后的偏心率在初始軌道半長徑上統(tǒng)計情況x軸是粒子的初始半長徑,y軸是該初始半長徑處對應(yīng)的粒子演化后的偏心率。第1個圖為所有測試粒子的結(jié)果,后面5圖分別為初始w=0°,72°,144°,218°,288°的統(tǒng)計結(jié)果中國煤化工CNMHG112天文學(xué)進展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°的粒子的偏心率都被激發(fā)離開初始值e=01之外,其他山值還有較多的測試粒子的偏心率仍保持在e=0.3以下。(3)在3:4共振處,初始值山=144時仍然有較多粒子的偏心率維持在e=0.3以下,其他值則大部分都被激發(fā)到很高的值。(4)非平運動共振處,2.2~35AU附近的偏心率較穩(wěn)定,而靠近木星的38~45AU附近的偏心率則被激發(fā)到很高的值,4~45AU附近的測試粒子全部逃逸(5)1.5~2AU區(qū)域內(nèi)的測試粒子的偏心率演化后可達0.4,有可能進入地球軌道44統(tǒng)計演化后的測試粒子的軌道傾角由圖6和圖7可以發(fā)現(xiàn):(1)在部分平運動共振處,測試粒子的軌道傾角很多有較大的增長,例如1:4、1:3、1:2共振處。(2)在2:3共振處,初始值ω=28的粒子的軌道傾角被激發(fā)很少維持在初始值8°附近,其他初始凵值還有較多的測試粒子的軌道傾角在8°附近。(3)在3:4共振處,初始值山=144時仍然有較多粒子的軌道傾角維持在8°附近,其他u值則大部分都被激發(fā)得離開這個區(qū)域(4)非平運動共振處,22~35AU處的軌道傾角較穩(wěn)定,而靠近木星3~45AU的軌道傾角則有很大的變化(5)1.5~2AU區(qū)域內(nèi)的測試粒子的軌道傾角出現(xiàn)不穩(wěn)定情況,相比于初始值有較大變化。綜合43與44與實際情況相比較,模擬結(jié)果(圖5()圖7(b)與實際情況(圖5(a)圖7(a))有一定符合。經(jīng)過10°a量級的動力學(xué)演化1中國煤化工AU附近,測試粒子由于木星的攝動幾乎全部逃逸(初始u=144CNMHG35~4AU附近的測試粒子盡管沒有全部逃逸,但是在木星的攝動下已經(jīng)有部分粒子的偏心率和軌道傾角1期李嘉:主帶小行星的動力學(xué)模擬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在初始軌道半長徑上統(tǒng)計情況z軸是粒子的初始半長徑,v軸是該初始半長徑處對應(yīng)的粒子演化后的軌道傾角第一個分圖為所有測試粒子的結(jié)果,5個子圖分別為初始w=0°,72°,144°,218°,288°的統(tǒng)計結(jié)果中國煤化工CNMHG114天文學(xué)進展29卷16留601:41:32:51:22:33:522.533.544.5半長徑a/AU演化后的半長徑a/AU圖7(a)實際觀測到的主帶小行星軌道傾角在軌道半長徑上的分布;(b)模擬之后的所有測試粒子軌道傾角在演化后軌道半長徑上的分布被激發(fā);主帶內(nèi)側(cè)(22~3AU附近,不包括1:3、2:5等共振點)的測試粒子的偏心率和軌道傾角以及軌道半長徑都保持穩(wěn)定;主帶1:3等平運動共振處,測試粒子發(fā)生逃逸。靠近木星的不穩(wěn)定區(qū)域由于平運動共振2:3、3:4卻仍然有測試粒子存在,同時不同的初始u值對于演化結(jié)果有一定影響:對于2:3共振,在=288°時測試粒子全部逃離共振區(qū),其他4個值則有較多的粒子仍然穩(wěn)定在2:3共振區(qū)域;對于3:4共振,只有w=144°時有較多粒子仍然留在共振區(qū)域,其他值則基本全部逃逸。0.250.2雲(yún)0.1593.913923.933943.953.963.973.983.99a /AU圖8實際觀測到的2:3平運動共振處附近的小行V凵中國煤化工分布CNMHG本文僅做了106a的演化模擬,測試粒子采用的初始模型也僅考慮了軌道半長徑的均勻李嘉:主帶小行星的動力學(xué)模擬分布以及特定的5個近點角距等特殊條件,盡管演化的結(jié)果在某些方面與實際統(tǒng)計結(jié)果有大致的符合,但是考慮到主帶小行星的演化時間是幾十億年,而且主帶小行星的形成原因目前尚無定論,所以結(jié)果只能作為一個參考。實際情況下,偏心率和軌道傾角在半長徑上的分布往往呈現(xiàn)某些規(guī)律,這種規(guī)律的分布下所隱藏的動力學(xué)成因值得我們進行深一步的研究以2:3共振處的偏心率為例(如圖8所示),實際觀測到的在2:3平運動共振處附近的小行星的偏心率與其軌道半長徑存在明顯的函數(shù)關(guān)系。5結(jié)論本文回顧了主帶小行星研究的一些歷史,并粗略介紹了 Orbit9軟件。利用 Orbit9軟件,我們模擬了不同u值和軌道半長徑的25000顆粒子在接近實際情況的太陽系模型下的演化模擬時間為10°a。通過對演化結(jié)果的各種統(tǒng)計分析我們重現(xiàn)了平運動共振和木星的長期攝動對小行星帶的影響,例如1:4、1:3共振處在百萬年量級即可激發(fā)測試粒子的偏心率e和軌道傾角i,使測試粒子進入火星軌道以內(nèi),4~45AU范圍內(nèi)的小行星在木星的攝動下全部逃逸;同時發(fā)現(xiàn)不同的初始u值對于2:3、3:4演化結(jié)果有一定影響,希望為解釋目前這兩個共振處仍然存在小行星提供一定的參考。致謝感謝中科院上海天文臺研究員廖新浩老師的細(xì)致指導(dǎo),也非常感謝兩位審稿老師提出的寶貴意見。參考文獻[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]戴文賽,胡中為,閻林山等太陽系演化學(xué)(下冊)上海:上?？茖W(xué)技術(shù)出版社,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天文學(xué)進展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