天文学中的数值模拟：GRAPE相关
因为种种原因拖来拖去，数值模拟的第3篇介绍文章一直拖到现在才写成。此番介绍的是N体引力作用相关的过程，确切地说是利用GRAPE完成的N体模拟，其他算法的相关结果还是留待日后介绍。至于GRAPE是什么，可参见本站之前的简介文章。简单说来，它是专门设计用于求解引力作用的计算机硬件设备，速度了得。本文主要的内容也就是GRAPE官方网站上列出的几项较重要的工作。
在正文之前，再多说两句相关模拟的思路。在许多情况下，天体间的相互作用可以近似为自引力粒子的作用，如太阳系系统、星团、星系乃至整个宇宙。由于引力是长程力，难以被忽略，因此在参与粒子数目较多时，会带来引力计算量以粒子数平方的速度增加，而其他部分只是以粒子数增加。GRAPE是提高速度的一条途径，它不是从算法（如粒子网格和Tree算法，将远场粒子做近似）上做改进来节约时间，而是通过专门设计的硬件来有效提高速度。GRAPE本身专门用于引力计算，其他部分是主机的任务。
再来看看GRAPE的外观，这一台是东京大学的GRAPE-4，后文介绍的工作中，至少有一部分是用它来完成的：
顺便提一句所谓的Plummer软化因子，这在之前的硬件介绍中也有所提及，具体表达式请参见前文。引进该长度量纲的因子的目的是为了避免近场奇点的存在。若形象地说，则是在引力表达式中考虑该因子后，可以把粒子看做是一种分布而非纯几何点。据一位主修宇宙学的学姐介绍，取法可以是采取网格长度的若干分之一，不过因为本人没有过实际操作经验，所以也就不敢乱讲了。
月球的起源
很不爽的一点是，这项工作对应的文献（Kokubo et al. 2000）没有找到全文，本人又实在是懒得去图书馆翻Icarus的旧刊，所以就姑且粗粗一写了。
其实相关工作还有很多，Kokubo et al. 2000在摘要里就提到了Ida et al. 1997的一些情况，说是他们的结果和Ida相符得很好。这项结果说穿了就是高精度N体模拟，基础是时下流行的月球起源撞击说。由于撞击，地球周围产生了碎片盘，接下来就要计算由此盘产生的月球的性质。
先是初始盘面，中央为地球：
然后在自引力的作用下，粒子碰撞，盘面收缩。洛希极限之外，物质团块逐渐积聚起来。由于开普勒剪切效应，逐渐产生旋臂状的结构：
旋臂的力矩将物质传输到洛希极限外，旋臂末端则坍缩成为较大的星子，并逐渐扫清周围物质：
最终结果是，行星盘的演化多半会导致洛希极限外近圆轨道上单一大型卫星的产生，而这正是现在月球的情况。一旦月球尺度足够大，在地月引力作用下，洛希极限内的物质残余会落回地球。具体有多少物质转化为卫星则取决于盘初始的角动量，不过与盘的初始质量关系并不大。
这一过程的动画可以在这里下载。实际上这次模拟描述的正是当下月球起源的主流理论，最早由William K. Hartmann与Donald R. Davis发表在1975年的Icarus上。撞击物可能的来源是行星过程的产物。不过Kokubo et al. 2000的摘要里有一句乍看起来很吓人的话：The formation time scale of the moon is of the order of a month。本人对月球起源的理论细节不很熟悉，难道真的是只需要个把月的时间，整个吸积过程就可以完成了？这一点还望指教……
关于星团
主要内容来自Zwart et al. 1999。此项模拟的蓝本是星团R136，也就是NGC 2070，位于大麦哲伦云中。说来这个星团有点邪门，成员星相当年轻，但彼此间隔极密，预计有演化成球状星团的趋势。正因为其中成员星演化与动力学演化的耦合，才让它成为研究的焦点。
致密星团R136（图片提供： N. Walborn et al., WFPC2, HST, NASA）
模拟的一个关键问题是恒星彼此的碰撞，由此会导致一系列后果，如X射线双星、毫秒脉冲星、蓝离散星、超大质量恒星等，此外还考虑了恒星的演化。在本人接触的一点模拟皮毛中，基本还是在“collisionless”的条件下打转转，也知道碰撞情况不是那么好处理的。为解决这一问题，Zwart等人使用了名为starlab的软件包，它专门为致密星团的模拟而设计，并备有GRAPE接口以及恒星演化模块SeBa，不过在实际使用时还要考虑大质量恒星质量损失的问题。至于恒星演化在模拟上的具体实现方法，还是等到以后有机会再介绍吧。
简要说来，Zwart et al. 19994认为，恒星碰撞并合后的产物是一颗新的恒星，且新恒星的等效年龄比并合前有所减少。而依照先前的模拟结果，并合过程中的质量损失可以忽略不计。实际模拟参与的恒星数为12000余颗或是6000余颗（实际上比R 136的恒星数少几倍），前者变换不同的质量损失率，后者的系统弛豫时标和密度分布也有所不同。
模拟得出的结果涉及恒星的并合与星团结构的演化，不过初始条件并未考虑双星。并合产生的大质量恒星质量可达太阳的180余倍，且由于过程中的refresh，这些恒星的寿命比单星来得长。与理论预期不同的是，模拟中，碰撞并合的次数要多上10倍左右，而且碰撞过程的主导是大质量而非中小质量恒星。这些差异的可能来源是先前对双星形成以及大质量恒星在动力学摩擦的作用下向星团中心转移过程的忽略以及。由于并合发生的次数基本上随恒星数目成正比，因此又有了可将结论用于更大规模星团的推断。至于星团的结构，核心区域是随着时间扩展的，原因在于双星加热以及星风损失等。
星团模拟一瞥，对应恒星数目为12288。碰撞产生的恒星以白色表示，壳层则是超新星爆发的产物。（Zwart et al. 1999）
考虑R 136的实际情况，其中质量最大的几颗蓝离散星可能正是过往碰撞史的写照。它们有着较小的年龄以及极丰富的氢元素含量。
冷暗物质晕问题
要说这篇讲述冷暗物质晕的Fukushige & Makino 1997的工作背景，是在1995年。哈勃太空望远镜的高分辨率观测表明，椭圆星系核心的密度并非如先前预期一般均匀，而是呈幂律分布，越靠近中心的地方密度越高，亮星系中是，一般星系更可达到
。按照标准的冷暗物质理论，这是难以解释的。
研究冷暗物质晕演化所经常采用的图景是所谓的“hierarchical clustering”，是由下到上由小到大的过程。主要问题在于，不论是什么时候，哪次研究，结论也都是建立在模拟的基础之上的。先前的介绍文章中已经讲过，数值模拟的一大不确定因素就是运算过程中引入的非物理解，这也算是本人对数值模拟的第一印象。先前的模拟中，除了均匀的中心密度分布，固然也出现过形式的分布，却难以找到确切的理论解释，而这些工作较低的质量分辨率以及相对大的软化也让结果屡遭质疑。
排除模拟假象的一条途径是增加精度。Fukushige & Makino 1997的粒子数目（70余万）是前辈工作的10倍以上，又考虑了较弱的软化，加之双精度下GRAPE运算的有效性，这一目的是达到了。单个质量400万倍于太阳的粒子分布在直径2兆秒差距的球状区域内，密度峰的选择则是在8兆秒差距范围内对密度丰度使用COBE数据作归一，并进行高斯平滑。之后，就是考虑结合哈勃流之后的N体演化。
此次所得的密度分布比略平，不过比
来得更陡。这可以归因于冷暗物质晕的非等温分布。而通过改变软化势的大小也说明，先前认为的平坦分布似乎确实应归罪于较大的软化势带来的假象。
模拟一瞥，依次为红移等于8.7、5.1和1.8时的对应情况。（图片来源：Fukushige & Makino 1997）
这项工作还得到了速度与温度的分布情况。不过由于较亮的星系中心密度分布较为平坦，并不能单纯用冷暗物质来解释。文章是将此解释为星系中央特大质量黑洞的作用。
冷暗物质晕的温度结构图（图片来源：Fukushige & Makino 1997）
关于星系晕的N体模拟思路还可以进一步推广到星系团演化上，只是需要扩大模拟的尺度，如这段动画。
其实在GRAPE之外，更多的引力相关N体模拟都是在算法上动脑筋，如Tree算法或是粒子网格一类，这些东西容我后文讲解吧，最近忙就一个字……
写作本文时顺便发现了一个有意思的东西，所谓的ACS (Art of Computational Science) Project。这是由GRAPE的两位主要开发者Piet Hut与Jun Makino发起的，目前主要包括开源知识（Open Knowledge）与Maya两大计划，前者关注代码开发过程的知识共享，后者致力于致密恒星系统的软件开发。这东西值得研究一下……
