自己重力多体系の物理

牧野淳一郎

前回は、Violent relaxation の結果として自己重力系の外側(virial radius のずっと外側)については、

さて、それでは、中心部の構造についてはなにもいえないのであろうか？これは実はまだ良くわかっていない問題である。10年前に、 Navarro たち（ApJ 1997, 490, 493) は、数値計算の結果をもとに以下のような主張をした

彼らはなぜそのようなことが起きるかについての解釈とか説明は特に与えていないが、例えば Syer and White (MN 1998, 293, 337)といった人達が説明を考えてはいる。

しかし、実は、Navarro たちの結果の解釈は割合すぐに異論がでており、 CDM と初期条件を制限しても、例えば Fukushige and Makino (ApJ 1997, 477 L9) とか Moore et al. (ApJ 1998, 499, 5L) を見ると、上の「ユニバーサル」な形になったのは数値誤差のせいという主張がなされている。これらの結果では Navarro たちのものより中心で等温に近くなっている。 Navarro らの結果は 1 万粒子程度であるが、 Fukushige ら、 Moore らは 100 万粒子程度であり、数値誤差の影響が小さくなっていることは間違いない。「真の」のスロープがどうなるかについては現在も活発な研究が続いているが、現時点では、CDM からの数値実験でできるハローの中心部は NFW よりもスロープの傾きが大きいということはほぼ万人の認めるところになったようである。

全く余談であるが、福重君はこの仕事とその続きで天文学会研究奨励賞をもらった。

さらにもっと余談であるが、福重君がこの計算をした動機は、1996年に GRAPE-4 を使って何か大きな計算で高い実効性能を出したい(Gordon Bell Prize に応募したい)、そのためにはとにかく粒子数が沢山必要で、普通の treecode とかではできないようなシミュレーションでよいネタがないか？とか色々考えたことである。

つまり、現状では、中心部（というか、 half mass radius より内側）の構造は、

という状況であるといえる。現在では、福重君の仕事からさらに 3 桁程度粒子数を増やした、

粒子程度の計算も行われているが、中心部の構造がどうなるか、それは何故か、というのは依然明確にはなっていない。

1.2 シミュレーションの結果

図は、現時点で最大粒子数、最高分解能の計算: Springel et al (2008) である。これは粒子数を3桁程度の範囲で変えて、5個のシミュレーションをしており、分解能が高い計算は低い計算の、計算が信用できそうな範囲は再現している。そういう意味では、計算結果は収束していて、求まっている答はそれなりに正しいと考えられる。

粒子数が少ないと中心のほうが信用できなくなるのは、熱力学的な緩和のため。これはまた後でもう少し詳しく。

こちらは密度勾配をグラフにしたもの。計算の範囲では中心にいくに従って浅くなっていて、いかなるべきでもない。

これは Navarro et al. 2008. 傾向は同じで、中心にいくに従って密度勾配が浅くなっている。

というような感じになる。数値計算の結果を疑うわけではないが、何故そうなっているかは今のところ全く理解されていない。感覚的には

つまり、数値計算の結果はあっても、それを理解できる理屈がついてない、というのが現状での大きな問題点である。

数値計算の結果が理論的になにかしら理解できない時には、それが「正しい」かどうかは判断のしようがない。

もっとも、こんなところの計算にはそもそも意味があるのか？という疑問はある。銀河サイズのダークマターハローの中心 1 kpc 以内の構造に、観測的になんか意味はあるのか？ということである。もちろん、これには2つの回答がありえる。一つは「ない」というもので、その理由は少なくとも実際にの銀河ではその辺はどうせバリオンの影響のほうが大きいからである。もうひとつは、「ある」というもので、もっと小さい、バリオンをトラップしてない本当にダークなハローならもちろんこのようなものは残っているはずである。

といった話があり、こういうことを考える時には小さい質量のダークマターハローの構造は重要である。原理的には、小さなものが生き残っていて太陽近傍にある、といった可能性も考慮する必要がある。

1.3 原理的な問題

小さなハローの構造を考える、という時には明らかに問題なことは、宇宙論的

体シミュレーションは、「無衝突系」の正しいシミュレーションになっていないことである。

宇宙論的 N 体シミュレーションでは、以下のようにして初期条件を作る。

というような感じ。実際のパワーは何十桁か先まで続いているが、粒子数の限界があるのでどこかで打ち切る。それよりちょっと低い波数でゆらぎも打ち切る。

で、コールドダークマターだと、小さいものが先にできる。このため、シミュレーションで最初にできる「ハロー」は粒子10個程度と非常に小さい。ここで、 (私を含む)業界の信念は、階層的な合体成長のプロセスが構造を決めているので、最小の構造がちゃんと計算できているかどうかは結果に影響しないはず、というものだが、これは検証されたことはない。また、ダークマター検出に関係するような最小スケールのハローではこの信念は明らかに間違っている。

1.4 最小ハロー

最小ハローがどれくらいの大きさかは、大体地球質量くらいと考えられている。これは、ダークマターがいつ他のものからデカップルして、自由に動くようになるかできまるが、現在好まれているダークマター粒子のモデルでは地球質量くらい。サイズは 100 AU くらいになる。そうすると、

1.5 これまでの研究

といってもほとんどこれだけ: Diemand et al. 2005, Nature 433, 389。基本的に普通の宇宙論的シミュレーションで、地球スケールのハローを

粒子くらいで分解している。密度構造を議論するには粒子数不足。

通常の宇宙論的シミュレーションでは、実際にはずっと先にあるカットオフとは別に、ずっと質量が大きいところで質量分解能の限界がくる。そのちょっと手前でカットオフを入れる。

Ishiyama et al. の計算では、物理的なカットオフを入れる。それをちゃんと表現するため、粒子の質量はそこに十分沢山の粒子が入る程度に小さくする。

実線がカットオフあり、波線がカットオフなしのもの。カットオフがあるほうが中心密度が高くなる。ほぼ

のべきで直線でのびているように見える。カットオフなしだと銀河スケールと同様に、中心にむかって浅くなる。

つまり、最小サイズのダークハローは $\rho \propto r^{-1.5}$ と書けるような深いべきの密度構造を持ち、中心密度が高いと考えられる。

これは、特に対消滅によるガンマ線輻射には大きな影響がある。ガンマ線輻射率が密度の二乗に比例するので、フラックスが中心にむけて対数的に発散するからである。

が問題である。とはいえ、なぜ

かはまだあまり理解が進んでいない。密度の上限については、いわゆる ``pha space argument'' を考えてみる。これは、無衝突ボルツマン方程式が6次元位相空間での非圧縮流れを与える、ということを使う議論である。

色々な仮定すると、初期の位相空間密度が $\sim 10^{15} {\rm M_{\odot} pc^{-3} (km/s)^{-3}}$ というのがでてくる。これくらいになったところで密度一定のコアができると思うと、

このようなハローが現在まで生き残っているかどうかを調べるため、主要な破壊メカニズムと考えられる恒星との衝突をシミュレーションしてみた。

非常に近くを通った場合でも中心の高密度領域は簡単には壊れない。完全に壊すためにはインパクトパラメータが $b=5\times 10^{-5}{\rm pc}$ 以下でないといけないとわかった。

これは宇宙年齢程度では普通起きない(太陽系もこんなのがあると破壊される)

現在、「京」を使って、この最小ハローからもっと大きなハローがどのようにできるか、といった辺りを計算しようとしている。

2 2体緩和

2.1 前置き

のうちの最初の３つまでが済んだということになる。もっとも、それぞれ大きなテーマであり、力学平衡モデルについては軸対称、あるいは３軸不等なモデル、あるいはディスク系とその不安定モードなどについては全く扱えていない。 (これは後半でなんとかしたい)力学平衡への緩和過程についても、例えば合体の場合、あるいは理想化された１次元系を例にとって具体的な話をする必要もあったかもしれない。が、しかし、世の中は無衝突系だけではないし、また、現実は無衝突系である場合でも、それを数値的にモデル化したものはそうでなくなってしまうことがある。従って、衝突系の進化がどのようなものかということは、実際になんらかの自己重力系を扱う場合には必ず理解しておく必要がある。というわけで、これからしばらくで衝突系の進化というものを考えることにする。

2.2 ２体緩和とはなにか？

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

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

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

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

が、良く考えると問題なのは、そもそも有限粒子数であるものを力学平衡とみなすとはどういうことかということである。このあたりを考えていると段々混乱してくるので、まず、理想化された状況から考えていくことにしよう。

2.3 理想化：一様等方な分布

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

もちろん、この場合エネルギー交換はないので速度は変わらず、単に散乱されるだけだが、しかし、この例は２体緩和のいくつかの重要な性質を示すのですこし詳しく見ていくことにする。分布している質点の質量を

、数密度を

とする。テスト粒子が一つの粒子から距離（インパクトパラメータ）

を速度

で通った時に曲がる角度は、実際にケプラー問題の解析解を使って

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

さて、散乱の方向はランダムであると思われるので、平均としては（一次の項は）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$ で発散しているのである。これについてはいくつかの考え方があった。例えば、初めて２体緩和の性質を理論的に調べた Chandrasekhar は、以下のように考えた。

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

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

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

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

(4)

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

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

(5)

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

ある大きさを持った多体系というものを考えてみよう。質量

、特徴的な半径（ビリアル半径か何か）

、粒子数

とすれば、ビリアル定理から

、力学的なタイムスケールが $t_d \sim \sqrt{R^3/GM}$ となる。これを使うと上の緩和のタイムスケールは

$\begin{displaymath} t_{\theta} \sim {N \over \log N} t_d \end{displaymath}$

(6)

2.4 どういうことを考えるかということ：流体との違い

２体緩和によって、最終的には系が熱力学的に進化するわけであるが、これが普通の流体（ガス）とは本質的に違うものであるということをここで再確認しておこう。

ガスの場合、粒子の平均自由行程はシステムサイズよりもはるかに小さい。液体であれば平均粒子間距離は粒子のサイズ程度であるし、気体であっても通常の状況では考えている現象の空間スケールに比べて平均自由行程は小さい。ちなみに、非常に希薄な気体とか、あるいは本当に空間スケールの小さい現象では平均自由行程が問題になる。これは例えば超高層での人工衛星の回りの気体の流れとか、あるいは最近の磁気ヘッドの回りの空気の流れとかいったものである。

とにかく、通常のガスの場合、平均自由行程がシステムサイズより小さく、システムサイズよりは小さく平均自由行程よりは大きいような空間スケールを考えると、そのなかでほぼ熱平衡になっていると思っていいことになる。いいかえれば、いわゆる Local thermal equilibrium (LTE) の仮定が使える。こうなると、温度とか圧力とかいった量が近似的（といっても実際上非常に高い近似精度で）に定義でき、そういったマクロな量で系の進化を扱う、特に熱の流れを拡散方程式で書くということが可能になる。

しかし、自己重力質点系では状況が全くことなる。まず、粒子数が無限大の極限では、平均自由行程も無限大であった。つまり、 LTE がなりたたないどころか、そもそも熱平衡に向かう（すなわちエントロピーを生成する）ようなメカニズムがなかったわけである。

粒子数が有限の場合も、依然として平均自由行程が長い、つまり、粒子数無限大の時の軌道から、他の粒子との相互作用によって段々ずれていくわけだが、そのずれる典型的なタイムスケールは $N t_d/\log N$ 程度であった。つまり、流体の場合とは全く逆に、ほとんど自由運動（というか、他の粒子全体が作るポテンシャルに沿った運動）をしていて、その場が有限の粒子で表現されるための揺らぎがあるので段々軌道が変わっていくということになるわけである。

従って、ローカルな熱平衡を仮定して拡散係数／輸送係数を求めるというのとは逆に、ある一つの粒子が系の中を動き回りながらどういうふうにエネルギー等を変化させていくかという観点で見ていくことになる。

これをすこし別ないい方をすれば、通常の空間のなかでの密度や温度の変化を考える代わりに、また6次元位相空間のなかでの分布関数の進化を考えるということに当たる。具体的には、これまで無視してきた「衝突項」というものをちゃんと評価して、どういうものかみてやろうということである。

2.5 バックグラウンドの分布のもとでの有限質量のテスト粒子の振舞い

さて、以下ではバックグラウンドの粒子分布のもとでの一つのテスト粒子の振舞いを考える。前と違うのは、バックグラウンドも動いていることと、テスト粒子も有限の質量を持つことである。バックグラウンドの粒子は一様に分布するものとし、ある速度分布に従うとする。さらに、バックグラウンドの粒子間の相互作用とかは考えないことにする。これで本当にいいかどうかはちょっと良くわからない問題であるが、まあ、とりあえずやってみることにしよう。

前と同じく、分布している質点の質量を

、数密度を

とする。テスト粒子が一つの粒子から距離（インパクトパラメータ）

を相対速度 $\mbox{\boldmath$V$}= \mbox{\boldmath$v$}_t - \mbox{\boldmath$v$}_f$ で通った時に曲がる角度は、実際にケプラー問題の解析解を使って

$\displaystyle \tan \theta$	$\textstyle =$	$\displaystyle {2 p/p_0 \over (p/p_0)^2-1}$
$\displaystyle p_0$	$\textstyle =$	$\displaystyle {G(m_t+m_f) \over V^2}$	(7)

$\displaystyle \Delta v_{垂直}$	$\textstyle =$	$\displaystyle {m_f \over m_t + m_f} V \sin \theta = 2V {m_f \over m_t + m_f} {p/p_0 \over 1 + (p/p_0)^2}$	(8)
$\displaystyle \Delta v_{平行}$	$\textstyle =$	$\displaystyle {m_f \over m_t + m_f} V (1-\cos \theta) = -2V {m_f \over m_t + m_f} {1 \over 1 + (p/p_0)^2}$	(9)

$\displaystyle <\Delta v_{垂直}^2>$	$\textstyle =$	$\displaystyle {2n_f\Gamma \over V}$	(10)
$\displaystyle <\Delta v_{平行}>$	$\textstyle =$	$\displaystyle - \left(1 + {m_t \over m_f}\right){n_f \Gamma \over V^2}$	(11)
$\displaystyle <\Delta v_{平行}^2>$	$\textstyle =$	$\displaystyle {n_f\Gamma \over V\ln \Lambda}$	(12)

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

(13)

上の式で、 $<\Delta v_{垂直}^2>$ の項は前にでてきた角度の曲がる項と同じものである。前の話と違うのは、ネットに速度が小さくなる成分がある、すなわち $<\Delta v_{平行}>$ が負で有限の値をもつということである。

これは、実は前にやった dynamical friction そのものである。つまり、回りが止まっているなかを粒子が走っていくと、それが回りを引っ張って動かすので、その分エネルギーを失って段々速度が落ちるわけである。これは、

が大きい（

が小さい）極限では

、つまり質量密度によっていて、バックグラウンドの粒子の質量に依存しないことに注意してほしい。これに対し、他の項は

に比例していて、質量密度が同じでも粒子の質量が大きいほうが値が大きくなるのは先週にやった通りである。

さて、ここではとりあえず1次と2次の項を求めたわけだが、それより先の項については考えなくてもいいのだろうか？ここでは粒子の軌道変化がたくさんの散乱のランダムな重ね合わせで書けるとした。この仮定が正しければ、たくさん散乱を受けた後の速度の分布は1次と2次のモーメントで決まるガウス分布になり、従って3次より高いモーメントの寄与は考えなくてもいいことになる。

問題はこの仮定が正しいかどうかであるが、実は理論的にはそれほど正確なわけではない。というのは、インパクトパラメータが例えば

の程度の散乱というのも現実におき、その効果はそれ以外の散乱すべての寄与に比べてせいぜい $\ln \Lambda$ 程度でしか小さくないからである。まあ、しかし、そんなことをいっていても高次の項があっては計算出来ないし、とりあえず $\ln \Lambda$ 程度で小さいということも確かなので、以下高次のモーメントは考えない。

この文書について...

コマンド行は以下の通りでした。:
latex2html -nomath_parsing -local_icons -show_section_numbers -split 0 note6-e.tex.

翻訳は Jun Makino によって平成24年9月30日に実行されました。