天文学中的数值模拟：超新星的引爆
最近读到马克斯·普朗克天体物理所的一篇新闻，称该所的研究者经过多年的努力之后，终于成功在计算机中自然引爆了源自大质量恒星死亡的核心坍缩型超新星。如果这一结果最终经受住了时间的考验，它必将成为天体物理研究史上的里程碑。超新星爆发，无论是核心坍缩型超新星还是起源于超重白矮星的Ia型超新星，从此都不再是纯粹的未知。当然，无论哪种超新星爆发都是很复杂的过程，其中牵扯的因素很多，相关的争议今后想必还会持续多年。
由于源自白矮星超过钱德拉塞卡极限的Ia型超新星的相关物理机制（相对）较为简单，只是以核聚变反应过程为主，其光变曲线的演化也完全由镍-56和钴-56同位素的衰变来驱动，解决得也相对比较早，大致可以追溯到10年之前，所以先介绍一下。无论什么样的模型，模拟Ia型超新星爆发的关键就是设法让局域的碳元素引爆扩散到整个星球上。不过早期的数值模拟已经证明，单纯考虑碳元素的暴震或者暴燃都不足以产生符合观测的外流物速度与超新星光谱。其中前者会在白矮星密度尚高的时候燃尽整个星体，但只能产生镍-56而缺乏中间元素，从而导致与实际不符的核合成过程；后者虽然需要经历星体膨胀，并在低密度环境下引燃涉及中间元素的核反应，但是总能量又低得不足以解释观测。所以为了实现真实再现Ia型超新星的目标，要么需要引入其他因素，要么需要将暴震与暴燃二者结合起来。
三维Ia型超新星的第一次自然引爆当然是FLASH的功劳。FLASH最初的开发目的正是去模拟以Ia型超新星为首的天体核燃烧。为了达到这个目标，它针对整个白矮星求解了描述可压缩流体的欧拉方程与描述自引力的泊松方程，并使用Khokhlov聚变俘获算法，通过假设核聚变火焰传播速度相当于最大层流火焰速度与湍流火焰速度，来确定碳元素向镁与其他元素的转化。芝加哥大学FLASH中心演化成功的第一种三维模型叫做引力束缚暴震（GCD）模型，之前在站里也提到过很多次了。简单地说，这种模型一开始是碳氧白矮星物质的积累导致星体核心附近引燃了核聚变反应，然后核聚变区域形成高温低密度泡状结构。泡体在浮力作用下一路上浮，最终速度超过音速，爆炸式地穿出星体。需要注意的是，泡体跨音速过程产生的弓形激波并不足以引燃撕碎星体的核聚变。相反，在泡体穿越星体表面的过程中，流出的物质在自引力束缚下从星球表面绕过，最终汇聚于突破点的对侧。这里因而受到物质流的挤压和加温，密度与温度上升，最终导致碳元素暴震的发生，并扩散到整个星球上去。当星体核心最终被暴震波扯碎后，超新星才进入同模膨胀阶段。模拟表明，暴震波可能并不是像早年的工作所认为的那样经由激波扩散，相反也可能是由于相对较缓和的密度梯度传播开去的。与传统的正中引爆相比，引力束缚暴震模型更为现实，因为白矮星核心是存在对流的，引燃核反应的位置不可能是严格的正中央。该模型可以产生与观测相符的爆发总能量以及超新星遗迹的特征结构，还成功给出了暴震发生前白矮星的核心密度与质量-密度分布之间的关系。之前在站里贴过核心引爆、高温泡上浮的模拟截图，这次就贴几张最终的暴震过程模拟图吧：
图片提供：DOE/NNSA/ASC/Alliance Flash Center/Univ. of Chicago
此外最近的研究表明，只要速度足够快，纯凭暴燃也不是不可以产生足够的能量来瓦解白矮星。而且在观测上，颜色偏红、光变曲线衰减速率较快的低光度Ia型超新星（也就是Iax型超新星）数量日益增加。虽然纯粹的暴燃难以产生正常Ia型超新星所需的镍-56同位素，不过它应付Iax型超新星应该是绰绰有余的。在此背景下，FLASH中心也对纯暴燃（PD）模型进行了三维模拟。它假设初始时超新星内部存在若干随机分布的“点火点”，再结合星体自转、碳氧组分比、对流形式等初始条件进行演化。所谓点火点，模拟中对应直径16千米的球形泡状区域，它们是白矮星质量接近或超过钱德拉塞卡极限，碳元素开始燃烧后的产物。点火点具体的形成数目存在争议，虽然前文所述的引力束缚暴震模型认为这样的区域只有一个，且略略偏离白矮星中心，但也有工作提出多个点火点是可能的。点火点内的物质达到了核统计平衡状态，温度是80亿开尔文左右。由于泡中高温燃尽物质表面包围着低温的反应燃料，由此产生的交换不稳定性会让反应区（也就是泡体表面）变得不稳定，由此诱发浮力驱动的湍动核燃烧。
通过改变点火点数目以及其他参数进行的一系列计算表明，随后核反应的速率是依赖初始点火点数目的。当存在100个左右的点火点时，白矮星初始的核反应速率较低，不过随后会迅速升高，对应星体膨胀速率的增长也存在一个由慢到快的演化过程；而点火点数量较多时，虽然初始核反应率很高，但随后并不会出现太大的增长，最终星体也难以全部燃尽。总的来说，初始点火点越少，则最终产生的能量与镍-56越多，爆发动能也越大；不同数目的初始点火点对应最终产物的分布也有所不同，数目较少则所得超新星的星等和光变曲线接近于Iax型超新星。
初始点火点数目不同的6种情况下纯暴燃模型Ia型超新星的演化过程，从上到下依次对应63、150、128、1100、1700与3500个点火点的情况。（图片来源：Long et al. 2014）
另一种Ia型超新星可能的爆发机制称作暴燃-暴震过渡（DDT）模型，很多人认为它是此类超新星最有可能的爆发途径之一，尤其对于单一白矮星的爆发（而非双白矮星并合）而言。在这一机制中，白矮星先是发生膨胀式的暴燃，随后在密度降到足够低，让湍流速度与层流速度相等的吉布斯长度小于核反应火焰的厚度时，转化为超音速的暴震燃烧，并释放出额外的能量。这种模型与纯粹的暴燃模型在初期演化没有任何本质区别，其星体最终膨胀的程度与核合成过程取决于早期的暴燃阶段。但是当前对关键的过渡过程的了解又很是匮乏，因此其中的不确定性很大。研究者也曾使用FLASH进行过相关的探讨，得出了与观测相符的结果，并给出了镍-56合成量与暴燃-暴震过程的过渡密度之间的关系。不过本人没能找到太多这方面的相关原始文献，只在此一提。
在FLASH中心的网站上找到了一张暴燃-暴震过渡模型的模拟图，不过不知出处。（图片提供：DOE/NNSA/ASC/Alliance Flash Center/Univ. of Chicago）
当然，在FLASH模拟Ia型超新星的前后，其他小组也使用不同的代码针对这一过程进行过模拟，并取得了不同程度的成功。不过它们一来名气没有那么大，二来本人也非其中的行家，因此还是暂且略过不谈。另外FLASH的结果也并非完全不存在争议，比如一开始在计算引力束缚暴震模型时就忽略了引燃核反应之前的星体对流，而纯暴燃模型也只考虑了无自转无对流的白矮星，这些出于简化计算的目的而引入的假设无疑是过度简化的，在一定程度上也削弱了结果的可靠性。Ia型超新星引爆的更多细节也还有待未来更真实的模拟来考察，不过现有的认识已经为白矮星爆发过程的了解带来不小的飞跃了。
与Ia型超新星相比，源自大质量恒星死亡的核心坍缩型超新星情况要更加麻烦一些。当质量大于太阳8倍的恒星核心有足够多的物质转化为铁，因而无法继续通过核聚变释放能量来对抗自身的引力的时候，核心坍缩就开始了。一旦星体核区密度达到原子核密度（也就是形成原中子星），中子与质子之间的排斥力足以对抗引力之后，坍缩骤停，星体核心由此发生反弹，强烈的反弹激波向外扫过外围仍在坍缩中的包层。按照最初的认知，这道反弹激波就应该足以逆转坍缩的过程，并将整个恒星瓦解，诱发超新星爆发。
但是自20世纪80年代初起，一系列数值模拟工作表明，如果只考虑流体过程，反弹激波的能量会被星周包层迅速吸收，最终在冲出星体之前，激波就会在核区内部停滞不前。由于无论如何改进分辨率、增加模拟的维度，所得结果都是类似的，因此为了引爆核心坍缩型超新星，激波冲出星周包层期间必然要由其他机制来补充能量，否则就无法应对真实观测了。又由于星体中心新生的中子星会散发出大量的高能基本粒子，又以中微子为甚；且超新星1987A爆发期间也确实有探测器接收到了成协的中微子，因此中微子相关过程成了探讨核心坍缩型超新星的首选，可惜相关的模拟也遭遇了重重阻碍。
核心坍缩型超新星爆发前，星体中心停止坍缩形成的反弹激波被中微子加热补充能量的过程示意。（图片提供：MPA）
在此期间还出现了一些旁门左道，比如之前在站里提过的快速自转强磁场恒星的爆发。它认为星体中心的物质在强烈的环向磁场的作用下射出了一对极向喷流，最终将恒星炸碎。但这一机制也只能在假设轴对称的二维情形下行得通，三维模拟中喷流初始的不对称扰动会迅速让其中的物质经由不稳定性扭曲扩散，无法突破受阻激波。看来引爆核心坍缩型超新星的正途还是要从中微子加热唤起激波复苏入手。
典型的超新星前身星在坍缩过程中，经由中微子释放的能量要比电磁辐射高上万倍。虽然中微子在通常环境下可以不受阻碍地从物质中穿过，但是在坍缩星这种极端的环境下，中微子光深还是不能忽略的，因此它们携带的部分能量理论上可以在星体深处堆积起来，并由受阻激波吸收，从而让激波重获冲出星周包层所需的动能。只要有1%的粒子能量可以转移到激波中去，超新星就可以爆发了。其实在很多年前，经由中微子过程诱发超新星爆发的计算机模拟就已经初步得以实现，还明确了星体物质的对流以及非径向运动在中微子加热过程中扮演的重要角色，但它们都是假设球对称的一维或假设轴对称的二维模拟。实际观测（尤其是对SN 1987A这样的近距离超新星进行的详尽观测）表明，超新星是高度不对称的。当对称性被破坏后，新的不稳定性浮出了水面，且原中子星在激波复苏后的吸积与外流过程也变得不可忽略，让模型复杂性大大增加。过去的三维中微子加热模拟并不能自然引爆核心坍缩型超新星，而是需要额外的干预，比如在达到某种条件后人为将星体炸开。
马普所针对20倍太阳质量恒星进行的坍缩与爆发模拟。（图片来源：Melson et al. 2015a）
在马普所之前，C. L. Fryer与M. S. Warren曾在10余年前实现过三维核心坍缩型超新星在中微子作用下的自然爆发，当时他们使用的是光滑粒子流体动力学模拟，再结合较为简单的“灰色”中微子输运模型。但是这种输运模型未免有过度简化之嫌，且从算法上就较易于引起爆发，因此不是太让人信服。马普所使用PROMETHEUS-VERTEX流体代码进行的新模拟在这一方面大大改进，考虑了所有3种中微子与能量相关的输运方程，还结合了较新的高密度物质物态方程。模拟考虑了9.6倍太阳质量与20倍太阳质量两种前身星，对后者还额外考虑了核子反弹、弱磁场、奇异夸克等复杂因素对中微子加热过程的影响，它们最多可以让中微子的不透明度发生10%左右的变化，最终成为爆发的关键。由于计算表明，无论在二维还是三维情况下，两种质量前身星都可以在激波反弹数百毫秒后成功引爆核心坍缩型超新星，且爆发能量与观测相符（只是三维情形对应的能量略高），它们第一次从根本上证明了中微子加热机制的可行性。
9.6倍太阳质量恒星的爆发过程，图中黑色实线表示各物质壳层随时间的演化，红色粗实线表示反弹激波波前的位置，红色区域表示星体内部物质强烈湍动的区域。可见激波在受阻120到130毫秒后再度复苏，最终在演化开始400毫秒后传播到6000千米外。（图片来源：Melson et al. 2015b）
与Ia型超新星相比，核心坍缩型超新星的模拟对计算资源的消耗也要大得多。当年FLASH中心在计算引力束缚暴震模型时，只在美国若干国家实验室的超级计算机上调用了不到800个CPU，总计花费数万个处理器机时。但是马普小组为了实现三维核心坍缩型超新星经由中微子机制的第一次自然引爆，在德国莱布尼兹计算中心的SuperMUC与西班牙巴塞罗那超级计算中心的地中海超算上动用了16000个CPU，前后共计消耗5亿个处理器机时，最终计算的只是0.5秒的演化过程。
第一次自然引爆三维核心坍缩型超新星期间使用的SuperMUC超级计算机（左，图片版权：LRZ 2012）与地中海（MareNostrum）超级计算机（右，图片版权：BSC 2013）。有趣的是，后者的机房是由教堂改造而来的。虽说乍看之下有些离谱，不过考虑教堂的防雷防静电性能一般都还不错，改作机房倒也恰当得很，否则如此结实的建筑简直是浪费。所以嘛，就好像教堂的最佳用途实为给超算中心做机房一样，必须把某些别有用心/不思进取的夏威夷土著心目中的所谓神山圣地改建成世界最大最先进的天文台才算是美观环保且高度尊重自然充分利用资源，什么原始文化都是扯淡，渎神最美科研最高！
接下来对于核心坍缩型超新星模拟来说，除了完善现有的物理过程（比如将这次忽略掉的奇异夸克对核子自旋的影响考虑在内）、进一步提高分辨率（这一点的前提当然是启用更强大的计算机，并调度更多的CPU）之外，更重要的是计算不同质量前身星的爆发过程，并寻找理论模型与观测量之间的联系。当然，如果好运来临，银河系或者邻近伴星系在不久之后发生超新星爆发，那就是验证理论再好不过的机会，因此随时做好观测的准备，让多波段多信使观测设备随时待命也是必需的。
末了是超新星模拟算法的一些小内幕，也是很多人认为数值模拟很“脏”的示例之一。本人并没有亲身参与过超新星的模拟，不过前阵在学习FLASH的过程中了解到了一些相关内容。FLASH也好，PROMETHEUS-VERTEX也罢，说是结合多种功能的流体模拟代码，其实如果单凭流体方程组的求解，是压根应付不了超新星爆发中心区域的情况的。因为这里的气体密度很高，相应的金斯长度很小；而为了避免出现气体团块的非物理分裂，每个网格的边长最长不能超过金斯长度的四分之一，这样一来，哪怕最密的自适应网格可能都难以应付中心区域的要求；另一方面，较高的密度意味着较短的自由下落时标，模拟的时间步长也需要随之缩短，最终可能会因步长过小而难以继续演化。PROMETHEUS-VERTEX对此问题的解决方案本人并不清楚，不过FLASH为此引入了具有特殊生灭规则的“下沉粒子（sink particle）”概念，用下沉粒子来描绘致密星核中心的连续介质演化。这些粒子具体的性质超出了本文的范畴，有兴趣者可以参见FLASH的使用手册的第19章；不过这样一来是不是又额外引入了数值算法方面的不确定性？本人不知，横竖近期也不会用到FLASH的这个功能嗯……
