N体シミュレーション

N体シミュレーション (エヌたいシミュレーション、N-body simulation) とは、天体物理学および天文学において、重力相互作用する"N"個の粒子の力学的な進化を数値的に計算するシミュレーションのことをいう。2体系つまりケプラー問題は可積分であるが、3体以上の系（重力多体系）は可積分ではなく、その力学的進化を定量的に予測するためには数値シミュレーションが必須である。太陽系や球状星団、銀河あるいは宇宙の大規模構造など、重力多体系は宇宙のあらゆる領域において重要な役割を果たすため、"N"体シミュレーションは宇宙に関する理論的研究において極めて重要な役割を果たしている。
ナイーブなN体シミュレーションの実装（直接総和法）は重力相互作用の計算に "O"("N"2) のコストを要するため、より大規模かつ長時間にわたるシミュレーションを実現することは計算機科学特に高性能計算 (HPC) の分野においても興味深い問題であり、複数のゴードン・ベル賞が"N"体シミュレーションの研究に対して与えられた。現在でも"N"体シミュレーションはコンピュータのベンチマークのためにしばしば用いられる。
概要.
太陽系における惑星の運動、原始惑星系円盤における微惑星の集積過程、球状星団の構造の力学的な進化、宇宙の大規模構造の形成といった問題は、系を構成する要素（惑星や恒星、暗黒物質など）が互いに重力相互作用を及ぼし合うものである。その進化を定量的に求めることは天文学および宇宙物理学の多くの分野で重要であるが、それは解析的な手法では困難であり、数値シミュレーションに頼らざるを得ない。formula_1体シミュレーションはこのような系の数値シミュレーション手法のひとつであり、系を構成する要素を多数の（formula_1個の）粒子とみなし、それらの粒子が互いにニュートン重力を及ぼし合うものとモデル化するものである。
典型的なformula_1体シミュレーションでは、formula_4 番目の粒子（質量 formula_5) の運動方程式は
であり、これを適切な初期条件のもとで数値的に積分することが主たる目標となる（宇宙の大規模構造など宇宙論的なスケールでの進化を扱う場合、後述のように方程式が修正されるが、本質的な差異はない）。その最も単純かつ非自明な問題が三体問題である。
シミュレーションの目的によって様々な数値計算上の困難が存在し、各分野毎に様々な手法が開発されてきた。
無衝突系と衝突系.
formula_1 体シミュレーションはその対象によって大きく無衝突系 (collisionless system) と衝突系 (collisional system) に分類される。これは二体緩和の効果が重要かどうかを意味し、それによってシミュレーションの性質が大きく変化する。
二体緩和とは、重力多体系において二体間の近接散乱による系の熱的な進化のことを言う。二体緩和による系の進化に要する時間スケール formula_9 は、粒子数 formula_1 と crossing-time formula_11 を用いて
と書ける。それ故に、極めて粒子数の大きなスケールでは二体緩和が効く時間は宇宙年齢よりも長くなり、その効果を無視することができる。例えば銀河は formula_13 個の恒星からなる系であり、その力学的な進化には二体緩和の効果は重要ではないと考えられている。一方、球状星団では formula_14 であり、二体緩和が重要である。
二体緩和が効かない無衝突系では粒子数無限大の極限に相当する。このとき重力場はある種の平均場として扱うことが可能であり、個々の粒子に起因する特異性を考慮する必要はない。そのため、シミュレーションに際してツリー法などのより効率的なスキームを使用することができる。また、シミュレーションで扱われる粒子は真の粒子ではなく、位相空間のある領域を代表する点であると解釈される。
時間積分.
積分子.
formula_1 体シミュレーションはエネルギーが保存するため、時間積分子としてリープ・フロッグ法などのシンプレクティック数値積分法がしばしば採用される。例えば太陽系の高精度シミュレーション、宇宙論的構造形成など。
衝突系では後述する可変時間刻みと相性の良い予測子修正子法やエルミート積分子も用いられる。
可変時間刻みと独立時間刻み.
自己重力系は一般に重力不安定性により密度ゆらぎが成長し、高密度領域を形成するように進化する。その結果、高密度領域の中心部では
が急速に短くなるため, 精度の良いシミュレーションを行うためには時間積分のタイムステップを動的に小さく調整することが望ましい。このような手法はと呼ばれる。
ところが、特に衝突系で連星形成が起こるような状況では、一部のformula_1体粒子は極めて短い時間で進化するものの、その他の大多数の粒子の軌道進化に小さな時間刻みが必要ないという可能性がある。この場合、最も小さな時間ステップに合わせて全体の時間刻みを調整すると、シミュレーションに多大な時間を要することになり、また不必要に小さな時間ステップに伴う数値積分誤差が累積する可能性がある。そこで必要な粒子のみ小さな時間刻み幅で時間積分を行う独立時間刻みという手法が開発された。この場合、その他の大多数の粒子については適当な補間を用いてその重力場を見積もり、必要な粒子のみ正確に時間積分を行うことになる。
重力相互作用.
特に無衝突系においてはシミュレーションの規模 (つまり粒子数 formula_1) を大きくすることが重要である。しかし粒子 formula_19 に作用する重力
をすべての粒子について愚直に計算する（これを直接総和法 direct summation という）ことは formula_21 という大きな計算時間を要するため、粒子数 formula_1 を大きくすると急速に計算時間が増大し、現実的な時間で計算を終えることができなくなる。このため、PM法とツリー法という重力計算の精度を下げてでもより効率的な相互作用の計算アルゴリズムが開発された。現在ではこれらの方法を組み合わせた P3M 法やtree-PM 法が大規模シミュレーションにおいて標準的な方法として採用されている。
Particle-Mesh 法.
Particle-Mesh (PM) 法は高速フーリエ変換に基づいて重力ポテンシャルの計算を行う。まず最初に計算領域にグリッドを生成し、各頂点での密度の値をその近傍の粒子分布に基づいて決定する。重力ポテンシャルは (フーリエ空間では) ラプラス方程式
により密度場に結びついているため、密度場のフーリエ係数 formula_24 を求め、そこから逆フーリエ変換することにより重力ポテンシャル formula_25 や重力場 formula_26 を求めることができる。この計算時間はグリッド数を formula_27 とすると formula_28 である。
なお近くの粒子からの重力は直接法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを Particle-Particle Particle-Mesh (P3M) 法と呼ぶ。
ツリー法.
Barnes & Hut (1986) により提案された、粒子分布をツリー構造という形で保持することにより重力相互作用の計算を formula_29 で行うアルゴリズムは現在 Barnes & Hut のツリー法として広く用いられている。これは計算領域を表す立方体を階層的により小さな立方体に分割し、最終的に各立方体がひとつ以下の粒子しか含まないようにすることにより、粒子分布の情報をツリーとして保持するものである。ツリーの深さは formula_30 であるため、ツリーの構成に要する計算量は formula_29 である。
ある粒子に作用する重力を計算する際には、遠方の粒子群からの寄与をまとめて計算することによりコストを削減する。この際に各立方体の重心および質量を用いるが、計算の精度を上げるために四重極モーメントなどを用いる場合もある。
なお近くの粒子からの重力はツリー法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを tree-PM 法と呼ぶ。
重力の近距離発散.
ふたつのformula_1体粒子間の距離が極めて小さくなると、両者の間に働く重力は任意に大きくなり得る。これは衝突系においては物理的に重要であり、その影響を正しくシミュレーションする必要がある。一方、無衝突系ではこの効果は物理的ではなく、カットオフにより除去される。
正則化.
衝突系において近距離発散に対処するために可変時間刻み幅を用いる場合、時間ステップが際限なく小さくなり、シミュレーションがほとんど進まなくなってしまう。しかしながら、二体問題における近距離重力に起因する特異性は見かけのものであり、適切な変数変換により除去することができる。この手続きは正則化 (regularization) として知られる。
Burdet-Heggie正則化は、時間座標 formula_33 を近接粒子の距離に応じて調整することで特異性を除去するもので、新しい時間座標 formula_34 を二体の粒子間距離 formula_35 と真の時間 formula_33 から
により定義するものである。このとき二体の相対位置ベクトル formula_38 の従う方程式は
へと帰着される。ここに formula_40 は二体の相対運動エネルギー、formula_41 は離心率ベクトル、formula_42 は他の天体による重力場である。この表式からわかるように、BH正則化により formula_43 での特異性が除去される。
1965年に提案されたクスターンハイモ・シュティーフェル変換 (Kustaanheimo-Stiefel transformation) は3次元直交座標 formula_44 を4次元のスピノルへと変換するもので、BH正則化よりも精度の良い結果が得られる。
近距離カットオフ.
無衝突系においては、formula_1体粒子は真の粒子ではなく、多数の粒子が占める位相空間上の領域を表す。そのため、formula_1体粒子間に働く重力が近距離で発散する効果は物理的ではなく、適当なカットオフにより除去される必要がある。最も簡単なカットオフとしては重力ポテンシャル formula_25 を
へと変更するものである。ここに formula_50 は距離の次元を持つ定数で、この距離スケールより近距離での重力の発散を抑える効果を持つ。このポテンシャルはformula_1体粒子がプラマーモデルであるような質量分布を持つと仮定した場合に得られるものに等しい。formula_50 の値は平均粒子間距離のオーダーに選ばれる。
宇宙論的formula_1体シミュレーション.
大規模構造の形成などの宇宙論的な問題はformula_1体シミュレーションが用いられる典型的かつ重要なセットアップであるが、宇宙膨張を考慮する必要があるという点で他の問題とは異なっている。この種のシミュレーションは非線型成長後の物質の密度ゆらぎのパワースペクトル、あるいはダークマターハローの密度プロファイルや質量関数を求めるために用いられる。これらの量は観測可能量であり、実際の観測データと比較することにより例えば宇宙論パラメータの制限を与えることができる。
運動方程式.
宇宙膨張は、宇宙の現在の大きさを1とする相対的な大きさを表すスケール因子 formula_55 により表され、その時間発展はフリードマン方程式
により与えられる。ここに formula_57 はハッブル定数、formula_58 は密度パラメータである。これにより、逆に独立変数として時刻 formula_33 の代わりにスケール因子 formula_60 を用いることができる。
粒子の座標としては宇宙膨張の効果を取り除いた共動座標 formula_61 が用いられる。これは固有座標 formula_38 と
という関係にある。粒子の速度はその微分 formula_64 であるが、初期宇宙での発散を回避するために formula_65 が用いられる。最終的な運動方程式は
である。
周期境界条件.
無限に広い計算領域を実現することは不可能であるため、宇宙論的シミュレーションでは周期境界条件が採用される。通常、計算領域は立方体であり、その一片の長さを formula_69 とするとき、座標 formula_70 と formula_71 の点は同一視される。
周期境界条件のもとでは隣接する立方体に自らの構造のコピーが存在するため、それが及ぼす重力を考慮する必要がある。これは分子動力学法においてクーロン電場に対して開発されたエバルトの方法を適用することで可能であり、付近のボックスからの重力は直接計算し、遠方のボックスからの重力をフーリエ級数の形で取り入れることにより効率的に精度よく計算される。例えば、座標原点にある質量 formula_72 の質点がつくる（規格化された）重力ポテンシャルは、周期境界条件のもとで
となる。ここに formula_74 は任意の正の定数、formula_75 は相補誤差関数である。
初期条件.
Λ-CDMモデルでは宇宙論的ゆらぎはインフレーション期にガウス分布に従って生成されたと考えられており、線型摂動の範囲では密度ゆらぎのパワースペクトルが指定されれば初期条件を確率的に生成することができる。パワースペクトルは与えられた宇宙論パラメータのもとで宇宙論的摂動論に基づいて計算することができる。なお宇宙論的formula_1体シミュレーションで最も広く用いられる初期条件は2次のラグランジュ摂動 (2LPT) に基づくものである。
プロジェクト.
ソフトウェア.
が中心になって開発されたは銀河や宇宙論的構造形成を主なターゲットとする無衝突系のシミュレーションコードであり、1998年にバージョン1が、2005年にバージョン2が公開された。スーパーコンピュータなどの大規模並列計算機で動かすために並列化されており、C言語によって実装されている。ライセンスはGNU GPL。
2010年代には牧野淳一郎らが中心となって汎用的な多体問題シミュレーションフレームワークであるFDPSが開発されている。
ハードウェア.
GRAPEは杉本大一郎、牧野淳一郎らによって東京大学で開発されたformula_1体シミュレーション専用計算機であり、重力相互作用の計算をパイプラインとして物理的に実装することにより効率的に計算を進めるものである。現在も国立天文台などで運用されている。
また GPUの活用（GPGPU）も進められている。
シミュレーション.
21世紀に入ってから、特に大規模なシミュレーションはスーパーコンピュータを長時間占有する必要があるため、大規模なシミュレーションを行いその出力を公開するプロジェクトが行われるようになった。その有名なものがであり、他に The Aquarius Project などがある。2010年代に実行されたはformula_1体シミュレーションだけでなく、星形成やAGNといったバリオン物理を考慮した流体シミュレーションを行っている。
歴史.
重力多体系の計算機を用いた研究すなわちformula_1体シミュレーションは1960年代から実際的な研究で採用されるようになった。例えば1963年のによる formula_80 体のシミュレーション、1964年の Hénon & Heiles による第三積分に関する数値的研究（これは formula_81 であるが）、1973年の Hénon による多体系の安定性の研究など。ジェームズ・ピーブルスは1970年に formula_82 体を用いて銀河団形成過程のシミュレーションを行った。その後も1979年には Efstathiou & Jones が formula_83 体による銀河回転の研究など、計算機の発達に伴ってより大規模なシミュレーションがなされるようになっていった。
より大規模なシミュレーションの要求は強く、1986年に Barnes & Hut はツリー法を導入し、同時期に PM 法も確立した。1989年にはGRAPEプロジェクトがスタートしている。
一方で積分スキームに関する研究も進められた。天体力学の分野からは1990年に吉田春夫によりシンプレクティック積分子の一般的な構成方法が示された。その翌年牧野はエルミート積分子を導入した。やがて対称型公式の有用性が認められるようになった。
2005年のでは formula_84 の宇宙論的構造形成シミュレーションが遂行されるに至った。