N体シミュレーション

N体シミュレーション (エヌたいシミュレーション、N-body simulation) とは、天体物理学および天文学において、重力相互作用するN個の粒子の力学的な進化を数値的に計算するシミュレーションのことをいう。2体系つまりケプラー問題は可積分であるが、3体以上の系（重力多体系）は可積分ではなく、その力学的進化を定量的に予測するためには数値シミュレーションが必須である。太陽系や球状星団、銀河あるいは宇宙の大規模構造など、重力多体系は宇宙のあらゆる領域において重要な役割を果たすため、N体シミュレーションは宇宙に関する理論的研究において極めて重要な役割を果たしている。

$N$体シミュレーションを用いて得られた銀河団における質量分布。

ナイーブなN体シミュレーションの実装（直接総和法）は重力相互作用の計算に O(N²) のコストを要するため、より大規模かつ長時間にわたるシミュレーションを実現することは計算機科学特に高性能計算 (HPC) の分野においても興味深い問題であり、複数のゴードン・ベル賞がN体シミュレーションの研究に対して与えられた。現在でもN体シミュレーションはコンピュータのベンチマークのためにしばしば用いられる。

概要

三体問題の数値積分は最も単純な$N$体シミュレーションと言える。

太陽系における惑星の運動、原始惑星系円盤における微惑星の集積過程、球状星団の構造の力学的な進化、宇宙の大規模構造の形成といった問題は、系を構成する要素（惑星や恒星、暗黒物質など）が互いに重力相互作用を及ぼし合うものである。その進化を定量的に求めることは天文学および宇宙物理学の多くの分野で重要であるが、それは解析的な手法では困難であり、数値シミュレーションに頼らざるを得ない。$N$体シミュレーションはこのような系の数値シミュレーション手法のひとつであり、系を構成する要素を多数の（$N$個の）粒子とみなし、それらの粒子が互いにニュートン重力を及ぼし合うものとモデル化するものである。

典型的な$N$体シミュレーションでは、$i=1,2,\cdots ,N$ 番目の粒子（質量 $m_{i}$) の運動方程式は

$m_{i}{\frac {d^{2}\mathbf {r} _{i}}{dt^{2}}}=-\sum _{j\neq i}Gm_{i}m_{j}{\frac {\mathbf {r} _{i}-\mathbf {r} _{j}}{\left|\mathbf {r} _{i}-\mathbf {r} _{j}\right|^{3}}}$

であり、これを適切な初期条件のもとで数値的に積分することが主たる目標となる（宇宙の大規模構造など宇宙論的なスケールでの進化を扱う場合、後述のように方程式が修正されるが、本質的な差異はない）。その最も単純かつ非自明な問題が三体問題である。

シミュレーションの目的によって様々な数値計算上の困難が存在し、各分野毎に様々な手法が開発されてきた。

• 太陽系シミュレーションの場合、長時間の進化の結果として最終的にどのような状態になるのかに興味がある。この場合、数値誤差の累積が最大の問題であり、高次のシンプレクティック積分を用いるなどの誤差を極力抑えるための工夫が必要となる。また、この場合中心天体（太陽）の重力が支配的で、それに対して摂動が加わった系であるとみなせるため、天体力学的手法など特有の手法が用いられる。
• 球状星団の場合、重力相互作用のコストが重い上に、近接散乱や連星形成などの効果を正しく計算することが必要である。初期のGRAPEは球状星団の力学進化の計算に用いられ、大きな成果を上げた。
• 銀河またはより大きなスケールを扱う問題の場合、これは無衝突系であるため $N$ 体粒子は真の粒子である必要はなく、ある空間領域に存在する暗黒物質あるいはバリオンを束ねたものである。一方で、より計算領域を大きく、またシミュレーションの分解能を小さくすることが必要である。このためスーパーコンピュータなどの大規模並列計算機が採用され、また重力相互作用の計算コストを下げるためにアルゴリズム上の改善が続けられている。

無衝突系と衝突系

球状星団は典型的な重力多体系（衝突系）であり、その力学的進化の理論的研究において$N$体シミュレーションが重要な役割を果たした。

$N$ 体シミュレーションはその対象によって大きく無衝突系 (collisionless system) と衝突系 (collisional system) に分類される。これは二体緩和の効果が重要かどうかを意味し、それによってシミュレーションの性質が大きく変化する。

二体緩和とは、重力多体系において二体間の近接散乱による系の熱的な進化のことを言う。二体緩和による系の進化に要する時間スケール $t_{\mathrm {relax} }$ は、粒子数 $N$ と crossing-time $t_{\mathrm {cross} }$ を用いて

$t_{\mathrm {relax} }={\frac {N}{\ln N}}t_{\mathrm {cross} }$

と書ける。それ故に、極めて粒子数の大きなスケールでは二体緩和が効く時間は宇宙年齢よりも長くなり、その効果を無視することができる。例えば銀河は $N\sim 10^{10}$ 個の恒星からなる系であり、その力学的な進化には二体緩和の効果は重要ではないと考えられている。一方、球状星団では $N\sim 10^{4-6}$ であり、二体緩和が重要である。

二体緩和が効かない無衝突系では粒子数無限大の極限に相当する。このとき重力場はある種の平均場として扱うことが可能であり、個々の粒子に起因する特異性を考慮する必要はない。そのため、シミュレーションに際してツリー法などのより効率的なスキームを使用することができる。また、シミュレーションで扱われる粒子は真の粒子ではなく、位相空間のある領域を代表する点であると解釈される。

時間積分

積分子

$N$ 体シミュレーションはエネルギーが保存するため、時間積分子としてリープ・フロッグ法などのシンプレクティック数値積分法がしばしば採用される。例えば太陽系の高精度シミュレーション、宇宙論的構造形成など。

衝突系では後述する可変時間刻みと相性の良い予測子修正子法やエルミート積分子も用いられる。

可変時間刻みと独立時間刻み

自己重力系は一般に重力不安定性により密度ゆらぎが成長し、高密度領域を形成するように進化する。その結果、高密度領域の中心部では自由落下時間

$t_{\mathrm {free} }={\frac {1}{\sqrt {G\rho }}}$

が急速に短くなるため, 精度の良いシミュレーションを行うためには時間積分のタイムステップを動的に小さく調整することが望ましい。このような手法は可変時間刻みと呼ばれる。

ところが、特に衝突系で連星形成が起こるような状況では、一部の$N$体粒子は極めて短い時間で進化するものの、その他の大多数の粒子の軌道進化に小さな時間刻みが必要ないという可能性がある。この場合、最も小さな時間ステップに合わせて全体の時間刻みを調整すると、シミュレーションに多大な時間を要することになり、また不必要に小さな時間ステップに伴う数値積分誤差が累積する可能性がある。そこで必要な粒子のみ小さな時間刻み幅で時間積分を行う独立時間刻みという手法が開発された。この場合、その他の大多数の粒子については適当な補間を用いてその重力場を見積もり、必要な粒子のみ正確に時間積分を行うことになる。

重力相互作用

特に無衝突系においてはシミュレーションの規模 (つまり粒子数 $N$) を大きくすることが重要である。しかし粒子 $i$ に作用する重力

$\mathbf {F} _{i}=-\sum _{j\neq i}Gm_{i}m_{j}{\frac {\mathbf {r} _{i}-\mathbf {r} _{j}}{\left|\mathbf {r} _{i}-\mathbf {r} _{j}\right|^{3}}}$

をすべての粒子について愚直に計算する（これを直接総和法 direct summation という）ことは ${\mathcal {O}}(N^{2})$ という大きな計算時間を要するため、粒子数 $N$ を大きくすると急速に計算時間が増大し、現実的な時間で計算を終えることができなくなる。このため、PM法とツリー法という重力計算の精度を下げてでもより効率的な相互作用の計算アルゴリズムが開発された。現在ではこれらの方法を組み合わせた P³M 法やtree-PM 法が大規模シミュレーションにおいて標準的な方法として採用されている。

Particle-Mesh 法

Particle-Mesh (PM) 法は高速フーリエ変換に基づいて重力ポテンシャルの計算を行う。まず最初に計算領域にグリッドを生成し、各頂点での密度の値をその近傍の粒子分布に基づいて決定する。重力ポテンシャルは (フーリエ空間では) ラプラス方程式

$-k^{2}\Phi _{\mathbf {k} }=4\pi G\rho _{\mathrm {k} }$

により密度場に結びついているため、密度場のフーリエ係数 $\rho _{\mathbf {k} }$ を求め、そこから逆フーリエ変換することにより重力ポテンシャル $\Phi$ や重力場 $-\nabla \Phi$ を求めることができる。この計算時間はグリッド数を $M$ とすると ${\mathcal {O}}(M\ln M)$ である。

なお近くの粒子からの重力は直接法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを Particle-Particle Particle-Mesh (P³M) 法と呼ぶ。

ツリー法

ツリー構築のアニメーション。

Barnes & Hut (1986) により提案された、粒子分布をツリー構造という形で保持することにより重力相互作用の計算を ${\mathcal {O}}(N\ln N)$ で行うアルゴリズムは現在 Barnes & Hut のツリー法として広く用いられている。これは計算領域を表す立方体を階層的により小さな立方体に分割し、最終的に各立方体がひとつ以下の粒子しか含まないようにすることにより、粒子分布の情報をツリーとして保持するものである。ツリーの深さは ${\mathcal {O}}(\ln N)$ であるため、ツリーの構成に要する計算量は ${\mathcal {O}}(N\ln N)$ である。

ある粒子に作用する重力を計算する際には、遠方の粒子群からの寄与をまとめて計算することによりコストを削減する。この際に各立方体の重心および質量を用いるが、計算の精度を上げるために四重極モーメントなどを用いる場合もある。

なお近くの粒子からの重力はツリー法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを tree-PM 法と呼ぶ。

重力の近距離発散

ふたつの$N$体粒子間の距離が極めて小さくなると、両者の間に働く重力は任意に大きくなり得る。これは衝突系においては物理的に重要であり、その影響を正しくシミュレーションする必要がある。一方、無衝突系ではこの効果は物理的ではなく、カットオフにより除去される。

正則化

衝突系において近距離発散に対処するために可変時間刻み幅を用いる場合、時間ステップが際限なく小さくなり、シミュレーションがほとんど進まなくなってしまう。しかしながら、二体問題における近距離重力に起因する特異性は見かけのものであり、適切な変数変換により除去することができる。この手続きは正則化 (regularization) として知られる。

Burdet-Heggie正則化は、時間座標 $t$ を近接粒子の距離に応じて調整することで特異性を除去するもので、新しい時間座標 $\tau$ を二体の粒子間距離 $r$ と真の時間 $t$ から

$d\tau ={\frac {dt}{r}}$

により定義するものである。このとき二体の相対位置ベクトル $\mathbf {r}$ の従う方程式は

${\frac {d^{2}\mathbf {r} }{d\tau ^{2}}}-2E_{2}\mathbf {r} =-\mathbf {e} +r^{2}\mathbf {g}$

へと帰着される。ここに $E_{2}$ は二体の相対運動エネルギー、$\mathbf {e}$ は離心率ベクトル、$\mathbf {g}$ は他の天体による重力場である。この表式からわかるように、BH正則化により $r\to 0$ での特異性が除去される。

1965年に提案されたクスターンハイモ・シュティーフェル変換 (Kustaanheimo-Stiefel transformation) は3次元直交座標 $x,y,z$ を4次元のスピノルへと変換するもので、BH正則化よりも精度の良い結果が得られる。

$u_{1}^{2}={\frac {x+r}{2}}\cos ^{2}\psi ,\ \ u_{2}={\frac {yu_{1}+zu_{4}}{x+r}},\ u_{4}^{2}={\frac {x+r}{2}}\sin ^{2}\psi ,\ u_{3}={\frac {zu_{1}-yu_{4}}{x+r}}$

近距離カットオフ

無衝突系においては、$N$体粒子は真の粒子ではなく、多数の粒子が占める位相空間上の領域を表す。そのため、$N$体粒子間に働く重力が近距離で発散する効果は物理的ではなく、適当なカットオフにより除去される必要がある。最も簡単なカットオフとしては重力ポテンシャル $\Phi$ を

$\Phi (r)=-{\frac {GM}{\sqrt {r^{2}+\epsilon ^{2}}}}$

へと変更するものである。ここに $\epsilon$ は距離の次元を持つ定数で、この距離スケールより近距離での重力の発散を抑える効果を持つ。このポテンシャルは$N$体粒子がプラマーモデルであるような質量分布を持つと仮定した場合に得られるものに等しい。$\epsilon$ の値は平均粒子間距離のオーダーに選ばれる。

宇宙論的$N$体シミュレーション

Gadgetを用いた宇宙論的構造形成シミュレーションのスナップショット。色は暗黒物質の密度分布を表す。フィラメント構造やダークマターハローが形成されていることが確認できる。

大規模構造の形成などの宇宙論的な問題は$N$体シミュレーションが用いられる典型的かつ重要なセットアップであるが、宇宙膨張を考慮する必要があるという点で他の問題とは異なっている。この種のシミュレーションは非線型成長後の物質の密度ゆらぎのパワースペクトル、あるいはダークマターハローの密度プロファイルや質量関数を求めるために用いられる。これらの量は観測可能量であり、実際の観測データと比較することにより例えば宇宙論パラメータの制限を与えることができる。

運動方程式

宇宙膨張は、宇宙の現在の大きさを1とする相対的な大きさを表すスケール因子 $a(t)$ により表され、その時間発展はフリードマン方程式

${\frac {1}{a}}{\frac {da}{dt}}=H_{0}{\sqrt {\Omega _{\mathrm {\Lambda } 0}+\Omega _{m0}a^{-3}+\Omega _{r0}a^{-4}}}$

により与えられる。ここに $H_{0}$ はハッブル定数、$\Omega _{x0}$ は密度パラメータである。これにより、逆に独立変数として時刻 $t$ の代わりにスケール因子 $a$ を用いることができる。

粒子の座標としては宇宙膨張の効果を取り除いた共動座標 $\mathbf {x}$ が用いられる。これは固有座標 $\mathbf {r}$ と

$\mathbf {r} =a(t)\mathbf {x}$

という関係にある。粒子の速度はその微分 ${\frac {d\mathbf {r} }{dt}}=H\mathbf {r} +a(t){\frac {d\mathbf {x} }{dt}}$ であるが、初期宇宙での発散を回避するために $\mathbf {w} :={\sqrt {a(t)}}{\frac {d\mathbf {x} }{dt}}$ が用いられる。最終的な運動方程式は

${\frac {d\mathbf {x} }{da}}={\frac {\mathbf {w} }{S(a)}}$

${\frac {d\mathbf {w} }{da}}=-{\frac {3}{2}}{\frac {\mathbf {w} }{a}}+{\frac {1}{a^{2}S(a)}}\nabla \Phi (\mathbf {x} )$

$S(a)=H_{0}{\sqrt {\Omega _{m0}+a^{3}\Omega _{\Lambda 0}}}$

である。

周期境界条件

無限に広い計算領域を実現することは不可能であるため、宇宙論的シミュレーションでは周期境界条件が採用される。通常、計算領域は立方体であり、その一片の長さを $L$ とするとき、座標 $x$ と $x\pm L$ の点は同一視される。

周期境界条件のもとでは隣接する立方体に自らの構造のコピーが存在するため、それが及ぼす重力を考慮する必要がある。これは分子動力学法においてクーロン電場に対して開発されたエバルトの方法を適用することで可能であり、付近のボックスからの重力は直接計算し、遠方のボックスからの重力をフーリエ級数の形で取り入れることにより効率的に精度よく計算される。例えば、座標原点にある質量 $m$ の質点がつくる（規格化された）重力ポテンシャルは、周期境界条件のもとで

${\hat {\Phi }}(t,{\mathbf {x} })=-{\frac {Gm}{a}}\sum _{{\mathbf {n} }\in \mathbf {Z} ^{3}}{\frac {{\mathrm {erfc} \,}(\alpha |{\mathbf {x} }-{\mathbf {n} }L|)}{|{\mathbf {x} }-{\mathbf {n} }L|}}-{\frac {Gm}{\pi aL}}\sum _{{\mathbf {h} }\neq {\mathbf {0} }}\exp \left(-{\frac {\pi ^{2}|{\mathbf {h} }|^{2}}{\alpha ^{2}L^{2}}}\right){\frac {1}{|{\mathbf {h} }|^{2}}}e^{i{\mathbf {k} }_{\mathbf {h} }\cdot {\mathbf {x} }}+{\frac {Gm}{a}}{\frac {\pi }{\alpha ^{2}L^{3}}}$

となる。ここに $\alpha$ は任意の正の定数、$\mathrm {erfc}$ は相補誤差関数である。

初期条件

Λ-CDMモデルでは宇宙論的ゆらぎはインフレーション期にガウス分布に従って生成されたと考えられており、線型摂動の範囲では密度ゆらぎのパワースペクトルが指定されれば初期条件を確率的に生成することができる。パワースペクトルは与えられた宇宙論パラメータのもとで宇宙論的摂動論に基づいて計算することができる。なお宇宙論的$N$体シミュレーションで最も広く用いられる初期条件は2次のラグランジュ摂動 (2LPT) に基づくものである。

プロジェクト

ソフトウェア

Volker Springelが中心になって開発されたen:GADGETは銀河や宇宙論的構造形成を主なターゲットとする無衝突系のシミュレーションコードであり、1998年にバージョン1が、2005年にバージョン2が公開された。スーパーコンピュータなどの大規模並列計算機で動かすために並列化されており、C言語によって実装されている。ライセンスはGNU GPL。

2010年代には牧野淳一郎らが中心となって汎用的な多体問題シミュレーションフレームワークであるFDPSが開発されている。

ハードウェア

GRAPEは杉本大一郎、牧野淳一郎らによって東京大学で開発された$N$体シミュレーション専用計算機であり、重力相互作用の計算をパイプラインとして物理的に実装することにより効率的に計算を進めるものである。現在も国立天文台などで運用されている。

また GPUの活用（GPGPU）も進められている。

シミュレーション

21世紀に入ってから、特に大規模なシミュレーションはスーパーコンピュータを長時間占有する必要があるため、大規模なシミュレーションを行いその出力を公開するプロジェクトが行われるようになった。その有名なものがミレニアム・シミュレーションであり、他に The Aquarius Project などがある。2010年代に実行されたen:Illustris projectは$N$体シミュレーションだけでなく、星形成やAGNといったバリオン物理を考慮した流体シミュレーションを行っている。

歴史

重力多体系の計算機を用いた研究すなわち$N$体シミュレーションは1960年代から実際的な研究で採用されるようになった。例えば1963年のen:Sverre Aarsethによる $N=100$ 体のシミュレーション、1964年の Hénon & Heiles による第三積分に関する数値的研究（これは $N=1$ であるが）、1973年の Hénon による多体系の安定性の研究など。ジェームズ・ピーブルスは1970年に $N=300$ 体を用いて銀河団形成過程のシミュレーションを行った。その後も1979年には Efstathiou & Jones が $N=500$ 体による銀河回転の研究など、計算機の発達に伴ってより大規模なシミュレーションがなされるようになっていった。

より大規模なシミュレーションの要求は強く、1986年に Barnes & Hut はツリー法を導入し、同時期に PM 法も確立した。1989年にはGRAPEプロジェクトがスタートしている。

一方で積分スキームに関する研究も進められた。天体力学の分野からは1990年に吉田春夫によりシンプレクティック積分子の一般的な構成方法が示された。その翌年牧野はエルミート積分子を導入した。やがて対称型公式の有用性が認められるようになった。

2005年のミレニアム・シミュレーションでは $N=1.0\times 10^{10}$ の宇宙論的構造形成シミュレーションが遂行されるに至った。

参考文献

• Binney, James; Tremaine, Scott (2008). Galactic Dynamics (Second ed.). Princeton University Press. ISBN 978-0-691-13027-9 
• Aarseth, Sverre J. (2010). Gravitational N-Body Simulations: Tools and Algorithms (Cambridge Monographs on Mathematical Physics). Cambridge University Press. ISBN 978-0521121538 
• Mo, Houjun; van den Bosch, Frank; White, Simon (2010). Galaxy Formation and Evolution. Cambridge University Press. ISBN 978-0-521-85793-2 
• 牧野淳一郎「重力多体系の数値計算(<シリーズ>物性研究者のための計算手法入門)」『物性研究』第76巻第3号、物性研究刊行会、2001年6月。 
• 牧野淳一郎, 福重俊幸, 小久保英一郎, 川井敦, 台坂博, 杉本大一郎 (2007年3月13日). “N体シミュレーション啓蟄の学校教科書” (PDF). 国立天文台. 2020年5月24日閲覧。

関連項目

• 計算物理学
• 多体問題
• 三体問題
• ピタゴラス三体問題

外部リンク

• 小松信義, 木綿隆弘, 木村繁男「重力エネルギーが支配的なN体系の速度緩和」『理論応用力学講演会 講演論文集』第61回理論応用力学講演会セッションID: GS06-01、日本学術会議 「機械工学委員会・土木工学・建築学委員会合同IUTAM分科会」、2012年、222-222頁、doi:10.11345/japannctam.61.0.222.0、NAID 130005020725。