广义相对论数值解探索

1915 年，爱因斯坦提出了广义相对论，将引力的相互作用与时空弯曲的几何效应相联系。通过求解广义相对论的爱因斯坦引力场方程，我们能够深入研究极端致密天体的复杂运动规律，并揭示宇宙的演化过程。然而，精确求解爱因斯坦引力场方程却是一项极其艰难的任务。甚至对爱因斯坦引力场方程的数值求解方法也一直是物理学、天文学以及计算数学领域中的重要研究课题。数十年来，几代物理学家和数学家对这一问题进行了不懈探索，直至本世纪初才取得重大突破。本文将重点回顾过去半个多世纪以来对爱因斯坦引力场方程数值求解的探索历程。

引力相互作用是自然界四种基本相互作用之一。历史上，最早出现并相对完备的引力理论是17世纪牛顿提出的万有引力理论，它系统地描述了引力与距离平方成反比的规律。在此基础上，经由欧拉、拉格朗日、拉普拉斯等人的发展与完善，牛顿引力理论在描述地球、月球以及太阳系其他行星的运动时取得了令人瞩目的成功。然而，牛顿引力理论在面对高速运动和极端强引力环境(如中子星、黑洞等)时，其局限性逐渐显现。1915 年，爱因斯坦提出的广义相对论，成功解决了这一难题，为现代引力理论奠定了坚实的基础。广义相对论不仅能够描述高速运动物体在极端强引力系统中的相互作用，还为我们提供了一个完整的引力理论框架。

在广义相对论中，引力相互作用被解释为弯曲时空的几何效应。具有大质量的天体会使周围的时空发生弯曲，从而产生一系列可观测的现象，包括引力红移、光线偏折以及雷达信号回波延迟(即夏皮罗延迟)等。更为重要的是，广义相对论预言了宇宙的膨胀、黑洞的存在以及引力波的形成与传播，这些预言无一例外都已通过现代天文观测得以验证。时至今日，广义相对论是最为成功和完备的引力理论。

在广义相对论中，描述引力场动态变化的方程是爱因斯坦引力场方程

R_μν– (1/2)R∙g_μν=(8πG/c⁴ )∙T_μν          (1)

该方程是张量型方程，包含了描述时空曲率的几何量和描述物质的能量动量张量。能量动量张量全面反映了物质的能量密度、能量流以及应力等物理特性。若已知物质的能量动量张量，即可通过求解爱因斯坦引力场方程得到相应的时空度规以及曲率，从而进一步得到弯曲时空所产生的可观测效应。物质的能量动量张量受制于时空几何，弯曲时空的几何效应又会反过来影响时空中物体的运动。正如物理学家约翰·惠勒所指出的“物质告诉时空如何弯曲，时空告诉物质如何运动”。这一理论在天文学中具有重要的应用价值：通过求解爱因斯坦引力场方程，人们能够深入研究任意天体的运动行为，乃至探究整个宇宙的演化过程。

爱因斯坦引力场方程是自然科学中最为复杂的方程之一，其求解具有较高的难度。在已知的解析解中，首先被发现的具有里程碑意义的解是1916年得到的史瓦西解，该解描述了球对称引力源外部的引力场特性。基于史瓦西解，科学家能够精确计算出太阳系因广义相对论效应导致的观测现象，例如水星近日点进动以及星光在太阳附近的偏折等。这些现象均已超出了牛顿理论的解释范围，并被观测数据检验。另一个有划时代意义的解析解是1963 年提出的克尔解，它描述了具有轴对称性的旋转引力源外部的引力场特性。克尔解在现代天文学中具有广泛的应用，特别是在研究星系中心超大质量黑洞的性质时，它为理解黑洞的吸积过程和黑洞影像的形成提供了理论基础。这些解析解的发现和应用，不仅深化了我们对广义相对论的理解，也为探索宇宙中的极端天体现象提供了重要的理论工具。

由于爱因斯坦引力场方程的高度复杂性，能够解析求解的天体系统极为有限，仅包含史瓦西解、赖斯纳-诺德斯图姆解(RN解)、克尔解以及克尔-纽曼解等。然而，这些解几乎仅适用于静态或稳态的引力系统，即引力场本身不随时间发生变化。对于那些随时间快速演化的动态引力系统，如致密星体在自身引力作用下的坍缩过程，或是由两个致密天体组成的双星系统，直接解析求解爱因斯坦引力场方程变得异常困难，甚至不现实。对于这些动态变化的引力场，必须依赖于计算机进行高精度的数值模拟。

20 世纪40 年代，电子计算机的发明为科学研究开辟了新的可能性，这一技术迅速在各个学科中得到广泛应用。在引力物理学领域，研究人员开始尝试利用计算机程序计算爱因斯坦引力场方程的数值解，并通过数值模拟研究引力场动态演化过程。这种研究方法起始于20 世纪60 年代。经过数十年的持续发展，逐渐发展成为一个独立学科——数值相对论。

数值相对论要解决的是动态变化引力场的数值计算。也就是说，给定某一个时刻t₀的引力场(称为初值)，通过数值计算得到其后时刻t > t₀的引力场。这一理解在其他经典物理理论中不存在任何困难。但在相对论中，由于时间和空间的概念是统一的，不可分割，这使问题变得比较复杂。通俗地说，在广义相对论中，要标定物理系统的“时间”，需要在时空中做一系列切片处理，将4 维时空切成一系列3 维超曲面。我们要求切出的3 维超曲面是类空的，即任意一张3 维面上的不同位置之间不存在因果联系。这样便可以在这一系列类空超曲面上标定时间，即对Σ₁、Σ₂、Σ₃、…、Σ_n标定时间t₁、t₂、t₃、…、t_n。这样的将4 维时空进行切片并标定时间的方式被称为时空的3+1 分解，如图1 所示。1959 年，Arnowitt、Deser 和Misner 首先提出了一种将时空3+1分解的方式，被称为ADM分解

ds²= g_μνdx^μdx^ν= -α²∙dt²+ h_ij∙( βⁱdt + dxⁱ) ∙( β^jdt + dx^j)  (2)

经过这样的3+1 分解后，将(2)式表示的时空度规g_μv带入到爱因斯坦引力场方程中，便可得到用于数值计算的方程式，即h_ij、α、βⁱ关于时间变量t 的偏微分方程。给定某一初始时刻t₀的初值g_μv，原则上我们可以计算得到其后时刻t >t₀的引力场，并研究随时间的变化。然而，在爱因斯坦引力场方程中，时空度规g_μv的各个分量并不是互相独立的，10 个分量中只有6 个是动力学变量，剩下的4 个与坐标系选取有关。在3+1 分解中的α 和βⁱ都是与坐标选取相关的变量，它们被称为规范变量。除此之外，在实际计算中，为了使数值计算更加方便，还会引入辅助变量(如超曲面Σ 的外曲率K_ij以及其他变量)，将爱因斯坦引力场方程从关于时间t 的二阶偏微分方程降阶为一阶偏微分方程。这样处理后，最终可以得到数值计算中真正求解的方程，它们表面上看和爱因斯坦引力场方程(1)式差异很大，但从数学本质上是等价的，这些方程称为数值相对论的计算方程形式。规范变量选取的不同，引入辅助变量的方式不同，会导致数值相对论中存在不同的计算方程形式。在现代数值相对论程序计算中，常见的计算方程形式有GH(Generalized Harmonic) 方程形式、BSSN(Baumgarte-Shapiro-Shibata-Nakamura) 方程形式、CCZ4(Conformal Covariant Z4)方程形式、Z4C 方程形式等。

图1 时空的3+1 分解

从理论层面上来讲，数值相对论中使用的各类计算方程形式与原本的爱因斯坦引力场方程(1)从数学上是完全等价的，因此计算结果应完全一致。规范变量的选取仅仅提供了不同坐标系(或不同观测者)来测量引力场，它们不应影响最终的物理结果。然而，在实际的数值计算中，由于计算机计算不可避免地存在数值误差，同时不同偏微分方程的形式在数值求解过程中采用不同的数值方法，其在数值计算中的行为也会有所差异(有些计算收敛迅速，有些计算收敛慢，甚至在有些情况下可能超出数值方法的收敛域而无法逼近真实解)。因此，在数值相对论中，不同计算方程的选取和规范变量的选取往往会导致数值计算结果出现显著差异。

事实上，只有少数特定的选择能够获得良好的计算结果，而大多数随意的选取可能导致数值不稳定，甚至使程序崩溃。在20 世纪中叶，数值相对论中最棘手的问题是数值稳定性问题。研究者们在那个时期几乎无暇顾及计算效率和计算精度问题，更难以得到有意义的物理结果。

1964 年，Hahn 和Lindquist 开始尝试数值求解双黑洞系统的爱因斯坦引力场方程，但开始计算后，程序很快就崩溃了，没有得到任何有物理意义的结果。这种情况被人们称为数值不稳定，这是由于数值计算的累计误差快速增长，最终导致数值计算中出现非数(NaN)。20 世纪70年代，Eppley 和Smarr 对于爱因斯坦引力场方程数值计算的坐标选择做了深入研究，发现坐标选择会严重影响数值计算的稳定性。他们针对一些特殊物理情形，充分发掘其内在对称性，发展了特殊和精致的坐标供数值计算使用，并取得了良好的效果。但缺点是这些坐标选择不具备通用性，无法对非对称性时空做一般性推广。

早期，人们认为数值不稳定性源自于计算使用的网格分辨率不够高，只要网格分辨率足够高，一定能够获得稳定的数值计算结果。20 世纪90 年代，受当时激光引力波干涉观测计划(LIGO)的驱使，美国多个研究所和高校联合组成“双黑洞大挑战联盟”，计划构建大型数值相对论软件以实现对双黑洞演化问题的数值模拟。一款名为Cactus 的软件就是在这样的背景下应运而生。但该联合团队的研究表明，数值不稳定性不能简单地归咎于计算使用的网格分辨率不够高。进一步的研究结果显示：计算网格分辨率、数值计算方法、爱因斯坦引力场方程的计算方程形式、边界条件等因素都会对计算的数值稳定性产生影响，但这些因素都不能唯一地决定数值稳定性。到2000 年，LIGO的硬件设施建设已基本完成，但数值相对论的数值稳定性问题并没有得到解决。在这种情况下，基普·索恩悲观地说：“很可能引力波的探测比数值相对论的突破还要更早实现”。

2005 年，数值相对论在稳定性问题上取得重大突破，来自加州理工学院的Pretorius 首次完成了双黑洞合并全过程的数值模拟，为这一领域的发展奠定了重要基础。随后的2006 年，德克萨斯大学布朗斯维尔分校和美国国家航空和航天局(NASA)的数值相对论团队也各自独立地突破了数值稳定性难题。三组研究人员虽然采取了不同的计算方法，但他们对双黑洞合并过程中引力波波形的模拟结果显示出极高的一致性。这一重要发现不仅验证了数值相对论在强场动态过程中的有效性，同时也为引力波天文学的发展提供了可靠的理论支撑。受到以上工作的启发，自2006 年以后，多个研究机构的数值相对论团队也分别陆续解决了数值稳定性问题，并发展出多个数值相对论计算程序。从此，数值相对论和引力波波形建模进入新时代。在数值稳定性问题得到突破之后，数值相对论的重大发现是双黑洞合并中的引力波辐射会导致最终形成的大黑洞出现反冲。

2006 年以后，数值相对论的发展集中在以下几个方面。第一个方面是对致密双星系统的合并过程进行高精度计算，以帮助构建引力波波形的理论模版，这些引力波模版对于LIGO和Virgo 的地面引力波探测至关重要，如图2 所示。

图2 数值相对论在引力波探测中的作用

在LIGO 对首例引力波事件GW150914 的观测中，数值相对论提供的引力波波形模版扮演了重要的角色。第二个方面是寻找新的计算方法，进一步提高数值计算精度、数值稳定性，并降低所需的计算资源。经过近20 年发展，目前已经有一些方式可以将计算精度提高1~2 个量级(如使用CCZ4，Z4C 计算方程)，并将所需的计算资源降低2 个数量级(如使用高效的弯曲坐标网格计算)。第三个方面是将数值相对论拓展到更多的物理系统(除双黑洞之外的其他致密双星系统、引力坍缩系统、宇宙学系统)以及更广泛的理论模型中(这涉及其他的对广义相对论的替代理论，例如f(R)引力、标量张量引力、爱因斯坦-高斯-博内引力、爱因斯坦-麦克斯韦-标量耦合引力、全息引力理论等)。

在数值相对论中，我们可以从几个层次来介绍数值相对论的算法。第一个层次是如何从爱因斯坦引力场方程出发，根据问题约化出数值计算所需的计算方程形式；第二个层次是针对具体计算方程形式，设置相应的数值方法和边界条件；第三个层次是选定适当的计算网格和并行算法，以完成数值计算。

在数值相对论中，通常有两种方式对爱因斯坦引力场方程进行约化，以得到数值计算所需的方程式。一种方式是将爱因斯坦引力场方程约化为柯西初值问题，另一种方式则是约化为特征初值问题。把爱因斯坦引力场方程约化为柯西初值问题的第一步是对时空进行3+1 分解；第二步是挑选动力学变量(如超曲面Σ 的度规h_ij、外曲率K_ij和其他变量)，得到它们关于时间t 的偏微分方程；第三步是指定规范变量(如α、βⁱ和其他规范变量)以及它们随时间的变化关系。将爱因斯坦引力场方程化为特征初值问题时，后面两步的顺序与柯西初值问题相反：即在第一步时空3+1 分解之后，第二步首先指定规范变量和规范自由度(它们包含时空度规g_μv的组合)，通过数学处理使得g_μv变量呈现特定形式，第三步是将特定规范代入爱因斯坦引力场方程中，获得所有g_μv变量关于时间t 的偏微分方程。在数值相对论中，将爱因斯坦引力场方程约化为柯西初值问题后可以得到ADM 方程形式、BSSN 方程形式、CCZ4 方程形式、Z4C方程形式等，约化为特征初值问题后可以得到GH方程形式。

在偏微分方程的计算中，主流的数值方法有三类：有限差分方法、谱方法和有限元方法。有限差分方法是用网格上的数值差分来替代偏微分，这种方法的优势是物理图像和编程都比较简单，且有较高的并行计算能力，但缺点是计算精度不够高。谱方法是利用基函数将动力学变量展开，并将偏微分方程的计算转化为关于展开系数的方程计算。谱方法具有很高的计算精度，但基函数的全局性导致计算中需要全局数据交换，这限制了其并行计算能力。有限元方法是将计算的区域划分为一系列单元格，在单元格内使用谱方法，单元格与单元格之间使用类似于差分方法的处理方式。有限元方法可以结合谱方法的高精度和差分方法的高并行可扩展性。理论上预期用有限元法求解爱因斯坦引力场方程可以达到很好的并行可扩展性。但爱因斯坦引力场方程的弱方程形式很难构造，目前有限元方法在数值相对论的应用还不够多，有待进一步发展。

爱因斯坦引力场方程的数值求解中，我们还需要处理边界条件。首先，我们必须处理计算区域网格的外边界数据。同时，如果涉及黑洞的计算，还可能需要处理黑洞内部的内边界数据。黑洞内部存在时空奇点时，在奇点处时空曲率等各类物理量变为无穷大，这样会使得程序很快崩溃并停止运行。数值相对论中有两类方式处理黑洞的奇点问题，一种是直接剪除法，另一种是移动穿刺法。直接剪除法是将黑洞事件视界内部的区域直接排除在计算区域之外，这样剪除的边界(也就是黑洞视界)构成一个人为边界，然而如何给定此边界处的边界条件是一个非常微妙的问题。移动穿刺法的思路是在黑洞事件视界内部人为填充某种度规场，使得时空奇点不可见。常见的处理是在史瓦西黑洞的最大延拓时空取某一时刻的类空超曲面(它包含两个渐近平直区)，人为将其中一个渐近平直区的无穷远对应到黑洞内部的原点。该思路可以推广到任意黑洞情形。这样的操作把无穷大区域压缩到有限区域，把无穷远点压缩成一个“奇点”。但这样的“奇点”不是时空奇点，曲率并没有奇异性，它只是个弱奇点。人们把这个弱奇点称为穿刺点。在移动穿刺法中，随着时间演化，穿刺点和填充区域也会跟随黑洞一起移动，并自动填充度规场。

爱因斯坦引力场方程的数值计算中，也需要对计算用的网格进行特殊处理。在计算网格的处理中，可以利用自适应网格细化算法，在有限的计算资源实现较高的精确度。在数值相对论程序中，自适应网格细化通常是使用一系列分辨率由低到高的网格组成，远离引力源(如黑洞和其他致密天体)的远场区域使用分辨率较低的网格，靠近引力源的区域使用分辨率更高的网格，如图3 所示。自适应网格细化在处理多尺度问题中非常有效，它可以在物理量变化更剧烈的区域使用分辨率更高的网格，而在物理量变化较为平滑的区域使用分辨率较低的网格，在保持计算精确度的同时大大节省计算资源。双星系统合并释放引力波的问题便是一个典型的多尺度问题，双星本身和它们合并产物的大小和线度都局限在很小的区域(可将其记为距离单位)，但研究引力波的信号则需要关注远离引力源的远场区域(通常在1000 个距离单位或更多)。对于这样的问题，如果在全部区域使用分辨率很高的网格，则对计算资源需求已经远远超出了现在计算机的能力；如果在全部区域使用分辨力很低的网格，则不但计算精度很低，还可能导致数值不稳定。

图3 自适应网格细化示意图

数值相对论求解爱因斯坦引力场方程是典型的计算密集型问题，可使用图形处理器GPU来实现硬件加速，以提高数值计算效率。目前已有研究团队使用GPU 进行数值相对论的计算，实现了大约10倍的加速比(GPU+CPU混合计算对比纯CPU计算)。

值得注意的是，在数值稳定性解决的过程中，有几个措施发挥了至关重要的作用。首先是自适应网格细化，其次是计算方程形式的选取，最后还有涉及黑洞视界内部和黑洞奇点处的数据处理和边界条件。在计算方程式的选取中，只有特定的计算方程形式才能够实现数值稳定。这些方程格式可分为两大类，一类是GH方程形式(这是Pretorius在2005 年使用的方程)，另外一类是BSSN 方程形式(这是当今数值相对论软件中最常用的计算方程形式)，以及以BSSN 方程形式为基础的改进(如CCZ4 和Z4C 方程形式)。在黑洞内部和黑洞奇点处的数据处理中，如果用BSSN方程形式(或CCZ4、Z4C方程形式)进行计算，必须使用移动穿刺法处理黑洞奇点和黑洞视界内部的数据，才能实现数值稳定；反之，如果使用GH方程形式计算，必须使用直接剪除法处理黑洞视界内部的数据，才能克服数值不稳定性。其他的尝试，例如在BSSN方程形式计算中使用直接剪除法，或是在GH方程形式中使用移动穿刺法，都将导致数值不稳定性，这些数值不稳定的原因还是未解决的公开问题。是否有其他方式得到准确和稳定的数值结果，目前也是一个公开的问题。

自从21 世纪初双星合并过程的数值稳定性被解决，国际上已经有20 多个研究团队开发了求解爱因斯坦引力场方程和引力场演化的数值计算程序，这些计算程序多使用有限差分方法或谱方法进行数值计算，能够完成各种情形下致密双星(双黑洞、双中子星、其他致密天体等)系统演化的数值计算。这里，我们选取几个典型的数值计算程序进行介绍，更多软件可参考表1。

表1 目前国际上的数值相对论计算程序汇总(按照字母顺序排序)

Einstein Toolkit 是目前最广泛使用的数值相对论计算程序，由美国路易斯安纳州立大学为首的研究团队进行开发。它由Cactus软件平台发展而来，经过多年的开发和维护，有众多的功能和分支(每一个分支对应相应的Thorn)，能够对各种情形下的致密双星系统进行数值模拟。它拥有众多教程和完善的网站资源，吸引了广大用户使用其进行研究。

NRPy/BlackHole@Home 由美国伊利诺伊大学的研究团队开发，其目的是利用弯曲网格，对具有对称性的引力系统演化进行高效的计算。对具有球对称性和轴对称性的系统，用该软件可大大减少计算资源消耗。相比传统的立方网格，NRPy 使用的高效弯曲网格可以让数值模拟的内存使用减少1~2 个数量级，同时计算效率也大大提升。对于具有对称性的系统，使用NRPy 可以在桌面级电脑上完成计算，而无需像Einstein Toolkit 那样借助大型服务器和超级计算机。

AMSS-NCKU 是中国首个数值相对论计算程序，能对双黑洞和多黑洞系统的演化进行数值计算。该软件开发于2008 年，由当时中国科学院数学与系统科学研究院的曹周键和刘润球、清华大学的都志辉以及中国台湾成功大学的游辉樟等人合作完成。2012 年，曹周键及其合作者发展出Z4C计算格式，并在AMSS-NCKU 软件中实现数值计算，提升了计算准确度。2013 年，AMSS-NCKU软件实现了GPU计算功能，使得该软件的计算速度加速约10 倍。目前，AMSS-NCKU 软件有较高的运行效率，可以作为引力波波形建模的重要工具。最近，AMSS-NCKU还增加了Python 操作接口，这使得该软件对用户具有更好的可操作性。

自1915 年爱因斯坦提出广义相对论，确立了引力相互作用的现代理论，至今已有整整110 年的历史。根据广义相对论，引力相互作用的本质是时空弯曲，这种弯曲的程度及其形态完全由物质分布决定。在广义相对论中，描述引力场的核心方程是爱因斯坦引力场方程。原则上，通过求解爱因斯坦引力场方程，可以系统地研究任意天体系统的复杂运动现象、探究致密星体的坍缩过程、乃至揭示宇宙本身的演化规律。

由于爱因斯坦引力场方程的复杂性，只有非常有限的系统可以解析求解，更多无对称性和随时间变化的引力系统无法解析求解，只能利用计算机数值计算。数值相对论便是利用计算机数值求解爱因斯坦引力场方程并得到引力场随时间演化的信息。数值相对论从上世纪60 年代开始发展，至今已有60 多年的历程。20 世纪中，几代科研人员受困于数值稳定性问题，无法得到有意义物理结果，他们在黑暗中苦苦摸索。直到21 世纪初，数值稳定性问题才得以解决，这使得我们可以通过数值计算来研究引力场的变化、引力波的释放，并提取相关物理结果。如今，多个研究团队已经开发出若干数值相对论计算程序，用来对致密双星系统的合并过程进行数值计算，帮助构建引力波波形的理论模版。这些努力在LIGO 和Virgo 的地面引力波探测中发挥了重要作用。

尽管数值相对论在数值稳定性方面取得了进展，我们已经能够计算出双黑洞和双中子星合并过程产生的引力波波形，但这并不意味着数值相对论的研究已经结束。下一代引力波探测计划(例如LISA、太极、天琴等)的实施，对更低频率的引力波波源(如涉及极端质量比的系统)以及更高精度的引力波波形模板的需求日益迫切。研究这类极端质量比系统释放的引力波需要耗费大量的计算资源。这凸显出开发更高效算法、提高数值相对论的计算效率与计算精度的重要性。相信未来的数十年里，爱因斯坦引力场方程的求解和数值相对论将在科学研究中发挥重要作用。

本文选自《现代物理知识》2025年3期YWA编辑

更多精彩

《现代物理知识》

搜索微信号

mpihep

长按扫码关注我们