本文轉載自 賽先生
原創： 李碩
天文學家中有這樣一群人，他們既不觀星，也不刷公式，而是透過模擬來研究浩瀚星空。在星團、星系這樣的恆星系統中，往往包含著上百萬顆恆星，規模如此驚人的恆星系統該如何處理？本期「賽先生天文」為你解讀天文學中恆星系統數值模擬的開端與發展。
球狀星團模擬的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年於北京大學物理學院天文學系獲得理學博士學位。研究領域是引力系統演化，集中在超大質量黑洞與星系的共同動力學演化。
在此十分感謝公眾號賽先生為我們提供眾多的優質文章！在下只做好文章的搬運工，希望好的東西可以被更多的人看到
