応用課題¶

著者:	高橋博之 (国立天文台)

モーメント法を用いた相対論的輻射流体の数値解法¶

近年のスーパーコンピュータの発達と数値計算コードの高精度化に伴い、より高解像度・長時間の磁気流体計算を行う事が可能になってきた。その一方で、磁気流体に自己重力や、特殊・一般相対性理論、電気抵抗、熱伝導、粘性、輻射輸送など、より多くの物理を含めた現実的な計算も行われている。この章では相対論的流体に輻射の効果を無矛盾に取り入れた相対論的輻射流体数値計算の紹介をする。輻射場は簡単のため、光子分布関数を運動量空間で積分したモーメント量(輻射エネルギー密度、フラックス)で記述する。モーメント方程式の導出は省略するが、詳しくは数値天文学マニュアル(http://www.astro.phys.s.chiba-u.ac.jp/hpci/ws2012/TechMan2012/)や Mihalas & Mihalas (1984), Kato et al. (2008) , Takahashi et al. (2013) を参照してほしい。

以下、メトリックは平坦とし( $\mathrm{diag}\ g_{\mu \nu}=(-1,1,1,1)$ )、添字のギリシャ文字とアルファベットはそれぞれ時空4次元( $0,1,2,3$ )と空間3次元 $(1,2,3)$ を表す。

相対論的輻射流体方程式¶

相対論的輻射流体方程式は以下で記述される:

(1) $\frac{\partial (\rho \gamma)}{\partial t} + \frac{\partial (\rho u^j)}{\partial x^j}= 0,$

(2) $\frac{1}{c}\frac{\partial}{\partial t} \left[(\epsilon + p_g)\gamma^2 - p_g\right] + \frac{\partial}{\partial x^j} \left[(\epsilon + p_g)\gamma \frac{u^j}{c}\right]= G^0,$

(3) $\frac{1}{c^2}\frac{\partial}{\partial t} \left[(\epsilon + p_g)\gamma u^i \right] + \frac{\partial}{\partial x^j} \left[(\epsilon + p_g)\frac{u^i u^j}{c^2} + \delta^{ij} p_g\right]= G^i,$

(4) $\frac{\partial E_r}{\partial t} + \frac{\partial F^j_r}{\partial x^j} =-cG^0,$

(5) $\frac{1}{c^2} \frac{\partial F^i_r}{\partial t} + \frac{\partial P^{ij}_r}{\partial x^j} =-G^i.$

上から順に質量保存、ガスのエネルギー保存、ガスの運動量保存、輻射の0次モーメント式、輻射の1次モーメントの式を表す。 $\rho, \ p_g$ はガス密度、ガス圧で $\delta^{ij}$ はクロネッカーデルタである。 $\epsilon$ は静止質量エネルギー密度を含めた流体のエネルギー密度を表しており、polytropic関係を仮定すると内部エネルギー密度 $e$ 、比熱比 $\Gamma$ を用いて

(6) $\epsilon = \rho c^2 + e = \rho c^2 + \frac{p_g}{\Gamma-1},$

とかける。 $u^\nu=(c \gamma, \gamma v^j)$ は流体４元速度で、3元速度 $v^i$ 、ローレンツ因子 $\gamma$ との間に $u^i = \gamma v^i = v^i/\sqrt{1 - \beta^2}$ の関係がある。ただし $\beta^i=v^i/c, \ \beta = \sqrt{\beta_k \beta^k}$ である。

$E_r [\mathrm{erg\ cm^{-3}}],\ F_r^i [\mathrm{erg\ cm^{-2}\ s^{-1}}],\ P_r^{ij} [\mathrm{erg\ cm^{-3}}]$ はそれぞれ輻射場のエネルギー密度、フラックス、ストレスで、輻射とガスは輻射4元力密度 $G^\mu$ を通してエネルギー運動量のやり取りを行う。

(7) $G^0 =-\frac{\rho\gamma \kappa_0}{c} \left(4\pi \mathrm{B} - c E_r +\beta_j F_r^j\right) -\frac{\rho \gamma \sigma_0}{c} \left[\frac{u^2 E_r}{c} + \frac{u_j u_k P_r^{jk}}{c} - \left(2\gamma^2 - 1\right)\beta_jF_r^j \right],$

(8) $G^i = &-\frac{\rho \gamma \kappa_0}{c} \left(4 \pi \mathrm{B} \beta^i - F^i_r +v_k P_r^{ik} \right) \\ &+\frac{\rho \gamma \sigma_0}{c} \left[F_r^i - \gamma E_r u^i - v_k P_r^{ik} +v^i\left(\frac{2 \gamma u_k F_r^k}{c^2}-\frac{u_j u_k P_r^{jk}}{c^2}\right) \right].$

ただし $u = \sqrt{u_k u^k}$ である。 $\kappa_0$ , $\sigma_0$ はそれぞれ共動座標系で評価した吸収係数、及び散乱係数( $\mathrm{cm^2 \ g^{-1}}$ )を表す。 $\mathrm{B} [\mathrm{erg\ cm^{-2}\ s^{-1}}]$ は黒体輻射強度でガス温度を用いて

(9) $\mathrm{B} = \frac{\sigma_\mathrm{SB}}{\pi}T_g^4$

となる。ここで $\sigma_\mathrm{SB}=5.67\times 10^{-5}~\mathrm{erg \ cm^{-2}\ s^{-1}\ deg^{-4}}$ はStefan-Boltzmann定数である。このガス温度は状態方程式によって決まる:

(10) $p_g = \frac{\rho k_B T_g}{\mu m_p}.$

$k_B=1.38\times 10^{-5}\ \mathrm{erg \ deg^{-1}}$ はボルツマン定数、 $m_p=1.672\times 10^{-24}\ \mathrm{g}$ はイオン質量、 $\mu$ は平均分子量である。

これらの式に加えて、クロージャー関係と呼ばれる $P_r^{ij}$ を与える式(輻射の状態方程式)が必要となる。ここでは2つのクロージャー関係式を紹介する。

共動座標系で等方的な放射を仮定した場合、 $P_r^{ij}$ は

(11) $P_r^{'ij} = \frac{\delta^{ij}}{3}E_r',$

と表される。このように等方的な輻射圧を与える近似を Eddington近似 と呼ぶ。ダッシュ $'$ は共動座標系の物理量を表しているが、式 (4) - (8) で必要なのは観測者系のストレス $P_r^{ij}$ である。そこで式 (11) にローレンツ変換を用いることで $P_r^{ij}$ に対する関係式が得られる:

(12) $P^{ij}_r &+ \left[ -\frac{\delta^{ij}}{3}\frac{u_k u_l }{c^2} +\frac{u^i u^j u_k u_l}{(1+\gamma)^2 c^4} \right] P_r^{kl} +\frac{u^i u_k P_r^{jk} + u^j u_k P_r^{ik}}{(\gamma+1)c^2} \\ &= \frac{\delta^{ij}}{3}\gamma^2 \left(E_r - 2 \frac{\beta_j F_r^j}{c}\right) -\frac{u^i u^j E_r}{c^2} +\frac{u^i F_r^j + u^j F^i_r}{c^2} +\frac{2 u^i u^j u_k F_r^k}{(\gamma+1)c^4}.$

$P_r^{ij}$ は $3\times 3$ の対称行列のため、 $E_r, F_r^i, u^i$ を与えて $6\times 6$ の線形方程式を解くことで $P_r^{ij}$ が得られる。

Eddiington近似は光学的に厚い状況では正しいが、光学的に薄い場合放射は等方的にならないため正しくない。非等方な放射を考慮したクロージャーとして M-1クロージャー が提案されている:

(13) $P_r^{ij} = \left[\frac{3(1 - \chi)}{2}\frac{\delta^{ij}}{3} + \frac{3\chi - 1}{2}n^i n^j\right]E_r$

(14) $\chi = \frac{3 + 4 f^2}{5 + 2\sqrt{4 - 3 f^2}}, \ n^i = \frac{F^i_r}{F_r}, \ f^i = \frac{F^i_r}{c E_r},\ f = \sqrt{f_k f^k}, \ F_r = \sqrt{F_{r,k} F_r^k}.$

課題¶

熱平衡状態への遷移 “1d_lte/”¶

ガス温度と輻射温度( $T_r=(E'_r/a_R)^{1/4}$ 、ただし $a_R=4\sigma_\mathrm{SB}/c$ )が等しくない場合、ガスと輻射は吸収・再放射を通してエネルギーのやり取りを行い、熱平衡状態( $T_g = T_r$ )へと近づいていく。この節では非熱平衡状態から熱平衡状態へと遷移していく過程を考える。

仮定としてガスは定常・静的とし、ガス密度( $\rho = 0.025\ \mathrm{gcm^{-3}}$ )、ガス温度( $T_g = 10^6\ \mathrm{K}$ )は一定とする。吸収係数は $\kappa_0=0.4\ \mathrm{cm^2\ g^{-1}}=\mathrm{const}$ とし、散乱係数は $\sigma_0=0$ とする。輻射温度は初期に $T_r=10^8\ \mathrm{K}$ とする。

吸収を通して輻射エネルギーが減少し、熱平衡状態に近づいて行くことを確認してください。特に$t-log E_r$でプロットした場合、どのように輻射エネルギーが減衰するかを確認してください。
輻射エネルギーが指数関数的に減衰する理由を考えてみてください。(ヒント: 式 (4) で $F_r = 0, \sigma_0 = 0, v^i=0, \mathrm{B}=\mathrm{const}, \frac{\partial}{\partial x^j}=0$ とする。)

1次元輻射の伝搬 “1d_transport/”¶

輻射場がガス中を伝搬すると、輻射場は吸収過程により減衰しながら伝搬していく。この節では1次元一様静的なガス中を伝搬する輻射場を考える。

計算領域のサイズは $L=100\ \mathrm{cm}$ で、ガス密度、ガス温度はそれぞれ $0.01\ \mathrm{g \ cm^{-3}}, \ 10^6\ \mathrm{K}$ とし、輻射場はガスと熱平衡にあるとする。吸収係数は $\kappa_0 = 10\ \mathrm{cm^2 \ g^{-1}}=\mathrm{const}$ とし、散乱係数は $\sigma_0=0$ とする。この一様媒質中に $x=0$ から輻射を $T_r =10^{8}\ \mathrm{K}$ で照射する。

$x-\log E_r$ 上で輻射が時間とともに伝播していく様子を確認してください。
吸収係数を変更した場合( initial.f90 内の変数 absorption )、どのように変化するか確認してください(ただし大きすぎると数値的に解けなくなります)。
クロージャー関係を変えた場合( initial.f90 内の変数 radiation を1から2に変更)、波面の伝播する速度が変わることを確かめ、なぜそのようになるか考えてください(ヒント: 光学的に薄い極限を考え、式 (4) - (5) で $G^\mu=0$ とし、Eddington近似の場合は $P_r^{ij} =\delta^{ij}/3$ 、M-1クロージャーの場合は $P_r^{xx}=E_r,$ その他 $\ =0$ として特性速度を求めてください)。
$x=0$ から注入された光の波面が $x=100\ \mathrm{cm}$ の境界を通過すると定常状態になるの定常状態になったときの、 $x-\log E_r$ の傾きをEddington近似を用いた場合とM-1クロージャーを用いた場合でそれぞれ求めてください。また、その傾きが何によって決まるかを考えてみてください。(ヒント：式 (4) (5) で定常・静的なガスを仮定し、 $E_r \gg 4\pi \mathrm{B}/c$ とする。)

2次元輻射場の伝播〜日陰問題〜 “2d_shadow/”¶

2次元輻射場伝播問題の初期条件。

輻射衝突の計算。

この節では光学的に厚い物体に光を照らした場合の計算を行う( 2次元輻射場伝播問題の初期条件図 )。密度 $\rho=10^{-5}\ \mathrm{g\ cm^{-3}}$ 、温度 $T_g=10^{6}\ \mathrm{K}$ の熱平衡にある一様媒質中に半径 $20\mathrm{cm}$ の重たいクランプ $\rho = 1\ \mathrm{g \ cm^{-3}}$ を置く。このクランプに対して $x=-100\ \mathrm{cm}$ の境界から温度 :math`T_r = 10^8mathrm{K}` の輻射を照射する。吸収係数は $\kappa_0=1\ \mathrm{cm^2 \ g^{-1}}$ 、散乱係数は $\sigma_0 =0$ とする。

Eddington近似を用いた場合とM-1クロージャーを用いた場合でどのように結果が変わるか確認してください。
$\kappa_0=0.0, \ \sigma_0=1\ \mathrm{cm^2 \ g^{-1}}$ とした場合( initial.f90 内 absorption = 0, scattering=1 )、どのように結果が変わるか確認してください。

2次元輻射場の伝播〜輻射の衝突〜 “2d_collision/”¶

光学的に薄いガス中を伝播する２本の光柱を考える( 輻射衝突の計算の設定図 )。ガスは一様静的とし、密度 $\rho=10^{-5}\ \mathrm{g\ cm^{-3}}$ 、温度 $T_g=10^{6}\ \mathrm{K}$ とする。以下の様に境界から光を注入し、二つの光が交差する様子を確かめる。

(15) $E_r = a_R T_r^4,\ F_r^x = F_r^y = \frac{c E_r}{\sqrt{2}}, & \ if\ x \le -100\ \mathrm{cm}\ and \ y = [-80\ \mathrm{cm}, -60 \ \mathrm{cm}],\\ E_r = a_R T_r^4,\ -F_r^x = F_r^y = \frac{c E_r}{\sqrt{2}}, & \ if\ x \ge 100\ \mathrm{cm}\ and \ y = [-80\ \mathrm{cm}, -60 \ \mathrm{cm}],\\$

ただし $T_r=10^8\ \mathrm{K}$ とする。吸収係数、散乱係数はともに $0$ とし、M-1クロージャーを用いる。

2本の光柱が真っ直ぐに進み、衝突することを確認してください。
光同士は相互作用をしないため、光が衝突することは非物理的である。なぜこのような事が起きるのか考えてみてください。
2本の光の強度が異なる場合にどうなるか調べてみてください。

(発展編)2次元輻射場の伝播〜ガスの加速〜 2d_acc/¶

ここでは輻射によってガスが加速されて行く様子を計算する。計算ボックスは2次元 $x=[-10^{4}, 3\times 10^4]\ \mathrm{cm}$ 、 $x=[-2\times 10^{4},2\times 10^4]\ \mathrm{cm}$ とし、一様なガス( $\rho = 10^{-6} \ \mathrm{g\ cm^{-3}}, T_g = 10^9 \ \mathrm{K}$ )で満たす。ここに円筒形のガスクランプ(半径 $3\times 10^{3}\ \mathrm{cm}$ 、密度 $10^{-2}\ \mathrm{g\ cm^{-3}}$ )をおき、 $x=-10^{4}\ \mathrm{cm}$ から輻射を照射する( $T_r = 3\times 10^{7}\ \mathrm{K}$ )。この問題ではガスの運動を追うため、静的なガスを仮定しない ( initial.f90 内、変数 hydro = .true. )。ガスクランプが加速されていく様子を確かめてください。