N体シミュレーション

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${\displaystyle 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}}}$

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

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

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

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

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

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

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

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

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

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

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

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

である[20]。

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

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

${\displaystyle {\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}}}}$

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

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

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

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

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

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

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