天文学中的数值模拟：编外篇之五·碰撞诱发磁重联过程（上）
先前本人研究过的一项课题是在解释伽玛暴瞬时辐射的ICMART(内部碰撞诱发磁重联与湍流）模型框架下探讨光变曲线的形态。ICMART模型是本人在国外的导师为了解决伽玛暴瞬时辐射经典内激波理论所面临的诸多问题而提出的，这些问题包括难以解释光球辐射的缺失、辐射效率过低、与统计关系相矛盾等等。该模型的大致思想是，在瞬时辐射期间，由于爆发中心能源抛出的一系列抛射物壳层速度不同，彼此会发生碰撞。但由于整个辐射区是磁场主导的，碰撞导致的内激波辐射很大程度上是被压制的。不过碰撞过程会导致磁力线扭曲，诱发雪崩式的级联磁重联，在此过程中粒子被众多小型重联区加速，最终产生了所见的伽玛光子。根据该模型的预言，伽玛暴瞬时辐射光变曲线分为快慢两种成分，前者对应单独的重联过程，而后者与壳层的碰撞有关。这一预言在很大程度上已经得到了理论和观测数据的支持，因此ICMART模型看上去颇有前景。
ICMART模型给出的伽玛暴辐射产生过程。（图片来源：Zhang & Yan 2011）
但是到现在为止，ICMART模型最大的不确定性就在于碰撞诱发磁重联过程的细节。这一过程高度复杂且非线性，无法用解析方法来描述，只能进行粗略的量级估计。至于这些估计可靠与否，甚至连碰撞诱发的雪崩式重联事件会不会发生，都还是没有人能说得清楚的。更加麻烦的是，为了确定粒子加速过程的细节，粒子模拟是必需的，具体说来是在相对论性磁流体/等离子体理论框架下进行三维PIC或SPH模拟，而有能力进行相关工作的小组，在全球范围内都是屈指可数的。
去年在美国洛斯阿拉莫斯国家实验室工作的李晖老师造访紫金山天文台，本人在此期间得知他的小组正在进行磁重联过程的PIC与磁流体模拟，而且与自己国外导师的一名学生合作，试图探讨伽玛暴中心能源区域可能发生的现象。近日相关论文已经在网络上公开（详见Deng et al. 2015以及Guo et al. 2015），于是就在此简要介绍一下他们的结果。由于相关内容较多，本文分上下两篇写成，上篇主要关注Deng et al. (2015)，重点是相对论性碰撞诱发磁重联，进而引起能量耗散的过程。
Deng et al. (2015)所用的程序是以洛斯阿拉莫斯实验室在10多年前开发的LA-COMPASS为基础改造的，其核心是通过分段抛物线方法（PPM）来求解流体的黎曼问题，并利用所谓的输运限制（constrained transport，实在不知道应该如何翻译）来确保磁场散度为零的苛刻条件。PPM法可以追溯到著名的Colella & Woodward (1984)，FLASH的前身PROMETHEUS也是与此系出同门。对于带有间断面的问题求解来说，这应该说是一种效果较好的算法，只消一两个格点就可以精确判明激波所在。不过对于刚刚啃完FLASH使用手册又继续抱着E. F. Toro的黎曼解器厚砖大啃特啃的本人来说，最近真是看到这几个字母就头疼不已啊……
与ICMART模型提出的壳层碰撞不同，Deng et al. (2015)考虑的是两个（最后又加到了多个）球状磁化团块彼此之间的碰撞，团块的尺度是1012到1013厘米。这一尺度远小于根据光深要求反推出的伽玛暴源区半径（1016厘米左右），而比当年本人根据观测到的光变时标选取的重联区尺度（3×1011厘米甚至更小）略小。考虑我们所见的伽玛暴辐射是集中在张角很窄的喷流中的，喷流截面的半径倒是应该远小于源区半径而远大于单个重联区，所以这样的尺度选取也还算合理。
在整个模拟区域内，假设初始物质密度分布均匀，且存在一个轴对称的大尺度磁场，其轴向分量保持通量为零，而环向磁场是基于黑洞吸积盘发电机理论构筑的。两个团块之间的相对位置用z方向的团块中心距离zd以及x方向的距离xs，而二者y轴坐标保持一致。由于团块只在z轴方向存在相对初始速度，因此前一参量决定了碰撞过程的延迟时间，而后一参量描述了两团块在碰撞方向上的彼此偏离程度。需要注意的一点是，这项工作主要选择了5/3而非极端相对论性的4/3作为气体的绝热指数。在大多数情况下这种处理问题不大，因为模拟是在团块的质心静止系中进行的，抛射物整体运动被无视，而重联过程只是中等相对论性的，速度还不会太极端。但在最猛烈的重联区，抛射物有可能达到极端相对论的速度，因此可能会出现某些偏差。
除了重联过程本身，模拟还特别关注了耗散效率。ICMART模型的一大卖点就是辐射效率超高。辐射效率在这里可以近似理解为初始系统能量最终耗散为抛射物内能的比例，由于是磁场主导的环境，最终可以近似为事件终止时的磁化因子σend加1的倒数，因此如果σend接近1，那么辐射效率可以轻易达到50%（内激波模型的效率一般只有百分之几），足以解释绝大多数的伽玛暴观测。但是实际的效率也需要通过数值模拟来确定。Deng et al. (2015)并没有计算具体的粒子加速过程以及辐射机制，这里的耗散效率是由非碰撞与碰撞过程耗散掉的电磁能量之差再比上非碰撞过程耗散的能量来求得的。
用于测试的双磁化团块系统结构图，左上：磁力线；右上：磁化系数分布；左下：x方向磁场；右下：y方向磁场。（图片来源：Deng et al. 2015）
最初测试用的系统构型如上图所示，两个磁化团块彼此有着相同的磁化因子（取为8）以及确定的相对速度和位置，然后发生碰撞。出于计算耗散效率的考虑，首先作者进行了非碰撞模拟，明确这一过程耗散掉的电磁能量总数。由于初始磁场并不是处在完全的平衡态的，两个团块会膨胀，并在磁场的驰豫过程中将部分电磁能转化为动能和热能，最终磁压和气压之间达到平衡，系统进入准静态状态。磁场的驰豫过程发生得很快，在此期间电磁能量迅速下降，随后几乎保持不变。
含有碰撞的过程就比较有趣了。首先通过数值模拟，有一个问题得以明确：磁化团块的碰撞是非弹性的，它们会黏结在一起。与伽玛暴早期余辉的光变曲线形态相似，磁化团块碰撞期间的电磁能量演化也可以分为陡降期、平台期、常规衰减期，最终是准静态阶段。当然，这二者之间是否存在物理关联还不敢妄言，毕竟是完全不同的过程，而且磁化团块的碰撞要解决的也是伽玛暴瞬时辐射而非余辉问题。从下图可见，考虑碰撞之后，电磁能量最初的下降要比非碰撞过程略缓慢一些，但由于碰撞本身也会消耗掉相当一部分能量，而且碰撞过程会压缩磁力线，因此最终剩余的电磁能更少。
电磁能量在两个磁化团块碰撞（实线）与非碰撞（虚线）事件中的演化，下方的曲线是耗散效率的演化。（图片来源：Deng et al. 2015）
与非碰撞过程一样，碰撞过程的陡降阶段也对应磁场的驰豫调整，随后团块的碰撞才会真正发生。这一过程最终可以达到的耗散效率是35%左右，比ICMART模型最初的估计（σend取为1，效率为50%）略低，但也远高于内激波模型。其实看过论文之后，本人也并不是很清楚判断终止模拟的条件，现在的理解是达到准静态阶段后，系统的参数不再存在明显变化，因此可以认为演化已经完成，具体操作还有待日后向这项工作的参与者请教了。这里非常有说服力的一点是，模拟过程中最终达到的σend略大于1，类比ICMART模型的理论框架推算出的辐射效率对应为35.7%，非常接近模拟结果，这反过来证明了提出ICMART模型的Zhang & Yan (2011)给出的量级估计是较为可信的。
看一看流场中各种物理量的分布图，可能有助于更好地理解这几个阶段：
系统演化的4个阶段中各物理量的分布，最左侧两列为电流，中列为速度，最右一列为密度。（图片来源：Deng et al. 2015）
图中可见在前两个阶段，两个团块接触面上是存在明显的电流片的，这是磁重联存在的明证。由于初始两个团块具有方向相反的x方向磁场（y方向磁场的差异与xs有关），在碰撞过程中二者相互接触并发生挤压，反向磁力线被扭曲，由此发生重联也是在情理之中。值得注意的是在第2个阶段（平台期）的x轴方向速度场，电流片似乎导致了大量三维指向不同的速度分量，这说明大尺度的多重重联事件有望诱发众多子喷流——这正是ICMART模型所期待的。
不过平台期为何存在？大致说来，可以认为是碰撞事件不仅诱发了磁重联，还导致了团块中物质和磁场的强压缩，后者可以增加电磁能量，所以在一定程度上总电磁能量大致不变。随后的常规衰减阶段就是前述过程进一步的延续了，只是压缩效应为之减弱，因此无法平衡磁重联，总电磁能量下降。在这一阶段，电流强度较先前有所衰减，外流物速度也从平台期的局域超阿尔芬与快磁声波速度下降到阿尔芬速度和快磁声波速度或者更低的水平。最后耗散逐渐停止，系统进入准静态阶段，能量演化曲线也趋于变为与非碰撞过程平行的形态。
不同阶段系统速度场的分布情况，图中对阿尔芬速度进行了归一化。（图片来源：Deng et al. 2015）
为了尽量减少数值过程离散化带来的误差，Deng et al. (2015)还额外分析了不同尺寸的格点（但保持模拟区域尺度不变）对结果产生的影响。总的来说，模拟结果与分辨率关系不大，不过随着分辨率的上升，耗散效率稍有增加，且最终电磁能量的演化会聚性更好。更加重要的趋势是，分辨率越高，最终外流物可以达到的最大速度也就越大。由于较高的分辨率意味着较小的数值阻抗与较窄的重联区，这样的结果都还在情理之中，而且也增强了主要结论的可靠性。如果进一步提高分辨率，并将背景环境设置为全面磁化的喷流，那么还是有可能形成大量相对论性子喷流的。
随后文中又讨论了不同的初始参数对模拟结果的影响。简单说来，团块的相对运动速度越大，其耗散效率越高；背景压强越强，效率也越高；团块在碰撞方向上的偏离量xs与效率关联比较复杂，大致是在xs取3时效率最高，更大的偏离量会减小碰撞接触面积，因此效率会下降，但较小的偏离量不能保证y方向磁场相反，因此也会减弱磁重联的强度。不过不同的zd以及背景密度与最终的效率关联不大。另外模拟过程中气体绝热指数无论是取非相对论性的5/3还是极端相对论性的4/3，求出的结果都非常相似，因此虽然严格说来最好能够逐格点根据速度和温度计算特定时间步骤上的局域绝热指数，但考虑对结果改善不大而且还会大大增加计算量，这样做似乎不是非常有必要。
不同绝热指数给出的电磁能量和磁化因子演化，图中实线表示碰撞过程，虚线表示，红色表示非相对论性绝热指数，绿色表示极端相对论性绝热指数。（图片来源：Deng et al. 2015）
文章的高潮是多团块碰撞模拟。前面提到过，多团块碰撞有望产生ICMART模型所要求的子喷流。但真正意义上的多团块碰撞模拟是要考虑大量位置、尺度、速度以及出现时间随机分布的团块的，费时费力，因此这项工作只是试探性地计算了四团块的碰撞。这一过程可以达到更高的耗散效率（约40%），而且起初的碰撞形成的纽结结构在后期会进一步碰撞，形成新的重联区，且原始重联区也会随时间推移发生进一步的演化，最终整个系统结构极其复杂，正佐证了ICMART模型的设想——雪崩式重联。
四团块碰撞过程的电流（左）与速度场（右）演化过程。（图片来源：Deng et al. 2015）
因此虽然Deng et al. (2015)只是模拟了最简单的双团块和四团块碰撞诱发磁重联的过程，其结果却为ICMART模型带来了好消息，其主要假设一一得到了证实。此外这一结果还可以应用到解释活动星系核TeV快速光变以及蟹状星云高能耀发现象上去。为了进一步细化碰撞诱发磁重联的过程，未来还需要进行更复杂的多团块碰撞模拟。
不过且慢，这里还有另一个有待解决的重要问题：碰撞诱发的磁重联过程如何加速粒子？这就是后文要涉及的内容了，我们下篇见。
