next up previous contents index
次へ: 非線型構造の統計モデル 上へ: 摂動の非線形成長 前へ: 非線型成長の近似モデル   目次   索引

N体シミュレーション

構造形成の非線形性を調べるのにもっとも直接的な方法は,コンピュータ・シ ミュレーションを行うことである.コンピュータ上に仮想的な宇宙を作り, その進化を数値計算によって追っていくのである.コンピュータは原理的 に0と1の組合せによってすべてが表されており,その自由度は有限である.一 方,密度場の自由度は空間の点の数,すなわち無限大である.そこで、シミュ レーションにおいては,無限大の自由度の系を有限の自由度の系で近似するこ とが必要になる.そこである程度小さいスケールの構造を無視して,興味あ るスケールの構造を調べるのである.どのくらい小さい構造まで調べられるか をシミュレーションの解像度という. 

ダークマターによる構造形成の場合にはシミュレーションの原理は比較的単純 である.ダークマターには重力以外の相互作用が働かない.そこで,物質を有 限個数$ N$ の粒子の集まりであると考えて,その粒子間に重力相互作用のみが働 くものとする.ある初期条件のもと,個々の粒子の位置と速度をコンピュータ 上に記憶させておく.ひとつひとつの粒子に働く重力を求めて加速度を計算す る.その加速度と粒子の持つ速度をもとにして時間をすこし進め,粒子が移動 する先の位置と速度を計算する.これを繰り返して,宇宙初期から現在にいた るダークマターの分布構造を具体的につくって,その進化を調べるのである. この手法を用いたシミュレーションを$ N$ 体シミュレーションという.

もう少し具体的に見てみよう. ある時刻$ t$ における粒子$ i$  ( $ i=1,2,\ldots,N$ )の共動座標における位置と 速度をそれぞれ $ {\mbox{\boldmath $x$}}_i$ , $ {\mbox{\boldmath $u$}}_i=\dot{{\mbox{\boldmath $x$}}}_i$ とする.これら $ 6\times N$ の数値 はコンピュータ上に記憶され,位相空間における物質分布を粒子で表現したも のとなっている.ある粒子に着目すると,その粒子の共動座標の速度 $ {\mbox{\boldmath $u$}}$ の時間変化率は次のラグランジュ微分

$\displaystyle \frac{d{\mbox{\boldmath$u$}}}{dt} = \frac{\partial{\mbox{\boldmat...
...mbox{\boldmath$u$}}\cdot{\mbox{\boldmath$\nabla$}}\right) {\mbox{\boldmath$u$}}$ (P.4.23)

で与えられる.ここで左辺の時間微分は着目する粒子の速度変化率を表したラ グランジュ的見方のもとでの微分を表し,右辺の時間に関する偏微分は空間の ある一点を固定してその場所での速度変化率を表したオイラー的見方のもとで の偏微分である.すると,オイラー的見方での運動方程式であるオイラー方程 式(6.1.18)をラグランジュ的見方に直すと,圧力のないダークマ ター($ p=0$ )に対して粒子$ i$ の運動方程式は

$\displaystyle \frac{d{\mbox{\boldmath$u$}}_i}{dt} + 2H{\mbox{\boldmath$u$}}_i = {\mbox{\boldmath$g$}}_i$ (P.4.24)

となる.ここで $ {\mbox{\boldmath $g$}}_i=-{\mbox{\boldmath $\nabla$}}\Phi({\mbox{\boldmath $x$}}_i,t)$ は共動座標における粒子$ i$ の位置での重力加速度 を表す量である.左辺第2項は宇宙膨張による引きずられの効果を表している. 重力加速度が求められれば,時間が少し進んだ $ t + \Delta t$ における時刻で の粒子の位置と速度は
    $\displaystyle {\mbox{\boldmath$x$}}_i(t+\Delta t) =
{\mbox{\boldmath$x$}}_i(t) + {\mbox{\boldmath$u$}}_i(t) \Delta t$ (P.4.25)
    $\displaystyle {\mbox{\boldmath$u$}}_i(t+\Delta t) =
{\mbox{\boldmath$u$}}_i(t)...
...\left[{\mbox{\boldmath$g$}}_i(t) - 2H{\mbox{\boldmath$u$}}_i(t)\right] \Delta t$ (P.4.26)

によって求めることができるP1.これを何度も繰り返すことにより,膨張宇宙の物質分布の時間進化 を非線形段階まで含めて追っていけるのである.

ここで,もっとも重要かつ計算時間のかかるのが重力加速度 $ {\mbox{\boldmath $g$}}_i$ の計算 である.重力は長距離力であることに対応して,重力加速度は局所的には定ま らない.式(6.1.21)で与えられる重力ポテンシャルの勾配から, 重力加速度は

$\displaystyle {\mbox{\boldmath$g$}}_i = G \int d^3\!x  \rho({\mbox{\boldmath$x...
...x$}}_i}{\left\vert{\mbox{\boldmath$x$}} - {\mbox{\boldmath$x$}}_i\right\vert^3}$ (P.4.27)

となる.ここで, $ \bar{\rho}$ の項の寄与は $ {\mbox{\boldmath $x$}} - {\mbox{\boldmath $x$}}_i$ の角度積分 により打ち消しあって消えている.$ N$ 体シミュレーションでは,密度は点粒 子で表現されているので,密度場は各粒子の位置におけるデルタ関数の和:

$\displaystyle \rho({\mbox{\boldmath$x$}},t) = \frac{m}{a^3} \sum_{j=1}^N \delta^3\left({\mbox{\boldmath$x$}} - {\mbox{\boldmath$x$}}_j\right)$ (P.4.28)

となる.ここで通常,各粒子の質量$ m$ はすべて同じものとする.したがって,

$\displaystyle {\mbox{\boldmath$g$}}_i = \frac{G m}{a^3} \sum_{j (j\neq i)} \fra...
...}}_i}{\left\vert{\mbox{\boldmath$x$}}_j - {\mbox{\boldmath$x$}}_i\right\vert^3}$ (P.4.29)

となる.

ここで,2つの粒子の位置が近付きすぎると,その粒子ペア$ i$ , $ j$ に対する 式(16.5.31)の項が非常に大きくなる.この場合,その2つの粒子の間の 引力がその他の粒子全体からの寄与に比べて大きくなり,2体相互作用によって 軌道がお互いに大きく曲げられることになる.ところが,このような2体相互作 用は連続な物質分布では起こらないものである.近似的に連続分布を粒子で置 き換えたために表れてしまうのであり,本来あってほしくない相互作用である. そこで,宇宙構造形成の$ N$ 体シミュレーションでは,このような2体相互作用 を起こらないように重力加速度を修正するということが行われている.すなわ ち,粒子の重力加速度を

$\displaystyle {\mbox{\boldmath$g$}}_i = \frac{G m}{a^3} \sum_{j=1}^N \frac{{\mb...
...dmath$x$}}_j - {\mbox{\boldmath$x$}}_i\right\vert^2 + \epsilon^2 \right)^{3/2}}$ (P.4.30)

と修正する.こうするとある距離$ \epsilon$ 以下に近付いた粒子にはもはや大 きな加速度は働かない.このパラメータ$ \epsilon$ はソフトニング長と呼ばれ ていて,粒子を点粒子でなくこの半径ぐらいに広がったものと考えることに対 応する.

この重力加速度は粒子すべてについて求める必要があり,これをまともに実行 しようとすると,各時間ステップごとに $ {\cal O}(N^2)$ の項を計算することに なる.粒子数が増えるとこの計算にかかる時間は非常に大きくなり,実行する のが難しくなる.そこで,この各粒子の重力加速度をいかに早く計算するか が$ N$ 体シミュレーションでは非常に重要になってくる.

重力加速度をより速く求めるひとつの方法は,粒子ペアの和を直接取らずに空 間を細かくグリッドに分けてそのグリッド上で重力ポテンシャルを計算し,そ の微分から重力加速度を求めるやり方である.これを粒子・メッシュ 法(Particle-Mesh method; PM法)という.この方法では,まず粒子の空間の位 置から,グリッド上の密度を計算する.その次にポアソン方程式をフーリエ変 換によって数値的に解く.そうして得たグリッド上の重力ポテンシャルの値の 勾配から,各粒子の位置での重力加速度を内挿により求めるのである.数値的 なフーリエ変換は一見,フーリエ空間の各点について空間の各点を足し合わせ るためにグリッド数$ N_{\rm g}$ の2乗の演算を必要とするように見える.十分 な解像度を保つためにはグリッド数は粒子数と同じぐらいにしなければならな いので,直接粒子のペアを足し合わせるのと計算時間の面で変わりがないよう に見える.ところが,フーリエ変換の計算には,高速フーリエ変換法(FFT)とい う数値アルゴリズムがよく知られており,実際には  $ {\cal O}(N_{\rm
g}\log_2\! N_{\rm g})$ のオーダーの演算で済んでしまうのである.これは たとえば,空間を一辺あたり1000分割した $ N_{\rm g}=10^9$ の3次元グリッドに 対し, $ {\cal O}(N_{\rm g}^2)$ のアルゴリズムでは30年かかる計算も30秒で終 わってしまうことになる!

このようにPM法は高速に重力加速度を求められる単純で強力な手法ではあるが, グリッドサイズ以下の解像度がないというのが欠点である.このため,比較的 近くにある粒子から受ける力を正確に評価できない.そこでこれを改良すべく 考え出されたのが粒子・粒子/粒子・メッシュ 法(Particle-Particle/Particle-Mesh method; P$ {}^3$ M法)である.この方法 では,各粒子について,その粒子から遠くにある粒子から受ける力をPM法によっ て評価し,近くにある粒子から受ける力は式(16.4.30)の直接 の和により評価するという方法である.この方法により,FFTによる高速性を保 持つつ,近距離の力も正確に評価できるようになる.

フーリエ変換を用いずに高速化する手法として,ツリー法(Tree method)という ものもある.加速度を計算するとき,直接和(16.4.30)ではあ らゆる粒子を別々に考えていたが,ツリー法ではある程度遠くにある粒子群を まとめて一つの点で代表させてしまうという近似をする.これにはいろいろな やり方が考えられるが,よく行われているのは立方体のシミュレーションボッ クスを,一辺の長さがもとの半分になる8個の立方体に分割し,さらに同様にそ れら各立方体をさらに8個に分割し,ということを繰り返していく.この分割の 連鎖は立方体の中に粒子がひとつしか含まれなくなったところで打ち切る.し たがって場所によって分割の回数は異なる.こうして階層的な立方体ができ, 各立方体には一般に親と子の立方体がある.ただしもともとのシミュレーショ ンボックスの立方体は親を持たず,粒子をひとつしか含まない立方体は子を持 たない.すべての立方体どうしの親子関係を図にすると,もとの立方体を木の 幹にして,ちょうど木が何度も枝わかれしている様子に似ている.枝の先端に 粒子が一つずつ対応する.この準備のもと,個々の粒子に働く力を計算すると きに,ある粒子から見てある程度小さな立体角以下におさまって見える立方体 の中にあるすべての粒子は,その立方体中の粒子の重心にすべての粒子がある ものとして,ひとつひとつの粒子からの力を計算せずに一度の計算で済ませる. こうすることで考えている粒子の近くにある粒子についてはずっと小さな立方 体までツリーをたどってより精度よい位置で力をを計算し,遠くにある粒子は 比較的大きな立方体で位置を近似する.この方法の演算数は粒子数に対 し, $ {\cal O}(N\log N)$ のようにスケールすることが知られていて,直接全て の和をとる $ {\cal O}(N^2)$ の演算に比べて計算時間が$ \log N/N$ 倍と,非常に 速くなる.

これら以外にも重力加速度を計算するさまざまな方法が開発されているが,い ずれにしても計算時間を短縮するために近似を伴う.すべての粒子のペアの引 力をそのような近似なしにすべて足しあげる直接法は,通常のコンピュータで 行うと時間がかかりすぎるので非現実的である.だが、この重力計算をするこ とだけに特化した専用の計算機を用いることによって直接法による$ N$ 体計算 を現実的に行うことができる.実際に,このための専用チップを開発して搭載 したGRAPEというハードウェアが開発されて実用化されている.

実際の$ N$ 体シミュレーションによって宇宙の構造形成の様子を表した例を図 16.1に示した.

図 16.1: $ N$ 体シミュレーションによるコールドダークマター分布の時間進化
\includegraphics[height=23pc]{figs/fig00.eps}



Footnotes

... によって求めることができるP1
ここでは原理を説明しているだけで, 必ずしもこのような変数が使われるとは限らない. 例えば,この例は時間 間隔について1次精度しかない.実際には数値的な安定性などの観点から異 なる変数を使ったり,時間間隔を最適に選ぶなどさまざまな工夫が施され る.

next up previous contents index
次へ: 非線型構造の統計モデル 上へ: 摂動の非線形成長 前へ: 非線型成長の近似モデル   目次   索引

All rights reserved © T.Matsubara 2004-2010
visitors, pageviews since 2007.5.11