呂振偉，李木子，郭越凡，張 帆

(北京師范大學天文系，北京100875)

1 引 言

物體的運動導致空間中質(zhì)量能量分布的改變，從而引起時空曲率的微小變化；而時空的曲率又決定著物體如何運動。時空曲率的擾動如漣漪一般在宇宙中以光速傳播且攜帶能量，這就是基于1916年愛因斯坦廣義相對論所預言的引力波輻射。2016年2月，LIGO(Laser Interferometer Gravitational Wave Observatory)科學合作組織宣布位于美國Hanford和Livingston 的兩臺激光干涉引力波探測器首次成功捕捉到來自雙黑洞并合的引力波信號GW150914[1]。這是引力波天文學史上具有里程碑意義的事件，為百年來物理學家們對于引力波是否為可測量的物理波的爭論畫上了句號，且再次驗證了廣義相對論的正確性。

目前高新LIGO[2]和高新VIRGO(Virgo interferometer)[3]已經(jīng)升級完畢并投入到聯(lián)合探測中，截止到2017年底(O1 和O2 已運行，2019年將開始O3 運行)，已有5 個雙黑洞系統(tǒng)并合引力波信號[1,4–7]及1 個雙中子星并合信號GW1710817[8]被HLV (Hanford-Livingston-Virgo)引力波探測器網(wǎng)絡(luò)成功探測到。對雙黑洞系統(tǒng)的質(zhì)量估計在表1 中列出。為降低低頻上的地震噪聲，日本3 km 臂長的KAGRA 引力波探測器建于地下，且應(yīng)用了冷凍鏡面的技術(shù)，以改善熱噪聲的干擾，于2016年3月開機并進行了為期1 個月的試運行，預計2018年將正式投入觀測[9]。2024年，LIGO-india 也會加入到全球地基引力波探測器網(wǎng)絡(luò)中[10]。空間引力波探測項目如激光干涉空間天線(Laser Interferometer Space Antenna, LISA[11])，DECIGO[12]，太極(TAIJI[13,14])和天琴(TIANQIN[15])的優(yōu)勢是可在無需考慮地球曲率的情況下增加臂長，且避免了地球上的地震噪聲和人為噪聲的干擾，使其在低頻上的靈敏度大大提升。LISA 的設(shè)計臂長為2.5×106km，工作頻段在10?4～10?1Hz，其科學目標為探測星系中心超大質(zhì)量雙黑洞系統(tǒng)并合、極端質(zhì)量比雙黑洞并合等低頻引力波信號。

表1 已觀測到的雙黑洞并合引力波事件(M⊙表示太陽質(zhì)量)

作為地基引力波探測器和空間引力波探測器的重點目標，雙黑洞系統(tǒng)附近的物質(zhì)通常比較稀薄，因此其繞轉(zhuǎn)-并合- 鈴宕的演化過程在傳統(tǒng)天文觀測的電磁波段上不可見。而在該過程中，系統(tǒng)能量以引力波的形式向外輻射，且由于引力波與物質(zhì)之間的相互作用比較微弱，引力波信號可從被極強引力效應(yīng)支配的近密雙星并合發(fā)生區(qū)域逃逸，并且攜帶黑洞系統(tǒng)的直接信息，如雙黑洞的質(zhì)量、自旋、軌道離心率等。當高新LIGO 和高新VIRGO 在其設(shè)計靈敏度上運行或者第三代引力波探測器建設(shè)完成時，我們將能探測到雙黑洞系統(tǒng)各階段引力波輻射的更多細節(jié)，而高探測率帶來的大量信號樣本既是進一步研究黑洞物理的寶貴資源，又是對引力波信號處理和波源參數(shù)估計方法的全新挑戰(zhàn)。

目前關(guān)于雙黑洞系統(tǒng)起源的研究還不深入，而我們特別關(guān)注恒星級雙黑洞系統(tǒng)，因該質(zhì)量范圍的黑洞在宇宙中大量存在，且有間接的光學觀測可相互驗證。普遍認為恒星級雙黑洞系統(tǒng)的起源主要有大質(zhì)量雙星演化起源和動力學起源[16–18]這兩種情況。不同起源的系統(tǒng)具有不同的黑洞自旋和質(zhì)量分布，大質(zhì)量雙星演化起源會使得雙黑洞自旋角動量方向和公轉(zhuǎn)軌道角動量方向具有高度同向性，而動力學起源則會導致自旋朝向的隨機分布。前身星核坍縮后的爆發(fā)一般是各向異性的，這可能會對自旋-軌道角動量方向夾角有重要影響[19,20]，因此大質(zhì)量雙星演化是否能如理論預言的一樣，使自旋-軌道角動量同向還需要驗證。我們?nèi)匀患僭O(shè)自旋-軌道角動量是同向的，并且有小角度彌散。如何利用引力波信號攜帶的波源信息進行參數(shù)估計，自旋-軌道角動量方向的同向性是否能用來推斷系統(tǒng)的起源，回答這些問題對進一步理解雙黑洞系統(tǒng)背后的物理至關(guān)重要。

本文主要面向剛?cè)腴T引力波研究領(lǐng)域的學者和學生，介紹如何通過分析引力波信號推斷波源參數(shù)來限制雙黑洞系統(tǒng)起源與演化模型的方法。我們假設(shè)黑洞形成的天體物理過程決定了雙黑洞質(zhì)量和自旋等參數(shù)分布的大致輪廓；而利用引力波信號推斷的波源質(zhì)量自旋信息和動力學方法測量的結(jié)果來驗證不同起源雙黑洞系統(tǒng)質(zhì)量、自旋分布的理論預言。綜合考慮理論模型和實驗觀測的結(jié)果，就可對于不同起源雙黑洞系統(tǒng)給出一個參數(shù)化的模型來描述黑洞質(zhì)量分布與自旋分布。隨后，我們對所構(gòu)建的模型進行采樣并模擬一個數(shù)據(jù)集，其中一部分為具有自旋-軌道同向性的大質(zhì)量雙星演化起源，另一部分為動力學起源。將該數(shù)據(jù)集作為LALInference[21]的輸入(該軟件提供了引力波參數(shù)估計的標準方法，見4.1 節(jié))，我們得到質(zhì)量和自旋等參數(shù)的后驗分布。進而計算貝葉斯系數(shù)(Bayes factors)和總質(zhì)量差等對模擬系統(tǒng)起源進行推斷，得到的結(jié)果可驗證利用自旋分布作為判據(jù)區(qū)分黑洞系統(tǒng)起源的有效性，回答數(shù)據(jù)集中不同起源所占的相對比例，以及不同起源總質(zhì)量分布的差異等問題。

本文的結(jié)構(gòu)如下，第2 章中介紹已有的各種理論演化模型，總結(jié)它們的異同；第3 章中，我們結(jié)合理論模型和觀測結(jié)果，重新構(gòu)建黑洞的參數(shù)模型；在第4 章中，進行采樣并得到模擬數(shù)據(jù)集，從而利用LALInference 軟件包將理論預言不同起源的雙黑洞系統(tǒng)質(zhì)量、自旋等參數(shù)的分布情況轉(zhuǎn)化成可以與引力波觀測結(jié)果直接比對的后驗分布，進而達成限定諸如“不同起源雙黑洞系統(tǒng)所占的比例”這類基本天體物理問題的目的。第5 章中，我們對本方法的結(jié)果進行討論，總結(jié)現(xiàn)階段工作中還需提升的部分及對未來工作的展望。

2 恒星級黑洞的觀測和理論模型

根據(jù)引力波探測的結(jié)果，人們發(fā)現(xiàn)了另外一種恒星級黑洞質(zhì)量分布規(guī)律。如首次探測到的雙黑洞質(zhì)量分別為29M⊙和36M⊙，并合后黑洞質(zhì)量更是高達62M⊙[1]，以前從未探測到類似結(jié)果。而且，多數(shù)探測到的雙黑洞并合事件的黑洞質(zhì)量在20M⊙～40M⊙之間[1,4–7]，這為黑洞以及雙黑洞系統(tǒng)存在多種起源的觀點提供了觀測支持。

恒星級黑洞的理論模型主要有兩種。一種是綜合了恒星模型和超新星流體動力學模擬的結(jié)果，由Belczynski 等人[23]在2012年提出的大質(zhì)量恒星演化的物理模型，該模型解釋了瞬態(tài)低質(zhì)量X 射線雙星中黑洞的質(zhì)量觀測結(jié)果，并且合理地解釋在中子星和黑洞之間的質(zhì)量缺失問題。但對于黑洞質(zhì)量大于15M⊙的情況，無法進行描述。另一種是脈沖對不穩(wěn)定型超新星(pulsational pair instability supernovae)理論[24,25]，該理論可以解釋引力波探測的大質(zhì)量黑洞現(xiàn)象，這主要是由于大量的物質(zhì)在坍縮之前被噴出。

Belczynski 等人[23]提出的大質(zhì)量恒星演化的物理過程，對超新星爆發(fā)內(nèi)部工作機制給出了嚴格的限制。特別是核坍縮型超新星是在恒星坍縮后100～200 ms 爆發(fā)的，這意味著爆發(fā)是由10～20 ms 內(nèi)生長起來的不穩(wěn)定性驅(qū)動的?；蛘哒f，如果將來的觀測發(fā)現(xiàn)了缺失的黑洞，那表明這種不穩(wěn)定性需要大于200 ms 的生長時間，這樣才能形成質(zhì)量為3M⊙～5M⊙的致密天體。

脈沖對不穩(wěn)定性超新星理論[24,25]可以解釋引力波觀測中的大質(zhì)量黑洞現(xiàn)象。由于大量的物質(zhì)在坍縮之前被噴出，這種不穩(wěn)定性發(fā)生在大質(zhì)量恒星(80M⊙～150M⊙)演化的后期，伴隨著高溫低密度區(qū)引起的正負電子對熱聚集，會對狀態(tài)方程有重要的影響，也只有質(zhì)量在80M⊙～150M⊙范圍內(nèi)的恒星有足夠高的熵可以產(chǎn)生這種不穩(wěn)定性。由此產(chǎn)生的黑洞質(zhì)量小于40M⊙，并且大多數(shù)黑洞的質(zhì)量約為34M⊙，這與引力波觀測的結(jié)果一致。對于質(zhì)量在150M⊙～250M⊙之間的原恒星，由于對不穩(wěn)定性(pair instability)而不能形成黑洞。

3 不同起源的雙黑洞系統(tǒng)的質(zhì)量分布建模

目前恒星級雙黑洞的起源主要有兩種：大質(zhì)量雙星演化起源和動力學起源[16–18,26]。大質(zhì)量雙星演化發(fā)生在星系場(galactic fields)中，雙星通過大質(zhì)量雙星演化形成雙黑洞，具有自旋和公轉(zhuǎn)軌道高度同向的特點。動力學起源發(fā)生在球狀星團(globular clusters)等恒星密度高的環(huán)境中，自旋和軌道一般是隨機的。下面我們先構(gòu)造雙黑洞系統(tǒng)中質(zhì)量較大的主黑洞(m1)的質(zhì)量分布模型，然后應(yīng)用冪律分布構(gòu)造質(zhì)量較小的次黑洞(m2)的質(zhì)量分布[27,28]。由于光學觀測上黑洞質(zhì)量在5M⊙～15M⊙范圍內(nèi)，存在2M⊙～5M⊙范圍內(nèi)的質(zhì)量缺失，并且大多數(shù)在(7.8±1.2)M⊙范圍內(nèi)[22]。結(jié)合Belczynski 的理論模擬結(jié)果，可以假設(shè)小質(zhì)量黑洞(5M⊙～20M⊙)服從冪律分布(power-law distribution)，并且M5M⊙。由引力波探測得到的黑洞質(zhì)量則大得多，主要分布在20M⊙～40M⊙之間[1,4–7]；同時，利用脈沖對不穩(wěn)定性超新星理論得出：大質(zhì)量恒星(80M⊙～150M⊙)演化得到的黑洞質(zhì)量小于40M⊙，并且大多數(shù)為34M⊙。因此，可以用高斯分布來描述大質(zhì)量黑洞(20M⊙～40M⊙)的分布，兩部分疊加得到[24]：

其中，λ是兩部分(主黑洞和次黑洞)所占比例，H(mmax?m1)是截止函數(shù)，mmax是冪律分布的截止質(zhì)量，mpp是高斯分布的平均質(zhì)量，方差為σpp，Λ={α,β,λ,mmin,mmax,mpp,σpp}，mmin是冪律分布的初始質(zhì)量，p(m1|Λ)是在Λ的條件下m1的條件分布。例如，可以取ΛA={1.5,3,0.2,3,35,38,1.5}和ΛB={1.5,2,0.1,3,35,35,1}，分別作為大質(zhì)量雙星演化起源(模型A)和動力學起源(模型B)的質(zhì)量分布，用來進行MCMC 采樣產(chǎn)生樣本，如圖1所示。

圖1 模型A 和模型B 的質(zhì)量分布

對于來自大質(zhì)量雙星演化的雙黑洞系統(tǒng)，我們假設(shè)黑洞誕生時爆發(fā)的各向異性對自旋傾角的影響很小，它們的自旋軌道傾角在[0?,10?]內(nèi)均勻分布或者符合高斯分布；而對于來自動力學起源的雙黑洞系統(tǒng)，其自旋和軌道朝向是隨機分布的。最后用馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo, MCMC)或者多級抽樣(nested sampling)[21]的方法進行采樣，得到一組引力波事件。把這些事件對應(yīng)的參數(shù)作為樣本輸入LALInference，可以得到相應(yīng)的后驗概率。

4 貝葉斯推理及結(jié)果

貝葉斯推理的特點是每次探測到的信號都可以用以前的后驗概率作為先驗概率，然后計算新的后驗概率，通過連續(xù)疊加得到最終的后驗概率。這是一種被天文學家廣泛應(yīng)用且行之有效的檢驗方式，可以有效地減少計算量，節(jié)約大量時間和計算成本，這對于計算量非常大的引力波數(shù)據(jù)分析工作至關(guān)重要。在接下來的數(shù)年內(nèi)，我們會得到越來越多的引力波事件，如何精確地推斷致密雙星的起源，不同起源所占的比例，以及質(zhì)量自旋等參數(shù)的分布情況，是需要解決的三個主要問題。

4.1 LALInference模擬引力波模板并進行貝葉斯推理

該軟件是由C 語言編寫，并且包含了用Python 寫的后期處理工具，是處理引力波探測數(shù)據(jù)的標準方法，可以用來計算功率譜和給定噪聲曲線的穩(wěn)定高斯噪聲，還可以使用所有的引力波模板；在給定模板和模型參數(shù)的條件下，得到數(shù)據(jù)的似然函數(shù)，以及參數(shù)空間的隨機采樣等。下面對LALInference 的基本原理進行簡單介紹。

本文主要采用其參數(shù)估計部分，當輸入一個引力波事件的數(shù)據(jù)時，LALIference 會比較輸入的探測信號和內(nèi)置的理論模板，得到關(guān)于該事件參數(shù)(質(zhì)量、自旋等)的似然函數(shù)(likelihood function)，也就可得到在探測到該事件的條件下參數(shù)的條件分布，即參數(shù)的后驗分布：

其中，θ包括了所有的模板參數(shù)，且∫p(θ)dθ= 1，d是探測到的數(shù)據(jù)集合。p(θ)是基于目前所有的信息假設(shè)的先驗分布。p(d|θ)是似然函數(shù)，可以由探測信號的統(tǒng)計性質(zhì)給出。p(d)被稱為統(tǒng)計證據(jù)，是似然函數(shù)和先驗相乘后對所有參數(shù)的求和，在參數(shù)估計時，p(d)被用來歸一化。當我們輸入不同的模板參數(shù)，LALInference 內(nèi)置模板函數(shù)可以計算相應(yīng)的波形模板，然后通過比較波形模板和探測數(shù)據(jù)得到似然函數(shù)。

一般我們假設(shè)探測器的噪聲是穩(wěn)定的高斯噪聲，因此，噪聲的似然函數(shù)可以用簡單的高斯分布來描述：

圖2 雙黑洞并合引力波事件GW151226 的噪聲功率譜

當探測數(shù)據(jù)中有引力波信號時，假設(shè)信號是h，那么n=d ?h，包含引力波信號的似然函數(shù)：

下面我們由LALInference 相關(guān)函數(shù)得到的不同質(zhì)量黑洞系統(tǒng)的引力波信號波形，進一步了解相應(yīng)的頻率和繞轉(zhuǎn)速度變化曲線。我們在計算時，設(shè)雙黑洞質(zhì)量相等，即m1=m2。圖3 和圖4 僅給出引力波的h+偏振。

圖3 不同質(zhì)量的雙黑洞系統(tǒng)引力波波形

圖4 不同質(zhì)量的雙黑洞系統(tǒng)引力波波形

雙黑洞質(zhì)量相等(m1=m2)時，在相互繞轉(zhuǎn)到并合時期的頻率和繞轉(zhuǎn)速度隨時間變化情況如圖5 和圖6 所示?？梢钥吹?，在并合前數(shù)毫秒，頻率和繞轉(zhuǎn)速度指數(shù)增加。

圖5 并合前引力波頻率隨時間變化

圖6 并合前雙黑洞繞轉(zhuǎn)速度隨時間變化

對于先驗概率分布，假設(shè)起源A 的質(zhì)量在10M⊙～80M⊙之間，為均勻分布；起源B的質(zhì)量在10M⊙～80M⊙之間，呈高斯分布，方差σ2= 400。由此可以進一步計算后驗分布，結(jié)果如圖7 所示。

圖7 質(zhì)量分布采樣

對于引力波源GW150914，假設(shè)質(zhì)量自旋等參數(shù)都是均勻的先驗分布，通過5 000 個樣本計算后驗概率?？梢缘玫劫|(zhì)量估計如果用更多的樣本，就能得到更加嚴格的限制。某些參數(shù)組合的散點分布如圖8 所示。

圖8 質(zhì)量m1 和m2 后驗分布以及赤經(jīng)和赤緯后驗分布的散點圖

4.2 推斷引力波波源

假設(shè)探測到一個新的事件，通過參數(shù)估計可以得到波源各個參數(shù)的分布(貝葉斯估計得到的都是概率描述，概率論的觀點得到的參數(shù)是確定值)。該并合事件的起源及其概率可以通過計算貝葉斯系數(shù)來進行推斷[24]。假設(shè)HA是大質(zhì)量雙星演化的模型假設(shè)，HB是動力學起源的模型假設(shè)，d是模擬的引力波波形數(shù)據(jù)，則在模型HA的條件下，得到事件d的概率：

由于樣本是服從p(θ|d)分布的，所以上式中的積分可以轉(zhuǎn)化為求和：

式中，θ是所有的模型參數(shù)，對這些參數(shù)積分后得到的統(tǒng)計證據(jù)p(d|HA)是與具體參數(shù)無關(guān)的量。p(θi)是假設(shè)的先驗分布，p(θi|H)的后驗分布皆可通過LALInference[21]軟件得到。模型A 和模型B 下，我們進一步比較引力波源的統(tǒng)計證據(jù)，得到：

其中，表示模型假設(shè)A 和B 的統(tǒng)計證據(jù)之比。當1 時，表示引力波信號來自模型A 的可能性比來自B 的可能性大很多。一般選取= 8 作為模型判斷的標準，當8 時，我們就認為信號來自模型A。

4.3 推斷不同起源所占比例

假如我們有大量探測事件，如何得到每種起源所占的相對比例？以及大概需要多少事件？我們可以模擬產(chǎn)生引力波信號樣本，其中給定比例的事件來自大質(zhì)量雙星演化，其自旋和軌道朝向高度一致；另一部分來自動力學起源。利用貝葉斯推理，來自大質(zhì)量雙星演化的樣本所占的比例，可以被有效地區(qū)分出來。也就是說，可以用來檢驗引力波是否可以用來區(qū)分自旋和軌道同向與否。這對于理解致密雙星演化動力學過程至關(guān)重要。

假設(shè)起源A 所占的比例為CI，B 是CII，CI+CII=1，則后驗概率為：

其中，

由于在某個起源下獲得引力波信號的概率與該起源所占的比例無關(guān)，所以可得p(dk|HICI)=p(dk|HI),p(dk|HIICII)=p(dk|HII)，并且p(HI|CI)=CI,p(HII|CII)=CII=1?CI。則式(10)可寫為：

為便于計算，等式兩邊同時取對數(shù)，得到：

其中，p(CI|d)是在數(shù)據(jù)d的條件下，得到起源A 在所有事件中所占的比例CI的分布。p(CI)是給定的先驗分布，p(dk|CI)是在CI條件下探測到事件k的概率。隨著事件數(shù)的增加，可以得到p(CI|d)隨事件個數(shù)變化的規(guī)律，即可估計出確定每種起源所占的相對比例的事件數(shù)，這對于雙黑洞系統(tǒng)的起源研究和雙黑洞并合率的研究至關(guān)重要。

4.4 質(zhì)量分布的差異

我們假設(shè)不同起源的雙黑洞具有不同的質(zhì)量分布，同時也具有不同的自旋分布，然后使用自旋分布信息來區(qū)分不同起源的黑洞，檢驗不同起源質(zhì)量分布的差異，例如，對兩種起源的總質(zhì)量進行比較。換言之，我們采用自旋表征的樣本來對不同起源的質(zhì)量分布進行估算。由金斯半徑估算可得，場星系里由大質(zhì)量雙星演化誕生的雙黑洞質(zhì)量比球狀星團里的要小一些，因此可假設(shè)大質(zhì)量雙星演化的黑洞比動力學起源的黑洞總質(zhì)量小。通過對分布π(m|Λ,H)∝p(m1|Λ,H)·p(m2|m1,Λ,H)·Vobs(m)進行采樣，從而得到不同起源的質(zhì)量分布[24]。q=MI?MII的分布，其中MI=mI1+mI2,MII=mII1+mII2，p(q)由下式表示：

其中，

上式中，p(MII|di)=p(MII|II,di)p(II|di)+p(MI|I,di)p(I|di)。

以上采用的是貝葉斯等級推理(Bayesian hierarchical inference)的方法，數(shù)據(jù)模擬可以通過LALInference 軟件來完成。

5 總結(jié)與展望

本文介紹了引力波信號的后期處理，主要是用LALInference 模擬引力波模型并進行貝葉斯推理；然后利用貝葉斯公式推斷引力波起源、不同起源的比例以及不同起源質(zhì)量分布的差異檢驗公式，對貝葉斯等級推理的公式進行系統(tǒng)的推導；最后從理論上討論了貝葉斯方法的可行性。

由致密雙星并合產(chǎn)生的引力波信號包含了豐富的波源信息，包括致密雙星的質(zhì)量、自旋、軌道傾角、離心率等。我們主要研究質(zhì)量和自旋兩個參數(shù)在推斷引力波源和內(nèi)稟參數(shù)分布上的應(yīng)用。事實上，引力波的波形中包含的參數(shù)高達15 個，當我們考慮更多的參數(shù)時，很難找到多參數(shù)分布的簡單描述，計算量也會大大增加，實際操作受到很大的影響。這時可以考慮采用機器學習的方法，通過對構(gòu)建的多層神經(jīng)網(wǎng)絡(luò)進行訓練，可能會給出更加精確的解。并且，深度神經(jīng)網(wǎng)絡(luò)的方法雖然在訓練深度網(wǎng)絡(luò)模型時會耗費大量的時間，但是對于已經(jīng)訓練好的神經(jīng)網(wǎng)絡(luò)模型，幾乎可以在瞬間對上述問題進行判斷。以上的問題是未來引力波天文學的主要問題之一，我們希望隨著觀測事件的增加，能提出更好的理論模型來解釋這些問題。