熱海王星系統(tǒng)HD 106315的近2:1平運(yùn)動(dòng)共振捕獲與軌道演化?

黃秀敏 季江徽 董 瑤

(1 中國科學(xué)院紫金山天文臺(tái)南京210023)(2 中國科學(xué)院行星科學(xué)重點(diǎn)實(shí)驗(yàn)室南京210023)(3 中國科學(xué)技術(shù)大學(xué)天文與空間科學(xué)學(xué)院合肥230026)

1 引言

1995年, Mayor和Queloz[1]利用視向速度測量法在類太陽型恒星51 Pegasi周圍發(fā)現(xiàn)第1顆熱木星, 他們也因此分享了2019年諾貝爾物理學(xué)獎(jiǎng). 隨著測量精度的不斷提高, 截至2019年12月天文學(xué)家通過視向速度法(Radial Velocity)、凌星法(Transit)、微引力透鏡法(Microlensing)、直接成像法(Direct Imaging)、天體測量法(Astrometry)等手段[2],發(fā)現(xiàn)并確認(rèn)4100多顆系外行星(exoplanet.eu). 這些行星大小不一形態(tài)各異, 包括熱木星(Hot Jupiters)、亞海王星、類地行星、超級地球(super-Earth)等類型, 它們怎樣形成以及軌道構(gòu)型如何演化, 無疑是當(dāng)前系外行星研究的前沿科學(xué)問題.

在發(fā)現(xiàn)的系外行星中, 有一類天體的質(zhì)量接近海王星且其離恒星的平均距離約為1 au (天文單位)以內(nèi), 通常我們稱之為熱海王星(Hot Neptunes). 目前觀測統(tǒng)計(jì)表明, 在已發(fā)現(xiàn)的300余顆熱海王星中偏心率較大(>0.2)的約占1/4, 分別分布于27個(gè)不同的行星系統(tǒng)[3].

作為Kepler空間探測任務(wù)的延續(xù), K2在2014年至2018年3月完成了16期觀測工作并發(fā)現(xiàn)數(shù)百顆系外行星[4–5]. 而本文的研究對象熱海王星系統(tǒng)HD 106315正由K2任務(wù)所發(fā)現(xiàn). 恒星HD 106315距離地球(107.3±3.9)pc, 它是一顆質(zhì)量和半徑與太陽相似的F光譜型恒星, 年齡約為40億年, 其質(zhì)量為(1.07±0.03)M⊙, 半徑為(1.18±0.11)R⊙, 有效溫度為(6290±60)K[6]. 該行星系統(tǒng)中包含兩顆大偏心率的熱海王星, 表1中列出了HD 106315系統(tǒng)兩顆行星的軌道觀測數(shù)據(jù), 其中,M為行星質(zhì)量,M⊕代表地球質(zhì)量,Rp為行星半徑,R⊕代表地球半徑,R?為中央恒星半徑,P為行星軌道周期,i和ω分別為軌道傾角和近星點(diǎn)幅角. HD 106315b和HD 106315c在密近軌道上(軌道半長徑a分別為0.09012 au和0.1526 au)以較大的軌道偏心率e(分別為0.3和0.23)圍繞明亮的中央恒星運(yùn)動(dòng). 這兩顆行星軌道周期比接近2, 故它們圍繞中央恒星的軌道呈近2 : 1平運(yùn)動(dòng)共振構(gòu)型. Yimathlmaz等[7]還指出行星HD 106315b的主要組成為50%的水和50%的巖石, 介于氣態(tài)行星和巖石行星之間. 結(jié)合視向速度測量結(jié)果, 這兩顆行星的質(zhì)量、內(nèi)部結(jié)構(gòu)及大氣特征值得關(guān)注. 其他研究表明, 該行星系統(tǒng)在更長周期的軌道上還可能存在第3顆行星[6].

表1 HD 106315系統(tǒng)行星參數(shù)[6]Table 1 Planetary parameters of HD 106315 system[6]

通過對Kepler多年觀測的數(shù)據(jù)分析[9], 科學(xué)家發(fā)現(xiàn)有很多系外行星系統(tǒng)呈現(xiàn)平運(yùn)動(dòng)共振(Mean Motion Resonances, MMRs)軌道構(gòu)型. 法國數(shù)學(xué)家和物理學(xué)家Laplace最早將軌道共振比率定義為兩顆行星圍繞同一中央恒星在相同的時(shí)間間隔內(nèi)完成軌道數(shù)的比例, 即兩星軌道平運(yùn)動(dòng)速率n′和n呈整數(shù)比:n′/n=p/(p+q)(p、q均為正整數(shù)). 此時(shí),行星軌道周期也呈現(xiàn)對應(yīng)的整數(shù)比(p+q)/p. 系外行星中常見的共振軌道構(gòu)型有2 : 1、3 : 2共振[10], 其中2 : 1平運(yùn)動(dòng)共振十分普遍[11–13], 凌星法發(fā)現(xiàn)的行星系統(tǒng)中有16%顯示為近似2 : 1平運(yùn)動(dòng)共振行星系統(tǒng), 2 : 1共振也被認(rèn)為是行星系統(tǒng)演化過程中比較穩(wěn)定的軌道構(gòu)型. 此外還有一類特殊的平運(yùn)動(dòng)共振—3體共振(也稱為Laplace共振)[14–15], 是指系統(tǒng)內(nèi)3顆行星兩兩處于2 : 1平運(yùn)動(dòng)共振, 例如太陽系內(nèi)木星的3顆衛(wèi)星Ganymede、Europa、Io就是4:2:1軌道構(gòu)型.

在系外行星系統(tǒng)中與HD 106315行星系統(tǒng)類似的近2:1平運(yùn)動(dòng)共振有很多, 例如GJ 876、GJ 1061、HD 136352等. GJ 876[16]系統(tǒng)中的GJ 876b和GJ 876c質(zhì)量接近木星質(zhì)量, 而HD 106315b和HD 106315c為海王星質(zhì)量行星. 在考慮行星與原行星盤相互作用時(shí), GJ 876行星系統(tǒng)由于質(zhì)量更大而可能形成原行星盤內(nèi)間隙環(huán)(gap), 從而使行星之間的盤內(nèi)物質(zhì)被快速清除, 抑制了外部軌道行星進(jìn)一步向內(nèi)遷移. 對于HD 106315行星系統(tǒng), 需要判斷行星質(zhì)量是否滿足產(chǎn)生間隙環(huán)的條件, 進(jìn)而影響后續(xù)的軌道遷移過程. 同時(shí), HD 106315行星系統(tǒng)構(gòu)型更緊湊, 恒星引力和行星間的攝動(dòng)作用對行星共振遷移的影響也會(huì)更顯著.

對于軌道半長徑小于0.1 au的熱海王星, 來自中央恒星的潮汐耗散效應(yīng)也會(huì)對軌道半長徑和偏心率產(chǎn)生影響. 與軌道遷移相比, 潮汐耗散作用時(shí)標(biāo)更長, 對于海王星質(zhì)量的行星, 該潮汐耗散時(shí)標(biāo)通常大于106yr. 通過對HD 106315軌道演化機(jī)制的研究, 我們可深入了解大偏心率熱海王星系統(tǒng)演化歷史及當(dāng)前系統(tǒng)穩(wěn)定性.

綜上所述, 本文工作將圍繞熱海王星系統(tǒng)HD 106315近平運(yùn)動(dòng)共振構(gòu)型的形成機(jī)制和軌道遷移模型展開. 第2節(jié)簡要介紹HD 106315系統(tǒng)的軌道演化動(dòng)力學(xué)模型, 包括I類軌道遷移和共振捕獲. 第3節(jié)給出本文使用的參數(shù)空間和數(shù)值模擬方法, 并據(jù)此得到觀測值對應(yīng)的近2:1共振軌道構(gòu)型及初始軌道參數(shù)范圍. 第4節(jié)考慮了潮汐耗散效應(yīng), 結(jié)合半分析方法和數(shù)值模擬, 討論HD 106315系統(tǒng)在不同軌道構(gòu)型下的潮汐演化結(jié)果. 第5節(jié)總結(jié)了HD 106315系統(tǒng)軌道遷移機(jī)制和潮汐效應(yīng)對軌道共振構(gòu)型的影響.

2 軌道遷移模型

2.1 HD 106315系統(tǒng)的I類軌道遷移

一般認(rèn)為, 行星起源于圍繞在年輕恒星周圍的氣體和塵埃盤. 在恒星引力主導(dǎo)的收縮過程中, 氣體和塵埃因具有不同的角動(dòng)量, 被分配到空間不同位置. 角動(dòng)量小的物質(zhì)被吸引到系統(tǒng)中心, 角動(dòng)量大的物質(zhì)形成原行星盤. 根據(jù)Hayashi[17]等提出的原行星盤經(jīng)典物理模型, 原行星盤直徑范圍為10–1000 au, 壽命約106–107yr.

根據(jù)系外行星原位起源(in-situ)模型和Lee等[18]的工作, 在距離中央恒星大約0.1 au處質(zhì)量為10M⊕的原行星可以達(dá)到雪崩吸積(Runaway accretion)的臨界核質(zhì)量, 在百萬年內(nèi)快速吸積氣體生長為氣態(tài)巨行星. Hansen等[19]和Chiang等[20]提出了類海王星行星的原位吸積模型, 行星在當(dāng)前軌道附近直接通過塵埃吸積和行星子(planetesimals)并合形成. 他們認(rèn)為, 目前觀測到的超級地球和類海王星行星形成于氣體耗散殆盡的過渡原行星盤中, 一般不會(huì)經(jīng)歷劇烈的角動(dòng)量交換過程, 從而易保持行星形成時(shí)刻的近圓軌道.

HD 106315系統(tǒng)內(nèi)兩顆行星均處于短周期密近軌道, 但是軌道偏心率較大. 這意味著該行星系統(tǒng)更可能通過軌道遷移機(jī)制形成, 且軌道遷移過程中可能存在不穩(wěn)定的軌道偏心率激發(fā)和衰減. Terquem等[21]通過對質(zhì)量為0.1～1M⊕的原行星軌道演化進(jìn)行數(shù)值模擬, 發(fā)現(xiàn)它們可以并合形成海王星質(zhì)量的行星, 且這些類海王星行星相伴出現(xiàn), 形成多行星系統(tǒng). Cresswell等[22]指出, 對于無法形成原行星盤環(huán)形空隙的類海王星行星, 其軌道偏心率和軌道傾角受盤內(nèi)氣體耗散作用的衰減時(shí)標(biāo)遠(yuǎn)小于軌道遷移時(shí)標(biāo).

原行星盤驅(qū)動(dòng)的軌道遷移模型中, 行星和氣體盤之間產(chǎn)生不可逆的角動(dòng)量交換,軌道遷移的速率是由氣體盤的物理參數(shù)和行星質(zhì)量決定的. 根據(jù)行星的質(zhì)量不同, 一般分為兩種類型的軌道遷移機(jī)制: 第I類遷移(type I migration)和第II類遷移(type II migration). 類地行星和類海王星往往啟動(dòng)速度較快的第I類遷移, 質(zhì)量較大的氣態(tài)巨行星則會(huì)發(fā)生第II類軌道遷移, 且在其軌道附近形成空隙環(huán), 使行星周圍塵埃和氣體密度下降, 共旋力矩(corotation torques)和Lindblad力矩[23]影響減小, 進(jìn)而減慢行星遷移速率.

對于HD 106315系統(tǒng), 考慮I類軌道遷移. I類軌道遷移速度較快, 地球質(zhì)量的原行星可在105yr內(nèi)到達(dá)中央恒星處. I類軌道遷移時(shí)標(biāo)具體還取決于行星盤的結(jié)構(gòu)、行星質(zhì)量比、行星盤表面密度等. 同時(shí), 兩顆行星軌道傾角差值較小, 根據(jù)行星盤內(nèi)的軌道遷移平面效應(yīng)理論, 可以假設(shè)軌道遷移初始時(shí)刻兩顆行星共面.

針對HD 106315行星系統(tǒng)從行星誕生后不久開始的軌道遷移過程, 我們首先考慮原行星盤內(nèi)殘余的氣體和塵埃對行星軌道遷移動(dòng)力學(xué)模型的影響. 原行星盤通過這些作用力對行星施加外部力矩, 為了保證這類力矩的持續(xù)性, 需要行星在原行星盤內(nèi)穩(wěn)定地遷移并且沒有較大質(zhì)量的吸積. 根據(jù)Masset等[24]早期工作中提出的觀點(diǎn), 行星受到的共旋力矩與其軌道遷移速率成比例, 且會(huì)對軌道遷移產(chǎn)生正反饋, 加快遷移速率.

在行星遷移的初始時(shí)刻, 通常假設(shè)兩顆行星位于相距較遠(yuǎn)的兩個(gè)軌道上, 以使行星軌道遷移的時(shí)標(biāo)盡量小于行星盤星云物質(zhì)完全耗散所需的時(shí)間. 一般可以假設(shè)行星b和行星c的初始軌道半長徑a之比為ab:ac≈1:2. 行星遷移的時(shí)標(biāo)是行星盤粘度系數(shù)的函數(shù), 即原行星盤內(nèi)軌道遷移的速率·a是行星盤粘滯系數(shù)v的函數(shù)[25]:

式中原行星盤粘滯系數(shù)[26]v=αH2?,α是用于表征原行星盤粘性程度的無量綱參數(shù),H為原行星盤標(biāo)高, 此處?是在行星軌道半長徑為a時(shí)的Kepler公轉(zhuǎn)角速度. 原行星盤內(nèi)可能的粘性來源主要有: 磁轉(zhuǎn)動(dòng)不穩(wěn)定性(Magnetorotational Instability, MRI)激發(fā)的磁流體動(dòng)力學(xué)擾動(dòng)以及由地球質(zhì)量的行星激發(fā)的密度衰減[16]等.

2.2 HD 106315系統(tǒng)的共振捕獲

在行星系統(tǒng)演化后期通常會(huì)形成穩(wěn)定的平運(yùn)動(dòng)共振軌道構(gòu)型, 目前系外行星共振捕獲理論包括: 原行星盤驅(qū)動(dòng)遷移模型、行星子散射遷移模型和潮汐演化模型[27]. 在不受外部擾動(dòng)作用下, 行星不再產(chǎn)生劇烈的軌道遷移. 由于此時(shí)行星盤內(nèi)剩余的氣體對行星產(chǎn)生作用很小可以忽略, 若行星質(zhì)量足夠大且軌道偏心率較小, 那么兩顆行星的軌道相對位置會(huì)長期保持穩(wěn)定而形成所謂的“共振鎖定”現(xiàn)象. 描述共振軌道構(gòu)型的共振角為:θ=(p+q)λ′?pλ ?q?′.λ′和λ分別為兩顆行星的平運(yùn)動(dòng)經(jīng)度:λ=n(t ?tσ)+?,其中,?和?′分別為兩顆行星的近星點(diǎn)經(jīng)度,? ≡?+ω, ?為軌道升交點(diǎn)經(jīng)度,ω為近星點(diǎn)幅角. HD 106315系統(tǒng)當(dāng)前軌道構(gòu)型處于近似2 : 1平運(yùn)動(dòng)共振, 下文將主要圍繞HD 106315系統(tǒng)進(jìn)入2:1共振捕獲的機(jī)制進(jìn)行討論.

在行星形成初期,軌道偏心率一般都較小. Batygin[27]提出,對于同樣的初始偏心率,隨著行星質(zhì)量比m1/m2以及行星總質(zhì)量與中央恒星質(zhì)量比(m1+m2)/m0同時(shí)增大(下標(biāo)0、1、2分別表示恒星和內(nèi)、外行星, 下文同), 確定形成共振捕獲的概率也越大. 他們還發(fā)現(xiàn)當(dāng)m1/m2= 1, (m1+m2)/m0= 10?5,e1和e2均小于0.03時(shí)對應(yīng)2 : 1平運(yùn)動(dòng)共振捕獲概率為100%. 對于HD 106315行星系統(tǒng),m1/m2=mb/mc= 0.829, 接近1, 且(mb+mc)/m0=7.87×10?5, 隨著兩顆行星偏心率增大, 2:1共振捕獲概率逐漸減小, 對應(yīng)的2:1共振捕獲臨界偏心率為emax=0.08.

同時(shí), Batygin[27]還在研究木星和土星共振構(gòu)型時(shí)得出2 : 1共振捕獲概率隨行星近星點(diǎn)經(jīng)度差值??的等高線圖, 當(dāng)??=π時(shí), 行星確定形成共振捕獲的數(shù)據(jù)覆蓋范圍最小, 當(dāng)??= 0時(shí), 行星確定形成共振捕獲覆蓋的參數(shù)范圍最大. 因此, 我們在設(shè)置HD 106315行星系統(tǒng)的軌道初始參數(shù)時(shí), 兩顆行星的初始軌道近星點(diǎn)經(jīng)度均設(shè)為0.

當(dāng)行星軌道形成穩(wěn)定共振構(gòu)型之后, 其偏心率演化是由近星點(diǎn)平均進(jìn)動(dòng)速率決定的. 兩顆行星近星點(diǎn)軌道進(jìn)動(dòng)速率達(dá)到同步時(shí), 軌道偏心率也會(huì)受到激發(fā). 且在實(shí)際遷移過程中, 存在一個(gè)偏心率衰減模型, 該衰減項(xiàng)與行星間的平運(yùn)動(dòng)共振作用對偏心率的激發(fā)項(xiàng)達(dá)到平衡, 此時(shí)偏心率會(huì)停止增長. 因此通過調(diào)整偏心率衰減的數(shù)學(xué)模型可以改變共振捕獲后的偏心率平衡值, 以此得到接近觀測值的軌道遷移機(jī)制. 該結(jié)論將應(yīng)用到數(shù)值模擬的參數(shù)空間選擇中.

3 數(shù)值模擬與結(jié)果分析

3.1 MERCURY6算法

傳統(tǒng)的N-Body動(dòng)力學(xué)積分算法在進(jìn)行長時(shí)間的積分時(shí)無法避免產(chǎn)生累積的能量誤差, 且積分速度較慢. 本工作在對HD 106315行星系統(tǒng)進(jìn)行數(shù)值模擬時(shí)將采用辛算法積分器[28], 它不會(huì)產(chǎn)生長時(shí)間的能量誤差累積, 并且當(dāng)系統(tǒng)內(nèi)行星數(shù)量較少、天體質(zhì)量較集中時(shí), 積分的速度將更快. 但是, 辛算法的缺點(diǎn)是當(dāng)兩顆行星發(fā)生近距離交會(huì)時(shí), 積分結(jié)果會(huì)變得不精確. 為了解決這個(gè)問題, 我們采用了適用于N-Body問題積分的MERCURY6算法包[29], 在處理天體近距離交會(huì)問題時(shí)使用傳統(tǒng)的積分方法, 其余部分采用混合辛算法(Hybrid Symplectic).

在對HD 106315行星系統(tǒng)進(jìn)行軌道遷移數(shù)值模擬時(shí), 除了考慮中央恒星的引力作用和行星之間的攝動(dòng)作用, 還將考慮加入原行星盤耗散效應(yīng)后的行星軌道遷移動(dòng)力學(xué)模型. 由于原行星盤的耗散作用會(huì)使軌道半長徑和偏心率產(chǎn)生衰減, 我們在HD 106315系統(tǒng)遷移模型中添加了準(zhǔn)確的半長徑a和偏心率e指數(shù)衰減項(xiàng), 所以可以將哈密頓函數(shù)對應(yīng)分解為3部分:Ha、He和Hgrav(即分別由軌道半長徑a和偏心率e的衰減以及引力作用項(xiàng)表示的系統(tǒng)能量函數(shù)), 并將半長徑和偏心率的衰減項(xiàng)寫成哈密頓函數(shù)具體積分項(xiàng)的形式. 行星b軌道遷移的啟動(dòng)時(shí)刻受到行星c遷移速率的影響, 我們固定行星c初始時(shí)刻的軌道遷移速率為·ac, 假設(shè)原行星盤中的行星偏心率衰減速率滿足Lee等[16]提出的由原行星盤耗散效應(yīng)產(chǎn)生的軌道遷移偏心率衰減模型, 我們僅考慮行星c的偏心率衰減,即e的衰減速率滿足: ·e2/e2=?K|·a2/a2|. 其中衰減系數(shù)K是無量綱的正常數(shù). 鑒于·a/a和·e/e實(shí)際上是半長徑和偏心率變化的時(shí)標(biāo)(倒數(shù)),K值實(shí)際上是偏心率變化時(shí)標(biāo)和半長徑變化時(shí)標(biāo)之間的比值, 偏心率最終平衡值只與K相關(guān).

3.2 參數(shù)空間與軌道遷移結(jié)果

常見的設(shè)定軌道初值的方法有: 固定兩顆行星的初始位置, 改變行星遷移的速率;或固定行星遷移的速率, 改變兩顆行星的初始軌道位置. 這兩種選取參數(shù)空間的方法均可以使行星在一定時(shí)間內(nèi)遷移至目標(biāo)軌道位置, 且可通過空間尺度與時(shí)間之比進(jìn)行等價(jià)變換. 根據(jù)2.1節(jié)介紹的第I類軌道遷移半長徑和偏心率演化模型, 可初步判斷行星軌道遷移初始階段在原行星盤內(nèi)的相對位置, 但由于初始時(shí)刻行星的具體軌道位置與遷移時(shí)標(biāo)均未知, 本文在一定初始軌道范圍內(nèi)固定兩顆行星半長徑比值及行星c初始時(shí)刻遷移的速率, 同時(shí)設(shè)定軌道遷移時(shí)標(biāo)為105yr. 假設(shè)初始軌道半長徑滿足ac,0/ab,0= 2或接近此值, 以避免兩顆行星因初始時(shí)刻相距過近而過早地產(chǎn)生密近交會(huì), 從而使積分過快停止. 假設(shè)HD 106315c的軌道遷移速率滿足: ·ac=?5×10?5au·yr?1, 再根據(jù)2.2節(jié), 當(dāng)兩顆行星偏心率均小于0.03時(shí), 將100%形成2 : 1平運(yùn)動(dòng)共振捕獲, 隨著偏心率增大, 形成2:1共振捕獲的概率將減小.

現(xiàn)設(shè)兩顆行星軌道初始偏心率均為e= 0.01, 其余軌道參數(shù)設(shè)置如表2所示. 表中所列的幾組軌道半長徑是經(jīng)過模型線性化假設(shè)得到的初始參數(shù)組合,Pc/Pb是數(shù)值模擬得到的兩行星軌道周期比. 其中, 初始值ab,0= 0.2–0.6 au、ac,0= 0.4–1.2 au得到的軌道周期比為2 : 1. 圖1中所示為ab,0= 0.2–1.0 au、ac,0= 0.4–2.0 au時(shí)分別對應(yīng)的軌道周期比變化圖像. 在初始偏心率等于0.01的情形下,ab,00.6 au、ac,01.2 au時(shí), 可形成2:1平運(yùn)動(dòng)共振捕獲. 當(dāng)ab,0=0.8 au、ac,0=1.6 au和ab,0=1.0 au、ac,0=2.0 au時(shí),HD 106315系統(tǒng)會(huì)分別進(jìn)入3 : 2和4 : 3平運(yùn)動(dòng)共振捕獲. 所以根據(jù)共振比確定軌道遷移的初始半長徑范圍為: 0.2 auab,00.6 au、0.4 auac,01.2 au.

表2 初始軌道參數(shù)及對應(yīng)軌道周期比Table 2 Initial orbital parameters and corresponding orbital period ratio

圖1中, case1和case2在形成2 : 1共振捕獲一段時(shí)間后, 對應(yīng)的軌道周期比均出現(xiàn)下降, 這是由于兩顆行星初始時(shí)刻軌道位置比較靠近, 第I類軌道遷移速率較快, 兩顆行星周期比迅速減小從而脫離2:1平運(yùn)動(dòng)共振.

除了滿足2 : 1平運(yùn)動(dòng)共振捕獲, 還需要兩顆行星軌道半長徑能遷移到觀測值附近,即接近ab= 0.09012 au、ac= 0.1526 au. 使用表2中的參數(shù)進(jìn)行數(shù)值模擬, 得到圖2所示的軌道半長徑演化圖像. 在case1–case4中, case1和case3會(huì)因發(fā)生近距離交會(huì)而使系統(tǒng)失穩(wěn), 只有在case2中, 當(dāng)ab,0～0.4 au、ac,0～0.8 au時(shí), 兩顆行星半長徑可以同時(shí)穩(wěn)定遷移到當(dāng)前觀測數(shù)值附近(行星b和c當(dāng)前軌道半長徑分別為0.09012 au和0.1526 au),且在65000 yr附近兩顆行星軌道半長徑之差約為0.06 au. 故ab,0～0.4 au、ac,0～0.8 au可能是符合HD 106315行星系統(tǒng)軌道遷移的初始軌道位置.

圖1 不同初始軌道半長徑下得到的軌道周期比Fig.1 Orbital period ratio of planets with different initial semi-major axes

圖2 不同初始軌道半長徑下的行星軌道半長徑演化Fig.2 Evolution process of planetary semi-major axis with different initial semi-major axes

在HD 106315初始半長徑分別為ab,0= 0.4 au、ac,0= 0.8 au、eb,0=ec,0= 0的前提下, 本節(jié)考慮了原行星盤對軌道遷移產(chǎn)生的耗散效應(yīng), 通過改變行星盤內(nèi)的偏心率衰減系數(shù)可以調(diào)整軌道遷移時(shí)標(biāo)及最終軌道偏心率數(shù)值. 令兩顆行星初始偏心率均為0,設(shè)置不同偏心率衰減系數(shù)K, 得到的偏心率平衡值見表3. 為了進(jìn)一步說明偏心率衰減項(xiàng)K與偏心率最終穩(wěn)定值的關(guān)系, 我們給出偏心率隨時(shí)間的變化圖像. 圖3對應(yīng)的分別是K= 0、100、200、300時(shí)的軌道偏心率演化.K= 0即沒有添加偏心率衰減效應(yīng)時(shí),兩顆行星的偏心率在平運(yùn)動(dòng)共振作用下不斷激發(fā), 直到超出觀測值.K= 300時(shí), HD 106315c的軌道偏心率受到抑制而無法與HD 106315b形成共振穩(wěn)定, 且HD 106315b的偏心率快速增大, 系統(tǒng)快速失穩(wěn). 根據(jù)觀測值ec= 0.23、eb= 0.3可知K= 200是當(dāng)ab,0= 0.4 au、ac,0= 0.8 au、eb,0=ec,0= 0時(shí)最接近HD 106315的實(shí)際軌道遷移過程的偏心率衰減系數(shù).

表3 偏心率衰減系數(shù)及對應(yīng)的平衡值Table 3 Attenuation coefficient and corresponding equilibrium value of eccentricity

圖3 不同偏心率衰減系數(shù)下的軌道偏心率演化Fig.3 Orbital eccentricity evolution with different attenuation coefficients of eccentricity

綜上, 我們的數(shù)值模擬結(jié)果給出了HD 106315行星系統(tǒng)的初始軌道半長徑和偏心率范圍. 假設(shè)初始軌道半長徑為ab,0= 0.4 au、ac,0= 0.8 au, 兩顆行星初始偏心率均為e=0. 隨著行星在行星盤內(nèi)受粘滯力作用,HD 106315c先開始向內(nèi)遷移,在最初10000 yr, 行星HD 106315b保持在初始軌道位置不變. 隨后兩顆行星軌道逐漸接近, 在HD 106315c的引力攝動(dòng)作用以及中央恒星的引力作用下,HD 106315b開始向內(nèi)遷移,并進(jìn)入2 : 1平運(yùn)動(dòng)共振捕獲, 兩顆行星的軌道偏心率被激發(fā). 該行星系統(tǒng)軌道半長徑和2 : 1平運(yùn)動(dòng)共振對應(yīng)的共振角演化圖像見圖4.

圖4 K = 200對應(yīng)的軌道半長徑和共振角演化圖像Fig.4 Evolution of semi-major axis and resonance angle with K = 200

4 HD 106315系統(tǒng)的潮汐演化

4.1 一般長期攝動(dòng)理論修正

根據(jù)Kepler空間望遠(yuǎn)鏡觀測的系外行星數(shù)據(jù), 軌道半長徑a >0.2 au的行星軌道偏心率均值在0.3左右, 最大可接近1, 且分布比較均勻. 對于軌道半長徑a <0.2 au的系外行星, 其軌道偏心率通常較小(e <0.2)或接近0, 密近軌道的偏心率降低主要是由于主星的潮汐耗散效應(yīng)[30–32]. 行星軌道偏心率在來自中央恒星的潮汐作用下衰減, 行星公轉(zhuǎn)角動(dòng)量轉(zhuǎn)化為中央恒星的自旋角動(dòng)量, 行星偏心率衰減時(shí)標(biāo)主要由行星潮汐耗散系數(shù)(tidal dissipation factor)Qp決定.

由第3節(jié)的軌道遷移結(jié)果可知, HD 106315行星系統(tǒng)中, 較大的軌道偏心率使得系統(tǒng)共振軌道構(gòu)型不穩(wěn)定, 考慮軌道遷移停止后的偏心率衰減效應(yīng)可得到較穩(wěn)定的共振軌道構(gòu)型. 且觀測時(shí)刻系統(tǒng)內(nèi)氣體已耗散殆盡, 無法通過氣體耗散作用使系統(tǒng)偏心率有效衰減. HD 106315b的軌道半長徑為0.09012 au, 此時(shí)主要考慮外軌道HD 106315c的引力攝動(dòng)作用及中央恒星對HD 106315b的潮汐作用. 系外行星潮汐耗散系數(shù)Qp與行星組成和內(nèi)部結(jié)構(gòu)有關(guān), 在沒有精確質(zhì)量和半徑測量數(shù)據(jù)作為限定條件時(shí), 通常使用估計(jì)值. HD 106315b是一顆海王星質(zhì)量的行星, 根據(jù)Gavrilov等[33]的估計(jì), 類海王星的潮汐耗散系數(shù)Qp～105.

與經(jīng)典單行星系統(tǒng)潮汐理論對比, 雙行星系統(tǒng)的潮汐演化模型中添加了中央恒星對內(nèi)外軌道行星的潮汐力[34?38],ft1和ft2分別是內(nèi)、外軌道行星受到的來自中央恒星的潮汐力:

其中G是引力常數(shù),m0為恒星質(zhì)量,R1,2、n1,2、Q′1,2分別為內(nèi)外行星的半徑、公轉(zhuǎn)速率及修正后的潮汐耗散系數(shù).Q′ ≡3Q/(2k), 其中勒夫數(shù)(Love number)k與天體潮汐有效剛性及密度分布有關(guān), 無法精確測定, 此處將實(shí)際的行星勒夫數(shù)k參數(shù)化為3/2,使Q′=Q, 下文統(tǒng)一使用潮汐耗散系數(shù)Q.v1= ·r1、v2= ·r2, 分別為內(nèi)外兩顆行星的速度, ?1,2是內(nèi)外軌道行星的自轉(zhuǎn)速率.

HD 106315行星系統(tǒng)內(nèi)兩顆行星軌道半長徑差值僅為0.06 au, 行星間的引力攝動(dòng)作用明顯, 故此處需考慮引力攝動(dòng)和潮汐耗散效應(yīng)耦合作用下的雙行星系統(tǒng)的潮汐演化.Laskar等[39]在一般長期攝動(dòng)理論中添加了相對論和潮汐效應(yīng), 修正后的系統(tǒng)長期攝動(dòng)方程可寫為:

其中,z是行星的位置向量, 行星的運(yùn)動(dòng)位置采用復(fù)變量表示法:zk=ekei?k,ek和?k分別是第k顆行星的偏心率和近星點(diǎn)經(jīng)度,A、δA、δB均為對角矩陣,A是修正前的攝動(dòng)方程系數(shù)矩陣,Atot是添加修正項(xiàng)δA之后的復(fù)數(shù)部分系數(shù)矩陣,δAkk是δA對角線上的元素值,δAkk為相對論效應(yīng)修正項(xiàng)δA(1)kk和潮汐效應(yīng)保守項(xiàng)δA(2)kk的疊加,δB為潮汐效應(yīng)的耗散項(xiàng). 系統(tǒng)長期攝動(dòng)方程的解為:

其中,uk(t)表示任意時(shí)刻t行星k的位置,γk為考慮耗散項(xiàng)修正定義的系數(shù),gk為矩陣A的本征值,uk(0)為一般長期攝動(dòng)方程的初解, 可通過行星偏心率初值求解. 復(fù)數(shù)形式的行星偏心率為uk(t)的線性組合, 線性組合的系數(shù)來自Atot與δB矩陣運(yùn)算得到的復(fù)數(shù)矩陣.

4.2 潮汐演化理論分析與數(shù)值模擬

盡管潮汐耗散時(shí)標(biāo)相比軌道遷移時(shí)標(biāo)長很多, 但在考慮相對論和潮汐效應(yīng)的長期演化過程中, 潮汐效應(yīng)保守項(xiàng)和相對論效應(yīng)會(huì)減弱內(nèi)軌道行星因外軌道行星引力攝動(dòng)產(chǎn)生的偏心率激發(fā), 從而使行星半長徑和偏心率均會(huì)緩慢衰減[40–41]. 潮汐演化數(shù)值模擬仍然使用MERCURY6算法, 并將4.1節(jié)的(2)、(3)式添加到潮汐演化動(dòng)力學(xué)模型中. 本章潮汐演化僅考慮遷移結(jié)束后的軌道演化, 軌道半長徑初始值設(shè)在觀測值附近, 且潮汐演化時(shí)標(biāo)遠(yuǎn)大于軌道遷移時(shí)標(biāo). 故潮汐演化結(jié)果與遷移階段的初始軌道構(gòu)型并不直接相關(guān).

首先, 固定初始軌道半長徑和軌道偏心率:ab= 0.09012 au、ac= 0.1526 au、eb= 0.4、ec= 0.3, 取行星潮汐耗散系數(shù)分別為: 102、103、104、105, 通過改變行星潮汐耗散系數(shù)Q觀察不同演化時(shí)標(biāo)下兩顆行星的潮汐演化結(jié)果, 得到的演化結(jié)果如圖5所示. 該系統(tǒng)的中央恒星HD 106315年齡約為40億年, 可推測該行星系統(tǒng)剩余壽命約為50–60億年. 由圖5可知, 兩顆行星的軌道偏心率衰減時(shí)標(biāo)隨潮汐耗散系數(shù)的增大而變長,當(dāng)Q=102、103時(shí),行星的偏心率衰減曲線越過了目前的觀測值,而當(dāng)Q104時(shí),在當(dāng)前軌道位置上無法通過潮汐效應(yīng)使偏心率在系統(tǒng)剩余壽命內(nèi)衰減至觀測值eb= 0.3、ec=0.23. 根據(jù)前期的研究結(jié)果[41–43], 潮汐耗散系數(shù)與行星軌道的偏心率衰減時(shí)標(biāo)成反比, 而不影響行星軌道半徑演化的最后值. 我們的研究結(jié)果也與這些結(jié)果一致.

圖5 不同潮汐耗散系數(shù)下的軌道偏心率演化Fig.5 Orbital eccentricity evolution with different tidal dissipation factors

此外, 為了探究兩顆行星的初始偏心率差值是否會(huì)對行星潮汐演化產(chǎn)生影響, 假設(shè)行星潮汐耗散系數(shù)Q=100, 改變初始時(shí)刻兩顆行星偏心率差值?e=eb,0?ec,0, 得到的偏心率演化結(jié)果如圖6所示, 其中ec,0= 0.2, 0.1

由圖6可知, 當(dāng)?e較小時(shí), 行星偏心率振蕩幅度較小, 隨著?e增大, 兩顆行星偏心率振蕩幅度均變大, 且初始時(shí)刻HD 106315c偏心率受HD 106315b引力攝動(dòng)激發(fā)較明顯, 隨后波動(dòng)下降, 反之HD 106315b偏心率不會(huì)被HD 106315c激發(fā), 而是直接振蕩衰減, 且行星偏心率衰減時(shí)標(biāo)不會(huì)隨?e發(fā)生變化. 該研究結(jié)果可以解釋為: 外行星偏心率一定時(shí),隨著內(nèi)行星偏心率變大(?e變大), 內(nèi)行星的潮汐衰減效應(yīng)增強(qiáng), 在其引力攝動(dòng)影響下,外行星的軌道衰減也增大. 這與Mardling等[36]的研究結(jié)果一致, 導(dǎo)致外行星偏心率衰減的主要原因在于它受到內(nèi)行星潮汐衰減效應(yīng), 而不是其自身受到的潮汐效應(yīng).

圖6 引力攝動(dòng)和潮汐效應(yīng)下的行星軌道偏心率演化(理論分析)Fig.6 Orbital eccentricity evolution with perturbation and tidal dissipation (theoretical analysis)

第3節(jié)中系統(tǒng)的軌道演化結(jié)果表明, HD 106315系統(tǒng)可能經(jīng)歷了2:1共振捕獲及脫離共振的過程. 本節(jié)給出了HD 106315系統(tǒng)脫離2 : 1共振的可能物理圖像: 行星系統(tǒng)在潮汐演化過程中會(huì)經(jīng)歷偏心率和軌道半長徑的衰減, HD 106315b由于受到更強(qiáng)的潮汐效應(yīng)會(huì)進(jìn)一步向內(nèi)遷移, 而HD 106315c受到的潮汐效應(yīng)較弱, 從而使兩顆行星軌道周期比發(fā)生變化, 使系統(tǒng)脫離2:1平運(yùn)動(dòng)共振構(gòu)型. 但是潮汐演化數(shù)值模擬表明, HD 106315系統(tǒng)內(nèi)兩顆行星潮汐效應(yīng)作用太弱, 并不會(huì)使2:1共振捕獲的兩顆行星軌道周期比增大至脫離共振.

行星軌道潮汐演化數(shù)值模擬結(jié)果見圖7, 初始軌道參數(shù)同圖5. 當(dāng)潮汐耗散系數(shù)Q= 100時(shí), 潮汐效應(yīng)造成的軌道半長徑衰減會(huì)使系統(tǒng)軌道周期比逐漸增加. HD 106315系統(tǒng)內(nèi)兩顆行星Q103, 潮汐演化過程中軌道周期比幾乎沒有改變, 驗(yàn)證了中央恒星產(chǎn)生的潮汐效應(yīng)不是HD 106315系統(tǒng)當(dāng)前脫離2 : 1平運(yùn)動(dòng)共振構(gòu)型的原因. HD 106315系統(tǒng)可能是由于在原行星盤內(nèi)軌道遷移初期即進(jìn)入共振捕獲, 且平運(yùn)動(dòng)共振激發(fā)偏心率時(shí)標(biāo)遠(yuǎn)小于盤內(nèi)氣體粘滯力對偏心率的衰減時(shí)標(biāo), 導(dǎo)致HD 106315b偏心率持續(xù)受到激發(fā)的同時(shí)兩顆行星軌道距離縮短, 從而脫離2:1平運(yùn)動(dòng)共振捕獲.

圖7 引力攝動(dòng)和潮汐效應(yīng)下的行星軌道半長徑和偏心率演化(數(shù)值模擬)Fig.7 Orbital semi-major axis and eccentricity evolution with perturbation and tidal dissipation(numerical simulation)

5 總結(jié)與展望

本文研究了具有大偏心率的熱海王星系統(tǒng)HD 106315, 該行星系統(tǒng)目前處于近2 : 1平運(yùn)動(dòng)共振軌道構(gòu)型. 我們研究了在I類軌道遷移和中央恒星的潮汐效應(yīng)作用下,HD 106315到達(dá)當(dāng)前觀測位置以及脫離2 : 1共振的軌道演化機(jī)制. 本文結(jié)合理論分析和數(shù)值模擬方法反演HD 106315行星系統(tǒng)軌道遷移過程, 確定符合觀測值的軌道遷移初始參數(shù)范圍為ab,0～0.4 au、ac,0～0.8 au、eb和ec0.03. 行星軌道遷移過程中因原行星盤氣體粘滯力產(chǎn)生的偏心率衰減系數(shù)K～200. 系統(tǒng)內(nèi)兩顆行星約在65000 yr演化至目前觀測的軌道上.

通過軌道遷移數(shù)值模擬發(fā)現(xiàn), HD 106315系統(tǒng)可能經(jīng)歷了2 : 1平運(yùn)動(dòng)共振捕獲和脫離的過程. 為了解釋行星系統(tǒng)脫離2 : 1平運(yùn)動(dòng)共振構(gòu)型的演化機(jī)制, 我們討論了該熱海王星系統(tǒng)的潮汐演化. 對比潮汐演化理論分析和數(shù)值模擬結(jié)果, 發(fā)現(xiàn)當(dāng)行星潮汐耗散系數(shù)Q= 100時(shí), 潮汐效應(yīng)造成的軌道半長徑衰減會(huì)使系統(tǒng)軌道周期比發(fā)生變化, 可能是系統(tǒng)脫離共振構(gòu)型的原因. HD 106315系統(tǒng)內(nèi)兩顆行星Q103, 來自中央恒星的潮汐效應(yīng)并不會(huì)使行星系統(tǒng)產(chǎn)生明顯的偏心率和軌道半長徑衰減, 中央恒星產(chǎn)生的潮汐效應(yīng)不是HD 106315系統(tǒng)當(dāng)前脫離2 : 1平運(yùn)動(dòng)共振構(gòu)型的原因. 由于Crossfield等[6]根據(jù)恒星HD 106315的視向速度數(shù)據(jù)推測該系統(tǒng)內(nèi)可能存在第3顆行星d, 且該行星質(zhì)量45M⊕, 而較大質(zhì)量的行星d的存在可能使內(nèi)部軌道的行星b和c進(jìn)入或脫離特定的共振構(gòu)型[44].

在對HD 106315行星系統(tǒng)進(jìn)行初始半長軸篩選時(shí)發(fā)現(xiàn),ab,00.8 au、ac,01.6 au時(shí)對應(yīng)的行星軌道遷移有形成3 : 2平運(yùn)動(dòng)共振和4 : 3平運(yùn)動(dòng)共振的趨勢, 這是由于隨著相鄰軌道上的行星距離縮小導(dǎo)致的軌道周期比下降. 但是潮汐演化過程中, 軌道周期比卻呈上升趨勢. 所以在未來的工作中, 我們將繼續(xù)關(guān)注在行星軌道遷移和潮汐共同作用下, 熱海王星行星系統(tǒng)由2:1平運(yùn)動(dòng)共振轉(zhuǎn)化為其他共振的關(guān)鍵影響因素和限制條件以及大偏心率密近軌道上的熱海王星系統(tǒng)的穩(wěn)定性和潮汐演化時(shí)標(biāo). 同時(shí), 密近軌道上的熱海王星也是系外行星大氣的重要研究目標(biāo), 隨著系外行星大氣觀測方法和演化理論的發(fā)展[45], 熱海王星大氣特性研究也將進(jìn)一步促進(jìn)系外行星的形成和演化理論的發(fā)展.