天文学中的数值模拟：Woosley与大质量坍缩星
终于决定了，简介过后，数值模拟问题就按一个一个小话题写下去，可能是介绍某一研究对象，也可能介绍某人或是某团队的工作。写作时间可能会拖得很久，毕竟这是一个艰深的领域，别说我一个小字辈了，就算是大牛也很难对其了解得面面俱到，因此倘或其中有错误，还请大家不吝赐教。此番先从Stanford Woosley的坍缩星相关工作说起，一来与本人的专业有关，二来他也算是我最早听说的数值模拟高人之一。
这位Woosley究竟是何方神圣呢？加州大学Santa Cruz分校教授兼天文系系主任，HETE-2的合作者，曾任美国物理学会天体物理分会执行委员会成员与马克斯·普朗克研究所的洪堡教授，International Journal of Modern Physics E编辑，2005年美国物理学会Hans A. Bethe奖得主，获奖理由是“在恒星演化、元素合成、核心坍缩型与Ia形超新星理论与伽玛暴解释等方面广泛而重要的贡献，特别是伽玛暴的坍缩星模型”（for his significant and wide-ranging contributions in the areas of stellar evolution, element synthesis, the theory of core collapse and type Ia supernovae, and the interpretation of gamma-ray bursts–most notably, the collapsar model of gamma-ray bursts）……不管怎么说，如今若谈到坍缩星或是伽玛暴的模拟研究，Woosley的工作是非说不可的。
在Woosley的主页上，赫然列着三百来篇他参与撰写的论文，还不包括会议文集，乍看能把来者吓得一楞一楞的。合作者也是牛人层出不穷，一开头就是著名的D. Clayton（这学期本人屡次被该人的砖头教科书折磨），后面更是有多篇重量级文章，让人感叹不已啊。
当然我不可能先通读这300余篇文章，再来写本文。所以就挑自己了解的一些内容零散讲一讲好了，主要是坍缩星相关的一些问题。首先是本人对Woosley的第一印象，还是在中学时代看到的动画。当时并不了解这张图片背后的真义，更不清楚伽玛暴究竟为何物，还以为霍金的黑洞蒸发理论早就把一切问题都解决了呢。只是现在已经知道，霍金的观点似乎从来就没有成为过伽玛暴的主流，就算是当年爆发距离悬而未决的时候，坍缩星与中子星起源才算得上是成了大气候的模型。
上图是Woosley的某项坍缩星模拟结果，时间跨度27秒。与此有关的大量计算核心都是二维流体模拟。既然是以数值模拟的名义进行介绍，模拟思路和方法也会是结论以外的另一个介绍重点。
大质量恒星的坍缩几十年以来都是令理论家头大的问题，甚至有什么“上千个理由表明核心坍缩型超新星不会爆发，但它还是爆发了”之类的说法。这个内容以后还会再专门撰文介绍。不过有一点值得注意，也就是超大质量恒星。倘或恒星质量过大的话，坍缩铁核之外的氧与硅包层就会有较大的密度和厚度，黑洞在核心坍缩阶段即可形成，因此在许多计算中往往无法得到预期的超新星爆发，人称所谓的“失败的超新星”（failed supernova，作为对比，正统的核心坍缩型超新星理论中，引发星体爆炸的是坍缩停止，产生中子星之后产生的反弹激波）。如果不考虑自转，事情倒也好说，星体本身直接以流体动力学时标坍缩成黑洞就是了。但真实恒星的外包层往往有着很大的角动量，黑洞形成后，吸积盘关联问题和随后的一系列过程那是相当的棘手。
Woosley正是研究上述“失败的超新星”的先行者之一，他在此方面的工作始于1983年与Bodenheimer的合作。发表的论文中探讨了一颗具有自转的25倍太阳质量恒星的坍缩过程。模拟的基础是流体力学基本方程，包括连续性方程以及能量动量守恒方程，并未考虑磁场因素。计算始于内区的8倍太阳质量，其中高速坍缩的核心质量1.5倍于太阳。为方便起见，系统采用的是对称构型，实际只求解了二维平面的一个象限。考虑核反应对温度的敏感性，径向网格的划分主要照顾了包层的内区，也就是温度较高的地方。包层外围（2至8倍太阳质量）的成分用纯氧近似，内区考虑氧参与的核反应。至于核心，初始状态干脆将中子星视为质点，与包层衔接处也只用边界条件带过；随着角动量与物质的积累，换作多方指数为3的多方球模型，相当于相对论性简并气体。
坍缩星演化后期密度（左）与速度场（右）的分布情况，可见存在氧燃烧产物的高速外流。（图片来源：Bodenheimer & Woosley 1983）
这项模拟的一大关键是角动量及其分布的选取。Bodenheimer与Woosley先是采用O型星的典型数值作为标准状态，之后为检验对参数的敏感性，又取了不同的角动量与其他参数（如压力梯度与引力之比）计算数次。由此得出的一大重要结果是，只要包层内区角动量与压力合适，超新星爆发仍会发生，光变曲线类似II型超新星。
Bodenheimer & Woosley 1983只是试探性的计算，首当其冲的毛病就是没有把考虑坍缩核心与包层同时考虑在内。接下来要讲的是Woosley 1993，也是一篇相当经典的文献，引文次数多达700余次，讨论的是源于“失败的超新星”的伽玛暴，爆发的成因是具有自转的沃尔夫—拉叶坍缩星中心产生的黑洞喷射出的火球。之所以说是失败，是因为这类超新星与一般低质量的超新星不同，由于中心形成的不是中子星而是黑洞，外层下落的物质无法反弹产生激波而使星体炸飞。这次除了角动量，恼人问题的还有粘滞产生的热量。文章本身并没有讲到太多的数值模拟，只是提出了伽玛暴成因的关键结论，还讨论了吸积盘的结构。当年的伽玛暴模型种类过百，正处筛选时期。随后Woosley的理论逐渐成为了焦点，虽然现在看来还是存在一些细节问题，比如辐射的方向性。
必然要关注的还有Woosley与他的高足Andrew MacFadyen、Weiqun Zhang等人在1990年代末的工作，正是MacFadyen给这些工作带来了漂亮的可视化。期间的不少计算涉及到一个叫做PROMETHEUS的二维流体力学代码，似与超新星爆发模拟有关，可以追溯到Fryxell, Mueller & Arnett 1989。这个PROMETHEUS也是个不得了的程序，被后人引用无数，只是原始文献的全文却是不大好找。可以确定的一点是，前面的简介中提过的FLASH这个牛物，就是以PROMETHEUS作为流体力学计算模块基础的。
PROMETHEUS应用的一个例子：超新星1987A的湍流模拟。湍流使超新星前身星的分层洋葱结构被搅乱，大大增加了过程的复杂性。（图片来源：Mueller, Fryxell & Arnett 1991）
就拿MacFadyen & Woosley 1999作为例子。这是一篇大红大紫的文章，引文数量也是高达600多。主要内容是对源自“失败超新星”的伽玛暴的详细模拟，前身星星体质量在太阳的10倍以上，氢包层业已损失殆尽。这里的情况其实和Bodenheimer & Woosley 1983有些类似，角动量再次成为关键因素：角动量过低，物质径直落入黑洞，无法产生爆发；过高的话，吸积盘形成的地方距离黑洞太远，中微子损能不明显，吸积过程受抑制。只有角动量适中的旋转氦星才能形成致密吸积盘，并由中微子过程或是喷流将能量辐射出去，从而形成爆发。为弄清参数敏感性，也同样是选取了多种组合。
由此得出的重要推论是，伽玛暴更容易在低金属丰度的恒星中形成，原因是这类恒星半径较小，质量损失也没有高金属丰度的大质量恒星那般严重，因此可以保持所需的角动量。所以爆发很可能随着红移有所演化，这也是当下方兴未艾的研究课题。
这项工作对PROMETHEUS的主要修改是物态方程方面，包含了正负电子组成的费米气体，对于高密区还进行了库仑修正，此外还用了又一个PPM代码。初始条件，也就是演化到终点的前身星，是由一个叫做KEPLER的程序计算而得的，这个程序也是Woosley早年参与的成果之一。下图是KEPLER的计算实例。
20世纪70年代利用KEPLER求出的不同质量的恒星演化晚期电子丰度与核心坍缩速度关系图。（图片来源：Weaver, Zimmerman & Woosley 1978）
模拟的过程持续时间为20秒，分为3大阶段，总计200万个时间步长。第一阶段持续约两秒，主要过程是吸积盘的形成。长度10余秒的第二阶段是吸积盘向黑洞传送物质的准静态过程，传送速率与吸积盘从外包层获取物质的速率相近，每秒不到太阳质量的0.1倍，比自由下落的速率约莫低一个数量级。期间的中微子与磁场作用有可能让能量在黑洞两极积累，为后面的爆发奠定基础。最后一个阶段是星体爆发，积累的能量（1051尔格的量级）沿喷流喷出，将恒星沿轴向瓦解，并在前行数百天文单位的距离之后引发伽玛暴。但模拟并没有将这一阶段完全覆盖。
坍缩核心周围的吸积盘在准稳定吸积状态的截面图，在此模型中并无喷流形成。左为中心区域，右为全局图景。（图片来源：MacFadyen & Woosley 1999）
大致的图景如此，当然其中还有大量细节问题，比如粘性的影响，中微子光度问题，黑洞的克尔参数演化，极区能量与动量的积累。最终模拟结果是相对论性强喷流的产生，张角大约可以覆盖1%的天区。喷流瓦解了吸积盘，并于穿越恒星表面之后演变为高度相对论性。依初始粘性、吸积率、角动量与视线方向等参数的不同，所得的爆发光度相差很大，正与据信为大质量恒星坍缩起源的长暴观测事实吻合。
喷流产生0.824秒后的能流分布图典型，洛伦兹因子是10左右，张角10度。此时喷流已经向外运动了7000千米。（图片来源：MacFadyen and Woosley 1999）
其实还需要再做做喷流的文章。喷流的产生是个持续的不稳定过程，正符合Rees & Mészáros 1994对的伽玛暴内激波推断。倘或喷流物质速度不均匀，分别具有洛伦兹因子与
、相隔
抛出的物质在前行
的距离后就会发生碰撞产生内激波，并释放大量能量，更详细的结果还要讨论散射光学厚薄的转换。喷流进动之前也为Woosley考虑过，不过是看得一头雾水，算了，先不说了。
后来的MacFadyen, Woosley & Heger 2001还讨论了另一种机制。之所以说MacFadyen & Woosley 1999的计算是以“失败超新星”为基础，那是因为黑洞与吸积盘形成在先，外行激波出现在后；这回干脆考虑形成于正统超新星爆发过程中的喷流，并无失败之说。这样的理论可以解释持续时间长达几百秒的伽玛暴，或是光度极高的超新星，同时还能说明某些超新星遗迹的不对称。
第二种机制产生的喷流结构。左上：轴向速度；右上：径向速度；左下：压力分布；右下：密度分布。左为全局图，右为中心区域放大图。（图片来源：MacFadyen, Woosley & Heger 2001）
有没有未决的问题呢？其中一个涉及流体方程。Woosley等人早期的计算基本只涉及了非相对论性的流体，后来才考虑了相对论因素；而且从恒星坍缩直至爆发形成还要用到半解析的结论，并未做到用数值方法再现全过程；关于伽玛暴前身星的金属丰度也存在不少争议，甚至有说爆发应形成于双星系统中，并非遵循Woosley提出的图景……不过一切也不能抹杀上述工作的意义，而针对问题的改进也在进行之中。
话说那些核合成相关的文献，本人是暂时无力全部理解，所以本文就先不涉及了。最后找几张MacFadyen的Scheides Cluster Simulations相关图片，不过只知道上排两图与坍缩星吸积盘和喷流有关，下排究竟是什么东西至今没搞明白……
