next up previous
Next: 3 数値計算の方法等 Up: 重力多体系 Previous: 1 はじめに

Subsections

2 少ない粒子数での計算

というわけで、問題は、少ない粒子数で計算するとはそもそもどういうことか ということだが、 これには実は少なくとも以下の3種類がある。第一は、銀河の 例のように、実質粒子数無限大の系を有限個の粒子で表すことである。第二は、 例えば100万体というふうに無限というほど粒子数が大きいわけではない系、 例えば球状星団を、10万とか1万といった、実際より少ない粒子数でモデル化 することである。最後に、例えば土星リングの一部だけをとり出してきて、周 期境界条件を入れて計算するような、系全体を扱うのではなく部分系を切り出 して扱う場合がある。

このそれぞれで、どのような問題が起きるかは非常に違う。なお、最初のと2 番目の違いがわからないかもしれないが、そのことを含めてこれから議論して いこう。

2.1 粒子数無限大の系

この例は例えば宇宙論での構造形成の計算や、空間スケールが違うだけともい えるが初期ゆらぎからの銀河形成の計算である。現在の標準的な理解では、宇 宙の全質量のほとんどは CDM (冷たい暗黒物質)が担っている。暗黒物質の正 体は依然明らかではないわけだが、とりあえず巨大なブラックホールがその辺 にゴロゴロしていると考えない限り2、一つの粒子としての質量は小さなものである。つまり、実効的な 粒子数は無限大と考えていい。 この場合、$N$体計算は 6 次元空間の中での分布関数を、離散近似しているこ とになる。

この場合、入ってくる誤差には2種類ある。一つは、重力不安定が線型ないし 弱い非線型段階で成長している時にはいってくる誤差であり、もうひとつは強 い非線型、というよりビリアライズした系で、粒子数が有限であるために入っ てくる誤差である。 ちゃんと初期条件を作れば(これはこれでいろいろ未解決 な問題があるが)、第一の誤差よりも第二の誤差が圧倒的に効くので、以下、第 二の誤差について述べる。

ここで、数値計算そのものは完全に正しいとすれば、入ってくる誤差は基本的には2体緩和である。従って、2体緩和がなにか理解できていればいい。簡単に復習しよう。

2.1.1 ２体緩和とはなにか？

まず、２体緩和とはいったいどういうものかというところから話を始め ることにする。原理的には、これがなにかというのは結構厄介な問題で ある。

有限粒子数の自己重力多体系を考えると、これは以下のような進化をす ると考えられる。まず、最初は力学平衡になかったとすると、とりあえ ず力学平衡に落ちつく。粒子数が無限大であれば、無限に細かく見れば 無限に時間がたっても真の力学平衡に到達するわけではないが、まあ、 漸近はしていく。この時、各粒子は与えられたポテンシャルの中を運動 するだけになり、それ以上進化することはなくなる。

さて、実際には有限粒子数であるので、そもそも真の力学平衡というも のはない。有限の質量をもった各粒子が系の中を運動するのに従って、ポ テンシャルは必ず変化するからである。この変化によって各粒子の軌道 も変化することになる。

それでは、粒子の軌道の変化を、粒子数が有限であることから来る成分 とそれ以外に分離することは可能であろうか？系が力学平衡に あるとみなすことができればそれは可能である。つまり、力学平衡にあれ ば、粒子のエネルギー変化は定義によりすべて粒子数が有限であること によるからである。

が、良く考えると問題なのは、そもそも有限粒子数であるものを力学平衡とみ なすとはどういうことかということである。このあたりを考えていると段々混 乱してくるので、以下、理想化された状況から順番に考えていくことにしよう。 なお、この節の内容は、 L. Spitzer の Dynamical Evolution of Globular Clusters[Spi87] にほぼ沿っている。2体緩和に関しては、重力多 体系での議論とプラズマ物理での議論は全くパラレルなものであり、上記の本 の内容自体、同じ著者の Physics of Fully Ionized Gases[Spi56] のものとほとんど同じである。


2.1.2 一様等方な分布

理想化といえば一様等方な分布を仮定することである。例えばマックス ウェル分布があって、その中の一つの粒子をとって考えるということを したいわけだが、これは結構厄介なのでさらに簡単な例を考える。すな わち、速度0で空間内に一様（ランダム）に分布した質点を考え、その中 を質量0のテスト粒子を飛ばして見る。

もちろん、この場合エネルギー交換はないので速度は変わらず、単に散 乱されるだけだが、しかし、この例は２体緩和のいくつかの重要な性質 を示すのですこし詳しく見ていくことにする。分布している質点の質量 を $m$、数密度を $n$ とする。図1のように、テスト粒子が一つの粒子から距離（イン パクトパラメータ） $b$ を速度 $v$で通った時に曲がる角度$\theta$は、実際に ケプラー問題の解析解を使って

$\displaystyle \tan \theta$ $\textstyle =$ $\displaystyle {2 b \over (b/b_0)^2-1}$  
$\displaystyle b_0$ $\textstyle =$ $\displaystyle {Gm \over v^2}$ (2)

で与えられる。単位時間当たり、インパクトパラメータが $(b,b+db)$の 範囲にある散乱の回数は $2 \pi nv b db $ である。

Figure 1: 2体散乱
\begin{figure}\begin{center} \leavevmode \epsfxsize =12cm \epsffile{2body.eps}\end{center}\end{figure}

さて、散乱の方向はランダムであると思われるので、平均としては（一 次の項は）0になる。しかし、 2次の項は0にならない。これは

\begin{displaymath} <\Delta \theta^2> = 2 \pi n v \int_0^{b_{max}} \delta \theta^2 bdb \end{displaymath} (3)

で与えられることになる。

この式から既にいろいろな性質がわかる。が、その前に理論的な困難を 解決しておく必要があるであろう。すなわち、この積分は $b\rightarrow \infty $ で発散しているのである。これについてはいく つかの考え方があった。例えば、初めて２体緩和の性質を理論的に調べ たチャンドラセカール[Cha43]は、以下のように考えた。

「平均粒子間距離よりもインパクトパラメータが大きいような散乱は、 多体の干渉によって効かなくなるのでそこで積分を打ち切ってよい」

しかし、多体の干渉というようなものが実際にあるかどうかはあきらか ではない。もっと素直な解釈は、実際に系にあるすべての粒子と常に同 時に相互作用しているのだから、システムサイズくらいまで全部いれる （系が構造を持つ場合はちょっとややこしいが、密度の空間依存も積分 のなかに入れて全空間で積分する）というものである。

数値実験の結果などから、後者の解釈すなわち全体が効くというほうが正しい ということはかなり昔から大体わかっていた。歴史的には、どちらの解釈が正 しいかについてはかなり最近まで論争があって、完全に決着がついたといえる のは 94-5年頃である。が、現在では後者の解釈が正しいということに疑いの 余地はない。

式(3) から、適当に近似すると

\begin{displaymath} <\Delta \theta^2> \sim Gn v^{-3} m^2 \log (R/r_0) \end{displaymath} (4)

となる。ここで $R$は先に述べたシステムの大きさ、 $r_0$ は「大きく 曲がる」ためのインパクトパラメータの値で、$b_0 = GM/v^2$の程度である。

さて、これからどんなことがわかるかというわけだが、これから、逆に 角度変化が1の程度になる時間というのを求めてみると、

\begin{displaymath} t_{\theta} \sim {v^3 \over Gnm^2 \log \Lambda} \end{displaymath} (5)

となる。ここで $\Lambda$は上の $R/r_0$を単に書き換えただけである。

今、$\log \Lambda$ の質量依存性といったものを無視すると、散乱のタ イムスケールは速度の3乗、数密度の逆数、質量の2乗の逆数に比例する ということがわかったことになる。特に、質量密度一定の場合というも のを考えてみると、タイムスケールが各粒子の質量に比例するというこ とがわかる。

ある大きさを持った多体系というものを考えてみよう。質量$M$、ビリアル 半径 $R$、粒子数 $N$ とすれば、ビリアル定 理から $<v^2>/2 = GM/R$、力学的なタイムスケール、つまり典型的な速度を持 つ粒子が系を横断するのにかかる時間が $t_d = R/v\sim \sqrt{R^3/GM}$ となる。これを使うと上の緩和のタイムスケールは

\begin{displaymath} t_{\theta} \sim {N \over \log N} t_d \end{displaymath} (6)

となる。つまり、力学的なタイムスケールに比べて、2体緩和のタイムスケー ルはほぼ $N$ 倍長いということになる。このために粒子数が大きいほど無衝 突系に近付くわけである。

2.1.3 バックグラウンドが速度分布を持つ場合

詳しく式を書くとそれだけでページリミットを超えるので、書くほうは楽でい いがそうもいかないので結果だけを書く。

速度分布を熱平衡、すなわち

\begin{displaymath} f_0(\mbox{\boldmath$v$}) = {n_f \over (2\pi \sigma^2)^{3/2}} \exp\left({- v^2/2\over \sigma^2}\right) \end{displaymath} (7)

とすると、単位時間当りの速度変化は
$\displaystyle <\Delta v_{平行}>$ $\textstyle =$ $\displaystyle -4 { n_f\Gamma\over \sigma^2}\left( 1 + {m \over m_f}\right)G(x)$ (8)
$\displaystyle <\Delta v_{平行}^2>$ $\textstyle =$ $\displaystyle 2\sqrt{2}{n_f\Gamma \over \sigma } G(x)/x$ (9)
$\displaystyle <\Delta v_{垂直}^2>$ $\textstyle =$ $\displaystyle 2\sqrt{2}{n_f \Gamma\over \sigma }{{\rm erf}(x) - G(x)\over x}$ (10)
$\displaystyle <\Delta E>$ $\textstyle =$ $\displaystyle \sqrt{2}{n_f \Gamma\over \sigma }\left[-{m \over m_f}{\rm erf}(x) + \left( 1 + {m \over m_f}\right)x {\rm erf}'(x)\right]$ (11)

ここで ${\rm erf}$は誤差関数であり、
\begin{displaymath} G(x) = {{\rm erf}(x) - x{\rm erf}'(x) \over 2x^2} \end{displaymath} (12)

また $x = v_t/(\sqrt{2}\sigma)$である。

ここで$\Gamma$は

\begin{displaymath} \Gamma = 4\pi G^2 m_f^2 \ln \Lambda \end{displaymath} (13)

である。ただし、 leading term でない項は適当に落ちてたりするので注意。

まず、速度の1次の項を見てみる。 $G$ の振舞いを考えると、$v$が大きい極 限では速度変化は速度の2乗に反比例す る。これに対して、$v$が小さい極限では、$f$を一定と見なすことが出来るの で、速度変化は $v$に比例する。

これは、タイムスケールを考えてみると、速度が大きい極限では減速のタイム スケールが $v^3$ であるのに対し、逆の極限では 一定になるということであ る。すなわち、非常に速度が大きい粒子が出来てしまうとこれはなかなか減速 しない。もちろん、自己重力系の場合には、そのようなものは系のなかに留ま るのは困難である。このような粒子（超熱的粒子）が問題になるのはプラズマ の場合である。

速度が小さいほうではタイムスケールがある一定値、つまりは $v\sim \sigma$ で決まる値あたりになる。

この1次の項は、 dynamical friction を表している。これ が問題になる場面は、例えば恒星系が質量の違う２つの成分から出来ているよ うな場合である。力学平衡状態で、分布関数に質量依存がないようなものを考える と、これは熱平衡から遠く離れている。従って、上の式で決まるタイムスケー ルで重いものがエネルギーを失い、軽いものがエネルギーを得る。これは、エ ネルギー等分配に向かう普通の熱力学的な進化である。

が、自己重力系ではこのエネルギー交換の結果熱平衡に近付くとは限らな い。つまり、重いものがエネルギーを失い、 軽いものがエネルギーを得るということは、それぞれの分布関数が変わり、空 間分布も変わるということである。具体的には、重いものは中心に集中した分 布になり、軽いものは外側に押し出される。その結果それぞれの成分の速度分 散がどうなるかの細かいところは初期条件に依存するが、普通は物が中心に集まれ ば重力も強くなり、力学平衡では結局速度分散も大きくなる。このために、熱 平衡からはかえって遠くなってしまう。

さて、次に、2次の項を見てみる。速度に平行な成分も垂直な成分も、 $v$が 大きい極限では0にいく。特に、 垂直な成分は$v$ に反比例する。これに対し、 速度が$0$の極限では、どちらも一定値に収束する。これはテスト粒子が停止 している極限でも、回りの粒子によって揺さぶられるということを表している わけである。

2.1.4 ２体緩和のタイムスケール

前節では、実際にバックグラウンドの粒子も動いている場合について、２体緩和 によって粒子の速度、エネルギーがどう変化するかの期待値（のモーメント） を求めた。ここでは、前回の結果を使って、まずいくつかの重要な場合について ２体緩和がどのように働くかを議論する。

今、フィールドに質量 $m_f$の粒子が一様に分布しており、テスト粒子として質量 $m_t$のものがこれもまた一様に分布しているとする。さらに、どちらも速度 分布はマックスウェルで与えられるとする。ここでは等分配を考えるので、そ れぞれの粒子１個当たりのエネルギーを $E_f$, $E_t$ と書く。今、テスト粒 子のエネルギー変化の平均を考えると、


\begin{displaymath} {d<E_t> \over dt} = 4\pi \int v_t^2f(v_t)<\Delta E_t>dv_t \end{displaymath} (14)

と書けることになる。これに前節で求めた $<\Delta E>$ を入れて実際に積分を 実行することができて、結果は


\begin{displaymath} {d<E_t> \over dt} = 2 \sqrt{\frac{6}{\pi}} {m_t n_f \Gamma \over m_f} {<E_f> - <E_t> \over (v_t^2 + v_f^2)^{3/2}} \end{displaymath} (15)

となる。

ここで、いくつかの極限的な場合を考えておくことは有益であろう。まず、 $m_t » m_f$で $v_t \sim v_f$ という状況を考えてみる。これはつまり非常 に重いものと軽いものが、同じような空間分布、速度分布で広がっている場合 である。この時は上の式で $E_t » E_f$ なので、

\begin{displaymath} {d\log <E_t> \over dt} = -\sqrt{\frac{3}{\pi}} {m_t n_f \Gamma \over m_f v^3} \end{displaymath} (16)

となる。なお、この時、変化率はテスト粒子の速度に分母の $v_t^2$ を通し てしか依存しないので、 $v_t\rightarrow 0$ の極限でエネルギー変化（減 速）のタイムスケールは一定値にいき、それは $v_t \sim v_f$ の時の値とそ れほど違わない。

次に、$m_t \sim m_f$で $v_t « v_f$という状況を考えてみる。この時は上 の式を

\begin{displaymath} {d(<E_t>/m_t) \over dt} = 2 \sqrt{6/\pi} { n_f \Gamma \over ... ...2 v_f^{3}} =8 \sqrt{6\pi}G^2\ln \Lambda n_f m_f <E_f> v_f^{-3} \end{displaymath} (17)

となる。 ここで
\begin{displaymath} \Gamma = 4\pi G^2 m_f^2 \ln \Lambda \end{displaymath} (18)

を使った。 さらに、 $E_t = m_tv_t^2/2$ などを使って書き直せば
\begin{displaymath} {d(v_t^2) \over dt} =4 \sqrt{6\pi}G^2\ln \Lambda n_f m_f^2 v_f^{-1} \end{displaymath} (19)

を得る。つまり、速度が小さい極限では、一定の率でエネルギーをもらうわけ である。言いかえれば、温度が2倍になるタイムスケールというものは、温度 に比例して小さくなるともいえる。

さて、通常「2体緩和のタイムスケール」という時には何を指してい るかというと、この等分配のタイムスケールのことではないのが普通である。 が、時と場合によっていろんなものが出てくるが、まあ同じようなものである。 普通に使われるのは、

\begin{displaymath} t_r = {1 \over 3} {v_m^2 \over <v_{平行}^2>_{v = v_m}} = {v... ...\over 1.22 n\Gamma} = {0.065 v_m^3 \over nm^2G^2 \log \Lambda} \end{displaymath} (20)

とするものである。 ここで $v_m$はr.m.s. 速度である。$1/3$になにか意味があるわけではなく、 こう定義したというだけである。

これは、ローカルな量で定義されていて、例えば系全体の緩和時間といったも のを考えるのにはちょっと不便なこともある。というわけで、いわゆる half-mass relaxation time $t_{rh}$ というものを導入しておく。これは、 半径$r_h$ の中に質量の 1/2があるとして、その中の密度は一様であるとし、 さらにビリアル定理から出てくる $T \sim 0.2 GM^2/r_h$といった関係 を使って出すことが出来る。ここでビリアル半径ではなく half-mass radius を使ったのは慣習に従っただけで深い意味はない。この2つは通常 30% 程度 の範囲で一致する。 これは

\begin{displaymath} t_{rh} = 0.138 {Nr_h^{3/2} \over M^{1/2} G^{1/2}\log \Lambda} \end{displaymath} (21)

となる。

ここで注意しないといけないことは、 $t_{rh}$はあくまでも球対称に近 い系の half mass radius のあたりでの緩和時間であるに過ぎないとい うことである。従って、球状星団全体の緩和時間とか、あるいは楕円銀 河、銀河団といったものには有効な概念であるが、球対称から大きくず れた銀河とか、あるいは half mass radiusのずっと外側やずっと内側で は全く違ったものになっていることに注意する必要がある。さらに、速 度分布が非等方であるとか、回転がメインであるとかでもまた話が全く 変わってくる。このような場合、ローカルな緩和時間、あるいはエネル ギー変化自体の式に戻って考えないと、タイムスケールについて全く間 違った推定をしてしまうことになる。

2.1.5 実例

2体緩和が比較的新しい研究でも問題を起こしていた例として、有名な NFW プ ロファイル[NFW96]がある。彼らは、ダークマターだけの構造形成計算 を、様々な宇宙モデルに対して行い、ダークマターハローの構造が宇宙モデル に依存するか、また依存するとすればどのように依存するかを調べた。その結 果が良く知られている NFW プロファイルであり、半径の小さい極限で $\rho \propto 1/r$、大きい極限で $\rho \propto 1/r^3$ となり、その間を滑らか につないだ形をしている。つまり、

\begin{displaymath} \rho \propto {1 \over r_*(1 + r_*)^2} \end{displaymath} (22)

ここで $r_*$ は適当にスケールした半径である。 しかも、この形はユニバーサルで、例えば一つの銀河でも銀河団でも同じ であるというのが彼らの結果であった。

ここでの問題は、この結果を信じていいかどうかである。彼らの使った粒子数 は $10^4$ 程度であり、中心部での二体緩和時間が宇宙年齢よりはるかに短い。 従って、中心部が信用できるかどうかはわからない。というわけで、 我々はとにかく粒子数を増やした計算をしてみた[FM97]。 これは約 80万粒子で、GRAPE-4 を使った直接計算であって計算精度も(おそらく必要 以上に)高いものである。

結果は、NFW がいっているような中心で $1/r$ になる分布にはならなず、む しろ $\rho \sim r^{-1.4}$程度になった。また、銀河サイズのハローでは、 信用できるのは半径 1-2kpc より外側だけであった。

結局、 NFW の結果が ``universal''になったのは数値計算の誤差のためであっ た。すべての計算で粒子数が同じ程度、ソフトニングパラメータも同じ程度だっ たので、結果がその2つで決まっていたのである。

振り返って見るとあまりに明白な誤りではあるが、しかし実際にこういう間違 いをした人はいるわけである。

2.2 有限の粒子数の系の場合

ここまで述べてきたことから、実際の系の粒子数が有限であるということの意 味は明白であろう。つまり、有限ということは、2 体緩和の 効果が現実の系の進化に影響する、あるいは進化を駆動する原動力になってい るということである。

このような系の典型は球状星団である。典型的な 2体緩和のタイムスケールが $10^9$ 年程度になり、系の進化は基本的には2体緩和によって進む。従って、 原理的にはタイムスケールを読み変えて、2体緩和のタイムスケールを合わせ てやれば、粒子数が大きい系も小さい系も「同じ」振舞いをするはずである。 細かいことをいいだすと、 $\log \Lambda$ の半径依存性が微妙に変わったり もするが、まあ普通は問題ではない。

もしも、2体緩和だけで話が済むなら、上のスケーリングで話は終わっていて、 別に難しく考えることはない。もちろん、話はそれではすまないからこう書く わけである。例えば現実の球状星団の進化であれば、2体緩和の他に少なくと も以下の2つは系の進化を大きく左右する。一つは構成要素である個々の恒星 の進化である。主系列から離れた星は巨星になって質量放出したり、さらには 超新星爆発を起こして質量を失ったりする。これにより星団全体の重力エネル ギーが小さくなり、星団は2体緩和がなくても膨張することになる。もう一つ は、親銀河の潮汐場である。大雑把にいって親銀河のポテンシャルと自分のポ テンシャルで決まるロシュ半径を超えて外側にいったものは星団から脱出して しまう。星団の質量が減るとロシュ半径も小さくなるので、 ますます星団か ら星は脱出しやすくなる。

今、恒星進化はちょっと置いて潮汐場の効果と2体緩和の関係を考えてみよう。 もしも2体緩和がなく、潮汐場も時間的に定常である(星団の銀河内の軌道が円 軌道である)なら、 一旦力学平衡に達した後は星団から星が逃げることはない。 つまり、星団から星が逃げるのは2体緩和の効果である。これは Spitzer の教 科書にものっているとおりで、理論的におかしいことはない。

しかし、問題は、それでは星団からの星の蒸発率は2体緩和時間だけで決まる かどうかである。1990 年頃まではそう考えられていた。この考え方の集大成 ともいえるのが、 Chernoff and Weinberg[CW90] によるフォッカープランク方程式 によるモデルに近似的に潮汐場の効果を入れ、さらに恒星進化の効果もとりい れたモデルである。 このモデルにより彼らは非常に広い範囲の初期モデルか らの球状星団の進化のサーベイを行い、 現在まで生き残る範囲等を求めた。

フォッカープランク方程式を解くのは $N$ 体計算ではなく、分布関数に対す る拡散方程式を時間積分するわけだが、ここでの基本的な仮定は緩和時間が力 学的なタイムスケールに比べて十分長い、つまり粒子数が十分に大きいという ことである。他にも速度分散が等方的とか分布が球対称とかいろんな仮定はあ るが、とりあえずこれらはそれほど悪さをしない。

力学的なタイムスケールと星団の寿命の関係に初めて着目したのは Fukushige and Heggie[FH95] である。彼らは、 Chernoff and Weinberg のモデルのい くつかについて同様な計算を、しかしフォッカープランク方程式を解くのでは なく GRAPE-3 を使った $N$ 体計算で行なった。$N$体計算の粒子数は実際の 球状星団の粒子数に比べて小さかったので、彼らは力学的なタイムスケールを 恒星進化のタイムスケールに合わせた。これでも、星団の寿命が非常に短い場 合には正しい結果が得られると考えられる。しかし、彼らの検討したいくつか の場合について、得られた寿命は Chernoff and Weinberg の結果の 10倍にも 及ぶものであった。つまり、力学的なタイムスケールが星団の寿命に大きな影 響を持つのである。

この状況をより系統的に調べたのは Portegies Zwart et al.[PHMM98] である。彼らは、 広い範囲で粒子数、つまり2体緩和の時間スケールと力学的タイムスケールの比 を変化させ、星団の寿命がどのように変わるかを調べた。恒星進化のタイムス ケールには2体緩和を合わせた。この結果、やはり寿命は粒子数に依存し、寿 命が宇宙年齢程度の場合でも、現実の球状星団の寿命は粒子数無限大の極限で あるフォッカープランク方程式を解いてでてくるものとはかなり違うことがわ かったのである。

このようなずれが生じる理由は、単にエネルギー的は系に束縛されなくなった 粒子でも、実際に系から離れていくには有限の時間がかかり、その時間は力学 的なタイムスケールであるということである。つまり、定性的には極めて当然 のことなのであるが、このような効果は最近まで考慮されていなかったのであ る。

潮汐場と2体緩和のカップリングによる進化については、最近になっても Baumgardt [Bau01]による新しい結果がでている。

こういった話になると$N$体計算の問題とはいいがたいかもしれないが、しか し、数値シミュレーションと現実の系の違いという意味ではこのような問題は 単純な数値計算の問題以上に重要であるといっても間違いではないであろう。

next up previous
Next: 3 数値計算の方法等 Up: 重力多体系 Previous: 1 はじめに
Jun Makino
平成15年3月19日