天文学中的数值模拟：PIC，兼说SPH
这篇文章应该算是烂尾中的烂尾了：从2009年春天开始动笔，中间写写停停一直拖了两年半还有余。最近将其赶完的原因倒也很实际，因为本人有心做一些相关的工作，赶紧先突击一下基本概念再说。
数值模拟往往需要求解偏微分方程组，传统的求解方法是为解域划分网格，可以用随流体运动的拉格朗日网格，也可以用固定于空间之上的欧拉网格。这次要说的Particle-in-Cell（PIC）方法则结合了以上二者，在表示流场的欧拉网格之上叠加以表示物质的拉格朗日粒子；而光滑粒子（SPH）方法更为干脆，索性屏弃了网格，直接以带有分布的粒子代替。
PIC的概念可以追溯到1955年，甚至早于第一款Fortran语言编译器的面世，后来John M. Dawson、Oscar Buneman、Charles K. Birdsall等人在20世纪60年代将其发扬光大。由于可以直接给出粒子本身的运动情况，PIC在使用中可以求解比网格差分更为微观的问题，是等离子体研究的重要工具。对于天体物理背景来说，粒子加速以及非热辐射都可以顺便给出，但是宏观大尺度现象并不适合用这一方法求解。作为对比，网格法是从描述流体（或磁流体）基本方程组入手，适合求解宏观现象。
很自然地，这种方法所应用的粒子疏密是物质密度的反映，粒子数量的选择不宜过少，且为了抑制非物理波动，一般以不太均匀的分布为佳。这里的粒子并不是真实的粒子，而是许多真实粒子的代表，至于各个网格的物理量则是其中所有模拟粒子相关参数的总和。使用模拟粒子也是有理论依据的，因为决定带电粒子在电磁场中运动的洛伦兹力效果只取决于电荷与质量之比，只要这一数字保持不变，模拟得到的结论也就是一样的。
计算过程分欧拉步与质点步，输运量根据质点运动给出。欧拉步还算好说，跟普通的网格没有本质区别。之前学习计算天文的时候，老师给出的例子是使用FTCS方法对动量差分，用交错法对能量差分，但是这个例子只针对普通流体，并没有涉及电磁场。如果考虑带有电磁场的等离子体，由于需要求解麦克斯韦场方程，需要换用其他方法，比如傅里叶分析一类，不过本人准备研究的那套代码暂时还没有拿到手，不知道具体操作方式。
质点步要求算每个模拟粒子的运动。首先根据欧拉步的结果，将求出的网格中央点速度加权平均，换算成粒子速度，然后根据运动方程（在等离子体的环境下就是洛伦兹方程）求出它们的新坐标，并用新的参数求出下一时刻的网格物理量，再重复上述步骤。为了抑制数值计算引入的非物理波动，一般粒子并非点粒子，而是带有分布。
可能是与模拟尺度受限有关（PIC算法要求格点尺度至少要与粒子的回旋运动半径同级，而不能过小，因此大尺度模拟所需的格点数量太多，很容易会超过计算机的能力。），PIC的天文应用并不多（多半PIC工作与实验等离子体有关），本人导师的合作者、阿拉巴马大学亨茨维尔分校的Ken-Ichi Nishikawa算是其中的先驱者之一，现在专事此业的研究小组哪怕在全球范围内也是屈指可数。下面就以Nishikawa的工作为例，展示一些PIC模拟结果。这些工作的基础程序是考虑狭义相对论效应的三维PIC并行代码，估计今后自己就要去啃它了。先是相对论性激波面附近的粒子密度演化，其对应的天文背景与喷流的传播有关，图中红、绿、蓝线分别表示密度跳变、接触间断面以及激波波前的运动情况：
同一背景下等离子体的Weibel不稳定性带来的电流密度分布：
有了粒子数密度与电磁场，从最基本的运动方程出发求解辐射是很自然的事情，图中深浅红线分别表示喷流电子在视线方向与法线夹0度与5度角时的辐射谱，深浅蓝线分别表示背景电子视线方向与法线夹0度与5度角时的辐射谱，由于不稳定性产生的磁场强度不是很大，这里的谱形更类似于轫致辐射：
当然如果磁场足够强，那么就可以得到类似同步辐射的光谱了，不过下图严格来说算不得是PIC模拟，只是在背景场中引入了两个检验粒子而已，不过使用的程序与前述工作相同，所以也在此一提。图中不同颜色的曲线则表示不同的视角：
图片来源：Nishikawa et al. 2009
虽说PIC的应用主要是针对小尺度现象的，不过也有例外，比如Nishikawa与合作者进行的相对论性PIC模拟，给出了黑洞吸积盘喷流的形成过程，速度不同的粒子可以沿喷流向黑洞两极喷出。不过这只是个toy model，黑洞尺度比天体物理中真正应付的要小很多，但不论如何这一工作尚属独家：
顺便一提的是，传统的（磁）流体模拟也可以引入检验粒子，那么它与PIC的区别在哪里？依个人理解，这主要在于算法是否自洽。PIC方法本身就可以得出粒子的运动情况，而检验粒子不能单凭网格给出，同时检验粒子也不参与流场的演化。
说完PIC，再提一提SPH，其实自己对后者的概念接触得还要多一些，不过它跟本人的课题倒是关系不大。前面说过，SPH方法是比PIC更进一步完全舍弃了网格，密度只以质点分布来表示。它起源于20世纪70年代末，最初就是为了天文应用而开发。由于取消了网格，质点又能自然地追踪流体本身的运动，SPH可以更容易地计算密度跳变或者离体的情况，实际经验也表明它在能量、动量、角动量等方面的守恒性也不错。
类似于PIC，SPH的粒子也是带有分布的（光滑粒子嘛）。模拟前首先要选择适当的分布做核，粒子密度分布形态即为粒子质量与核函数之积，而场量则是所有粒子的累加。至于将流体方程，大可参考《光滑粒子流体动力学》一书，本文不赘述。
SPH这一方法其实跟本人当下的课题关系不大，原因是其用途至今多半限于流体，对带磁场的等离子体的应用仍处起步阶段，据说是由于磁场散度为0这个特性很难在该算法下很好地保持的缘故。据本人所知，SPH的主要应用包括宇宙大尺度结构的形成与演化、星系并合过程、恒星形成等，比如本系列文章第一篇中的星系并合模拟应该就是这一算法的产物。至于具体工作本文就不细说了，留给后文。
BTW，Nishikawa先生的造访也真是不迟不早，于是Miramar航展+中途岛号航母参观计划就此泡汤……
