J. Plasma Fusion Res. Vol.88, No.9 (2012) 467‐474

講座 輻射流体シミュレーション

１．輻射流体シミュレーションの基礎
大西直文
東北大学大学院工学研究科航空宇宙工学専攻（原稿受付：２
０
１
２年５月２
５日）

輻射流体シミュレーションの基本的な構造と，解くべき方程式系の性質について概説する．また，現在利用されている代表的な輻射輸送計算手法について紹介し，それらの特徴と利点，欠点についていくつかの観点から比較を行う．現在の輻射流体シミュレーションにおいては，輻射モーメント方程式に基づく数値解法が最もよく輸送計算に用いられているため，その背景を説明した後に，古典的な流束制限拡散法に始まり非等方性を考慮できる手法までエディントン因子の評価法という側面からまとめて紹介する．
Keywords:
radiation transfer, laser-produced plasma, astrophysical flow, flux-limited diffusion, Eddington factor

１．
１ はじめに体中の原子過程という３つの要素で構成されるが，本講座

高温プラズマでは，電磁波や中性子，ニュートリノなど，

ではそのうち後者２つに焦点を当て，それらの計算手法の輻射と呼ばれる比較的平均自由行程の長いエネルギーキャ

概要を紹介することで，これから輻射流体コード開発を試リアが，しばしばプラズマの流体力学的挙動を大きく左右

みる読者の一助としたい．

する．特に，エネルギー密度の高い爆発的現象において輻

輻射輸送方程式の解法や流体計算との結合方法は，プラ射の寄与が顕著となることが多く，それゆえ衝撃波の物理

ズマの光学的厚さや輻射強度に応じて適切に選択する必要とともに研究されてきた［１］．そのような輻射流体力学的がある．次節に続く本章では，代表的な輻射輸送計算手法

現象は，宇宙物理においては古くから研究されており，実を挙げながらどのような手法を選択すべきかを述べる．典験室においてもレーザー生成プラズマの基本的挙動を決定

型的な輻射流体シミュレーションでは，輻射輸送計算に全する．例えば重力崩壊型の超新星爆発では，原始中性子星

計算時間の９割以上を費やすため，なるべく計算負荷の低から放射されるニュートリノがバウンス時に生成される衝い手法を選択することが求められる．また，レーザープラ

撃波を加熱し，爆発現象に大きく寄与すると考えられていズマのように輻射場が等方に近い場合に用いられる輻射

る［２］．また，慣性核融合ではそもそも輻射によって爆縮モーメント方程式の解法について詳説し，そこに非等方効を駆動する方式もあるため，現象の予測に輻射流体力学の果を取り入れることのできる可変エディントン因子法を紹

理解が不可欠である［３］．

介する．第２章では，超新星爆発におけるニュートリノ輸

ところが，そこでは連続体と飛程の異なる粒子（光子を送を例に，非等方性の強い輻射場における輸送計算手法に含む）が混在した輸送現象や，さらにその粒子の生成，消ついて紹介する．角度方向の情報を保持する輸送計算手法

失過程が複雑に絡んでいる．光と流体という異なる時間・は一般に計算負荷が高くなるが，大型計算機の発展に伴

空間スケールの現象を含む物理の解析的な取り扱いは難し

い，そのような手法を用いた輻射流体計算が実現可能にく，数値シミュレーションの担う役割が極めて重要になっなってきている．代表的手法である !! 法について，ニューてくる．現在，世界的に見ると，輻射流体コードと呼ばれトリノ輸送計算への応用を紹介する．最後に第３章では，るものは宇宙物理と慣性核融合のコミュニティを中心とし輸送計算に必要となる放射係数，吸収係数がどのように算て多数存在する．しかし，輻射流体力学は，流体力学，熱出されるかについて述べる．流体の輻射特性は，それを構

力学，統計力学，運動論，原子核物理など，専門性の異な成する原子が輻射の放出・吸収によって状態遷移すること

るいくつかの要素から構成されるため，コード開発には分で決まる．この過程を原子過程と呼び，適切なモデル化が野を横断する複数の研究者によるプロジェクト的な枠組み必要となる．代表的なモデル化とその計算手法，および輻

が必要である．残念ながらそのような組織的研究が不得意

射スペクトルの表現方法について解説する．

な我が国では，成熟した輻射流体コードが存在していると

１．
２ 輻射流体コードの構造は言い難く，したがって輻射流体コードに関する指南書も十分でないのが現状である．輻射流体コードは大きく分け

輻射流体力学の方程式は，輻射と相互作用する物質て流体，輻射輸送，およびそれらの相互作用を決定する流（以降，流体と表現する）を支配する流体方程式をベースに

1. Fundamentals of Radiation Hydrodynamics Simulation
OHNISHI Naofumi

author’s e-mail: ohnishi@rhd.mech.tohoku.ac.jp

４６７!2012 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

し，輻射との相互作用項を付加するという形で表すことがまま輻射流体コードの構造を表す．輻射流体コードを構成できる．例えば流体が非粘性圧縮性である場合，輻射流体

する３つの要素に対しては，それぞれ独立に数値計算手法

方程式は以下のようになる［４］．

が発展している．多くの輻射流体コードでは，各要素に対応した支配方程式を別々の手法で解き，要素間で物理量の

(
'
'#(
!&"
'
#!
(
'

（１）

('
'#!!"
$#(
#!
'##!&!P(
!&"
'
(
'

（２）

している．

('
$!"(
P(
#!!*
$!&(
!&"
)
'
##"
'
&"
(
'

（３）

性質，特に解を得るにあたっての困難さの要因となるもの

やりとりをすることで結合計算を行う方式を採用している．図にはそのやりとりに用いられる代表的な物理量を記数値計算法の紹介に入る前に，もう少し支配方程式系のについて述べておく必要がある．まずはこの方程式系が極

ここで '，#，&，$はそれぞれ密度，流速，圧力，流体のめて多次元の問題であることに注意されたい．それは，後全エネルギーを表す．流体に対する付加項は，運動方程式

述するように輻射輸送方程式が輻射を担う粒子に対するボに輻射エネルギー流束 !の時間微分項と輻射圧力テンソル

ルツマン方程式であること，さらには，原子過程が多電子

P の発散項，エネルギー保存式には輻射エネルギー密度 "系の量子力学に支配されることに由来する．したがって，の時間微分項と輻射エネルギー流束の発散項，および輻射

輻射流体力学の方程式系を「真面目に」解くことは断念せ圧による仕事の項 #"P(が入る．ただし，流れ場が相対'&"ざるを得ない．さらに方程式系には硬直性がある．例えば論的である場合，相互作用を流体の共動系で評価する必要

多くの場合，流体の伝搬速度は音速であり，輻射の伝搬速があるため，方程式の形が若干異なってくる．また，宇宙度は光速であるため，それらの違いが硬直性を招くことは流体現象ではしばしば流体の圧力に比べて輻射圧が卓越す

容易に想像できる．このため，収束解を得るには多大な演るが，レーザー生成プラズマでは無視できる場合が多いた

算量を費やすことになる．

め，式（２）の輻射に関連する項，および式（３）の右辺は省以上のことから，輻射流体シミュレーションには如何に略されることが多い．輻射流体力学方程式の詳しい導出に「サボる」かという観点が必要となる．実のところ，輻射流

ついては，参考文献を参照されたい［４‐６］．

体力学の解を正確に知ることなど，ほとんどの場合不可能

輻射に関連する量 "，!，P は，周波数 %の輻射強度なのである．つまり，輻射流体力学の妥当な解を得ること（specific intensity）#%に対する伝播方向 Ω のモーメントとは，「適当な」近似をして演算量を減らし，必要な計算精して定義される．

度と使用できる計算資源のトレードオフの中でベストな解%%
Ω%%#! "%%%"$! ! "#
%%
! $&$
!
%

%

Ω%%#! !%%%!$! ! #
%Ω%
!

$&

%

!

%

P $! ! "#Ω%%#! P %%%
%ΩΩ%
! $&$
!

を探す作業になる．したがって輻射流体シミュレーション（４）の研究者には，様々な近似解法とその適用範囲を知り，場合によって解法を使い分けることが求められる．次節では（５）

輻射輸送方程式について概要を説明した後，代表的な近似解法をいくつか紹介する．

（６）

１．
３ 輻射輸送方程式

ここで，$，"%，!%，P %はそれぞれ光速，周波数 %におけ以降本章では輻射は光子（電磁波）であるとし，また，相る輻射エネルギー密度，輻射エネルギー流束，輻射圧力テ対論的な効果は無視する．輻射輸送方程式は，光子の分布

ンソルである．そして位置 "と方向ベクトル Ω の関数であ関数 # についてのボルツマン方程式から導出でき，輻射強

Ω!
る#
"
!
'
(は，後述する輻射輸送方程式によって支
%##
%'

度#
%%
# について以下のように記述される．
%#$

配される．ただし，輻射流体力学と呼ばれる現象においては，流体と輻射の方程式が独立に成立するのではなく，流体中の原子過程に従って輻射が流体から発せられ，また同時に流体に吸収，散乱されるため，それらは相互に依存している．さらに，流体の内部エネルギーは電子とイオンの熱運動だけでは決まらず，特に部分電離プラズマではイオ
ン化ポテンシャルの部分が比熱 !( の大部分を占めることになるため，原子過程を適切に考慮することで状態方程式という形でこの効果を流体へ取り入れる必要がある．すな
わち輻射流体力学は複合系の力学であり，方程式系は高度に非線形となる．
図１

長い距離を伝搬する輻射，中距離の輸送を担う流体，そして微視的スケールの原子過程と，輻射流体力学は多階層構造の問題とみることもできる．このような輻射流体力学の現象論的構造を示したのが，図１である．この図はその４６８

輻射流体コードの構造．輻射流体コードは大きく３つの要
素に分けられ，それぞれの要素間で物理量のやり取りが行われる．図中の括弧で囲った物理量は，レーザープラズマ
では通常考慮されないが，状況によっては（慣性系の変換や輻射による原子の励起および電離）必要になる情報を表す．

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishiが必要ないことがわかる．この事実が，後に述べるモーメ

#
!,
'
"Ω#
'#
'$&
'!#
'#
'
%,
*

ント方程式を用いた輻射流体方程式の解法に繋がる．
輻射輸送方程式の解は，方程式中の相互作用項がどの程

Ω!
Ω!
Ω ,)
Ω ,!
Ω)
Ω ,!
"
!
*
%(
*
%(
#Ω ,
)
+
)
!#
"
#
!
"*'! *
'(
'(
"(
（７）

度，あるいはどのように寄与するかに大きく依存する．本

ここで，&
'，#
'はそれぞれ放射係数，吸収係数を表す．まを基にした解の性質の分類について述べる．また，最後にΩ)は方向 Ω ,から方向 Ω への
Ω ,!
た，*'は散乱係数で，%(輻射輸送方程式の近似解法を数例紹介し，その特徴につい節の残りでは，いくつかの重要となる概念を挙げ，それら

散乱カーネルを表す．#'は定義より明らかなようにエネルてまとめる．

ギー流束の単位を持つが，独立変数として位置３次元と方１．
３．
１ 光学厚さ向２次元を持つため，ある点におけるある方向のみの，つ輻射輸送方程式の解の性質を分類するために最もよく使まり単位立体角中のエネルギー流束を表しており，ベクト

われる指標が光学厚さである．光学厚さ +'は以下のようにル量としてのエネルギー流束 !'とは異なる．

定義される．

光子はその性質から伝搬速度が変化しない．したがっ#'"*')+#)('$!

て，式（７）の左辺には通常のボルツマン方程式にみられる（１４）

ような加速度による速度空間の移流項がない．右辺は流体との相互作用項であり，第１項と第２項はそれぞれ輻射のここで，)は輻射の伝搬距離を表す．ほぼ同じ概念に光学

放出，吸収を表す．これらの項によって，流体から放出さ深さがあるが，光学深さは特定の区間に対して使われるのれたエネルギーは輻射に，輻射が吸収された場合は流体へ

に対し，光学厚さは系全体の値として語られる場合が多とエネルギーが移動する．式（７）の右辺第３項は散乱項でい．ここでは，特定の区間について議論するわけではないある．ここでは，コンプトン散乱等のエネルギー変化を生ので，光学厚さという言葉を使う．

じる散乱は考えていない．そのため，輻射は流体によってさて，この光学厚さであるが，吸収係数や散乱係数がそ散乱するが，流体へのエネルギーフィードバックはない．

れぞれに対応する平均自由行程 &$!
"'$の逆数で与えられ散乱によって光子の伝搬方向は変化するため，積分中の第ることから，積分区間を通して光子がどれだけ流体と衝突

１項 は 別 の 方 向 か ら 考 え て い る 方 向 に 入 っ て く る

するかを表す．ここで，'は衝突相手の数密度，$は衝突断（scattered-in）光子，第２項は今の方向から別の方向に移面積である．+'&!であれば光子にとって流体場は不透明

動する（scattered-out）光子の寄与を表している．

（opaque）であり光 学 的 に 厚 い（optically thick），+'%!

ここで後の議論のために，輻射輸送方程式に対し方向ベであれば流体場は透明（transparent）で光学的に薄い（op-クトル Ω のモーメントを取ってみる．簡単のため散乱の寄tically thin）という．

与を無視すると，０次および１次のモーメント方程式はそ

まず，光学的に厚い場合を考えよう．光学的に厚い場合，れぞれ以下のように表される．

光子は流体中ですぐに吸収あるいは散乱されてしまう．とこ ろ で，局 所 熱 平 衡（local thermodynamic equilibrium:

,
"'
!'$"
(&
#'"'
"'#
'!%
,
*

（８）

!,!'
P'$!#'!'
"%'#
%,
*

（９）

LTE）状態にある系ではキルヒホッフの法則から放射係数が吸収係数に比例する．
&
'$#
'!'

これらの両辺をさらに周波数 'で積分すると，（１５）

ここで，!'はプランク関数で温度の関数である．そのため，吸収が顕著な場合は放射係数も高くなっており，多数

,
"
!$!
(&
#'"')
"
#'
"'#
(
'!%
,
*

（１０）

!,!
P $!!#'!'#'
"%'#
%,
*

（１１）

の光子が放出されている．放出された多数の光子がその周辺で吸収されると，そこで無視できない加熱が起こり再び放出が起こる．結果として光子は吸収，放出を繰り返し，放出は等方的に起こるため，光子は実効的にランダムとなり，式（２），（３）に現れる輻射に関連した項は，これウォークする状況になる．散乱が卓越する場合も多数の散らの方程式を用いて書き換えることができる．

乱を繰り返すため，やはりランダムウォークする．した,()#))##"()"'#($!!#'!'#',*

がって，光学的に厚い場合，輻射輸送は拡散的な現象にな（１２）り，後述するように輻射輸送方程式は方向２次元を落として輻射エネルギー密度の拡散方程式へと帰着することがで,$P （１３）
#+
$"()
(&
#'"')
"
#'"##
"'#
*
(
$!
(
(
'#)
'!%
,
*

きる．
では，光学的に薄い場合はどうであろうか．光学的に薄例えば，式（１３）の右辺第２項を無視すると，流体のエネル

い場合，光子はほとんど流体と衝突することなく系をすりギー方程式のソース項がその時刻での &'，#'，"'の積分だ

抜けていく．したがって，光子がただ移流する自由流（freeけで表されることを意味している．つまり，式（８），（９）-stream）の状況となる．この場合，問題は非常に簡単にから "'や !'が求まるのであれば，輻射流体計算には #'なったように感じられる．確かに，まったく放出以外の相
４６９

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

互作用がなければ，単に輻射冷却が起こるだけあり，輻射

ルは量子力学的な準位間の遷移過程によって決まるので，輸送方程式を解く必要すらなくなる．つまり，図１の上と電子の束縛状態間の遷移はいわゆる線スペクトルとなっ右の要素だけで系を表現することができる．ところが，流て，エネルギー空間における変動が激しい．エネルギー空体の挙動に無視できないぐらいの吸収が起こるとなると，間の粗視化を行わず，特定の光子エネルギーへのすべてのやはり何らかの形で輸送方程式を解く必要が生じる．光学輻射遷移過程を逐次評価する手法を line-by-line 法とよぶ的に薄い場合，流体の時間スケールでは伝搬速度の大きいが，この手法は輻射スペクトルの詳細を知りたいときには輻射はほぼ定常になるため，輸送方程式の時間微分項を省有効であるものの，膨大な計算時間を要するため輻射流体

略できる場合が多い．ただし，非等方性が重要となるため，計算には不向きである．
そこで，エネルギー空間をある幅を持った領域で分割方向の次元を落とすことはできない．
１．
３．
２ LTE と non-LTE

し，その領域ごとに輻射輸送方程式を解くといった方針が図１に示したように，輻射流体計算において，流体と輻取られる．このようにエネルギー空間を有限のエネルギー

射輸送の橋渡しをするのが原子過程で決まる物質の輻射特群（bin ともいう）に分割することを多群近似，あるいはマ性である．物質が局所的に熱平衡にあるとする LTE の場

ルチバンドモデルとよぶ［８］．分割数が多いほど実際のス

合，物質の熱力学的状態は密度と温度に従い，その状態にペクトルに近づけることができるが，それに比例して計算

よって輻射特性も決定づけられる．これは，プラズマが比量も増大するため，不必要に分割すべきではない．また，

較的低温で高密度のときに成立する．ところが，レーザー

輻射流体シミュレーションは，通常，流体力学的な挙動にプラズマのように高温で低密度な領域が含まれる場合や，注目している場合が多く，輻射スペクトルの詳細まで再現

物質の状態が瞬時に変化するような 状 況 で は 必 ずしもすることは求められない．実際，エネルギー空間の解像度

LTE が成立しない．このような状態を非局所熱平衡（nonが低くても流体の挙動は十分予測できるため，分割数は通-LTE）という．

常数十のオーダーで，多くても百前後であることが多い．

LTE でなければ non-LTE な の で，一 口 に non-LTE と

ただし，特定のスペクトルがエネルギー輸送に大きく貢献

いっても様々な状態が考えられる．最もLTEから遠い状況しており，かつそのスペクトル構造が輸送に大きく影響す下では輻射特性に寄与する束縛電子の量子力学的状態が統る場合はその限りではない．

計的分布に従わないため，時々刻々状態遷移過程を追跡す

ちなみに，まったくエネルギー空間の分割を行わず，平る必要がある．しかし，多くの場合，ある熱力学的条件下均化された放射係数や吸収係数を用いて輸送計算を行うこでの定常状態を仮定できるため，non-LTE といっても時間とを灰色（gray）近似とよぶ．特に，光学的に薄い場合に依存のないある種の平衡状態を表す．これを準定常状態

はプランク平均，光学的に厚い場合にはロッセランド平均（quasi-steady state: QSS）という．最もよく仮定されるQSSとよばれる平均化が施される［１］．

は衝突輻射平衡（collisional radiative equilibrium: CRE）あ１．
３．
４ 様々な解法
さて，上述のような輻射輸送方程式の性質を考えながるいは衝突輻射定常状態（collisional radiative steady state:ら，いくつかの代表的な解法について概観してみる．

CRSS）と呼ばれる状態で，電子衝突過程として電離・再結合および励起・脱励起が起こると同時に，輻射過程とし

まず，モンテカルロ法はおそらく最も古くから輻射輸送ては再結合・脱励起だけが起こると仮定し，それらが全体方程式に用いられている解法の一つであろう．モンテカル

としてつり合っている状態である［７］．比較的高温である

ロ法では，サンプル粒子となる光子をランダムに発生さ程度密度が高いときにこの仮定は成り立つ．一方，アブせ，それに確率論的な衝突イベントを繰り返し与えること

レーションプラズマのように高温かつ低密度な場合は衝突で，輻射場の発展を追う［９］．光子発生の試行回数が多い

過程がほとんど起きないので，基底状態からの電子衝突とほど正確な解を得ることができるが，逆に少ないと統計的輻射遷移のみによって QSS になっており，これをコロナ平な誤差が生じるため，ノイズの多い解となる．モンテカル

衡とよぶ．

ロ法は，多次元問題の解の概形を得るには適しているが，流体計算や輻射輸送計算に使用する放射係数や吸収係数

高精度の解を得るには膨大な数のサンプル粒子が必要となは，密度や温度といった熱力学的諸量のみによって決定さる．また，光学的に厚い場合，衝突が頻繁に起こるため，自れなければ，膨大な数の状態遷移を記述する大規模な連立

由行程を衝突区間にとると光子がほとんど空間を移動でき

レート方程式をインラインで解く必要がある．そのため，ないという問題も生じる．そのため，どちらかというと光

通常は上記のような QSS を仮定し，輻射流体計算に必要な学的に薄い場合に適した解法である．また，散乱が容易に物理量を密度と温度の関数として予めテーブル化してお導入できるのも特徴の一つである．
非等方性の強い輻射場，つまり光学的に薄い系の解法とく．それらの具体的な計算手法については，本講座の第３章を参照されたい．

して，現在最もよく用いられるのは通称 "! 法とよばれる

１．
３．
３ 多群近似離散座標（discrete ordinate）法であろう［１０］．これは輻，つまり光子エネルギー #!の関数
輻射強度 !
!は周波数 !

射輸送方程式を伝搬方向も含めて離散化して解く手法であである．したがって !
!は光子エネルギー空間にも次元をる．ただし，当然ながら他の手法に比べて演算量が多い．

もっており，数値計算するためには，自然な発想としてこまた，方向次元の離散化が不十分であると，遠くの発光源の空間についても離散化することを考える．輻射スペクトが捉え難くなるいわゆる ray effect と呼ばれる問題が生じ
４７０

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishi

る［１１］．ただし，散乱も考慮することができ，モンテカル
ロ法のような統計ノイズもないため，計算資源さえ許せば

f %%広範な問題に適用することができる．この手法の応用例に

P%
$
!%

ついては，第２章で詳しく説明する．

%"
%ΩΩdΩ

（１６）

%"
dΩ
%

もしこの f %を何らかの方法で決めることができるなら，P %解の定常性を仮定することができれば，最も簡単な手法
は光線追跡法であろう．この手法では，解を必要とする場を!%で表すことができるので，式（８），（９）からなる方程所に向かって，外部境界から輸送方程式を光路に沿って積式系を閉じることができる．これは，流体の方程式におけ分するだけである．ただし，この手法でも #% 法と同様にる状態方程式の役割に似ており，モーメント方程式に基づ

ray effect が生じるので，なるべく多くの方向に対して積分

く数値解の信頼性はエディントン因子の決定方法に大きくを実行する必要がある．局所的な発光が重要になる場合，依存する．

それを遠くで捉えるためには膨大な光線数を必要とする．

エディントン因子の決定方法の詳細については後の議論

また，そのままでは散乱を扱うことが難しいため，定常かに譲ることとして，まずは f %が何らかの方法で得られた場

つ散乱の無視できる系で有効である．

合に，モーメント方程式を如何に離散化するかを考える．
P %＝f %!%を式（９）に代入すると，以下の式を得る．

最後に，本章の後半で詳しく説明するモーメント方程式に基づく解法は，光学的に厚い系で用いられる手法であ

!'!%
!
"$&#
'
$!#%!%%!%($'(

る．輻射輸送方程式の方向モーメントを取った方程式を扱うため，方向の次元がない分，演算量を大幅に削減できる．

（１７）

こ の と き，!%の 時 間 微 分 が 無 視 で き る と す る と，

ただし，非等方性を捉えることが苦手であり，それを克服するためには後述する付加的な手続きを必要とする．輻射

（８）から以下の
!
!%$!'
$!
#%(
&#
'
%!%(と表せるので，式と流体との相互作用が強い領域は光学的に厚い場合が多い

式を得る．

ため，輻射流体シミュレーションでは安定で演算量も少な

$
'
!%
!
&$
#%!%
!&#
'
! &#
"$$
%!%(
%!$
#%
'
(

いモーメント方程式に基づく手法が最もよく用いられる．
以上，これらの手法をまとめると表１のようになる．光（１８）

学厚さの程度，散乱の重要性に加え，何らかの解を得よう

この式は!%に対する拡散方程式になっており，一般的な拡とするならば，計算精度の重要度，計算資源の限界等を考

散方程式の解法を用いて数値解を得ることが可能となる．

慮した上で計算手法を選定する必要がある．

また，前節で説明したように光学的に厚い系では輻射場が拡散的になるため，上式はそのような状況下で成り立つ式

１．
４ モーメント方程式の解法

で あ る と 考 え ら れ る．特 に，光 学 的 に 十 分 厚 け れ ば，

ここから本章の残りを使って，モーメント方程式（８），

I（I は単位テンソル）となって，問題は最も簡単化
f %$'
!
!
#(

（９）の解法について解説する．前節で述べたように，輻射される．例えば，平板一次元であれば有限体積法とオイ

流体シミュレーションにおける輸送計算で最もよく用いら

ラー陰解法によって次のように離散化できるであろう．

れる手法である．モーメント方程式は元の輻射輸送方程式'' '"!
!'"!!%!
!&'"!!!&' $ %& &! #&"! &"!
'")"(#&"!!"")

よりも方向ベクトルの分だけ次元が少ないが，その代わり
未知変数が "
%から !%，!%，P %の３つに増えている．ところがモーメント方程式は２つであるから，このままでは解

' '"!
'
%
! !%
!'"!
'
' '"!
!& （１９）#&&$!$$$! & & ' &!! &!! $&#&!!!"")を得るための式が足りない．式の数を増やすためにさらに高次のモーメントを取ったとしても，また新しく未知変数が追加されるため，方程式系を閉じることができない．そここでは添字の %は省略し，等間隔格子を仮定した．簡単

こでいわゆるクロージャーとして，次のエディントン因子のために $，#，%
はすべて 'ステップでの値を用いたので，（テンソル）を導入する．

上式は !&'"!に関する連立一次方程式になっており，それを SOR や ICCG などの連立一次方程式の解法を用いて解く

表１

ことで !&の時間発展を得ることができる［１２］．ただし，

輻射輸送計算手法の比較．これらの評価はあくまで筆者の経験によるものであり，数値的裏付けはない．それぞれに提案されている改良法やプログラミング技法によっても相対的評価は異なってくるが，各々の手法の基本的なアイデ
アにそうような評価にした．

採用する解法によって解を得るまでの計算時間が大きく変わるため，適切な解法の選択が必要である．
式（１６）からわかるとおり，エディントン因子は輻射強度"%で定義されるため，定義どおりに求めようとすると "%moment eqs.
ray trace Monte Carlo#%computational cost
low
high
high
mediumrequired memory
small
large
large
mediumanisotopic field
poor
fine
fine
fine
opaque fieldfine
medium medium
poor
scattering field
N/A
available complicated available
time-dependence available available
N/Aavailableのために元の輻射輸送方程式を解かなければならない．したがって，何らかの仮定の元に近似的に決定すべきである．特に，流れ場の状態に従ってエディントン因子を変化させる手法を可変エディントン因子（variable Eddington
factor: VEF）法という．エディントン因子の決定方法は大きく分けて２つある．１つは局所的な物理量のみによって
４７１

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

決める方法であり，もう１つは非局所的な影響も取り入れ

光学的に薄い場合であり，*!'*"'となっていると考え%'る方法である．光学的に厚い状況では，光子が遠くまで衝られる．したがって，式（１７）は実質的に移流方程式になる．

突なしに進むことができないため，輻射場は局所的な情報

このときの移流速度は光速で決まり，陽解法で安定に計算で決まり，前者の取り扱いが適用できる．局所的な取り扱するためには #+を小さく設定する必要がある．そのため

いでは，あまり多くの演算を必要としないため，効率的な離散化は陰的に行われるべきであり，例えば平板一次元で計算が期待できる．しかし，光学的に薄い場合は輻射場があれば以下のような離散式が考えられる［１６］．

遠くの光源によって決まる可能性があるため，後者のよう

""
*""
!#)*
*
* *"" #
#
#
"""
!""
)
"
(%
$
")*""#")*"#+!%
!'
!
"
)
) )
#,
（
２１）
#
#)*
'
#
+
*""
#
"""
#)"""##!
**$)$)#+ ""'
#+
""'
#
#
"""
"""

な取り扱いが必要となる．
１．
４．
１ エディントン因子の局所近似法
エディントン因子を局所的に決定する方法は，これまでいくつも試みられている．基本的には，光学的に厚い限界

*
* *""
(
"*""!("!)"" )"" ) ) "（２２）#,で拡散近似に，薄い限界では自由流になるように構築される．ただし，光学的に薄い場合にエディントン因子を局所的に決定可能な根拠は原理的に存在しないため，決定的な

ここで添字 'は省略し，")，#)"""#はスタッガード格子で定手法はない［１３］．したがって，光学的に薄い領域での解が義している．このような離散化であれ ば，結 果 と し てそれほど重要でないが，厚い領域と混在しているような場")*""の連立一次方程式に帰着させることができる．

合に，安定かつ効率的に両方の領域をカバーするための手法とみるべきであり，以下に挙げる様々な方法から実際に

ただし，%" 法では光学的に薄い限界で伝搬速度が '"$+に漸近することが指摘されており，そのような領域では 'どれを使用するかは，そのような観点で選択されるべきでに正しく漸近する %""$法が好ましいとされている［１７］．

ある．

最大エントロピー法

流束制限拡散法

ボーズ−アインシュタイン統計に従う光子がエントロ

流束制限拡散（flux-limited diffusion: FLD）法は，エディ

ピーを最大にするように分布関数を求め，そこから決まる

I としたときのントン因子の決定法というよりは f '#(""$)

エディントン因子を用いる方法が Minerbo の最大エントロ式（１８）の近似解法である．エディントン因子が一定の値をピー法である［１８］．輻射場の軸対称性を仮定すると，エ

取 る と き，式（１８）に お け る エ ネ ル ギ ー 流 束 はディントン因子は以下のように表される．

!'#!(
'"
$')
$
'"'で表されるため，光学的に薄くなって$'が小さくなったり，"'の勾配が急峻になると，流束が光

""
"$
""
(
!(
)
!")
f '# (
I" (##速から決まる限界値（free-stream limit）である '"'を超え

ただし，"#!'"!'*である．Minerbo の原論文では， (*る可能性がある．その場合，非物理的な流束となるため，

+#*
は#
!'"
'
"'*の関数として与えられており，そのようなこれを超えないように流束を制限する必要がある．

表現は後述の $" 法にも使用することができる．

制限方法はいくつか考えられるが，例えば以下のような

一方，一種のFLDとして使用されることも多く，式（２０）ものが挙げられる．
'!'#!
'"'#!!)*('"'
$'"*
"
"'
$
'"'*

（２３）

と 同 様 に !'が "'の 勾 配 で 表 現 さ れ る と す る と， (は"$'"')をパラメータとして
&#*
(
'"'*

（２０）

(
#&"&#&#

ここで，!)*( は流束制限された拡散係数である．一見，こ（２４）

れはエディントン因子を操作しているようにもみえるが，と表される［１９］．このとき，"#'"'"*'"'*である．また，それとは若干思想が異なる．本来，一次元で考えた場合，&は流束制限関数であり，最大エントロピー法では以下のエディントン因子は光学的に厚い領域から薄い領域に向けようになる．
!
"&#) if !$&$$#$"#'""#
(
"
#
&(
&)
#$
"
!"
&
"
"
)
(
"&
"
"#
+
if
"
!&!&
$
#
$
#

て 1/3 か ら１に 変 化 す る よ うに 構 築 さ れ る べ き で あ る［１４］．拡散近似では輻射場はあくまで等方なままなので，いくら流束を制限したところで，非等方な自由流を表現す

Levermore-Pomraning 法と M1 法

ることは難しい．また，上記の表現では光学的に厚い領域
でも０次精度しか達成できないという指摘もあり，他の流輻射場の角度依存性を示す規格化した輻射強度を定義

束制限法が提案されている［１５］．

P 1 法と P 1/3 法%" および %""$法では，エディントン因子を一定とする
が，代わりにエネルギー流束の時間発展方程式も解くこと
で，エネルギー流束が流束限界を超えないようにする．例
I，%""$法では f '#I とする［１３］．
えば，%" 法では f '#(""$)
エネルギー流束は式（１７）の時間微分を残したまま解くことによって求めるが，時間微分項が重要になってくるのは（２５）

し，それが時間的，空間的に十分緩やかに変化するという仮定に基づいてエディントン因子を決定することもできる［２０］．この Levermore-Pomraning（L-P）法も FLD として使用されることが多く，&#*"$'"')をパラメータと('"'*すると，以下のような &の近似式を使って式（２３），（２４）からエディントン因子を求めることができる［１４］．
&(
&)
#
４７２

#"&
&"&#
&"$

（２６）

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishi

リッドと捉えることができる．
最大エントロピー法と L-P 法で評価した )を図２にそれ

光線追跡法によるエディントン因子の評価

ぞれ示した．'は光学厚さの程度を示しており，光学的にエディントン因子は流体の時間スケールにおいて定常で厚いとき 'は小さくなり， )は 1/3 に近づく．光学的に薄あると仮定すると，エディントン因子を評価するための輸い場合は 'が大きくなり，)が１に漸近することがわかる．

送方程式の時間微分項を無視できる．そこで，以下の定常

ここで，)は式（２３）中の変数であるが，一次元のときにはな輸送方程式を光線追跡法で解くことを考える．

これがそのままエディントン因子になる．

&#
Ω#
%#
$$
$$#
$!"
$#
$
&*

Levermore はこれとは別に，現在 $! 法のクロージャー
としてよく用いられるエディントン因子の決定法も提案し（２８）ている［１９］．$! 法は主に宇宙物理における輻射流体計算

ただし，*は光路に沿った距離を表し，簡単のため散乱項で採用されており，式（８），（１７）を !$と !$から求まるエは省略した．光線追跡のアルゴリズムは単純であり，#$

デ ィ ン ト ン 因 子 を 用 い て 解 く 手 法 で あ る．こ の 場 合，
'$&
!$!
"
(
!$&を用いて，は場所と方向の関数であるから，エディントン因子を得た)$

い場所 &からある方向に遡る光路と外部境界との交点 %を求め，点 % から点 &に向かって上式を積分すれば点 &

'
"
#"$

"

'
"
%""!$!#

"

での #
$や "
$はセルで定義される値$が求まる．このとき，#（２７）なので，通過するセルをどれだけ横切るかを逐次評価する必要がある．

という式からエディントン因子が求められ，光学的に薄い

上 記 の よ う な 計 算 は long characteristics 法 と よ ば れ

状況下の検証も行われている［２１］．
１．
４．
２ 非局所的に決定するエディントン因子［２２］，手続きとしては単純であるが，エディントン因子が局所的にエディントン因子を決定する場合，遠くの光源

必要なすべての格子点で，ある程度の方向離散点の数だけの影響をうまく取り入れることができない．モーメント方

光線追跡を行わなければならず，演算量はたちまち膨大な

程式に基づく解法は本来そのような状況に適していないものとなる．そこで，図３に示すように，光線追跡は隣りが，これを克服するために非局所的なエディントン因子ののセル境界からのみ行い，初期値となる境界での輻射強度評価が試みられている．

は周辺の値を補間することによって評価する short characteristics 法が提案されている
［１６，
２３］．

式（１６）の定義に立ち戻ると，エディントン因子を評価するためには #$が求まっていればよい．ところが，そもそもshort characteristics 法によって演算量を抑えることは#$は輻射輸送方程式を解かなければ決められないので，こできるが，光源から遠くにいくほど非等方性が失われてしれを用いてエディントン因子を定義するのは本末転倒であまう問題がある．一方 long characteristics 法では，演算量ろう．しかし，計算負荷を小さく抑えながらある程度の輸を厭わなければ精度よくエディントン因子を得ることがで

送方程式の解が得られるのならば，その解からエディント

きるが，有限の方向離散点ではいずれ遠くの光源を捉えらン因子だけを評価し，時間発展はモーメント方程式を用いれなくなるという光線追跡法の根本的な問題は避けることて評価するという手続きが可能となる．これは，エディン

ができない
［２４］．

トン因子を評価するための輻射輸送計算にはそれほど精度を求めず，精度が重要となるエネルギー輸送の評価に関してはモーメント方程式の解法に委ねるという発想から導かれており，連続体的なアプローチと運動論的手法のハイブ

図３

図２

エディントン因子の近似法の比較．実線は最大エントロ
ピー法で破線は L-P 法．一次元のときのエディントン因子となる )は光学厚さを表すパラメータ R = | %E | / ( "E ) に対して 1/3 から 1 まで変化する．

４７３

光線追跡法によるエディントン因子評価法の概念図．long
characteristics 法では，ある方向の輻射強度が点 P で必要な場合，そこから光路を遡って外部境界との交点 O を求める．そして，点 O から P に向かって，輸送方程式を積分する．一方，short characteristics 法で点 P´の輻射強度を評価する場合，隣のセル境界との交点 O´から P´の間のみ積分する．積分の初期値となる点 O ’における輻射強度は，周辺の輻射強度から補間して求める．

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012［７］D. Salzmann, Atomic Physics in Hot Plasmas (Oxford University Press, Oxford, 1998).
［８］D.E. Cullen and G.C. Pomraning, J. Quant. Spectrosc. Radiat. Transfer 24, 97 (1980).
［９］Y. Wu, M.F. Modest, and D.C. Haworth, J. Comput. Phys.
223, 898 (2007).
［１０］小林啓祐：原子炉物理（コロナ社，東京，１
９
９
６）
．
［１
１］P.J. Coelho, J. Quant. Spectrosc. Radiat. Transfer 73, 231(2002).
［１
２］小国力 編著，村田健郎，三好敏郎，J.J.ドンガラ，長谷川秀彦著：行列計算ソフトウェア−WS，スーパーコ
ン，並列計算機（丸善，東京，１
９
９
１）
．
［１
３］G.L. Olson, L.H. Auer and M.L. Hall, J. Quant. Spectrosc.
Radiat. Transfer 64, 619 (2000).
［１４］N.J. Turner and J.M. Stone, Astrophys. J. Suppl. Ser. 135,
95 (2001).
［１５］J.E. Morel, J. Quant. Spectrosc. Radiat. Transfer 65, 769(2000).
［１６］J.C. Hayes and M.L. Norman, Astrophys. J. Suppl. Ser. 147,
197 (2003).

モンテカルロ法とのハイブリッド計算
連続体に対する方程式と運動論的方程式のハイブリッド
解法を考えるとき，計算格子ベースの手法と粒子的手法の結合計算が最も現実的な選択である．実際，光線追跡法も一種の粒子的手法であり，そのような発想の範疇を脱しない．ところが，光線追跡法は粒子の飛行経路を固定した計算手法であるため，予め用意した飛行経路にない情報を取り入れることができないという意味では，計算格子ベースの手法と解釈することもできる．すなわち，計算空間中の点光源を隈なく捉えるためには網羅的に飛行経路，つまり離散方向を設定する必要があり，それは計算資源の面から考えて到底現実的ではない．
そこで，もう一つの粒子的手法としてモンテカルロ法を考えることは有益であろう［２５］．前節で述べたように，モンテカルロ法は解の統計誤差を除去するために多くのサンプル粒子を必要とするが，解の概形を得るだけであれば，比較的少ない計算負荷で目的を達成することができる．エディントン因子を評価するだけであれば，多くの計算負荷［１７］K.H. Simmons and D.Mihalas, J. Quant. Spectrosc. Radiat.
Transfer 66, 263 (2000).
［１８］G.N. Minerbo, J. Quant. Spectrosc. Radiat. Transfer 20,
541 (1978).
［１９］C.D. Levermore, J. Quant. Spectrosc. Radiat. Transfer 31,
149 (1984).
［２０］C.D. Levermore and G.C. Pomraning, Astrophys. J. 248,
321 (1981).
［２１］D. Aubert and R. Teyssier, Mon. Not. R. Astron. Soc. 387,
295 (2008).
［２２］梅村雅之，中本泰史：ながれ 15, 457 (1996).
［２３］J.M. Stone, D. Mihalas and M.L. Norman, Astrophys. J.
Suppl. Ser. 80, 819 (1992).
［２４］N. Ohnishi, K. Sugai and Y. Ogino, J. Phys. Conf. Ser. 244,
022078 (2010).
［２５］N. Ohnishi, Toward an accurate numerical simulation ofradiation hydrodynamics in laser ablation plasmas,
Proc.9th International Conference on High Energy Density
Laboratory Astrophysics, 2012.

を必要としない可能性が高い．発生確率を放射係数で重み付けすることで，光源からより多くのサンプル粒子を生成し，点光源の寄与をなるべく全域に行き渡らせることも可能である［９，２５］．

参考文献［１］Y.B. Zel'dovich and Y.P. Raizer, Physics of Shock Waves andHigh-Temperature Hydrodynamic Phenomena (Dover, New
York, 2002).
［２］K. Kotake, K. Sato and K. Takahashi, Rev. Prog. Phys. 69,
971 (2006).
［３］S. Atzeni and J. Meyer-ter-Vehn, The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot Dense
Matter (Oxford University Press, Oxford, 2004).
［４］D. Mihalas and B.W. Mihalas, Foundations of Radiation Hydrodynamics (Dover, New York, 1999).
［５］G.C. Pomraning, The Equations of Radiation Hydrodynamics(Pergamon Press, Oxford, 1973).
［６］J.I. Castor, Radiation Hydrodynamics (Cambridge University Press, Cambridge, 2004).

!!!!!!!!!!!!!!!!!

おお

にし

なお

ふみ

大 西 直 文東北大学大学院工学研究科航空宇宙工学専攻准教授．２０
０
１年３月大阪大学大学院工学研究科電気工学専攻にて学位を取得（博士（工学））．専門は輻射流体力学，高温気体!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

力学．出身は島根県で大学は大阪，就職は関東を飛び越して
仙台，と東にシフトし，おかげで大きな地震に２度遭いました．帰省が大変で，毎年盆には家族を乗せて仙台から島根を走破します．関西より西の道路で関東ナンバーを見ると，勝ったと思ってしまいます．

!!!!!!!!!!!!!!!!!!!!!!!!!!

４７４

J. Plasma Fusion Res. Vol.88, No.10 (2012)６
１０‐６１７

講座 輻射流体シミュレーション

２．非等方性の強い輻射場における輸送計算：
超新星爆発におけるニュートリノ輻射輸送の例
住吉光介
沼津工業高等専門学校・教養科物理（原稿受付：２
０
１
２年８月２
５日）

輻射流体シミュレーションの記述において，輻射のエネルギー・角度分布が重要となる場合について，宇宙物理における重力崩壊型超新星の問題を例として概説する．重い星の進化の最期には，超新星爆発という華々しい天体現象が起こる．その爆発メカニズムの解明では，高温高密度の極限状態におけるニュートリノ粒子の反応・輸送過程が重要な役割を担っている．超新星においてニュートリノ輻射輸送の緻密な記述が果たす役割を解説して，ニュートリノ輻射輸送計算による発展について概観し，多次元計算における近似手法から６次元ボルツマン方程式の直接解法による最新の計算手法までを紹介する．
Keywords:
radiation transfer, Boltzmann equation, neutrino, supernova

２．
１ 宇宙物理分野における輻射輸送

しまう．こうした天体現象は，重力崩壊型超新星と呼ばれ，

輻射流体計算は，宇宙分野においても多く現れ，様々な

古くから観測例があり，年間を通じて数多くの観測がなさ

天体のダイナミクスにおいて重要な鍵を握っている．銀河れている．なかでも有名なものとして，例えば１０５４年に観や星などの形成過程においては，輻射によるエネルギー・

測された記録が存在する超新星爆発があり，その残骸であ

運動量の輸送が流体への影響を与えて，天体の構造や進化るカニ 星 雲（中 心 部 に パ ル サ ー）が 知 ら れ て い る．まが変わってくることになる．一方，星が進化して，やがてた，１９８７年に観測された超新星爆発は，宇宙からのニュー

最期にどうなるのか，というのも根源的な疑問である．重トリノ粒子を初めて検出したことが日本人のノーベル物理

たい星はやがて潰れてしまい，大爆発（超新星爆発）を起学賞の受賞へと繋がり，一般社会へも広く知られることとこすことが知られている．この爆発現象の鍵を握るのが，なった．

輻射輸送である．ただし，この場合は，ニュートリノとい爆発のあと，中心には高密度天体（中性子星あるいはブう素粒子がダイナミクスにおける輸送過程を担っている．

ラックホール）が残されて，外層においては爆発的な元素

このニュートリノ輻射輸送と呼ばれる過程を扱うことは，合成が行われて外部へまき散らされる．これらの物質放出

宇宙物理分野における重要な課題の一つであり，様々な計がもとで，新たな星が生まれる材料となり，星の集まりで算手法が開発されて，数値シミュレーションにより爆発メある銀河の進化にも影響を及ぼしている．我々の経済・生カニズムの解明が行なわれてきた．特に，現象記述には，

活に関わりのある貴金属（プラチナ・金など）や原子力エ

エネルギー・角度分布について緻密な取り扱いが要求されネルギーの燃料であるウランなどは超新星爆発のような爆るため，様々な近似を越えて，より厳密な計算手法の開発

発的な現象により作られたとされているが，まだ起源はよが行われてきた．本章では，超新星爆発におけるニュート

くわかっていない．

リノ輻射輸送流体計算を例に取って，ボルツマン方程式を

このように宇宙・物質の発展の源である重力崩壊型超新扱う計算手法の最近の発展について解説を行う．

星であるが，その爆発メカニズムは，４０年以上にわたる精力的な研究にも関わらず，未だ爆発の引き金となる根源的２．
２ 超新星爆発とニュートリノ

な部分が解明されていない［１，２］．実際のところ，超新星爆２．
２．
１ 重力崩壊型超新星とは発の数値シミュレーションにおいて，星の重力崩壊から爆

太陽の質量の１０倍を越える質量をもつ重たい星は，進化発を再現することは難しく，爆発した場合であっても起源の最期に自らを支えることができず，重力崩壊を起こしてとなるメカニズムを確定できずにいるのが現状である．こ潰れてしまう．高温高密度になった中心コアで跳ね返りの難しさの原因は，超新星コアにおける極限状態での物質（バウンス）が起こり，発生した衝撃波がやがて星全体を吹の性質とニュートリノ輻射輸送の様相が絡みあって，流体

き飛ばし，きわめて明るい輝きを放ち，やがて暗くなって

ダイナミクスが複雑になっていることにある．中でも，

2. Computations of Radiation Transfer with Anisotropic Distributions: Examples of the Neutrino Transfer in Supernova Explosions
SUMIYOSHI Kohsuke
author’s e-mail: sumi@numazu−ct.ac.jp

６１０!2012 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

2. Computations of Radiation Transfer with Anisotropic Distributions: Examples of the Neutrino Transfer in Supernova Explosions

K. Sumiyoshi

ニュートリノ粒子が高温高密度の物質において，どのよう考えられている．大まかにエネルギーの移り変わりを見るに相互作用して伝搬したり輸送されるのかが問題の鍵となと，重力崩壊によって鉄コア（太陽質量の 1.4 倍程度の質る．

量）の半径が小さくなる（〜103 km から〜10 km へ）こと

２．
２．
２ 超新星爆発メカニズム

により重力エネルギーが解放される．このエネルギーは，まず，爆発に至るまでの超新星ダイナミクスを簡単に説明

1053 erg 程度であり，圧縮による温度上昇に伴い熱エネル

したい（図１）．星の中心にある鉄コアは重力崩壊により圧ギーに転換されて，その後はニュートリノ対の生成へと使縮されて，すぐに中心での密度は3×1014 g/cm3（原子核物

われる．ニュートリノのほとんどは外に逃げてしまい，超11質密度）を越え，温度は 10 K にも達する．この密度は，通新星ニュートリノとなる．地上の検出器により測定された常の原子核内部で陽子・中性子が安定に存在する状態を越ニュートリノのエネルギーと個数から，ニュートリノ放出

えている．密度上昇により，核力による斥力が急激に働きの合計エネルギーは 1053 erg であることが判明したため，だすため，内部コアはこれ以上圧縮できず，堅くなり反発星の重力崩壊によるエネルギー解放という大枠のシナリオ

してコアバウンスする．引き続き外部コアは重力により落は正しいことが明らかとなった．

下してくるので，この両者の間に衝撃波が発生する．バウ

一方，外層が吹き飛ばされる際の爆発エネルギーはンスによる外向きの衝撃波が，鉄コア表面まで突き抜ける

1051 erg であることが観測からわかっており，このエネル

ことができれば，星の最外層に達して，超新星爆発となる．

ギーをどのように説明するのかが課題である．最近の数値

ここで重要な役割をするのがニュートリノである．重力シミュレーションによる結果では，コアバウンスによる衝崩壊の際に，原子核が電子捕獲反応を起こしてニュートリ

撃波は，伝搬途中のエネルギー損失のため，途中で停滞しノが発生するほか，高温部分では熱的にニュートリノ対がてしまうことがわかっており，何らかのエネルギー的な手大量に発生する．これらのニュートリノは，中心密度が高助けが必要である．放出されるニュートリノのうち，一部

いために，すぐさま逃げることはできず，中心コアに閉込

は物質を加熱することに使われて，衝撃波を後押しするこ

められた状態にある（図１ニュートリノ閉じ込め）．物質と（ニュートリノ加熱）もわかっており，この量が十分かと頻繁に反応・散乱しながら，やがて外へ向けて飛んでいを明らかにすることが重要な課題である．つまりニュート

く．こうして放出された大量のニュートリノの一部を，地リノ放出分のうち，約１％のエネルギーを物質に与えて，球上で検出したものが，先に述べた超新星ニュートリノで爆発エネルギーに寄与するか否かを明らかにするため，ある．

ニュートリノ輻射エネルギーの一部が移行される過程を扱う，緻密な輻射輸送計算が必要不可欠ということである．

このニュートリノ発生・放出過程は，エネルギー移行の観点からも重要であり，その一部が爆発に寄与していると

２．
３ 超新星爆発の数値シミュレーション
２．
３．
１ 必要な物理過程
爆発メカニズムを解明すべく，重力崩壊型超新星のダイ
ナミクスを記述するためには，大枠となる流体力学とニュートリノ輻射輸送の記述，高温高密度における物質の熱力学的性質（状態方程式）とニュートリノ・物質間の相互作用（反応率）の詳細なミクロ物理データが必要である．
超新星の構造・ダイナミクスは重力の支配下にあり，自然界の４つの基本的な力（重力，弱い相互作用，強い相互作
用，電磁相互作用）のすべてが関与する系となっている．
この中で，流体計算自体については，計算法や計算機資源の発展により，３次元空間であっても高い解像度での計算が可能になってきた．しかし，ニュートリノ輻射による流体ダイナミクスへの影響（圧力の寄与，加熱・冷却）が大きく，流体の状態変化によりニュートリノ反応・輻射が変動するため，流体とニュートリノ輻射の両者を同時に解くことが必要である．さらに，ニュートリノ輻射を求めるにあたっては，後述のように単純な近似を用いることができない．特に，ニュートリノ反応率は，ニュートリノの種
類・エネルギー・角度に依存しており，個々の反応過程の詳細を計算に取り入れる必要がある．また，実験室では到
図１

鉄コアの重力崩壊から超新星爆発に至るまでの概略．中心
密度が上がるにつれて，ニュートリノは逃げ出すことができなくなり，物質中に閉じ込められる．中心コアの跳ね返りで衝撃波が発生して，外層部分を吹き飛ばす．中心には高密度天体が残される．

達できない高温高密度領域での物質の状態方程式，ニュートリノ・物質の相互作用は，原子核物理の手法を駆使して理論的に求めて，広範囲をカバーする核データとして整備しなければならない．
６１１

Journal of Plasma and Fusion Research Vol.88, No.10 October 2012２．
３．
２ ニュートリノ輻射の役割
ニュートリノ輻射は，超新星のダイナミクスの各ステー
ジで重要な役割を果たしている．なかでも爆発か否かを決めるのに重大な局面となるのはコアバウンスで打ち上げられた衝撃波が，途中で停滞してしまった状況である．コア
バウンス後 100 ms 程度で衝撃波は 100 km を越えたあたりで停滞して，外へは進まなくなってしまう．伝搬の際にエネルギー損失（鉄を核子に分解するためにエネルギーが消費される）があり，コアバウンス自体による初期の衝撃波
エネルギーは使い果たしてしまうためである．この時点
で，中心には（のちに中性子星となる）原始中性子星が誕
生し，密度・温度が高くニュートリノを多く含んだ状態にあり，表面からニュートリノが徐々に外へ流れ出している．ニュートリノは物質からエネルギーを持ち出して伝搬

図２

コアバウンス後 100 ms における速度（停滞衝撃波）と加熱率（負では冷却）の分布を半径の関数として示した．球対称のニュートリノ輻射流体計算では衝撃波が 100 km 付近
で停滞してしまう．原始中性子星の表面ではニュートリノ
が放出されて冷却領域となり，ニュートリノの一部は衝撃波付近の物質に吸収されて，その付近において加熱領域を
構成する．

図３

コアバウンス後に衝撃波が停滞した頃のニュートリノ輻射
輸送の様子．中心にはニュートリノを多く含み高温高密度の原始中性子星が誕生している．拡散領域で徐々に放出されたニュートリノは，低密度の外層では自由伝搬となる．
停滞衝撃波付近ではニュートリノ加熱が重要となるが，この辺りは拡散でも自由伝搬でもない中間領域である．光学
厚さは桁で大きく変わり，光学的に厚い領域から薄い領域まで存在する．

する の で，物 質 に と っ て は 冷 却 過 程 と な る．こ れ ら のニュートリノが，外部へ逃げ出す途中で物質に吸収されると，ニュートリノの持つエネルギーが物質へ与えられることになり，物質にとっての加熱過程となる．これらの冷却・加熱領域がどのような配置となるのか，そしてエネル
ギー移行量はどれくらいなのかを定量的に求めなければならない．図２には，バウンス後の速度および加熱率の典型的な分布例を示した．150 km付近で停滞した衝撃波の付近にニュートリノ加熱領域が存在する．その内側はニュートリノを放出する冷却領域である．
このエネルギー移行過程を探る問題が難しい理由の一つ
は，超新星コアの中心から表面までの密度・温度の範囲が非常に広く，輻射輸送の様相が大きく変動していることにあ る．図３に 示 す よ う に，超 新 星 コ ア の 中 心 密 度 は1014 g/cm3 以上，鉄コア表面密度は105 g/cm3 程度であるので，シミュレーション計算でカバーすべき密度範囲は１０桁に及ぶ．中心温度は 1011 K，表面温度は 109 K 程度である．
こうした高温高密度物質の中でニュートリノが相互作用しながら，どのように分布して，輸送されるのかを調べなけ
ればいけない．ここで物質とは，ニュートリノ以外の構成
粒子で，陽子・中性子・原子核・電子・陽電子・光子の集まりである．これらの粒子群は強い（核力）相互作用と電磁気相互作用により，十分短い時間スケールで頻繁に反応を し て 常 に 熱 統 計 平 衡 状 態 に 保 た れ て い る．一 方 で，ニュートリノは弱い相互作用によってのみ物質と相互作用をしており，その反応時間スケールが長い場合には平衡が成り立たないため，物質とは分けて，ニュートリノのエネ
ルギー・角度分布を記述する必要がある．つまり，流体力
学で追う部分は物質の部分であり，ニュートリノは輻射の収・放出・散乱により，熱・化学平衡状態に達している．

部分として扱い，両者が結合した輻射流体力学を解くこと

ニュートリノのエネルギー分布は，物質の温度と化学ポテとなる．

ンシャルで決まる，Fermi-Dirac分布と一致しており，角度

超新星コアにおけるニュートリノ輻射輸送は，中心にお分布は等方である．もちろん物質分布は一様ではなく，密ける準平衡状態から外部へ向けて自由伝搬するまで，非常度・温度・組成の空間分布に従ってニュートリノ分布にもに幅広い過程を経由するものとなっている．図４のように勾配ができる．この勾配に従って拡散現象により，ニュー

重力崩壊時には，ニュートリノ閉じ込めにより徐々に数密

トリノの輸送が進む．この領域では，角度分布が等方から度が増えていき，Fermi-Dirac 分布・等方分布に近づいて若干ずれている程度である．外へ向かって密度・温度が下いく．密度温度が高くなった中心部ではニュートリノ・物

がっていくと，ニュートリノ反応率が下がり，もはや平衡質間の相互作用は頻繁に起こっており，ニュートリノ吸を保つことはできなくなる．密度が 1011 g/cm3 程度になる
６１２

Lecture Note

2. Computations of Radiation Transfer with Anisotropic Distributions: Examples of the Neutrino Transfer in Supernova Explosions

K. Sumiyoshi

２．
４ 輻射流体計算シミュレーションの発展
２．
４．
１ 計算手法の発展と爆発メカニズム解明
超新星爆発の研究は，流体力学とニュートリノ輻射輸送の両者を解く計算手法の発展とともに進んできており，新しく近似や次元数をあげた計算により新たな知見が得られる，という歴史をたどってきたといっても過言ではない．
ニュートリノ輻射輸送の記述によっては，爆発に有利／不利に働くこともあり，注意深く近似の影響を評価しながら進めなければいけない．例えば，重力崩壊時になるべく
ニュートリノをたくさん溜め込んでおき，バウンス後にそ
図４

重力崩壊における中心部でのニュートリノエネルギー分布の例．中心で発生したニュートリノは，初期（Initial collapse）に は わ ず か に 溜 ま る 程 度 だ が，中 心 密 度 が1012 g/cm３程度（Neutrino trapping）になると物質内に閉込められ，分布値が大きくなり始める．この時点では，角度分布は等方ではなく，同じエネルギーでも伝搬角度によって値が異なっている．コアバウンス（Core bounce）までにはニュートリノ閉じ込めは十分となり，熱・化学平衡に達して，等方なフェルミ分布に従っている．

のエネルギーを多く放出するとともに，衝撃波付近で物質加熱に十分に寄与することが起きれば，爆発には有利である．これらの要因は，ニュートリノ輻射輸送の記述により左右される事柄であり，近似方法によっては効果を過大／過小評価してしまう恐れがある．
一方，空間次元数も爆発の可否に大きな影響を与える要因である．計算規模の拡大に応じて制約を外していくと，新たな流体力学的不安定性に伴う現象が現れて，爆発メカと，ニュートリノ反応が稀になり，ニュートリノは放出さニズムの本質ともいえる効果が明らかになってきた．球対れて，外へ伝搬して出ていくようになる．これは，光子で称から軸対称計算へ拡張した際には対流の効果が現れて，言うところの光球（photosphere）に対応しており，ニューさらに赤道面対称性の制限を外した計算では，上下非対称

トリノ光球（neutrinosphere：図３における光学厚さ !!!

となる定在降着衝撃波の不安定性が発見されるなど，多次付近）という．これより十分に外部へ到達すれば，非常に元的な発展に応じてニュートリノ加熱効率が変わり，爆発

薄い物質分布においてニュートリノは自由伝搬しており，へ転ずる後押しになりうることが判明した［３］．それでは，

その角度分布は，伝搬方向の前方ピークへ集中している．

軸対称性を外した３次元計算ではどうなっているのか，

しかし，拡散・ニュートリノ光球から，自由伝搬になるま

ニュートリノ輻射の影響とともに調べることが，現在の研での途中では，一部のニュートリノが物質と相互作用し究の焦点となっている［４］．

て，吸収・散乱を受ける．このように超新星コアでは，エ２．
４．
２ ニュートリノ輻射輸送計算手法の発展

ネルギー分布において平衡分布ではなく，角度分布におい超新星におけるニュートリノ輻射輸送の計算は，計算資て等方でも前方ピークでもない，中間領域が存在している．

源の規模に応じて近似計算から厳密計算へと発展してい爆発メカニズム解明において鍵となっているのは，このる．現在までに，１次元（球対称）においてはニュートリ

中間領域である．ニュートリノを放出する領域（原始中性ノ輻射輸送を厳密に解く手法が確立しており，第一原理計

子星の表面）では物質の冷却が進むが，その放出率は物質算的な意味合いで，ニュートリノ輻射流体計算が行われるの組成やニュートリノエネルギーに依って大きく異なるたようになった．Multi-energy group に対する輻射輸送を，め，ニュートリノ輻射の詳細計算が必要である．例えば，ボルツマン方程式により直接解く !" 法あるいは Edding-

一口にニュートリノ光球といっても，相互作用のエネル

ton factorを介したモーメント法により解いており，ニュー

ギー依存性が強いため，その位置はエネルギーごとに違っトリノ輻射輸送方程式と流体力学方程式を組み合わせて，ている．つまり，ある半径から平衡分布でニュートリノを

重力崩壊からコアバウンス，そして衝撃波伝搬が爆発に至放出している，という様な単純化をすることはできず，放るか，までを系統的に調べることが行われている．多くの出ニュートリノのエネルギー・角度分布（スペクトル）おグループによる数値シミュレーションの結果，球対称計算

よび空間分布が必要である．一方，ニュートリノを吸収すにおいては爆発を再現できないことが確実となった［５］．

る領域では物質の加熱が進むが，その吸収率は，ニュート

また，球対称におけるニュートリノ輻射流体の第一原理計

リノ放出スペクトルの詳細を元にして得られるニュートリ

算は，ニュートリノ・原子核物理がダイナミクスへ及ぼす

ノ輻射に基づいて，エネルギーごとに物質分布との反応率

影響の検証や超新星ニュートリノのスペクトル・時間変動を計算して評価しなければならない．この加熱領域においを予測することにも用いられている［５，６］．

ては，ニュートリノは自由伝搬にはなっておらず，ニュー

球対称における近似計算から厳密計算へ至る過程では，

トリノ流束を解析的あるいは近似的に扱うことはできない．

近似法の影響が詳細に調べられて，のちの多次元計算での

つまり，ニュートリノ放出・吸収および爆発への影響を調応用・近似評価に活かされている．近似手法としては，初

べるには，ニュートリノ輻射の全領域における振る舞いを期の light bulb，leakage 法のほか，flux limiter を用いた拡解くことが不可欠である．

散近似が広く用いられてきた．light bulb 法では，あらかじ
め決めたニュートリノ光球から，温度・化学ポテンシャル
により指定した平衡分布を持ってニュートリノが外へ流れ
６１３

Journal of Plasma and Fusion Research Vol.88, No.10 October 2012

出してくる，という単純化を行う．これらは流体力学的不安定性や元素合成過程などを系統的に調べる際に今も用いられている．leakage 法では，ある決められた密度（ニュートリノ光球）より高密度領域には，ニュートリノが閉じ込
められており，（拡散・伝搬時間スケールで決める）Escape 時間スケールに従ってニュートリノが自由伝搬して流 出 す る，と い う 近 似 を 行 う．flux limited diffusion 法
は，大枠の計算方法として拡散近似を取り，自由伝搬領域
においてニュートリノ伝搬速度が光速を超えないように，拡散係数に対して，関数形（flux limiter）により制限を与

図５

えておくものである（第１章を参照）．
拡散近似においては，中心部分での記述には問題がない
が，原始中性子星表面から衝撃波付近までの中間領域では，flux limiter によって与えられる関数形に頼る方法となっている．この拡散近似による結果を，モンテカルロ法

コアバウンス後 100 ms におけるニュートリノ角度因子の分
布：停滞衝撃波（図２）での流体力学量分布の状況でニュートリノ輻射輸送計算を行って得られたニュートリノ分布関2数により角度変数 #$= cos %$の１，
２次モーメント‹#$›，‹#$›を計算した．２つのモーメント量は，等方分布ではそれぞれ 0，1/3，前方ピーク集中では１の極限値を取る．中間の値となる衝撃波付近（〜100 km）の分布が重要である．

や輻射を厳密に解く方法での結果と比較することにより，近似による影響を定量的に評価することができる．例えル化から脱却して，ニュートリノ輻射輸送を組み込む方向

ば，半径 &における（核子一個あたり）ニュートリノ加熱にある．３次元における流体不安定性の探索においては，

率 #$を求めるには，全核子数に対する陽子・中性子の個

light bulb 法などの簡単化によりニュートリノ加熱効果を，角度因
数比 $%，ニュートリノ光度・エネルギー $"$!
!$%取り入れて計算が行われている．拡散近似をベースとして子が関わっており，

ray-by-ray 法と組み合わせた３次元計算が計算機資源の限

"$!
!#"$'
$%
$#&
&)'
#$#""!! # $!
!
"
+
.(
*
),+
!
&
#
&
'
$
%

界まで使って行われて，最初の爆発シミュレーション結果（１）が得られるようになった［４］．
こうした２次元・３次元計算においては，ある程度の爆のように表され，ニュートリノ分布における角度変数＊１の発過程が見られるようになっているが，複数の爆発メカニ

１次モーメント &#$'で測られる前方集中度の値に反比例

ズム候補の是非について結論に至るところまでは達していしている．ここで，&'はニュートリノ分布による平均を表

ない．爆発エネルギーが十分でないことや，研究グループ

す．光度 "$!
%は，それぞれ，10$#，エネルギー !$"$，半径 &!

間において計算手法と結果が異なることも多く，いまだ議

erg/s，15 MeV，107 cm の値で規格化してある＊２．この時，論が残っている．なかでも主流となっているのは，多次元拡散近似では角度因子&#$'が過大に評価されており，加熱における流体力学的不安定性とニュートリノ加熱の組み合率を過小評価してしまうことが明らかになっている［７］．

わせである．衝撃波面が大きく歪んで非対称になり，一部

コアバウンス後の衝撃波付近（〜150 km）では，図５のよが大きな半径に達した際に，ニュートリノ加熱を起こす時

うに角度因子（#$の１，２次モーメント）は，拡散と自由間と領域が増えて，爆発へ達するだけのエネルギーを得

伝搬での極限値の間にあり，これらを精密に求めなければる．というのが大筋のシナリオである．
（図６左）ここでいけない．

問題となるのが，ニュートリノ輻射輸送である．多次元的52

空間２次元（軸対称）においては，ニュートリノ輻射輸

な流体分布におけるニュートリノ輻射を信頼度の高い方法

送計算は近似手法ながら，ミクロ物理も詳細に取り入れたで計算して，ニュートリノ加熱率を緻密に評価しなけれ

緻密な数値シミュレーションが行われている．拡散近似をば，近似への疑いを拭いさった最終的な答えを見いだすこ

用いた計算では，動径方向以外のニュートリノ輻射輸送のとはできない．つまり，３次元ニュートリノ輻射輸送を近記述において有利であるが，中間領域での角度分布推移に似なく解くことが課題である．

は問題点が残されている．他方，中間領域での記述におけ２．
５ ３次元ニュートリノ輻射輸送計算

る信頼度の高い方法としては，ray-by-ray 法が用いられている．これは，球対称輻射輸送計算コードを活用するもの

２次元・３次元になると，輻射輸送問題の難易度は一気で，中心から様々な空間角度方向について，独立に球対称に上がってしまう．ニュートリノ分布の記述には，空間３

ニュートリノ輻射輸送問題を解いて，流体計算と結合する次元とニュートリノ運動量３次元を合わせた６次元位相空ものである．この場合，動径方向以外の輻射輸送について間が必要である＊３．空間２，３次元の超新星計算は，位相は記述することができず，空間角度依存性が強すぎるとい空間５，６次元での扱いとなり，メモリ・計算量ともに巨う欠点が残っている．

大となってしまうため，ニュートリノ輻射輸送計算につい空間３次元におけるニュートリノの扱いは，簡単なモデては近似を行うのが現実的である．拡散近似ではニュート＊１ ニュートリノ伝搬方向を動径（&）
方向から測った角度 %$により決まる変数 #$"(,-%$
＊２ 核子１個あたりの加熱率100 MeV/s/nucleonは，換算すると1 gあたりの加熱率およそ1020 erg/g/sに対応する．
（図２を参照）＊３ 空間３次元およびニュートリノ運動量３次元（ニュートリノエネルギー１つ !$と伝搬方向を指定する角度２つ %$，"$）

６１４

Lecture Note

2. Computations of Radiation Transfer with Anisotropic Distributions: Examples of the Neutrino Transfer in Supernova Explosions

K. Sumiyoshi

図６ （左）
３次元超新星爆発におけるニュートリノ輻射輸送の状況を模式図で示した．流体不安定性により歪んだ衝撃波面に対して，引き続き外層から物質が降り積もっている．中心の原始中性子星から放出されたニュートリノが多次元の流体分布においてどのように輸送されるのかが，ニュートリノ加熱の大きさおよび爆発の可否を決める鍵である．（右）６次元ボルツマン方程式計算コードにより得られるニュートリノ密度分布の例．３次元超新星計算での中心コア分布を与えて，ニュートリノ輻射輸送計算を行った．

リノエネルギー１つのみ，ray-by-ray 近似では，エネ ル

$
$ #
$ &
!&
!& ! "$&"$"$"&$ !#$%$'
$%
#&
&&
&
#$
( &" &

ギー１つ・角度１つを方向ごとに扱って，計算を軽量化している．現在のところ，多次元ニュートリノ輻射輸送を近

"
"
!!#$
!!#$
&
,-"$ &&,-% & $
(
+%
+"$$
"#
$
!#
$(
$%
$%
&& &
%
"$
(
+%
(
+% &

似せずに解いた例は２次元における数例ほどしかなく，３次元では例がない．筆者らは，こうした状況を打破すべく，開発して，３次元ニュートリノ輻射輸送を解くことに成功

"
!!#$
(
+"$ &
$
$
$
$
!!
#!
"#
"
&&
" #!
( &,))(-(,+(+%

した［８］．

の 形 と な る［１０］．空 間３次 元 に よ る 微 分 項 に 加 え て，

２．
５．
１ ６次元ボルツマン方程式

ニュートリノの進行方向を指定する角度空間２次元による

空間３次元におけるボルツマン方程式を解く計算コードを

ここで "# 法による計算の概略を簡単に説明しておきた（３）

微分項が含まれている．ニュートリノの伝搬方向を記述す

い．角度依存性を積分したモーメント法とは異なり，"#る角度変数（#$#&,-%$，"$）については，グリットの配置

法では，輻射の角度分布もグリッド上で定義して輻射輸送に条件があり，ガウス積分で用いられる重み係数を用いて方程式を離散化して扱う方針を取る．これにより，角度分

グリットを定義する等の工夫が行われる．これは，ボルツ布について仮定することなく，拡散領域から自由伝搬領域

マン方程式を離散化する際に，複数の微分項による差分表までを扱うことができる．中間領域における任意の角度分

現の間で同時に満たすべき関係式が存在するためである

布を扱えることがメリットであるが，計算量が膨大になる［８，９］．また，定式化には座標系の指定による相対論の取りというデメリットもある．また，光源から遠いところでの扱いも問題となる．式（３）による保存形のボルツマン方程

自由伝搬においては，角度グリッドが有限であるために前式に従い，式を差分化して数値計算を行う．この時，計算方ピークの分布を記述するのが難しいという数値的な問題

負荷としては重くなるが，時間差分に関しては陰解法を採も あ る［９］．た だ し，超 新 星 爆 発 の 場 合 に お い て は，用する．ニュートリノ反応のエネルギー依存性により反応

ニュートリノ源となる中心天体の半径が 50 km 程度である

率が何桁も違っており方程式が堅いこと，陽解法ではのに対して，衝撃波およびニュートリノ加熱領域の位置はニュートリノが光速で伝搬することによる時間ステップの数 100 km 程度であるので，ある程度の角度グリット数を制限が厳しいこと，ニュートリノ分布が熱・化学平衡解へ設定することで対応できる．

到達するための安定性を保証すること，などが理由であニュートリノ輻射輸送では，フェルミ粒子を扱うこともる．実際に，球対称計算においては陰解法による取り扱い

あり，分布関数自体を解く手法が用いられている．ニューが用いられており，爆発に至るかどうかの長時間計算で役

トリノ分布関数に対するボルツマン方程式は，立っている．移流項の差分化については，セル境界での

$
$
$
$ &
$
$
!!
!&
" #!
"
#&
( &
' #!
( &,))(-(,+

ニュートリノ分布値を用いて表している．この際，中央差分と風上差分を採用しており，拡散・自由伝搬の間で重み（２）をつけて切り換える方式を取っている．
衝突項には，ニュートリノ反応による生成・消滅過程である．ここでニュートリノ分布関数は６次元位相空間での関数である．これを３次元球座標を用いて表して，さら

$
$
!!
%%
$
%%
!!$
#!!$%-$$!
"$"!'*($% （４）-$#!
( '*(-!$%-に保存形に直すと

６１５

Journal of Plasma and Fusion Research Vol.88, No.10 October 2012

選択も重要な点である．現在の計算ではpoint Jacobi前処理

および，散乱過程と Bi-CGSTAB 法を組み合わせて用いているが，反復法が

(#+)
)
%
!#
)
+
!
+)&
)
!!)
!
" $!#
%%
+
)
3%
%&
#%2
'#
, 2+)3&&"%

収束しにくい場合にも対応できる新たな解法［１２］を開発して採用することも行っている．

(#+)
)
（５）
+)!+&
)
!!)
"#
%&
+
)
3%
%%
#%2
&&
"
%

２．
５．
２ ３次元超新星コアでのニュートリノ輻射輸送計算を組み入れる［８］．それぞれ陽子・中性子・原子核・電子

実際の超新星コアにおけるニュートリノ輻射を解く問題へ完成した６次元ボルツマン方程式を解く計算コードは，等をターゲッ ト と す る 過 程 が あ り，反 応 率 を 代 表的にの適用が行われて，これまでの近似的な方法では扱えな，%,.-2%，%2+)3%（以下同様）などと表した．
+&
+&
+
!
+)&
%)*2%

かった新たな様相が明らかにな り つ つ あ る．図６右 に

これらはニュートリノのエネルギー・角度に依存しておは，３次元超新星計算での流体分布をもとに６次元ボルツり，詳細釣り合いを加味して，反応過程ごとに具体的な表

マン方程式計算コードによりニュートリノ分布の時間発展

式を組み込む必要がある．また，輻射における光子（ボーを追った輻射輸送計算でのニュートリノ密度分布の例を示

ズ粒子）の扱いと異なり，ニュートリノはフェルミ粒子でした．２次元・３次元に変形した超新星コアにおいて，ray

あり，分布関数 )%は０から１の間に制限されている．フェ-by-ray 近似によるニュートリノ輻射輸送計算では動径方ルミ統計に従って，ニュートリノが生成される際に，行き向のニュートリノ伝搬のみ記述していたが，６次元ボルツ

先の量子状態が満たされている場合には抑止される効果，

マン方程式計算コードによる結果では，動径方向だけでな

ブロッキングファクターが必要である．)%が１の時は反応

く，非動径方向（３次元球座標での '，$の変わる方向）への形で衝突項に組み込
レートがゼロとなるように，%!!)%&のニュートリノ伝搬（流束）を記述することができる．まむ．またニュートリノ・反ニュートリノで対生成・消滅す

た，拡散領域の外でも多次元的なニュートリノ輻射輸送のる過程記述が可能である．実際の超新星ダイナミクスへ適用するには，流体計算とニュートリノ輻射輸送計算を結合して，(#+) 0)-1
)
%
!#
+
!
+)&
)
%
!
" $!#
%)
%
#%)
//'#
, 0)-1&&"%

時間発展を解かなければならない．現在，多次元流体計算
コードと６次元ボルツマン方程式計算コードを組み合わせ(#+) 0)-1
+
!
+)&
!!)
!!)
"#
%
%
%
%&
%&（６）
#%,.2
&&
"
%

た，ニュートリノ輻射流体計算コードの開発テストを行っている．また並列化により大規模化へ向けた対応も進んで
いる．

もある．)
%は反ニュートリノの分布関数である．この場合には，２つのニュートリノ分布が積として現れるため，方２．
６ 大規模計算シミュレーション

程式が非線形となるが，これらの過程は高温高密度の平衡&( 法による３次元ニュートリノ輻射輸送計算が可能に状態に近い領域でのみ起こるため，片方の分布を固定して線形化することができる．

なってきた．しかし，現状の計算機資源では不十分なため，反応率の計算では，高温高密度物質での組成（陽子・中いくつかの制約をおいて進まざるを得ない．ここで３次元

性子・原子核・電子・陽電子・光子）の詳細が必要とな超新星爆発でニュートリノ輻射輸送を解く際に必要な計算

る．高温高密度や中性子過剰な領域では，核子・原子核の機資源について述べて，将来のスーパーコンピュータに向性質も変化しており，ニュートリノ反応率自体へも影響がけた展望を記しておきたい．

及ぶ．

６次元ボルツマン方程式において，ある時間ステップでここまでの式を陰解法の形のニュートリノ分布を保存するには，６次元配列が必要となる．空間３次元のグリット数を$20)+,，ニュートリノ空間"*"!!"*
"*"!(
$"'
",

３次元のグリット数を $%とすると，$20)+,$%がニュートリ（７）ノ分布に必要な配列の並び数である．現状の計算資源での*"!

$!の形とな実 行 規 模 と し て，$20)+,$"%&##"#&$，$%$'#!"#!$

る．ここで "は，６次元空間グリット上でのニュートリノ

を想定すると，ニュートリノ３種の分布情報だけで 20 GB

分布関数 )%を順番に並べたベクトルである．計算の手順と程度の規模となる．

で表現して，整理すると線形方程式 !"しては，時間ステップ *におけるニュートリノ分布のもと

ボルツマン方程式を陰解法で解くために，大規模疎行列

で，移流・衝突項による行列要素を計算して，大規模疎行による線形方程式の解法が主な計算負荷であり，行列要素

列による線形方程式を解いて，次の時間ステップ*"!におを納めるためのメモリが大規模となる．大規模疎行列の対

けるニュートリノ分布を求める．

角に並ぶブロック行列のサイズ # は，ニュートリノ反応の行列 !は特定のパターンを持った大規模疎行列 で あ

寄与でサイズが異なり，エネルギーが変わる散乱過程を含り，時間発展の対角要素，移流項による非対角要素，衝突

めない計算では，# $$'%$$% であり，必要なメモリは項によるニュートリノ状態変化を表す対角ブロック行列か# "$(%$20)+,，演算量は # #$(%$20)+,でスケールする＊４．上らなる．この疎行列による方程式を解くためには，反復法

記のメッシュ数の場合，行列要素の格納に２TB 程度が必による数値解法［１１］を用いており，前処理と反復法の手法要であり，浮動小数点演算数は，ニュートリノ分布の時間＊４ 反復法の演算量をブロック行列の前処理（反転を想定）による演算回数で見積もった．

６１６

Lecture Note

2. Computations of Radiation Transfer with Anisotropic Distributions: Examples of the Neutrino Transfer in Supernova Explosions

１ステップ更新で 6×10１２程度となる．すべての反応（およ謝K. Sumiyoshi

辞

本解説は，山田章一，長倉洋樹，固武慶，滝脇知也，松び相対論）を考慮する場合には，ブロック行列がさらに大古栄夫，今倉暁，櫻井鉄也氏らとの多次元超新星爆発の研きくなり，メモリと演算量が巨大となる［２］．
超新星爆発の計算では，鉄コアの重力崩壊を初期条件と究プロジェクト［新学術領域研究（２０１０５００４，
２０１０５
００５），科

して，コアバウンスを経て，衝撃波が停滞・復活するまで研費（２２５４０２９６，２４２４４０３６），HPCI 戦略プログラム分野５］の長時間発展（約１秒）を追う必要がある．一方，ニューに基づいたものである．３次元ニュートリノ輻射輸送計算

トリノ輻射輸送計算における時間ステップは，陰解法ではは KEK，京都大学 YITP，東京大学，大阪大学 RCNP にお−5（空間解像度に依るが典型的に）105s 程度であり，10 時

けるスーパーコンピュータを用いて行われた．

間ステップの計算が必要である．系統的な計算を行うことを考えると現在の計算資源では，次元を落とすか，グリッ

参考文献

ド数を抑えて，長時間発展計算を行うのが精一杯である．

［１］H.-T. Janka, K. Langanke, A. Marek, G. Martínez-Pinedo
and B. Muller, Phys. Rep. 442, 38 (2007).
［２］K. Kotake, K. Sumiyoshi, S. Yamada, T. Takiwaki, T.
Kuroda, Y. Suwa and H. Nagakura, Prog. Theor. Exp.
Phys. 2012, 01A301 (2012).
［３］N. Ohnishi, K. Kotake and S. Yamada, Astrophys. J. 641,
1018 (2006).
［４］T. Takiwaki, K. Kotake and Y. Suwa, Astrophys. J. 749,
98 (2012).
［５］K. Sumiyoshi, S. Yamada, H. Suzuki, H. Shen, S. Chiba andH. Toki, Astrophys. J. 629, 922 (2005).
［６］K. Nakazato, K. Sumiyoshi, H. Suzuki and S. Yamada,
Phys. Rev. D 81, 083009 (2010).
［７］S. Yamada, H.-T. Janka and H. Suzuki, Astron. Astrophys.
344, 533 (1999).
［８］K. Sumiyoshi and S. Yamada, Astrophys. J. Suppl. 199, 17(2012).
［９］J.I. Castor, Radiation Hydrodynamics (Cambridge University Press, Cambridge, 2004).
［１０］G.C. Pomraning, The Equations of Radiation Hydrodynamics(Pergamon Press, Oxford, 1973).
［１１］R. Barrett et al., 長谷川里美，長谷 川 秀 彦，藤 野 清 次（訳）：反復法 Templates（朝倉書店，１
９
９
６）
．
［１
２］A. Imakura, T. Sakurai, K. Sumiyoshi and H. Matsufuru,
JSIAM (2012) in press.

京コンピュータでは，空間２次元での解像度が高い系統計
算が可 能 と な る が，空 間３次 元 で の 長 時間 発 展 計 算 はExaflops スケール以上のスーパーコンピュータが必要である．さらに，一般相対論での検証，重元素合成過程の探求が控えている．計算規模を順次スケールアップしていき，究極の爆発計算により超新星現象の全貌を明らかにしていくことが将来にわたる課題である．

２．
７ まとめ
星の進化の最期におきる超新星爆発という華々しい現象
では，ニュートリノ輻射輸送が重要な役割を果たしている．超新星コアでは，拡散から自由伝搬まで幅広い環境が存在しており，輻射輸送を一貫して記述する必要がある．
途中で停滞してしまう衝撃波を後押しするニュートリノ加
熱過程が働くのは，ニュートリノ分布が非等方な中間領域である．爆発か否かを探索するには，緻密なニュートリノ
輻射輸送計算が不可欠である．
本章では，超新星爆発におけるニュートリノ輻射輸送の役割と最近の発展について解説した．球対称計算では，第一原理的なニュートリノ輻射流体計算が可能となっており，系統的な研究が行われている．球対称計算では爆発を再現することができないことが判明しており，多次元での数値シミュレーションが盛んに行われている．この時，ニュートリノ輻射輸送については主に近似手法が用いられており，爆発メカニズムの新たな様相が明らかになってき
た．最後に，ニュートリノ輻射輸送の近似をせずに扱う方向として，!!法に基づいて空間３次元でのボルツマン方程式を解く計算手法について紹介した．６次元位相空間での時間発展を追うために，計算資源・時間ともに膨大なもの

こうした研究では，計算機・並列アルゴリズム等の研究者との連携が欠かせない．シミュレーション計算の規模が大きくなるにつれて，一つの分野だけではなく，総合的に技術を蓄積して応用していくことが不可欠である．超新星
爆発は宇宙物理の話ではあるが，工学分野との情報交換や協力を通じて，極限状態での物質・反応・流体の複合３次元問題へのチャレンジを，広く他分野での問題へも活かしていきたいと考えている．

!!!!!!!!!!!!!!!!!

なった．

すみ

よし

こう

すけ

住 吉 光 介東京都立大学大学院出身，理化学研究所等を経て，現在は沼津工業高等専門学校教養科教授．専門は原子核と宇宙にまたがる物
理で，超新星爆発が主な研究テーマ．スパ

コンに関わる仕事も多く，原子核データ整備・天文数値シ
ミュレーション・計算科学と広く共同研究を展開中．東京都
郊外の育ちだが，両親の影響で江戸っ子の職人気質を受け継ぐ．長年の Macintosh ユーザ．こだわりの趣味を我慢しつつ，端末に向かう日々．

!!!!!!!!!!!!!!!!!!!!!!!!!!

り，ボルツマン方程式を直接解く計算手法が現実のものと!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!

となるが，近年のスーパーコンピュータの急速な進展によ

６１７

J. Plasma Fusion Res. Vol.88, No.11 (2012)６
３９‐６４８講座 輻射流体シミュレーション

３．輻射輸送における原子過程
西川亘

岡山大学大学院自然科学研究科（原稿受付：２
０
１
２年７月２
５日）

プラズマ中での輻射輸送を計算する上で必要な原子過程の解析と具体的な計算方法を紹介する．特に実験室で生成されるプラズマの場合，プラズマのサイズが小さく，温度・密度が急激に変化するので，個々の場所で原子過程の熱平衡を仮定する局所熱平衡近似が成り立たない場合も多い．また扱うイオンの束縛電子数が多い場合には個々の状態を個別に取り扱うことが難しくなり，様々な近似計算が必要となる．ここでは，筆者がこれまでに携わったモデルを中心に輻射輸送を解く上で必要な原子過程のモデルの詳細と扱う際の問題点等を解説する．
Keywords:
atomic process, emissivity, opacity, LTE, TE, CRE, UTA

３．
１ プラズマ中原子過程と輻射輸送吸収係数を計算することを前提においているが，輻射場，

プラズマ中で輻射によるエネルギー輸送が無視できない

特にラインスペクトル状の強い輻射場がある場合について

場合，流体方程式と同時に輻射輸送方程式を解く必要があ

簡便に輻射励起の寄与を取り扱う方法を紹介する．この章

る．一般に原子番号の小さい原子からなるプラズマの場合を読めばひととおりの放射係数，吸収係数の計算方法と問は，比較的，流体の運動エネルギーが十分大きいため輻射

題点が理解できると期待している．

によるエネルギー輸送の寄与は小さいが，原子番号が大き「輻射輸送における原子過程」というタイトルからすれ

くなるにつれて輻射によるエネルギー輸送の寄与が大きくば輻射輸送方程式の放射係数，吸収係数を構成する物理かなり，流体方程式と輻射輸送方程式を同時に解くことが必ら説き起こすべきなのかもしれないが，原子過程研究の究須となる．輻射輸送方程式を解く上で必要な係数は，輻射

極は原子，イオンのまわりの電子状態など，量子力学で厳の放射係数（Emissivity）と吸収係数（Opacity）である．こ格に記述される素過程をいかに統計的に処理するか，といの章では，放射係数と吸収係数を計算する上で必要なプラ

う点に帰着できる．本章では，輻射輸送を解く場合に必要

ズマ中原子過程の主要な物理とモデル，問題点などを基礎なプラズマ中原子過程モデルがどのように作られ，計算さ方程式とともに解説する．また，筆者が経験したいくつかれているのかを示すため，あえて，サハ・ボルツマンの関の特に注意すべき点については個別に説明する．

係の導出より解説を始める．具体的なイメージが沸くよう

プラズマ中原子過程の基本的なモデルの説明と問題点をに，筆者が実際に計算した経験のある水素プラズマの電離

解説した後，プラズマ中原子過程の代表的な近似モデルで度の計算を具体例として示しながら説明を進める．

ある熱平衡，衝突輻射平衡を仮定して得られたポピュレー

なお，輻射流体シミュレーションと原子過程の関わり

ションの特徴を説明する．その上で，放射係数と吸収係数は，プラズマの輻射の放射，吸収係数を介して輻射輸送方の計算を構成する，自由自由，自由束縛（束縛自由），束縛

程式と，状態方程式を介して流体方程式と関係する．この束縛を説明し，実際にポピュレーションからスペクトルを章では輻射輸送と原子過程の関係だけを述べ，状態方程計算する際，原子番号が小さく線スペクトルが主となる場

式，特に流体方程式を解く上で必要になる，流体（プラズ合に重要となるラインプロファイル，プラズママイクロ

マ）の比熱と圧力については触れない．

フィールド，原子番号の大きなイオンの部分電離状態の場
UTA（unresolved transition array），密度が固体密度程

３．
２ 水素プラズマから始めるプラズマ中原子過程モデル度以上に大きくなった場合など，個別の場合について注意

プラズマ中の自由電子とイオン・原子に捕まっている束合など，個々の線スペクトルが重なり合い連続的に見えるすべき点を挙げる．複雑なイオンの場合，すべての状態を縛電子との間に熱平衡が成り立っている時，&"!回電離し考慮したレート方程式を解くことが難しいため，平均イオ

たイオンの数密度 #&"!とエネルギー !を持った自由電子

ンモデルに代表される簡略化を行いスペクトルの計算を行の数密度 #'#!$の積と，それらが結合した &回電離したイ

う場合もあるが，この際の注意点についても説明する．基オンの数密度 #& との比は解離前後の束縛電子の状態数の本的にこの章では光学的に薄いプラズマ近似で放射係数と，およびエ比 %&"!!
%&，ボルツマン因子 #%$#
!
("$#$
!""
!#
$$

3. Atomic Process Modeling and Radiation Transport
NISHIKAWA Takeshi

author’s e-mail: nishikawa.takeshi@okayama-u.ac.jp

６３９!2012 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.88, No.11 November 2012

ネルギー "を持った自由空間の電子の状態数 #&
*$'
$.#!
).
%

いたものが図１である．図１より考慮する束縛状態の最大の積

主量子数が変わると電離度が 0.5 になる温度・密度の場所$.#%%(""%*&"' '(""
""$
+
*-+!!
#
"# $ ).
%(
'(
+'&*
*

が大きく移動することがわかる．逆にいえば，考慮する束（１）

縛状態の最大主量子数をうまく選べば電離度が 0.5 になる
温度・密度を特定の場所に移動させることができる．しかで表せる．ここで $+は (回電離したイオンから (""回電しながら実際の局所熱平衡下の水素プラズマの電離度は一離したイオンになるための電離ポテンシャル，+! と &*は意的に決まっているはずである．

ボルツマン定数と電子温度である．自由電子状態数の式中束縛状態をどこまでとるのかについては古くから議論がの *はプランク定数である．,*を電子の質量とした時，自あり，大きく分けると物理的描像と化学的描像がある．荷由電子のエネルギー "とその運動量 .との関係は，自由空電粒子間の相互作用に立ち入る物理的描像に対し，化学的間 で の 関 係 式 "#.#!
,*を 用 い て，.を "に 変 換 し，#描像は相互作用の詳細は解かず，反応前後の状態数とエネ""!で積分するとルギーで解く．先に説明したサハの式も化学的描像に基づ
く．水素プラズマの精密な分光結果との比較をみる限り，%(""%* '(""$,$*-+!! + "&%#$*
#
%(
'(
+'&*
*

プラズマ中での相互作用から束縛状態の占有確率を計算
し，化学的描像を用いて占有密度（ポピュレーション）を%$(" *-+!! " ")"+'&*!

計算し，スペクトルを計算した方が実験で得られた分光ス
ペクトルをより正確に再現する［４］．

$!
#$,*+'&*'#$（２）# (""*-+!! + "#!
" ．
'(
+'&*
*#

ここでは著者が考案した化学的描像に基づき，考慮する
最大主量子数を与える方法［５］を簡単に紹介し，どこまでイオン化状態が一つだけ異なった基底状態のイオン数密度

束縛状態を考慮すべきかについて基礎的な指針を与える．

の比を与えるこの式は一般にサハの式と呼ばれている．イ

まず，プラズマ中の原子・イオンは統計的に分布するとマオン化状態が同じ時，基底状態と励起状態の関係は束縛状

イクロフィールド中に存在すると仮定する．水素プラズマ態の状態数の比とボルツマン因子だけになるので，平衡状を考えるとプラズマ中の原子のまわりのポテンシャルは，態におけるプラズマ中の励起状態を含む原子・イオンの数

原子核のポテンシャル (()!
$#%!/と均一な外場 #の重ね合

密度を与えるこれらの式はまとめてサハ・ボルツマンの関わせとなる．ここで ((，)，および #!はそれぞれ原子の核係と呼ばれている．ここで，電離度 ((,はすべての原子・イ

電荷，電気素量，真空の誘電率である．水素原子においてオンの数密度の和%!## %(とすると，((,#%*!
%!で定義

主 量 子 数 -の 状 態 の 電 子 の エ ネ ル ギ ー は

される．

# # #
&#
!,*)%((#!
! * - で与えられ，外場中に存在しても変化

さて，このサハ・ボルツマンの関係を使えばプラズマのしないとする．この場合，原子のまわりのポテンシャル分

電離度を計算できるであろうか？実際に，例えば水素プラ

布はサドルポイントを持ち，また，このサドルポイントのズマの電離度を計算しようとするとこれだけでは不十分で

高さは外場 #の関数となる．エネルギーがこのサドルポイ

ある．束縛状態をどこまで取り入れるべきなのか，というントより大きい電子を自由電子，小さい電子を束縛電子とことを決めないと計算できない．レーザー生成プラズマのシミュレーションでは主量子数 -が１０の状態まで取り入れることが多かった．一方，放電プラズマや磁場閉じ込め核融合のプラズマでは -#$!$"!!が使用されてきた．おそらくレーザー生成プラズマは初期状態が固体から生成される一方，放電，磁場閉じ込め核融合のプラズマは気体から生成されるということが理由の一つと考えられる．従来の多くのレーザー生成プラズマのシミュレーションが遮蔽水
素モデル（Screened Hydrogenic Model）をベースとしたものが多く，エネルギー準位を計算するために Mayer より導入され［１］，More によって改良された遮蔽水素モデルの計算に必要な遮蔽定数［２］が -#"!までしか提供されな
かったという点も理由の一つと考えられる．従来，それほど問題視されてこなかった取り入れる束縛状態数の問題であるが，実は状態方程式への応用に伴って古くから問題と認識されていた［３］．実際問題，電離度の等高線がそのプ
ラズマの温度・密度のどこにくるかは考慮する束縛状態で大きく変化する．この一つの例として，水素プラズマの電図１

離度が 0.5 になる等高線がそのプラズマの温度・密度のどこに来るか，考慮する束縛状態の最大主量子数依存性を描
６４０

水素プラズマの電離度が 0.5 の等高線．考慮する状態の最大主量子数によって温度，密度が変化する．

Lecture Note

3. Atomic Process Modeling and Radiation Transport

すると，主量子数 .の状態が存在する閾値電場がT. Nishikawa

上の準位でも比較的長い寿命を持つ準位もある上，多励起状態の状態数は各励起電子の取り得る状態数の積で状態数

# & $
*(
!.)# %--$ % +%
'%#
! +.

が増えるので，ボルツマン因子で減る分と状態数の増加の（３）関係を考慮する必要があり一筋縄ではいかない．特に原子と計算でき，外場 !の存在確率を %&!'とすると主量子数

番号が大きくて束縛電子が多いイオンの場合は内部状態をが .の状態が存在する確率は多く持ち，全エネルギーはイオン化ポテンシャルよりもは

!.)

/. ##

!

るかに上の状態であるにも関わらず，束縛電子分布という!',!%&観点からみたイオンの大きさという点では十分小さく，例（４）

えばイオン球半径内に収まる電子分布を持つイオンに制限と表すことができる．例えば，電場分布として，イオン間しても束縛電子が多くなるにつれて考慮すべき束縛状態がのクーロン相関を無視したホルツマーク場

どんどん増える．現実問題としても，4f が開殻のイオンのすべての電子状態を解くことのできる原子コードが存在し

%

#
"&
"'
!'
"0'
1
.
/&
,0
#!!'&
# "# 0-20&!0$!#'（５）% !

ないという問題もある．計算精度としては生真面目なモデ
ルから一歩下がるものの，それらを解析するために平均イを用いれば，式（４）の積分が実行でき，主量子数が .の状オンモデルやそれをさらに一般化した STA（Super Transi-態が存在す る 確 率 を 具 体 的 に 計 算 で き る．こ こ で "はtion Array）モデルなどが開発された．

"#!!
!!で定義され，!!はすべてのイオン化状態を分離して詳細な電子配置の状態を考慮して原子過程を解くモデルをDetailed Configuration

#!
$

(0*
$
%
!!##
%! 0"
%#
%
"&
!

（６）

Accounting（DCA）モデルと呼ぶ．多重項の違いによる状で表されるホルツマーク基準電場である．ここで $0と (0（DTA）モデルと呼ばれる．一方，各イオン化状態を別々はそれぞれ，外場を作り出すまわりのイオンの数密度と核に解かず，仮想的な平均イオンについてのみ電子状態を解電荷である．おおまかな束縛状態数を計算したい場合，

くモデルを平均イオンモデルという．プラズマ中イオンの態をすべて考慮した場合は，Detailed Term Accounting#"#!での近傍で "&"'"!
%であることを用いると，式$$%

電子状態は局所場近似で自己無撞 着 場 を 解 く Hartree-（４）の積分がFock-Slater 法などで解析する場合もあるが，多くは計算の

!.)

/. ##

!

&#()
!'
,!%&
# "( % "#+ $ ## %. (0 $0

簡便さと非平衡プラズマの応用，状態方程式など熱力学量の計算の容易さから，電子間の相互作用を遮蔽定数として（７）取り入れた遮蔽水素モデル（Screened Hydrogenic Model）と解析的に実行できる．

がよく用いられてきた．古くは吸収係数におけるラインス

ここでは束縛状態しか紹介しないが，自由状態の取り扱

ペクトルの重要性を鉄イオンの解析から示した H. Mayer

いに興味のある方は拙著
［６］等をご参考いただきたい．束［１］，非平衡原子過程の取り扱いを確立した Lokke［７］，縛状態数は電離度の等高線を大きく変えるが自由状態につ

今日でもなお参照される磁場閉じ込めプラズマでのクーリ

いては自由空間の電子の状態数で計算しても，プラズマ中

ングレートを計算した Post の計算［８］など広く使用されてのポテンシャルを考慮してもそれほど大きな変化はない．

きた．以上の結果は状態数として主量子数のみを考慮した

なお，ここで計算した束縛状態の消失やその結果として自平均イオンモデルであるが，方位量子数を含んだ状態を考由状態の変化を考慮することで連続準位の低下（Contin-

慮したモデルなどへの応用も行われ，金など原子番号が非

uum lowering）や圧力電離（Pressure ionization）を見積も常に大きいイオンからなる部分電離プラズマからの X 線ス

ることもできる．

ペクトルの解析などで実績を上げてきた［９］．
非平衡プラズマへの応用は次の節にゆずるとして，平均

３．
３ プラズマ中の原子過程モデルの簡素化と平均イオンモデル

イオンモデルを用いて熱平衡状態のポピュレーションを計

水素プラズマの場合など，比較的束縛電子数が少ない場

フェルミモデルで自由電子は連続体近似で計算し，束縛状合は３．
２で説明したような考慮する束縛状態をどこまで取態は遮蔽水素モデル等で計算した離散的なエネルギー準位算する方法を説明する．基本的には有限温度のトーマス・

ればよいか？という問題だけであるが，束縛電子数が多いで計算する．化学ポテンシャルを $とすると，エネルギー

複雑な電子状態が可能なイオンの場合はもう少し丁寧な取が０より大きい電子を自由電子とすれば自由電子数はり扱いが要求される．すなわち，束縛電子数が多くなると，$!
#$--,*&#
$-##!
" #
"
"!
+#
*&
# ,

最外殻の一電子励起状態だけではなく，二電子励起状態なども含む必要があり，考慮すべき束縛状態数は爆発的に増（８）で与えられる．ここでえる．どこまで考慮すべきかは悩ましい問題で，単にイオン化ポテンシャルよりも高い準位は寿命が短く，すぐに自

1"!#
,1
! "
1!0'
"-20&
%

#
0'
##
"&

動電離すると考えることができれば，それ以下の束縛状態#を考慮するだけではあるが，イオン化ポテンシャルよりも
６４１

（９）

Journal of Plasma and Fusion Research Vol.88, No.11 November 2012

で定義される．同じ化学ポテンシャルを用いて状態 &の占有確率（ポピュレーション）"& をフェルミ・ディラック分布
"& #

%
&
! !!
!"$'%! &
"
&"#$

（１０）

を用いて計算する．ここで %& と !& は状態 &の縮重度（状態数）とエネルギー準位である．"& は非整数の各準位の束縛電子数となるが電離度は$#&#$!# "&（１１）で計算できる．ここで $はプラズマを構成するイオンの核電荷である．ここで状態を示す添え字として一般的な &を使用したが，そのモデルで考慮している状態が主量子数のみの場合は (，方位量子数も含む場合は (，'で置き換えればよい．熱平衡状態の金プラズマの平均電離度の計算結果を図２（a）に示す．熱平衡の場合，(#!の状態のイオン化ポテンシャルが 100 keV 程度ある金でも電離した方が状
態数的にエネルギーを得するので 10 keV 程度の電子温度で完全電離する．図より同じ温度であれば密度が高くなる
ほど電離度が小さくなっている．これは密度が大きくなると式（８）より１原子あたりの自由電子の状態数が小さくな
り，結果として電離度が小さくなる．温度が低くて固体密度を超える辺りからは，こんどは束縛状態が減る影響で，
いわゆる圧力電離が起こり，逆に電離度が増える．
なお，ここで説明した平均イオンモデルで計算した熱平衡のポピュレーションと３．
２で説明したサハ・ボルツマンの関係を用いた計算の結果は一致しない．電子数が少ない
場合，特に水素イオンの場合は大きくずれることが知られている．計算機資源の発達した現在では，平均イオンモデ
図２

ルで放射係数や吸収係数を計算する必要性は少ないと思わ

熱平衡と衝突輻射平衡での金プラズマの電離度の比較．

れるが，固体密度付近からそれ以上の密度のデータが必要な場合，トーマス・フェルミモデルやそれを基づく平均イ

ラズマ中原子過程を解いた結果としての電離度がどうなる

オンモデルの結果を援用する場合があるかもしれない．幅

かがきれいに説明できるので，平均イオンモデルを用いた広い温度・密度のデータが必要な場合は別々のモデルで計古い結果ではあるが，金プラズマの平均電離度を例に取り算したデータを使わざるを得ない場合があるが，接続する

説明する．
プラズマ中の原子，イオンの各電子状態は大きく分けて

温度・密度で電離度，および輻射の放射，吸収係数が連続

電子もしくは輻射を伴い状態を変化させる．具体的にいう的に変化するかどうか確認が必要である．
平均イオンモデルでは水素様イオンのエネルギー準位をと，電子衝突電離とその逆過程の三体衝突再結合，電子衝

記述する量子数，例えば (，'で状態を決めている．一方，突励起とその逆過程の電子衝突脱励起，輻射電離とその逆すべてのイオン化状態の詳細なスペクトルのうち，同じ程過程の輻射再結合，および輻射励起とその逆過程の輻射脱

度のエネルギーのラインスペクトルを出す状態をスーパー励起である．三体衝突再結合が相対的に大きくなりつつあ

シェルとしてまとめて扱い，なるべく少ない計算で DCA/

る密度では二電子性再結合も重要になる．
電子衝突励起，電子衝突電離過程については入射電子のDTA モデルの結果を再現しようとするのが STA であるエネルギーの関数として与えられる衝突断面積を熱平衡電［１０］．

子の速度分布関数，マクスウェル分布で積分することによ３．
４ プラズマ中原子素過程と非平衡近似

り速度係数が得られるが，イオン化の場合は Lotz の式［１１］３．
３までで説明したモデルは熱平衡の場合で，量子統計など，励起の場合は振動子強度で表される簡便な公式を使の問題であるが，実際のプラズマ中では，原子過程の平衡

用することも多い．入射する電子によって起こるので当然状態を仮定できない場合もある．ここではプラズマ中原子のことながら，衝突電離，励起，脱励起のレートは電子密素過程とその速度（レート）係数の説明，典型的な平衡，準度に比例する一方，三体衝突再結合のレートは電離する電平衡モデルの結果と素過程との関係などの説明を行う．プ子にもよるので電子密度の２乗に比例する．またイオン化
６４２

Lecture Note

3. Atomic Process Modeling and Radiation Transport

T. Nishikawa

エネルギーや励起エネルギーのボルツマン因子を含むのでが大きい場合は衝突過程でポピュレーションは決まるのでイオン化エネルギーや励起エネルギーが小さいほど大きく

電離度の等高線は熱平衡モデルの場合と同じになる．密度なる．輻射電離も衝突過程と同様に衝突断面積を熱平衡電が小さくなると相対的に輻射脱励起の寄与が現れる．衝突子のマクスウェル分布で積分することにより速度係数が得

輻射平衡の場合，遷移エネルギーの２乗に比例する輻射脱られるが，Kramers の公式［１２］から計算される式を使うこ励起の影響が大きく電子温度が 10 keV 程度であってもとも多い．輻射電離は入射する輻射と電離する電子による!!!の状態まで電子はつまっている．図中の温度が高くので電子密度に比例する．輻射励起，脱励起は準位の寿命て密度が低い領域では密度が変わっても電離度はほとんどを表す A 係数で決まり，振動子強度に比例し，遷移エネル

変化していない．これはコロナ平衡で，電離度は電子衝突

ギーの２乗に比例する．

電離と輻射再結合で決まり，どちらも電子密度の１乗に比衝突過程だけで電子状態が決まっている場合，原子過程

例するので密度依存性がなくなる．密度が高くなるにつれは状態数と状態間のボルツマン因子だけできまり，前述のて電子密度の２乗に比例する三体衝突再結合の影響が出てサハ・ボルツマン分布となる．輻射輸送を考えるプラズマ

来て密度依存性が現れる．二電子性再結合を考慮したモデの温度・密度の変化する空間スケールが十分大きい場合，

ルでは，三体衝突再結合の影響が出る前に影響が出る．二原子過程と輻射場は個々の場所で熱平衡を仮定できる局所

電子性再結合の素過程は，束縛電子を電離させるほどに十熱平衡近似が使用できる．一方，原子過程の非平衡を考え分大きなエネルギーを持たない電子がイオンと衝突，束縛なければならないのは局所的な輻射の強度が強くて電子衝

電子を励起し，自身もエネルギーを失い再結合することに

突電離，励起よりも輻射電離，励起の方が大きい時と，局

より二電子励起状態ができ，その後，電子密度が比較的小所的な輻射の強度が小さくて輻射電離，輻射励起が無視でさい場合は別の電子が衝突する前にどちらかの励起電子がきる場合となる．実験室で局所的な輻射場が大きい場合は

輻射脱励起することにより再結合が完了する．電子密度が光励起の X 線レーザーなどの研究などに限られる．宇宙に大きくなると輻射脱励起が起こる前に次の電子が衝突し励目を向けると中性子星やブラックホールなどのコンパクト

起電子のエネルギーを奪い，一電子励起状態になる三体衝星の周りで観測される光電離プラズマなどが具体例として

突再結合になってしまう．通常の二電子性再結合の速度係挙げられる．一方，磁場閉じ込め核融合のプラズマ，レー

数は輻射脱励起の速度係数と二電子励起状態を生成する衝

ザーや放電で生成されるプラズマなど，実験室で得られる突過程の速度係数の比で与えられ，多くの計算では電子密

おおよそすべてのプラズマは，そのプラズマサイズが小さ度が増えると三体衝突再結合速度係数に上乗せされる形で

い，もしくは大きくても密度が小さいためにプラズマ中で導入されている．一方，三体衝突再結合の速度係数は電子生成された輻射はすぐにプラズマ外へ出てしまい，局所的衝突電離の速度係数との詳細釣り合いから計算されているな輻射場が小さいと近似することができる．この様な場

ので，単に二電子性再結合の速度係数を単に足すだけで原合，良く用いられるのが衝突輻射モデル（Collisional Radia-

子過程への寄与を取り込んだモデルでは，輻射再結合から

tive Model）である．最近は，後で述 べ る 衝 突 輻 射平衡

三体衝突再結合が効き始める領域では再結合レートが過大（Collisional Radiative Equilibrium）を衝突輻射モデルと決評価になっている可能性がある．

め打ちしている事例が多いが，少なくともオリジナルの論

３．
５ 放射係数と吸収係数

文［１３］を読む限りは平衡状態とは限らないし，彼らが衝突輻射と呼んだのは，励起状態の占有密度は基底状態から直熱平衡の場合はサハ・ボルツマンの関係，非平衡の場合

接励起，直接再結合ですべてできたわけではなくて，輻射

はレート方程式を解き，各励起状態，各イオン化状態の占

脱励起や衝突脱励起によるカスケードの影響が大きいこと

有密度が求まれば，輻射輸送方程式で用いる放射係数と吸をみいだし，励起状態の占有密度の生成過程として衝突輻収係数を計算することができる．束縛準位間の遷移による射モデルと命名したようである．よく使われるモデルとしラインスペクトルは速度方程式で用いた係数に３．
６で後述ては電離先の占有密度を０として解く純粋な電離プラズ

するラインプロファイルを考慮することにより計算でき

マ，基底状態の占有密度を０として解く純粋な再結合プラる．束縛自由，自由束縛のスペクトルはレート係数を計算

ズマで，特に磁場閉じ込め核融合のプラズマなど比較的密するのに用いた輻射の衝突による電離断面積から計算でき度が薄いプラズマの分光でのラインスペクトルの強度比なる．

どの解析によく用いられている．一方，流体コードの中で速度方程式の計算で出てこない放射係数と吸収係数の成の輻射輸送の計算で使用するモデルとして良く使用される

分は，自由電子同士の衝突過程による輻射の放射，吸収でのが衝突輻射平衡である．このモデルは，プラズマ中で生あり，制動放射，逆制動放射［１４］として計算される．なお，成した輻射はすべてプラズマ中の外に出てしまい，原子過制動放射，逆制動放射の計算では自由電子の速度分布関数

程に影響を及ぼさない一方，イオンはある一定の電子温度としてマクスウェル分布を仮定しているが，固体密度程度の熱浴に存在するとする．自由電子のエネルギーは輻射とになると自由電子の状態空間が制限され速度分布関数はなってプラズマの外に出ていくのだが，電子温度は変化し

フェルミ・ディラック分布となることを考慮する必要があないとすると平衡解が得られ，広くレーザー生成プラズマる．さらに温度が低い場合は，低いエネルギー準位に電子などの輻射流体シミュレーションで用いられている．図２がつまってしまうので束縛自由遷移閾値付近への遷移がで（b）に衝突輻射平衡での金プラズマの電離度を示す．密度

きなくなり，吸収係数のスペクトルにおける吸収エッジが６４３

Journal of Plasma and Fusion Research Vol.88, No.11 November 2012

だマイクロフィールドが計算できる［１７］．水素や水素様イ

なまってしまう．
特に高温プラズマ中で束縛状態がないイオンが多い場

オンの場合は電場があると縮退が解けてラインスペクトル合，自由自由過程しか輻射を放射，吸収することができなが分裂し，比較的大きなラインスペクトルシフトをもたら

いが，さらに大きなエネルギーの輻射を扱う場合は，トム

す．アルカリ金属などは縮退していないので弱い電場の場

ソン散乱の影響が無視できなくなる．通常の実験室で生成合は２次のシュタルク効果でラインスペクトルがシフトすされるプラズマ中の輻射によるエネルギー輸送の解析にトるが，さらに強い電場になると角運動量が適切な量子数でムソン散乱を考慮する必要はない．

なくなり，水素様イオンのようにシフト，分裂する．単純な水素やアルカリ金属などのシュタルク効果は比較的，計３．
６ ラインプロファイルと多群近似

算や実験データなどが得られるが［１８］，複雑なイオン，特平均イオンモデルのラインプロファイルの問題は大きくに原子番号が大きい部分電離イオンになるとそもそも多重２つに分けられる．一つは，特に原子番号の小さい原子か項に起因する分裂が元々プラズマ中のイオンからのライン

らなるプラズマの場合である．平均イオンモデルで計算す

スペクトルが備えているドップラー広がりや自然広がりでるイオンは非整数の束縛電子を持つ仮想的なイオンについマージしてしまっている場合も多く，輻射輸送を計算するて原子過程を計算する．そのため，実際のプラズマ中での上では個々のラインを考慮する必要がない場合が多い．

ラインスペクトルはイオン化状態の異なった複数のライン

X 線レーザーなどの利得計算を行う場合は，非常に幅のスペクトルとして現れるはずなのだが，平均イオンモデル

狭いラインスペクトルの輸送を計算する必要があるため，で計算されるラインスペクトルは１本で，温度・密度に対ラインスペクトルの詳細なプロファイルを計算することがしてラインの位置が連続的に変化する．後に説明する平均要求される．また速度勾配をもつレーザー生成プラズマを

イオンモデルの占有密度から各イオン化状態の占有密度を媒質とする場合は，プラズマの局所的な速度がドップラー

予測するなどして現実のイオンのラインスペクトルの位置効果としてラインスペクトルのシフトを起こす．

で放射係数，吸収係数を計算する必要がある［１５］．とは

輻射輸送を多群近似で解く場合には取り得る輻射のエネ

言っても現在では，特に原子番号の小さい原子からなるプルギー群の数に限りがあるため注意が必要である．たとえ

ラズマの場合，各イオン化状態を別々にして解くことが可輻射のエネルギーメッシュをたくさんとれるとしても輻射能であるので，平均イオンモデルを使用する理由はない．

スペクトルに特徴的なラインスペクトルが現れる場合は，もう一つは金などの原子番号が大きい場合である．当エネルギー輸送として考慮すべき輻射エネルギーの範囲で然，この場合でもラインスペクトルの位置の問題はある

その特徴的なラインスペクトルのプロファイルを記述できが，比較的密度が高い場合は，各イオン化状態のラインがるだけのエネルギー群を取ることは難しい．ラインスペク

その他のライン広がりによりマージして一つの UTA にみトルが特徴的な閉殻構造に近いイオン種からなるプラズマ

える．この場合，この広がり幅をイオン化状態の異なるイの場合，比較的エネルギー輸送における輻射の割合は小さ

オンの占有密度分布より求めればよい［１６］．

いので各輻射のエネルギー群を代表する放射係数と吸収係

STA のように考慮する状態数を最適化する場合も各ラ数で輻射輸送を解き，詳細なラインの輸送はポストプロセ

インスペクトルの位置は，ライン群を代表する位置にあスとしてそのラインのある輻射エネルギーの部分だけ詳し

り，通常は多重項による広がり幅をガウス分布で近似してく光線追跡などの手法で解く．

いる［１０］．

具体的に多群近似で計算するのに使用する放射係数と吸DCA モデルで多重項もフルに計算した DTA モデルの場

収係数の各エネルギー群での値はどう計算すればよいか？

合，ラインスペクトルの広がり幅は，イオンの熱運動に伴スペクトルがエネルギー群内で一様に増えたり減ったりしうドップラー効果，準位の寿命に起因する自然広がり幅がている場合は中央値の値をとればよいが，ある一つのエネ

あるが，ある程度密度の高いレーザー生成プラズマではプ

ルギー群の代表する輻射エネルギー帯域に，例えばライン

ラズマ中の統計的に分布する電場に起因するシュタルク効

スペクトルのように局所的に強い強度を持つ場合は注意が果に起因するラインスペクトルの分裂，シフトの方が一般必要である．放射係数についてはパワー（エネルギー）で的に寄与は大きい．また磁場閉じ込め核融合のプラズマであるのでエネルギー群内の放射係数スペクトルを積分しそはゼーマン効果によるラインスペクトルの分裂，シフトがのエネルギー群内に入るパワー（エネルギー）とすればよ重要となる場合も多いが，磁場自体は，外場やプラズマ中

い．一方，吸収係数の場合は，エネルギー群内の輻射エネの電流による磁場の大きさで評価される．一方，電場に関ルギーのスペクトル吸収係数の逆数（とプランク分布の強してはプラズマ中のイオン，電子によるためプラズマ中で度）で平均をとるロスランド平均をとったものをそのエネは必ず存在する．最も簡単なモデルはホルツマーク場であ

ルギー群の吸収係数とする．そのエネルギー群の吸収係数

り式（６）で定義されるホルツマーク基準電場程度の電場がスペクトルのまま平均をとるとプランク平均に相当するもプラズマ中には必ず存在している．イオンの運動エネルのが出てくるが，一般的な実験室プラズマではラインスペ

ギーとイオン間のクーロンポテンシャルの比である !が１

クトルの割合にもよるがプランク平均の方が１桁以上大きより十分小さい場合は，デバイ・ヒュッケルの２体分布関く，輻射輸送方程式での吸収係数を過大評価してしまう．

数を用いた APEX（Adjustable Parameter EXponential ap-

素性のよいモデル，例えば STA などで計算したプランク

proximation）を用いてクーロン相関を精度良く取り込ん

平均やロスランド平均吸収係数の値は DTA で計算した結６４４

Lecture Note

3. Atomic Process Modeling and Radiation Transport

T. Nishikawa

占有密度を計算している．

果とはおおむね一致する［１９］．

ここではさらに複雑な状態を含むモデルの場合は安易に

３．
７ UTA
（Unresolved Transition Array）と占有密度（ポピュレーション）の予測簡略化した占有密度から詳細な占有密度を予測する難しさ

部分電離金などの複雑なイオンからなるプラズマでは電モデルの結果を示す［２０］．従来の平均イオンモデルでは主子数が多く，取り得る状態数が非常に多いので，先に説明

量子数を考慮した状態しか考慮してこなかった．レーザー

した平均イオンモデルなど，何らかのモデルを使って計算

生成金プラズマの実験スペクトルを再現するためには，方を減らすことが行われてきた．その一方で，その減らした位量子数を考慮した状態を含めた平均イオンモデルを開発

各イオン化状態のイオン数の情報など，簡略化で失ってしする必要があった．図５で破線のシミュレーションは主量

まった情報を取り戻そうとする試みも数多く行われてき子数のみを考慮した原子過程モデルを用いて計算した結果

た．熱平衡状態のプラズマの場合は，平均イオンモデルので，スペクトルは 500 eV から 1 keV あたりの!"!!"の遷の例として，著者らが行った方位量子数を含む平均イオン

与える占有密度（ポピュレーション）をその準位に電子が移に伴うライン群と200 eV から 400 eV あたりの !""!#の存在する確率とみなし，各イオン化状態の占有密度を計算

遷移に伴うライン群からなっている．一方，実験で得られすることが行われている．例えば図３はアルミニウムのプているスペクトルは，300 eV 付近の比較的狭い範囲にまとラズマの計算結果である．実線は DCA モデルによる結果

まって出ている 4d-4f 遷移に伴うライン群と !"!!"遷移

で破線は平均イオンモデルから予測した結果である．熱平が副殻を考慮することにより，400 eV から 1 keV 程度まで衡モデルの場合，この手法はイオン化ポテンシャルがイオ

比較的広い範囲に分布するライン群からなるスペクトルでン化状態によらず一定とする近似で正しいことが証明され構成されている．

ている［１６］．

研究の当初は，遮蔽水素モデルで方位量子数を含んだエ

一方，衝突輻射平衡の場合も，主量子数のみを考量したモデルで，励起状態の占有密度から各イオン化状態の占有密度を計算する方法が提案されている［１５］．衝突輻射平衡における様々な密度の金プラズマの !"!から１０まで準位の占有密度をそれらの状態数と空虚度（Vacancy）で割った値を図４に示す．図４で密度が1.0（g/cm3）の場合は，各準位の値は片対数で直線上にのり，熱平衡状態の占有密度と一致し，傾きの逆数が電子温度になる．一方，衝突輻射平衡の場合，密度が小さくなると輻射脱励起の影響が大きくなる．輻射脱励起の速度係数は遷移エネルギーの２乗に比例するので下準位ほどずれが大きくなり，励起状態占有密度から計算した実効的な温度は実際の電子温度よりも小さくなる．一方，密度が低くても高い励起状態の占有密度
は準位間のエネルギー差が小さくなるため衝突過程の寄与が 大 き く な り ボ ル ツ マ ン 分 布 に 近 づ く．Hybrid atom

図４

平均イオンモデルを用いて衝突輻射平衡を解いた金プラズ
マのポピュレーションを状態数と空虚度で割った値のエネ
ルギー準位依存性．平均イオンモデルでは，熱平衡状態の場合，この値の傾きの逆数が電子温度となる．

図５

レーザー生成金プラズマからの軟 X 線スペクトルの実験結果とシミュレーション．

model では励起状態の占有密度がイオン化状態が異なっても同じ傾き，すなわち実効的な温度を持つとして励起状態の占有密度を計算し，その占有密度から各イオン化状態の

図３

アルミニウムプラズマのイオンアバンダンスのDCAモデルと平均イオンモデルからの予測の比較．

６４５

Journal of Plasma and Fusion Research Vol.88, No.11 November 2012

ネルギー準位は計算できたが，方位量子数を考慮した状態

クトルとして観測される UTA のエネルギーが準位の状態を含む平均イオンモデルの速度（レート）方程式の平衡解数で平均をとったラインの位置ではなく A 係数の重みをを解くことができなかった．そのため，次善策として，とって計算したライン群の位置に現れることを示し，A

図４で得られた衝突輻射平衡での励起状態の占有密度から係数の重みを取ったライン群の位置と分散を計算する手法

得られる実効的な温度を用いて，占有密度の方位量子数分であった［２１］．図８に 4d の電子数の関数としてスズの 4d布を予測することを思いついた．ところが，主量子数が異-4f遷移エネルギーの例を示す．4dと4fの各状態の状態数でなる準位間から得られた実効的な電子温度で方位量子数分

重みをとったエネルギー準位差をプロットした物が黒三角布を計算すると 4d-4f に対応するライン群の強度が小さす

で 4d-4f の遷移エネルギーを遷移の A 係数の重みで平均をぎ，実際の電子温度を使用すると大きすぎた．$##の状態取ったものが白丸である．図からわかるように，4d 殻の電間は $が異なる状態間よりエネルギー差 が 小 さ いので子数が変化しても 4d-4f の遷移エネルギーはあまり変化し$##の状態間の実効温度は実際の電子温度と衝突輻射平

ないため，比較的広い温度・密度領域で同じエネルギーの衡の実効的な温度の間にあると考え，$##間の占有密度

輻射を効率よく生成することができ，EUV 光源などの研を計算する温度として究・開発が精力的に行われている．図から各状態の状態数!#"!&"$!&''"!"%

で重みを取ったエネルギー準位差とA係数の重みで遷移エ（１２）

ネルギーを計算した値の差は，開殻4dの電子数が状態数のとし，"!
%とすると図６のスペクトルのように4d-4fに起#!

半分程度，入った時が一番大きく，10 eV 程度もある．レー

因するライン群と$##!$の遷移に起因するライン群との強度比を再現することができた．しかしながら得られたス
ペクトルには 400 eV，800 eV，1 keV あたり不自然な吸収スペクトルが現れた．後で方位量子数分布を含んだ平均イ
オンモデルの占有密度を解いた後でわかったことは，これらの吸収スペクトルは，占有密度の方位量子数分布を実効的な温度から予測で計算したために現れた，実際より過大評価された吸収線群であった．図７に示したフルに方位量子数を含んだ平均イオンモデルを解いた吸収係数，放射係数で計算した結果では，当然のことながら不自然な吸収は現れていない．実際，のちの計算で，"は温度と密度の関数
で，一定の値を取るわけではないことがわかったのだが，その依存性を手間をかけて探すよりは，可能ならば素直に複雑な式をそのまま解いた方がいい．

図７

方位量子数分布を考慮した平均イオンモデルを解いて計算
した結果．

図８

スズの 4d-4f 遷移エネルギー．状態数の重みで平均を取った準位間のエネルギー差と振動子強度の重みで平均をとったエネルギー差の比較．

さらに DTA から例えば $，#の量子数で記述する UTAを考えた時は，さらに悩ましい問題があることを示そう．
Bauche らによって有名になった UTA は，元々は分光スペ

図６

主量子数のみを考慮した平均イオンモデルを用いて計算したポピュレーションから方位量子数分布を励起状態の傾きから求めた実効的な電子温度の関数として計算した結果．

６４６

Lecture Note

3. Atomic Process Modeling and Radiation Transport

T. Nishikawa

ザー生成プラズマの場合，4d が開殻となる電子温度はせい

まう．通常は，電離度が変わると強い UTA の位置（エネル

ぜい 50 eV 程度である．さて，この場合，レート方程式をギー）が移動するため，比較的低温で密度が大きい部分で解く際のエネルギー準位としてはどちらを使うべきだろう

生成された輻射は温度が高く，密度が小さい領域を通り抜

か．状態数で各準位の平均を取った物をエネルギー準位とけるが，4d-4f 遷移の場合は UTA のエネルギーが大きく変してレート方程式を解くのが最初に考えられる．しかしな化しないため問題が起こる．このような場合に，Novikov

がら，このエネルギー準位間の差として計算したラインス

らは衝突輻射平衡と熱平衡時の放射係数と吸収係数の輻射ペクトルの遷移エネルギーは実際に実験で得られるライン場を特徴づける変数 &を用いて線形補間することにより簡スペクトルの位置とはずれている．それでは，A 係数を重易的に光励起の効果を取り入れる方法を提案した［２３］．次

みとして平均を取った遷移エネルギーを元にエネルギー準

世代リソグラフィー用の EUV 光源の高効率化のためのシ位を構成し，レート方程式を解いたらよいと思うかもしれミュレーションで実績を上げている［２４］．熱平衡時の放射ない．ところが，厄介なことに，A 係数を重みとして平均係数と吸収係数，衝突輻射平衡で計算した放射係数と吸収を取った遷移エネルギーでは，例えば，4d-4f の遷移から計

%#
"$#
%#
"$#
と し，&を係 数 を そ れ ぞ れ，$
% ，$
% ，#
% ，#
%

算した 4f のエネルギー準位と，もう一つの重要な遷移であ

&#!%$
%
!
$
"
#%$
"&$
%
!
$
%
%
%で表される，輻射輸送を解いて得る4d-5pから5p-5d，5d-4fと計算したエネルギー準位が一致られる各場所での光子の周波数 %における輻射エネルギー

しないのでレート方程式を解くべきエネルギー準位を構成密度とプランク分布の輻射エネルギー密度の比とした時，できない．実際，平均イオンモデルを用いて計算した際は

%#
"$#
$
$
$
!!&%
%#&
% "$
%

［２２］，4d-4f，4d-5p 等，実験のスペクトルを構成する重要

%#
"$#
#%
#%
#%#&
!!&%
"$

な遷移間のエネルギーがなるべく一致するようにエネル（１３）

ギー準位を構成する遮蔽定数を計算し，平均イオンモデルとして輻射輸送を解く．ただし，ここで熱平衡時の値はイでレート方程式を解いた．逆に，状態数で平均をとったエオン数密度が同じで電離度が衝突輻射平衡の値に一致す

ネルギー準位でレート方程式を解き，その占有密度でエネ

る，少し温度の低い熱平衡時の値を用いる．同じ電離度にルギーがずれた遷移エネルギーを使用して放射係数，吸収なる場合の値を用いるため，これらの放射係数と吸収係数係数を計算すれば，原子過程と輻射過程のつじつまが合わを用いて計算したプランク分布は少し温度が小さい場合のず，図６でみたように不自然な吸収や放射が起こってしま

ものに漸近するが，大きな問題にはならないようである．

うことは容易に想像できる．

同じ温度の熱平衡の占有密度は図２で示したように電離度

占有密度，特に衝突輻射平衡などの非平衡状態の占有密が大きく違うのでうまくいかない．あくまで輻射励起によ度を計算するために簡略化した状態モデルを使用した場り変化するのは励起状態数の割合であり，プラズマの電離合，その占有密度から元の簡略化する前の詳細な占有密度度が違うと励起状態の占有密度自体が大きく変わってしまを出すことは原理的には可能であろうが，十分，注意深くい，スペクトル強度がまったく異なったものになってしま検討が必要で，現在のように計算機資源が豊富な時代でうため，わざわざ電離度の一致する熱平衡の値を使用すは，もとの大規模な状態モデルで占有密度を解くことをおる．通常，流体シミュレーションコードでは輻射輸送を解すすめする．

く際に使用する放射係数と吸収係数はあらかじめ計算した数表から補完するなりして特定の温度・密度の値を取り出

３．
８ 光励起過程している．そのため熱平衡と衝突輻射平衡の２つの数表を

ここでは，輻射励起過程を簡易的に取り入れた原子過程

用意しておき，輻射場から計算した &と（１３）式で放射係数を流体運動と共に解くことができるモデルを紹介する．あと吸収係数を計算することはそれほど大きな計算負担とはる程度光学的厚さのプラズマスケール長を伝搬し，UTAならない．

スペクトルの強度がプランク分布における輻射エネルギー

参考文献

密度程度になる，典型的には，レーザー生成プラズマなどで応用できる．

［１］H. Mayer, Los Alamos Scientific Laboratory Report, LA647 (1947).
［２］R.M. More, J. Quant. Spectrosc. Radiat. Transfer 27, 345(1982).
［３］例えば，D.G. Hummer and D. Mihalas, Astrophys. J. 331,
794 (1988).
［４］W. Däppen, L. Anderson and D. Mihalas, Astrophys. J. 319,
195 (1987).
［５］T. Nishikawa, Astrophys. J. 533, 670 (2000).
［６］T. Nishikawa, Plasma Fusion Sci. 7, 1401142 (2012).
［７］W.A. Lokke, and W.H. Grasberger, Lawrence Livermore
National Laboratory Report No. UCRL-52276 (1977).
［８］D.E. Post, R.V. Jensen, C.B. Tarter, W.H. Grasberger and
W.A. Lokke, At. Data Nucl. Data Tables 20, 397 (1977).

次世代半導体リソグラフィーの光源として注目されているスズプラズマからの 4d-4f 遷移からなる EUV 光は図８に示すようにイオン化状態が異なっても同じ位置に現れる．
そのため，輻射流体シミュレーションにおいてもプラズマ
中での EUV 光の自己吸収が大きく，光励起を何らかの形で取り扱う必要がある．しかしながら，輻射場を考慮した原子過程をインラインで輻射流体計算と同時に解く事は，
まだ計算機能力の点で難しい．従来からよく用いられている衝突輻射平衡の放射係数と吸収係数をそのまま用いると特にレーザーが照射されているプラズマ表面では密度が小
さく，温度が高いので輻射脱励起の寄与が大きく，下準位の占有密度が大きくなり，結果，吸収が過大評価されてし
６４７

Journal of Plasma and Fusion Research Vol.88, No.11 November 2012［９］T. Nishikawa, H. Takabe and K. Mima, Laser Part. Beams
11, 81 (1993).
［１０］A.Bar-Shalom, J. Oreg, W.H. Goldstein, D. Shvarts andA. Zigler, Phys. Rev. A 40, 3183 (1989).
［１１］W. Lotz, Z. Phys. 216, 241 (1968).
［１２］H.A. Kramers, Philos. Mag. 46, 836 (1923).
［１３］D.R. Bates, A.E. Kingston and W.P. McWhirter, Proc. R.
Soc. London, Ser. A 267, 297 (1962).
［１４］例 え ば，Ya. B. Zel'dovich and Yu. P. Raizer, Physics ofShockwaves and High Temperature Hydrodynamic Phenomena (Academic Press, New York, 1966) Vol. 1 Chap. V.
［１５］M. Itoh, T. Yabe and S. Kiyokawa, Phys. Rev. A 35, 233(1987).
［１６］西川 亘，中村真嗣，高部英明，三間圀興：核融合研究
６
８別冊，145 (1992).
［１７］C.A. Igleisia, J.L. Lebowitz and D. MacGowan, Phys. Rev.

A 28, 1667 (1983).
［１８］例 え ば,M.G. Littman, M.L. Zimmerman, T.W. Ducas, R.
R. Freeman and D. Kleppner, Phys. Rev. Lett. 36, 788(1976).
［１９］3rd Int. Opacity Workshop & Code Comparison Study
WorkOp-III:94 Final Report, MPI für Quantenoptik,
Garching, MPQ 204 (1995).
［２０］T. Nishikawa, H. Takabe and K. Mima, Tech. Rep. Osaka
Univ. 41, 253 (1991).
［２１］C. Bauche-Arnoult, J. Bauche and M. Klapisch, Phys. Rev.
A 20, 2424 (1979).
［２２］K. Nishihara, T. Nishikawa et al., Proc. 3rd IFSA (2003).
［２３］A.F. Nikiforov, V.G. Novikov, V.B. Uvarov and A. Iacob,
Quantum-Statistical Models of Hot Dense Matter: Methods for
Computation and Equation of State (Birkhauser, 2005).
［２４］K. Nishihara et al., Phys. Plasmas 15, 056708 (2008).

!!!!!!!!!!!!!!!!!

かわ

たけし

西 川

にし

亘

!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

岡山大学大学院自然科学研究科，助教．大学院から一貫してプラズマ中原子過程の理論・シミュレーション研究に従事．学部は電気電子工学コースに所属するが，「電気」
とか「電子」のつく学会には一つも所属していない不心得者．
とはいえ昨今，世間で何かと話題に上る電気については大学で関わる身として多少なりとも知識は身につけるべく，出先で見つけた電気設備をいろいろ写真に収めています．
趣味はマラソン．

!!!!!!!!!!!!!!!!!!!!!!!!!!

６４８

J. Plasma Fusion Res. Vol.88, No.11 (2012)６
４９

講座 輻射流体シミュレーション

４．おわりに
長友英夫
大阪大学レーザーエネルギー学研究センター（原稿受付：２
０
１
２年１
０月３
１日）

輻射流体シミュレーションは輻射輸送と流体運動という時定数の異なる現象を結合して解くということから両者を厳密に解くことはむずかしく，その近似解法を適応すると多様になる．この分野の先駆的な研究分野の一つであるレーザープラズマ・核融合における輻射流体に関しては，これまでに本学会誌の小特集［１］などで紹介する機会があったので，今回は違った視点からの講座とした．また，流体解法の部分についても，数値（計算）流
体力学としての研究分野が幅広く認知されていることから割愛し，輻射輸送とその流体との結合に絞った．特に，
輻射輸送に主眼をおいたので，輻射輸送方程式を直接 SN 法で，あるいは拡散近似して解くなどの流体に起因する
光学的な特性を反映させた手法について，それぞれ第１章，２章で詳しく記述しているので理解していただけたのではないかと思う．
原子過程は輻射と流体を結びつける放射・吸収係数を決める重要なプロセスである．輻射流体シミュレーショ
ンでは，通常，予め用意された数表を参照するが，数字が羅列されているだけの数表からその作成するための労力は測りにくい．第３章ではその順を追った計算方法だけでなく適応範囲，注意点など流れがわかりやすく示してあり，大いに参考になると思う．実際に必要なデータベースをすべて自分で揃えるのは難しいし，いくつかのデータベースは文献［２］，Web サイト［３］等で公開されている．実際の輻射流体シミュレーションで利用するには少しデータを加工する必要があるが，公開されているデータは整理されていて検索しやすいだけでなく，広く使われていることから信頼性も高いと思われる．
本講座では，現象を支配する主たる物理の概要とその支配方程式からのシミュレーション方法を紹介したが，実際のシミュレーションを実施するまでにはまだステップがある．要求される計算精度，計算機資源に応じて適切な解法を採用しなければならない．例えば，拡散近似法，SN 法などは，最終的な求解部分は，大規模疎行列からなる連立方程式を解くことに行き着くので，これら解法の高速化，近年は特に並列化などのコーディングが必要になってきている．また，数値計算ライブラリや計算機技術に依存することから，その進展に追随していなけ
ればならない．この点は第１章でも述べられているとおり，分野横断的な組織による効率的な研究開発が求められており，実際，米国ではエネルギー省のプロジェクトとして位置づけられ開発研究が進められている．特に，Advanced Simulation and Computing（ASC）の中の Academic Strategic Alliance Program（ASAP）の枠組みで開発された Flash Code［４］は有名で，学術研究分野にも開放されている．一方，国内ではそのようなプロジェクト
開発は難しいのが現状だが，当分野以外でも類似した研究をされている方の協力，情報共有が重要だと思う．そのような観点からもこの講座をそのきっかけにしていただけたら幸いである．

参考文献［１］長友英夫他：プラズマ・核融合学会誌 82, 141 (2006).
［２］浜口智志 他 編：プラズマ原子分子過程ハンドブック（大阪大学出版会，２０１１）．
［３］村上 泉：プラズマ・核融合学会誌 88, 180 (2012).
［４］http://flash.uchicago.edu/site/flashcode/

4. Conclusion
NAGATOMO Hideo

author’s e-mail: naga@ile.osaka-u.ac.jp

６４９!2012 The Japan Society of Plasma
Science and Nuclear Fusion Research