数理解析研究所講究録 1413 巻 2005 年 186-192

186

長距離相互作用粒子系における不安定性、及び自発的構造形成
Instability and Selforganization in
Particle system with Long-range Interaction
東京大学・工学系研究科

鈴木将 (Masaru Suzuki)

論 (Satoshi Yukawa)伊藤伸泰 (Nobuyasu Ito)

湯川

School of Engineering,
The University of Tokyo

我々は、 長距離相互作用粒子系、 特に回転を持つディスク状重力系について、 その安定性及び線形波の分散関係を数値的に求めることを試みた。 数値計算には通常の粒子の代わりに、 平行に並べられた無限に長い棒状粒子を用いた 体計算を行$\mathrm{N}$うことで、 波の伝搬方向を 1 次元的に制約しつつ計算を効率化した。安定性に関する結果として、 このモデル化による特殊性が現われることなく、 速度が Gauss 分布である場合について導出されている厳密解と、 あるパラメータ空間における中立安定性が良く一致することが確かめられた。 しかし異なる速度分布の場合について調べてみると、 通常の近似的解釈で扱われるように安定性は速度分散のみで決定されす、 分布の 4 次のモーメントの大きさに依存して安定性の境界値と不安定化されるモードが大きく変化することが示された。 このことは、 例えば宇宙論における銀河の渦巻きパターンが密度波理論により解釈される際の、 特徴的なパターンについて新しい解釈を与える可能性を示唆している。

1

はじめ [こディスク状に分布した重力相互作用多体系は天文分野ではしばしば ‘stellar disk’ などとも呼ばれ、 平面状に降着した天体としての銀河の一つの理想的なモデルと見なされる。実際の宇宙論的な銀河は通常 stellar disk の構成粒子としての恒星成分と、 ガス成分からなる 2 成分系として扱われる。 後者の成分比が大きい系では局所平衡の成り立つ流体系として取り扱うことが可能であるのに対して、 恒星系の場合には粒子同士が無衝突であるために様々な奇妙な性質を持ち長く興味を集めてきた。特にディスク銀河にしばしば観測される渦巻状模様に関してらが流体モデルを用いて線形波の分散関係を解析的に求めてより [1]、 密度波理論としての解釈が進められてきた。 近年では
Kondoh らの非線形解析によりソリトン状パターンの形成をも予言されている [2]。 ただしこれらの解析が厳密に可能なのは、 特殊な問題を除いて速度分散が 0 の’cold disk’$\mathrm{L}\mathrm{i}\mathrm{n},\mathrm{S}\mathrm{h}\mathrm{u}$187-方で定常な構造をもつ実際の天体では有限の速と呼ばれる極限に限定されている度分散をもつことが本質的に不可欠な要素であると考えられる。 これは重力のみを考えた場合常に微小な揺らぎに対して不安定で一方的に密度の集中が進むのに対し、 速$\mathrm{c}$’’

度分散による拡散的な振る舞いとの間で 静水圧 的な力学平衡が成り立つためと説明できるが、 この効果を入れると後述する理由で、 流体的な記述が途端に困難になる。
ダイ
従ってこの系の安定性については、 それ自体が未解決な問題であるだけでなナミクスに関する理解とも密接に関わっている。の時間発展方程式とし連続極限においてこの系は、 一粒子分布で代表されるての無衝突 Bortzmann 方程式 (Vlasv 方程式)$\text{く_{}\urcorner}$$f(\mathrm{x}, \mathrm{v})$(1),$\frac{\partial}{\partial t}f+\mathrm{v}\cdot\nabla f-\nabla\Phi\cdot\frac{\partial}{\partial \mathrm{v}}f=0$として記述されることが一般的に知られている。 ここで は自己重力により作り出すポテンシャルなので Poisson 方程式により 2\Phi (x) $=4 \pi G\int fd\mathrm{v}$ と書ける。 このよう$\Phi$に位相空間における分布の時間発展が分かつてはいるもののダイナミクスに関する理解は必すしも容易ではない。
例えば、定常な分布に fig. 1 左のような摂動を与えた初期条件からの時間発展は右図.のように分布が引き伸ばされ一方的に位相が混ざる方向に進む (phase mixing)。 このことは、 熱力学的な緩和によりローカルな速度分布が Maxwell 分布に回復するといった他の系のようなメカニズムがはたらかない無衝突系としての特徴を端的に示していて、このため 空間について平均した上で実空間でのダイナミクスが非自明なものになって$v$いる。 実際に van Kampen の証明によると、 三次元等方的な系における波は$v$

空間でのmixing 効果により常に減衰し、 逆に音波的に減衰無く純実数の周波数で伝搬するモードを与えようとすると速度分布に特異点ができてしまうことが示されている [3].

Figure 1: Vlasov 系における phase mixing の模式図

188

同様の困難さは、 (1) 式から流体方程式系のモデルを構築する際にも現われる。 phasemixing 効果は、 言い換えれば の各次数のモーメントがダイナミクスに関して自由度を持つことを示しているため、 閉じた方程式系として記述できない (closure problem)。
$v$

これを回避する手法として、 回転の強い系など平均運動に比べて速度分散が小さいような状況を仮定した上で高次のモーメントの効果を切るという近似を用いて、 よりよい精度で記述する閉じた流体方程式系のモデルを考案する試みは多くなされてきてい6$[4, 5]_{\text{。}}$

回転ディスク系の安定性に関する解析としては、 Toomre により波長 及び系のローは順に平均回転運カルな状態を示す無次元量 $Q\equiv 2\Omega c/\pi G\sigma(G$ は重力定数、動の角周波数、 速度分散、 面密度) の 2 変数空間における中立安定線を非摂動の速度分布が Gaussian である場合に限定して (1) から直接、厳密な解として導出された。 しかし改めて熱力学的系ではないことを考慮すると (1) の定常解として Gaussian からかけ離れた分布を与えることも可能であり、 –般的に解が示されているとは言い難い。そこで、 上述の流体モデル方程式との対応としても安定性及び波のモードが議論され、 最も簡略には断熱流体近似 ( が定数) による 2 次形式の中立線が Toomre と比較的良い一致を見せるため今日でもしばしば引用される。 また Hunter は自身の導出による流体モデルをもとに、 4 次以上のモーメントが無視できるような分布では Toomre の解析値から有意にすれるとし、 その場合の安定線を示している。 このように、 速度分布を変化させた場合それによる安定性 波の分散関係等への寄与は未だ定量的には理解されて$\lambda$$\Omega,$$c,$$\sigma$$c$$=$

いない。

同様の問題に関する別の手法として、 Ka 石 により導入された 1 粒子軌道を基にしこれを用いた後続研究により、 van kampen が指摘した非減た行列法が挙げられる
衰モードが回転ディスクにおいては存在し得ることを示すなど部分的な解明が進んだ
ものの [7]、 全体像を明らかにするには至っていない。
このように解析的手法ではこの系特有の問題が付きまとうため数値計算による解明を試みることが効率的だとも考えられるが、 長距離相互作用系のシミュレーションにオーダーの計算が必要になるため一はナイーブには前粒子間の相互作用を考慮した般的には膨大な時間を要する。 近年の計算機の能力の進展と専用計算機等の出現によ} 万体レベルの計算も可能になっているものの、 ここで問題にしている
り今日では分散関係や安定性について系統的にサーチするという使い方は容易な計算量とはいえない。 そこで我々は線形問題に関する一般性を壊さないようなモデルとして、 棒状粒子を用いたシミュレーションにより大幅に計算量を効率化することで、 一般の分布に対する線形波のモード・中立安定線等を定量的に求めることを試みる。
$\mathrm{a}\mathrm{j}\mathrm{s}$$[6]_{0}$$N^{2}$$1\alpha$

1892 計算手法
ここでは、 任意の空間的な分布をもつ回転デイスク上の局所的な状態に対する分散関係・安定性等を議論するためのモデルを考える。 このため、 扱う系の定常状態は–様な空間分布で剛体回転しているとする。 さらに密度揺らぎの自由度を 軸の 1 次元的に制約し摂動状態 $f_{1}(x, u, v)$ の時間発展を追うことで計算の効率化を図る。 今興味のある対象の様な線形問題として捉えられる領域では 1 、これに直行する方向の揺らぎの時間発展に影響することが無いため、 この制約による物理的な影響は無の有無が$x$$f_{1}$

いと考えられる。

この目的に合うモデルを実装するために通常の質点系の代わりに、 fig 2 に示したような大さを持たず無限に長い棒状粒子を平行に並べた系を用いる。 さらにディスク面と平行な $x,$ 軸は 軸回りの系の回転に固定されていて, このため回転の効果を各粒子はコ$z$$y$

リオリカとして感じる。 -従っ$\circ$てこのコリオリカを伝えるために、 1 粒子の自由度としてFigure 2: 棒状粒子モデル
方向の速度 以外に潜在的に 方向の速度 を持つ必要がある。 これらを踏まえると、 番目の粒子が受ける力は粒子の線密度を $m$ と座標 x、

$u$

$x$

$v$

$y$

$i$

しで、
$F_{ix}$$=$$F_{iy}$$=$$\sum_{j\neq i}Gm\frac{x_{J}-x_{i}}{(x_{j}-x_{i})^{2}}+2\Omega v_{\mathrm{i}}$$-2\Omega u_{i}$,(2)(3)と書き表される。 系の 方向については周期的境界条件を適用した。 粒子数は $N=$$960,000$ を用いて計算を行ったが、 この粒子の多さはおもに速度分布を局所的にも正確に再現するためなので力の計算は系の 方向を 2, 400 に区切ったメツシュの分布をフー$x$$x$リエ変換することで近似的に求めた。

安定領域においては、 十分微小な摂動を与えて分散関係を計算する。 また、 不安定領域においても揺らきの成長の立ち上がりに注目する。
$\iota$1903 計算結果
ます-」 安定な領域で任意の微小の摂動を与えた初期条件からの時間発展を観察する。

この摂動分布に対して時間・空間 $(t-x)$ の二次元に関するフーリエ変換を行うことで、
波数.‘ 周波数 $(k-\omega)$ 空間におけるスペクトルとして分散関係を数値的に得る (fig 3)。
以前の章で述べたとおり、この系で一般的には phasemixing による振幅の減衰が強く現れることが考えられる。
その場合上記のような実数の$k,$$\omega$を用いて $\exp(ikx-i\omega t)$により展開すると、 スペクト
ルは連続的になってしまいこの手法による解析は困難であることが考えられる。 しかしこの系の今必要とするパラメータ領域では図に示した通り十分に鋭いピークが現れてぃて、 振動数 [。比、て減衰の$\mathrm{B}\frac{\simeq}{\mathrm{a}}3\Phi$$\in \mathrm{a}$$<\in$$<$omega

Figure 3: 周波数空間における振動モードのピーク $(Q=1.2)$

時定数の十分長いモードであることが伺える。

次に系の状態 $Q$ を減少させて、 徐々に最も周波数の低いモードを 0 点に近づけることで、 不安定化を起こす境界を見出すことを試みる。 具体的には同一の波数のモードと不安定側での成長率を内挿して 0点を決定する。 fig 3 は安定・不安定領域を分ける中立安定曲線を、 無次元化した波数$(c^{2}/G\rho)k$ 及び非摂動の系の状態を表すパラメータ $Q\equiv 2\Omega c/\pi G\sigma$ に対して描いてある。 曲線より下が不安定領域であり、 ある系の状態に対応する値における水平線が全ての波数においてこの曲線を横切ることがなければ安定であることを示している。離散的なプロットが我々の計算値を示す。 先す速度分布が Gaussian のケースでは Toomre
について、 安定領域側での周波数の二乗$\omega^{2}$$(i\gamma)^{2}$$Q$の解析による線$\frac{1}{2}\kappa=1-\exp(-\kappa^{2}Q^{-2})I_{0}(-\kappa^{2}Q^{-2})$($I_{0}$(4)は変形ベツセル関数) と精緻に一致している。 これにより我々の棒状粒子による計算が 2 次元的ディスク系の線形波に関する正確なモデルとなっていることも検証された。併記した点線は速度分散を定数と仮定した最も簡略な近似 (断熱近似) による 2 次式の分散関係式$\omega^{2}=\kappa^{2}-2\kappa+Q^{2}$の左辺を 0 とおいた曲線を示す。

(5)

131

$[mathring]_{6}^{1}$$\mathrm{k}$Figure 4: 中立安定線の解析値をよび数値計算値

次に本題となる速度分布に対する依存性を見るために、 一例として Gaussian にカットオフをいれてモーメント比 $<v^{4}>/<v^{2}>^{2}$ を 30 から 23 に減少させた分布について同様に中立安定性を求めた (同図) 。これが Gaussian のそれに比べて低い側にずれたことは、 高次のモーメントを切った近似による Hunter の導出式
$Q^{2}=0.5(-\kappa^{2}+2\kappa+\sqrt{-2\kappa^{4}-2\kappa^{3}+4\kappa^{2}})$.

(6)も分布によっては妥当であることを示唆している。

4

まとめ

本研究では効率的にモデル化することによって、 前例の少ない重力多体回転デイスク系の分散関係および安定性を数値的に求めることに成功した。 これにより理解された点は、$\bullet$Gauss 分布に限定すると中立安定曲線は Toomre のよる解析値と精緻に一致する。
この解析解は Vlasov 方程式からダイレクトに導出されたものなのでこの結果は妥当と考えられる。

上記の曲線と、 広く用いられている最低次の近似による 2 次曲線は比較的近い値を示している。 これは偶然であると考えられるものの、 この事実のために安定性は速度分布の違いに対して強くは依存しないと思われている節もある。
$\circ$モーメント比を Gaussian からすらした分布では中立安定線が有意に変化する。

192

・上の結果のより低次のモーメントから展開する方法による Hunter の近似は Gaussianに対しては誤っているように見えるものの、 高次のモーメントが小さな分布に対しては良い近似になっていることが示唆された。
これらを踏まえると、 この系の安定性の速度分布に関する情報として単に速度分散$c$のみで–義的に決まらす- 分布の形状に依存する。 ただし特に Gauss 分布の周辺では、個別の分布の詳細までは立ち入らすにいくつかの低次のモーメントの関数として表せることが期待される。 我々は追加的な計算によりこの関数形を数値的に求める事を引き続き試みている。 一方で観測論との対応としては大きな変動を経た結果として定常に存在するディスクが持つ分布にどの程度の幅があるかという点が重要になるかもしれない。

参考文献. Lin and. Shu. Astrophys. J., 140:646, 1964.

[11$\mathrm{C}.\mathrm{C}$[21S. Kondoh. R. Teramoto, and Z. Yoshida. Phys. Rev.

$\mathrm{F}.\mathrm{H}$, 61:5710, 2000.

$E$[31 N. G. van Kampen. Physica. 21:949, 1955.
[41C. Hunter. Stgtdies Appl. Math., 49:59, 1970.

[5] P. Amendt and P. Cuddeford. Astrophys. J,. 368:79, 1991.
[61 A. J. Kalnajs. Astrophys. J., 212:637, 1977.
$\lfloor 7]$C. Hunter. In Disks of Galaxies: Kinemtics, Dynamics and Perturbations, $p\prime vceedings$ of the 4thGuillermo Haro Conference, 2001, page 293. Astro. Soc. Pacific, 2002.