重力多体系の進化

牧野淳一郎

この講義では、天体物理学の扱う主要な対象の一つである自己重力多体系の進化とその理解について、その概要を扱うことにしたい。

「自己重力多体系」というと、難しそうに聞こえるが、基本的には「たくさんの質点が重力で引き合って、まとまった天体を作っている」というものである。例えば方程式で書けばすべての粒子が他のすべての粒子からの重力を受けて運動するというだけのことである。

しかし、実際に宇宙にある自己重力多体系を見ると、非常に多様な進化を遂げているということはすぐにわかる。この多様性がもっとも顕著なのは銀河であり、基本的には重力で星が集まっているだけなのに、円盤（渦巻）銀河、棒渦巻銀河、楕円銀河、あるいは不規則銀河などじつに多様な形をしたものがある。また、最近の特に HST による High-Z の観測から、このような銀河の形は進化してきたものである、例えば昔に遡って見ると「晩期型」銀河のほうが「早期型」銀河より多い（方向が逆だけど、これは定義の問題）ということがわかってきた。

これに対し、例えば「球状星団」と呼ばれる銀河内の星団は、比較的単純な丸い形をしていて、回転もほとんどしていない。また、半径方向の光度分布についても共通の特徴がある。

さらに小さいスケールで、惑星系を考えると、20年くらい前までは惑星系というと我々の太陽系しか知られていなかったので、多様かどうかも良くわかっていなかったが、系外惑星が見つかってくると極めて多様であることがわかってきた。

ここで銀河と惑星系には、「どちらもほぼ円盤である」という共通性があり、これはもちろん、どちらもガス雲が収縮していく過程でできたもので、初めから質点系であったわけではないからである。つまり、自己重力的なガスにはまたそれなりの共通な性質がある。

逆に大きいスケールに移って、銀河群、銀河団といったものを考えると、これらは（銀河の数が少ないこともあるが）良くわからない形をしているものが多い。結構丸くまとまっているものもあるが、縦に伸びたもの、二つの銀河団がお互いの回りを回っているものなど、多様な銀河団がある。また、比較的大きな銀河団のなかには、 cD と呼ばれる異常に大きな楕円（ほとんど球に近い）銀河を持つものが多い。

さらに、ここ 30 年ほどの、大規模な銀河分布についての研究から、宇宙には銀河、あるいは銀河団が一様に分布しているわけでは全くなく、大きなスケールで構造があるということがわかっている。

以上のようなさまざまな現象は、「重力多体系」という描像で（基本的には）統一的に理解することができる。その統一的な描像の概要を与えるのがこの講義の目標ということになる。

2 本題

まず、無衝突ボルツマン方程式を導いてその振舞いについて少し考える。無衝突ボルツマン方程式は、恒星系力学の基礎方程式といえる。

原理的には、恒星系力学の基礎方程式は何かというと、恒星系は自己重力多体系と考えていいわけで、各粒子(恒星)の運動方程式

$\begin{displaymath} {d^2 {\bf x}_i \over dt^2} = \sum_{j\ne i} Gm_j {{\bf x}_j - {\bf x}_i \over \vert{\bf x}_j - {\bf x}_i\vert^3}, \end{displaymath}$

(1)

が基礎方程式ということになる。数値計算にはもちろんこれを使うわけだが、理論的な扱いには不便である。というわけで、しばらくは(1粒子)分布関数 $f({\bf x},{\bf v},t)$ で話をする。ここでは、粒子数が「無限に大きい」と思って、6次元の位相空間 $({\bf x},{\bf v})$ の中の粒子の密度分布の時間進化を考える。この時の基礎方程式が(無衝突)ボルツマン方程式である。

以下、方程式を導く。位相空間での座標を $\mathbf{w}= ({\bf x},{\bf v},t)$ と書くことにする。また、重力ポテンシャルを $\Phi=\Phi({\bf x},t)$ とおくと、位相空間の中での粒子の流れは

$\begin{displaymath} \mathbf{\dot w}= (\mathbf{\dot x}, \mathbf{\dot v}) = ({\bf v}, -\nabla \Phi) \end{displaymath}$

(2)

$\begin{displaymath} {\partial f \over \partial t} + \sum_{i=1}^6 {\partial f \dot w_i \over \partial w_i} = 0, \end{displaymath}$

(3)

$\begin{displaymath} \sum_{i=1}^6 {\partial \dot w_i \over \partial w_i} = \su... ...rtial v_i} \left(\frac{\partial \Phi}{\partial x_i} \right)= 0 \end{displaymath}$

(4)

$\begin{displaymath} {\partial f \over \partial t} + {\bf v}\cdot \nabla f - \nabla \Phi \cdot {\partial f \over \partial {\bf v}} = 0, \end{displaymath}$

(5)

$\begin{displaymath} \nabla ^2 \phi = 4\pi G \rho. \end{displaymath}$

(6)

$\begin{displaymath} \rho = m\int d{\bf v}f, \end{displaymath}$

(7)

なお、以下の議論では（当分）

のことは忘れて、その代わり

が個数密度ではなくて質量分布であるということにしておく。

$\begin{displaymath} {\partial f \over \partial t} + {\bf v}\cdot \nabla f - \nabla \Phi \cdot {\partial f \over \partial {\bf v}}\end{displaymath}$

3 力学平衡とジーンズの定理

これらから、まず、「力学平衡状態」とはどう定義され、どういう性質があるかということを考え、それから具体的な平衡状態の例を見ていくことにする。

まず、「力学平衡」とは何かということだが、これは、上の無衝突ボルツマン方程式とポアソン方程式を連立させたものの定常解、すなわち、時間的に変化しない解ということになる。従って、ある分布関数

が力学平衡にあるということは、それによって決まるポテンシャル $\Phi$ を固定して考えた時に、

の時間微分が

になるということである。

直観的には、粒子数が無限に多い重力多体系で、時間がたっても形が変わらない、というのが力学平衡にある、ということである。

これは楕円銀河とかは実際にそういう状態にあるが、どうやってそういう状態が可能になっているかはちゃんと理解するのは意外に面倒である。以下、簡単な場合から考えていく。

ここで、わざわざ「力学」とつけるのは、もちろん平衡状態にはほかにもいろいろあるからである。もっとも重要なのは熱平衡の概念であるが、これはまた後で。

あと、平衡状態だけ考えでもしょうがなくて、私達は系がどう進化するかを理解したかったのでは？という疑問は当然あると思うが、これは

4 運動の積分

平衡状態というものを考える上で基本になるのは、「運動の積分」という概念である。ポテンシャル $\Phi$ のもとで、ある ${\bf x}, {\bf v}$ の関数

が運動の積分であるとは、その上で

$\begin{displaymath} {d \over dt} I({\bf x}, {\bf v}) = 0, \end{displaymath}$

(8)

$\begin{displaymath} {\bf v}\cdot \nabla I - \nabla \Phi \cdot {\partial I \over \partial {\bf v}} = 0 \end{displaymath}$

(9)

なお、「運動の積分」というときの流儀は２通りあって、一般に運動の保存量のことを「運動の積分」という流儀もあるが、ここでは位相空間の座標だけの関数であって同時に保存量であるものをさす。具体的には、たとえば１次元調和振動子で「初期の位相」というのは保存量だが運動の積分ではない。これは、時間が入ってくるからである。

4.1 例

エネルギー $v^2/2 + \Phi$ や、ポテンシャルが球対称（

だけの関数）の場合の角運動量ベクトル ${\bf L} = {\bf r}\times {\bf v}$ は運動の積分である。

5 ジーンズの定理

さて、上のように

を定義すると、以下の「ジーンズの定理」がなり立つことがわかる。

ジーンズの定理 任意の無衝突ボルツマン方程式の定常解は、運動の積分を通してのみ位相空間座標に依存する。逆に、任意の運動の積分の関数は定常解を与える。

いいかえると、分布関数

が定常であるためには、運動の積分

があって

の形で書けることが必要十分ということ。

証明だが、まず「定常ならば運動の積分で書ける」というほうを考えてみる。これは、

自体が運動の積分の定義を満たしているので、 OK。

逆のほうは、実際に

の全微分を

で書き下せば、それぞれが 0 になるということからいえる。

というわけで、これはなかなか強力な定理だが、一般の場合にはそれほど役に立つわけではない。というのは、ポテンシャルを与えた時に一般に運動の積分というのは 5 個あるはずだが、それらをすべて知っているということはないからである。

ただし、球対称とか軸対称とか条件をつけると、いろいろちゃんと決まるようになる。以下、まず球対称の場合を考える。

6 球対称の場合

球対称の場合、運動の積分はエネルギーと角運動量の3成分で4つある。一般にはもう一つあるが、これは特別な場合を除いてあまり意味がないので、定常な分布関数はエネルギーと角運動量だけで書けると思っていい。

いちおう、ここで、意味がないというのはどういうことかということを説明しておく。そのためには、意味がある特別な場合というのを考えるのがよい。これは、ケプラー軌道のような、軌道が閉じる場合である。この時には、エネルギーと角運動量の他に、軌道全体の向きを表す量（近点経度）が保存する。これはちゃんと保存量になっている。

しかし、一般には軌道が閉じない。このときでも、近点経度に対応するような保存量が実は存在しているが、それにも関わらず、ある軌道がエネルギーと角運動量で決まる部分空間を覆ってしまう（数学的には、もちろんすべての点を覆えるのではなく、任意の点について、いくらでも近くにいけるというだけだが）。こうなっていると、その積分に分布関数が依存すると、連続性とか微分可能性とかに困難を生じることになる。

さて、

は

と ${\bf J}$ によるということにしたわけだが、いま球対称な場合ということなので ${\bf J}$ の方向にではなく、絶対値だけに依存するのでないといけない。したがって、実は球対称の分布関数は一般に

と書けるということになる。

我々が扱いたいのは自己重力系なので、実際にこれを球対称の場合に書き下してみると

$\begin{displaymath} {1 \over r^2} {d \over dr} \left(r^2{d \Phi \over dr}\right)... ...2}v^2 + \Phi, \vert{\bf r}\times {\bf v}\vert\right) d{\bf v}, \end{displaymath}$

(10)

の場合

上の場合でもまだちょっと大変なので、さらに単純化してとりあえず

にもよらない場合というのを考えてみる。これには、なかなか特別な、空間上の各点で速度分散が等方的であるという性質がある。これはどういうことかというと、一般にある方向の速度分散というのは

$\begin{displaymath} <v_e^2> = {1 \over \rho} \int v_e^2 f(v^2/2 + \Phi)d{\bf v} \end{displaymath}$

(11)

$\begin{displaymath} \Psi = -\Phi + \Phi_0, \quad \quad \quad {\cal E} = -E + \Phi_0 = \Psi - v^2/2 \end{displaymath}$

(12)

$\displaystyle {1 \over r^2} {d \over dr} \left(r^2{d \Psi \over dr}\right)$	$\textstyle =$	$\displaystyle -16\pi^2 G \int_0^{\sqrt{2\Psi}} f(\Psi - {1 \over 2}v^2)v^2 d{\bf v}$
	$\textstyle =$	$\displaystyle -16\pi^2 G \int_0^\Psi f({\cal E})\sqrt{2(\Psi - {\cal E})} d\cal E.$	(13)

これで、一般に

を与えて $\Psi$ を求めるとか、あるいはその逆とかが出来る。

ただし、 $\Psi$ を与えて

を求めようってときには、求まった

が $f\ge 0$ の条件を満たすという保証はないので、そういうのは物理的には意味がない解ということになる。

8 球対称な分布関数の例

ここであげるのはあくまでも例であるが、さまざまな理由からその性質がよく調べられているものである。

8.1 ポリトロープとプラマーモデル

ある意味でもっとも簡単な分布関数の例は、 ${\cal E}$ の冪乗（パワー）で書けるものである。例えば

$\begin{displaymath} f({\cal E}) = \cases{ F{\cal E}^{n-3/2}& (${\cal E}>0$)\cr 0 & otherwize} \end{displaymath}$

(14)

$\begin{displaymath} \rho = c_n \Psi^n \quad (\Psi > 0) \end{displaymath}$

(15)

$\begin{displaymath} {1 \over r^2} {d \over dr} \left(r^2{d \Psi \over dr}\right) + 4\pi G c_n \Psi^n = 0 \end{displaymath}$

(16)

$\begin{displaymath} {1 \over s^2} {d \over ds} \left(s^2{d \psi \over ds}\right) + \psi^n = 0 \end{displaymath}$

(17)

実際には、上の Lane-Emden 方程式を解かないとポテンシャルや密度がどうなっているかはよくわからない。で、一般の

ではこの方程式には初等的な解はないが、

の場合には解があることが古くから知られている。これは

$\begin{displaymath} \phi = {1 \over \sqrt{1 + {1 \over 3} s^2}} \end{displaymath}$

(18)

密度が

で有限で、 $r\rightarrow \infty$ で

より速く落ちるので、質量は有限である。

これは、天文学的になにか素晴らしいものであるというわけではないが、球状星団のうち中心密度が低いものにはまあまあ似ていなくもない。とりあえず、これの意味は、解析関数で簡単に書ける自己重力系の self-consistent なモデルであるということである。

プラマーモデルは、いろんなシミュレーションの初期条件として使われることが多い。

8.2 Hernquist Model

プラマーモデルはその存在が100年くらい前から知られているが、こちらは論文が発表されたのが 1990 年（というわけで、 Binney & Tremaine のときにはまだ知られていなかった）という、非常に新しいモデルである(Hernquist, L., 1990, ApJ 356, 359)。これは、ポテンシャルを

$\begin{displaymath} \Phi = -{1 \over r + a} \end{displaymath}$

(19)

$\begin{displaymath} \rho = C{a^4 \over r (r+a)^3} \end{displaymath}$

(20)

$\begin{displaymath} {d \Phi \over dr} = GM_r/r^2 \end{displaymath}$

(21)

Hernquist Model には、「 $r^{1/4}$ 則をかなり良く再現する」という著しい特徴がある。

$r^{1/4}$ 則については後でその物理的解釈も含めて詳しく議論することにしたいが、要するに、観測される楕円銀河の表面輝度の対数（要するに等級ですね）が、半径の

乗に対して直線にのって見えるというものである。この性質と、一応解析関数で分布関数が書けるということのために、楕円銀河やダークハローのモデルとして広く使われるようになってきている。

ただし、このモデルにはいくつか妙な性質もあり、それについてもまた後で触れることになるはずである。

8.3 等温モデル

前回、無衝突ボルツマン方程式の定常解は熱平衡とは限らないし、そもそもエントロピーが発生しないのだから系が熱平衡に向かって進化するとも限らないという話をした。が、後で出てくるようないくつかの理由から、熱平衡状態について良く理解しておくことは結構大事である。熱平衡状態では（古典統計なので）分布関数はマックスウェル-ボルツマン分布、すなわち

$\begin{displaymath} f({\cal E}) = {\rho_1 \over (2\pi \sigma^2)^{3/2}} e^{{\cal ... ...\sigma^2)^{3/2}} \exp\left({\Psi - v^2/2\over \sigma^2}\right) \end{displaymath}$

(22)

$\begin{displaymath} {2 \over \sqrt{\pi}}\int_0^{\infty} e^{-x^2} dx = 1 \end{displaymath}$

(23)

$\begin{displaymath} \rho = \rho_1 e^{\Psi/\sigma^2} \end{displaymath}$

(24)

$\begin{displaymath} {1 \over r^2} {d \over dr} \left(r^2 {d \Psi \over dr}\right) = -4\pi G \rho \end{displaymath}$

(25)

$\begin{displaymath} {1 \over r^2} {d \over dr} \left(r^2 {d \log \rho \over dr}\right) = -{ 4\pi G \sigma^2}\rho \end{displaymath}$

(26)

$\begin{displaymath} \rho = {\sigma^2 \over 2\pi Gr^2} \end{displaymath}$

(27)

特に、渦巻銀河については、「回転速度が中心からの距離に（あまり）依存しない」（いわゆる flat rotation curve）という性質が知られていて、これを説明するためには上のような $\rho \sim 1/r^2$ のダークハローが必要であるということになっている。

特別ではない解は、中心密度を有限にして中心から外側に向かって解いていけばいい。この時でも、 $r\rightarrow \infty$ の極限では singular isothermal に近付く。

8.3.1 流体との関係

等温モデルは、エントロピー極大であり、分布関数がボルツマン分布になっているという特別な性質がある。このため、等温ガス球と実は同じ構造をとる。以下、ガス球について方程式を導いておく。静水圧平衡の式は

$\begin{displaymath} {dP \over dr} = -\rho {GM_r \over r^2} \end{displaymath}$

(28)

$\begin{displaymath} P = {k_B T \over m} \rho \end{displaymath}$

(29)

$\begin{displaymath} C {d \over dr} \left(r^2 {d \log \rho \over dr}\right) = -{ 4\pi G }\rho r^2 \end{displaymath}$

(30)

なお、ポリトロープでも、ポリトロピックな状態方程式を持つガス球の密度分布と stellar system のそれとは一致する。が、等温モデルの場合とは実は本質的な違いがある。等温モデルの場合は、分布関数そのものが一致する（ボルツマン分布であり、局所的にも大局的にもエントロピー最大）が、一般のポリトロープではそんなことはない（そもそもガス球ではジーンズの定理が成り立たないし、局所的にはボルツマン分布であるから）。

9 ここまでのまとめ

10 球対称ではない場合

10.1 軸対称の場合

球対称ではない場合の例として、まず軸対称の場合を考える。

この場合、運動の積分として確実に存在するのはエネルギーと、対称軸回りの角運動量だけである。但し、では他に運動の積分はないのかというと、話は既に簡単ではなくなる。これは以下のように理解できる。

軸対称ポテンシャルなので、座標を $(R, z, \phi$ ) の円筒座標にしてみる。

を決めると角運動量保存から $V_{\phi}$ が決まるので、こちらは解く必要がなく、

平面の中の運動を考えればよい。そうすると自由度が 2なので、運動の積分がエネルギーの他にもうひとつあると可積分になって、「解析的に解ける」。

$\begin{displaymath} V = aR^2 + bz^2 \end{displaymath}$

(31)

ところが、調和振動子に少し余計な項を加えると、もうよくわからないことが起こる。この有名な例が Hénon-Heiles （エノン-ハイレス）系である。 Hénon も Heiles も天文学者であり、この系は元々円盤銀河のような軸対称ポテンシャルの中での恒星の運動のモデル方程式として導かれたものである。

Hénon と Heiles (1964 AJ 69, 73) は、

として以下のようなものを考えてみた。

$\begin{displaymath} H = \frac{1}{2}(x^2 + y^2) +x^2y - \frac{1}{3}y^3 \end{displaymath}$

(32)

何故この形を考えたかというのは、論文には以下のように書いてある

つまり、物理的あるいは天文学的になにかの意味があるというよりは、このポテンシャルの中での軌道の振る舞いが妙であるということで選ばれたかなり人工的な例ではある。

もちろん、最初の項は調和振動であり明確な意味がある。非線型の第2項の意味はポテンシャルの等高線を書いて見ると分かる。

このように、非線型項によりポテンシャルが三角形になっている。特に

の等高線に注意すると、これは直線になっている。このことから、この三角形の頂点はポテンシャルの鞍点であり、力が 0 になることが分かる。

鞍点であるので、ポテンシャルが下がるのは三角形の内側を向いた方向とその反対側の方向だけであり、それに直交する方向ではポテンシャルは増えている。また、頂点以外では力は 0 でない。

軌道は例えばこんなふうである。これはエネルギー(ハミルトニアン)

$\begin{displaymath} H=\frac{v^2}{2} + U(x,y) \end{displaymath}$

(33)

これだけを見ていると別にどうということはないが、ここでいわゆるポアンカレ断面というものを書いて見ると非常に面白いことになっていると分かる。ポアンカレ断面とは、この場合のような自由度が2のハミルトン力学系で、一方の変数、例えば

が 0 になる時に

の値を書いたものである。

ここで、つながって線になっているのは基本的には一つの軌道であり、分かれた曲線は別の初期条件からの軌道に対応している。エネルギー保存を考えると、

で $v_x^2 \ge 0$ なので、

) はある閉領域の中にくる。保存量がエネルギーだけであれば、任意の初期条件から出発した軌道はこの領域の全ての点をいつかは通ると考えられるが、このようにそうならないということは、なんだかわからない保存量のようなものがあるということを示している。元々の軸対称ポテンシャルに戻ると軸回り角運動量は保存しているが、これは有効ポテンシャルにして回転方向の速度を消去するのに使っているので、2次元問題になった後では無関係である。

軸対称ポテンシャル内での運動にエネルギーと角運動量以外の保存量があるかどうかという問題は、「第3積分問題」と呼ばれる恒星系力学上の現在でも完全には解決されているとはいいがたい問題の一つである。とりあえずこの場合には、第3積分が存在しているように見える。

但し、妙なのは、ポアンカレ断面上の軌跡が交点をもつように見えることである。実際には、この「交点」のところを精密に調べると、奇妙なことがわかる。交点近くの十分に狭い領域では運動の積分が無くなって、ある領域内を軌道が埋めているのである。

この領域はエネルギーを大きくするとどんどん大きくなる。次は

の場合である。

保存量を持つ、規則的な運動をする領域よりも、規則的ではない運動をする領域のほうが広くなってしまっていることが分かる。少なくとも恒星系力学では、このような場合に規則的な軌道を regular orbit, そうでないものを irregular orbit または chaotic orbit という。

キャプションには

よりも大きいようなことが書いてあるが、気の迷いであろう。それはともかく、この場合にはほとんどの領域が規則的でない、つまり、第3積分をもたないような軌道になっていることが分かる。

ここではカオスかどうかを、エネルギー以外の積分がある(ように見える)かどうかで判断したが、実質的に同じことになるもうひとつの定義がある。それは、「十分に近い初期条件から出発した2つの軌道がどのように離れていくか」である。より厳密にはこれはリヤプノフ指数が正かどうかということになる。定義の詳細はここでは省くが、大雑把にはリヤプノフ指数が正なら十分近い2つの軌道が指数関数的に離れていく。このことは上の Hénon-Heiles 系でも観察され、周期的に見える領域では軌道は時間の1次で離れていくが、カオス的な領域では指数関数的に離れる。

というわけで、軸対称でも話は球対称の場合よりはるかに難しくなる。実際に理論的に軸対称な系、例えば楕円銀河とか回転がある球状星団とかを扱う時には第三積分はないものとしてジーンズの定理を使って分布関数を構成することはできるが、実際には多くの場合に第三積分的なものはあるのでそれ以外の解がないというわけではない。

なお、2次元円盤の場合は話が簡単で、全ての軌道は可積分である。なので、薄い円盤の場合にはほぼ可積分として扱うことができる。

10.2 軸対称でもない場合

多くの楕円銀河は実は軸対称でもなく、楕円の3軸の長さが全て違う。そうすると、角運動量の3成分が全て保存しなくなり、自明な保存量はエネルギーだけになる。

しかし実際にはその中で regular orbit と呼ばれる、エネルギー以外の保存量を持つ軌道があり、それによって3軸不等な形状を維持していると考えられている。

というわけで、今回の講義では球対称な系を中心に扱う。後半では円盤の話もする。

この文書について...

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

翻訳は Jun Makino によって平成25年7月17日に実行されました。