pt 称

重力多体系の進化

牧野淳一郎

1 ２体緩和（続き）

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

というわけで、いよいよバックグラウンドが動いている場合ということになる。 この時でも、相対速度の変化自体は前節に述べたもので正しいが、相対速度に フィールド粒子の速度成分が入ってくることになる。

以下、二種類の単位ベクトル系をとって、その上で考える。一つは $(\mbox{\boldmath$e$}_1, \mbox{\boldmath$e$}_2, \mbox{\boldmath$e$}_3)$であり、元の空間に固定されている。もう一つは $(\mbox{\boldmath$e$}_1', \mbox{\boldmath$e$}_2', \mbox{\boldmath$e$}_3')$であり、最初の成分を相対速度 $\mbox{\boldmath$V$}$に平行にとる。従って、 後者は相手の粒子によって違うわけである。この２つを考えることで、相対速 度の変化をもとの静止系でのテスト粒子の変化に焼き直す。

まず、1次の項は相対速度に平行な成分だけであった。このことから、ある方 向の速度変化は

$$<\Delta v_i> = <\mbox{\boldmath$e$}_i \cdot \mbox{\boldmath$e$}_1' \Delta v_{平行}>'$$ (1)

ということになる。もうすこし精密に書くと、右辺はインパクトパラメータと 相手の速度の積分なので、以下のように書けることになる。

$$<\Delta v_i> = \int dp2\pi pV\int d\mbox{\boldmath$v$}_f f(\mbox{\boldmath$v$}_f) \Delta v_i$$ (2)

$$= -\Gamma(1+m/m_f)\int {f(\mbox{\boldmath$v$}_f) \over V^2}(\mbox{\boldmath$e$}_i\cdot \mbox{\boldmath$e$}_1') d\mbox{\boldmath$v$}_f$$ (3)

ここで、 $\mbox{\boldmath$v$}_f = \mbox{\boldmath$v$}_t - \mbox{\boldmath$V$}$ はフィールド粒子の速度、 $f$は速度分布 関数である。 $p$についての積分を先にやったことに注意して欲しい。この積 分は $v_f = 0$ であった時の結果をそのまま使っている。

さて、2次の項についてであるが、前節で見たように $\mbox{\boldmath$v$}$ に垂直な成分、す なわち $\mbox{\boldmath$e$}_2'$と $\mbox{\boldmath$e$}_3'$の成分を考えればいい。従って、一つの方向から くるフィールド粒子との散乱を考えた時には

$$<\Delta v_i \Delta v_j> = [(\mbox{\boldmath$e$}_i\cdot\mbox{\boldmath$e$}_2')\mbox{\boldmat... ...h$e$}_j\cdot\mbox{\boldmath$e$}_3')\mbox{\boldmath$e$}_3']<\Delta v_{垂直}^2>/2$$

$$= [(\mbox{\boldmath$e$}_i\cdot\mbox{\boldmath$e$}_2')(\mbox{\boldma... ...e$}_3')(\mbox{\boldmath$e$}_j\cdot\mbox{\boldmath$e$}_3')]<\Delta v_{垂直}^2>/2$$

$$= Q_{ij}<\Delta v_{垂直}^2>/2$$ (4)

これもまた分布関数を掛けて積分すると、結局

$$<\Delta v_i \Delta v_j> = \Gamma \int {f(\mbox{\boldmath$v$}_f) \over V}Q_{ij} d\mbox{\boldmath$v$}_f$$ (5)

ということになる。これで一応必要な2次までの係数はすべて書けたわけだが、 あまり計算するのに使い易い形ではない。というのは、 $\mbox{\boldmath$e$}_i'$ とか $V$ と かいったものがまだややこしい形ではいったままであるからである。しかし、 もうちょっと簡単な形に書き直せることが知られている。まず、1次の項だが、

$${ \partial \over \partial v_i} \left({1 \over V}\right) = -{... ...-{\mbox{\boldmath$e$}_i\cdot \mbox{\boldmath$e$}_1' \over V^2}$$ (6)

という都合のよい関係がある。変形はとくにややこしいところはないと思う。 $\mbox{\boldmath$e$}_1'$をそもそも $\mbox{\boldmath$V$}$に平行にとったから上のように出来るわけである。 このため、

$$h(\mbox{\boldmath$v$}) = \int {f(\mbox{\boldmath$v$}_f) \ove... ...ldmath$v$}- \mbox{\boldmath$v$}_f\vert} d\mbox{\boldmath$v$}_f$$ (7)

なる関数 $h(\mbox{\boldmath$v$})$を導入して、

$$<\Delta v_i> = -\Gamma(1+m/m_f){\partial h \over \partial v_i}$$ (8)

ということになる。

2次の項についても同様な整理が可能である。$Q_{ij}$は $\mbox{\boldmath$e$}_i$, $\mbox{\boldmath$e$}_j$ の $\mbox{\boldmath$e$}_2'$と $\mbox{\boldmath$e$}_3'$によって張られる平面への写像の内積なので、 $\mbox{\boldmath$e$}_1'$と の内積の分をつけてやれば元の単位ベクトル同士の内積になる。つまり、

$$Q_{ij} = \delta_{ij} - (\mbox{\boldmath$e$}_i\cdot\mbox{\bol... ...dot\mbox{\boldmath$e$}_1') = \delta_{ij} - {V_i V_j \over V^2}$$ (9)

したがって、

$${\partial^2 V \over \partial v_i \partial v_j } = {1 \over V}\left(\delta_{ij} - {V_i V_j \over V^2}\right) = Q_{ij} / V$$ (10)

というわけで、

$$g(v) = \int f(\mbox{\boldmath$v$}_f)\vert\mbox{\boldmath$v$}- \mbox{\boldmath$v$}_f\vert d\mbox{\boldmath$v$}_f$$ (11)

とおけば、

$$<\Delta v_i \Delta v_j> = \Gamma {\partial^2 g \over \partial v_i\partial v_j }$$ (12)

1.2 バックグラウンド速度分布が熱平衡の場合

前節では、バックグラウンドの速度分布が任意のものについて、実際に計算可 能な式を導いた。ここでは、速度分布が等方的な場合につ いて式を単純化してみる。

速度分布が等方的な場合、$h$や$g$の積分を、 $\mbox{\boldmath$v$}_f$の絶対値方向と角度方 向に分けることができる。角度方向の積分については、 $\mbox{\boldmath$v$}$ と $\mbox{\boldmath$v$}_f$の なす角度を $\theta$とし、 $\mu = \cos \theta$ とすれば、球面上での積分 が、まず$h$については

$$\int {d\mbox{\boldmath$v$}_f \over \vert\mbox{\boldmath$v$}-... ...{1/2}} = 4\pi \cases{1/v, &($v>v_f$)\cr 1/v_f, &($v<v_f$)\cr}$$ (13)

となる。これは、球面上に分布する電荷の作るポテンシャルと同じ式である。 $g$についても同様に計算できて

$$\int { \vert\mbox{\boldmath$v$}- \mbox{\boldmath$v$}_f\vert}... ... \cases{3v+v_f^2/v, &($v>v_f$)\cr 3v_f+v^2/vf, &($v<v_f$)\cr}$$ (14)

となる。これから、

$$F_n(v) = \int_0^v\left({v_f \over v}\right)^n f(v_f)dv_f$$

$$E_n(v) = \int_v^{\infty}\left({v_f \over v}\right)^n f(v_f)dv_f$$ (15)

というものを考えると、

$$h(v) = 4\pi v[F_2(v) + E_1(v)]$$

$$g(v) = {4\pi v^3 \over 3} [3F_2(v) + F_4(v) + 3E_3(v) +E_1(v)]$$ (16)

これらから、最終的な結果、すなわち、バックグラウンドが動いている時の、 速度に平行な速度変化と垂直なそれを書き下せることになる。それらは、結局、

$$<\Delta v_{平行}> = -4 \pi \Gamma\left( 1 + {m \over m_f}\right)F_2(v)$$ (17)

$$<\Delta v_{平行}^2> = {8 \pi \Gamma v \over 3} [F_4(v) + E_1(v)]$$ (18)

$$<\Delta v_{垂直}^2> = {8 \pi \Gamma v \over 3} [3F_2(v)-F_4(v) + 2E_1(v)]$$ (19)

これらから、粒子のエネルギーの変化 $\Delta E$ を出すことができる。

$$\Delta E = v\Delta v_{平行} + <\Delta v_{平行}^2>/2+ <\Delta v_{垂直} ^2> /2$$ (20)

と書けるので、1次の項は

$$<\Delta E> = 4 \pi \Gamma v\left[E_1(v) -{m \over m_f}F_2(v)\right]$$ (21)

となる。2次の項については、 $(v\Delta v_{平行})^2$以外の項は小さいので 無視すると

$$<\Delta E^2> = {8 \pi \Gamma v^3 \over 3} \left[F_4(v) +E_1(v)\right]$$ (22)

となる。

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

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

とすると、上の係数等を具体的に計算できることになって、その形は

$$<\Delta v_{平行}> = -4 { n_f\Gamma\over \sigma^2}\left( 1 + {m \over m_f}\right)G(x)$$ (24)

$$<\Delta v_{平行}^2> = 2\sqrt{2}{n_f\Gamma \over \sigma } G(x)/x$$ (25)

$$<\Delta v_{垂直}^2> = 2\sqrt{2}{n_f \Gamma\over \sigma }{{\rm erf}(x) - G(x)\over x}$$ (26)

$$<\Delta E> = \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]$$ (27)

ここで ${\rm erf}$は誤差関数であり、

$$G(x) = {{\rm erf}(x) - x{\rm erf}'(x) \over 2x^2}$$ (28)

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

山ほど式はでてきたものの、全然なんだかわからないという気分になった人も まあいるのではないかと思うので、以下、上の式の意味についてちょっと考え てみよう。

まず、速度の1次の項を見てみる。これは、速度分布には$F_2$だけ を通して依存しているということに注目して欲しい。例えば、マックスウェル 分布のようなものを考えた時、 $v$が大きい極限では $F_2\sim 1/(2\pi v^2)$となるので、回りが止まっているときと同じく速度変化は速度の2乗に 反比例する。これに対して、$v$が小さい極限では、$f$を一定と見なすことが 出来るので $F_2 \propto v$ となる。

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

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

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

なお、、自己重力系ではこのエネルギー交換の結果熱平衡に向かうとは限らな いということに注意する必要がある。つまり、重いものがエネルギーを失い、 軽いものがエネルギーを得るということは、それぞれの分布関数が変わり、空 間分布も変わるということである。具体的には、重いものは中心に落ちるし、 軽いものは外側に押し出される。

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

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

「タイムスケール」というのは、普通ある量 $x$ の時間変化が

$${dx \over dt} = (-) {x \over T}$$ (29)

なる形で書ける時の $T$ のことである。もちろん、一般には$T$は$x$を含む いろんなものの関数である。

そういった意味で考えやすいのは、温度（平均運動エネルギ）が違う２つの空間一様な 分布が重なりあっている時に、どのようにして２つが近付いていくかというの ものである。これは、もちろん時間が立てば熱平衡に近付くわけである。以下、 実際に計算してみる。

2.1 等分配のタイムスケール：理論

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

$${d<E_t> \over dt} = 4\pi \int v_t^2f(v_t)<\Delta E_t>dv_t$$ (30)

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

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

となる。

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

$${d\log <E_t> \over dt} = -\sqrt{3/\pi} {m_t n_f \Gamma \over m_f v^3}$$ (32)

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

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

$${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}$$ (33)

となる。 ここで

$$\Gamma = 4\pi G^2 m_f^2 \ln \Lambda$$ (34)

を使った。 さらに、 $E_t = m_tv_t^2/2$ などを使って書き直せば

$${d(v_t^2) \over dt} =4 \sqrt{6\pi}G^2\ln \Lambda n_f m_f^2 v_f^{-1}$$ (35)

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

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

$$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}$$ (36)

とするものである。 ここで $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$といった関係 を使えばでてくる。これは

$$t_{rh} = 0.138 {Nr_h^{3/2} \over M^{1/2} G^{1/2}\log \Lambda}$$ (37)

となる。

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

これがもっともクリティカルに効いて来るのは、実際の天体においてというよ りはむしろシミュレーションにおいてである。これについては後で具体例で議 論しよう。

2.2 等分配のタイムスケール：実験

さて、実際に数値実験でエネルギー等分配に近付く過程を見て、これまで理論 的に考えてきたもの、特に積分の上限がどうなっているかみようというわけだ が、これはそれほど簡単ではない。というのは、「空間分布が無限一様」とか 「等温」とかいう初期条件が設定できないためである。

一例として、 Farouki and Salpeter (1982 APJ 253, 512, 1994 APJ 427, 676)の実験を取り上げてみる。彼らは、実際に２種類の違う質量の粒子からなる、 初期に力学平衡にある系を考えた。初期には単位質量あたりの運動エネルギー が同じであるようにした（つまり、空間分布、速度分布ともに同じということ）。 したがって、重い方が熱平衡に比べて余計に運動エネルギーを持っていること になる。

正確にいうと、実際に力学平衡にある系を作ったわけではなく、適当に球とか 立方体のなかに一様に粒子をばらまいて、落ちつくまで待ってから使うという方法を とっている。速度は、全体としてビリアル平衡になるようにマックスウェル分 布を与えた。

さて、この系を進化させると何が起こるかをまず考えてみる。重い方の温度が 高いので、平均的には重い方から軽い方にエネルギーが流れるであろう。しか し、その結果温度が近付くであろうか？これはとりあえず何ともいえない。と いうのは、ポテンシャルが一様ではないのでエネルギーの変化は粒子の速度分 布だけでなく密度分布、すなわち粒子の位置と重力エネルギにも影響するから である。エネルギーの変化が温度にどう影響するかは、もともとの粒子の分布 にも依存するわけである。

さらにややこしいことに、 $\log \Lambda$ がなんであるかという問題もある。 通常、$N$体計算においては、2粒子間の重力ポテンシャル $\phi_{ij}= -Gm_i m_j /r_{ij}$をそのまま使わないで、典型的には

$$\phi_{ij} = -G{ m_i m_j \over {r_{ij}^2 + \epsilon^2}}$$ (38)

という形に「ソフトニング」したものを使う。これは、数値計算上の困難を避 けるのが第一義的な理由である。つまり、純粋な $1/r$ポテンシャルを使うと、 2つの粒子が無限に近付き得る。無限に近付くこと自体は問題ではないが、こ の時速度と重力エネルギーも無限大になる。従って、非常に高い精度で計算し ても、エネルギーの誤差が非常に大きくなってしまう。

もちろん、近付いたら解析解に切替える、あるいは摂動法にする、または座標 変換をして特異性を消すといった方法があるが、どれもかなり面倒である。そ こで、多くの数値実験で、ポテンシャルをすこし変えて原点で発散しないよう にする。この場合、インパクトパラメータが $\epsilon$ よりも小さい時には ほとんど曲がらなくなるので、速度変化の積分の下のほうを修正する必要が起 きる。大雑把にいって $p<\epsilon$ の寄与を無視すればよいことになる。

図1に典型的な結果を示す。Farouki and Salpeter では、温度と かエネルギー変化を直接測定するのは断念し、その代わりに軽い粒子の half mass radius をとった。彼らは緩和時間として

$$T_R = {R_0 \over dR/dt}$$ (39)

つまり、半径の変化のタイムスケールをもって緩和時間であるとした。これは もちろん定数を別にすれば上で考えた定義と一致するものになっている。

図: 緩和過程

この結果から、彼らは $\log \Lambda$ と粒子数 $N$ について図2のような 関係を得た

図: $\log \Lambda$ と粒子数 $N$

上と下は $\epsilon$ の大きさが違い、上は$R/N$の程度、下は$R/N^{1/3}$の 程度である。ここで$F$は緩和時間を

$$t_{rh} = F {N \over 26\log \Lambda} t_d$$ (40)

と書くことにしたために出てくる定数である。26という定数は、前に出てきた $t_{rh}$ を $t_d$ を使って書き直したので、定数の $0.138$が違う値になっ ているだけのものである。

彼らの主な主張は、「この振舞いは上のカットオフがシステムサイズであると いう理論と一致している」というものである。それはまあかなり確かといって もよいと思われる。が、 $F$の値を推定すると 10 程度になり、かなり大き過 ぎるものになっている。これは、定義の違いもあるので評価はちょっと難しい。

もうちょっと良くあっている例を紹介しておこう。

図: 2体緩和による系の構造進化

図 3は Spurzem and Takahashi (1995, MNRAS 272, 772) によるものである。 ここでは、緩和時間といった定義のはっきりしないものを比べているのではな く、系の構造そのものの変化を、$N$体計算の結果とモーメント方程式から導 いた（球対称、等方の）フォッカー・プランク方程式の数値解とで比べている。 結果の一致は素晴らしいものである。

3 フォッカープランク近似

さて、数値計算の精度がどうという話を別にすれば、2体緩和を考える理由はそ れにより系がどう進化するかを理解するということである。そのためには、粒 子の速度変化のモーメントの式から、分布関数の変化についての方程式を導く ことが有用であろう。

粒子の物理量の変化の1次と2次のモーメントから（高次の項の寄与を無視して） 分布関数に対する移流拡散方程式を導くことができる。この操作を通常フォッ カープランク近似といい、でてきた方程式をフォッ カープランク方程式という。以下、その導出を行なってみる。 これは、ボルツマン方程式の衝突項を求めることに対応する。

今、 $P(\mbox{\boldmath$v$}, \Delta \mbox{\boldmath$v$})d\Delta v$ が、速度 $\mbox{\boldmath$v$}$の粒子が、時間$\Delta t$の間に速度変化 $\Delta \mbox{\boldmath$v$}$ （の近傍の $d\Delta \mbox{\boldmath$v$}$の範囲）を受ける確率であるとする。 すると、 $\Delta t$ たった後の分布関数は以下のように書ける

$$f(\mbox{\boldmath$v$}, t + \Delta t) = \int f(\mbox{\boldmat... ...th$v$}, \Delta \mbox{\boldmath$v$})d\Delta \mbox{\boldmath$v$}$$ (41)

これは $\mbox{\boldmath$v$}$ のところにやってくる可能性があるものすべてについての和をとっ ただけである。ここで、両辺を展開することを考える。左辺は

$$f + {\partial f \over \partial t} \Delta t$$ (42)

である。右辺は

$$\int\left[ fP - \sum_{i=1}^3 {\partial (fP)\over \partial v_... ... \Delta v_i\Delta v_j + ...\right]d\Delta \mbox{\boldmath$v$}$$ (43)

ここで、

$$<\Delta v_i> \Delta t = \int P\Delta v_i d\Delta \mbox{\boldmath$v$}$$ (44)

$$<\Delta v_i\Delta v_j> \Delta t = \int P\Delta v_i \Delta v_j d\Delta \mbox{\boldmath$v$}$$ (45)

とおけることを使えば、微分と積分の順序を入れ換えて

$${\partial f \over \partial t} = - \sum_{i=1}^3 {\partial (f ... ...^2 (f<\Delta v_i \Delta v_j>) \over \partial v_i\partial v_j}$$ (46)

これで、理屈の上では分布関数の変化が計算できるということになる。もちろ ん、実際にこれを解いて自己重力系の進化を調べるのは、必ずしも容易ではな い。その理由は、分布関数が6次元位相空間上で定義されること、それがジー ンズの定理を満たすように進化しなければいけないことである。具体的には、

• そもそも、球対称とかの仮定なしでは運動の積分だけの関数であるよう な分布関数が求まらない

• 実際に運動の積分の関数として分布関数が書ける場合でも、拡散係数を 求めるのはそれほど容易ではない。

• 実際に拡散係数が求まっても、自己重力系なので分布関数が変化すると 同時にポテンシャルも変化することになり、それによって粒子の軌道が断熱的 に変化することになる。そのへんを考慮してつじつまがあうような計算プログ ラムを作るのは容易ではない。

というような困難があり、従来は、分布関数をエネルギーだけの関数と近似す る計算しか行なわれていなかった。1994年頃に、Takahashi が初めて $f(E,J)$の場合に信頼できる結果を得ることに成功した。しかし、これも、拡 散係数を求める時に、フィールドの分布は等方的であるという近似を行なって いる。

なお、モンテカルロ法など、もうちょっといい加減な方法ではある程度のこ とは出来ている。この辺についてはまた後でもう少し詳しく触れる。

4 熱力学的進化

ここからは熱力学的な進化ということで、いくつかの理想化された 場合を考えておく。なお、以下はまた流体近似で話をする。自己重力質点系の 熱力学的な進化をそのまま考えるのはほとんど不可能だからである。原理的に は、今やったような Fokker-Planck 方程式の平衡解をつくって、その安定性 を調べればできるわけだが、これまでそのような研究はあまり行なわれていな い。また、熱平衡状態とその安定性に関するかぎり、ガスか質点系かという違 いには意味がない。まあ、そういうわけで、ガスで考えるということにもそれ なりの意味はある。

もちろん、平衡状態からのずれの時間発展はガスか質点系かで大きく違うわけ だが、おもな違いは熱伝導の係数の密度、温度への依存の違いとして理解でき る。

5 等温状態の安定性

まず、もっとも単純な場合ということで、等温の平衡状態を考える。話を簡単 にするために、球対称の 断熱壁の中のガスということにする。平衡状態は、静水圧平衡の式

$${dp\over dM } = -{M \over 4 \pi r^4},$$ (47)

$${dr\over dM} = {1 \over 4\pi r^2\rho},$$ (48)

を考えればいい。例によって球対称で、$M(r)$は半径$r$の中の質量、$p$と $\rho$は圧力と密度である。面倒なので、単位系として $G=M=R=1$となるよう にとる。$G$は重力定数、$M$は壁（断熱壁）のなかの全質量、$R$は断熱壁の 半径である。

温度は、状態方程式

$$p = \rho T$$ (49)

で決まる。単位質量当たりのエントロピーは

$$s = \ln(T^{3/2}\rho^{-1})$$ (50)

であり、境界条件は

$$r = 0\quad {\rm for}\quad M=0,$$

$$r = 1\quad {\rm for}\quad M=1.$$ (51)

である。この解自体は、すでに何度か扱ったように数値的に求めることが出来 る。

以下、熱力学的安定性について議論するわけだが、これにはいろんな流儀があっ て、なんだかよく理屈がわからないものもある。等温の自己重力ガスの安定性 について初めて議論したのは Antonov (1961) であり、もうちょっと詳しい議 論が Lynden-Bell & Wood (1968) によってなされた。しかし、これらはいず れも真の意味での安定性解析、すなわち、平衡解に対する摂動を考え、それが 成長するかどうかをしらべたというものではない。

そのような意味での安定性解析を初めて適用したのは、 Hachisu & Sugimoto (1978) である。彼らの方法は、大雑把にいうと以下のようなものである。

• 平衡解の回りの線形化した摂動に対する方程式を作る

• 適当な制約条件のもとで、全エントロピーの2次の変分$\Delta S^2$を最大にするよ うな摂動を求める。これは、具体的な手続きとしては固有値問題を解くことに 帰着される。

• 仮に求まった$\Delta S^2$が正であれば、これは、平衡状態がエントロ ピー極大の状態ではない、したがって、不安定点であるということを意味する。

等温状態なのでエントロピーは通常ならば極大値である。これは、任意のエン トロピーの再分配に対して、 $\Delta S=0, \Delta S^2<0$となっているとい うことを意味する。

これは、熱をちょっとどこかからとって別のところに与えると、それによる温 度変化を考えなければ（一次の変分）エントロピーは変わらない。また、温度 変化を考えると（二次の変分）、熱をもらった方は温度が上がっているのでも らうエントロピーは少なく、出したほうは逆に温度が下がるので出ていくエン トロピーが多い、従って、系全体としては普通は摂動を与えるとエントロピー が減る、すなわち、平衡状態はエントロピー極大に相当している。

以上から、もし、熱を取り去った時に温度が上がるようなことがあればエント ロピー極大ではないかもしれないということが想像できよう。もちろん、常識 的な熱力学の対象ではそんなことはあり得ないわけだが、自己重力系ではそう ではないというのはすでにビリアル定理のところでやった通りである。

つまり、自己重力系を全体として考えると、熱を奪うと系が小さくなり、単位 質量あたりの運動エネルギー、すなわち温度が大きくなるわけである。

断熱壁で囲んだ系では話はもう少しややこしいが、実際に、十分温度が低い、 重力の影響が大きいような系では$\Delta S^2$が正になるということを示した のが Hachisu & Sugimoto である。

この解析は非常に見事なものなので、是非元論文 (PTP 60, 13) を読んで見て 欲しい。まあ、それはそれとして、ここではもう少し違った解析方法をとって みよう。

5.1 等温状態からの時間発展

Hachisu & Sugimoto の方法では、摂動に対して$\Delta S^2$を求め、その符 号から安定か不安定かを決めている。この方法では、もちろん、熱力学的に安 定かどうかをきめることは出来るが、不安定性がどのように発展するかを調べ ることはできない。というのは、そのためには熱伝導の式もカップルさせて線 形応答を求めないといけないのに、そのような解析は行なっていないからであ る。というわけで、しばらく前にそういう解析をやってみた (Makino and Hut 1991、APJ 383, 181) ので、今日はその結果に基づいて話す。

熱伝導の式は

$$K{\partial T\over \partial r} = - {L \over 4\pi r^2}$$ (52)

と書ける。ここで、 $L(r)$は半径 $r$のところでの熱流束であり、 $K$ は熱 伝導の係数である。$K$は温度、密度の関数だが、ここでは等温に近いので密 度だけの関数として

$$K = \rho^{\alpha}$$ (53)

という形を仮定する。放射伝達であれば $\alpha=-1$である。

自己重力質点系の場合は、密度が高いほうが緩和が 速かった。このことを熱伝導係数でむりやりに表現すると、$\alpha=1$となる。

これは以下のように考えたことになっている。

速度分散が同じなら緩和時間 $T$ は単純に密度 $\rho$ に反比例する。自己重力系を考えると、 やはり速度分散が同じなら系の特徴的な大きさ $R$ は質量 $M\sim \rho R^3$ に比例するので、 $\rho \sim R^{-2}$ なる関係がある。

緩和時間を温度勾配と単位面積当りの熱流束の関係に直してみると、温度勾配 は $1/R \sim \rho^{1/2}$ の程度、 単位面積当りの熱流束は $\rho R/T \sim \rho^{3/2}$ の程度であ る。従って $K\sim \rho$ ということになる。

エントロピーについての式は

$${\partial L\over \partial r} = - 4\pi r^2\rho T {\partial s \over \partial t}\big{\vert}_{M}.$$ (54)

で与えられ、境界条件は

$$L=0\quad {\rm for}\quad M=0\quad{\rm and}\quad M=1.$$ (55)

ということになる。

5.2 線形化した方程式

微小な摂動に $\delta$ をつけることにして、線形化した方程式は

$${d\delta\ln p \over dM} = {M\over 4 \pi pr^4} (\delta\ln p + 4\delta \ln r),$$ (56)

$${d\delta\ln r \over dM} = -{1\over 4\pi r^3\rho}(3\delta\ln r + \delta \ln \rho),$$ (57)

$$\delta\ln p = \delta \ln \rho + \delta\ln T,$$ (58)

$$\delta s = {3\over 2}\delta\ln T - \delta \ln \rho,$$ (59)

境界条件は

$$3\delta \ln r + \delta \ln \rho =0\quad {\rm for}\quad M=0$$ (60)

$$\delta\ln r = 0 \quad\quad\quad\quad\quad {\rm for}\quad M=1$$ (61)

ということになる。

熱流束とエントロピーの変化については、もとの式が線形なのでそのまま使え る。つまり、 $L=\delta L$ なので、 ${\partial T\over \partial r} = {\partial \delta T\over \partial r}$ であり、また ${\partial s \over \partial t} = {\partial \delta s \over \partial t}$となる。これらは始 めから一次の微小量しか含んでいない。

この方程式系に対して、

$$\delta \ln p = \delta \ln p_0e^{\lambda t},$$

$$\delta \ln r = \delta \ln r_0e^{\lambda t},$$

$$\delta \ln \rho = \delta \ln \rho_0e^{\lambda t},$$ (62)

$$\delta \ln T = \delta \ln T_0e^{\lambda t},$$

$$\delta L = \delta L_0e^{\lambda t},$$

$$\delta s = \delta s_0e^{\lambda t}.$$

という形をした解を捜すわけである。ここで、添字 $0$がついたものは時間発 展解の空間依存性を表す。線形化したので固有モードが出てきて欲しいなとい う解析をしているわけである。これらを線形化した方程式に入れると、固有値 問題

$${d\delta\ln p_0 \over dM} = {M\over 4 \pi pr^4}(\delta\ln p_0 + 4\delta \ln r_0)\quad ,$$ (63)

$${d\delta\ln r_0 \over dM} = -{1\over 4\pi r^3\rho}(3\delta\ln r_0 + \delta \ln \rho_0)\quad ,$$ (64)

$$KT{d \delta \ln T_0\over d M} = - {\delta L_0 \over (4\pi r^2)^2\rho_0} \quad , {d \delta L_0\over d M} = - \lambda T \delta s_0\quad ,$$ (65)

$$\delta\ln p_0 = \delta \ln \rho_0 + \delta\ln T_0\quad ,$$ (66)

$$\delta s_0 = {3\over 2}\delta\ln T_0 - \delta \ln \rho_0\quad .$$ (67)

が出てくる。 ただし、境界条件は

$$3\delta \ln r_0 + \delta \ln \rho_0 =0\quad {\rm for}\quad M=0$$ (68)

$$\delta\ln r_0 = 0 \quad\quad\quad\quad\quad {\rm for}\quad M=1$$ (69)

$$\delta L=0\quad\quad\quad\quad\quad\quad {\rm for}\quad M=0$$ (70)

$$\delta L=0\quad\quad\quad\quad\quad\quad {\rm for}\quad M=1$$ (71)

ということになる。

さて、これを解かないといけないわけだが、数値解法まで触れている余裕がな いので詳細は省く。我々はいわゆる shooting method を使ったが、本当は緩 和法のほうが安全であったかもしれない。

以下、どういう答が求まり、それはどういうものかということを簡単にまとめ よう。

5.3 安定領域

面倒なので以下 $\alpha=1$の場合だけを考える。以下に示すのは第一固有値 （ここではすべての固有値が負なので、最も0に近いもの）に対応する固有関 数である。

=5cm =5cm =5cm

ここで、 $D$は中心の密度と壁のすぐ内側での密度の比である。$D=1$という のは、温度が無限に高くて重力エネルギーが相対的に小さい極限である。これ に対し、 $D=\infty$ は singular isothermal に対応する。

$D=1.05$ は、要するに重力が無視できる場合である。この時はもちろん応答 はベッセル関数かなにかで書ける。注意して欲しいことは、圧力の変化がない こと、エントロピーと温度がちゃんと比例関係にあることである。これは、重 力が無視できるので普通の振舞いをしているわけである。つまり、密度、温度 の変化が圧力変化がなくなるように働く。これは、静水圧平衡を保つためであ る。

なお、ここでは、中心から熱を奪って外に与えるようなものを考えているが、 その逆も固有関数であることに注意してほしい。これは、線形化した方程式の 解だからである。

さて、少し中心密度を上げると、摂動に対する圧力の応答が変わって、中心で 圧力が上がるようになる。これは、熱を奪われることに応答して縮むと、重力 も強くなるので、つじつまをあわせるにはもうすこし縮んで圧力を上げる必要 が起きるからである。このために、温度の応答は与えたエントロピーからはず れてくる。もっとどんどん温度をさげて、$D$を大きくすると、ついには、熱 を奪ったにもかかわらず、温度が中心でも上昇するようになる。

もちろん、この解は負の固有値に対応するものであり、いぜんとして安定であ る。それは、温度勾配としては依然として中心に向かって下がっていて、ちゃ んとエントロピー変化を打ち消す向きに熱がながれるからである。

5.4 中立安定

さて、もっと $D$を大きくすると、ついには固有値が $0$、すなわち与えられ た摂動が減衰しなくなる。この状況を以下に示す。

=8cm

与えられた摂動が減衰しないということは、温度勾配ができないということで ある。実際、応答は $\delta T$ が定数になっている。

5.5 重力熱力学的安定

さらにもっと温度を下げ、 $D$を大きくすると、ついには固有値が正になる。 以下にいくつかの例を示す

=5cm =5cm =5cm

=5cm =5cm

たくさん例を出したが、どの場合でも中心でエントロピーが減っているのに温 度が大きく上がり、それが外側の温度上昇を追い越している。その結果中心か ら外に向かう熱流ができるのである。

なお、ちょっと注意して欲しいのは、$\alpha$ の値によって応答が大きく違 うことである。$D$が同じ時、$\alpha=1$の方が$\alpha=-1$に比べて中心に集 まったような応答になっている。これは、熱伝導が密度の高いところで速いた めと考えて良い。

なお、以下で関係するのでここで述べておくが、$D$の値と系の全エネルギー の間の関係は単純ではない。安定領域では$D$を大きくするためには系を冷や せばよく、したがってエネルギーと$D$は一対一に対応する。しかし、$D=709$ の中立安定点はエネルギーの極小値になっていて、これよりエネルギーの低い 熱平衡解は存在しない。言い換えれば、$D>709$ の解には、それと同じエネル ギーをもった安定解が常に存在する。

5.6 Hachisu and Sugimoto の解析の意味

HS は、安定性の問題を $\delta s$ に対する変分問題として定式化した。具 体的には、局所的なエントロピーの2次の変分(1次の変分は定義により0なので) を計算し、それを最大化する $\delta s$ を固有値問題を解く形で求めている。 つまり、

$$\int_0^1 \delta s dM = 0,$$ (72)

$$\int_0^1 \delta s^2 dM = 1.$$ (73)

という制約を加えた上で、系全体のエントロピー変化を最大化しているわけで ある。

これに対応する固有値問題は

$$\Lambda \left[{d\over dM}\left({5p^2r^7 \over TM}{d\over dM}\right)-{3\over 8}\right]\delta s + \delta s + {3 \over 8} = 0,$$ (74)

ちょっと書き換えると

$$\xi{d\over dM}\left({5p^2r^7 \over TM}{d\over dM}\right)\delta\ln T = \delta s.$$ (75)

熱伝導からの我々の式を同じ形にすると

$${d\over dM}\left(4\pi r^4\rho K {d\over dM}\right)\delta\ln T = \lambda\delta s.$$ (76)

つまり、 HS の解析は

$$K = {\rho r^3 \over M}.$$ (77)

という形の熱伝導を考えたことに相当する。これはまあ大雑把にいって熱伝導 が密度に依存しないわけで、 星の場合と恒星系の場合の中間的なものになっ ている。

この文書について...

重力多体系の進化

この文書はLaTeX2HTML 翻訳プログラム Version 2002-2-1 (1.71)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds,
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

を日本語化したもの( 2002-2-1 (1.71) JA patch-1.9 版)

Copyright © 1998, 1999, Kenshi Muto, Debian Project.

Copyright © 2000, Jun Nishii, Project Vine.

Copyright © 2001, 2002, Shige TAKENO, Niigata Inst.Tech.

Copyright © 2002, KOBAYASHI R. Taizo, Project Vine.

を用いて生成されました。

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

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