数值相对论

阿尔伯特·爱因斯坦在1915年系统地阐释了广义相对论。[8]与狭义相对论中的情形类似，时间与空间在广义相对论中统一为时空，而演化则服从爱因斯坦方程。爱因斯坦方程是一组耦合非线性偏微分方程组。在这个方程组诞生后的近一个世纪的时间里，物理学家仅仅解出了少数几个精确解。并且，其中许多的精确解是在高度对称的条件下，也就是说在方程被很大程度简化时，推导出来的，比如均一且各向同性宇宙的弗里德曼解。[9]

为了在更普遍的条件下研究爱因斯坦方程，数值相对论应运而生。数值求解爱因斯坦方程的必要前提是将四维时空还原为各自独立的三维空间与一维时间，也就是时空的“3＋1解构”。这个过程可以通过几种或简或繁的途径实现。这个问题的首个突破是由理查德·阿诺维特、斯坦利·德塞尔以及查尔斯·W·米斯纳于1950年代后期利用哈密顿力学沿着保罗·狄拉克给出的途径实现的。他们对于爱因斯坦方程的推演结果称作ADM形式。[10]出于技术原因，他们的结果并不适用于数值计算。它的双曲性很弱，因而在实际计算中很难使用。不过数值相对论中许多用到3＋1解构的求解方法还是与ADM形式紧密相关。这种解构会将爱因斯坦方程整理为具有一定初始条件的柯西问题。这种形式已经是计算机数值求解惯常解决的问题。[11]

时空的坐标也不能唯一确定。即使固定了起始超曲面上的坐标，但在穿越到邻近的超曲面时，时空坐标也会因不同的取法而不同。[lower-alpha 4]这是数值相对论的一个特征。这种度规上的自由度虽然并不会影响实际上的物理过程，但对过程的描述可能会因其改变。假设从起始的超曲面到邻近的超曲面，时间与空间坐标上分别移动了${\displaystyle \alpha }$与${\displaystyle \beta ^{i}}$。位移描述方式上的不同，也就是说${\displaystyle \alpha }$与${\displaystyle \beta ^{i}}$的取法不同，会导致推导出的方程组形式发生变化。适当的取法会为数值方程的数值求解带来潜在的好处，但坐标和度规许多看来“自然”的取法会导致数值解不稳定，并导致模拟失败。[12]

茨维·皮兰

在ADM形式发表时，计算机的算力并不足以应对这一计算。1964年，哈恩和林德奎斯特首次尝试数值求解爱因斯坦方程。[13]随后在1970年代，斯马[14][15]与埃普利[16]继续这方面的尝试。这些早期的尝试主要针对轴对称（“2＋1维”）时空的情况。茨维·皮兰写出了首个可以追踪辐射引力波的柱对称系统演化的程序。[17]皮兰在此项工作中为数值相对论中许多概念奠基，比如自由演化与约束演化。它们是推导起始数据随时间演化的两种途径。[18][19]对称性的引入降低了内存与运算速度的门槛，使得科学家可以通过当时的超级计算机解决一些问题。[17]

理查德·斯塔克与茨维·皮兰在1980年代早期首次对于真实的天体物理学问题做了数值计算。[20]他们首次计算了旋转黑洞在产生阶段辐射出的引力波。在这个结果发表后的近二十年的时间里，数值相对论少有结果发表。这可能还是因为计算机算力不够。1990年代，美国“大挑战”中的双黑洞项目组成功地模拟了两个黑洞迎头相撞的过程。他们利用了问题轴对称的特点简化了求解过程。项目组还在后期处理阶段通过这个数值解计算了事件视界。[21]

而在三维时空中数值求解爱因斯坦方程方面最早的尝试专注于不自转的史瓦西黑洞。这种黑洞静止且具有球对称性。对于验证数值相对论，这一问题十分理想：一是由于科学家非常了解这个问题的解析解，二是由于不自转的黑洞的数值解必然对时间收敛，三是由于这个问题还包含了数值模拟中一种非常复杂的模拟对象——黑洞中心的引力奇点。1995年，安尼诺斯等人发表了这一问题最早的结果之一。[22]在文章中，他们写道：

目前的计算机缺少足够的内存与算力以模拟解析的三维时空。这一点阻碍了三维数值相对论的发展。

在之后的几年里，随着计算机功能更为强大，几组研究人员发展了几种可以提升运算效率的技术。拉撒路组使用了早前通过非线性ADM方程组对黑洞融合模拟的结果，基于单黑洞线性摄动理论，得到了更为稳定的程序。[23]而从黑洞模型入手，又有两种可以绕过奇点问题的技术：切除法与穿刺法。[24]2005年，比勒陀利乌斯综合上述方法，采用了适当的坐标条件在模拟双黑洞方面取得了重大突破。[25]这几种新方法的稳定性在之后的几年里进一步提升，使得人们可以模拟两个黑洞在并合前彼此旋转几十以至几百转的情形。除此之外，计算流体力学中的自适应计算网格细化方法也被引入到数值相对论中来。[26]

继“大挑战”后，科学家又在拉撒路计划（1998－2005）中利用天体物理学结果对双黑洞融合做短期数值模拟。不过当时，在超级计算机上模拟这个系统演变的过程中，对于爱因斯坦方程的数值求解无论如何都得不到稳定的结果。研究者前后结合了多种近似方法（比如后牛顿方法和单黑洞扰动）来对全过程做数值模拟。[23]

拉撒路计划标志着双黑洞问题当时最为前沿的进展，给出了许多可以用于天体物理学研究的相当精准结果，比如引力波带走的能量与角动量，[27][28]不同质量的黑洞融合时发出的脉冲[29]以及融合后的黑洞的最终质量、动量和角动量。[30]科学家还在这项计划中计算了融合过程中释出的引力波的具体形式。这项结果对于引力波探测器的研制非常重要。同时，他们还预言黑洞撞击会瞬间产生宇宙中最为强大的能暴。在仅仅一秒内以引力波形式放出的能量足以大过银河系自产生以来其中恒星放出的能量总和。而引力波辐射只占系统总质量减量的很小一部分。[31]

切除法在1990年代后期提出，[32]是指在模拟演化过程中将奇点周围的事件视界“切除”掉。理论上，由于因果律以及视界本身的性质，[lower-alpha 5]这种方法不会影响到对于切除部分之外部分的演化的模拟。也就是说，即使没有解出黑洞内的情况也不会妨碍到外部的精确求解。人们可以通过设置适当的边界条件实现“切除”的过程。[33]

尽管切除法在实际应用中非常成功，但它也存在两个小缺陷。首先，在使用切除法时，人们必须审慎地选用坐标条件。尽管物理效应可能不会传播到视界外，但坐标效应可以。比如，在使用椭圆坐标条件时，黑洞内的网格变化会立即传播到视界外。[34]这意味着，人们需要采用特征速度小于光速的双曲型坐标条件（比如谐和坐标条件）。[35]其次，在黑洞运动时，所切除的部分也必须随黑洞一起运动。[33]

切除法在随后几年里随着校准条件的改进稳定性逐步提高。并且科学家也解决了切除部分随运算区域运动的问题。[36][37][38][39][40][35]借由切除法，双黑洞运动轨迹及融合过程的首个长期稳定的计算结果于2005年发表。[25]

在“穿刺法”中，解分为包含黑洞奇点的解析部分（也就是“刺”）[41]和奇点之外数值构造的部分。这种方法是由布里尔－林德奎斯特算法（针对起始时静止的黑洞）推广而得。[42]而如果再做一般化处理则可得到鲍恩－约克算法（针对起始时旋转且运动的黑洞）。[43]直到2005年，所有采用穿刺法的结果都需要刺在模拟过程中坐标位置不变。不过，彼此密联的黑洞自然会受到引力的作用运动。这意味着坐标系本身或会“伸展”或会“扭曲”。这会导致模拟的部分阶段出现数值不稳定。另外那种避开奇点的方法在模拟物质坍缩形成黑洞的过程中，由于也会采用这种坐标条件，同样也会造成类似效应。[44]

2005年，研究人员首次展示了刺在坐标系中运动的可行性，从而解决了这种方法早前存在的部分问题，使其可以用于精确追踪黑洞的长期演化。[25][45][46]通过选择合适的坐标条件并对奇点周围的物理场做粗略的估计，[lower-alpha 6]人们可以得到绕着彼此旋转的双黑洞的数值解，并能精确算出它们释放的引力辐射。[47]

自适应网格细化是一种在数值相对论萌芽前即已发展的数值算法。它是在1980年代由乔普提克在研究标量场坍缩的临界现象时引入的。在那种情形中，场参数处于黑洞形成与空间扩张的临界点。[48][49]由于使用球对称，这项工作针对的是一维情形。之后这种算法被推广到二维情形中。[50]二维网格法也用于研究非均匀宇宙[51][52]和史瓦西黑洞。[53]目前，自适应网格已经成为数值相对论的标准工具，用于研究黑洞与其他稠密系统的融合，以及这种事件产生的引力波。[54][55]