N体シミュレーションのソースを表示

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

[[File:Galaxy_cluster_sim.png|thumb|320px|<math>N</math>体シミュレーションを用いて得られた[[銀河団]]における質量分布。]]

ナイーブなN体シミュレーションの実装（直接総和法）は重力相互作用の計算に ''O''(''N''<sup>2</sup>) のコストを要するため、より大規模かつ長時間にわたるシミュレーションを実現することは[[計算機科学]]特に[[高性能計算]] (HPC) の分野においても興味深い問題であり、複数の[[ゴードン・ベル賞]]が''N''体シミュレーションの研究に対して与えられた<ref>{{Cite journal |last1=Karp |first1=A. H. |last2=Geist |first2=A. |last3=Bailey |first3=D. |title=1996 Gordon Bell Prize Winners |journal=Computer |volume=30 |issue=1 |date=1997 |page=80-85 |doi=10.1109/MC.1997.562930}}</ref><ref>{{Cite web|和書|url=https://www2.ccs.tsukuba.ac.jp/Astro/news/ja/2012/11/16/ishiyama-gordon-bell/ |title=石山智明 ゴードン・ベル賞を受賞 |accessdate=2020-05-24}}</ref>。現在でも''N''体シミュレーションはコンピュータの[[ベンチマーク]]のためにしばしば用いられる<ref>{{Cite web |date= |url=https://www.hpc-technologies.co.jp/gpu-nbody-benchmark |title=GPU Benchmark - nbody |accessdate=2020-05-24}}</ref>。

== 概要 ==
[[File:Three-body Problem Animation.gif|thumb|320px|三体問題の数値積分は最も単純な<math>N</math>体シミュレーションと言える。]]

[[太陽系]]における[[惑星]]の運動<ref>{{Cite journal |last1=Kinoshita |first1=Hiroshi |last2=Nakai |first2=Hiroshi |title=Long-Term Behavior of the Motion of Pluto over 5.5 Billion Years |journal=Earth Moon and Planets |date=1996 |volume=74 |pages=165-173 |doi=10.1007/BF00117514}}</ref>、[[原始惑星系円盤]]における[[微惑星]]の集積過程<ref>{{Cite journal |last1=Kokubo |first1=E. |last2=Ida |first2=S. |title=On Runaway Growth of Planetesimals |journal=Icarus |volume=123 |pages=180-191 |date=1996}}</ref>、[[球状星団]]の構造の力学的な進化<ref>{{Cite journal |last=Makino |first=Junichiro |title=Postcollapse Evolution of Globular Clusters |journal=The Astrophysical Journal |date=1996 |volume=471 |page=796 |doi=10.1086/178007 |arxiv=astro-ph/9608160}}</ref>、[[宇宙の大規模構造]]の形成<ref>{{Cite journal |author=Springel, V., et al. |title=The Aquarius Project: the subhaloes of galactic haloes |journal=Monthly Notices of the Royal Astronomical Society |date=2008 |volume=391 |pages=1685-1711 |doi=10.1111/j.1365-2966.2008.14066.x |arxiv=0809.0898}}</ref>といった問題は、系を構成する要素（惑星や恒星、[[暗黒物質]]など）が互いに[[重力]]相互作用を及ぼし合うものである。その進化を定量的に求めることは[[天文学]]および[[宇宙物理学]]の多くの分野で重要であるが、それは解析的な手法では困難であり、[[数値解析|数値シミュレーション]]に頼らざるを得ない<ref name="astro-dic"/>。<math>N</math>体シミュレーションはこのような系の数値シミュレーション手法のひとつであり、系を構成する要素を多数の（<math>N</math>個の）粒子とみなし、それらの粒子が互いにニュートン重力を及ぼし合うものとモデル化するものである。

典型的な<math>N</math>体シミュレーションでは、<math>i = 1, 2, \cdots, N</math> 番目の粒子（質量 <math>m_i</math>) の[[運動方程式]]は
:<math>m_i \frac{ d^2 \mathbf{r}_i }{ d t^2 } = - \sum_{j \neq i} G m_i m_j \frac{ \mathbf{r}_i - \mathbf{r}_j }{ \left| \mathbf{r}_i - \mathbf{r}_j \right|^3 }</math>
であり、これを適切な初期条件のもとで数値的に積分することが主たる目標となる{{sfn|牧野|2001|p=399-400}}（[[宇宙の大規模構造]]など[[宇宙論]]的なスケールでの進化を扱う場合、後述のように方程式が修正されるが、本質的な差異はない）。その最も単純かつ非自明な問題が[[三体問題]]である。

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

* 太陽系シミュレーションの場合、長時間の進化の結果として最終的にどのような状態になるのかに興味がある。この場合、数値誤差の累積が最大の問題であり、高次のシンプレクティック積分を用いるなどの誤差を極力抑えるための工夫が必要となる。また、この場合中心天体（[[太陽]]）の重力が支配的で、それに対して摂動が加わった系であるとみなせるため、[[天体力学]]的手法など特有の手法が用いられる。<ref>{{Cite journal |last1=Saha |first1=Prasenjit |last2=Tremaine |first2=Scott |title=Symplectic Integrators for Solar System Dynamics |journal=Astronomical Journal |date=1992 |volume=104 |oage=1633 |doi=10.1086/116347}}</ref>
* 球状星団の場合、重力相互作用のコストが重い上に、近接散乱や連星形成などの効果を正しく計算することが必要である{{sfn|牧野|2001|p=399-400,413-414}}。初期の[[GRAPE]]は球状星団の力学進化の計算に用いられ、大きな成果を上げた。
* 銀河またはより大きなスケールを扱う問題の場合、これは無衝突系であるため <math>N</math> 体粒子は真の粒子である必要はなく、ある空間領域に存在する[[暗黒物質]]あるいは[[バリオン]]を束ねたものである。一方で、より計算領域を大きく、またシミュレーションの分解能を小さくすることが必要である。このため[[スーパーコンピュータ]]などの大規模並列計算機が採用され、また重力相互作用の計算コストを下げるためにアルゴリズム上の改善が続けられている。

=== 無衝突系と衝突系 ===
[[File:The globular cluster NGC 6388 observed by the European Southern Observatory.jpg|thumb|320px|[[球状星団]]は典型的な重力多体系（衝突系）であり、その力学的進化の理論的研究において<math>N</math>体シミュレーションが重要な役割を果たした。]]
<math>N</math> 体シミュレーションはその対象によって大きく無衝突系 (collisionless system) と衝突系 (collisional system) に分類される{{sfn|牧野,福重ほか|2007|p=24-25}}<ref name="BinneyTremaine122">{{harvnb|Binney & Tremaine|p=122-123}}</ref>。これは[[二体緩和]]の効果が重要かどうかを意味し、それによってシミュレーションの性質が大きく変化する。

[[二体緩和]]とは、重力多体系において二体間の近接散乱による系の熱的な進化のことを言う<ref>{{天文学辞典 |urlname=two-body-relaxation |title=2体緩和}}</ref>。二体緩和による系の進化に要する時間スケール <math>t_\mathrm{relax}</math> は、粒子数 <math>N</math> と crossing-time <math>t_\mathrm{cross}</math> を用いて
:<math>t_\mathrm{relax} = \frac{ N }{ \ln N } t_\mathrm{cross}</math>
と書ける{{sfn|牧野,福重ほか|2007|p=27}}。それ故に、極めて粒子数の大きなスケールでは二体緩和が効く時間は[[宇宙年齢]]よりも長くなり、その効果を無視することができる。例えば[[銀河]]は <math>N \sim 10^{10} </math> 個の恒星からなる系であり、その力学的な進化には二体緩和の効果は重要ではないと考えられている。一方、[[球状星団]]では <math>N \sim 10^{4-6}</math> であり、二体緩和が重要である<ref>{{Cite book |last=Spitzer |first=Lyman S. |year=1988 |title=Dynamical Evolution of Globular Clusters |publisher=Princeton University Press |isbn=978-0691084602}}</ref>。

二体緩和が効かない無衝突系では粒子数無限大の極限に相当する<ref name="BinneyTremaine122"/>。このとき重力場はある種の平均場として扱うことが可能であり、個々の粒子に起因する特異性を考慮する必要はない{{sfn|Binney & Tremaine|p=124}}。そのため、シミュレーションに際してツリー法などのより効率的なスキームを使用することができる<ref>Aarseth, p.94.</ref>。また、シミュレーションで扱われる粒子は真の粒子ではなく、[[位相空間 (物理学)|位相空間]]のある領域を代表する点であると解釈される<ref name="BinneyTremaine122"/>。

== 時間積分 ==
=== 積分子 ===
<math>N</math> 体シミュレーションはエネルギーが保存するため、時間積分子として[[リープ・フロッグ法]]などの[[シンプレクティック数値積分法]]がしばしば採用される。例えば太陽系の高精度シミュレーション<ref>{{Cite journal |last1=Wisdom |first1=Jack |last2=Holman |first2=Matthew |title=Symplectic maps for the N-body problem |journal=Astronomical Journal |volume=102 |page=1528-1538 |date=1991 |doi=10.1086/115978}}</ref>、宇宙論的構造形成<ref name="Springel2001"/>など。

衝突系では後述する可変時間刻みと相性の良い[[予測子修正子法]]{{sfn|牧野|2001|p=402-403}}や[[エルミート積分子]]{{sfn|牧野|2001|p=404-406}}も用いられる。

=== 可変時間刻みと独立時間刻み ===
自己重力系は一般に[[重力不安定性]]により[[密度]]ゆらぎが成長し、高密度領域を形成するように進化する。その結果、高密度領域の中心部では{{仮リンク|自由落下時間|en|Free-fall time}}
:<math>t_\mathrm{free} = \frac{ 1 }{ \sqrt{ G \rho } }</math>
が急速に短くなるため, 精度の良いシミュレーションを行うためには時間積分のタイムステップを動的に小さく調整することが望ましい{{sfn|牧野|2001|p=400}}。このような手法は{{仮リンク|可変時間刻み|en|Adaptive step size}}と呼ばれる{{sfn|Binney & Tremaine|p=205}}。

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

== 重力相互作用 ==
特に無衝突系においてはシミュレーションの規模 (つまり粒子数 <math>N</math>) を大きくすることが重要である。しかし粒子 <math>i</math> に作用する重力
:<math>\mathbf{F}_i = - \sum_{j \neq i} G m_i m_j \frac{ \mathbf{r}_i - \mathbf{r}_j }{ \left| \mathbf{r}_i - \mathbf{r}_j \right|^3 }</math>
をすべての粒子について愚直に計算する（これを直接総和法 direct summation という）ことは <math>\mathcal{O} ( N^2 )</math> という大きな計算時間を要するため、粒子数 <math>N</math> を大きくすると急速に計算時間が増大し、現実的な時間で計算を終えることができなくなる{{sfn|Mo, van den Bosch & White|p=766}}。このため、PM法とツリー法という重力計算の精度を下げてでもより効率的な相互作用の計算アルゴリズムが開発された。現在ではこれらの方法を組み合わせた P<sup>3</sup>M 法やtree-PM 法が大規模シミュレーションにおいて標準的な方法として採用されている<ref name="TreePM"/>。

=== Particle-Mesh 法 ===
Particle-Mesh (PM) 法は[[高速フーリエ変換]]に基づいて[[重力ポテンシャル]]の計算を行う{{sfn|Mo, van den Bosch & White|p=766-767}}。まず最初に計算領域に[[格子|グリッド]]を生成し、各頂点での密度の値をその近傍の粒子分布に基づいて決定する。重力ポテンシャルは (フーリエ空間では) [[ラプラス方程式]]
:<math>- k^2 \Phi_\mathbf{k} = 4 \pi G \rho_\mathrm{k}</math>
により密度場に結びついているため、密度場のフーリエ係数 <math>\rho_\mathbf{k}</math> を求め、そこから逆フーリエ変換することにより重力ポテンシャル <math>\Phi</math> や重力場 <math>- \nabla \Phi</math> を求めることができる。この計算時間はグリッド数を <math>M</math> とすると <math>\mathcal{O} ( M \ln M )</math> である。

なお近くの粒子からの重力は直接法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを Particle-Particle Particle-Mesh (P<sup>3</sup>M) 法と呼ぶ{{sfn|Binney & Tremaine|p=135}}{{sfn|Mo, van den Bosch & White|p=767}}。

=== ツリー法 ===
[[File:Barnes-Hut tree construction.png|frame|ツリー構築のアニメーション。]]
Barnes & Hut (1986) により提案された、粒子分布をツリー構造という形で保持することにより重力相互作用の計算を <math>\mathcal{O} ( N \ln N )</math> で行うアルゴリズムは現在 Barnes & Hut のツリー法として広く用いられている<ref name="BarnesHut1986">{{Cite journal |last1=Barnes |first1=Josh |last2=Hut |first2=Piet |title=A hierarchical O(N log N) force-calculation algorithm |journal=Nature |volume=324 |issue=6096 |pages=446-449 |year=1986 |doi=10.1038/324446a0}}</ref>。これは計算領域を表す立方体を階層的により小さな立方体に分割し、最終的に各立方体がひとつ以下の粒子しか含まないようにすることにより、粒子分布の情報を[[木構造 (データ構造)|ツリー]]として保持するものである。ツリーの深さは <math>\mathcal{O} ( \ln N )</math> であるため、ツリーの構成に要する計算量は <math>\mathcal{O} ( N \ln N )</math> である<ref>Aarseth, p.94-96.</ref>。

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

なお近くの粒子からの重力はツリー法で、遠方の粒子からの寄与は PM 法で計算する複合的な手法のことを tree-PM 法と呼ぶ<ref name="TreePM">{{天文学辞典 |urlname=tree-pm-method |title=ツリーPM法}}</ref>。

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

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

Burdet-Heggie正則化<ref>{{Cite journal |last1=Burdet |first1=Claude A. |title=Regularization of the two body problem |journal=Zeitschrift für angewandte Mathematik und Physik ZAMP |volume=18 |issue=3 |page=434-438 |date=1967 |doi=10.1007/BF01601283 }}</ref><ref>{{Cite journal |last1=Burdet |first1=Claude Alain |title=Theory of Kepler motion: the general perturbed two body problem |journal=Zeitschrift für angewandte Mathematik und Physik ZAMP |volume=19 |issue=2 |page=345-368 |date=1968 |doi=10.1007/BF01601478}}</ref><ref>{{Cite journal |last1=Heggie |first1=D. C. |title=Regularization Using a Time-Transformation Only |journal=Astrophysics and Space Science Library |volume=39 |page=34 |date=1973 |doi=10.1007/978-94-010-2611-6_3}}</ref>は、時間座標 <math>t</math> を近接粒子の距離に応じて調整することで特異性を除去するもので、新しい時間座標 <math>\tau</math> を二体の粒子間距離 <math>r</math> と真の時間 <math>t</math> から
:<math>d\tau = \frac{ d t }{ r }</math>
により定義するものである<ref name="BinneyTremaine208">{{harvnb|Binney & Tremaine|p=208-210}}</ref>。このとき二体の相対位置ベクトル <math>\mathbf{r}</math> の従う方程式は
:<math>\frac{ d^2 \mathbf{r} }{ d \tau^2 } - 2 E_2 \mathbf{r} = - \mathbf{e} + r^2 \mathbf{g}</math>
へと帰着される。ここに <math>E_2</math> は二体の相対運動エネルギー、<math>\mathbf{e}</math> は[[ルンゲ＝レンツベクトル|離心率ベクトル]]、<math>\mathbf{g}</math> は他の天体による重力場である。この表式からわかるように、BH正則化により <math>r \to 0</math> での特異性が除去される<ref name="BinneyTremaine208"/>。

[[1965年]]に提案された[[クスターンハイモ・シュティーフェル変換]] (Kustaanheimo-Stiefel transformation)<ref>{{Cite journal |last1=Kustaanheimo |first1=P. |last2=Stiefel |first2=E. |title=Perturbation theory of Kepler motion based on spinor regularization |journal=J. Reine Angew. Math. |volume=218 |date=1965 |page=204-219}}</ref> は3次元直交座標 <math>x, y, z</math> を4次元のスピノルへと変換するもので、BH正則化よりも精度の良い結果が得られる<ref name="BinneyTremaine210">{{harvnb|Binney & Tremaine|p=210-211}}</ref>。
:<math>u_1^2 = \frac{ x + r }{ 2 } \cos^2 \psi , \ \ u_2 = \frac{ y u_1 + z u_4 }{ x + r } , \ u_4^2 = \frac{ x + r }{ 2 } \sin^2 \psi , \ u_3 = \frac{ z u_1 - y u_4 }{ x + r }</math>

=== 近距離カットオフ ===
無衝突系においては、<math>N</math>体粒子は真の粒子ではなく、多数の粒子が占める[[位相空間 (物理学)|位相空間]]上の領域を表す。そのため、<math>N</math>体粒子間に働く重力が近距離で発散する効果は物理的ではなく、適当なカットオフにより除去される必要がある<ref name="BinneyTremaine123">{{harvnb|Binney & Tremaine|p=123-125}}</ref>。最も簡単なカットオフとしては重力ポテンシャル <math>\Phi</math> を
:<math>\Phi ( r ) = - \frac{ G M }{ \sqrt{ r^2 + \epsilon^2 } }</math>
へと変更するものである<ref name="BinneyTremaine123"/>。ここに <math>\epsilon</math> は距離の次元を持つ定数で、この距離スケールより近距離での重力の発散を抑える効果を持つ。このポテンシャルは<math>N</math>体粒子が[[プラマーモデル]]であるような質量分布を持つと仮定した場合に得られるものに等しい。<math>\epsilon</math> の値は平均粒子間距離のオーダーに選ばれる<ref name="BinneyTremaine123"/>。

== 宇宙論的<math>N</math>体シミュレーション ==
[[File:LCDM.jpg|thumb|320px|Gadgetを用いた宇宙論的構造形成シミュレーションのスナップショット。色は[[暗黒物質]]の密度分布を表す。[[銀河フィラメント|フィラメント]]構造や[[ダークマターハロー]]が形成されていることが確認できる。]]

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

=== 運動方程式 ===
宇宙膨張は、宇宙の現在の大きさを1とする相対的な大きさを表す[[スケール因子]] <math>a ( t )</math> により表され、その時間発展は[[フリードマン方程式]]
:<math>\frac{ 1 }{ a } \frac{ d a }{ d t } = H_0 \sqrt{ \Omega_{\mathrm{\Lambda}0} + \Omega_{m0} a^{-3} + \Omega_{r 0} a^{-4} }</math>
により与えられる。ここに <math>H_0</math> は[[ハッブル定数]]、<math>\Omega_{x 0}</math> は[[宇宙論パラメータ|密度パラメータ]]である。これにより、逆に独立変数として時刻 <math>t</math> の代わりにスケール因子 <math>a</math> を用いることができる。

粒子の座標としては宇宙膨張の効果を取り除いた[[共動座標]] <math>\mathbf{x}</math> が用いられる。これは固有座標 <math>\mathbf{r}</math> と
:<math>\mathbf{r} = a ( t ) \mathbf{x}</math>
という関係にある{{sfn|Mo, van den Bosch & White|p=163}}。粒子の速度はその微分 <math>\frac{ d \mathbf{r} }{ d t } = H \mathbf{r} + a ( t ) \frac{ d \mathbf{x} }{ d t }</math> であるが、初期宇宙での発散を回避するために <math>\mathbf{w} := \sqrt{ a ( t ) } \frac{ d \mathbf{x} }{ d t }</math> が用いられる<ref name="Springel2001"/>。最終的な運動方程式は
:<math>\frac{ d \mathbf{x} }{ d a } = \frac{ \mathbf{w} }{ S ( a ) }</math>
:<math>\frac{ d \mathbf{w} }{ d a } = - \frac32 \frac{ \mathbf{w} }{ a } + \frac{ 1 }{ a^2 S ( a ) } \nabla \Phi ( \mathbf{x} )</math>
::<math>S ( a ) = H_0 \sqrt{ \Omega_{m 0} + a^3 \Omega_{\Lambda 0} }</math>
である<ref name="Springel2001"/>。

=== 周期境界条件 ===
無限に広い計算領域を実現することは不可能であるため、宇宙論的シミュレーションでは[[周期境界条件]]が採用される{{sfn|Mo, van den Bosch & White|p=769}}。通常、計算領域は立方体であり、その一片の長さを <math>L</math> とするとき、座標 <math>x</math> と <math>x \pm L</math> の点は同一視される。

周期境界条件のもとでは隣接する立方体に自らの構造のコピーが存在するため、それが及ぼす重力を考慮する必要がある。これは[[分子動力学法]]において[[クーロンの法則|クーロン電場]]に対して開発された[[エバルトの方法]]を適用することで可能であり、付近のボックスからの重力は直接計算し、遠方のボックスからの重力を[[フーリエ級数]]の形で取り入れることにより効率的に精度よく計算される<ref>{{Cite journal |last1=Hernquist |firtst1=Lars |last2=Bouchet |first2=Francois R. |last3=Suto |first3=Yasushi |title=Application of the Ewald Method to Cosmological N-Body Simulations |journal=Astrophysical Journal Supplement |volume=75 |pages=231 |doi=10.1086/191530}}</ref>。例えば、座標原点にある質量 <math>m</math> の質点がつくる（規格化された）重力ポテンシャルは、周期境界条件のもとで
:<math>\hat{\Phi} ( t, {\mathbf{x}}) = - \frac{ G m }{ a } \sum_{{\mathbf{n}}\in \mathbf{Z}^3} \frac{ {{\mathrm{erfc}\,}}( \alpha | {\mathbf{x}}- {\mathbf{n}}L | ) }{ | {\mathbf{x}}- {\mathbf{n}}L | } - \frac{ G m }{ \pi a L } \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{ G m }{ a } \frac{ \pi }{ \alpha^2 L^3 }
</math>
となる。ここに <math>\alpha</math> は任意の正の定数、<math>\mathrm{erfc}</math> は相補[[誤差関数]]である。

=== 初期条件 ===
[[Λ-CDMモデル]]では宇宙論的ゆらぎは[[宇宙のインフレーション|インフレーション期]]に[[ガウス過程|ガウス分布]]に従って生成されたと考えられており、線型摂動の範囲では密度ゆらぎのパワースペクトルが指定されれば初期条件を確率的に生成することができる。パワースペクトルは与えられた宇宙論パラメータのもとで[[宇宙論的摂動論]]に基づいて計算することができる。なお宇宙論的<math>N</math>体シミュレーションで最も広く用いられる初期条件は2次の[[ラグランジュ摂動]] (2LPT) に基づくものである。<ref>{{Cite journal |last1=Hahn |first1=Oliver |last2=Abel |first2=Tom |title=Multi-scale initial conditions for cosmological simulations |journal=Monthly Notices of the Royal Astronomical Society |date=2011| volume=415 |pages=2101-2121 |doi=10.1111/j.1365-2966.2011.18820.x |arxiv=1103.6031}}</ref>

== プロジェクト ==
=== ソフトウェア ===
{{仮リンク|Volker Springel|de|Volker Springel}}が中心になって開発された[[:en:GADGET]]は銀河や宇宙論的構造形成を主なターゲットとする無衝突系のシミュレーションコードであり<ref name="gadgetWeb">{{Cite web |url=https://wwwmpa.mpa-garching.mpg.de/gadget/ |title=Cosmological simulations with GADGET |publisher=Volker Springel |accessdate=2020-05-24}}</ref>、[[1998年]]にバージョン1<ref name="Springel2001">{{cite journal |last1=Springel |first1=Volker |last2=Yoshida |first2=Naoki |last3=White |first3=Simon D.M. |date=April 2001 |title=GADGET: a code for collisionless and gasdynamical cosmological simulations |journal=New Astronomy |volume=6 |issue=2 |pages=79-117 |url=http://www.mpa-garching.mpg.de/gadget/gadget1-paper.pdf |accessdate=January 21, 2018 |doi=10.1016/S1384-1076(01)00042-2 |arxiv=astro-ph/0003162 |bibcode=2001NewA....6...79S }}</ref>が、[[2005年]]にバージョン2<ref name="Springel2005">{{cite journal |last=Springel |first=Volker |date=21 December 2005 |title=The cosmological simulation code GADGET-2 |journal=Monthly Notices of the Royal Astronomical Society |volume=364 |issue=4 |pages=1105-1134 |url=http://www.mpa-garching.mpg.de/gadget/gadget2-paper.pdf |accessdate=January 21, 2018 |doi=10.1111/j.1365-2966.2005.09655.x |arxiv=astro-ph/0505010 |bibcode=2005MNRAS.364.1105S }}</ref>が公開された。[[スーパーコンピュータ]]などの大規模並列計算機で動かすために[[並列計算|並列化]]されており<ref name="gadgetWeb"/>、[[C言語]]によって実装されている。ライセンスは[[GNU General Public License|GNU GPL]]。

2010年代には[[牧野淳一郎]]らが中心となって汎用的な多体問題シミュレーションフレームワークである[[FDPS]]が開発されている<ref>{{Cite journal |和書|publisher= |分子シミュレーション学会 |url=https://doi.org/10.11436/mssj.19.105 |title=連載「FDPS入門(1)」 |author=牧野淳一郎 |journal=アンサンブル |volume=19 |issue=2 |pages=105-110 |year=2017 |doi=10.11436/mssj.19.105 |accessdate=2020-06-23}}</ref>。

=== ハードウェア ===
[[GRAPE]]は[[杉本大一郎]]、[[牧野淳一郎]]らによって[[東京大学]]で開発された<math>N</math>体シミュレーション専用計算機であり、重力相互作用の計算を[[パイプライン処理|パイプライン]]として物理的に実装することにより効率的に計算を進めるものである<ref>{{Cite journal |last1=Sugimoto |first1=Daiichiro |last2=Chikada |first2=Yoshihiro |last3=Makino |first3=Junichiro |last4=Ito |first4=Tomoyoshi |last5=Ebisuzaki |first5=Toshikazu |last6=Umemura |first6=Masayuki |title=A special-purpose computer for gravitational many-body problems |journal=Nature |volume=345 |issue=6270 |page=33-35 |date=1990 |doi=10.1038/345033a0}}</ref>{{sfn|牧野,福重ほか|2007|p=57-61}}。現在も[[国立天文台]]などで運用されている<ref>{{Cite web|和書|date= |url=https://www.cfca.nao.ac.jp/content/grapegpu%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0 |title=GRAPE/GPUシステム - CfCA - Center for Computational Astrophysics |publisher=国立天文台 |accessdate=2020-05-24}}</ref>。

また [[Graphics Processing Unit|GPU]]の活用（[[GPGPU]]）も進められている<ref>{{Cite web|和書|url=http://www.sim.gsic.titech.ac.jp/Japanese/Research/gpu.html |title=青木研究室（東工大）：Research：ＧＰＵコンピューティング |accessdate=2020-05-24}}</ref><ref>{{Cite journal |last1=Miki |first1=Yohei |last2=Takahashi |first2=Daisuke |last3=Mori |first3=Masao |title=A Fast Implementation and Performance Analysis of Collisionless N-body Code Based on GPGPU |journal=Procedia Computer Science |volume=9 |pages=96-105 |date=2012 |doi=10.1016/j.procs.2012.04.011}}</ref>。

===シミュレーション ===
21世紀に入ってから、特に大規模なシミュレーションはスーパーコンピュータを長時間占有する必要があるため、大規模なシミュレーションを行いその出力を公開するプロジェクトが行われるようになった。その有名なものが{{仮リンク|ミレニアム・シミュレーション|en|Millennium Run}}<ref>{{Cite web |url=https://wwwmpa.mpa-garching.mpg.de/galform/virgo/millennium/ |title=The Millennium Simulation Project |accessdate=2020-05-24}}</ref>であり、他に The Aquarius Project<ref>{{Cite web |url=https://wwwmpa.mpa-garching.mpg.de/aquarius/ |title=The Aquarius Project |accessdate=2020-05-24}}</ref> などがある。2010年代に実行された[[:en:Illustris project]]<ref>{{Cite web |url=https://www.illustris-project.org/ |title=The Illustris Simulation |accessdate=2020-05-24}}</ref>は<math>N</math>体シミュレーションだけでなく、[[星形成]]や[[活動銀河|AGN]]といったバリオン物理を考慮した流体シミュレーションを行っている。

== 歴史 ==
重力多体系の[[計算機]]を用いた研究すなわち<math>N</math>体シミュレーションは1960年代から実際的な研究で採用されるようになった<ref>{{Cite web|和書
|url  = http://th.nao.ac.jp/MEMBER/tanikawa/21seiki/cm21j.pdf
|author  = 伊藤孝士, 谷川清隆
|title  = 21世紀の天体力学
|accessdate = 2020-05-24}} p.7.</ref>。例えば[[1963年]]の[[:en:Sverre Aarseth]]による <math>N = 100</math> 体のシミュレーション<ref>{{Cite journal |last=Aarseth |first=S. J. |title=Dynamical evolution of clusters of galaxies, I |journal=Monthly Notices of the Royal Astronomical Society |volume=126 |page=223 |date=1963 |doi=10.1093/mnras/126.3.223}}</ref>、[[1964年]]の Hénon & Heiles による[[第三積分]]に関する数値的研究<ref>{{cite journal |author1=Hénon, M. |author2=Heiles, C. |title=The applicability of the third integral of motion: Some numerical experiments |journal=The Astronomical Journal |volume=69 |pages=73-79 |year=1964 |bibcode=1964AJ.....69...73H |doi = 10.1086/109234 }}</ref>（これは <math>N = 1</math> であるが）、[[1973年]]の Hénon による多体系の安定性の研究<ref>{{Cite journal |last=Henon |first=M. |title=Numerical Experiments on the Stability of Spherical Stellar Systems |journal=Astronomy and Astrophysics |volume=24 |page=229 |date=1973}}</ref>など。[[ジェームズ・ピーブルス]]は[[1970年]]に <math>N = 300</math> 体を用いて銀河団形成過程のシミュレーションを行った<ref>{{Cite journal |last=Peebles |first=P. J. E. |title=Structure of the Coma Cluster of Galaxies |journal=The Astronomical Journal |date=1970 |doi=10.1086/110933}}</ref>。その後も[[1979年]]には Efstathiou & Jones が <math>N = 500</math> 体による[[銀河回転]]の研究<ref>{{Cite journal |last1=Efstathiou |first1=G. |last2=Jones |first2=B. J. T. |title=The rotation of galaxies: numerical investigations of the tidal torque theory |journal=Monthly Notices of the Royal Astronomical Society |volume=186 |pages=133-144 |date=1979 |doi=10.1093/mnras/186.2.133}}</ref>など、計算機の発達に伴ってより大規模なシミュレーションがなされるようになっていった。

より大規模なシミュレーションの要求は強く、[[1986年]]に Barnes & Hut<ref name="BarnesHut1986"/> はツリー法を導入し、同時期に PM 法も確立した<ref>{{cite book|title=Computer simulation using particles |author1=Roger W. Hockney|author2=James W. Eastwood |chapter=Particle-Particle-Particle-Mesh (P<sup>3</sup>M) Algorithms |pages=[https://archive.org/details/computersimulati0000hock/page/267 267&ndash;304] |publisher=CRC Press |year=1988 |isbn=9780852743928 |url=https://archive.org/details/computersimulati0000hock/page/267}}</ref>。[[1989年]]には[[GRAPE]]プロジェクトがスタートしている。

一方で積分スキームに関する研究も進められた。[[天体力学]]の分野からは[[1990年]]に[[吉田春夫]]によりシンプレクティック積分子の一般的な構成方法が示された<ref>{{cite journal |last = Yoshida |first = H. |title = Construction of higher order symplectic integrators |journal=Phys. Lett. A |year=1990 |volume=150 |pages=262 |doi = 10.1016/0375-9601(90)90092-3|bibcode = 1990PhLA..150..262Y |issue = 5-7 }}</ref>。その翌年牧野はエルミート積分子を導入した<ref>{{Cite journal |last=Makino |first=Junichiro |title=A Modified Aarseth Code for GRAPE and Vector Processors |journal=Publications of the Astronomical Society of Japan |date=1991 |volume=43 |pages=859-876}}</ref>。やがて[[対称型公式]]の有用性が認められるようになった<ref>{{Cite journal
|last1=Kokubo |first1=Eiichiro |last2=Yoshinaga |first2=Keiko |last3=Makino |first3=Junichiro |title=On a time-symmetric Hermite integrator for planetary N-body simulatio |journal=Monthly Notices of the Royal Astronomical Society |date=1998 |volume=297 |pages=1067-1072 |doi=10.1046/j.1365-8711.1998.01581.x}}</ref>。

[[2005年]]の{{仮リンク|ミレニアム・シミュレーション|en|Millennium Run}}では <math>N = 1.0 \times 10^{10}</math> の宇宙論的構造形成シミュレーションが遂行されるに至った<ref>{{Cite journal | last = Springel |first = Volker |display-authors=etal | date = 2005 | title = Simulations of the formation, evolution, and clustering of galaxies and quasars | journal = [[Nature (journal)|Nature]] | volume = 435 |issue = 7042 | pages = 629-636 | arxiv = astro-ph/0504097 | bibcode=2005Natur.435..629S | doi= 10.1038/nature03597 | pmid = 15931216 |hdl = 2027.42/62586 |url=https://deepblue.lib.umich.edu/bitstream/2027.42/62586/1/nature03597.pdf}}</ref>。

== 脚注 ==
{{Reflist|30em}}

== 参考文献 ==
* {{Cite book
|last  = Binney
|first  = James
|last2  = Tremaine
|first2 = Scott
|year  = 2008
|title  = Galactic Dynamics
|edition = Second
|publisher = Princeton University Press
|isbn  = 978-0-691-13027-9 |ref={{harvid|Binney & Tremaine}}}}
* {{Cite book
|last  = Aarseth
|first  = Sverre J.
|year  = 2010
|title  = Gravitational N-Body Simulations: Tools and Algorithms (Cambridge Monographs on Mathematical Physics)
|publisher = Cambridge University Press
|isbn  = 978-0521121538
}}
* {{Cite book
|last1  = Mo
|first1 = Houjun
|last2  = van den Bosch
|first2 = Frank
|last3  = White
|first3 = Simon
|title  = Galaxy Formation and Evolution
|year  = 2010
|publisher = Cambridge University Press
|isbn  = 978-0-521-85793-2 |ref={{harvid|Mo, van den Bosch & White}}
}}
* {{Cite journal |和書|author=牧野淳一郎 |title=重力多体系の数値計算(<シリーズ>物性研究者のための計算手法入門) |journal=物性研究 |volume=76 |issue=3 |year=2001 |month=6 |publisher=物性研究刊行会 |url=https://hdl.handle.net/2433/97012|2020-06-23 |ref={{harvid|牧野|2001}}}}
* {{Cite web|和書 |url=https://www.cfca.nao.ac.jp/~cfca/hpc/muv/text/note.pdf |format=PDF |title=N体シミュレーション啓蟄の学校教科書 |author=牧野淳一郎, 福重俊幸, 小久保英一郎, 川井敦, 台坂博, 杉本大一郎 |publisher=国立天文台 |date=2007-03-13 |accessdate=2020-05-24 |ref={{harvid|牧野,福重ほか|2007}}}}

== 関連項目 ==
* [[計算物理学]]
* [[多体問題]]
* [[三体問題]]
* [[ピタゴラス三体問題]]

== 外部リンク ==
* {{Cite journal |和書|author=小松信義, 木綿隆弘, 木村繁男 |title=重力エネルギーが支配的なN体系の速度緩和 |https://doi.org/10.11345/japannctam.61.0.222.0 |journal=理論応用力学講演会 講演論文集 |publisher=日本学術会議 「機械工学委員会・土木工学・建築学委員会合同IUTAM分科会」 |year=2012 |volume=第61回理論応用力学講演会 |issue=セッションID: GS06-01 |pages=222-222 |naid=130005020725 |doi=10.11345/japannctam.61.0.222.0}}

{{DEFAULTSORT:Nたいしみゆれいしよん}}
[[Category:計算物理学]]
[[Category:天文学に関する記事]]