N体シミュレーション

N体シミュレーション (エヌたいシミュレーション、N-body simulation) とは、天体物理学および天文学において、重力相互作用するN個の粒子の力学的な進化を数値的に計算するシミュレーションのことをいう。2体系つまりケプラー問題は可積分であるが、3体以上の系（重力多体系）は可積分ではなく、その力学的進化を定量的に予測するためには数値シミュレーションが必須である。太陽系や球状星団、銀河あるいは宇宙の大規模構造など、重力多体系は宇宙のあらゆる領域において重要な役割を果たすため、N体シミュレーションは宇宙に関する理論的研究において極めて重要な役割を果たしている。 ナイーブなN体シミュレーションの実装（直接総和法）は重力相互作用の計算に O(N2) のコストを要するため、より大規模かつ長時間にわたるシミュレーションを実現することは計算機科学特に高性能計算 (HPC) の分野においても興味深い問題であり、複数のゴードン・ベル賞がN体シミュレーションの研究に対して与えられた。現在でもN体シミュレーションはコンピュータのベンチマークのためにしばしば用いられる。 太陽系における惑星の運動、原始惑星系円盤における微惑星の集積過程、球状星団の構造の力学的な進化、宇宙の大規模構造の形成といった問題は、系を構成する要素（惑星や恒星、暗黒物質など）が互いに重力相互作用を及ぼし合うものである。その進化を定量的に求めることは天文学および宇宙物理学の多くの分野で重要であるが、それは解析的な手法では困難であり、数値シミュレーションに頼らざるを得ない。 N {\displaystyle N} 体シミュレーションはこのような系の数値シミュレーション手法のひとつであり、系を構成する要素を多数の（ N {\displaystyle N} 個の）粒子とみなし、それらの粒子が互いにニュートン重力を及ぼし合うものとモデル化するものである。 典型的な N {\displaystyle N} 体シミュレーションでは、 i = 1 , 2 , ⋯ , N {\displaystyle i=1,2,\cdots ,N} 番目の粒子（質量 m i {\displaystyle m_{i}} ) の運動方程式は m i d 2 r i d t 2 = − ∑ j ≠ i G m i m j r i − r j | r i − r j | 3 {\displaystyle 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 {\displaystyle N} 体粒子は真の粒子である必要はなく、ある空間領域に存在する暗黒物質あるいはバリオンを束ねたものである。一方で、より計算領域を大きく、またシミュレーションの分解能を小さくすることが必要である。このためスーパーコンピュータなどの大規模並列計算機が採用され、また重力相互作用の計算コストを下げるためにアルゴリズム上の改善が続けられている。 N {\displaystyle N} 体シミュレーションはその対象によって大きく無衝突系 (collisionless system) と衝突系 (collisional system) に分類される。これは二体緩和の効果が重要かどうかを意味し、それによってシミュレーションの性質が大きく変化する。 二体緩和とは、重力多体系において二体間の近接散乱による系の熱的な進化のことを言う。二体緩和による系の進化に要する時間スケール t r e l a x {\displaystyle t_{\mathrm {relax} }} は、粒子数 N {\displaystyle N} と crossing-time t c r o s s {\displaystyle t_{\mathrm {cross} }} を用いて t r e l a x = N ln ⁡ N t c r o s s {\displaystyle t_{\mathrm {relax} }={\frac {N}{\ln N}}t_{\mathrm {cross} }} と書ける。それ故に、極めて粒子数の大きなスケールでは二体緩和が効く時間は宇宙年齢よりも長くなり、その効果を無視することができる。例えば銀河は N ∼ 10 10 {\displaystyle N\sim 10^{10}} 個の恒星からなる系であり、その力学的な進化には二体緩和の効果は重要ではないと考えられている。一方、球状星団では N ∼ 10 4 − 6 {\displaystyle N\sim 10^{4-6}} であり、二体緩和が重要である。 二体緩和が効かない無衝突系では粒子数無限大の極限に相当する。このとき重力場はある種の平均場として扱うことが可能であり、個々の粒子に起因する特異性を考慮する必要はない。そのため、シミュレーションに際してツリー法などのより効率的なスキームを使用することができる。また、シミュレーションで扱われる粒子は真の粒子ではなく、位相空間のある領域を代表する点であると解釈される。 N {\displaystyle N} 体シミュレーションはエネルギーが保存するため、時間積分子としてリープ・フロッグ法などのシンプレクティック数値積分法がしばしば採用される。例えば太陽系の高精度シミュレーション、宇宙論的構造形成など。 衝突系では後述する可変時間刻みと相性の良い予測子修正子法やエルミート積分子も用いられる。 自己重力系は一般に重力不安定性により密度ゆらぎが成長し、高密度領域を形成するように進化する。その結果、高密度領域の中心部では自由落下時間（英語版） t f r e e = 1 G ρ {\displaystyle t_{\mathrm {free} }={\frac {1}{\sqrt {G\rho }}}} が急速に短くなるため, 精度の良いシミュレーションを行うためには時間積分のタイムステップを動的に小さく調整することが望ましい。このような手法は可変時間刻み（英語版）と呼ばれる。 ところが、特に衝突系で連星形成が起こるような状況では、一部の N {\displaystyle N} 体粒子は極めて短い時間で進化するものの、その他の大多数の粒子の軌道進化に小さな時間刻みが必要ないという可能性がある。この場合、最も小さな時間ステップに合わせて全体の時間刻みを調整すると、シミュレーションに多大な時間を要することになり、また不必要に小さな時間ステップに伴う数値積分誤差が累積する可能性がある。そこで必要な粒子のみ小さな時間刻み幅で時間積分を行う独立時間刻みという手法が開発された。この場合、その他の大多数の粒子については適当な補間を用いてその重力場を見積もり、必要な粒子のみ正確に時間積分を行うことになる。 特に無衝突系においてはシミュレーションの規模 (つまり粒子数 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法とツリー法という重力計算の精度を下げてでもより効率的な相互作用の計算アルゴリズムが開発された。現在ではこれらの方法を組み合わせた P3M 法やtree-PM 法が大規模シミュレーションにおいて標準的な方法として採用されている。 Particle-Mesh (PM) 法は高速フーリエ変換に基づいて重力ポテンシャルの計算を行う。まず最初に計算領域にグリッドを生成し、各頂点での密度の値をその近傍の粒子分布に基づいて決定する。重力ポテンシャルは (フーリエ空間では) ラプラス方程式 − k 2 Φ k = 4 π G ρ k {\displaystyle -k^{2}\Phi _{\mathbf {k} }=4\pi G\rho _{\mathrm {k} }} により密度場に結びついているため、密度場のフーリエ係数 ρ k {\displaystyle \rho _{\mathbf {k} }} を求め、そこから逆フーリエ変換することにより重力ポテンシャル Φ {\displaystyle \Phi } や重力場 − ∇ Φ {\displaystyle -\nabla \Phi } を求めることができる。この計算時間はグリッド数を M {\displaystyle M} とすると O ( M ln ⁡ M ) {\displaystyle {\mathcal {O}}(M\ln M)} である。 なお近くの粒子からの重力は直接法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを Particle-Particle Particle-Mesh (P3M) 法と呼ぶ。 Barnes & Hut (1986) により提案された、粒子分布をツリー構造という形で保持することにより重力相互作用の計算を O ( N ln ⁡ N ) {\displaystyle {\mathcal {O}}(N\ln N)} で行うアルゴリズムは現在 Barnes & Hut のツリー法として広く用いられている。これは計算領域を表す立方体を階層的により小さな立方体に分割し、最終的に各立方体がひとつ以下の粒子しか含まないようにすることにより、粒子分布の情報をツリーとして保持するものである。ツリーの深さは O ( ln ⁡ N ) {\displaystyle {\mathcal {O}}(\ln N)} であるため、ツリーの構成に要する計算量は O ( N ln ⁡ N ) {\displaystyle {\mathcal {O}}(N\ln N)} である。 ある粒子に作用する重力を計算する際には、遠方の粒子群からの寄与をまとめて計算することによりコストを削減する。この際に各立方体の重心および質量を用いるが、計算の精度を上げるために四重極モーメントなどを用いる場合もある。 なお近くの粒子からの重力はツリー法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを tree-PM 法と呼ぶ。 ふたつの N {\displaystyle N} 体粒子間の距離が極めて小さくなると、両者の間に働く重力は任意に大きくなり得る。これは衝突系においては物理的に重要であり、その影響を正しくシミュレーションする必要がある。一方、無衝突系ではこの効果は物理的ではなく、カットオフにより除去される。 衝突系において近距離発散に対処するために可変時間刻み幅を用いる場合、時間ステップが際限なく小さくなり、シミュレーションがほとんど進まなくなってしまう。しかしながら、二体問題における近距離重力に起因する特異性は見かけのものであり、適切な変数変換により除去することができる。この手続きは正則化 (regularization) として知られる。 Burdet-Heggie正則化は、時間座標 t {\displaystyle t} を近接粒子の距離に応じて調整することで特異性を除去するもので、新しい時間座標 τ {\displaystyle \tau } を二体の粒子間距離 r {\displaystyle r} と真の時間 t {\displaystyle t} から d τ = d t r {\displaystyle d\tau ={\frac {dt}{r}}} により定義するものである。このとき二体の相対位置ベクトル r {\displaystyle \mathbf {r} } の従う方程式は d 2 r d τ 2 − 2 E 2 r = − e + r 2 g {\displaystyle {\frac {d^{2}\mathbf {r} }{d\tau ^{2}}}-2E_{2}\mathbf {r} =-\mathbf {e} +r^{2}\mathbf {g} } へと帰着される。ここに 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} 体粒子は真の粒子ではなく、多数の粒子が占める位相空間上の領域を表す。そのため、 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 } の値は平均粒子間距離のオーダーに選ばれる。 大規模構造の形成などの宇宙論的な問題は N {\displaystyle N} 体シミュレーションが用いられる典型的かつ重要なセットアップであるが、宇宙膨張を考慮する必要があるという点で他の問題とは異なっている。この種のシミュレーションは非線型成長後の物質の密度ゆらぎのパワースペクトル、あるいはダークマターハローの密度プロファイルや質量関数を求めるために用いられる。これらの量は観測可能量であり、実際の観測データと比較することにより例えば宇宙論パラメータの制限を与えることができる。 宇宙膨張は、宇宙の現在の大きさを1とする相対的な大きさを表すスケール因子 a ( t ) {\displaystyle a(t)} により表され、その時間発展はフリードマン方程式 1 a d a d t = H 0 Ω Λ 0 + Ω m 0 a − 3 + Ω r 0 a − 4 {\displaystyle {\frac {1}{a}}{\frac {da}{dt}}=H_{0}{\sqrt {\Omega _{\mathrm {\Lambda } 0}+\Omega _{m0}a^{-3}+\Omega _{r0}a^{-4}}}} により与えられる。ここに H 0 {\displaystyle H_{0}} はハッブル定数、 Ω x 0 {\displaystyle \Omega _{x0}} は密度パラメータである。これにより、逆に独立変数として時刻 t {\displaystyle t} の代わりにスケール因子 a {\displaystyle a} を用いることができる。 粒子の座標としては宇宙膨張の効果を取り除いた共動座標 x {\displaystyle \mathbf {x} } が用いられる。これは固有座標 r {\displaystyle \mathbf {r} } と r = a ( t ) x {\displaystyle \mathbf {r} =a(t)\mathbf {x} } という関係にある。粒子の速度はその微分 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}}}} である。 無限に広い計算領域を実現することは不可能であるため、宇宙論的シミュレーションでは周期境界条件が採用される。通常、計算領域は立方体であり、その一片の長さを L {\displaystyle L} とするとき、座標 x {\displaystyle x} と x ± L {\displaystyle x\pm L} の点は同一視される。 周期境界条件のもとでは隣接する立方体に自らの構造のコピーが存在するため、それが及ぼす重力を考慮する必要がある。これは分子動力学法においてクーロン電場に対して開発されたエバルトの方法を適用することで可能であり、付近のボックスからの重力は直接計算し、遠方のボックスからの重力をフーリエ級数の形で取り入れることにより効率的に精度よく計算される。例えば、座標原点にある質量 m {\displaystyle m} の質点がつくる（規格化された）重力ポテンシャルは、周期境界条件のもとで Φ