基于天體測量法探測系外行星的仿真模擬與軌道求解算法?

黃平輝 季江徽

(1中國科學(xué)院紫金山天文臺 南京 210008) (2中國科學(xué)院行星科學(xué)重點實驗室 南京 210008) (3中國科學(xué)院大學(xué) 北京 100049)

基于天體測量法探測系外行星的仿真模擬與軌道求解算法?

黃平輝1,2,3?季江徽1,2?

(1中國科學(xué)院紫金山天文臺 南京 210008) (2中國科學(xué)院行星科學(xué)重點實驗室 南京 210008) (3中國科學(xué)院大學(xué) 北京 100049)

天體測量法是探測系外行星一種有效的手段,具有其他探測方法所不具備的獨特優(yōu)勢,如可以探測到行星的3維軌道參數(shù)及確定行星的質(zhì)量等,這些將為現(xiàn)有的行星樣本提供重要補充.高精度的空間天體測量衛(wèi)星Gaia(G lobal Astrometric Interferometer for Astrophysics)在2013年成功發(fā)射升空,可以預(yù)期未來基于其探測數(shù)據(jù),將可能發(fā)現(xiàn)大量的長周期類木行星.針對αCentauri A、HD 62509和GJ 876系統(tǒng),根據(jù)Gaia單次天體測量精度生成對應(yīng)的仿真數(shù)據(jù),使用Lomb-Scargle周期譜分析行星軌道周期信號,并采用M arkov Chain Monte Carlo(MCMC)算法對行星系統(tǒng)進行軌道反演,得到的結(jié)果與行星的初始參數(shù)相吻合.

天體測量學(xué),行星與衛(wèi)星:演化,天體測量衛(wèi)星:Hipparcos,Gaia,方法:Lomb-Scargle周期譜,MCMC

1 引言

1995年,Mayor等人通過地面視向速度觀測發(fā)現(xiàn)第1顆圍繞主序星公轉(zhuǎn)的系外行星——飛馬座51 b[1],迄今已有3000多顆系外行星通過不同的方法被發(fā)現(xiàn)(見圖1).目前探測系外行星主要有以下方法:視向速度法(Radial Velocity),凌星法(Transit), TTV(Transit Tim ing Variation)法,直接成像法(Direct Imaging),微引力透鏡法(M icrolensing),脈沖星計時法(Pulsar Tim ing)和天體測量法(Astrometry).由圖1可知,目前發(fā)現(xiàn)大多數(shù)系外行星通過視向速度法和凌星法被發(fā)現(xiàn).在早期,系外行星主要通過地面上的視向速度觀測發(fā)現(xiàn),而在2009年美國航空航天局(NASA)發(fā)射的Kelper衛(wèi)星通過凌星法發(fā)現(xiàn)了4600多顆有凌星信號的系外行星候選體,其中大部分據(jù)信是行星.2010年,M uterspaugh等人基于天體測量法利用Palomar天文臺的PHASE(Palomar High-precision Astrometric Search for Exoplanet System s)巡天計劃發(fā)現(xiàn)了一顆疑似質(zhì)量為1.5MJupiter的行星HD 176051 b[2].2013年,歐空局(ESA)成功發(fā)射了新一代高精度天體測量衛(wèi)星Gaia,這將會極大推動系外行星的發(fā)現(xiàn),進而揭示行星系統(tǒng)的多樣性與復(fù)雜性.

圖1 已被發(fā)現(xiàn)的系外行星的半長徑與行星質(zhì)量.橫坐標(biāo)為行星的軌道半長徑,縱坐標(biāo)為行星質(zhì)量,太陽系內(nèi)的行星由黑色六邊形表示,數(shù)據(jù)來源于exop lanet.eu.Fig.1 The sem i-m a jor axes and m asses of the detected exop lanets.The x-coord inate ind icates the sem i-m a jor ax is o f p lanets,and the y-coord inate ind icates the p lanetary m ass.The p lanets of so lar system are show n as b lack hexagons.The data cam e from exop lanet.eu.

與視向速度法和凌星法比較,天體測量法有如下特點:(1)可以給出行星的軌道傾角i,解除行星質(zhì)量與軌道傾角的簡并,從而導(dǎo)出行星的質(zhì)量,結(jié)合視向速度觀測可求解出行星的全部動力學(xué)參數(shù),從而建立行星的3維軌道[3];(2)視向速度法和凌星法更容易發(fā)現(xiàn)偏心率大的行星系統(tǒng),而天體測量法則容易觀測到偏心率小的行星系統(tǒng)[4];(3)視向速度法和凌星法易于探測短周期和edge on(傾角i=90°)的行星系統(tǒng),而天體測量法則易于探測到長周期和face on(傾角i=0°)的行星系統(tǒng),可豐富目前的行星樣本;(4)視向速度法要求觀測恒星是中晚型、低自轉(zhuǎn)速度、色球活動性弱的恒星[2],而天體測量法對恒星類型、恒星自轉(zhuǎn)速度要求較低、受恒星表面黑子、耀斑、非徑向?qū)α鞯挠绊戄^視向速度法要小,由此更容易探測圍繞小質(zhì)量恒星公轉(zhuǎn)的行星[5];(5)天體測量法廣泛應(yīng)用于雙星觀測,因此易于從雙星觀測中發(fā)現(xiàn)系外行星[2].

天體測量學(xué)是古老的天文學(xué)分支,是關(guān)于位置的天文學(xué),通過測量目標(biāo)天體在天球上的位移與速度,給出天體的坐標(biāo)、自行和視差等信息.地面上的天體測量歷史悠久,最早完善的星表可追溯到第谷(Tycho)時期.1989年,歐空局發(fā)射了Hipparcos(High precision parallax collecting satellite),首次在太空中以高精度的天體測量對太陽附近的恒星進行觀測,并給出了高精度的Hipparcos星表[6].作為Hipparcos的繼任者,經(jīng)過數(shù)據(jù)優(yōu)化處理后,Gaia將有能力最終以最優(yōu)至10μas的精度進行天體測量觀測[7]. Gaia的核心科學(xué)目標(biāo)是測量銀河系約10億顆恒星(銀河系恒星數(shù)的百分之一)的自行、視差、視向速度和相關(guān)的恒星參數(shù),建立3維的銀河系模型,同時Gaia有能力在5 yr的工作壽命內(nèi)探測約20000顆的長周期類木行星[7].而國內(nèi)正積極推進的天體測量空間衛(wèi)星STEP(Search for Terrestrial Extra-solar Planet),期望以最優(yōu)至1μas的窄角天體測量精度來探測發(fā)現(xiàn)距離太陽15 pc內(nèi)的類地行星.

本文第2節(jié)介紹天體測量法探測系外行星的原理；第3節(jié)介紹尋找周期信號及進行軌道擬合中所采用的方法；第4節(jié)將針對3個具體的算例(αCentauri A、HD 62509和GJ 876),開展行星軌道參數(shù)反演研究；第5節(jié)將對天體測量法探測系外行星中的問題進行討論；第6節(jié)為總結(jié)與展望.

2 天體測量法探測系外行星

2.1 天體測量法探測系外行星原理

天體測量可給出恒星的5個天體測量參量:赤經(jīng)(α),赤緯(δ),赤經(jīng)自行(μ?α),赤緯自行(μδ)和周年視差(?).當(dāng)恒星存在行星圍繞其公轉(zhuǎn),行星會對恒星造成額外的一個小幅度的周期性橢圓擺動.因此可通過測量恒星在天球上的位置變化來探測系外行星,在恒星的運動中減去恒星自行與周年視差,然后通過對殘差進行細致的傅立葉分析后,確認最強的周期信號為行星引力擾動造成,最后再對行星進行軌道擬合.

行星對恒星造成的天體測量信號α(Astrom etric Signature)[8]:

a為行星的軌道半長徑,d為恒星離太陽的距離,M?為恒星的質(zhì)量,Mp為行星的質(zhì)量.木星質(zhì)量行星在1 pc距離、5 au處圍繞太陽質(zhì)量恒星公轉(zhuǎn)造成的天體測量信號為5mas;地球質(zhì)量行星在1 pc距離、1 au處圍繞太陽質(zhì)量恒星公轉(zhuǎn)造成的天體測量信號為3μas.

由圖2可知,目前探測到的大多數(shù)系外行星的天體測量信號都小于1mas,而通過凌星法發(fā)現(xiàn)的系外行星的天體測量信號大部分小于1μas[7].洋紅色框表示Gaia在5 yr的工作期限內(nèi)可以觀測到軌道周期大于0.2 yr小于6 yr、天體測量信號大于10μas的行星[7]；靛藍色框表示STEP預(yù)計能觀測到軌道周期大于0.2 yr小于5 yr、天體測量信號大于1 μas的行星.

圖2 已探測到的系外行星的天體測量信號,橫坐標(biāo)為行星的軌道周期,縱坐標(biāo)為行星的天體測量信號,點的大小正比于行星的質(zhì)量[7].HST-FGS(F ine G u idance Sensor)、G aia和STEP的探測精度用點虛線標(biāo)注.Fig.2 The astrom etric signatu re o f detected exop lanets.The x-coord inate ind icates the period o f p lanets,and the y-coord inate ind icates the astrom etric signatu re.The size o f p oints is p rop ortiona l to the p lanetary m ass[7].The detective ab ilities of HST-FGS(Fine Gu idance Sensor),Gaia,and STEP are p resented by the dot-dashed lines.

系外行星圍繞恒星運動的軌道參數(shù)定義見圖3：軌道半長徑a,偏心率e,軌道傾角i,升交點經(jīng)度?,近星點角距ω和過近星點時刻τ.選取天球切平面作為系外行星軌道的參考平面,天球切平面與軌道平面的夾角定義為軌道傾角i(0°–90°),升交點經(jīng)度?為兩平面交線北點的位置角,過近星點時刻τ由真近點角ν(t)表示[8].

圖3 系外行星軌道參數(shù)的定義[8]Fig.3 The definition o f orbita l elem en ts o f an exop lanet[8]

將圖3軌道根數(shù)(a,e,i,?,ω,τ)和軌道周期P、行星質(zhì)量mp轉(zhuǎn)化為Thiele-Innes根數(shù)[9]:

在(2)式中A、B、F、G為Thiele-Innes根數(shù),E(t)為行星的偏近點角.

對于存在n顆行星的系統(tǒng),如果行星之間的相互引力作用在觀測時間內(nèi)可以忽略(線性化),則恒星在天球切平面上赤經(jīng)、赤緯方向的變化分別為

Πα,t和Πδ,t分別為赤經(jīng)、赤緯方向上的視差因子,˙μ?α和˙μδ為恒星自行的加速項.對于絕大多數(shù)恒星,在短時間的觀測時間段內(nèi),自行的加速項是可以忽略的.因此采取上述的線性化模型,共需要求解7n+5個參數(shù),含5個天體測量參數(shù)(α,δ,μ?α,μδ,?)和n顆行星各7個參數(shù)(P,e,i,?,ω,τ,mp).

2.2 天體測量法與視向速度法結(jié)合

對于已由視向速度法發(fā)現(xiàn)的系外行星系統(tǒng),可通過下式來聯(lián)合視向速度數(shù)據(jù)和天體測量數(shù)據(jù)進行同時擬合[2]:

其中K為行星運動對恒星引起的視向速度半幅度(Radial Velocity Sem i-Amp litude),定義如下[10]:

例如木星以周期為11.86 yr圍繞太陽公轉(zhuǎn),對太陽造成的視向速度半幅度為12.5ms?1.

3 算法

3.1 Lom b-Scarg le周期譜

視向速度數(shù)據(jù)和天體測量數(shù)據(jù)是時序信號數(shù)據(jù),對于非均勻采樣的時序信號,可以使用改進的快速傅立葉變換算法Lomb-Scargle周期譜,找出周期信號[11?12]

tj為時間,Xj為時序信號,PX(ω)為Lomb-Scargle周期譜功率,τ由下式計算所得:

值得注意的是,Lomb-Scargle周期譜與使用正弦函數(shù)作最小二乘擬合得到的結(jié)果并無差別,同時PX(ω)并不代表真實的正弦信號擬合幅度,而是代表信號的統(tǒng)計顯著性[12].

由于恒星自轉(zhuǎn)周期、恒星活動性周期、行星間相互攝動、觀測采樣方案的周期性及行星軌道的偏心率,甚至隨機噪聲也會引起周期譜中異常的信號峰(alias),需要對周期譜中的alias進行辨認.Lomb-Scargle周期譜的FAP(False A larm Probability)功率值z0可由下式計算[12]:

p0為出現(xiàn)假信號的可能性,即FAP,N為選取的頻率計算數(shù)目.

3.2 M arkov Chain M on te Car lo(M CM C)

從觀測數(shù)據(jù)中確認行星周期信號后,需要對行星進行軌道擬合,一般使用LM (Levenberg-Marquardt)算法[13]或MCMC(Markov Chain M onte Carlo)算法[14?17]進行軌道擬合,求出行星的軌道參數(shù).LM算法特點是計算速度快,但容易得到局部最優(yōu)解,而MCMC算法能得到全局最優(yōu)解,但耗費時間長,而本文的所舉算例均使用MCMC算法進行軌道擬合.MCMC算法的原理可簡單表示為

其中

χ為歸一化的行星過近星點時刻,M1表示恒星周圍存在行星,p(M1|I)為M1的先驗概率分布,p(D|M1,I)為M1的似然函數(shù),p(D|I)為全局似然函數(shù).p(P,α,e,i,?,ω,χ|M1,I)為求解參數(shù)的先驗概率之乘積,積分的上下限由所求解的參數(shù)范圍確定.對于信噪比高的情況,可以使用均勻分布作為求解參數(shù)的先驗分布.而對于信噪比低的情況,則使用Jeffrey分布作為求解參數(shù)的先驗分布[15].

4 算例

為了說明擬合方法,本文中我們選取3個算例:αCentauri A、HD 62509和GJ 876.αCentauri A目前并沒有發(fā)現(xiàn)有系外行星圍繞,HD 62509則伴有1顆木星大小的行星,GJ 876則發(fā)現(xiàn)存在4顆行星.由文獻[7]的表2可知Gaia在未來的5 yr任務(wù)中有可能分別對它們進行107次、52次和52次的觀測(其單次觀測精度為34.2μas).對于Hipparcos和Gaia這類廣角巡天天體測量衛(wèi)星,對某一天區(qū)的觀測依賴于衛(wèi)星本身的姿態(tài)與軌道[18].在Gaia正式釋放數(shù)據(jù)之前,其軌道與姿態(tài)及對某一天區(qū)的觀測方案是不可知的,因此模擬數(shù)據(jù)中暫不考慮這些因素所帶來的影響,而觀測時刻是均勻分布隨機生成的.對于窄角天體測量觀測,通常需要有背景星定標(biāo),但在本文中我們不考慮這些數(shù)據(jù)的預(yù)處理過程,同時本文作為行星軌道擬合的示例,采用的恒星的視差與自行數(shù)據(jù)在擬合過程中視為恒值,不參與擬合,也不考慮星表本身所帶來的誤差.因此對于n顆行星的系統(tǒng),這里只擬合7n個參數(shù).

4.1 αCen tau ri A

αCentauri是距離太陽最近的雙星,距離僅為1.34 pc,同時αCentauri還與比鄰星(Proxima Centauri)組成三合星系統(tǒng).2012年Dumusque等人通過視向速度發(fā)現(xiàn)αCentauri B存在一顆地球大小、周期為3.2357 d的行星[19],但在2016年Ra jpau l等人卻證實α Centauri B b是不存在的,之前的行星信號是由于視向速度觀測的窗函數(shù)所造成[20],因此在這里我們不選取αCentauri B b作為例子.

αCentauri A的光譜型與太陽均為G2V型[21],質(zhì)量為1.1 M⊙[22],我們模擬產(chǎn)生了一顆海王星質(zhì)量的行星b在3 au處的橢圓軌道圍繞αCentauri A公轉(zhuǎn)(其他軌道根數(shù)預(yù)先假定或隨機產(chǎn)生),加入Gaia單次天體測量觀測的精度34.2μas的正態(tài)分布誤差[7],生成107組天體測量的模擬數(shù)據(jù).因為αCentauri雙星的半長徑為35.6 au,軌道周期為79.91 yr[19],對于5 yr的觀測來說影響較小,因此模擬數(shù)據(jù)中沒有考慮另一顆恒星引力作用對模擬數(shù)據(jù)的影響.

由模擬數(shù)據(jù)得到圖4的Lomb-Scargle周期譜,可以看到存在一個周期約為1750 d的周期信號,信號的FAP低于0.01%,因此可以認為是真實的周期信號,而并非由于隨機噪聲引起.我們對這個信號進行MCMC模擬,得到表1中模擬行星b的擬合結(jié)果及擬合誤差.值得注意的是,行星b的周期真值約1810 d,擬合值約1769 d,與Lomb-Scargle周期譜給出的值略有不同,這是由于Lomb-Scargle周期譜的低頻率部分的分辨率不足和信噪比不高所造成.單獨的天體測量觀測存在簡并性,升交點經(jīng)度和近星點角距各有±180°的不確定性,因此表1給出的是近星點經(jīng)度(longitude of pericentre,?=?+ω).同時擬合的是開普勒運動橢圓,過近星點時刻可由過特定時刻的平近點角M0表示,采取的時刻是JD 2457023.5,即2015年1月1日0時0分(UTC),圖5為歸一化天體測量擬合曲線及擬合殘差.

表1 αCentauri A b的擬合參數(shù)Tab le 1 The fitting p aram eters o fαCen tau ri A b

圖4 αCentau ri A的模擬天體測量數(shù)據(jù)的Lom b-Scargle周期譜,其中存在一個1750 d左右的周期信號.虛線表示1%的FAP,點虛線表示0.1%的FAP,點線表示0.01%的FAP.Fig.4 The Lom b-Sarg le p eriodogram o fm ock astrom etric data o fαCen tau ri A.A 1750-day p eriod ex ists in the p eriodogram.The 1%,0.1%,and 0.01%o f FAP are ind icated by the dashed line,dot-dashed line, and dotted line,respectively.

圖5 左:αCentau ri A的歸一化天體測量擬合曲線和擬合行星b后的殘差;右:αCentau ri A的天體測量法擬合 曲線Fig.5 Left:The norm a lized astrom etric cu rve ofαCen tau ri A and the residua l a fter fitting p lanet b; R ight:The astrom etric cu rve o fαCen tau ri A

4.2 HD 62509

HD 62509 b是一顆在2006年由Reffert等人通過視向速度法發(fā)現(xiàn)、圍繞K 0III型主星HD 62509的系外行星[23].HD 62509 b的周期約590 d,半長徑1.69 au,最小行星質(zhì)量Mpsin i為2.9 MJupiter.生成模擬的天體測量數(shù)據(jù),并結(jié)合Lick的視向速度實測數(shù)據(jù)[23],計算Lomb-Scargle周期譜,找到最強的周期信號,對周期信號進行軌道擬合.

由圖6和圖7的Lomb-Scargle周期譜可以看到Lick視向速度數(shù)據(jù)和模擬的天體測量數(shù)據(jù)中存在一個約590 d的周期信號,不同F(xiàn)AP對應(yīng)的功率在圖中表示.對HD 62509的Lick視向速度數(shù)據(jù),結(jié)合模擬的天體測量數(shù)據(jù)通過(4)式,使用MCMC算法進行擬合,得到表2中擬合結(jié)果及擬合誤差.由于結(jié)合了視向速度數(shù)據(jù),天體測量觀測中的升交點經(jīng)度?和近星點角距ω可以唯一確定,不再存在±180°的不確定性.近星點時刻用特定時刻的平近點角M0表示,采取時刻為Lick首次視向速度觀測時刻JD 2451808.39.圖8和圖9分別為歸一化視向速度法擬合曲線及擬合殘差、歸一化天體測量法擬合曲線及擬合殘差.

表2 HD 62509 b的擬合參數(shù)Tab le 2 The fitting param eters o f HD 62509 b

圖6 HD 62509的Lick視向速度觀測數(shù)據(jù)的Lom b-Scargle周期譜,其中存在一個590 d左右的周期信號.虛線表示1%的FAP,點虛線表示0.1%的FAP,點線表示0.01%的FAP.Fig.6 The Lom b-Sarg le p eriodogram of L ick rad ia l velocity(RV)data of HD 62509.A 590-day p eriod exists in the periodogram.The 1%,0.1%,and 0.01%of FAP are indicated by the dashed line,dot-dashed line,and dotted line,respectively.

圖7 HD 62509的模擬天體測量數(shù)據(jù)的Lom b-Scargle周期譜,其中存在一個590 d左右的周期信號.虛線表示1%的FAP,點虛線表示0.1%的FAP,點線表示0.01%的FAP.Fig.7 The Lom b-Sarg le p eriodogram o f m ock astrom etric data of HD 62509.A 590-day period exists in the p eriodogram.The 1%,0.1%,and 0.01%o f FAP are ind icated by the dashed line,dot-dashed line,and dotted line,respectively.

圖8 HD 62509的歸一化視向速度擬合曲線和擬合行星b后的殘差Fig.8 The norm a lized RV cu rve of HD 62509 and the residua l a fter fitting p lanet b

圖9 左：HD 62509的歸一化天體測量法擬合曲線和擬合行星b后的殘差;右：HD 62509的天體測量法擬合曲線Fig.9 Left:The norm a lized astrom etric cu rve of HD 62509 and the residua l a fter fitting p lanet b;R igh t: The astrom etric cu rve of HD 62509

4.3 G J 876

GJ 876是距離地球4.689 pc[24]的M 4型恒星[25],通過長達20 yr的視向速度觀測發(fā)現(xiàn), GJ 876存在著4顆行星圍繞其公轉(zhuǎn)[26?30].GJ 876行星系統(tǒng)是距離地球最近的多行星系統(tǒng),是首個發(fā)現(xiàn)存在軌道共振(Mean Motion Resonance,MMR)的行星系統(tǒng)[28],同時也是首個發(fā)現(xiàn)存在拉普拉斯共振的系外行星系統(tǒng)[30?32].

在GJ 876系統(tǒng)中,除去以約2 d周期圍繞GJ 876公轉(zhuǎn)的行星d外,行星c(周期約30 d)、行星b(周期約61 d)和行星e(周期約124 d)存在著4:2:1的拉普拉斯共振.由于GJ 876系統(tǒng)中行星與行星間的相互攝動非常強烈,傳統(tǒng)的開普勒橢圓擬合無法很好擬合長達20 yr的視向速度觀測數(shù)據(jù).Laughlin等[33]和Rivera等[34]在2001年先后分別提出對GJ 876系統(tǒng)的N體積分牛頓動力學(xué)數(shù)值解,在傳統(tǒng)的LM算法上加入N體積分算法,求出行星在特定時刻的軌道解,從而降低視向速度數(shù)據(jù)擬合的χ2.

而本文采取的是改進的開普勒橢圓擬合,采用文獻[30]中表3的N體積分牛頓動力學(xué)數(shù)值求解得到的軌道根數(shù),恒星質(zhì)量為0.32 M⊙.為了與文獻[30]中的結(jié)果比較,我們選取在JD 2450602.093至JD 2452428.093內(nèi)生成52組模擬的天體測量數(shù)據(jù),再對其進行改進的開普勒橢圓擬合.

對于使用開普勒橢圓擬合天體測量數(shù)據(jù),需要擬合7個參數(shù)(α,e,i,?,ω,M0,mp),而對于存在2:1軌道共振的行星系統(tǒng),近星點角距存在著快速的線性變化:˙?= |j4|Cre|j4|?2cos?[35],其中?為共振角,j4為共振角中的系數(shù),Cr為含拉普拉斯系數(shù)的參數(shù).對于處于平運動共振的行星,其共振角總是在一定范圍內(nèi)擺動,等式右邊可以視為恒值,即位于2:1軌道共振的行星的近星點存在進動,因此可以額外增加變量˙ω共8個參數(shù)作為單行星的改進開普勒橢圓擬合,改進的開普勒橢圓擬合本質(zhì)上是一個不斷進動的開普勒橢圓擬合.根據(jù)文獻[30]的表3,行星b和行星c的天體測量信號分別為298μas和58 μas,而行星d和行星e的天體測量信號分別只有0.28μas和9.64μas,遠低于Gaia單次天體測量精度34.2μas[7],因此這里只對行星b和c進行共平面的改進開普勒橢圓擬合.

通過使用MCMC算法,最終得到表3和表4中行星b與行星c的擬合結(jié)果,可以看到使用改進的開普勒橢圓擬合有較好的擬合結(jié)果.由于我們并沒有加入GJ 876的視向速度數(shù)據(jù)進行擬合,所以升交點經(jīng)度?與近星點角距ω由近星點經(jīng)度?表示.圖10為行星b和c的天體測量法時序擬合曲線及擬合殘差,圖11為行星b和c的天體測量法曲線及擬合殘差.

表3 GJ 876 b的擬合參數(shù)Tab le 3 The fitting p aram eters o f G J 876 b

表4 GJ 876 c的擬合參數(shù)Tab le 4 The fitting param eters o f G J 876 c

圖10 改進的開普勒擬合行星GJ 876 b和G J 876 c天體測量法時序曲線及擬合殘差Fig.10 The m od ified K ep lerian astrom etric cu rve of GJ 876 b and G J 876 c and the fitting residua l

圖11 行星G J 876 b和G J 876 c改進的開普勒擬合曲線及擬合殘差Fig.11 The m od ified K ep lerian fitting cu rve o f G J 876 b and G J 876 c and the fitting residual

5 討論

5.1 觀測采樣方案的影響

Ford曾針對SIM(Space Interferometer M ission)的觀測方案討論到,對于定點的天體測量觀測,不同的觀測采樣方案會對行星探測效率有不同的影響[4].同時由于天體測量法探測系外行星首先要擬合出恒星的周年視差,而周年視差的周期為1 yr,天體測量法會降低軌道周期在1 yr附近行星的探測效率,因為圍繞類太陽恒星的宜居帶行星的軌道周期在1 yr左右,這對探測宜居帶行星帶來一定的困難[4].

5.2 天體物理噪聲

由于恒星表面的耀斑、黑子、米粒組織和非徑向?qū)α鞯?會造成恒星星像的抖動,造成圖像光心的偏移.這類影響是恒星自身所固有的,因此稱為天體物理噪聲.天體物理噪聲對視向速度和天體測量法均有影響,對于平靜穩(wěn)定的低自轉(zhuǎn)恒星,天體物理噪聲對其視向速度造成3–5m/s的抖動,而對天體測量觀測所造成的抖動的影響可由下式進行估計[36?41]:

σpos、σvR和σm分別是天體測量法、視向速度和星等均方根值.R?為恒星的半徑, v sin i為恒星自轉(zhuǎn)速度徑向投影,Prot為恒星的自轉(zhuǎn)周期.

5.3 視差和自行的加速項

恒星在銀河系中圍繞銀心不斷運動,太陽也在以約2.5億年的周期圍繞銀心運動,視差和自行是由于恒星相對太陽的位置和運動引起的.由于銀河系恒星運動的復(fù)雜性,視差和自行并不是不變的.恒星的視差與自行有如下變化[42]:

μ為恒星的總自行,vr為恒星的視向速度,A為天文單位.對于大多數(shù)恒星,短期內(nèi)的天體測量觀測中自行變化可以忽略,但對一些高自行的恒星,如Barnard星,必須考慮自行的加速項.

5.4 行星間引力作用

對于一些行星與行星間引力作用強烈的系統(tǒng),行星間攝動在觀測時段內(nèi)造成的天體測量法變化較大,是不適合用(3)式進行軌道擬合的,例如GJ 876系統(tǒng)[43?46]、55 Cnc系統(tǒng)[47?49]等.對于這類行星間引力作用強烈的行星系統(tǒng),傳統(tǒng)的開普勒軌道解不能很好地擬合觀測數(shù)據(jù).對于存在軌道共振的系統(tǒng),對于算例所示的GJ 876系統(tǒng),改進的開普勒橢圓擬合有一定實用意義,但此方法不能普適于非共振的行星系統(tǒng),因此求解N體積分的牛頓動力學(xué)數(shù)值軌道解有其必要性.

6 總結(jié)與展望

在本文中,我們針對αCentauri A、HD 62509和GJ 876系統(tǒng),基于Gaia單次天體測量法精度生成多組對應(yīng)的仿真數(shù)據(jù),使用Lomb-Scargle周期譜分析行星軌道周期信號,進而采用MCMC算法對3個行星系統(tǒng)進行軌道反演,解算結(jié)果與行星的初始參數(shù)相吻合.

目前發(fā)現(xiàn)的大多數(shù)系外行星所造成的天體測量信號小于1 mas,而地面的天體測量法精度難以達到1 m as以下,但隨著高精度的空間天體測量衛(wèi)星的升空,未來將有大量通過天體測量法被發(fā)現(xiàn)的系外行星.哈勃空間望遠鏡(HST)上的FGS能以1 m as的天體測量法精度進行窄角天體測量法[50],詹姆斯韋伯空間望遠鏡(JWST)作為哈勃望遠鏡的繼任者,也會安裝上FGS.Gaia在釋放第1批數(shù)據(jù)后,預(yù)期可發(fā)現(xiàn)大量長周期的類木行星.STEP預(yù)期有能力以1μas的窄角天體測量法精度發(fā)現(xiàn)距離太陽15 pc內(nèi)的類地行星.在2016年Gaia公布其觀測數(shù)據(jù)后,我們可以利用本文的行星3維軌道擬合方法來對Gaia的數(shù)據(jù)進行處理,以發(fā)現(xiàn)更多的系外行星并開展科學(xué)研究.另外在將來的工作中,我們還將對行星軌道的擬合方法做進一步改進,在行星開普勒解的基礎(chǔ)上發(fā)展N體積分動力學(xué)數(shù)值軌道解,這將有助于我們更好地了解系外行星系統(tǒng).

[1]M ayor M,Q ueloz D.Natu re,1995,378:355

[2]M u terspaugh M W,Lane B F,K u lkarni S R,et a l.A J,2010,140:1657

[3]Pou rbaix D,Jorissen A.A&AS,2000,145:161

[4]Ford E B.PASP,2004,116:1083

[5]M akarov V V,Beichm an C A,Catanzarite J H,et a l.A p JL,2009,707:L 73

[6]Perrym an M A C,Lindegren L,K ovalevsky J,et al.A&A,1997,323:L 49

[7]Perrym an M,Hartm an J,Lindegren L.A p J,2014,797:14

[8]Perrym an M.The Exop lanet Handbook.New York:Cam b ridge University P ress,2011

[9]K ovalevsky J,Seidelm ann P K.Fundam en tals o f A strom etry.New York:Cam b ridge Un iversity P ress, 2004

[10]Cumm ing A,M arcy G W,Bu tler R P.A p J,1999,526:890

[11]Lom b N R.A&SS,1976,39:447

[12]Scarg le J D.A p J,1982,263:835

[13]P ressW H,Teukolsky S A,Vetterling W T,et al.Num erical Recipes.New York:Cam b ridge University P ress,1992

[14]Ford E B.A p J,2006,642:505

[15]G regory P.Bayesian Logical Data Analysis for the Physical Sciences:A Com parative App roach w ith M athem atica?Supp ort.New York:Cam b ridge Un iversity P ress,2005

[16]G regory P C.M NRAS,2007,381:1607

[17]G regory P C.M NRAS,2007,374:1321

[18]Lattanzi M G,Spagna A,Sozzetti A,et a l.M NRAS,2000,317:211

[19]Dum usque X,Pep e F,Lov is C,et a l.Natu re,2012,491:207

[20]Ra jpau l V,A igrain S,Roberts S.M NRAS,2016,456:L 6

[21]Hoffleit D,Jaschek C.The B right Star Catalogue.5th ed.New Haven:Yale University Observatory, 1991

[22]K ervella P,Th′even in F,S′egransan D,et a l.A&A,2003,404:1087

[23]Reffert S,Quirrenbach A,M itchell D S,et al.Ap J,2006,652:661

[24]Van Leeuw en F.A&A,2007,474:653

[25]von B raun K,Boya jian T S,van Belle G T,et a l.M NRAS,2014,438:2413

[26]M arcy G W,Butler R P,Vogt S S,et al.Ap JL,1998,505:147

[27]Delfosse X,Forveille T,Beuzit J L,et a l.arX iv:astro-ph/9812008,1998

[28]M arcy G W,Butler R P,Fischer D,et al.A p J,2001,556:296

[29]R ivera E J,Lissauer J J,Bu tler R P,et a l.A p J,2005,634:625

[30]Rivera E J,Laugh lin G,Bu tler R P,et al.Ap J,2010,719:890

[31]W ang S,Ji J,Zhou J L.A p J,2012,753:170

[32]W ang S,Ji J.A p J,2014,795:85

[33]Laugh lin G,Cham bers J E.A p JL,2001,551:109

[34]R ivera E J,Lissauer J J.A p J,2001,558:392

[35]M u rray C D,Derm ott S F.So lar System D ynam ics.New York:Cam b ridge Un iversity P ress,1999

[36]Saar S H,Donahue R A.A p J,1997,485:319

[37]Queloz D.11th Cam b ridge W orkshop on Coo l Stars,Stellar System s and the Sun,2001,223:59

[38]Hen ry G W,Donahue R A,Ba liunas S L.A p JL,2002,577:111

[39]Saar S,Seager S.arX iv:astro-ph/0305436,2003

[40]Hatzes A P.AN,2002,323:392

[41]Pau lson D B,Coch ran W D,Hatzes A P.AJ,2004,127:3579

[42]D rav ins D,G u llberg D,Lindegren L,et a l.ASPC,1999,185:41

[43]Ji J,Li G,Liu L.A p J,2002,572:1041

[44]Ji J,Liu L,K inosh ita H,et a l.A p JL,2003,591:57

[45]季江徽,劉林.天文學(xué)報,2006,47:402

[46]Ji J H,Liu L.ChA&A,2007,31:187

[47]Ji J H,K inoshita H,Liu L,et al.A p JL,2003,585:139

[48]Ji J H,K inoshita H,Liu L,et al.RAA,2009,9:703

[49]Ji J H,Lin L,Zhou J L,et a l.ChA&A,2003,27:127

[50]Nelan E P,Lup ie O L,M cA rthu r B,et a l.SP IE,1998:237

A nalogue Sim u lation and O rb ital Solving A lgorithm of A strom etric Exop lanet Detection

HUANG Ping-hui1,2,3JIJiang-hui1,2

(1 Pu rp le M oun tain O bserva to ry,Chinese A cadem y o f Scien ces,Nan jing 210008) (2 K ey Labo ra to ry o f P laneta ry Scien ces,Chinese A cadem y o f Scien ces,Nan jing 210008) (3 Un iversity o f Chinese A cadem y o f Scien ces,Beijing 100049)

Astrometry is an effectivemethod to detect exoplanets.It hasmany advantages that other detectionmethods do not bear,such as providing three dimensional planetary orbit and determ ining the planetarymass.Astrometry w illenrich the sam ple of exop lanets.As the high-precision astrometric satellite Gaia(G lobal Astrom etry interferometer for Astrophysics)was launched in 2013,therew illbe abundant long-period Jupiter-size p lanets to be discovered by Gaia.In this paper,we specify theαCentauri A,HD 62509,and GJ 876 systems,and generate the synthetic astrometric data w ith the single astrom etric precision of Gaia.Then we use the Lomb-Scargle periodogram to analyse the signature of planets and the Markov Chain Monte Carlo(MCMC)algorithm to fit the orbit of p lanets.The simulation results are well coincide w ith the initial solutions.

astrometry,planets and satellites:evolution,astrometric satellites:Hipparcos,Gaia,m ethods:Lomb-Scargle periodogram,M arkov Chain M onte Carlo(MCMC)

P144;

A

10.15940/j.cnki.0001-5245.2016.05.007

2016-01-29收到原稿,2016-02-27收到修改稿

?國家自然科學(xué)基金項目(11273068,11473073)、中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(XDB09000000)、中國科學(xué)院空間科學(xué)戰(zhàn)略性先導(dǎo)科技專項(XDA 04060900)、中國科學(xué)院新興與交叉學(xué)科布局試點項目(KJZD-EW-Z001)及紫金山天文臺小行星基金會資助

?phhuang@pmo.ac.cn

?jijh@pm o.ac.cn