多次元輻射輸送・輻射流体力学シミュレーション
中本泰史（筑波大学・計算物理学研究センター）
概要：輻射輸送問題は，一般的には，配位空間３次元・運動量空間３次元の合わせて６次元の位相空間上の関数（輻射強度）を求める問題である．対象に高い空間対称性が特にない一般的な場合には取り扱う位相空間の次元が大きいため，大規模な計算はほとんど行われてこなかった．本稿では，空間次元が２次元・３次元といった多次元問題において輻射輸送を計算するための試みを簡単に紹介する．現時点では，多次元を扱えるようにするために何らかの近似を用いる必要があるが，そのやり方は，研究対象や求めたい計算の精度などに応じて，臨機応変に対処すべきである．
内容：
1. はじめに：輻射とは
2. 輻射輸送計算：定常問題
3. 輻射流体力学計算
4. 計算の実例
5. まとめと今後
1. はじめに：輻射とは「輻射」は宇宙物理学において，(1) 観測手段，(2) 「力」を及ぼす，(3) エネルギー輸送を担う，という意味において，重要な役割を持っている．大規模な観測装置により，天体の詳細な観測データが得られるようになればなるほど，その観測データからより良く真の天体の姿を導き出すためには，多次元空間内の輻射輸送のことを正確に考慮しなければならなくなる．
輻射は，輻射強度と呼ばれる量によって記述される（右図参照）．輻射強度は，光子の分布関数と，数係数のみの違いがある．すなわち輻
射強度は本質的には，光子の分布関数と理解しても良い．したがって，輻射強度は，配位空間３次元と光子の運動量空間３次元の合わせて６次元の位相空間上の関数となる．通常，運動量の３次元は，「方向２次元×振動数１次元」という呼び方で呼ばれることが多い．
輻射強度の変化を支配する方程式は，輻射輸送方程式と呼ばれる．輻射強度は光子の分布関数であると言うことから，Boltzmann 方程式との類推からも，輻射輸送方程式を導くことが出来る．

輻射輸送方程式で注意すべき点の一つは，散乱現象を考慮した場合，この方程式は積分・微分方程式になっているということである．このため，一般的に解を求めることは難しい．特に多次元の場合では，ほとんど不可能である．よって，数値計算を用いて解を求めることになる．

2. 輻射輸送計算：定常問題の解法
宇宙物理学において興味のある多くの問
題では，輻射が伝播する時間尺度よりも，流体が流れたり物質が輻射の影響で変化したりする時間尺度の方が短い．そうでない問題ももちろん数多く存在するが，以下では，上記の条件が満たされている場合を考える．その場合，輻射輸送方程式の時間微分項を無視して議論を進めることが出来る．
すなわち，問題は輻射場の定常解を求める問題に帰着される．
このときの一般的な解法は，次のようなものである．
(1)まず最初に，一本の光線上の輻射輸送を考える．定常輻射輸送方程式は，空間に関する微分と積分を含む方程式となっているので，これを形式的に積分する．すると，右のような形式解が得られる．これを数値的に評価することにより，光線上の輻射強度が求められる．
(2) ただし，散乱光を含む部分は積分項となっていて，このままでは調和的な解になっていない．そこで，反復法により，調和的な解を得る．反復法には，種々の方法が考案されている．代表的なものは，Feautrier法，Λ-Iteration 法， Accelerated Λ-Iteration法，Operator Perturbation 法，などである．

以上のような方法を用いることにより，一般的に，多次元空間内の定常輻射場を計算することが出来る．ただし，実際には計算量が膨大となるために，少しでも計算量を減らす工夫をしつつ，大型計算機を利用することが必要となる．計算量を減らす工夫には問題毎にいろいろあるが，代表的な一つは，重複している光線上の計算を何度もしないようにすることである．これを実現する方法として，
Short Characteristics 法 や

Accelerated Ray

Tracing 法などがある．大型計算機を効率よく利用するためには，並列化を考慮することが重要であるが，輻射輸送問題は，方向・振動数に対して輻射が独立なので，本質的に高い並列性を内在している．唯一，上流から下流に向かって情報が流れるという性質があるので，空間に対する計算順序に気をつける必要が出てくるが，この点も，克服することは可能である．実際，これらの工夫をこらし，大型計算機を効率よく利用すれば，大きな規模の３次元輻射輸送問題も数値計算することが可能になっている．

3. 輻射流体力学計算
流体の運動と輻射輸送を同時に考慮する「輻射流体力学」の計算法を簡単に紹介する．
これも一般的な方程式を書き下すことは出来るが，それをそのまま数値計算することは難しい．よって，問題毎に，よく見極めながら，近似を適用して解くべき式を導出すると言うことが行われる．
以下では，数値計算に用いられる基礎方程式の導出の一例を紹介する．ここでは，星形成問題など，非相対論的な問題を想定している．用いられる仮定は，右の通りである．

説明の簡単のため，１次元平行平板に対する方程式を右に書き下す．ここでは，一般的基礎方程式のうち，v/c の１次の次数の項までを残した．この式に対し，いろいろな状況を想定し，いずれの状況においても他の項よりも小さくなっている項を無視することにする．
そうすると，v/c の 0 次の方程式が得られ．
その方程式には，見かけ上は，v/c の 1 次の項も含まれているが，実際上は，v/c の 0 次の方程式と理解して良い．
得られた方程式系（右図参照）のうち，上の３つは，ふつうの流体力学の基礎方程式系に輻射の効果を追加した形になっている．輻
射の効果としては，輻射圧と輻射エネルギー輸送の効果が入ってきている．なお，輻射粘性の効果は，ここでは現れない．それは，v/cの1次の項まで正確に考慮して初めて現れる．
得られた方程式の４つ目と５つ目は，輻射のエネルギー保存・運動量保存に相当するものである．以上の５つの式は，３次元空間上の量を求める方程式系になっている．しかしこれだけでは方程式系は閉じていない（情報が足りない）ので，この段階では，６次元の輻射輸送方程式も解くことが要請される．
以上で数値計算で解くべき基礎方程式（の一例）が得られたが，これをどうやって得かという数値解法にも，いくつかのやり方がある．ここでは，３つほどのアルゴリズムを簡単に紹介する．
3.1 Flux Limited Diffusion 法
一般に，光学的に厚い状況に対して有効な計算手法として，拡散近似法がよく知られている．
これは，輻射エネルギーフラックスを物質温度の勾配と結びつけて計算するものである．この場合，輻射場は各点でほとんど等方的になることもあって輻射輸送方程式を解く必要がなくなり，３次元空間上の流体力学の方程式＋輻射フラックスの効果，のみを計算すればよいことになる．これは，計算量の観点から見ると非常に有利な点であり，このために，拡散近似が成り立つ状況に対しては，すでに多次元の数値計算もなされてきている．
しかし，拡散近似は光学的に薄い状況に対しては適用できない．そのような状況においては，輻射エネルギーが光速以上で伝播することになり，物理的に正しい解が得られなくなってしまう．この欠点を補うために，flux limitter という概念が導入された．これは，光学的に厚い状況に対しては，従来の拡散近似で得られるフラックスを与え，光学的に薄い状況においては，エネルギー輸送速度が光速になるように制限を与えるような関数を，輻射エネルギーフラックス

にかけておくことによって，物理的に不自然な解を除外するようにした方法である．この方法では実際，安定して数値計算を実行することが出来る．この計算法を，Flux LimitedDiffusion (FLD)法と呼ぶ．
右図には，flux limitter の一例を示した．FLDの長所は，拡散近似と同じく，輻射輸送方程式を直接解く必要がないので，計算量が３次元流体力学計算のものと同程度の（後の方法に比べて相対的に）小さいものとどまる，という点があげられる．短所としては，limitter の与え方が任意であるため，光学的厚さが 1 程度のあたりでの輻射エネルギーフラックスの評価が怪しいこと，物質の温度勾配と輻射エネルギーフラックスが常に平行と仮定しているが，一般的にはその仮定は成り立たないこと，などが挙げられる．しかし，一般的には短所もあるが，この方法の長所（計算が軽い）は，従来，あまり強調されていなかったように思われる．もっと見直しても良いのではないかと思う．
3.2 Simple Method
正式な名称は筆者は知らないので，便宜的に「Simple Method」と呼んでおく．後で述べる方法に比べて，より simple だからである．この方法では，拡散近似をやめ，きちんと輻射輸送方程式を解くことにより，そこから得られる情報を利用することを考える．このとき，輻射エネルギー密度や輻射エネルギーフラックスが輻射強度の方向積分で得られることを利用し，３次元方程式系の中に現れるこれらの量を，輻射輸送方程式から得られた輻射強度を用いて計算する．
この方法はアルゴリズムが（相対的に）簡単であるが，数値安定性に関しては，やや問題があると思われる．
3.3 Variable Eddington Factor Method
Simple Method の数値安定性を向上させるため，３次元方程式系と輻射輸送方程式の 間 の 繋が り を や や 弱 く した も の が ，Variable Eddington Factor (VEF) Method である．３次元方程式系の中に現れる輻射圧テンソルと輻射エネルギー密度の間の関係は，一般に，あるテンソル倍になっている．この係数を輻射輸送方程式から得られる輻射強度を用いて計算し，この係数が（微小時間の間は）一定であることを利用して，３次元方程式系を時間積分する．なお，３次元方程式系の中には，輻射エネルギー密度と輻射エネルギーフラックスに対する時間微分の式も含まれる．
この方法は一般に，Simple Method よりも数値安定性が高いと思われる．ただし，アルゴリズムが複雑になることが欠点である．
3.4 多次元輻射流体力学計算法の比較
ここで，いままで紹介した３つの方法の比較を簡単にしておく（右図）．参考のために，輻射輸送を考慮しない，単純な流体力学計算法の場合も示した．精度・計算量・複雑さ，などいくつかの観点があるが，それぞれの問題に応じて，適宜使い分ければ良いだろう．
ここでは，FLD 法のコストパフォーマンスの良さを指摘しておく．この方法は，計算量は単純流体力学計算の場合とそれほど変わらないにもかかわらず，輻射輸送を考慮することが出来るものである．確かに輻射輸送に対する精度に限界はあるが，より高精度の方法に比べてその不正確性はそれほど大きくなく，かつ，計算量は劇的に小さくてすむ．すなわち，そこそこの精度が軽い計算量で実現できるのである．したがってこの FLD 法は，もっとその利用価値が認識されても良いと思われる．

4. 計算の実例
ここでは，多次元輻射輸送計算・輻射流体力学計算を用いた研究の実例をいくつか簡単に紹介する．
4.1 １次元軸対称：フィラメン
ト状星間雲の重力収縮
無限に長い，１次元軸対称フィラメント状星間雲の重力収縮過程を，VEF 法を用いた輻射
流体力学計算により調べた．注
目点は，星間雲が等温から非等温に変わるところである．輻射
輸送過程を詳細に考慮した結果，この遷移の様子を正確に知ることが出来た．
この計算の結果明らかになったことは，「等温から非等温への遷移」と「星間雲の光学的厚さが１を超えること」は，独立の現象であると言うことである．これはよく勘違いされているが，一般的には別の現象であることを認識しておく必要がある．
この計算の応用として，形成される「星」の最小質量を見積もることが出来る．それは現在の星形成領域においては，約 0.007 太陽質量である．これは，褐色矮星の最小質量がこの程度であることを示唆している．実際にこれまでに観測されている褐色矮星の最小質量(約 0.01 太陽質量)は，この計算結果と調和的である．
4.2 大 質量形 成過 程：２ 次元
FLD による輻射流体力学計算

4.3 T Tauri 型星の構造の推定：２次元輻射輸送計算
特 異 な ス ペ ク ト ル (SpectralEnergy Distributioin; SED)を持つ
T Tauri 型星の構造を推定するため，モデル計算を行った．この天体群は，星周円盤の外側が，従来のモデルでは理解できないくらいに高温になっている．こ
こではその理由として，薄い密度で存在するハローによって中心星からの輻射が反射され，それが円盤外側を加熱しているという「ディスク・ハローモデル」を考えた．
このモデルに基づき，実際に観測されるような SED が得られるのかどうかを，２次元輻射
輸送計算を行うことによって確かめた．

計算の結果得られた物質の密度分布と温度分布が右図に示されている．色で示されているのが密度分布，線で示されているのが温度分布である．これで見られるように，温度分布は決して球対称的なものではない．また，直接は示されていないが，温度分布の勾配と輻射エネルギーフラックスが平行でない場所もあることも見て取れるだろう．

このような温度・密度分布を持つ天体からの SED を計算したものが，右図である．線で示してあるのは，数値計算で求められたもの．点で示してあるのは，観測結果である．このモデ
ルでは，観測から得られる SEDを良く再現していることがわかる．

4.4 原始星の構造の推定
4.3 と同様にして，T Tauri 型星よりも，より若い段階の天体である「原始星」の構造を推定することも出来る．
「原始星」を記述するために多くのパラメータが必要になるが，適切な手順を踏めば，これらを観測データからひとつひとつ求めて決定していくことが出来る．

実際，そのような手順を踏んで「原始星」の構造を推定した結果得られた SED と観測結果を比較したものが，右図である．
計算結果（実線）と観測結果（点）が良く一致していること，すなわち，構造の推定がうまく出来たことがわかる．

4.5 宇宙の再電離過程：３次元輻射輸送計算
非一様な密度場を持つ宇宙が紫外線によって再電離される状
況を，３次元輻射輸送計算により求めた．
この計算では，密度場はZeldovich 近似により与え，周囲から一定の強度の等方的な紫外線を入射させ，その結果実現される電離構造を計算した．
右図は，計算の結果得られた電離度分布の一例である．これは，Z=15 の宇宙におけるものである．
赤方偏移が変わると密度が変わるので，同じ強度の紫外線に対しても，電離の仕方が変わってくる．それを示したのが，右図である．Z が小さいと平均密度は小さくなるので，ガスはより電離しやすくなることがわかる．

なお，３次元空間内の輻射輸送を正確に計算した結果として，３次元輻射輸送に特有な効果である，遮蔽効果が実際に計算の中でも観測された．

4.6 宇宙の再電離過程：３次元輻射流体力学計算
宇宙の再電離過程を，輻射流体力学的に解いて調べようとす
る Gnedin らの試みを紹介する．
Gnedin らは流体力学計算コードをベースにして，それに輻射輸送の効果を取り込んで計算しようとしている．ただし，正確に輻射輸送を考慮したのでは計算量は莫大になって実行が不可能なので，輻射輸送の効果は近似的にしか扱えない．
Gnedin が試みた方法の一つ
は Local Optical Depth 近似と呼ばれるもので，注目する点と光源の点の２点における光学的厚さのみから注目する点に到達する輻射量を評価するものである．これは，輻射が伝播する途中の媒質の情報を一切無視しており，明らかにおかしい近似であるが，計算を実行可能にするために，敢えて彼らはこのような近似法をとった．
もう一つの方法は，Optically Thin Variable Eddington Tensor 近似法と呼ばれるもので，実際の光学的厚さに関係なく，光学的に薄い場合の Variable Eddington Tensor (Factor) の値を用いて計算を進めるものである．これも先と同様，途中の媒質の情報を一切無視するもので，物理的には明らかに正しくない近似法である．
上の図は，この２つの方法を用いて同じ問題を計算した結果を比較したものである．このような簡単な近似法の間でも，結果は大きく異なっている．もちろん，物理的には両者とも正しい近似とは思えない．したがって，現在発表されているような３次元輻射流体力学計算による結果は，十分注意して見る必要があると思われる．

4.7 銀河形成：３次元 SPH による輻射流体力学計算
QSO の近傍におけるガス雲から銀河が形成される際には，そのガス雲のガス密度と，QSOからの紫外線の強度によって，形成される銀河の姿が大きくことなくと考えられる．このことが，SPH と輻射輸送計算を結合させた輻射流体力学計算により，確かめられている．

右図は，紫外線が（相対的に）弱く，ガス密度も高くて遮蔽効果が有効に効いている場合のガス雲の進化を示している．この場合には，ガスは十分冷えることが出来て，大量の星（青色で表示された中心部の密集物）が形成される．銀河の形状としては，楕円銀河になっている．

一方，紫外線が（相対的に）強い場合には，ガスはなかなか冷却できず，星はあまり作られない（右図）．この場合には，星が少なく，ガスの残りが多い円盤銀河に近い銀河が形成される．

5. まとめと今後

参考文献：「Foundations of Radiation Hydrodynamics」Mihalas & Mihalas, 1999, Dover
輻射流体力学の基礎が非常に詳しい．主に連続光を扱うことを想定して書かれている．
「Radiative Transfer in Moving Media」Sen & Wilson, 1998, Springer
１次元球対称問題に特化しているが，ラインの輻射輸送などについても一通り書かれている．