J. Plasma Fusion Res. Vol.90, No.8 (2014)５００‐５０４

講座 高密度相対論プラズマの粒子シミュレーション技法

６．粒子シミュレーションにおける輻射輸送と核反応モデル
6. Radiation Transport and Nuclear Reaction Models for Particle Simulations
城 " 知 至，千 徳 靖 彦１）
JOHZAKI Tomoyuki and SENTOKU Yasuhiko１）
広島大学大学院工学研究院，１）ネバダ州立大学リノ校物理学科（原稿受付：２０１
４年６月２３日）

高密度相対論プラズマでは，熱輻射や高エネルギー電子からの放出 X 線・&線によるエネルギー輸送が重要となる場合がある．また，高エネルギー粒子や光子による核反応過程も物理過程として非常に興味のある現象である．本章では，これらの過程を粒子コードに組み込む手法について紹介するとともに，いくつかの研究成果について紹介する．
Keywords:
radiation transport, nuclear reaction, particle simulation

６．
１ はじめに生成過程）を計算することを目的とし，輻射強度に対する輸送方程式を数値的に解く輸送計算法を紹介する．
++
輻射場は光子の分布関数 #(!
&!
$!
(&'で記述される．ここ

本講座ではレーザープラズマ相互作用領域から高密度固体密度領域までの幅広い密度領域を取り扱うためのシミュレーション技法を紹介してきた．最終章となる本章では，+，(はそれぞれ，光子の周波数，位置ベクトル，+で，(，&，$高エネルギー密度科学や相対論レーザープラズマで重要と方向ベクトル，時刻である．輻射の輸送方程式は，この光なる X 線や &線によるエネルギー輸送過程，および高密度物質中での核反応過程を粒子コードに組み込む手法とその子の分布関数に対するボルツマン方程式から導出され，輻
++++
射強度 !
（添え字 (は周波数，"，&!
$!
(%(##
(!
&!
$!
('$"&'(&応用例について紹介する．

%は光速とプランク定数を表す）に対して，次式となる．
! $! ,+#%!
! (""("$'(!%(!
((" ,６．
２ 輻射輸送相対論的強度のレーザーと物質の相互作用では，相対論（１）(((+++++
$'
$!
$'
!
$'
#$ $．
!&#(&!*'(#!
(&(&的高エネルギー電子が生成され，物質中にエネルギーを輸送することで，高エネルギー密度プラズマが生成される．

")

高エネルギー電子により加熱された物質内では，バルクか

ここで，'(，%(，*'(は放出，吸収，散乱係数を表す．また，ら輻射される X 線によるエネルギー輸送が重要となる場合

右辺第三項は散乱項で，括弧内の第一項は着目している方+(から+からの散乱による流出，第二項は任意の方向 $向$がある．また加熱物質からの特性 X 線等の放出・輸送特性
介されたように，物質中を伝播する相対論的電子からの制+への流入を表し，#(&++('は %++('+$!とな$$!
$$!
$#$")#(&(+ から $+への散乱確率である．以下では簡単化のため，る$

動放射 &線や，1022 W/cm2 を超えるような超高強度レー散乱の効果は無視する．

はプラズマ状態を知るうえで重要である．一方，５章で紹

ザーと電子の相互作用における輻射減衰による &線放射，

輻射輸送方程式の解法としては，拡散近似や光線追跡法

さらに電子・陽電子対消滅による &線放射等は，高エネル等があるが，輻射場の非等方性が強く，かつ非定常な場をギー密度科学において重要な物理現象の一つである．これ解くために，ここでは輸送方程式を直接離散化して解く方

ら X 線や &線は，その波長がレーザー波長に比べ３ケタ以法を用いる．輻射強度は，実空間方向の次元数に加え，周上短くなるため，粒子シミュレーションにおいて電磁波と波数と伝播方向の次元をもつ．周波数については，計算対して Maxwell 方程式に基づいて解くことは現実的ではな象とする周波数領域を有限数の周波数群に分割し，輸送方い．また X 線輸送を解くためには，放出・吸収過程等の原程式を各周波数群の幅 "($で積分し，周波数群毎の輻射強子過程を含めた計算が必要となる．ここでは，粒子シミュ度!
(に対する輸送方程式とする（多群化）．この$ $%(#"( !

レーションにおいて，X 線や &線放射によるエネルギー輸際，群の切り方は，等間隔である必要はなく，問題に応じ送や，放射と粒子の相互作用（例えば電子陽電子対消滅・て適切に切る．例えば，バルクからの放射で連続成分と線$Graduate School of Engineering, Hiroshima University, Higashi-Hiroshima, HIROSHIMA 739-8527, Japan
Corresponding author’s e-mail:tjohzaki@hiroshima-u.ac.jp

５００!2014 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

6. Radiation Transport and Nuclear Reaction Models for Particle Simulations

T. Johzaki and Y. Sentoku

成分とが混在する場合は，線成分領域については細かく切対して各超粒子の重みに応じて分配する．結合計算の際，り，連続成分領域については粗く切る等の工夫をすること

輻射輸送における格子構造を粒子コードと同じものとしで，少ない群数で精度のよい計算が可能となる．伝播方向)については，## 法［１］等の離散化法により，適切な離散"て，毎ステップ輸送計算を行うことも可能である．しかし，方向を取り，各離散化方向毎に輸送方程式を解く．この時，)に対する積分は各離散化方向毎の重みをかけて伝播方向"ため，現象に応じて輻射輸送の計算格子を粗くし，計算負総和を取ることで表す．##法による離散化の詳細について能である．各周波数群における放出・吸収係数は，バルクは参考文献［２］を参照されたし．また，外部からの輻射ソーの組成に対して温度・密度の関数として予め作成しておス（例えば XFEL による X 線と物質の相互作用など）の場

く．

輸送計算は計算負荷，特にメモリーの使用量が大きくなる荷を低減する（粒子計算の数ステップ毎に行う）ことも可

合は，その入射方向に離散方向を取る（もしくは集中させ衝突・電離過程を含む相対論電磁粒子コード PICLS［４］る）ことで，ビーム状の輻射場を解くことも可能である．に輻射輸送を結合して行った二次元結合計算の一例を示

周波数と伝播方向について離散化した×

す．固体密度の Cu ターゲット（初期電離度３を仮定）に高離散化方向数 個の輸送方程式は，散乱を無視した場合に周波数群数

強度レーザー（強度 1020 W/cm2，波長 1 μm，P 偏光，スは各離散化方程式ごとに独立なものとなり，個々の離散化

ポットサイズ 10 μm$，パルス長 330 fs フラットパルス）を方程式を実空間方向について解けばよい．方程式は左辺の照射した．本計算では，FLYCHK コード［５］を用いて，Cu移流項と右辺の生成・吸収項（非移流項）からなる．空間に対する非局所熱平衡（non-LTE）オパシティテーブルを変数については，電場や磁場と同様に空間格子点上に離散作成した．このときのレーザーの吸収率は４０％であった．

点 を 取 り，移 流 項 に つ い て は CIP 法［３］等 を 用 い て 解

図２にレーザー照射終了間際のターゲット内での輻射によき，非移流項は，各格子点における吸収・放出係数を用いるエネルギー輸送の空間プロファイル（a）と，レーザー照て輻射強度の変化量を求める．

射終了後に計測したターゲット内の輻射スペクトル（b）を示す．図２（a）においてレーザーは左から $"32 μm を集光

次に，粒子コードとの結合について考える．粒子コードからはバルクからの輻射や高エネルギー粒子による %線放点として垂直入射している．スペクトルは，（１）$"32 μm

出（制動放射や対消滅）が計算され，これが輻射輸送方程（レーザー照射軸上），（２）40 μm，（３）50 μm，（４）60 μm の式の生成項となる．一方，輻射場は空間を伝播する際に各各位置で x 軸方向に積分したものである．相対論レーザー

格子点におけるプラズマ状態により決まる吸収係数に応じプラズマ相互作用により生じた高速電子がターゲット内をて，輻射強度が減衰する（吸収項）．吸収により失われた輻

伝播・加熱するため，レーザー照射軸近傍で高温となる．

射エネルギーはプラズマに付与され，電離や温度上昇に費その領域では，輻射によるエネルギー損失が顕著である．

やされる．また，対生成の場合は，各素過程の断面積と %線並びに電磁場や原子核の分布に応じて生成率を計算する．ここではバルクプラズマとの結合について述べる．
図１に結合計算の模式図を示す．
粒子コードにて各格子点での温度・密度・組成を評価し，それを基に，予め作成しておいたオパシティテーブルデータから各格子点での放出・吸収係数$&#'%を求める．
'!
これを用いて，離散化した輻射輸送方程式を解く．求まった輻射強度と放出・吸収係数 $&#'%を用いて各格子点で'!
#)!&の輻射エネルギー損失率 "&$%"# !
#"%'を求'!
'%'""(#!

める．この値が各格子点におけるバルクプラズマの加熱率（負値の場合はプラズマからのエネルギー損失率）となる．
この値を粒子コードに返し，各格子内のバルク電子成分に

図２
図１

粒子コードと輻射輸送計算の結合模式図．

５０１

レーザー照射終了間際の輻射によるエネルギー輸送の空間
プロファイル（a）とターゲット内の輻射スペクトルのy方向依存性（b）．

Journal of Plasma and Fusion Research Vol.90, No.8 August 2014

その周辺領域は低温のため吸収係数が大きいので，中心加ミュレーションで行うために，核反応過程の導入が行われ熱領域で放出された輻射は周囲に伝播するものの，放出領ている．本節では，これらの核反応過程のうち，レーザー域から 10 μm 程度の領域でほとんどの輻射が吸収されてい加速イオンとターゲット内原子核との核融合反応過程を取ることがわかる．ターゲット外への輻射放出はレーザー照り扱う方法とともに，研究の一例を紹介する．

射面からの放出がほとんどで，その割合は照射レーザーエ相対論レーザープラズマ相互作用で加速されたイオンネルギーの１％であった．一方，ターゲット内のスペクトが，固体ターゲット内に入射し，ターゲット内の標的核種

ル（図２（b））を見ると，光子エネルギー 2〜8 keV におけとビーム核融合反応を起こす場合を考え，核融合反応の生る連続線がレーザー照射軸上の最も温度の高い領域で強起確率と生成中性子の放出角分布を評価する方法を述べ

く，軸から離れるにつれて弱くなっている様子がわかる．

る．ターゲット内に入射する高エネルギーイオン（質量また，1〜2 keV 並びに 8〜9 keV 領域に線スペクトルが見)-，速度 ,-，エネルギー #-）とターゲット内の静止した標えており，レーザー軸中心近傍の高温領域では，比較的電的核種（質量 )5）との反応を考える．この時の核融合反応離が進んでLi様Cuからの放出線スペクトルが見えている．の生起確率 $-' 5$#-%は，標的核種の数密度を *5，反応断

また，軸から離れるにつれて温度が低下するため高電離状，時間ステップ幅を #+面積を *-' 5$（必ずしも粒子計#-%

態からの線スペクトル強度が低くなっていくという，特徴算の時間幅に合わせる必要はなく，計算負荷を考えその数倍にしてもよい）とすると，$-' 5$#-%#-%#+"*5,-*-' 5$的な線スペクトルが示されていることがわかる．
このように，粒子コードに輻射輸送を結合することで，で計算される．標的核種の各格子内の密度は，本講座２章

ターゲット内での輻射によるエネルギー輸送を考慮できるで紹介した高次補間法により，電流密度等の計算と同じ手と共に，ターゲットから放出される輻射スペクトルや放出法で評価する．核融合反応断面積は，例えば Bosch と Hale

分布を評価し，実験結果と比較することで，相対論レー［６］に よ っ て 重 心 系 の 相 対 運 動 の 運 動 エ ネ ル ギ ー

ザー物質相互作用の理解がより進むことが期待される．ま#), "#-)5")-!)5%に対して，次式で与えられている．
$

た，相対論レーザープラズマ実験において，+%線計測によ*$#), %"り高速電子の伝播特性等が評価されているが，シミュレーションに高速電子による+%線放出過程を組み込むことで，実験との直接比較が可能となる．更には高速電子由来の '%$#), %．
#), 062$"*"#! ), %（２）

ここで，"* は Gamov 因子で，線の輸送や X 線自由電子レーザーの X 線輸送等も，適切な"* "(%&-&5&#)3'#，周波数領域，離散化方向を設定し，放出・吸収係数を準備（３）!'%であり，%は微細構造定数（(#""!
"$(!$'!%）で，&-，""$することで可能となる．このように輻射輸送過程を粒子コードに組み込むことで，高輻射場における非平衡プラズ&5は反応核種の原子番号，)3は反応核種間の換算質量で，マの研究など，様々な展開が期待される．

)3'#は keV 単位である（"*#を Gamov 因子として定義する場合もある）．%$#), %は宇宙物理学的 %ファクターで，文６．
３ 核反応過程

献［６］に次式で与えられている．

レーザー強度が 1018 W/cm2 を超える相対論レーザー領域では，レーザーの光子圧が Gbar を超え，電子のみならず"%$#), %イオンも加速され，MeVを超える高エネルギーイオンが生%%!"!#), $!#!#), $!$!#), $!%!#),!&%．
"!#), $%%""!#), $"#!#), $"$!#),"%%（４）成される．したがって，固体ターゲットの加熱や高速点火
レーザー核融合における爆縮コア加熱において，高エネル

係数 !"#!&，""#"%については代表的な反応について文

ギー電子だけでなく，高エネルギーイオンも寄与する．ま献［６］のTable IVを参照されたし．また，他の反応の反応断た，このような高エネルギーイオンが高密度ターゲット内面積については，NRL の Plasma Formulary［７］や日本原子に入射すると，ターゲット内の原子核との間でビーム核融力研究開発機構の核データ研究グループが提供している核合反応を起こす．実験的には，この核反応から生成する中データライブラリーの Web site［８］等から入手可能である．

性子を計測することで，ターゲット表面におけるレーザー次に，反応生成物の放出角を評価する．核融合反応で生加速イオン特性を評価することが可能であり，ターゲット成する粒子の放出角は等方ではない．重心系における微分

加熱に対する高エネルギーイオンの寄与を調べることがで断面積（/*"/$）は，Brookhaven 国立研究所の National Nu-きる．応用としては，レーザー加速イオンを標的物質に照

clear Data Center の Webサイトの評価済み核データセット

射した際に生じる核反応を利用した中性子ビーム生成など（Evaluated Nuclear Data File; ENDF）［９］から入手可能であがある．また，高Z材に超高強度レーザーを照射すると，発3る．例えば，入射・標的核共に重水素で，D（d，n）
He 反応で生した高エネルギー電子がターゲット内で原子核と相互作生成する中性子の場合，実験室系のエネルギー #/"100 eV

用して，もしくは制動放射により '線を放出し，その '〜10 MeV の領域に対して，反応断面積 *$#/%で規格化し線が原子核と相互作用することで電子・陽電子対を生成す

た /*"/$（すなわち，$ 方向への放出の確率分布関数）をる．この過程により電子・陽電子プラズマを生成し，宇宙放出角の極角の余弦（/$"/&/.
14)）に対して14)の .

物理の模擬実験に適用することなどが考えられている．こLegendre 多項式展開で近似した際の係数が与えられていれらの過程をレーザー照射から核反応まで一貫して粒子シる．実験室系での入射重水素エネルギー #/に対して，生成５０２

Lecture Note

6. Radiation Transport and Nuclear Reaction Models for Particle Simulations

T. Johzaki and Y. Sentoku

中性子の放出角余弦が (#+/0*の周りの ,(となる規格化

核融合反応エネルギーを %，着目放出粒子の重心系での放"確率分布 '("",%,(（$$("",%,(#"）が，$!"'

出角余弦を (&)$#+/0*，入射粒子の運動エネルギーを"("とし，&"#)")$## #，&##)#)%## #である．

) + ("#'("",%$#$ ",%",%+!!*$*(ビーム核融合などの生起確率が非常に小さく，反応によ（５）る粒子数の変化が無視できる場合は，反応による粒子の増(#""& $，"%$$(%# ( , (減は考えず，単に放出核種のスペクトル評価等を行うだけ!
!で，(#１〜６にの形で近似され，係数 &($",%は，&!#"数の変化を考慮すべきであろう．また，生成核種の飛程が対してテーブルとして与えられている．ここで，+''.$",%

計算体系と同程度の場合は，生成核種の輸送を解く必要がは反応断面積，+$("",%は放出中性子に対する微分断面積ある．反応による粒子の減少については，対象粒子の重みである．図３に D（d，n）3He 反応で生成される中性子の放を反応率に応じて減ずることで考慮できる．一方，粒子生出角の規格化確率分布を示す．

成の場合，予め生成核種用に計算配列を準備しておくか，##(#!

でよい．一方，反応確率が大きい場合は，反応による粒子

実際の微分断面積は，式（５）と反応断面積 +''.$",%を別途，輻射輸送の場合と同様に輸送方程式を解く等，工夫用いて，次式で評価される．
+ $"%("",%+$("",%$# ''. , ')#" %*( #+ $(""& $"%$$(%# ''. , #) (#! # ( , (#が必要となる．
上記の手法は，入射粒子とターゲット内標的核種との反応について述べたが，衝突計算において行う粒子間ペアリ．（６）ングを利用することで，高エネルギー粒子間，バルク粒子間の反応も取り扱うことができ，より一般的なモデルとなる．
3
相対論電磁粒子コード PICLS［４］に D（d，n）
He ビーム核一方，方位角方向'については 0〜2)で等方であるとする．
実際に計算を行う上では，まず実験室系から重心系への座融合反応を組み込んで行った二次元粒子シミュレーション

標変換を行い，重心系での放出角 *を，上述した確率分布の一例［１０］を示す．対象は阪大レーザー研で行われた PWM

関数に応じた乱数で 0〜)の範囲で決定し，方位角'についレーザー［１１］（波 長 1.054 μm，パ ル ス 長 500 fs［半 値 全ては０〜2)の一様乱数により決定する．最後に放出粒子の幅］，ス ポ ッ ト 直 径 5.4 μm［半 値 全 幅］
，集 光 強 度運動を実験室系に変換し，実験室系での放出角度 %，$2×1019 W/cm2）を CD フォイルに垂直照射した実験［１２］

並びにエネルギーを求め，放出角分布やエネルギースペクで，レーザー軸から 54.7°ならびに 90°方向で中性子スペクトルを評価する．実験室系−重心系間の変換については，トルが計測されている．図４に実験にて観測された中性子本講座第３章の２体衝突モデルと同じである．放出粒子の実験室系でのエネルギー "-$は，反応の運動学から次式で求められる．
"-$#&""-""&#!"-""# %")##"# &"&#"-"!"-"" %"(&) ．（７）)#%$

ここで，反応の前後で粒子の質量和の変化（質量欠損）は小さいものとして無視し，反応前の入射粒子を添え字１，標的粒子を２で，反応後の生成粒子を添え字３，４で表し，着目放出粒 子の添え字を３としている．また，粒子の質 量を ) で表し，反応前後の粒子の質量和を # #)"")##)$")%，

図３

3He 反応で放出される中 性 子 の 規 格 化 確 率 分 布D（d，n）．0°が入射粒子の運動方向．放出確率は入射粒子の（(f ，Ed）運動方向軸に沿って前後方に高くなっており，この傾向は入射粒子エネルギーの増加とともに強くなる．

５０３

図４CD フォイル（5.4 μm 厚）に PWM レーザー［１１］を垂直照射した実験にて，レーザー軸から 54.7°ならびに 90°方向で観
3He 生成中性子
測された中性子スペクトル．破線は D（d，n）の基準エネルギー 2.45 MeV．

図５3He 生成中
2D PIC シミュレーションより得られた D（d，n）性子スペクトル．

Journal of Plasma and Fusion Research Vol.90, No.8 August 2014
3スペクトルを示す．ピーク値は D（d，n）
He 反応で生成されクロ的に扱うことが困難で，多階層的シミュレーションがる中性子の基準エネルギー 2.45 MeV よりも高エネルギー行われている現象も，今後の計算機開発や計算手法の進展側にシフトし，かつ非常にブロードなスペクトルとなってによっては，すべてをミクロ的に取り扱うシミュレーショいて，ビーム核融合が支配的であることがわかる．また，ンが可能となるかもしれない．本講座において紹介した手［個／sr］の 中 性 子 が 観 測 さ れ て い
各 観 測 方 向 で７×１０４法が，今後，解析を行っていこうとする読者の一助となれる．一方，二次元粒子シミュレーションで得られたスペクば幸いである．

トルを図５に示す．シミュレーションで５５°方向で観測された中性子数は５×１０７［個／（sr・cm）］であった．単位長さ

参考文献

あたりとなっているのは #方向の一様性を仮定した二次元［１］C.E. Lee, Los Alamos Scientific Laboratory Report LA2595, (1962).
［２］小林啓祐：原子炉物理（コロナ社，１９９６）第６章．
［３］T. Yabe et al., Comput. Phys. Commun. 66, 233 (1991).
他，CIP 法に関しては，さまざまなバリエーションがあり，
矢部孝他：CIP 法―原子から宇宙までを解くマルチスケール解法（森北出版，２００３）等を参照されたし．
［４］Y. Sentoku and A.J. Kemp, J. Comput. Phys. 227, 6846
(2008).
［５］H.-K. Chung et al., High Energy Density Phys. 1, 3 (2005).
［６］H.-S. Bosch and G.M. Hale, Nucl. Fusion 32, 611 (1992).
［７］J.D. Huba, NRL plasma formula (NRL, Washington DC,2011). NRL web サ イ ト か ら 入 手 可 能 http://www.nrl.
navy.mil/media/publications/plasma-formulary/［８］日本原子力研究開発機構核データ研究グループの核データ web サ イ ト http://wwwndc.jaea.go.jp/index̲J.
html［９］Brookhaven 国立研究所の National Nuclear Data Centerの Web サイトの評価済み核データセッ ト(EvaluatedNuclear Data File; ENDF)http://www.nndc.bnl.gov/exfor/endf.htm［１０］千徳靖彦：超高強度レーザープラズマの計算機シミュレーションによる研究，大阪大学博士学位論文（１９９９）．
［１１］H. Fujita et al., Proc. SPIE 3047, 501 (1997).
［１２］N. Izumi et al., Phys. Rev. E 65, 036413 (2002).
［１３］T. Moritaka et al., J. Phys.: Conf. Series 454, 012016 (2013).

（!
-"平面）シミュレーションのためである．#方向について
4
レーザー集光サイズの２倍を仮定すると，5.4×10［個／sr］となる．z 方向サイズの任意性はあるが，実験結果とよい一致を得た．また，ピークのシフトやエネルギー広がりについてもよい一致を示しており，実験で観測された中性子の大半がレーザー加速重水素とターゲット内の重水素原子核とのビーム核融合反応由来であることが示された．

６．
４ まとめ
本章では粒子シミュレーションに輻射輸送と核反応過程を導入する手法について概説し，その適用例を個別に示した．最近では，QED 過程（制動放射と対生成）と生成 !
線輸送過程を粒子コードに組み込み，高 Z 材における相対論レーザープラズマ相互作用による陽電子生成過程のシミュレーションも行われている［１３］．このように，計算機の急速な性能向上に伴い，粒子シミュレーションの適用領域は，レーザープラズマ相互作用領域のみならず，高密度固体や爆縮プラズマ領域にまで及ぶようになってきた．それに伴い，衝突過程，電離過程に加え，輻射輸送や核反応過程などを含む物理モデルの深化した粒子コードが開発され，様々な研究領域に用いられるようになってきている．
また，現状では，時空間スケールが大きくて，すべてをミ

５０４