N体シミュレーション
Volker Springel(ドイツ語版)が中心になって開発されたen:GADGETは銀河や宇宙論的構造形成を主なターゲットとする無衝突系のシミュレーションコードであり、1998年にバージョン1が、2005年にバージョン2が公開された。スーパーコンピュータなどの大規模並列計算機で動かすために並列化されており、C言語によって実装されている。ライセンスはGNU GPL。
2010年代には牧野淳一郎らが中心となって汎用的な多体問題シミュレーションフレームワークであるFDPSが開発されている。
無衝突系においては、 N {\displaystyle N} 体粒子は真の粒子ではなく、多数の粒子が占める位相空間上の領域を表す。そのため、 N {\displaystyle N} 体粒子間に働く重力が近距離で発散する効果は物理的ではなく、適当なカットオフにより除去される必要がある。最も簡単なカットオフとしては重力ポテンシャル Φ {\displaystyle \Phi } を
Φ ( r ) = − G M r 2 + ε 2 {\displaystyle \Phi (r)=-{\frac {GM}{\sqrt {r^{2}+\epsilon ^{2}}}}}
へと変更するものである。ここに ε {\displaystyle \epsilon } は距離の次元を持つ定数で、この距離スケールより近距離での重力の発散を抑える効果を持つ。このポテンシャルは N {\displaystyle N} 体粒子がプラマーモデルであるような質量分布を持つと仮定した場合に得られるものに等しい。 ε {\displaystyle \epsilon } の値は平均粒子間距離のオーダーに選ばれる。
ここに E 2 {\displaystyle E_{2}} は二体の相対運動エネルギー、 e {\displaystyle \mathbf {e} } は離心率ベクトル、 g {\displaystyle \mathbf {g} } は他の天体による重力場である。この表式からわかるように、BH正則化により r → 0 {\displaystyle r\to 0} での特異性が除去される。
1965年に提案されたクスターンハイモ・シュティーフェル変換 (Kustaanheimo-Stiefel transformation) は3次元直交座標 x , y , z {\displaystyle x,y,z} を4次元のスピノルへと変換するもので、BH正則化よりも精度の良い結果が得られる。
u 1 2 = x + r 2 cos 2  ψ , u 2 = y u 1 + z u 4 x + r , u 4 2 = x + r 2 sin 2  ψ , u 3 = z u 1 − y u 4 x + r {\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}}}
ふたつの N {\displaystyle N} 体粒子間の距離が極めて小さくなると、両者の間に働く重力は任意に大きくなり得る。これは衝突系においては物理的に重要であり、その影響を正しくシミュレーションする必要がある。一方、無衝突系ではこの効果は物理的ではなく、カットオフにより除去される。
粒子の速度はその微分 d r d t = H r + a ( t ) d x d t {\displaystyle {\frac {d\mathbf {r} }{dt}}=H\mathbf {r} +a(t){\frac {d\mathbf {x} }{dt}}} であるが、初期宇宙での発散を回避するために w := a ( t ) d x d t {\displaystyle \mathbf {w} :={\sqrt {a(t)}}{\frac {d\mathbf {x} }{dt}}} が用いられる。最終的な運動方程式は
d x d a = w S ( a ) {\displaystyle {\frac {d\mathbf {x} }{da}}={\frac {\mathbf {w} }{S(a)}}}
d w d a = − 3 2 w a + 1 a 2 S ( a ) ∇ Φ ( x ) {\displaystyle {\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 Ω m 0 + a 3 Ω Λ 0 {\displaystyle S(a)=H_{0}{\sqrt {\Omega _{m0}+a^{3}\Omega _{\Lambda 0}}}}
である。
特に無衝突系においてはシミュレーションの規模 (つまり粒子数 N {\displaystyle N} ) を大きくすることが重要である。しかし粒子 i {\displaystyle i} に作用する重力
F i = − ∑ j ≠ i G m i m j r i − r j | r i − r j | 3 {\displaystyle \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 という)ことは O ( N 2 ) {\displaystyle {\mathcal {O}}(N^{2})} という大きな計算時間を要するため、粒子数 N {\displaystyle N} を大きくすると急速に計算時間が増大し、現実的な時間で計算を終えることができなくなる。このため、PM法とツリー法という重力計算の精度を下げてでもより効率的な相互作用の計算アルゴリズムが開発された。現在ではこれらの方法を組み合わせた PM 法やtree-PM 法が大規模シミュレーションにおいて標準的な方法として採用されている。