天文学家中有这样一群人，他们既不观星，也不刷公式，而是通过模拟来研究浩瀚星空。在星团、星系这样的恒星系统中，往往包含着上百万颗恒星，规模如此惊人的恒星系统该如何处理？本期“赛先生天文”为你解读天文学中恒星系统数值模拟的开端与发展。

球状星团模拟的 RGB 图像（Credit: MPA）

　　撰文：李硕（中国科学院国家天文台）
　　编辑：韩越扬　　

　　提起天文学家，我们脑海里最先想到的可能就是一群夜猫子。他们昼伏夜出，没事喜欢守在望远镜前，不知道在捣鼓些啥，总之看上去很厉害就是了。或者，像有些科幻电影里的角色一样，手里抱本厚厚的《星系动力学》，办公室黑板上写满了连符号都看不懂的公式。当然后者看上去更拉风些。毕竟守着望远镜“发呆”的难度系数并不算高，可没事在黑板上刷公式就不是谁都玩得起的了。

　　其实，还有那么一小撮天文学家，他们觉得天天盯着望远镜和黑板太没创意，而有了计算机就可以当上帝，想要什么都可以用计算机来模拟。这群喜欢没事趴在显示器前，十指噼里啪啦敲个不停的家伙管自己叫：搞数值模拟的。

　　今天，我们就来看看这个行当是怎么发家的。篇幅有限，我们只聊最“简单”也是最早开始发展的数值模拟——恒星系统直接动力学模拟中的直接N体数值模拟。

　　预测恒星的运动

　　图 1. 著名的球状星团杜鹃座 47，含有数百万颗恒星（图源：哈勃拍摄，NASA and Ron Gilliland (Space Telescope Science Institute)；地面拍摄，David Malin, © Anglo-Australian Observatory）

　　我们的银河系是由数千亿颗恒星组成的。其中还有大量由众多恒星聚集而成的星团。既然这些系统都是由恒星构成的，而支配恒星运动的又是万有引力，那么我们是不是能够通过计算每颗恒星受到的引力来预测他们的运动呢？

　　想法挺好，但问题的关键是计算量。以最简单的三体系统为例，要想近似地在纸面上计算，需要把系统演化的过程分成若干间隔尽可能小的时间段（这就是计算机模拟中所谓的时间步长），然后在每个时间点计算每个天体受到其它天体的引力之和，再根据已有的位置速度“预测”下一时间每个天体的位置速度，循环反复。

　　聪明的你可能已经意识到了，这类计算的复杂度随着天体数目是以平方指数增加的。而一个恒星系统，不说像银河系那样有几千亿颗恒星，也不说像球状星团那样动辄几十万颗恒星（不信？可以试着在图 1 里面数数），就算只有十几颗，计算量也足以让人吐血了。所以说，没有计算机，这事也就是想想，真要去做那简直是疯了。

　　“疯狂”天文学家

　　图 2. 最早的星系交汇模拟与合并星系 NGC2207 照片的比较（图源：Erik Holmberg, 1941, The Astrophysical Journal, 94, 385; NGC2207: ESO）

　　现在有请“疯狂”天文学家 Erik Holmberg 登场。这位仁兄在没有计算机的年代完成了两个恒星系统交汇的模拟。图 2 左边是他模拟的恒星系统交汇过程中的形态演化，右图则是我们拍到的合并星系 NGC 2207。

　　这项工作可是在 1941 年发表的，还要再等五年，人类第一台计算机才会诞生。前面提过，没有计算机帮助这样的计算可是要出人命的。而 Holmberg 活到了 92 岁，显然并没在 33 岁的时候就累死在稿纸堆里。那他是怎么做到的呢？

　　Holmberg 找了 74 个灯泡来代表两个星系。对，你没看错，74 个灯泡！他给不同的灯泡通上了不同的电压来代表星系中恒星的密度分布。越靠近中心电压越高，灯泡也就越亮，他通过移动灯泡来模拟星系在交汇过程中的演化。可怎么计算每个灯泡下一步应该挪到哪里去呢？Holmberg 找来了光电管测量各处的光强。因为光与引力一样是遵循距离平方反比衰减的，他巧妙地利用光强的测量代替了引力的计算。

　　就这样，史上第一个天体系统动力学演化的模拟就算是完成了。当然，这个模型的分辨率很低——用 37 个灯泡代替整个星系是有点勉强的。那多加灯泡不行吗？醒醒吧，刚才我们提到的球状星团有多少恒星来着？几十万颗。几十万个灯泡，要开烧烤店么？这画面太美，实在不敢想。

　　图 3. Hoerner 用来完成首例星团数值模拟的晶体管计算机西门子 2002（图源：Wikipedia）

　　计算机诞生后，德国天文学家 Sebastian von Hoerner 在五十年代末开始了在晶体管计算机上的模拟。当时的计算机没有显示器和键盘（图3）。所有程序都要借助穿孔的纸带输入计算机，得到的输出也是一条打满了孔的纸带。那时候程序员们每天的工作就是扎眼，代码不按行来计算而是按捆算。Hoerner 在 1959 年就是用这样的计算机实现了 16 颗星的数值模拟。梦想总是要有的，万一实现了呢？Hoerner 同学估算了一下，发现按照当时的计算速度，“只要”花 1400 亿年就能模拟一个球状星团了。

　　“头号玩家”入局

　　几乎是同时，另一位重要的玩家参与到了这项游戏中，这便是挪威人 Sverre Aarseth。这位奇人开创了 Nbody 系列模拟程序、获得过国际象棋通讯赛国际大师称号、炒股挣钱买过法拉利、50 岁时迷上登山并因此冻掉过好几根脚趾、80 岁时玩野外漂流，为了骗取工作人员许可谎称自己才 60。最有个性的是，在剑桥大学工作了几十年却一直是普通研究员，相当于现在的博士后。要说学术贡献不够吧，有一颗小行星可是以他命名的，用来表彰其贡献的。他的人生经历简直可以写一本书了。呃，其实书也已经写了。

　　言归正传，Aarseth 在A. Schlüter 的建议下为每个天体赋予了独有的计算步长。这就避免了两个密近天体因为步长太小而拖慢整个计算的情况。这一方法在 1963 年提出，于 1986 年才最终完善为等级阻塞时间步长方案。从 80 年代初开始，Aarseth 用 FORTRAN 语言编写了著名的 Nbody 程序，并在之后与合作者一道维护更新至今。

　　为了更真实地模拟恒星系统，他们又把目标对准了双星。银河系中双星非常普遍。而由于其轨道过于靠近，需要非常小的时间步长才能保证计算精度。这就导致包含双星的计算会非常缓慢。借助数学的帮助，Aarseth 等人在六七十年代通过一种叫做 KS（Kustaanheimo & Stiefel）正则化的坐标变换，将双星的轨道积分转换成了一种简单的谐振子计算。后来 Seppo Mikkola 与 Aarseth 又在 1992 年将这一工具应用在了多个天体上，即所谓的 KS 链。

　　同时，为了进一步提高精度，日本天文学家 Junichiro Makino 与 Aarseth 在 91-92 年引入了厄米积分方法，通过泰勒展开的数学变换巧妙地用一阶导数取代了模拟中所需的二三阶导数计算，从而在不增加计算复杂性的前提下有效地提高了计算精度。

　　随着计算机硬件的飞速发展，1996 年，德国天文学家 Rainer Spurzem 与 Aarseth 一起实现了 10000 个天体的数值模拟。而在同一年， Makino 借助专门设计的计算设备 GRAPE-4 实现了 32768 个天体的模拟。

　　GRAPE 是一种专门用来进行引力计算的特殊硬件，很容易进行并行化计算。在 GPU 异构加速计算兴起之前，GRAPE 是最高效的引力系统模拟硬件。但是，GRAPE 专精于引力计算，很难被大规模应用于其它领域，因此制造成本高昂。而 GPU 在这方面拥有与 GRAPE 类似的能力，且成本更低。在 GPU 计算兴起后，GRAPE 也就慢慢地退出了历史舞台。

　　图 4. 天龙数值模拟使用的超级计算机 Hydra（图源：马克思普朗克计算与数据中心）

　　悬赏一瓶威士忌

　　九十年后期，Aarseth 的 Nbody 模拟程序已经发展为 Nbody6 并加入了各种真实的物理过程。随着高性能计算机群的发展，人们开始考虑将模拟工作交给多个 CPU 并行完成。1999 年，Spurzem 在 Nbody6 的基础上开发出了第一个并行版本 Nbody6++。几年之后，一些简单的动力学模拟（包括 GRAPE）已经可以处理百万天体量级的计算。

　　当然这些计算还无法兼顾恒星演化等星团中重要的物理过程。为此，英国数学与天文学家 Douglas C. Heggie 专门悬赏，将为第一个完成真实球状星团模拟的工作赠送一瓶他珍藏的威士忌。

　　2000 年后，GPU 加速计算迅速兴起。Spurzem 尝试借助德国基金会帮助建造 GPU 机群未果。他得到的答复是：“我们已经有了非常强大的超级计算机，为什么还要 GPU 机群呢？”而中国此时为了鼓励 GPU 计算正准备资助建造若干个 GPU 超级计算机。就这样 Spurzem 来到了中国，开始将 Nbody++ 与 GPU 加速相结合。

　　与此同时，Aarseth 与日本天文学家 Keigo Nitadori 合作，为 Nbody6 加入了 GPU 计算支持，使模拟工作可以在一台有 GPU 加速卡的工作站上完成。Nbody 数值模拟由此进入了 GPU 加速时代。

　　图 5. 天龙数值模拟星团“快照”与球状星团 M13 照片（图源：模拟结果，Wang, L.， Monthly Notices of the Royal Astronomical Society， 2016, 458, 1450；M13 照片，NASA, ESA, and STScI/AURA）

　　终于，在 2015-2016 年，经过与 Spurzem、Aarseth 等人多年的合作，来自北京大学的博士生王龙完成了 Nbody 程序的并行 GPU 版本——Nbody6++GPU，并在先进的通用 GPU 异构机群上（图4，可以和图 3 比较一下）一举实现了百万恒星量级的球状星团模拟——天龙星团数值模拟，毫无悬念地赢得了 Heggie 的威士忌。

　　图 5 中可以看到经过可视化处理的数值模拟结果。左图是模拟结果按星团中不同种类天体呈现的“观测”效果，中间的圆图是所有天体叠加的结果，与右图中 Hubble 望远镜拍摄的球状星团 M13 非常接近。而与观测不同的是，数值模拟可以轻松地区分各种天体，甚至连观测中基本无法发现的黑洞都能轻松标识。左图中最右边的白色背景小图就是恒星质量级黑洞在星团中的分布。

　　故事讲到这里该告一段落了。但还要多啰嗦两句，我们只介绍了直接N体方法数值模拟的成长，实际上随着这一领域的发展，出现了很多方法用以研究恒星系统的演化。即便是N体模拟，在九十年代以后也出现了大量的替代程序，比如 Starlab、phi-GPU、Hi-GPU 等。篇幅所限就不一一介绍了。

　　图 6. 银河系中心恒星轨道，所有恒星都围绕着在中心不可见的超大质量黑洞（图源：UCLA Galactic Center Group – W.M. Keck Observatory Laser Team）

　　最后想说的是，Nbody 的故事还远未到落幕的时刻。我们虽然已经能够实现星团的模拟，但仍然有很多物理过程无法很好地实现。比如说星系中心的超大质量黑洞（图6）或球状星团中心可能存在的中等质量黑洞，现有的N-body 程序还无法很好地处理这些极大质量天体与普通天体的相互作用。要想模拟一个有着几千亿颗恒星的星系，我们恐怕还有相当长的路要走，好在我们总能找到几个喜欢天天坐在电脑前面发呆的家伙来干这事。

　　致谢：Rainer Spurzem 教授为本文的撰写提供了大量资料，在此深表谢意！

　　作者简介

　　李硕，中国科学院国家天文台助理研究员。2012 年于北京大学物理学院天文学系获得理学博士学位。研究领域是引力系统演化，集中在超大质量黑洞与星系的共同动力学演化。

　　　　参考资料

　　[1] Holmberg, Erik“On the Clustering Tendencies among the Nebulae. II. a Study of Encounters Between Laboratory Models of Stellar Systems by a New Integration Procedure”， The Astrophysical Journal, 1941, 94, 385

　　[2] von Hoerner, S.， “Die numerische Integration des n-Körper-Problemes für Sternhaufen. I”， Zeitschrift fur Astrophysik, 1960, 50 

　　[3] Aarseth, S. J.， “Dynamical evolution of clusters of galaxies, I”， Monthly Notices of the Royal Astronomical Society， 1963, 126, 223

　　[4] Aarseth, S. J. & Zare, K.， “A regularization of the three-body problem”， Celestial Mechanics, 1974, 10, 185

　　[5] Makino, J. & Aarseth, S. J.， “On a Hermite integrator with Ahmad-Cohen scheme for gravitational many-body problems”， Publications of the Astronomical Society of Japan, 1992, 44, 141

　　[6] Mikkola, S. & Aarseth, S. J.， “An implementation of N-body chain regularization”， Celestial Mechanics and Dynamical Astronomy， 1993, 57, 439

　　[7] Spurzem, R. & Aarseth, S. J.， “Direct collisional simulation of 10000 particles past core collapse”， Monthly Notices of the Royal Astronomical Society， 1996, 282, 19

　　[8] Makino, J.， “Postcollapse Evolution of Globular Clusters”， The Astrophysical Journal, 1996, 471, 796

　　[9] von Hoerner, S.， “How it all started”， Dynamics of Star Clusters and the Milky Way， ASP Conference Series, Vol. 228. Edited by S. Deiters, B. Fuchs, R. Spurzem, A. Just, and R. Wielen. San Francisco: Astronomical Society of the Pacific. ISBN: 1-58381-060-9, 2001.， p.11

　　[10] Aarseth, S. J.， “Gravitational N-Body Simulations”， pp. 430. ISBN 0521432723. Cambridge, UK: Cambridge University Press, November 2003. (New Edition 2010)

　　[11] Wang, L.， Spurzem, R.， Aarseth, S.， Giersz， M.， Askar, A.， Berczik, P.， Naab, T.， Schadow， R.， Kouwenhoven, M. B. N.， “The DRAGON simulations: globular cluster evolution with a million stars”， Monthly Notices of the Royal Astronomical Society， 2016, 458, 1450