黃秀敏 季江徽 董 瑤

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

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

作為Kepler空間探測(cè)任務(wù)的延續(xù), K2在2014年至2018年3月完成了16期觀測(cè)工作并發(fā)現(xiàn)數(shù)百顆系外行星[4–5]. 而本文的研究對(duì)象熱海王星系統(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)兩顆行星的軌道觀測(cè)數(shù)據(jù), 其中,M為行星質(zhì)量,M⊕代表地球質(zhì)量,Rp為行星半徑,R⊕代表地球半徑,R?為中央恒星半徑,P為行星軌道周期,i和ω分別為軌道傾角和近星點(diǎn)幅角. HD 106315b和HD 106315c在密近軌道上(軌道半長(zhǎng)徑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é)合視向速度測(cè)量結(jié)果, 這兩顆行星的質(zhì)量、內(nèi)部結(jié)構(gòu)及大氣特征值得關(guān)注. 其他研究表明, 該行星系統(tǒng)在更長(zhǎng)周期的軌道上還可能存在第3顆行星[6].

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

通過對(duì)Kepler多年觀測(cè)的數(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)對(duì)應(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體共振(也稱為L(zhǎng)aplace共振)[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)遷移. 對(duì)于HD 106315行星系統(tǒng), 需要判斷行星質(zhì)量是否滿足產(chǎn)生間隙環(huán)的條件, 進(jìn)而影響后續(xù)的軌道遷移過程. 同時(shí), HD 106315行星系統(tǒng)構(gòu)型更緊湊, 恒星引力和行星間的攝動(dòng)作用對(duì)行星共振遷移的影響也會(huì)更顯著.

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

綜上所述, 本文工作將圍繞熱海王星系統(tǒng)HD 106315近平運(yùn)動(dòng)共振構(gòu)型的形成機(jī)制和軌道遷移模型展開. 第2節(jié)簡(jiǎn)要介紹HD 106315系統(tǒng)的軌道演化動(dòng)力學(xué)模型, 包括I類軌道遷移和共振捕獲. 第3節(jié)給出本文使用的參數(shù)空間和數(shù)值模擬方法, 并據(jù)此得到觀測(cè)值對(duì)應(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)對(duì)軌道共振構(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)快速吸積氣體生長(zhǎng)為氣態(tài)巨行星. Hansen等[19]和Chiang等[20]提出了類海王星行星的原位吸積模型, 行星在當(dāng)前軌道附近直接通過塵埃吸積和行星子(planetesimals)并合形成. 他們認(rèn)為, 目前觀測(cè)到的超級(jí)地球和類海王星行星形成于氣體耗散殆盡的過渡原行星盤中, 一般不會(huì)經(jīng)歷劇烈的角動(dòng)量交換過程, 從而易保持行星形成時(shí)刻的近圓軌道.

HD 106315系統(tǒng)內(nèi)兩顆行星均處于短周期密近軌道, 但是軌道偏心率較大. 這意味著該行星系統(tǒng)更可能通過軌道遷移機(jī)制形成, 且軌道遷移過程中可能存在不穩(wěn)定的軌道偏心率激發(fā)和衰減. Terquem等[21]通過對(duì)質(zhì)量為0.1～1M⊕的原行星軌道演化進(jìn)行數(shù)值模擬, 發(fā)現(xiàn)它們可以并合形成海王星質(zhì)量的行星, 且這些類海王星行星相伴出現(xiàn), 形成多行星系統(tǒng). Cresswell等[22]指出, 對(duì)于無法形成原行星盤環(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)而減慢行星遷移速率.

對(duì)于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í)刻兩顆行星共面.

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

在行星遷移的初始時(shí)刻, 通常假設(shè)兩顆行星位于相距較遠(yuǎn)的兩個(gè)軌道上, 以使行星軌道遷移的時(shí)標(biāo)盡量小于行星盤星云物質(zhì)完全耗散所需的時(shí)間. 一般可以假設(shè)行星b和行星c的初始軌道半長(zhǎng)徑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)高, 此處?是在行星軌道半長(zhǎng)徑為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)剩余的氣體對(duì)行星產(chǎn)生作用很小可以忽略, 若行星質(zhì)量足夠大且軌道偏心率較小, 那么兩顆行星的軌道相對(duì)位置會(huì)長(zhǎng)期保持穩(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]提出,對(duì)于同樣的初始偏心率,隨著行星質(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í)對(duì)應(yīng)2 : 1平運(yùn)動(dòng)共振捕獲概率為100%. 對(duì)于HD 106315行星系統(tǒng),m1/m2=mb/mc= 0.829, 接近1, 且(mb+mc)/m0=7.87×10?5, 隨著兩顆行星偏心率增大, 2:1共振捕獲概率逐漸減小, 對(duì)應(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ù)范圍最大. 因此, 我們?cè)谠O(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)共振作用對(duì)偏心率的激發(fā)項(xiàng)達(dá)到平衡, 此時(shí)偏心率會(huì)停止增長(zhǎng). 因此通過調(diào)整偏心率衰減的數(shù)學(xué)模型可以改變共振捕獲后的偏心率平衡值, 以此得到接近觀測(cè)值的軌道遷移機(jī)制. 該結(jié)論將應(yīng)用到數(shù)值模擬的參數(shù)空間選擇中.

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

3.1 MERCURY6算法

傳統(tǒng)的N-Body動(dòng)力學(xué)積分算法在進(jìn)行長(zhǎng)時(shí)間的積分時(shí)無法避免產(chǎn)生累積的能量誤差, 且積分速度較慢. 本工作在對(duì)HD 106315行星系統(tǒng)進(jìn)行數(shù)值模擬時(shí)將采用辛算法積分器[28], 它不會(huì)產(chǎn)生長(zhǎng)時(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).

在對(duì)HD 106315行星系統(tǒng)進(jìn)行軌道遷移數(shù)值模擬時(shí), 除了考慮中央恒星的引力作用和行星之間的攝動(dòng)作用, 還將考慮加入原行星盤耗散效應(yīng)后的行星軌道遷移動(dòng)力學(xué)模型. 由于原行星盤的耗散作用會(huì)使軌道半長(zhǎng)徑和偏心率產(chǎn)生衰減, 我們?cè)贖D 106315系統(tǒng)遷移模型中添加了準(zhǔn)確的半長(zhǎng)徑a和偏心率e指數(shù)衰減項(xiàng), 所以可以將哈密頓函數(shù)對(duì)應(yīng)分解為3部分:Ha、He和Hgrav(即分別由軌道半長(zhǎng)徑a和偏心率e的衰減以及引力作用項(xiàng)表示的系統(tǒng)能量函數(shù)), 并將半長(zhǎng)徑和偏心率的衰減項(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í)際上是半長(zhǎng)徑和偏心率變化的時(shí)標(biāo)(倒數(shù)),K值實(shí)際上是偏心率變化時(shí)標(biāo)和半長(zhǎng)徑變化時(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類軌道遷移半長(zhǎng)徑和偏心率演化模型, 可初步判斷行星軌道遷移初始階段在原行星盤內(nèi)的相對(duì)位置, 但由于初始時(shí)刻行星的具體軌道位置與遷移時(shí)標(biāo)均未知, 本文在一定初始軌道范圍內(nèi)固定兩顆行星半長(zhǎng)徑比值及行星c初始時(shí)刻遷移的速率, 同時(shí)設(shè)定軌道遷移時(shí)標(biāo)為105yr. 假設(shè)初始軌道半長(zhǎng)徑滿足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所示. 表中所列的幾組軌道半長(zhǎng)徑是經(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í)分別對(duì)應(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ù)共振比確定軌道遷移的初始半長(zhǎng)徑范圍為: 0.2 auab,00.6 au、0.4 auac,01.2 au.

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

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

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

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

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

在HD 106315初始半長(zhǎng)徑分別為ab,0= 0.4 au、ac,0= 0.8 au、eb,0=ec,0= 0的前提下, 本節(jié)考慮了原行星盤對(duì)軌道遷移產(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對(duì)應(yīng)的分別是K= 0、100、200、300時(shí)的軌道偏心率演化.K= 0即沒有添加偏心率衰減效應(yīng)時(shí),兩顆行星的偏心率在平運(yùn)動(dòng)共振作用下不斷激發(fā), 直到超出觀測(cè)值.K= 300時(shí), HD 106315c的軌道偏心率受到抑制而無法與HD 106315b形成共振穩(wěn)定, 且HD 106315b的偏心率快速增大, 系統(tǒng)快速失穩(wěn). 根據(jù)觀測(cè)值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ù)及對(duì)應(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)的初始軌道半長(zhǎng)徑和偏心率范圍. 假設(shè)初始軌道半長(zhǎng)徑為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)軌道半長(zhǎng)徑和2 : 1平運(yùn)動(dòng)共振對(duì)應(yīng)的共振角演化圖像見圖4.

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

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

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

根據(jù)Kepler空間望遠(yuǎn)鏡觀測(cè)的系外行星數(shù)據(jù), 軌道半長(zhǎng)徑a >0.2 au的行星軌道偏心率均值在0.3左右, 最大可接近1, 且分布比較均勻. 對(duì)于軌道半長(zhǎng)徑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)型. 且觀測(cè)時(shí)刻系統(tǒng)內(nèi)氣體已耗散殆盡, 無法通過氣體耗散作用使系統(tǒng)偏心率有效衰減. HD 106315b的軌道半長(zhǎng)徑為0.09012 au, 此時(shí)主要考慮外軌道HD 106315c的引力攝動(dòng)作用及中央恒星對(duì)HD 106315b的潮汐作用. 系外行星潮汐耗散系數(shù)Qp與行星組成和內(nèi)部結(jié)構(gòu)有關(guān), 在沒有精確質(zhì)量和半徑測(cè)量數(shù)據(jù)作為限定條件時(shí), 通常使用估計(jì)值. HD 106315b是一顆海王星質(zhì)量的行星, 根據(jù)Gavrilov等[33]的估計(jì), 類海王星的潮汐耗散系數(shù)Qp～105.

與經(jīng)典單行星系統(tǒng)潮汐理論對(duì)比, 雙行星系統(tǒng)的潮汐演化模型中添加了中央恒星對(duì)內(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), 無法精確測(cè)定, 此處將實(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)兩顆行星軌道半長(zhǎng)徑差值僅為0.06 au, 行星間的引力攝動(dòng)作用明顯, 故此處需考慮引力攝動(dòng)和潮汐耗散效應(yīng)耦合作用下的雙行星系統(tǒng)的潮汐演化.Laskar等[39]在一般長(zhǎng)期攝動(dòng)理論中添加了相對(duì)論和潮汐效應(yīng), 修正后的系統(tǒng)長(zhǎng)期攝動(dòng)方程可寫為:

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

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

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

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

首先, 固定初始軌道半長(zhǎng)徑和軌道偏心率: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億年, 可推測(cè)該行星系統(tǒng)剩余壽命約為50–60億年. 由圖5可知, 兩顆行星的軌道偏心率衰減時(shí)標(biāo)隨潮汐耗散系數(shù)的增大而變長(zhǎng),當(dāng)Q=102、103時(shí),行星的偏心率衰減曲線越過了目前的觀測(cè)值,而當(dāng)Q104時(shí),在當(dāng)前軌道位置上無法通過潮汐效應(yīng)使偏心率在系統(tǒng)剩余壽命內(nèi)衰減至觀測(cè)值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

此外, 為了探究?jī)深w行星的初始偏心率差值是否會(huì)對(duì)行星潮汐演化產(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)歷偏心率和軌道半長(zhǎ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)造成的軌道半長(zhǎ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)氣體粘滯力對(duì)偏心率的衰減時(shí)標(biāo), 導(dǎo)致HD 106315b偏心率持續(xù)受到激發(fā)的同時(shí)兩顆行星軌道距離縮短, 從而脫離2:1平運(yùn)動(dòng)共振捕獲.

圖7 引力攝動(dòng)和潮汐效應(yīng)下的行星軌道半長(zhǎ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)前觀測(cè)位置以及脫離2 : 1共振的軌道演化機(jī)制. 本文結(jié)合理論分析和數(shù)值模擬方法反演HD 106315行星系統(tǒng)軌道遷移過程, 確定符合觀測(cè)值的軌道遷移初始參數(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演化至目前觀測(cè)的軌道上.

通過軌道遷移數(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)的潮汐演化. 對(duì)比潮汐演化理論分析和數(shù)值模擬結(jié)果, 發(fā)現(xiàn)當(dāng)行星潮汐耗散系數(shù)Q= 100時(shí), 潮汐效應(yīng)造成的軌道半長(zhǎng)徑衰減會(huì)使系統(tǒng)軌道周期比發(fā)生變化, 可能是系統(tǒng)脫離共振構(gòu)型的原因. HD 106315系統(tǒng)內(nèi)兩顆行星Q103, 來自中央恒星的潮汐效應(yīng)并不會(huì)使行星系統(tǒng)產(chǎn)生明顯的偏心率和軌道半長(zhǎng)徑衰減, 中央恒星產(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ù)推測(cè)該系統(tǒng)內(nèi)可能存在第3顆行星d, 且該行星質(zhì)量45M⊕, 而較大質(zhì)量的行星d的存在可能使內(nèi)部軌道的行星b和c進(jìn)入或脫離特定的共振構(gòu)型[44].

在對(duì)HD 106315行星系統(tǒng)進(jìn)行初始半長(zhǎng)軸篩選時(shí)發(fā)現(xiàn),ab,00.8 au、ac,01.6 au時(shí)對(duì)應(yīng)的行星軌道遷移有形成3 : 2平運(yùn)動(dòng)共振和4 : 3平運(yùn)動(dòng)共振的趨勢(shì), 這是由于隨著相鄰軌道上的行星距離縮小導(dǎo)致的軌道周期比下降. 但是潮汐演化過程中, 軌道周期比卻呈上升趨勢(shì). 所以在未來的工作中, 我們將繼續(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), 隨著系外行星大氣觀測(cè)方法和演化理論的發(fā)展[45], 熱海王星大氣特性研究也將進(jìn)一步促進(jìn)系外行星的形成和演化理論的發(fā)展.