5.回転円盤の不安定性

Note

2023/05/19は、前半で乱流粘性拡散方程式の導出、後半で円盤の自己重力不安定性の話を行います。 授業ノートはこのリンクのgoogle jamboardを使います。 zoomは通常と同じものを使います。

この章でやること

• 宇宙における回転円盤としての銀河、原始惑星系円盤を知る。
• 降着円盤の観測的制約を知る
• 降着円盤と粘性の関係や磁気回転不安定性(MRI)について概要を知る
• いくつかの不安定性(自己重力不安定)について学ぶ

降着円盤を学ぶモチベーションのおさらい

• X線連星のように、定常的にX線を出すような高エネルギー現象を説明したい (円盤のエネルギー放出機構)
• 原始惑星系円盤のように、ある程度の時間で消失する現象を説明したい （円盤の散逸）
• 増減光を説明したい (円盤の不安定性)

前半の問題設定

質量\(M_{\rm star}\)の中心星の周りの回転円盤の時間発展を記述する式を導出することを目標にする。座標系は円筒座標\((r,\theta,z)\)を取る。面密度\(\Sigma\)で、粘性応力が働いている円盤を仮定している。流体の基礎方程式を前提知識として始める。

降着円盤の基礎方程式 - 特に角運動量輸送

ここは前回授業でやっていた場合は飛ばします

さて、降着円盤の基礎方程式を考えてみよう。 まず、流体力学のナビエ・ストークス方程式を思い出しておこう（どうしても思い出せない場合は、A1.ナビエ・ストークス方程式をクリック）。

連続の式

\[\frac{\partial \Sigma}{\partial t}+\frac{1}{r} \frac{\partial}{\partial r}\left(r \Sigma v_r\right)=0\]

導出はナビエ・ストークス方程式のページを参照

角運動量進化の式

円盤内にある薄い円環を考えよう。その円環の単位質量あたりの角運動量は、\(r^2\Omega\)となる。 円環の面積は\(2\pi r \Delta r\)、面密度は\(\Sigma\)であることに注意する。 円環の角運動量のラグランジュ微分は、円筒座標のラグランジュ微分に注意して,

\[\frac{\partial}{\partial t}\left(\Sigma r^2 \Omega\right)+\frac{1}{r} \frac{\partial}{\partial r}\left(\Sigma r^3 \Omega v_r\right)=\mathscr{G}\]

ここで、トルクに関わる項は\(\mathscr{G}\)とおいた(花文字のGです)。 トルクとは、なんらかの力を介してよその円環等と角運動量の遣り取りをする力だと思おう。 この右辺の\(\mathscr{G}\)がどういう物理的意味を持つかどうかを考え、その中身を書いてみよう。

式をよく見てみると、左辺第一項は、角運動量の時間進化を表していそうである（角運動量は、(質量)×(距離)\(^2\)×(角速度))。 そこで、とある細い円環を考えて、その角運動量輸送を考えよう。 \(r\)から\(r+dr\)の範囲内の円環の面積は\(2\pi rdr\)なので、この円環の角運動量は\(\Sigma r^2 \Omega2\pi rdr\)である。 つまり、先の式の両辺に\(2\pi rdr\)をかけたら、細い円環の角運動量輸送の式になりそうである。 その時の右辺は、粘性トルク\(G(r)\)を導入することで

\[\mathscr{G} \cdot 2 \pi r d r=G(r+d r)-G(r)\]

と書けるだろう（だいぶずるい感じがするが、こうすると物理の意味を考えながら立式できて良い）。 すなわち、いま書き下したかった\(\mathscr{G}\)は、粘性トルク\(G(r)\)を用いて

\[\mathscr{G}=\frac{1}{2 \pi r} \frac{d G}{d r}\]

とかけそうである。そう思って\(G(r)\)を書いてみよう。

回転円盤における粘性トルクの形を予想してみる

ここでは、速度差があったときに、その速度差を小さくしようと働く力である。 粘性係数\(\mu\)とシアストレスの関係は、ミクロなダストの抵抗則 7.ガス円盤中の固体微粒子の運動/#-stokes の導出時に出てくるが、

\[T_{x}=\mu \left|\frac{dv_x(y)}{dy}\right|\]

という形で関係していることを思い出しておこう。

降着円盤で円環に働く速度シアはどんなものだろうか？ おそらく最も強いものは、円環の回転運動の速度シアだろう。 （隣の円環と回転速度に差があり、かつ粘性があれば、回転速度のシアに応じて粘性トルクという力が発生しそう、ということ）

すなわち、ちょっと回転速度の\(r\)微分を考えてみれば良さそうである。\(v_\theta\)の\(r\)微分は

\[\frac{d v_\theta}{d r}=\frac{d}{d r}(r \Omega)=\Omega+r \frac{d \Omega}{d r}\]

となる。速度シアは第二項で効いていそうである。 ここから、粘性ストレステンソルを予測してみると、その形は

\[(\text{粘性ストレステンソル})=\mu \times \text{(速度シア)}=\mu r \frac{d \Omega}{d r}\]

だと予想される。トルクとは\(r\times F\)だったのを思い出して一発\(r\)をかけてから、\(\theta\)方向及び\(z\)方向に積分してみると

\[G(r)=\int r d\theta \int dz \mu r^2 \frac{d\Omega}{dr} = 2 \pi \nu \Sigma r^3 \frac{d\Omega}{dr}\]

これで、目的だった\(\mathscr{G}\)が導出できたので、最初の式に代入して

\[\frac{\partial}{\partial t}\left(\Sigma r^2 \Omega\right)+\frac{1}{r} \frac{\partial}{\partial r}\left(\Sigma r^3 \Omega v_r\right)= \frac{1}{r}\frac{\partial}{\partial r}\left(\nu \Sigma r^3 \frac{d\Omega}{dr}\right)\]

これが、降着円盤において、角運動量が速度シアに起因する粘性ストレスで輸送されるときの方程式である。

ケプラー回転円盤における角運動量進化の式

このままだとまだわからんので、とりあえずケプラー回転を仮定して、もう少し式を整理してみよう。 ケプラー回転の場合、角速度は

\[\Omega = \left( \frac{GM}{r^3}\right)^{1/2}\]

だったので、例えばこれの\(r\)微分は

\[\frac{\partial \Omega}{\partial r} = -\frac{3}{2} (GM)^{1/2} r^{-5/2}= -\frac{3}{2} \frac{\Omega}{r} \]

などとかけそうなので、これと連続の式を用いて、\(v_r\)を消すように努力すると、

\[\frac{\partial \Sigma}{\partial t}=\frac{3}{r}\frac{\partial}{\partial r}\left[r^{1/2}\frac{\partial}{\partial r}(\nu\Sigma r^{1/2})\right]\]

と、ケプラー回転を仮定してもいい場合の降着円盤の面密度の時間進化を粘性を用いて書き表すことができた。

なお、この式の形は、拡散方程式と同じ形になっていることを確認しておこう。拡散方程式は、物理量の時間微分が、物理量の空間二階微分に比例するという形であった。

時間発展の可視化

Lynden-Bell & Pringle 1974

降着のタイムスケール

ここで、導出した式を用いて、ざっくりと円盤の降着のタイムスケールを見積もってみよう。 両辺のオーダー見積もりをすると

\[\frac{\Sigma}{\tau} \sim \frac{\nu \Sigma}{\Delta r^2}\]

よって

\[\tau \sim \frac{\Delta r^2}{\nu}\]

すなわち、粘性が大きいほど短い時間で円盤が拡散するという結果が得られた。

ところで粘性トルクがなかったらどうなるの

ここまでだるい計算をして粘性降着円盤の基礎を見てきたが、そもそも粘性トルクが仮になかったら、どうなるのか。 最後の式に\(\nu=0\)を代入して

\[\frac{\partial \Sigma}{\partial t}=0\]

すなわち、面密度は全く時間発展しない。 これは、もし粘性がなければ、すなわち速度シアがあっても互いに力のやり取りをしなければ、角運動量は輸送されず、結果として降着円盤は初期状態をひたすら保つ、ということを示している。

つまり、粘性が必要だよね、という結論に至る。

降着円盤におけるガスの粘性

粘性とは、(長さ)×(速度)の次元を持つ、流体のねばりけを表す量である。

分子粘性

分子粘性とは、流体の分子の運動による粘性で、ガスの平均自由行程\(\lambda_{\rm mfp}\)とガスの音速\(c_s\)を用いて下記のように表される。

\[\nu_m= \lambda_{\rm mfp} c_s\]

分子粘性の値を自分でざっくり見積もってください。更に、分子粘性によって円盤の角運動量輸送が決まっている場合の、円盤の降着タイムスケールを計算してください（答えなし）

円盤のグローバルな粘性

降着円盤においてガス降着に必要な粘性は分子粘性では足りない。 円盤において乱流渦があったとする。乱流渦の大きさを\(L\)、乱流の速度を\(V\)とすれば、粘性は

\[\nu= LV\]

のオーダーを持つはずである。

乱流の大きさや速度はもちろん円盤の性質や場所によって変わるだろうがなにかうまくスケールする方法はないだろうか？ 乱流は渦の最大長はせいぜいスケールハイトを超えないだろう。 また、乱流の速度は音速に比例すると思う。 すると、円盤のグローバルな粘性は円盤のスケールハイト\(h\)と音速\(c_s\)を用いて、無次元量\(\alpha\)を導入し、

\[\nu= \alpha h c_s\]

とかけそうである。

乱流粘性の値を自分でざっくり見積もってください。更に、乱流粘性によって円盤の角運動量輸送が決まっている場合の、円盤の降着タイムスケールを計算してください（答えなし）

ここで、\(\alpha\)はどう分配してもいいのだが、よく使われるのは、長さを\(\sqrt{\alpha}h\)、速度を\(\sqrt{\alpha}c_s\)とする記述である。特にガスの乱流速度を　\(\sqrt{\alpha}c_s\) と記述するのは今後も出てくるので覚えておいてほしい。

(2023-05-19 前半ここまで)

後半の問題設定

回転円盤には不安定性が存在する。ここでは特に、円盤の重力による不安定性について紹介し、どのような条件で不安定になるかを議論する。

自己重力不安定

別のノートに書いてしまったので、A2.線形解析を参照。

回転円盤の重力不安定(Toomre's Q)

無限に薄い回転円盤が重力的に不安定になる条件を、線形解析から見積もってみよう。 線形解析が何かわからない人は、このページ(A2.線形解析)に重力不安定性（ジーンズ不安定性）の例をあげて解説しておいたので見てみてほしい。 知ってる人は見なくて大丈夫。

本解析から出てくる条件はToomreのQ値と呼ばれ、例えば回転する銀河において自己重力によって星ができるかを判定したり、重い原始惑星系円盤において自己重力によってガス惑星ができるかどうかを判定することができる。

基礎方程式をたてよう。 円筒座標系\((r,\theta,z)\)をセットし、無限に薄い円盤を考えて、密度\(\rho\)は\(z\)積分して面密度\(\Sigma\)を使おう。 状態方程式は\(p\propto \rho^\gamma\)を使おう。定数はどうでもいいから\(p=\Sigma^\gamma\)でいい。

全ての基礎方程式を書くのは大変なので、線形化した方程式を書こう。 \(x=x_0+x_1\)と書き、\(x_1\)は微小量とする。 面密度\(\Sigma\),動径方向の速度\(v_r\),回転方向の速度\(v_\theta\)についてそれぞれ微小量の二次以降を無視すると

連続の式

\[\frac{\partial \Sigma_1}{\partial t}+\frac{1}{r}\frac{\partial \Sigma}{\partial r}(r \Sigma_0 v_{r,1}) + \Omega \frac{\partial \Sigma_1}{\partial \theta} + \frac{\Sigma_0}{r} \frac{\partial v_{\theta,1}}{\partial \theta} =0\]

導出

円筒座標系における連続の式をまず思い浮かべて

\[\frac{\partial \rho}{\partial t}+\frac{1}{r}\frac{\partial}{\partial r}(r \rho v_r) + \frac{1}{r}\frac{\partial}{\partial \theta}(\rho v_\theta) + \frac{\partial}{\partial z} (\rho v_z)=0\]

両辺を\(z\)積分して\(\Sigma\)に変換しつつ眺める。 無限に薄い回転円盤なので\(v_z=0\)は良さそう。 第3項だけちょっとだるそうなので中を見ておくと、\(v_{\theta}=\Omega r\)であることに注意して

\[ \begin{align} \text{(第3項)}&=&\frac{1}{r}\frac{\partial}{\partial \theta}(\rho v_\theta)\\ &=&\Omega \frac{\partial \rho}{\partial \theta} + \frac{\rho}{r}\frac{\partial v_\theta}{\partial \theta} \end{align} \]

あたりに注意しながら、代入して二次の項を消していけばでそう

\(r\)方向の運動方程式

ここで、\(dh=dp/\Sigma\) と導入して $$ \frac{\partial v_{r,1}}{\partial t}+\Omega\frac{\partial v_{r,1}}{\partial\theta}+2\Omega v_{\theta,1} =- \frac{\partial h_1}{\partial r} - \frac{\partial \Phi_1}{\partial r} $$

導出

驚きのだるさ。A1.ナビエ・ストークス方程式を眺めながら1タームずつ確認しよう。粘性は\(\nu=０\)としよう。 $$ \frac{\partial v_r}{\partial t}+v_r \frac{\partial v_r}{\partial r}+\frac{v_\theta}{r} \frac{\partial v_r}{\partial\theta}+v_z \frac{\partial v_r}{\partial z}-\frac{v_\theta^2}{r} =-\frac{1}{\rho} \frac{\partial p}{\partial r} +F_r, $$ 外力は重力を入れておけばよかろう。\(F_r=-\nabla \Phi\)だと思うと、 $$ \frac{\partial v_{r,1}}{\partial t}+\Omega\frac{\partial v_{r,1}}{\partial\theta}+2\Omega v_{\theta,1} =-\frac{1}{\rho} \frac{\partial p_1}{\partial r} -\frac{\partial \Phi_1}{\partial r} $$ ここでちょっと面倒なのが\(\rho\)の取り扱い。いま、\(z\)方向を積分して使いたいので、\(dh=dp/\rho\)という量(エンタルピー)を導入した。

\(\theta\)方向の運動方程式

エピサイクル振動数を $$ \kappa=\sqrt{ r \frac{\partial \Omega^2}{\partial r} + 4\Omega^2} $$ と定義すると

\[ \frac{\partial v_{\theta,1}}{\partial t}+\Omega \frac{\partial v_{\theta,1}}{\partial \theta}+ \frac{\kappa^2}{2\Omega}v_{r,1} =- \frac{1}{r}\frac{\partial h_1}{\partial \theta} - \frac{1}{r}\frac{\partial \Phi_1}{\partial \theta}, \]

導出

とりあえずA1.ナビエ・ストークス方程式　の\(\theta\) 方向について粘性を無視して見てみると $$ \frac{\partial v_\theta}{\partial t}+v_r \frac{\partial v_\theta}{\partial r}+\frac{v_\theta}{r} \frac{\partial v_\theta}{\partial \theta}+v_z \frac{\partial v_\theta}{\partial z}+\frac{v_r v_\theta}{r}=-\frac{1}{\rho r} \frac{\partial p}{\partial \theta} +F_\theta $$ 注意すべきはたぶん第二項で

\[ \begin{aligned} v_r \frac{\partial v_\theta}{\partial r}&=& v_r \frac{\partial (\Omega r)}{\partial r}\\ &=&r v_r \frac{\partial \Omega}{\partial r} + v_r \Omega \end{aligned} \]

すると、第５項も$ v_r \Omega$なので

\[(第2項)+(第5項)=r v_r \frac{\partial \Omega}{\partial r} + 2v_r \Omega\]

ところで、急に天下り的に\(\frac{d \Omega^2}{dr}\)を計算すると

\[\frac{d \Omega^2}{dr}=2\Omega \frac{d \Omega}{dr}\]

なので

\[ \begin{aligned} (第2項)+(第5項)&=&\frac{r v_r}{2\Omega} \frac{\partial \Omega^2}{\partial r} + 2v_r \Omega\\ &=&\frac{v_r}{2\Omega} \left( r \frac{\partial \Omega^2}{\partial r} + 4\Omega^2\right) \end{aligned} \]

で、このカッコの中がエピサイクル振動数の二乗。という無理やりな式変形を頭に留めておく。 まあ、ご想像の通り\(r\)方向の運動方程式とにた形にしたかっただけ。

分散関係の導出

さて、摂動を与える方向として、方位角方向\(\theta\)を考える。 そこで、\(x=x_a \exp[i(m\theta-\omega t)]\)という波の形を考えよう。

（まだここ導出途中だが、もう代入するだけ）

分散関係式は

\[(m\Omega-\omega)^2 = k^2 c_s^2 - 2\pi G\Sigma_0|k| + \kappa^2\]

ここで、\(m=0\)、すなわち軸対称モードを考えよう。すると、

\[\omega^2 = k^2 c_s^2 - 2\pi G\Sigma_0|k| + \kappa^2\]

ちょっと式を整理すると

\[\omega^2 = c_s^2(k - \frac{\pi G\Sigma_0}{c_s^2})^2 + \kappa^2 - \left(\frac{\pi G\Sigma_0}{c_s}\right)^2\]

さて、とある\(k\)にたいして\(\omega^2<0\)になれば、不安定な解が存在するということだった。 下に凸な二次曲線がマイナス側にいくには、最後の定数項が負になれば不安定。すなわち

\[\kappa^2 - \left(\frac{\pi G\Sigma_0}{c_s}\right)^2 < 0\]

なら不安定となる。整理すると、

\[\frac{\kappa c_s}{\pi G\Sigma_0} < 1\]

であれば不安定な解が存在する。

この左辺を通常\(Q\)と置き、Toomreの\(Q\)値と呼ぶ。

分散関係をプロット

Toomre's Qの偉いところ

線形解析は大変ではあるが、結果を非常にシンプルなパラメータQに落とし込み、誰でも重力不安定性かどうかを判定可能にしたところが偉い。 みなさんが理論計算をするときは、ぜひ観測の人や分野の違う人でもシンプルに使える条件式に落とし込めると良いと思う。

Toomre's Qの具体例 - 原始惑星系円盤の場合

原始惑星系円盤は、分子雲コアが重力収縮する中で、初期に持っていた角運動量によって円盤状の天体が形成されることで誕生すると考えられている。 原始惑星系円盤は形成初期の段階では、おそらく最も質量が重いので、重力不安定が作用してもおかしくないので、形成初期の円盤の進化を考えてみよう。

初期に非常に重たい円盤が形成されたとしよう。 その円盤が\(Q<1\)を満たすのであれば、おそらく巨大ガス惑星が随所に形成されると考えられる。 ただし、これらの惑星は早いタイムスケールで中心性に落下してしまうと考えられている。

そこまで重くない円盤で\(1<Q<2\)程度であった場合は、自己重力によってクランプ状の構造はできないまでも、渦状腕（スパイラルアーム）が形成されることがわかっている。 この渦状腕を通して角運動量輸送が発生するため、このように重たい円盤を降着円盤に見立てて、粘性降着の\(\alpha\)を見積もる試みもされてきた。

参考: Takahashi et al. 2013, Zhu et al. 2010

以上のように、自己重力が効くほど重たい円盤は、自己重力不安定によってクランプや渦状腕を作って質量を内側に落とす。 すなわち、円盤形成初期においては重たすぎる円盤をself-regulateして質量を輸送するメカニズムとして機能していると考えられる。 この妥当性を検討してみるため、原始惑星系円盤の質量を考えてみよう。

観測結果から考える原始惑星系円盤の重力不安定

円盤全体の質量から重力不安定を考えてみる。 まず、重力不安定を円盤質量から考えやすくするため、Toomre's Qの式を少しいじってみよう。 円盤の質量は、仮に面密度が一定であれば\(M_{\rm disk}=\pi r^2 \Sigma_0\)くらいである。 また、円盤のスケールハイトは\(h=c_s/\Omega\)である。 ここで、\(\Sigma_0\)もスケールハイトも円盤内の適当な場所での値であることに注意。 エピサイクル振動数は\(\Omega\)くらいだと思い、\(\Omega^2=GM/r^3\)であることに注意すると、分母分子に\(1/\Omega^2\)をかけて

\[Q=\frac{\kappa c_s}{\pi G\Sigma_0} \approx \frac{\Omega c_s /\Omega^2}{\pi G\Sigma_0 \times r^3/GM_{\rm star}} \approx \frac{M_{\rm star}}{M_{\rm disk}} \frac{h}{r}\]

よって不安定になる条件\(Q<1\)は

\[ \frac{M_{\rm disk}}{M_{\rm star}}>\frac{h}{r}\]

と書き換えられる。典型的なアスペクト比\((h/r) \sim 0.1\)くらいなので、円盤質量が中心星質量の0.1倍くらいより大きければ重力不安定になりそうな気がしてくる。

適当な円盤の適当な場所を想定し、スケールハイトを見積もってみよう。例えば、中心星が太陽質量、中心から30 auで、温度が20 Kのときはアスペクト比はいくつになるだろうか（答えなし）

円盤質量の観測的制限は案外難しい。 円盤の質量の大半を占めるのは水素ガスであるが、中性水素は高温領域でしか効率的に輝線を出さない。 そのため、COなどの分子輝線を用いたガス面密度の推定がなされるが、COとH2の存在量比の不定性が大きい。

ここでは、ダストと呼ばれる固体成分のミリ波観測を用いる。 円盤内の固体成分はガスの1/100程度であると考えられているが、不定性は大きい。 また、固体成分が微惑星のように大きな天体になった場合、それらはミリ波では観測できない。 その一方で、ミリ波で見ると原始惑星系円盤の質量を担う外側領域は光学的に薄いことが多く、質量測定に使える。 更に、ミリ波の連続波は輝線より感度高く観測できるため、サーベイ観測が進んでおり、全体像を掴むには良い観測手法となっている。 これらの不定性を頭に入れつつ、ダストのミリ波観測を見てみよう。

Tychoniec et al. 2018

上の図は、複数の若い星形成領域におけるダスト質量の統計分布である。 Class 0から1,2と進むにつれてダスト質量は減っていることがわかる。 Class 0 天体は中央値250地球質量とわかっている。

原始惑星系円盤の質量はどれくらいより重いと重力不安定になりうるだろうか。 先程の見積もりから、ガス質量が0.1太陽質量を超えたあたりから不安定になると考えられる。 仮にガスがダストの100倍の質量あるとすると、ダストが1e-3太陽質量=300地球質量より重いと不安定になりそうな感じがする。

先程の図を見ると、Class 0の中央値(250地球質量)よりも少し重い天体はすべて重力不安定になっていそうに見える。 だが、ここまでの見積もりには多くの不定性が含まれており、例えばこれらの天体は固体が濃集してダスト/ガス比が0.1くらいなのかもしれない。

参考: Kratter and Lodato 2016, arXiv, ARAA

磁気回転不安定性(Magneto-Rotational Instability, MRI)

ここは次回授業範囲になるそうなのでやりません

参考 Youtube video by Xuening Bai

粘性降着円盤を考える上でのMRIの必要性

降着円盤のShakura-Sunyaevモデルでは、速度シアに対して粘性が働く必要があった。 この粘性の物理的理解は脇において、Shakura-Sunyaevモデルでは粘性の大きさを無次元パラメータ\(\alpha\)をもちい、スケールハイト\(h\)と音速\(c_s\)を用いて

\[ \nu = \alpha h c_s \]

と表す。

ではここで、この粘性の大きさはどのくらい必要かを考えてみよう。 まず最初に思いつく粘性の起源は分子粘性である（が、これだと全然足りないことを今から見ていく）。 \(\(\nu_m= \lambda c_s\)\)

\[\lambda = \frac{1}{n \sigma_{mol}}\]

ここで、\(n=10^{12}{\rm~cm^{-3}}\)、\(\sigma_{mol}= 2\times10^{-15}{\rm~cm^2}\)とおくと

\[\nu_m \sim 2.5 \times 10^7 {\rm [cm^2~s^{-1}]}\]

粘性タイムスケールは

\[t_{\nu}\sim\frac{r^2}{\nu}=3\times 10^{13} {\rm~[years]}\]

原始惑星系円盤の散逸に宇宙年齢を待つわけには行かない。 すなわち、分子粘性は降着円盤の粘性を説明するのには不適当である。

円盤の粘性起源としての磁気回転不安定性による乱流

磁気回転不安定性についてかんたんに触れておこう。 磁気回転不安定性(Magneto-rotational instability, MRI)はBalbus and Hawley 1991によって提唱された不安定性で、差動回転する円盤を鉛直磁場が貫いていた場合、不安定性が発展して乱流状態になるものである。

*Velikov (1959) によって発見されたもので、Chandrasekhar (1961)の本にも記載されているらしい

下の図は、鉛直方向に貫いた磁場が時間発展する様子を表す。

Balbus and Hawley 1998

MRIの分散関係(磁気回転不安定性が起こる条件)

線形解析は、Balbus and Hawley 1998 の"IV. THE LINEAR STABILITY OF MAGNETIZED ACCRETION DISKS B. Weak-field shearing instability: a simple treatment" というところを追うとわかる。

内容を大幅に省略すると、円筒座標系において、z軸対称なガス円盤を考える。 磁場は非常に弱いとする（円盤の力学的平衡状態に影響を与えない）。 摂動はz方向に\(e^{ikz}\)で与える。 非圧縮性流体を考える。 すると、分散関係は、式(111)のとおり

\[ \omega^4-\omega^2\left[\kappa^2+2\left(\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}\right)^2\right]+\left(\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}\right)^2\left(\left(\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}\right)^2+\frac{d \Omega^2}{d \ln R}\right)=0 \]

と与えられる。

ただし、\(\kappa\)はエピサイクル振動数、\(\mathbf{u}_{\mathbf{A}}=B/\sqrt{4\pi \rho}\)はアルフベン速度を表す。

分散関係の可視化

ちょっと式をケプラー円盤の場合を想定して簡略化してみる。 ケプラー円盤なのでエピサイクル振動数は角振動数に等しいとして、また\(\frac{d\Omega}{dR}=-\frac{3\Omega}{2R}\)になりそうなので、

\[ \left(\frac{\omega}{\Omega}\right)^4-\left(\frac{\omega}{\Omega}\right)^2\left[\left(\frac{\kappa}{\Omega}\right)^2+2\left(\frac{\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}}{\Omega}\right)^2\right]+\left(\frac{\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}}{\Omega}\right)^2\left(\left(\frac{\mathbf{k} \cdot \mathbf{u}_{\mathbf{A}}}{\Omega}\right)^2-3\right)=0 \]

となる。これは、横軸\(k u_A/\Omega\)にとって、縦軸を\((\omega/\Omega)^2\)にとれば、このグラフが負になれば\(\omega^2<0\)となって不安定が起きるところが可視化されそうである。やってみると

(どうでもいいけどBingAIにコード書いてもらった。plt.contour(X, Y, Z, levels=[0]みたいな感じでおすすめされた)

これは、\(k u_A/\Omega<\sqrt{3}\)なら必ず\(\omega^2\)が負になる領域が存在することを意味している。 つまり、波数が\(k<\sqrt{3}\Omega/u_A\)を満たす摂動であれば、必ず不安定化する。 さらに最大成長率は、図の一番下に来てる点なので（導出は端折るが）、\(k u_A/\Omega=\sqrt{15}/4\)のとき、\(|\omega/\Omega|=3/4\)で達成される。

出典：Balbus and Hawley 1998, Balbus and Hawley 1991

驚くべき点？は以下のとおりである。

• 不安定になる必要条件は多くの天体で容易に達成される
  • criticalな波数、\(\omega^2\)が負から正に変わるときの波数は、ケプラー回転を仮定しない場合は\(ku_A=-\frac{d\Omega^2}{d \ln R}\)に対応している。小さい波数を取れれば不安定になる。じゃあどうやっても不安定が起きないためには、\(\frac{d\Omega^2}{d \ln R}>0\)なら良い。が、そんな天体は早々存在しなさそうである。
• 磁場がどんなに弱くても、小さい波数に対して不安定。
  • 磁場が弱いとき、すなわちアルフベン速度がめっちゃ小さくても、不安定が起こる。磁場が弱くても起こるので磁場が全然無視できないのである。
• タイムスケール
  • 最大成長率は、\(\omega\sim \Omega\)くらいで、すなわち回転のタイムスケールくらいで十分発展する。めっちゃはやい。

非理想MHD

ここまでは、理想MHDにおける磁気回転不安定性の発展の話をしてきた。 しかし、原始惑星系円盤では非理想MHDの効果が効いてくる。

Ohmic dissipation

誘導方程式を書くと、

\[ \frac{\partial \boldsymbol{B}}{\partial t}=\nabla \times(\boldsymbol{v} \times \boldsymbol{B})+\eta \nabla^2 \boldsymbol{B} \]

ここで、\(\eta\)はmagnetic diffusivityで、

\[ \eta= \frac{c^2}{4\pi \sigma} \]

である。

まず、誘導方程式を書き下すと

\[ \frac{\partial \boldsymbol{B}}{\partial t}=\nabla \times(\boldsymbol{v} \times \boldsymbol{B})-\nabla \times\left[\frac{4 \pi \eta}{c} \boldsymbol{J}+\frac{\boldsymbol{J} \times \boldsymbol{B}}{e n_e}-\frac{(\boldsymbol{J} \times \boldsymbol{B}) \times \boldsymbol{B}}{c \gamma \rho \rho_i}\right] \]

非理想MHDにおける分散関係

Sano and Miyama 1999