沒有它們，給你引力波你也不認識！ | 賽先生天文
● ● ●
封面圖：目前LIGO-Virgo官方證認的50個引力波信號，包含了GWTC-1和GWTC-2，圖片來源：LIGO-Virgo | Frank Elavsky, Aaron Geller | Northwestern
撰文 | 吳仕超、趙天宇（北京師范大學）
責編 | 韓越揚、呂浩然
通常來講，飛機和雷達之間是“貓和老鼠”般的競爭關系：雷達發(fā)射電磁波，并利用回波尋找飛機，而回波會因為信號衰減而變得微弱。于是，雷達探測飛機則會變成一個典型的弱信號、強噪聲的問題。人們由此發(fā)展出了相應的匹配濾波（matched filtering）數據處理方法，以提高信號的信噪比[1]，使之更容易被探測到。雷達接收到的回波和其發(fā)射的電磁波具有基本相同的波形，所以雷達的匹配濾波可以直接實現(xiàn)。
而引力波的探測也存在相似的地方：從數據分析的角度來講，引力波探測也是典型的弱信號、強噪聲問題。相應地，人們自然也想到了通過使用匹配濾波的方法來增強引力波的探測能力。
為了實現(xiàn)匹配濾波，我們需要預先知道信號的形狀，也就是波形（waveform）。由于引力波探測器是被動接收引力波信號，因此它的波形不能預先得知。對此，人們的應對辦法是從理論上，針對各種引力波波源構造引力波的波形庫（template bank），也就是所謂的引力波模板（template）。
因為實際的引力波信號，波形受波源參數影響且參數未知，所以引力波模板庫通常需要覆蓋整個可能的波源參數空間，一個模板庫通常包含數百萬個引力波模板。不準確的引力波模板會與實際引力波信號波形不吻合，損失匹配濾波信噪比，從而導致引力波信號丟失。所以引力波波源建模至關重要。
圖1：GW150914在引力波探測器中的波形及軌道演化示意圖，圖片來源：文獻[2]。
在引力波被實際探測到之前，雙星并合系統(tǒng)是科學家預測的最有可能的波源，同時也是引力波天文學研究中最重要的引力波波源之一。因此，它也是人們首先想到要構建引力波模板的系統(tǒng)。
雙星系統(tǒng)中的雙黑洞系統(tǒng)不涉及復雜的潮汐相互作用，是理論描述上最“干凈”的引力波波源系統(tǒng)。為了獲得雙黑洞系統(tǒng)的引力波波形，人們需要也只需要求解愛因斯坦場方程即可。但愛因斯坦場方程的復雜性是出了名的！一兩千頁的專著《愛因斯坦場方程的精確解》不干別的，只為了列舉到現(xiàn)在為止已知的愛因斯坦場方程的解。
圖2：非常“簡單”的愛因斯坦場方程
更“可惜”的是，一兩千頁中與天文系統(tǒng)有對應關系的解并不多。換句話說，愛因斯坦場方程的解析求解非常困難。愛因斯坦本人也深諳其引力場方程的困難程度，他本人以及后來非常多的科學家們發(fā)展出弱場、低速、微擾等近似愛因斯坦場方程的計算分析辦法。但對于引力波輻射最強的雙星并合過程，雙黑洞的相對速度在臨近并合前會越來越快，最后接近光速（如圖1），引力場強超過地球附近場強6個量級，廣義相對論效應強到一切近似計算方法都失效。
于是，人們想到了數值方法求解愛因斯坦場方程，也就是通常人們說到的數值相對論（numerical relativity）?；谶@個設想，基普·索恩（Kip Thorne）早年畫出了雙黑洞并合引力波波形的示意圖（如圖3）。
圖3：索恩設想的雙黑洞引力波波形（下方坐標軸內），圖片來源：Kip Thorne。
既然都“淪落”到數值求解愛因斯坦場方程了，無非就是寫寫程序，花費些電力資源，喝喝咖啡等著波形出來就可以了。這也是人們最早對數值相對論的預期。帶著這個想法，Susan G. Hahn和Richard W. Lindquist最早在1964年就開始嘗試數值求解雙黑洞的愛因斯坦場方程。
Hahn是計算數學大家Peter Lax的學生，當時他在IBM上班。Lindquist是有名的廣義相對論學家。等他們寫好程序，啟動計算機，坐下來準備喝咖啡的時候，預期以外的事情發(fā)生了：程序運行十多步就崩潰了, 輸出了一大堆NaN（編者注：NaN在計算機科學中代表一堆意義不明的值），沒有得到任何有物理意義的結果。兩人瞪直眼睛，為什么？這件事開啟了數值相對論學家們長達近半個世紀的苦思。這個問題被人們稱之為數值相對論的穩(wěn)定性問題。
廣義協(xié)變性原理是廣義相對論的理論基礎之一。該原理從著名的馬赫原理發(fā)展而來，愛因斯坦也以之為傲。該原理使得愛因斯坦場方程以張量形式出現(xiàn)，與數值計算所需的偏微分方程僅一紙之隔，人們需要選擇坐標系（圖4）把愛因斯坦場方程寫成偏微分方程的形式。
20世紀70年代，普林斯頓大學的Kenneth Eppley和德克薩斯大學奧斯汀分校的Larry Smarr對于愛因斯坦場方程數值計算的坐標選擇做了深入研究，發(fā)現(xiàn)坐標選擇會嚴重影響數值計算的穩(wěn)定性。他們針對一些特殊物理情形, 充分發(fā)掘其內在對稱性, 發(fā)展了特殊、精致的坐標以供數值計算使用, 取得了不錯的效果。但缺點是這些坐標選擇不具有通用性, 無法對非對稱性時空做一般性推廣。他們還提出了自動搜索對稱性的坐標條件等, 但發(fā)現(xiàn)還是不能保證計算的穩(wěn)定性。
20世紀80年代，Tsvi Piran對一些簡單的模型化問題給出了初步的數值計算結果。這種利用時空對稱性選擇坐標系的方式在引力/規(guī)范對偶數值相對論問題中也被經常使用。但可惜的是，我們面對的真實引力波波源系統(tǒng)往往什么對稱性都沒有，比如雙黑洞并合系統(tǒng)。
圖4：Kruskal-Szekeres坐標下的史瓦西黑洞，圖片來源：文獻[3]。
直到20世紀90年代初，LIGO的初步硬件（initial LIGO）建設正式啟動[4]。眼看引力波探測的數據處理很快就要開始，人們對引力波模板的需求迫在眉睫。人們急紅了眼，猜測數值相對論所遇到的不穩(wěn)定性問題會不會是由于計算所用網格不夠精細導致的？為了把計算網格鋪得足夠細，就需要用非常多的計算節(jié)點并行計算。為了把復雜的愛因斯坦場方程和復雜的并行計算綜合到一起，就需要許多人的協(xié)作。
綜上，美國多個研究所和高校聯(lián)合組成“雙黑洞大挑戰(zhàn)聯(lián)盟”（The Binary Black Hole Grand Challenge Alliance，圖5）, 計劃構建大型數值相對論軟件以實現(xiàn)對雙黑洞問題的數值模擬。一款名為Cactus的軟件[5]就是在這個背景下誕生的，而如今數值相對論領域最成熟的開源軟件Einstein Toolkit[6]就是基于Cactus發(fā)展起來的。但該聯(lián)合團隊的研究結果卻表明，數值相對論的不穩(wěn)定性問題并非只是計算分辨率的問題。雙黑洞大挑戰(zhàn)聯(lián)盟后來也就不了了之了。
圖5：上世紀90年代雙黑洞大挑戰(zhàn)聯(lián)盟的介紹網頁（http://www.crpc.rice.edu/CRPC/demos/BlackHole/）
數值計算不穩(wěn)定是不是因為所計算的偏微分方程本身就不穩(wěn)定呢？針對這個想法，若干的數值相對論工作者從偏微分方程的角度出發(fā)分析計算所用偏微分方程的雙曲性質。人們發(fā)現(xiàn)，雖然各種各樣的用來數值計算的偏微分方程形式都來自于原來的愛因斯坦場方程，但各自的偏微分方程性質卻千差萬別。
基于對稱雙曲性質的要求，加州理工學院和康奈爾大學聯(lián)合數值相對論小組提出了廣義調和坐標形式的計算用愛因斯坦偏微分方程。這些方程形式在2003年之前就被提出來了，但數值相對論的不穩(wěn)定問題依然進展緩慢。
Bernd Brügmann教授1996年博士畢業(yè)，論文是圈量子引力。圈量子引力界的著名物理學家Martin Bojowald和Thomas Thiemann都是他的師弟?？釔壑袊鴩宓腂rügmann教授喜歡通過簡單的組合得到復雜事物，這促成他從圈量子引力到數值相對論的研究興趣轉變。
不同于雙黑洞大挑戰(zhàn)聯(lián)盟關于計算網格問題的想法，Brügmann教授把雙黑洞問題當作多尺度物理問題來看待，他最早把科學計算領域的自適應網格細化（adaptive mesh refinement）技術引入到數值相對論（圖6），建立了現(xiàn)在稱之為BAM的數值相對論軟件。Brügmann教授喜歡合作，但只限于十人以下的合作。這就是為什么他在德國馬普引力所工作期間拒絕了把BAM引入到Cactus當中。所以Cactus的自適應網格細化大科學計算平臺部分是2003年之后才由Eric Schnetter發(fā)展起來的，現(xiàn)在被稱為Carpet。
1999年，Brügmann教授基于他發(fā)展的自適應網格細化技術把雙黑洞的數值計算時間推進到了35個時間單位（數值計算把雙黑洞問題無量綱化后的時間單位）。以現(xiàn)在計算雙黑洞約1000個時間單位的標準來看，35個時間單位微不足道，但在當時這已經是不小的進展了。到2000年，Brügmann教授和合作者一起把計算時間推進到了50個時間單位。
圖6：雙黑洞計算中的自適應網格細化，越靠近黑洞的區(qū)域網格越密集。圖片來源：文獻[7]
2000年對引力波發(fā)展是個特殊的年份，一方面LIGO初步硬件基本建成，眼看就可以開展數據采集了；另一方面，Kip Thorne這一年60歲大壽。對引力波探測的極度樂觀和對數值相對論發(fā)展的極度悲觀使得索恩在其生日宴會上說“很可能引力波的成功探測比數值相對論的突破還要早日實現(xiàn)”。當然，后來的發(fā)展和索恩的這句話正好相反。
2000年以后，數值相對論繼續(xù)發(fā)展。Brügmann教授與其合作者相繼提出了黑洞穿刺的處理方式和針對BSSN計算方程形式的特殊坐標選擇條件。這些技術把雙黑洞計算時間推進到150個時間單位。也許做事精細是德國人的特點，Brügmann教授一直“固執(zhí)”地認為黑洞穿刺點存在坐標奇點，處理坐標奇點的特殊函數使得穿刺點不能動。
在人們還在順著Brügmann教授的思路往下思考時，由加州理工學院和康奈爾大學聯(lián)合的數值相對論小組正在廣義調和坐標方程形式和譜方法兩個問題上苦苦掙扎。數值相對論界有名的“獨狼”——Frans Pretorius于2003年從Matthew Choptuik處學成畢業(yè)到加州理工學院進行博士后研究。
獨狼就是獨狼，一個人把調和坐標方程形式套用到他在博士期間發(fā)展的自適應網格細化軟件平臺上。2005年秋，數值相對論學家們從睡夢中被人叫醒，“雙黑洞計算穩(wěn)定性問題解決了！”。人們都覺得難以置信，Pretorius卻報告說：“我就這么一做，沒問題??！”。
當人們還在思考到底發(fā)生了什么的時候，2006年，美國德克薩斯大學布朗斯維爾分校的數值相對論小組（現(xiàn)在的羅切斯特理工學院數值相對論小組）放棄Brügmann教授固定黑洞穿刺點的想法獲得了成功。同年美國國家航空航天局NASA的數值相對論小組不小心忘記了固定穿刺的黑洞，計算出與德克薩斯大學布朗斯維爾分校類似的結果。所以他們兩家同時宣布，“我們也可以穩(wěn)定計算雙黑洞了！”。和Pretorius一起，他們三個小組計算所得的引力波波形對比結果如圖7所示。明眼人皆可看出，異常吻合！
圖7：最早實現(xiàn)雙黑洞穩(wěn)定計算的三個數值相對論小組所得引力波波形對比。圖片來源：文獻[8]
后來人們認識到，方程形式、坐標條件、網格精細程度、邊界條件、黑洞奇點處理方式和數值計算方法等都對愛因斯坦場方程數值計算的穩(wěn)定性起到了關鍵性作用。而且這些問題的處理方式需要合適地搭配才能保證數值相對論計算的穩(wěn)定性。
直到現(xiàn)在，也沒有嚴格的數學理論給出數值相對論穩(wěn)定計算的充分必要條件，更別提數值相對論的計算收斂性理論。所以有人會說數學意義的數值相對論還沒有被建立起來，甚至可以說還沒有開始。但無論如何，世界上已經有十多個數值相對論小組能夠穩(wěn)定計算雙黑洞和其他雙星問題，并給出相應的引力波波形。
實際上，從數值相對論的“能算”到引力波波形模板的建立，中間還隔著遍歷參數空間的問題，僅帶有軌道進動的圓軌道雙黑洞模板由15個參數來描述，如果是帶有軌道偏心率或者包含中子星的模板則需要更多參數來描述，人們把這個問題稱為數值相對論問題中的維數災難。
同時，圖3中索恩理想化的三段相接波形的想法在實際數據處理中也會帶來若干問題?，F(xiàn)在，人們對這個問題的處理方式是，在深度理解數值相對論結果的基礎上，建立半解析波形模型，類似一種“模糊”處理，將精度降低。目前LIGO-Virgo科學合作組織主要使用的模板EOBNR和IMRPhenom系列都是這樣的波形模型。
由于數值相對論波形的計算成本過于高昂（一個波形需要約3萬核心小時的計算資源，一次引力波信號搜索需要約1000萬個波形，合計約10的11次方核心小時的計算資源），很難直接用于實際數據分析，生成速度更快且精度尚可的半解析模型就成為了目前引力波數據分析的首選。在建立這些模型的過程中，數值相對論的波形會被當做標準波形來作比對。更具體的建模細節(jié)參見文獻[9]。
從數值相對論的計算、“求穩(wěn)”到“模糊”處理，目前，人們還在繼續(xù)對雙星并合系統(tǒng)引力波波形做更深入的理解，力圖建立完備的雙星并合系統(tǒng)引力波模型和模板庫，為未來的空間引力波探測和下一代的地面引力波探測提供堅實的理論支持。
除了雙星系統(tǒng)，我們還可能探測到其它引力波波源，也許廣義相對論存在適用邊界，我們還需要其它的引力波波源和引力理論所描述的波形模板。數值相對論和波形分析建模的故事還會伴隨引力波天文學很長一段時間。
作者簡介：
· 吳仕超，北京師范大學天文系研究助理，KAGRA科學合作組織成員。研究方向是引力波波形模板構建和引力波數據分析。
· 趙天宇，北京師范大學天文系博士研究生。研究方向是數值相對論。
參考資料：
1.Matched filter https://en.wikipedia.org/wiki/Matched_filter
2.Abbott B P, Abbott R, Abbott T D, et al. GW150914: Implications for the stochastic gravitational-wave background from binary black holes[J]. Physical review letters, 2016, 116(13): 131102.
3.Misner C W, Thorne K S, Wheeler J A. Gravitation [M]. Macmillan, 1973.
4.Abramovici A, Althouse W E, Drever R W P, et al. LIGO: The laser interferometer gravitational-wave observatory [J]. Science, 1992, 256(5055): 325-333.
5.Cactus http://cactuscode.org/
6.Einstein Toolkit http://einsteintoolkit.org/
7.Brügmann B. Adaptive mesh and geodesically sliced Schwarzschild spacetime in 3+ 1 dimensions[J]. Physical Review D, 1996, 54(12): 7361.
8.Baker J G, Campanelli M, Pretorius F, et al. Comparisons of binary black hole merger waveforms[J]. Classical and Quantum Gravity, 2007, 24(12): S25.
9.蔡榮根, 曹周鍵, 韓文標. 并合雙星系統(tǒng)的引力波理論模型[J]. 科學通報, 2016 (14): 1525-1535.
10.Abbott R, Abbott T D, Abraham S, et al. GW190412: Observation of a binary-black-hole coalescence with asymmetric masses [J]. arXiv preprint arXiv:2004.08342, 2020.
11.Abbott R, Abbott T D, Abraham S, et al. GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object [J]. The Astrophysical Journal Letters, 2020, 896(2): L44.
12.Cao Z, Han W B. Waveform model for an eccentric binary black hole based on the effective-one-body-numerical-relativity formalism [J]. Physical Review D, 2017, 96(4): 044028.
13.曹周鍵, 都志輝. 數值相對論與引力波天文學[J]. 中國科學: 物理學, 力學, 天文學, 2017, 47(001): 55-72.
14.曹周鍵，愛因斯坦方程和數值相對論，《數學與人文》叢書第二十五輯，68頁，高等教育出版社 (2018).
15.Baumgarte T W, Shapiro S L. Numerical relativity: solving Einstein's equations on the computer [M]. Cambridge University Press, 2010.
16.Cai R G, Cao Z, Guo Z K, et al. The gravitational-wave physics [J]. National Science Review, 2017, 4(5): 687-706.
制版編輯 | Livan
