天文学中的数值模拟：数值相对论天体物理概述
数值模拟系列文章写到现在，其实一直还是在牛顿力学的框架下打转转。当然，之前介绍过的FLASH出于精确模拟天体爆发现象的考虑，包含了相对论性流体的计算模块，否则本人最近也不会出于工作需要去跟它一再较劲。但FLASH以及其他一些类似的程序涉及的只是狭义相对论，并没有将需要用到广义相对论的强引力系统涵盖进来。若要全面涵盖极端天文现象，对广义相对论进行数值求解是非做不可的。
说来本人在读研究生之初还是选修过Maurice van Putten开设的广义相对论数值方法相关课程的，不过因为那个学期的课程负担太重，最终没能抽出时间参加考试。实在对不起van Putten教授（连带对不起广义相对论课程的主讲教师外加某两册教材的作者）的是，因为平时的研究不大会涉及相关内容，课程中所学到的知识到现在基本上是原封不动地还了回去。因此撰写本文期间，自己也实在是心虚得很哪……
本文主要参考的是Lehner & Pretorius (2014)这篇综述。在介绍数值相对论天体物理的求解思路之前，首先要明确这一领域的适用范围。强引力系统分为两类，一是像致密星吸积系统这样的，虽然主天体引力场很强，但与它存在相互作用的吸积盘、喷流等结构自身的引力却微不足道，因此处理起来要简单得多。只要不去探讨致密星周边的细节，大可使用牛顿近似，只是在计算物质运动时要考虑狭义相对论的尺缩钟慢等效应。就算是双致密星系统，当两星间隔较远时，后牛顿近似也足矣。后一类系统就比较棘手了，比如并合中的双致密星，其中参与相互作用的天体都会显著影响整个系统的引力场构型，因此就必须用到广义相对论理论。
数值相对论的核心是爱因斯坦场方程，乍看形式很简单，Gμν = 8πGTμν。但你要知道，上式用张量形式表示的，把μ和ν的指标轮换一遍那就是16个方程，而且方程之间彼此还存在非线性的耦合，这可是比牛顿引力场方程或者流体的欧拉方程要抓狂多了。更要命的是，强引力场周围必须要考虑时空的弯曲，对于黑洞来说还存在视界这种让人头大的东西，于是数值求解广义相对论的过程中，像非相对论（包括狭义相对论）模拟那样一套网格用到底是万万不可能的，而必须一再地重新定义坐标系……
在实际计算中，用得最为广泛的是爱因斯坦场方程的两种变体——谐和方程或BSSN方程（好吧，后一种究竟是啥玩意本人也不清楚）。它们都属于双曲型偏微分方程，特征速度就是光速。由于广义相对论时空耦合的特性，时间和空间方面的解析都非常重要。空间离散化通常可以使用高阶差分或者所谓“谱分解（Spectral Decomposition）”（嗯，本人再一次不知道后者是啥意思）来实现，时间离散化倒要直观得多。为了保证求算结果满足场方程的耦合条件，同时便于数值上的具体实现，通常会采用“约束阻尼”法来修改运动方程，让修改后的方程产生合适的数值“阻尼”，不致于偏离真实解太远。至于什么自适应网格、并行处理之类，就都属于数值模拟界的国际惯例了，没有太多的特别之处。
如前所述，数值相对论天体物理主要用于密近双致密星系统的研究。这类系统是在引力辐射的驱动下演化的，通过辐射引力波损失轨道动能，星体间距逐渐缩小，最终走向并合，形成单一的克尔黑洞（对于低质量双中子星系统来说也可能会形成大质量中子星）。很自然地，随着引力波探测器的发展，现在这类研究很大程度上是为引力波研究服务的，首先要通过数值模拟来计算出并合产生的引力波信号形态，再以此为参考，在探测器的监测数据中大海捞针一般地进行搜索。另一方面，由于这类过程往往会导致短伽玛射线暴等极端过程，数值相对论对于天文学的理论研究来说也有着重要的意义。
通过求解广义相对论方程模拟出的双中子星并合形成短暴的过程，两星的质量分别相当于太阳的1.2与1.6倍。并合过程中，较小的中子星被撕碎，其上的物质累加到较大的中子星上。该过程最终的产物是一个内部温度高达数千亿度的致密天体（图中的红黄色区域），最终它将在数十毫秒内坍缩为黑洞，周围较冷的物质（蓝色和紫色）则形成吸积盘系统，产生短暴。（图片提供：R. Oechslin & H.-Th. Janka, MPA）
整个并合过程大致可以分为三段，分别是绕转（inspiral）期、非线性期以及余波（ring-down）期。绕转期的黑洞行为可以用后牛顿近似很好地描述，对于中子星来说再加上潮汐力即可；黑洞形成后的余波期用黑洞扰动理论也能应付。最棘手的是描述并合本身的非线性期，引力作用以及时空演化最为复杂，也是数值相对论发挥作用的时候。由于场方程存在以指数增长的非物理解，在计算中若规避不善，往往会导致程序无法继续运行，因此非线性期的数值模拟方法直到差不多10年前才被彻底解决，双星系统得以在计算机中自然演化下去。
非线性期可以继续细分为晚期绕转、冲入以及并合后早期三段。晚期绕转期间，双星仍旧保持轨道运动，但绕转周期陡降，轨道速度极高。在并合后早期，黑洞或超重中子星刚刚形成，由于变形严重，还不能用扰动理论应付，只能直接求解场方程。冲入阶段是非线性期最麻烦的部分，双星速度的径向分量剧增，甚至可能超过了轨道运动速度，最终带来了并合。这一过程与黑洞最内稳定圆轨道处发生的过程直接相关。
左：两个等质量黑洞并合期间发出的引力波辐射（图片来源：Buonanno et al. 2007）。图中的颜色表示Newman-Penrose标量与径向坐标之积，与引力波h+分量的二阶时间导数成正比，绿色、紫色与红色分别表示零值、正向与反向，上方一排是并合前的情况，左下是并合期间引力波辐射接近峰值的情况，右下是并合后的余波期。右：质量不等的并合黑洞产生的引力波辐射波形（图片来源：Taracchini et al. 2013）。
由于黑洞无毛（或者说只有三毛），黑洞系统的并合模拟要简单一些，具体计算只需两个黑洞的质量比、初始自转矢量、轨道偏心率与半长轴即可。这样的过程也不会产生电磁波辐射，因此模拟的目标也只是探讨引力波辐射的细节（参见下图），弄清并合期间损失的总能量（相当于黑洞最内稳定圆轨道的束缚能，依较大黑洞角动量不同，合较小黑洞静能的3.8%到42%）与角动量，以及最终形成的黑洞遭受的kick（无自转黑洞最终可达每秒175千米，超过了大型星系的逃逸速度）。除了服务引力波探测计划之外，黑洞并合相关的模拟工作还促使研究者去搜索并合黑洞的电磁场信号、黑洞kick相关的激波辐射、巨型黑洞相对星系核心的偏移等现象，并帮助人们进一步了解特大质量黑洞与星系并合和演化过程的关系。
由于参数的增加、物态方程的复杂性以及磁化物质的存在，涉及中子星的并合系统模拟起来比较麻烦。另外一旦双星轨道靠得足够近，后牛顿项变得不可忽略之后，潮汐力导致的中子星形变也是必须要考虑的。这方面的研究者除了关注引力波辐射的波形演化之外，又额外涉及了中子星坍缩的时标——如果初始系统由两颗中子星组成，并合过程一开始的产物是作为中间态的超重中子星，其最终坍缩为黑洞的过程所需时间除了与轨道构型有关，还取决于物态方程，不同物态引发的引力波辐射也不尽相同，因此搜索相关信号有助于澄清目前尚不明朗的中子星物质。可惜双中子星的并合产生不了太强的引力波信号，哪怕升级后的LIGO也只能对近距离的并合中子星进行，这种方法的运用还有待时日。另外与黑洞不同的是，并合期间的中子星潮汐瓦解事件也会为并合产物带来共振，这一点也能反映在引力波信号中。在引力波之外，与并合相关的短暴连带抛射物r过程同位素衰变驱动的kilonova研究不必多说，中微子辐射的探讨也是方兴未艾。
等质量双中子星并合过后反电子中微子的辐射分布图。（图片来源：Sekiguchi et al. 2011）。
如果是黑洞+中子星的系统，情况又有所不同。这里主要问题在于中子星星体在并合前即已瓦解的可能性（是否瓦解则取决于黑洞的最内稳定圆轨道半径和潮汐半径之间的关系，进而取决于黑洞自转与质量）。若系统中黑洞的最内稳定圆轨道大于潮汐半径，那么结果与双黑洞的并合相差不大；否则整个过程更接近于吸积，中子星先瓦解，再逐渐落入黑洞，并且会有相当一部分物质留在最内稳定圆轨道之外形成吸积盘。这种情况下的引力波辐射会大大减弱，其具体表现也能成为透露星体物态的线索。
除了致密星的并合，数值相对论的另一个用武之地是坍缩形成致密星的过程，如核心坍缩型超新星的爆发。爆发过程至今在数值上没能很好地解决，不过这并不妨碍人们研究致密星本身的形成以及相关的引力波辐射，也没有妨碍研究者辨认出星体爆发的若干元凶，如中微子加热、流体不稳定性、核燃烧等。由于方程复杂且涉及诸多过程，相关模拟是相当耗时的。这一问题的全面解决还有待时日，不过现有的结果已经说明，超新星相关的引力波辐射频率完全落在了现有地面探测器的观测范围内，若有河内超新星事件发生，找到它们的引力踪迹还是力所能及的。
未来的数值相对论天体物理将进一步扩展自身的研究领域，去探索含有白矮星或主序星的双星演化；同时进一步改善代码，提升计算速度；再有就是更好地衔接数值相对论与后牛顿近似和微扰法得出的结果。由于当前不同小组对类似问题给出的模拟结果比较接近，相对论的研究者对数值方法还是很有信心的。而随着LIGO升级工作的即将完成，数值方法也将在相对论的研究领域发挥更大的作用。
末了自己要加上几句对具体数值算法的看法。首先是黑洞的处理。据说现有程序大抵是用已知的黑洞解直接植入初始条件中，由此来回避数值计算过程中难以应付的奇点，但是对于坍缩期间黑洞的形成应该如何处理，本人还未见详细描述。这一点应该是数值模拟的一大难题，更有说是由于方程过于复杂，由超重中子星坍缩形成黑洞的过程计算起来非常不经济，当前甚少考虑。二是“下沉粒子（sink particle）”的使用。这个概念是在学习FLASH期间接触到的，不过本人应该暂时不会用到它，只是在此一提。这种粒子有着独有的生灭规则，引入它的目的是为了避免坍缩后期出现数值上的物质团块分裂（为保证数值稳定性，金斯长度必须用不小于32个格点来解析，低于此值则结果不可靠，但后期金斯长度很短，过度加密网格不现实）以及程序因自由下落时标短于时间步长而导致的中断。这种粒子的具体实现方法和特性并非本人当前关注的话题，只是想说，在坍缩星的模拟工作中，要经历流体到粒子再到奇点的过程，看上去人为性很强啊，总感觉还是没有那么的自然……
末了的末了，尽管综述中没有涉及，但《星际穿越》中的这个场景也应该属于数值相对论天体物理的应用吧，虽然该片本人没有看过，但这不就是黑洞吸积盘的模拟图像嘛，考虑了强引力场中光线的弯曲：
图片来源：Interstellar
嗯嗯，更加接近真实的图像应该是这样的，左为侧视右为俯视。图中考虑了相对论性速度运动的气体流的多普勒增亮效应，但没有考虑红蓝移。不过对于宽波段连续谱辐射来说，这红蓝移考虑与否也差别不大吧？
图片提供：Phil Armitage/Christopher Reynolds
