没有它们，给你引力波你也不认识！| 赛先生天文
引言：
自2015年LIGO探测到首个引力波事件GW150914以来，引力波天文学已取得长足进步。从寻找信号，到推断波源和族群性质，再到检验广义相对论，几乎所有引力波数据分析过程都离不开精确的引力波模板。那么，引力波模板是什么？在引力波模板构建中的过程中，数值相对论又扮演了至关重要的角色，它又是怎么发展起来的？
相对论和引力波大家一定有所耳闻，但数值相对论和引力波模板可能就触及到大家的知识盲区了。没有它们，就算引力波探测器探测到了引力波，我们也无法精确分析其携带的有用信息！本期赛先生天文，让我们认识一下引力波模板和数值相对论。
撰文 | 吴仕超、赵天宇（北京师范大学）
责编 | 韩越扬、吕浩然
何谓“引力波模板”
通常来讲，飞机和雷达之间是“猫和老鼠”般的竞争关系：雷达发射电磁波，并利用回波寻找飞机，而回波会因为信号衰减而变得微弱。于是，雷达探测飞机则会变成一个典型的弱信号、强噪声的问题。人们由此发展出了相应的匹配滤波（matched filtering）数据处理方法，以提高信号的信噪比[1]，使之更容易被探测到。雷达接收到的回波和其发射的电磁波具有基本相同的波形，所以雷达的匹配滤波可以直接实现。
而引力波的探测也存在相似的地方：从数据分析的角度来讲，引力波探测也是典型的弱信号、强噪声问题。相应地，人们自然也想到了通过使用匹配滤波的方法来增强引力波的探测能力。
为了实现匹配滤波，我们需要预先知道信号的形状，也就是波形（waveform）。由于引力波探测器是被动接收引力波信号，因此它的波形不能预先得知。对此，人们的应对办法是从理论上，针对各种引力波波源构造引力波的波形库（template bank），也就是所谓的引力波模板（template）。
因为实际的引力波信号，波形受波源参数影响且参数未知，所以引力波模板库通常需要覆盖整个可能的波源参数空间，一个模板库通常包含数百万个引力波模板。不准确的引力波模板会与实际引力波信号波形不吻合，损失匹配滤波信噪比，从而导致引力波信号丢失。所以引力波波源建模至关重要。
数值相对论又是什么
图1：GW150914在引力波探测器中的波形及轨道演化示意图，图片来源：文献[2]。
在引力波被实际探测到之前，双星并合系统是科学家预测的最有可能的波源，同时也是引力波天文学研究中最重要的引力波波源之一。因此，它也是人们首先想到要构建引力波模板的系统。
双星系统中的双黑洞系统不涉及复杂的潮汐相互作用，是理论描述上最“干净”的引力波波源系统。为了获得双黑洞系统的引力波波形，人们需要也只需要求解爱因斯坦场方程即可。但爱因斯坦场方程的复杂性是出了名的！一两千页的专著《爱因斯坦场方程的精确解》不干别的，只为了列举到现在为止已知的爱因斯坦场方程的解。
图2：非常“简单”的爱因斯坦场方程
更“可惜”的是，一两千页中与天文系统有对应关系的解并不多。换句话说，爱因斯坦场方程的解析求解非常困难。爱因斯坦本人也深谙其引力场方程的困难程度，他本人以及后来非常多的科学家们发展出弱场、低速、微扰等近似爱因斯坦场方程的计算分析办法。但对于引力波辐射最强的双星并合过程，双黑洞的相对速度在临近并合前会越来越快，最后接近光速（如图1），引力场强超过地球附近场强6个量级，广义相对论效应强到一切近似计算方法都失效。
于是，人们想到了数值方法求解爱因斯坦场方程，也就是通常人们说到的数值相对论（numerical relativity）。基于这个设想，基普·索恩（Kip Thorne）早年画出了双黑洞并合引力波波形的示意图（如图3）。
图3：索恩设想的双黑洞引力波波形（下方坐标轴内），图片来源：Kip Thorne。
既然都“沦落”到数值求解爱因斯坦场方程了，无非就是写写程序，花费些电力资源，喝喝咖啡等着波形出来就可以了。这也是人们最早对数值相对论的预期。带着这个想法，Susan G. Hahn和Richard W. Lindquist最早在1964年就开始尝试数值求解双黑洞的爱因斯坦场方程。
Hahn是计算数学大家Peter Lax的学生，当时他在IBM上班。Lindquist是有名的广义相对论学家。等他们写好程序，启动计算机，坐下来准备喝咖啡的时候，预期以外的事情发生了：程序运行十多步就崩溃了, 输出了一大堆NaN（编者注：NaN在计算机科学中代表一堆意义不明的值），没有得到任何有物理意义的结果。两人瞪直眼睛，为什么？这件事开启了数值相对论学家们长达近半个世纪的苦思。这个问题被人们称之为数值相对论的稳定性问题。
如何“求稳”
广义协变性原理是广义相对论的理论基础之一。该原理从著名的马赫原理发展而来，爱因斯坦也以之为傲。该原理使得爱因斯坦场方程以张量形式出现，与数值计算所需的偏微分方程仅一纸之隔，人们需要选择坐标系（图4）把爱因斯坦场方程写成偏微分方程的形式。
20世纪70年代，普林斯顿大学的Kenneth Eppley和德克萨斯大学奥斯汀分校的Larry Smarr对于爱因斯坦场方程数值计算的坐标选择做了深入研究，发现坐标选择会严重影响数值计算的稳定性。他们针对一些特殊物理情形, 充分发掘其内在对称性, 发展了特殊、精致的坐标以供数值计算使用, 取得了不错的效果。但缺点是这些坐标选择不具有通用性, 无法对非对称性时空做一般性推广。他们还提出了自动搜索对称性的坐标条件等, 但发现还是不能保证计算的稳定性。
20世纪80年代，Tsvi Piran对一些简单的模型化问题给出了初步的数值计算结果。这种利用时空对称性选择坐标系的方式在引力/规范对偶数值相对论问题中也被经常使用。但可惜的是，我们面对的真实引力波波源系统往往什么对称性都没有，比如双黑洞并合系统。
图4：Kruskal-Szekeres坐标下的史瓦西黑洞，图片来源：文献[3]。
直到20世纪90年代初，LIGO的初步硬件（initial LIGO）建设正式启动[4]。眼看引力波探测的数据处理很快就要开始，人们对引力波模板的需求迫在眉睫。人们急红了眼，猜测数值相对论所遇到的不稳定性问题会不会是由于计算所用网格不够精细导致的？为了把计算网格铺得足够细，就需要用非常多的计算节点并行计算。为了把复杂的爱因斯坦场方程和复杂的并行计算综合到一起，就需要许多人的协作。
综上，美国多个研究所和高校联合组成“双黑洞大挑战联盟”（The Binary Black Hole Grand Challenge Alliance，图5）, 计划构建大型数值相对论软件以实现对双黑洞问题的数值模拟。一款名为Cactus的软件[5]就是在这个背景下诞生的，而如今数值相对论领域最成熟的开源软件Einstein Toolkit[6]就是基于Cactus发展起来的。但该联合团队的研究结果却表明，数值相对论的不稳定性问题并非只是计算分辨率的问题。双黑洞大挑战联盟后来也就不了了之了。
图5：上世纪90年代双黑洞大挑战联盟的介绍网页（http://www.crpc.rice.edu/CRPC/demos/BlackHole/）
数值计算不稳定是不是因为所计算的偏微分方程本身就不稳定呢？针对这个想法，若干的数值相对论工作者从偏微分方程的角度出发分析计算所用偏微分方程的双曲性质。人们发现，虽然各种各样的用来数值计算的偏微分方程形式都来自于原来的爱因斯坦场方程，但各自的偏微分方程性质却千差万别。
基于对称双曲性质的要求，加州理工学院和康奈尔大学联合数值相对论小组提出了广义调和坐标形式的计算用爱因斯坦偏微分方程。这些方程形式在2003年之前就被提出来了，但数值相对论的不稳定问题依然进展缓慢。
Bernd Brügmann教授1996年博士毕业，论文是圈量子引力。圈量子引力界的著名物理学家Martin Bojowald和Thomas Thiemann都是他的师弟。酷爱中国围棋的Brügmann教授喜欢通过简单的组合得到复杂事物，这促成他从圈量子引力到数值相对论的研究兴趣转变。
不同于双黑洞大挑战联盟关于计算网格问题的想法，Brügmann教授把双黑洞问题当作多尺度物理问题来看待，他最早把科学计算领域的自适应网格细化（adaptive mesh refinement）技术引入到数值相对论（图6），建立了现在称之为BAM的数值相对论软件。Brügmann教授喜欢合作，但只限于十人以下的合作。这就是为什么他在德国马普引力所工作期间拒绝了把BAM引入到Cactus当中。所以Cactus的自适应网格细化大科学计算平台部分是2003年之后才由Eric Schnetter发展起来的，现在被称为Carpet。
1999年，Brügmann教授基于他发展的自适应网格细化技术把双黑洞的数值计算时间推进到了35个时间单位（数值计算把双黑洞问题无量纲化后的时间单位）。以现在计算双黑洞约1000个时间单位的标准来看，35个时间单位微不足道，但在当时这已经是不小的进展了。到2000年，Brügmann教授和合作者一起把计算时间推进到了50个时间单位。
图6：双黑洞计算中的自适应网格细化，越靠近黑洞的区域网格越密集。图片来源：文献[7]
2000年对引力波发展是个特殊的年份，一方面LIGO初步硬件基本建成，眼看就可以开展数据采集了；另一方面，Kip Thorne这一年60岁大寿。对引力波探测的极度乐观和对数值相对论发展的极度悲观使得索恩在其生日宴会上说“很可能引力波的成功探测比数值相对论的突破还要早日实现”。当然，后来的发展和索恩的这句话正好相反。
2000年以后，数值相对论继续发展。Brügmann教授与其合作者相继提出了黑洞穿刺的处理方式和针对BSSN计算方程形式的特殊坐标选择条件。这些技术把双黑洞计算时间推进到150个时间单位。也许做事精细是德国人的特点，Brügmann教授一直“固执”地认为黑洞穿刺点存在坐标奇点，处理坐标奇点的特殊函数使得穿刺点不能动。
在人们还在顺着Brügmann教授的思路往下思考时，由加州理工学院和康奈尔大学联合的数值相对论小组正在广义调和坐标方程形式和谱方法两个问题上苦苦挣扎。数值相对论界有名的“独狼”——Frans Pretorius于2003年从Matthew Choptuik处学成毕业到加州理工学院进行博士后研究。
独狼就是独狼，一个人把调和坐标方程形式套用到他在博士期间发展的自适应网格细化软件平台上。2005年秋，数值相对论学家们从睡梦中被人叫醒，“双黑洞计算稳定性问题解决了！”。人们都觉得难以置信，Pretorius却报告说：“我就这么一做，没问题啊！”。
当人们还在思考到底发生了什么的时候，2006年，美国德克萨斯大学布朗斯维尔分校的数值相对论小组（现在的罗切斯特理工学院数值相对论小组）放弃Brügmann教授固定黑洞穿刺点的想法获得了成功。同年美国国家航空航天局NASA的数值相对论小组不小心忘记了固定穿刺的黑洞，计算出与德克萨斯大学布朗斯维尔分校类似的结果。所以他们两家同时宣布，“我们也可以稳定计算双黑洞了！”。和Pretorius一起，他们三个小组计算所得的引力波波形对比结果如图7所示。明眼人皆可看出，异常吻合！
图7：最早实现双黑洞稳定计算的三个数值相对论小组所得引力波波形对比。图片来源：文献[8]
后来人们认识到，方程形式、坐标条件、网格精细程度、边界条件、黑洞奇点处理方式和数值计算方法等都对爱因斯坦场方程数值计算的稳定性起到了关键性作用。而且这些问题的处理方式需要合适地搭配才能保证数值相对论计算的稳定性。
直到现在，也没有严格的数学理论给出数值相对论稳定计算的充分必要条件，更别提数值相对论的计算收敛性理论。所以有人会说数学意义的数值相对论还没有被建立起来，甚至可以说还没有开始。但无论如何，世界上已经有十多个数值相对论小组能够稳定计算双黑洞和其他双星问题，并给出相应的引力波波形。
稳定之后的“模糊”处理
实际上，从数值相对论的“能算”到引力波波形模板的建立，中间还隔着遍历参数空间的问题，仅带有轨道进动的圆轨道双黑洞模板由15个参数来描述，如果是带有轨道偏心率或者包含中子星的模板则需要更多参数来描述，人们把这个问题称为数值相对论问题中的维数灾难。
同时，图3中索恩理想化的三段相接波形的想法在实际数据处理中也会带来若干问题。现在，人们对这个问题的处理方式是，在深度理解数值相对论结果的基础上，建立半解析波形模型，类似一种“模糊”处理，将精度降低。目前LIGO-Virgo科学合作组织主要使用的模板EOBNR和IMRPhenom系列都是这样的波形模型。
由于数值相对论波形的计算成本过于高昂（一个波形需要约3万核心小时的计算资源，一次引力波信号搜索需要约1000万个波形，合计约10的11次方核心小时的计算资源），很难直接用于实际数据分析，生成速度更快且精度尚可的半解析模型就成为了目前引力波数据分析的首选。在建立这些模型的过程中，数值相对论的波形会被当做标准波形来作比对。更具体的建模细节参见文献[9]。
结 语
从数值相对论的计算、“求稳”到“模糊”处理，目前，人们还在继续对双星并合系统引力波波形做更深入的理解，力图建立完备的双星并合系统引力波模型和模板库，为未来的空间引力波探测和下一代的地面引力波探测提供坚实的理论支持。
除了双星系统，我们还可能探测到其它引力波波源，也许广义相对论存在适用边界，我们还需要其它的引力波波源和引力理论所描述的波形模板。数值相对论和波形分析建模的故事还会伴随引力波天文学很长一段时间。
作者简介：
· 吴仕超，北京师范大学天文系研究助理，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.蔡荣根, 曹周键, 韩文标. 并合双星系统的引力波理论模型[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.