核データニュース，No.120 (2018)

部会・連絡会セッション

核データ部会 [「シグマ」特別専門委員会共催]「我が国における核データ計算コード開発の現状と将来ビジョン」(4) PHITS コードにおける核反応モデルの役割と高度化

日本原子力研究開発機構
橋本 慎太郎hashimoto.shintaro@jaea.go.jp

佐藤 達彦

sato.tatsuhiko@jaea.go.jp

岩元 洋介

iwamoto.yosuke@jaea.go.jp

小川 達彦

ogawa.tatsuhiko@jaea.go.jp

古田 琢哉

furuta.takuya@jaea.go.jp

安部 晋一郎

abe.shinichiro@jaea.go.jp

高度情報科学技術研究機構
仁井田 浩二1.

niita@rist.or.jp

はじめに粒子・重イオン輸送計算コード PHITS (Particle and Heavy Ion Transport code System)[1]は、中性子をはじめとする様々な放射線の振る舞いを模擬することができるため、加速器施設の遮へい設計や宇宙線・地球科学、放射線防護研究等の幅広い分野で利用されている。その応用利用の多くで中性子・陽子・重イオンが重要な役割を果たすことから、これらの粒子が起こす核反応を正確に記述するため、核反応モデルや核データライブラリの整備及びこれらの高度化は必要不可欠となる。PHITS の開発は、日本原子力研究開発機構を中心に多数の研究機関の協力を得て進められており、新規モデルの導入や既存モデルの改良を通じて、計算結果の信頼性の向上に努めている。
本稿では、PHITS で取り扱う物理現象について示し、それらを再現するために導入している各モデルやデータライブラリの役割についてまとめる。次に、我々が行ってきた核反応モデルの高度化について紹介する。さらに、様々な入射粒子やエネルギーの条件で行ったベンチマーク[2]の結果をふまえて、今後のモデル開発に関する展望について述べる。

2.

PHITS で取り扱う物理現象と核反応モデルの役割
PHITS は、核反応が起こるスケール（数 fm）と比べて、数m から数 1,000km のマク－ 26 －

ロなスケールにおける放射線の挙動を解析する計算コードであり、放射線の振る舞いは、輸送過程と衝突過程の 2 つを交互に繰り返す形式で記述している。輸送過程では放射線による物質の電離作用（荷電粒子の場合）を考慮しており、衝突過程では放射線と物質中の原子や原子核との散乱・反応を取り扱っている。
放射線が荷電粒子の場合、輸送過程における電離作用を連続エネルギー損失仮定に基づいた阻止能計算によって評価している。ここでは、ATIMA[3]や SPAR[4]モデルを利用しており、これらの計算結果から物質中を放射線が通過した際に物質に与えるエネルギー量（発熱量）が求まる。放射線が中性子やガンマ線（光子）といった中性粒子の場合は、物質中の原子核や原子との散乱・反応によって発生する陽子や電子等の 2 次荷電粒子を通じて、間接的に物質にエネルギーを与える。なお、PHITS では、2 次荷電粒子の生成を無視し、中性粒子が “平均的に” 物質に与えるエネルギーを求めるカーマ近似も選択できる。この場合は 2 次粒子の輸送を省略するため、輸送計算全体の時間が短くなる。
ただし、これは局所的にエネルギーを付与する近似となっており、発生する 2 次荷電粒子の飛程が短い低エネルギー放射線の場合に有効である。
衝突過程は放射線が物質中の原子や原子核と衝突して起こす様々な散乱・反応を指しており、弾性散乱や非弾性散乱、吸収反応や核分裂反応といった様々な反応が含まれる。
PHITS では、次の 3 つの stage を経由することにより衝突過程を模擬している。
1.

Initial stage: 核反応の発生

2.

Dynamic stage: 前平衡過程（＋直接過程）

3.

Static stage: 複合核過程

最初の initial stage では、断面積モデルを利用することで放射線の平均自由行程を計算し、その値を参照しながら乱数によって放射線の移動距離を決定する。この計算は輸送過程の前に行い、決定した距離だけ移動した後、核反応イベントを発生させる。核反応イベントは Dynamic stage と Static stage の 2 段階で構成されている。前者では 10-22 秒程度の短時間で推移する前平衡過程や直接過程の計算を行い、後者では 10-18 秒程度の比較的長い時間で進展する複合核の崩壊過程の計算を行う。Dynamic stage における前平衡過程や直接過程は核内カスケード等の反応モデルを用いて計算しており、多数の核子が関与する複雑な反応を記述した後、ある 1 つの終状態（放出粒子の数と種類及びそれらのエネルギーと角度）を決定する。Static stage では、Dynamic stage で生成された複合核の崩壊を統計崩壊モデルにより計算する。励起エネルギーをもった残留核を複合核とみなし、ガンマ線や中性子・陽子等の粒子放出の他、核分裂による崩壊を記述する。以上の計算では乱数を用いており、放出粒子の種類やエネルギーが具体的に決まるため、モンテカル
ロシミュレーションで利用し易い event-by-event な結果を得ることができる。

－ 27 －

表 1 核反応を模擬する際の stage と役割及び各 stage で使用できるモデル
Stage

役割モデル

Initial stage

核反応の発生Pearlstein-Niita, Kurotama, Nasa(Tripathi), Shen,MWO, JENDL-4.0

Dynamic stage

前平衡過程

INCL4.6, JAM, INC-ELF, BERTINI, JQMD,（＋直接過程）JAMQMD, 光 核 反 応 , ミ ュ ー オ ン 核 反 応 ,JENDL-4.0, （DWBA との組み合わせモデル）

Static stage

複合核過程

GEM, SMM, JENDL-4.0+Event Generator Mode

表 1 に、各 stage で使用できるモデルについてまとめた。Initial stage では、初期設定で動作するデフォルトモデルとして、核子入射に関しては Pearlstein-Niita の式[5-7]を使用しており（ただし、後述する 20MeV 以下の中性子は除く）、原子核入射については Kurotama
モデル[8,9]を使用している。Kurotama モデルは核子入射にも使用でき、他にも、Nasa(Tripathi)[10-12]や Shen[13]を原子核入射の断面積モデルとして選択することができる。また、重陽子入射に特化したモデルとして MWO の式[14]も利用可能である。
次に、Dynamic stage で利用できるモデルとして、核内カスケードモデルである
INCL4.6[15] 、 JAM[16] 、 INC-ELF[17] 、 BERTINI[18] の 他 、 量 子 分 子 動 力 学 モ デ ル
JQMD[19,20]及び JAMQMD が PHITS には組み込まれている。デフォルトモデルは、3GeV未満の核子入射の場合は INCL4.6、3GeV 以上の核子入射の場合は JAM、3GeV/u 未満の原子核入射の場合は JQMD、3GeV/u 以上の原子核入射の場合は JAMQMD である。ただし、3GeV/u 未満の軽イオン（重陽子、三重陽子、3He、4He）を入射粒子とする場合は、INCL4.6 をデフォルトで使用する。加えて、陽子と重陽子入射における直接過程を考慮するための DWBA（歪曲波ボルン近似）[21]や光核反応[22]及びミューオン核反応[23]についても、簡易な設定変更により利用可能である。
最後の Static stage では、GEM[24]による統計崩壊計算のみをデフォルトで実行する。
標的が重い原子核の場合は、SMM[25]を利用することにより、3 個以上に残留核が分解する反応も考慮することができる。
以上に加えて、20MeV 以下の中性子入射については、デフォルトで評価済み核データ
ライブラリ JENDL-4.0[26]を用いた独自の計算手法を実行する。Initial stage と Dynamicstage に相当する部分は、ライブラリに含まれた断面積データを参照しながら、乱数により反応イベントを決定する。ただし、この反応イベントは、ライブラリが持っている中性子とガンマ線のスペクトル分布から放出粒子を決定するため、粒子相関が考慮されず、エネルギーや運動量の保存則が満たされない。これらの粒子の平均的な振る舞いにのみ－ 28 －

着目する場合は問題ないが、検出器応答を模擬する際など、核反応における粒子相関が重要となる場合は、Static stage において Event Generator Mode[27,28]を使用する必要がある。

3.

核反応モデルの高度化
本章では、最近我々が行った核反応モデルの高度化について、幾つかを抜粋して簡単に紹介する。
まず、断面積モデルの高度化の例として Hybrid-Kurotama モデル[9]について述べる。
元々Kurotama モデル[8] は電磁波を吸収するいわゆる完全黒体に着想を得て開発されたモデルであり、一定の半径をもつ原子核同士がぶつかると必ず反応を起こす、という描像に基づいて核子-原子核間或いは原子核-原子核間の反応断面積を計算する。Kurotama
モデルは、様々な原子核同士の反応断面積に関してその入射エネルギー依存性を再現するが、100MeV/u 以下の低エネルギー領域では有効に機能しなかった。そこで、高エネルギー側に Kurotama モデル、低エネルギー側に Nasa (Tripathi)モデル[10-12]を採用し、これらを 100-400MeV/u のエネルギー領域で滑らかに接続した Hybrid-Kurotama モデルを開発した。この Hybrid モデルには全てのエネルギー領域で有効に機能する汎用性があるため、PHITS version 3.02 より原子核入射に関するデフォルトの断面積モデルとして利用している。
量子分子動力学 JQMD[19]は原子核入射反応を記述するデフォルトモデルであるが、これを相対論的に記述することで JQMD-2.0[20]へと改良を行った。JQMD は、原子核中の核子 1 つ 1 つの広がりをガウス波束で表現し、核子同士の核力を相互作用とする運動方程式にしたがって、原子核中の核子の時間発展を記述する。JQMD は原子核同士の反応によって生成される中性子や陽子の放出スペクトルを再現可能なモデルであり、核子衝突におけるハドロン生成の素過程断面積を記述できる JAM[16]と組み合わせた JAMQMDは、3.0GeV/u 以上の高エネルギー領域でも有効に機能する。しかし、標的原子核が数個の核子を失った target-like な残留核を生成する反応については再現性が十分でなかったため、この反応が主として発生する周辺衝突（衝突径数の大きな衝突）の記述を改善した。
具体的には、周辺衝突の計算に際に発生する原子核の非物理的な崩壊が問題であったた
め、従うべき運動方程式と核子間の相互作用を相対論の形式に拡張することによりこれを解決した。相対論的に反応を記述することで、座標変換や時間経過による数値計算が適切に実行され、非物理的な反応イベントが排除可能となった。この改善により、target-like な残留核を生成する反応の再現性が大きく向上した。ただし、JQMD-2.0 を使用
すると計算時間が増加するため、デフォルトモデルとはせず、関連するパラメータの設定により動作可能とした。
DWBA（歪曲波ボルン近似）は、量子力学的に核反応を記述する手法であり、原子核－ 29 －の殻構造を考慮していわゆる直接過程の反応断面積を計算することができる。例えば、重陽子(d)入射における 6Li(d,n)7Be 反応では、7Be が基底状態や離散的な励起状態に遷移するため、放出される中性子(n)は離散的なスペクトル構造を示す。通常の核内カスケードモデルは古典的な描像に基づいているため、このような離散的なスペクトルを再現できないという問題があった。そこで、DWBA によって計算した微分断面積を PHITS に組み込み、INCL の結果と組み合わせることで、重陽子入射反応における中性子放出スペクトルの再現性を向上させることに成功した[21]。Li や Be 標的の重陽子入射反応における中性子スペクトルは加速器を用いた小型の中性子源として利用されており、その研究開発に PHITS を応用する際に、本機能は重要な役割を果たす。
ミューオン核反応[23]に関しては、仮想光子を経由して核反応を起こすという描像に基づきモデル開発を行った。その際、既に実装していた光核反応モデル[22]を利用することで、仮想光子吸収後の巨大共鳴や準重陽子共鳴の挙動を再現可能とした。また同時に、ミューオンが原子との衝突で引き起こす制動放射・電子対生成についてもモデルを導入した。更に、負ミューオンが原子軌道に捕獲された際に起こる特性 X 線の放出や、原子核に捕獲された際に生じる核反応も記述可能とすることで、ミューオンが関与する様々な反応を PHITS で模擬できるようになった。近年、ミューオンはピラミッドや火山等の巨大な構造物の内部をイメージングするために利用されており[29-31]、シミュレーションによる予測精度の向上が、これらの工学的な応用研究を大きく進展させることが期待される。
最後に、Event Generator Mode ver. 2[28]について紹介する。元々Event Generator mode[27]は、JENDL-4.0[26]等の評価済み核データライブラリを利用する際に、粒子相関を考慮した反応イベントを生成するために開発された。PHITS で核データライブラリを利用する
場合、ライブラリに含まれる中性子とガンマ線の放出スペクトルを参照しながら、乱数を用いて放出粒子のエネルギーや角度を決定する。この手法は、平均的に中性子とガンマ線の振る舞いを記述する場合は問題がないが、本来放出されるべき陽子等の荷電粒子の寄与を無視するため、粒子相関が重要となる検出器応答を模擬することが困難となる。
Event Generator mode は、上記の計算で粒子を放出したエネルギーや運動量の “残り” を残留核に与えて陽子や α 粒子を放出させることで、エネルギーと運動量保存則を満たしかつ粒子相関を考慮した反応イベントを生成する。ただし、従来の手法はライブラリを参照して放出する粒子が 1 つの場合を想定していたため、(n,2n)反応のように複数の中性子を放出する反応を再現することが難しかった。そこで、中性子を放出する段階でエネルギーや運動量の相関を考慮する手法を採用したのが Event Generator Mode ver. 2 であり、この改良を通じて、(n,2n)等の複数中性子を放出する反応が再現可能となった。

－ 30 －

4.

PHITS のベンチマーク
我々は、PHITS に組み込まれた各モデルの信頼性を検証するために、様々な条件で測定された実験値との比較を通じてベンチマークを行った[2]。原子核反応に関しては、
20MeV から 1GeV までの入射エネルギーの範囲において、中性子、陽子、軽イオン生成、target-like な放射性核種の生成反応の断面積を対象とした。加えてコンクリートや鉄材における中性子の遮へい計算も行った。その他、アクリル樹脂に 12C 線を照射した際の線量の深さ分布や電子・光子を金属材に照射した際の陽子の生成スペクトルに関しても、実験値との比較を行った。これらの多くの結果において PHITS の結果は実験値を再現し、高い信頼性をもつことが示された。
しかし、幾つかの条件では計算値と実験値に大きな違いが示され、モデルの信頼性が十分でない反応があることがわかった。例えば、単色中性子線源として利用されている7Li(p,xn)反応については、DWBA によって鋭いピーク構造は再現できるものの[32]、それ以外の成分が実験値と比べて過小評価であった。他にも、4He 入射反応における軽イオン生成反応や陽子を U 標的に入射した際の核分裂生成物についても、再現性が十分ではなかった。これらに関しては、INCL4.6 や GEM の課題として今後改良を進める予定である。

5.

まとめと今後の展望
本稿では、PHITS で取り扱う物理現象とそれを記述するために組み込まれている様々な物理モデルについて示し、近年行ってきた高度化について紹介した。我々が行ったベンチマーク結果[2]を通じてその信頼性を確認できたが、一部の反応については精度が十分でないことがわかった。加えて、PHITS を利用した応用研究[33,34]においても実験値とのずれが報告されており、今後の課題としてモデル改良に努めたい。また、核反応モデルや核データがもつ不定性が PHITS の計算結果に与える影響を定量的に見積もるために、系統的不確かさ評価機能の開発を進めている。
PHITS の利用者は年々増加しており、それに伴って得られるフィードバックも大きく、その要望に応えながら各モデルの信頼性を高めてきた。フィードバックには核データに関するものも多く、これを核データコミュニティへの要望として伝え、またその成果をPHITS 利用者へと展開してきた。差し迫っては、最近の重大な成果の 1 つである
JENDL-4.0/HE[35]を PHITS で利用できるよう整備を行いたい。また、上記の系統的不確かさ評価の際に必要となる共分散データについて、重要度の高い標的核種の選別を行い、要望として挙げる予定である。これらは一例であるが、今後も核データコミュニティの協力を得ながら核反応モデルと PHITS の開発に取り組んでいきたい。

－ 31 －参考文献
[1] T. Sato et al., “Features of Particle and Heavy Ion Transport code System (PHITS) version
3.02,” J. Nucl. Sci. Technol. 55(6), 684-690 (2018).
[2] Y. Iwamoto et al., “Benchmark study of the recent version of the PHITS code,” J. Nucl. Sci.
Technol. 54(5), 617-635 (2017).
[3] ATIMA, https://web-docs.gsi.de/~weick/atima/[4] T.M. Armstrong and K.C. Chandler, “SPAR, a FORTRAN program for computing stopping
powers and ranges for muons, charged pions, protons, and heavy ions,” ORNL Report,
ORNL-4869 (1973).
[5] S. Pearlstein, “Systematics of neutron emission spectra from high energy proton
bombardment,” Nucl. Sci. Eng. 95(2), 116-127 (1987).
[6] S. Pearlstein, “MEDIUM-ENERGY NUCLEAR DATA LIBRARIES: A CASE STUDY,NEUTRON- AND PROTON-INDUCED REACTIONS in 56Fe,” Astrophys. J. 346, 1049-1060(1989).
[7] K. Niita et al., “High-energy particle transport code NMTC/JAM,” Nucl. Instrum. Meth. B
184, 406-420 (2001).
[8] K. Iida, A. Kohama, and K. Oyamatsu, “Formula for Proton–Nucleus Reaction Cross Section
at Intermediate Energies and Its Application,” J. Phys. Soc. Japan 76, 044201 (2007).
[9] L. Sihver et al., “Current status of the ‘‘Hybrid Kurotama model’’ for total reaction crossSections,” Nucl. Instr. Meth. B 334, 34-39 (2014).
[10] R.K. Tripathi, F.A. Cucinotta, J.W. Wilson, “Accurate universal parameterization ofabsorption cross sections,” Nucl. Instr. Meth. B 117, 347-349 (1996).
[11] R.K. Tripathi, J.W. Wilson, F.A. Cucinotta, “Accurate universal parameterization ofabsorption cross sections II –neutron absorption cross sections,” Nucl. Instr. Meth. B 129,11-15 (1997).
[12] R.K. Tripathi, F.A. Cucinotta, J.W. Wilson, “Accurate universal parameterization ofabsorption cross sections III –light systems,” Nucl. Instr. Meth. B 155, 349-356 (1999).
[13] W.-Q. Shen et al., “Total reaction cross section for heavy-ion collisions and its relation to theneutron excess degree of freedom,” Nucl. Phys. A 491(1), 130-146 (198).
[14] K. Minomo, K. Washiyama, and K. Ogata, “Deuteron-nucleus total reaction cross sections upto 1 GeV,” J. Nucl. Sci. Technol. 54, 127-130 (2017).
[15] A. Boudard et al., “New potentialities of the Liège intranuclear cascade model for reactions
induced by nucleons and light charged particles,” Phys. Rev. C 87, 014606 (2013).
[16] Y. Nara et al., “Relativistic nuclear collisions at 10A GeV energies from p+Be to Au+Au withthe hadronic cascade model,” Phys. Rev. C 61, 024901 (2000).

－ 32 －[17] Y. Sawada et al., “Intranuclear cascade with emission of light fragment code implemented in
the transport code system PHITS,” Nucl. Instrum. Meth. B 291, 38-44 (2012).
[18] H. Takada et al., “An upgraded version of the nucleon meson transport code,”
NMTC/JAERI97, JAERI-Data/Code 98-005, Japan Atomic Energy Research Institute (1998).
[19] K. Niita et al., “Analysis of the (N,Xn’) reactions by quantum molecular-dynamics plus
statistical decay model,” Phys. Rev. C 52, 2620-2635 (1995).
[20] T. Ogawa et al., “Energy-dependent fragmentation cross sections of relativistic C-12,” Phys.
Rev. C 92, 024614(2015).
[21] S. Hashimoto et al., “New approach to description of (d,xn) spectra at energies below 50
MeV in Monte Carlo simulation by intra nuclear cascade code with Distorted Wave Born
Approximation,” Nucl. Instrum. Meth. B 333, 27-41 (2014).
[22] S. Noda et al., “Improvement of photonuclear reaction model below 140 MeV in the PHITS
code,” J. Nucl. Sci. Technol. 52, 57-62 (2015).
[23] S. Abe and T. Sato, “Implementation of muon interaction models in PHITS,” J. Nucl. Sci.
Technol. 54, 101-110 (2016).
[24] S. Furihata, “Statistical analysis of light fragment production from medium energy
proton-induced reactions,” Nucl. Instrum. Meth. B 171, 251-258 (2000).
[25] T. Ogawa et al., “Analysis of multi-fragmentation reactions induced by relativistic heavy ionsusing the statistical multi-fragmentation model,” Nucl. Instrum. Meth. A 723, 36-46 (2013).
[26] K. Shibata et al., “JENDL-4.0: a new library for nuclear science and engineering,” J. Nucl.
Sci. Technol. 48, 1-30 (2011).
[27] K. Niita et al., “Event generator models in the particle and heavy ion transport code system;PHITS,” J. Korean Phys. Soc. 59, 827-832 (2011).
[28] T. Ogawa et al., “Development of a reaction ejectile sampling algorithm to recover kinematic
correlations from inclusive cross-section data in Monte-Carlo particle transport simulations,”
Nucl. Instrum. Meth. A 763, 575-590 (2014).
[29] L.W. Alvarez et al., “Search for hidden chambers in the pyramids,” Science 167, 832-839(1970).
[30] K. Nagamine et al., “Method of probing inner-structure of geophysical substance with thehorizontal cosmic-ray muons and possible application to volcanic eruption prediction,” Nucl.
Instr. Meth. A 356, 585-595 (1995).
[31] D. Carbone et al. “An experiment of muon radiography at Mt Etna (Italy),” Geophys. J. Int.
196, 633-643 (2013).
[32] S. Hashimoto et al., “PHITS Simulation of Quasi-Monoenergetic Neutron Sources from7Li(p,n) Reactions,” Energy Procedia 71, 191-196 (2015).

－ 33 －[33] H. Iwamoto et al., “Impact of PHITS spallation models on the neutronics design of anaccelerator-driven system,” J. Nucl. Sci. Technol. 53(10), 1585-1594 (2016).
[34] H. Wang et al., “Spallation reaction study for the long-lived fission product 107Pd,” Prog.
Theor. Exp. Phys. 2017(2), 021D01 (2017).
[35] S. Kunieda et al., “Overview of JENDL-4.0/HE and benchmark calculation,” JAEA-Conf
2016-004, 41-46 (2016).

－ 34 －