# Fileset

[JSA VI Mtot^2 IMFP_rev2.pdf](https://mdr.nims.go.jp/filesets/7d341135-55ee-4f10-9adf-e0c715b403a1/download)

## Creator

[田沼繁夫](https://orcid.org/0000-0003-2628-9941), [篠塚寛志](https://orcid.org/0000-0001-5147-1396)

## Rights

[Creative Commons BY-NC Attribution-NonCommercial 4.0 International](https://creativecommons.org/licenses/by-nc/4.0/)

## Other metadata

[表面電子分光法における信号の減衰は如何に記述されるか？ VI.　元素における双極子行列要素平方の計算とPenn algorithmによるIMFPの計算](https://mdr.nims.go.jp/datasets/29e4ff80-3428-42bf-937e-7e22fbc6f798)

## Fulltext

名称未設定 1 表⾯電⼦分光法における信号の減衰は如何に記述されるか？ VI. 元素における双極⼦⾏列要素平⽅の計算と Penn algorithm による IMFP の計算  ⽥沼繁夫*, 篠塚寛志  物質・材料研究機構 技術開発・共⽤部⾨  〒 305-0044 つくば市並⽊ 1-1 tanuma-sh@tbd.t-com.ne.jp, tanuma.shigeo@nims.go.jp  元素固体における双極⼦⾏列要素平⽅を光学的エネルギー損失関数(ELF)から計算し，単体原⼦の結果と⽐較した。その結果，固体元素における双極⼦⾏列要素平⽅(𝑀tot2 )は対応する原⼦のそれよりも⼤幅に⼩さく，さらに原⼦番号依存性は両者で⼤きく異なっていた。また，IMFP 値から Bethe の式を⽤いて𝑀tot2 の計算を⾏い，ELF から直接計算される𝑀tot2 と⽐較した。その結果，元素固体ではスカンジウムを除いて Bethe の式から計算された𝑀tot2 は，ELFから直接計算した値よりも⼤きく，その差はおよそ 5%以内であった。さらに，実測したエネルギー損失関数と Lindhard の誘電関数から電⼦の⾮弾性平均⾃由⾏程（IMFP）を計算する Penn のアルゴリズムの詳細について，Mathematica を⽤いて具体的に計算⽅法を詳細に述べた。この⽅法の利点は，damping factor を含む Lindhard 関数を使⽤することにより，⾮常に簡単に IMFP 計算がプログラミングできることである。しかし，未確定の damping factor がインプットパラメータに⼊ること，さらに計算時間が⻑いと⾔う⽋点がある。IMFP計算を⾼速に⾏うためには damping factor を含まない⻑波⻑極限における Lindhard 関数を使⽤して FORTRAN や Julia による計算が必要となる。そこで，それらに代表される⼿続き型プログラミングにおいて，Penn のアルゴリズムで実際に数値計算を⾏う場合に必須となる項⽬について詳細に解説した。     2 How is the Signal Attenuation in Surface Electron Spectroscopy Described? VI.   The calculation of the squares of the dipole matrix elements in elemental solids and the IMFPs by the Penn algorithm  Shigeo Tanuma* and Hiroshi Shinotsuka Research Network and Facility Services Division, National Institute for Materials Science 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan * tanuma-sh@tbd.t-com.ne.jp, tanuma.shigeo@nims.go.jp  In elemental solids, the squares of the dipole matrix elements（𝑀tot2 ) were calculated from optical energy loss functions (ELFs) and compared with those of isolated atoms. It was found that these squares were significantly smaller in the elemental solid than in the corresponding atoms, and there was a notable difference in atomic number dependence between the two. Moreover, 𝑀tot2  calculated using the Bethe formula from the inelastic mean free path (IMFP) values, when compared with those directly derived from the ELF, were larger for elemental solids, except for Sc, with the difference being approximately 5%. This lecture also detailed the use of the Penn algorithm for deriving the IMFP from the measured ELF and the Lindhard dielectric function, using Mathematica for explanation. Mathematica's advantage is its facilitation of IMFP calculation using the Lindhard function, which incorporates a damping factor. However, it also has disadvantages, such as requiring an undetermined damping factor and extensive computation time. For rapid IMFP calculations, it is advisable to use the Lindhard function in the long-wavelength limit without the damping factor and to perform calculations with FORTRAN or Julia. Therefore, the lecture provided an in-depth explanation of the essential points for numerical computation using the Penn algorithm within procedural programming, as exemplified by FORTRAN.     3 １．はじめに  前回までに導いた微分⾮弾性散乱断⾯積 𝑑2𝜎/𝑑𝜔𝑑𝑞 から IMFP を実際に計算する⽅法の詳細について実践的に解説する。今回は Penn のアルゴリズを⽤いて対象物質のエネルギー損失関数（ELF）から IMFP を 10 keV以下，すなわち⾮相対論が使⽤できるエネルギー範囲の計算について解説する。  Penn の原論⽂[1]の IMFP 計算は，その簡略版である SPA (single pole approximation)の計算法に関する記述は⼤変詳しく，実践的に記載されている。しかし，今回に紹介する full Penn algorithm(FPA)についての記載は⾮常に簡略で，原論⽂を指針として計算するのは⼤変に難しい。そこで，Penn の⽅法をベースとして，篠塚らが開発した⾼速 FPA 計算法（Fortran による）をここでは紹介する。さらに，Penn の原論⽂に忠実に計算できる⽅法として Mathematica による直感的に理解できる計算⽅法についても紹介する。これらの両者の⽐較は⼤変に興味深い。講義は前回の宿題である𝑀tot2から始める。  2. 全⾮弾性散乱の双極⼦⾏列要素平⽅ 𝑀tot2  元素(孤⽴原⼦）における𝑀tot2と元素固体の𝑀tot2 を⽐較してみよう。図１に Salvat ら[2]が計算した単体元素における𝑀tot#iと我々が元素固体の ELF から下式（シリーズ V[3]の式 28）で計算した𝑀tot#を原⼦番号に対してプロットした。  𝑀tot2 =  𝑆(−1) =⎝⎜⎜⎛ 2𝑅πΩ𝑝2⎠⎟⎟⎞ ∫ Im [ −1𝜀(𝐸)] 𝑑(𝐸)∞0                              (1)  ここで，R は Rydberg 定数，Ω𝑝 = 28.816(𝜌/𝐴w)1/2 eV，Im[ −1𝜀(𝐸)] は光学的エネルギー損失関数である。また，𝑆(−1)は光学的振動⼦強度（ここでは ELF）から得られる𝑀tot2を限定的に表すものとしてここに導⼊する。ii  図１より，固体元素の𝑀tot2 は概して原⼦におけるそれと⽐べて⼩さい。また，原⼦における𝑀tot2は明確に周期構造を⽰している。すなわち原⼦の殻構造を反映し，原⼦番号の連続する希ガス元素とアルカリ⾦属の間に⼤きながギャップが⾒られる。そして，𝑀tot2の極 i 相対論的平⾯波ボルン近似により，self-consistent Dirac-Hartree-Fock-Slater potential から計算された原⼦波動関数を⽤いて計算した。 ii 𝑀tot2は式(2)が⽰すように，σ(𝐸)からも計算できるので，ELF からのみ求められる物理量として S(-1)を導⼊した。  4 ⼤値は，アルカリ⾦属より⼀つ⼤きな原⼦番号を持つアルカリ⼟類⾦属で⽰している。⼀⽅，元素固体では価電⼦数の増加に対しての変化は単体原⼦に⽐較して⼩さい。原⼦番号が増加するに従って，すなわち光吸収に関与する電⼦数が増加するので，𝑀tot2の値は増加する傾向は，単体原⼦よりも明瞭に⾒られる。  前回の講義で述べたように，𝑀tot2は式(1)からだけでなく，電⼦の⾮弾性散乱断⾯積 σ からも求めることができる。[3] すなわち，𝑀tot2を使って全⾮弾性散乱断⾯積を表現すると  σ(𝐸) = 4π𝑎02(𝐸/𝑅) [𝑀tot 2 ln (4𝑐tot 𝐸𝑅 ) + γtot (𝐸/𝑅) + 𝑂 (𝑅2𝐸2)].               (2)  𝑁𝜎𝜆 =  1より，この式から TPP-2M で使われる Modified Bethe (M. Bethe)式が導かれる。[4]  𝜆(𝐸)−1 = 𝑁  𝜎(𝐸) = 1𝐸 𝐸p2{𝛽[ln (𝛾𝐸)] − (𝐶/𝐸) + (𝐷/𝐸2)}(nm−1).       (3)  すると，𝑀tot2は下式で求めることができる。  𝑀tot2 =𝐸𝑝24𝜋𝑎02𝑅𝑁 𝛽 = 2.880𝑁𝑣𝛽                                                     (4)    ここで，式(1)における ELF から計算された𝑆(−1)と式(4)から求められる𝑀tot2 を⽐較してみる。固体元素における結果を図２に⽰す。図２(a)から両者には，当然であるが，⾮常に強い相関があることがわかる。しかし，𝑀tot2が対応する𝑆(−1)よりも少し⼤きいように⾒える。そこで，図 2(b)には𝑀tot2と𝑆(−1)の⽐を𝑆(−1) の関数として⽰す。⽐の範囲は 0.99 ~ 1.07であり，平均は 1.03 (中央値 1.03)である。この傾向は以下のように我々は考えている。𝑆(−1)は式（１）が⽰すように，𝑞 =  0である光学的 ELF から計算されるので，𝑞 >  0における ELF の寄与は無視されている。⼀⽅，式(2)-(4)から得られる𝑀tot2では IMFP データから計算されているが，q > 0 領域における⾮弾性散乱の寄与を含んでいる。したがって，𝑀tot2は 𝑞 > 0 における⾮弾性散乱の分だけ𝑆(−1)より⼤きいと考えられる。  ３．Penn algorithm による IMFP の計算  5  講義 IV[5]で述べたように，IMFP は電⼦の平均速度もしくは⾮弾性散乱断⾯積 σ から計算することができる。ここでは後者を⽤いることにする。iii前回の講義 V [3]における式(21)から始めよう。この式はHartree 単位系で表した⾮相対論的微分⾮弾性散乱断⾯積である。また，今回もHartree 単位系を使⽤する。  𝑑2𝜎𝑑𝜔𝑑𝑞 = 1𝜋𝑁𝐸 Im [ −1𝜀(𝑞, 𝜔)]1𝑞                                            (5) 定義 𝑁σλ = 1 より，  𝜆(𝐸)−1 = 1𝜋𝐸 ∬ 𝐷𝑑𝑞𝑑𝜔Im [ −1𝜀(𝑞, 𝜔)]1𝑞                               (6) 積分領域𝐷は下式で与えられる。                       𝐷 = {(𝜔, 𝑞)|0 ≤ 𝜔 ≤ (𝐸 − 𝐸F), 𝑞− ≤ 𝑞 ≤ 𝑞+}                    (7)  ここで，𝐸Fはフェルミレベル，𝑞−, 𝑞+は動⼒学的に許容される運動量移送の最⼩値と最⼤値である。また，エネルギー𝐸は⾦属における伝導帯の底を基準にしている。iv まず，𝑞−, 𝑞+を計算する。    講義 IV[5],V[3]の図 1 にあるような，散乱を考える。すると，運動量移送 𝑞 は  𝑞2 = 𝑘02 + 𝑘𝑛2 − 2𝑘0𝑘𝑛cos𝜃                                                   (8)  が得られる。このとき，q の最⼤値と最⼩値は cos 𝜃 =  −1, 1に対応することは明⽩である。 したがって 𝑞± = (k0 ± 𝑘𝑛)  となる。これをエネルギーに変換する。⼊射エネルギーを𝐸，エネルギー損失をωとすると，散乱後の電⼦ネルギー𝐸’ は (𝐸 −  𝜔 )となる。すると，𝐸と𝑞の関係からv式（８）は                  𝑞±=√2𝐸 ± √2(𝐸 − ω).                               (9)  iii Penn は電⼦の平均速度から IMFP を計算したが，ここでは Penn のアルゴリズムを,(𝑞, ω)平⾯における ELF を optical ELF から計算する⽅法と限定して話を始める。 iv 数値計算を簡単にするため。 v Hartree 単位系では 𝐸 = 𝑞2/2(= 𝑘2/2 = ħ2𝑘2/2𝑚). 𝑝 = (ℏ𝑘) = 𝑘, ここで kは波数。  6    次は𝜔のみの光学的 ELF Im[−1/𝜀(𝜔)]から 𝑞, 𝜔 の関数である完全な ELF Im[−1/𝜀(𝑞, 𝜔)]を推定する Penn のアルゴリズムについて解説しよう。vi このアルゴリズムは運動量依存性 qが統計的近似を⽤いて決定されるモデル誘電関数に基づいている。誘電関数のω依存性は光学的 ELF Im[−1/𝜀(𝜔)]によって決まる。そこで，運動量移送がゼロのときのモデル誘電関数を測定された光学誘電関数と等しくすれば(𝑞, ω)平⾯全体の ELF Im[−1/𝜀(𝑞, 𝜔)]を得ることができる。この実験的に求められる光学データを IMFP 計算に⽤いるという基本的なアイデアは，Howie と Stern，そして Powell によって開発された。[6,7]   Penn のアルゴリズムでは光学的 ELF Im[−1/𝜀(𝜔)]と式(6)におけるIm[−1/𝜀(𝜔, 𝑞)]は以下の式で関係づけられる。  Im [ −1ε(𝑞, ω)] = ∫ 𝑑ω𝑝g(ω𝑝)Im [-1𝜀𝐿(q,ω;ωp)]                                      (10)∞0  このときε𝐿は⾃由電⼦ガスに対する Lindhard のモデル複素誘電関数，ω𝑝はプラズモンエネルギー (= √4π𝑛)，n は電⼦密度，𝑔(ω𝑝)は Im[−1𝜀(𝑞=0,𝜔)] = Im[−1𝜀(𝜔)] を満たす係数で，次式であたえられる。  𝑔(𝜔) = 2𝜋𝜔Im [−1𝜀(𝜔)]                                                                     (11)  Lindhard 誘電関数は，様々な表現があるが，複素数項（damping factor)を含む場合は，次式で表される。vii  vi オリジナルよりは Tanuma-Powell-Penn の論⽂が簡潔にまとまっていて読みやすい。https://doi.org/10.48505/nims.4168 (ここからフリーで⼊⼿できる。) vii Mathematica では複素数計算が楽にできるので，この式は便利である。また，Merminの誘電関数にも damping factor を導⼊するために，この式が使われている。  7 𝜀𝐿(𝑧, 𝜇) = 1 + 𝜒2𝑧2 (12 + 18𝑧 [1 − (𝑧 − 𝜇)2]× ln [𝑧 − 𝜇 + 1𝑧 − 𝜇 − 1] + 18𝑧 [1 − (𝑧 + 𝜇)2]× ln [𝑧 + 𝜇 + 1𝑧 + 𝜇 − 1])                                                     (12)  ここで, 𝑧 = 𝑞/(2𝑘F), 𝜇 = (𝜔 + i𝛾)/(𝑞𝑘𝑉 ), 𝜒2 = 1/(𝜋𝑘F) である。kF, kVはそれぞれ Fermi 波数，Fermi 速度でありHartree 単位系では， 𝑘F = 𝑘V = (3𝜋4 )1/3𝜔𝑝2/3                                                             (13)  である。  ⼀般的に，よく使われる Lindhard 誘電関数ε𝐿(= ε1𝐿 + 𝑖ε1𝐿)は γ →  0 における次式である。（⻑波⻑極限の式） 𝜀1𝐿 = 1 + 1𝜋𝑘𝐹 𝑧2 [12 + 18𝑧 {𝐹 (𝑧 − 𝑥4𝑧) + 𝐹 (𝑧 + 𝑥4𝑧)}]                           (14𝑎)  𝜀2𝐿 = 18𝑘𝐹 𝑧3 ×⎩{⎨{⎧𝑥     for 0 < 𝑥 < 4z(1 − 𝑧)1 − (𝑧 − (𝑥/4𝑧))2     for |4𝑧(1 − 𝑧)| < 𝑥 < 4𝑧(1 + 𝑧)0     otherwise          (14𝑏)   ここで，x = ω/EF，𝐸𝐹 = 𝑘𝐹2 /2である。式(14b)における𝜀2𝐿 の領域区分の詳細は Ritchie の論⽂8に詳しいので，これを参考のために図＊３に⽰した。  基本式である式(6), (7)に Penn のアルゴリズムによる式(10)を代⼊すると以下の式が得られる。  𝜆(𝐸)−1 = 1𝜋𝐸 ∭ 𝐷𝑑𝜔𝑝𝑑𝑞𝑑𝜔 × 1𝑞 {2𝜋𝜔𝑝Im [−1𝜀(𝜔𝑝)]} Im [−1𝜀𝐿(𝑞, 𝜔, 𝜔𝑝)]             (15𝑎)  𝐷 = {(𝜔, 𝑞, 𝜔𝑝)|0 ≤ 𝜔 ≤ (𝐸 − 𝐸𝐹 ),√2𝐸 − √2(𝐸 − 𝜔) ≤ 𝑞 ≤√2𝐸 + √2(𝐸 − 𝜔), 0 < 𝜔𝑝 < ∞} (15𝑏)  これで IMFP は実測した光学的 ELF と Lindhard の ELF から計算することができる。実際の計算ではwpは⽤いる光学的 ELFデータのエネルギー損失量に相当し，その上限，下限が𝜔𝑝積分の範囲に相当する。   8 3.1 Mathematica を⽤いた FPA アルゴリズムによる IMFP の計算  Mathematica を⽤いると，今までに述べた計算式に忠実にしたがって IMFP を計算することができるので，最初に取り上げる。ここで注意すべき点は Lindhard の誘電関数としてdamping factor γを含む式(12)を使う点である。𝛾項は Penn の原論⽂には存在しない。これは𝛾 → 0 における Lindhard 誘電関数（式 14)を使うことを前提としているためである。我々は，𝛾  ≤  0.1 eVとすれば，𝛾 = 0 における計算によく⼀致することを確認している。⼀⽅，𝛾 → 0における式(14)を単純に式(15)にいれて，Mathematica で計算すると ELF が発散して計算は不可能となることを付記しておく。式(14)を⽤いる計算の詳細は FORTRAN によるIMFP 計算で紹介する。  Mathematica 13.3 を⽤いて，実際に IMFP を計算しよう。Mathematica では計算式は基本的に関数として，数式をそのまま記述すれば良い。複素関数の積分もそのまま⾏うことができ，デフフォルトで条件の悪い場合は，その点を回避して計算してくれる。実際に以下にFPA による IMFP 計算コードを⽰すが，とても簡単に計算できることがわかると思う。  １）メモリークリアー    Remove["Global`@*"]  Mathematica は基本的にインタープリター⾔語なので，最初にすべての定数，変数などをクリアするのが良い。  ２）関数の定義  Mathematica では計算は関数を基本とする。そこで，式(10) -(14)を幾つかの関数に分けて記述する。関数は遅延型で記述すると，関数が呼び出された時点で評価される。したがって，記述する順番には制限はない。  （１）Lindhard 関数と Penn の ELF の計算コード  式(12)では Lindhard 誘電関数は z, 𝜇 の関数となっているが，実質的には q (= q), ω (=w) , ω𝑝 (= wp)の関数であり，次式で記述できるviii。 epsL[q_, w_, wp_]: = 1 +4kf[wp]q2πf[z[q, wp], μ[q, w, wp]]                                   (16𝑎) ここで，   viii q: 運動量移送，ω: エネルギー損失，ω𝑝: 電⼦密度, γ ∶ plasmonの寿命  9 𝑓 [z_, u_]: = 12+ 18z((1 − (z − u)^2)Log[z − u + 1z − u − 1] + (1 − (z + u)^2)Log[z + u + 1z + u − 1])     (16𝑏) z[q_, wp_]: =q2kf[wp]                                                                                                       (16𝑐) μ[q_, w_, wp_]: =w + Iγqkf[wp]                                                                                               (16𝑑)  ここで I は虚数単位, 𝑘Fは ω𝑝の関数であり，式(13)より，  kf[wp_]: = (3 Pi 4⁄ )^(1 3⁄ )wp^(2 3⁄ ).                                                                                (17) 式(16),(17)より式(14)の Lindhard の ELF と𝑔(𝜔) = 2𝜋𝜔Im[−1𝜀(𝜔)]は，それぞれ次式となる。 elfL[q_, w_, wp_]: = Im[ −1epsL[q, w, wp]]                                                                        (18)  optical ELF の実測値を dat２: ( energy, ELF; 図４(a)) のファイルとして，spline補間を使い，最⼩ energyと最⼤ energy 間の ELF 値を求めるための関数は次式で作ることができる。  elfgs = Interpolation[dat2, InterpolationOrder → js];                                               (18𝑎)  ここで，InterpolationOrderは spline 関数の次元をあたえる。任意のエネルギーＥにおけるELF は elfgs[E]として求められる。本計算では，直線近似を採⽤し，js=1 としている。すると式(11)の𝑔(𝜔)は 𝑔[wp_]: = ( 2𝜋wp)elfgs[wp]                                                                                         (19) となる。  （２）式(15b)で与える積分範囲Dの計算 q 積分の下限値 qmin(=√2𝐸 − √2(𝐸 − 𝜔),)と上限値，qmax (=√2𝐸 + √2(𝐸 − 𝜔),)は  qmin[ek_, w_]: = Sqrt[2ek] − Sqrt[2(ek − w)]                                                            (20𝑎)  qmax[ek_, w_]: = Sqrt[2ek] + Sqrt[2(ek − w)]                                                           (20𝑏)  10  である。w積分は，上限，下限を原⼦単位で与えればそのままの式が使える。wp は計算対象となる物質の optical ELF の下限値 emin と上限値 emaxとすれば良い。 すると，式(14)の IMFP λ (= imfp3D)は, エネルギー ekの関数として，下記のように q, w, wp の三重積分として記述できる。  imfp3D[ek_]: = (π ek) NIntegrate[g[wp] elfL[q, w, wp] q⁄ , {wp, emin, emax}, {w, 0, ek − ef},⁄  {q, qmin[ek, w], qmax[ek, w]}, Method → {GlobalAdaptive, MaxErrorIncreases→ maxe}, PrecisionGoal → ms]                                                                           (21)   ここで，NIntegrate は数値積分コマンド，積分の順番は q → w → wp の順で計算される。式(16)-(21)では ef とγが定数として含まれているので，式(21)を使⽤する際は Hartree 単位で ef とγに値を⼊⼒しておくことが必要である。  （３）IMFP の数値計算の実⾏  計算に必要なパラメータ，光学的 ELF（の⼀部）の例を付録 A に⽰す。 ここで，光学的 ELF は式(19)で述べたように elfgsで, Fermi エネルギーを ef (a.u.), 減衰係数をγ (a. u. )とする。電⼦のエネルギーEx における IMFP はimfp3D[Ex]で与えられる。  実際には エネルギーはシリーズで計算するので，繰り返し計算が必要となるが，Mathematica では Do ループを使うよりは，Map 演算⼦を使うのが⼀般的である。そのイメージは関数 f[x] に{a, b, c, d, e}を作⽤させると，  In: Map[f , {a, b, c, d, e}]                                                                        (22𝑎) Out: {f[a], f[b], f[c], f[d], f[e]}                                                             (22b)  となる，計算するエネルギーを eV単位で eklistixとして与えると，Hartree 単位系で伝導帯の底から測ったエネルギーリストが得られる。  ix 例えば eklist={10.0, 20.0, ....., 2000.0} など。  11  eklistau = (eklist au2eV⁄ + ef)/. constcondition;                                  (23) ここで， constcondition = {au2eV−> 27.211396, au2nm−> 0.0529177};        (24)  するとエネルギーeklist に対応する IMFP は  mfp1 = Parallelize[Map[imfp3D, eklistau]] ∗ au2nm/. constcondition//AbsoluteTiming;   (25)  で計算することができる。Parallelizeは⾃動並列化処理にコマンド, AbsoluteTimingは，計算時間を測定するコマンドである。mfp1 には１⾏⽬に計算時間，２⾏⽬に IMFP 値が格納される。そこで，(エネルギー, IMFP)をペアとする。すなわち  mfpc = Transpose[{eklist, mfp1[[2]]}];                                             (26)  計算に使⽤した Al の ELF と計算した IMFP 値をエネルギーの関数として図４に⽰す。  ListPlot[mfpc, ScalingFunctions → {"Log", "Log"}, AxesLabel → {"Ek −Ef (eV)", "IMFP (nm)"}, LabelStyle → 12, GridLines → Automatic]                      (27)  図４より，Al の IMFP は 40 eV付近で最⼩値を取り，このエネルギー位置からエネルギーが減少,もしくは増加するに従って IMFP が増⼤する傾向があることがわかる。  このように，Mathematica を使⽤すれば Penn のアルゴリズムによる IMFP 計算は簡単に⾏うことができる。しかし， 3 eVから 1096.6 eVまでの，計 60個の計算に要した時間は，γ = 0.2 eVで 2900秒,およそ 48 分である。したがって，理想的と推察するγ ≤ 0.1 eVの計算や，10 keV 以上の⾼エネルギー領域での計算は多⼤な時間を要することになる。この問題を克服するのは⼤変に難しい。⾼ネルギー領域計算では Fortran や Julia など⾼速処理が可能な⾔語を⽤いるのが実⽤的である。  Mathematica で数値積分の精度が要求されないときは，モンテカルロ法による積分が有効である。著者はプログラムのチェックによく⽤いている。この実⾏には式(21)の Methodを以下に変更すれば良い。  12  Method → "AdaptiveMonteCarlo", PrecisionGoal → ms, MaxRecursion → 500         (28)  計算結果を図４(b)に⻘の■で⽰す。数値積分精度が悪く(被積分関数の性質が良くない)， IMFP 値のばらつきは⼤きい。しかし，エネルギー依存性は⼗分に把握することができる。このときの計算時間は 148秒であり，通常計算に⽐べて，所要時間はおよそ 1/20 となっている。  3.2 FORTRAN を⽤いた IMFP の計算（γ → 0 における Lindhard 関数を⽤いた FPA アルゴリズム）  式(12)に⽰した Lindhard 関数は虚数項を含むために，インプットパラメータとして不定（未知）の damping factor を⼊⼒する必要がある。これを避けるためにはγ → 0 における式(14)に⽰した Lindhard 関数を⽤いればよい。⼀般的にはこの形式がいろいろな応⽤に⽤いられている。ここでは，⼿続き型プログラミングの典型的な Fortran の使⽤を前提として話を進めることにする。  先に述べたように，式(14)，式(15)を⽤いて IMFP を計算すると，𝜀2𝐿 = 0 ，かつ𝜀1𝐿 = 0 の点のところで発散する。すなわちプラズモン極が存在する点である。これを避けるために，式(15a)もしくは式(12)にける ELF を𝜀2𝐿 = 0と𝜀2𝐿 ≠ 0の２つに領域に分割する。  Im [1𝜀(𝑞, 𝜔)]= Im [1𝜀(𝑞, 𝜔)]𝑝𝑙+ Im [1𝜀(𝑞, 𝜔)]𝑠𝑒                               (29)  ここで，pl は𝜀2𝐿 = 0が成り⽴つ領域で，図３に⿊の実線⽰すプラズモン分散（𝜀1𝐿 = 0）に沿った (𝑥, 𝑧)平⾯上の領域のみが有効となり，式(10)に寄与するプラズモンを表す。 se は𝜀2𝐿 ≠0 の領域における⼀電⼦励起の寄与からなる ELF を表す。 (1) Im[1𝜀(𝑞,𝜔)]𝑝𝑙パートの計算 ELF を誘電関数の実部ε1と虚数部ε2を⽤いて表記すると Im  [−1𝜀 ] = 𝜀𝜀12 + 𝜀22                                                                                (30)  ここで，”pl”の領域では𝜀2𝐿 = 0 であるので，δ関数をもちいた ELF を考えることができる。 13 すなわち，δ関数を近似式を利⽤して，𝛿(𝑥) = 1/𝜋[𝛾/(𝑥2 + 𝛾2)], (𝛾 → +0) と表現すると，容易に ELF はδ関数で記述できる。すなわち，  Im  [−1𝜀(𝑞, 𝜔, 𝜔𝑝)]= 𝜋𝛿(𝜀1),              (𝜀2 → +0)                                      (31)  ここで，𝜀1𝐿は実質的には𝑞, 𝜔, 𝜔𝑝 の関数であるので，この内の⼀つの変数で式(31)を表現することができる。xここでは Shinotsuka らの論⽂[9]に従い，𝜔𝑝の関数とする。 ここで，δ関数の公式 𝛿(𝑓 (𝑥)) = ∑  𝑖1|𝑓 ′(𝑎𝑖)|𝛿(𝑥 − 𝑎𝑖) を⽤いると，  𝛿(𝜀1(𝑞, 𝜔, 𝜔𝑝)) = ∑  𝑖1|∂𝜀1(𝑞, 𝜔, 𝜔𝑝)/ ∂𝜔𝑝|𝜔𝑝=𝜔𝑝,𝑖𝛿(𝜔𝑝 − 𝜔𝑝,𝑖)                         (32)  すると，式(10),(31),(32)より，xi  Im[−1𝜀(𝑞,𝜔)]𝑝𝑙= ∫ 𝑑ω𝑝g(ω𝑝)Im [−1𝜀𝐿(𝑞,𝜔;𝜔𝑝)]𝑝𝑙= g(ω0)∞0𝜋|∂𝜀1(𝑞,𝜔,𝜔𝑝)/∂𝜔𝑝|𝜔𝑝=𝜔0 𝜃(𝑞−(𝜔; 𝜔𝑝) − 𝑞)    (33)   ここで，𝜔0は 𝜀1𝐿(𝑞, 𝜔, 𝜔0) = 0を満たす解, 𝜃はヘヴィサイド関数（ステップ関数），𝑞−(ω; ω𝑝) =−𝑘𝐹 + √𝑘𝐹2 + 2ωである。これは図３におけるプラズモンの分散曲線が，x = 4z(𝑧 + 1)と交差する点に相当し，物理的にはプラズモンが消滅する位置である。この𝜔0値を実際に数値として求めるのはかなりやっかいである。⼀つの⽅法としては２分法により，式(14a) 𝜀1𝐿(𝑞, 𝜔, 𝜔𝑝) = 0を解くことが考えられる。また，解が存在しないときは，式(33)はゼロとなる。   式(33)にある，𝜀1𝐿(𝑞, 𝜔, 𝜔𝑝)の𝜔𝑝の微分に，Shinotsuka らは以下の式を得ている。xii  x ⽤いる変数で𝜀1𝐿の微分操作に⼤きな影響が出る。式(14b)からは𝑞もしくは𝜔 での微分が楽そうであるが，この⽅法だと ELF を(𝑞, ω)平⾯の図,すなわち Bethe Surface として⽰すことが難しくなる。 xi デルタ関数𝛿(𝑥)とヘヴィサイド関数𝜃(𝑥)にはdθ(𝑥)/dx = δ(𝑥) の関係がある。 xii この微分操作は結構厄介で，論⽂により異なっていることがあるので，計算に⽤いる微分式を選ぶ際には注意が必要である。  14   ∂ε1𝐿#𝑞, ω, ω𝑝$∂ω𝑝=13πω𝑝𝑞𝑧2 %In &𝑌− + 1𝑌− − 1& + In &𝑌+ + 1𝑌+ − 1&'                                (34)  ここで 𝑌± ≡ 𝑧 ± 14(𝑧/𝑥)−1 である。数値計算の誤差を減らすために，Shinotsuka らは 𝑧/𝑥 が⼩さいときと，⼤きいときは以下の式を⽤いることを推奨している。xiii  の場合。 ln  |𝑌− + 1𝑌− − 1| + ln  |𝑌+ + 1𝑌+ − 1| ≈ − 643𝑧𝑎2{3 + 48(1 + 𝑧2)𝑎2 + 256(3 + 𝑧2)(1 + 3𝑧2)𝑎4}       (35) ここで，𝑎 ≡ 𝑥/[𝑧(𝑧2 − 1)]である。  の場合,     (36) ここで， 𝑏 ≡ 𝑥/[𝑧(𝑧2 − 1)]である。  (2) Im[1𝜀(𝑞,𝜔)]𝑠𝑒;⼀電⼦励起部の計算 ⼀電⼦励起による寄与は図３における領域 I,II，すなわち𝜀2𝐿  ≠ 0 の領域であり，式(10),(14)から計算される。  Im  [−1𝜀(𝑞, 𝜔)]𝑠𝑒= ∫  ∞0𝑑𝜔𝑝𝑔(𝜔𝑝)Im [−1𝜀𝐿(𝑞, 𝜔; 𝜔𝑝)]𝜃(𝑞+(𝜔, 𝜔𝑝) − 𝑞)𝜃(𝑞 − 𝑞−(𝜔, 𝜔𝑝)).    (37)  ここで，𝜃はヘヴィサイド関数で，図３の領域，I,II の左右両端を表し，以下の式で与えられる。 𝑞±(ω; ω𝑝) = ± 𝑘𝐹 + √𝑘𝐹2 + 2ω                                                            (38)   xiii これらの式は式(34)をマクローリン展開して得られたものである。 z / x < 0.01z / x >100ln Y− +1Y− −1+ ln Y+ +1Y+ −1≈ ln z+1z−1#$%&'(2+ 4zb2 1+ (1+ z2 )b2 + 13(3+ z2 )(1+ 3z2 )b4)*+,-. 15 𝑢 ≡ 𝜔/𝑞𝑘𝐹 ≪ 1と𝑢 ≫ 𝑧 + 1の場合には，計算誤差を抑えるために以下の式を使⽤することを推奨する。[4] 𝑢 ≡ 𝜔/𝑞𝑘𝐹 <  0.01のとき，    ε1𝐿 = 1 + 2π𝑞𝑧 (12+ 14𝑧 [(1 − 𝑧2 − 𝑢2) ln|𝑧+1𝑧−1| + (𝑧2 − 𝑢2 − 1)2𝑢2𝑧(𝑧2−1)2]) ,   (39𝑎) ε2𝐿 = 𝑢𝑞𝑧                                                                                                     (39𝑏)            𝑢/(𝑧 + 1) > 100 のとき， 𝜀1𝐿 = 1 −𝜔𝑝2𝜔2 {1 + (𝑧2 + 35)1𝑢2}                                                     (40𝑎)) ε2𝐿 = 0                                                                                              (40𝑏)  (3) (𝑞, ω)積分計算  以上から式(10)における ELF， Im[ −1ε(𝑞,ω)] が求められたので，式(6), (7)を⽤いて IMFPが(𝑞, 𝜔) の２重積分として計算できる。繰り返しになるが，式を記述すると  𝜆(𝐸)−1 = 1𝜋𝐸 ∬ 𝐷𝑑𝑞𝑑𝜔Im [ −1𝜀(𝜔, 𝑞)]1𝑞                        (41𝑎)  𝐷 = {(𝜔, 𝑞)|0 ≤ 𝜔 ≤ (𝐸 − 𝐸F), 𝑞− ≤ 𝑞 ≤ 𝑞+}           (41𝑏)  𝑞±=√2𝐸 ± √2(𝐸 − ω)                                         (41𝑐)  ここで，E は conduction band の底から測った電⼦のエネルギー, EF はターゲット試料における実測したフェルミエネルギーである。  ここで，式(41a)を q について積分すれば，いわゆる DIIMFP, differential inelastic mean free path, が得ら れ る 。 xiv  篠 塚 ら は ω𝑝，𝑞積分を様々な ω, 𝐸 に つ て ⾏ い ，DIIMFP(𝜔, 𝐸)のテーブルを作成し，𝜔につてスプライン積分を⾏い，IMFP(E)を決定してい xiv  𝑑𝜆(𝐸)−1/𝑑𝐸 = 𝑑(𝑛σ(𝐸))/𝑑𝐸，すなわちエネルギーE の電⼦が⾮弾性散乱したときのエネルギー損失の確率分布を与える。  16 る。[4]  3.3 IMFP 値の⽐較  図５に Mathematica および Fortran よって計算された IMFP 値の⽐較を⽰す。そのエネルギー範囲は 3 eV (もしくは 10 eV)から 2 keV である。Mathematica の計算は𝛾 =  1, 0.1,0.005, 0.002 eV の計算結果を⽰す。  左側の図(A,C,E)は IMFP 値をエネルギーの関数としてプロット⽰したものである。これらより，Fortran と Mathematica で計算された IMFP 値は，およそ 10 - 30 eV のエネルギー位置で値の⼤⼩が逆転しているが，概して 30 eV 以上では両者の計算結果はよく⼀致していることが分かる。 両者の違いを詳細に知るために，図(B, D, F)では両者の相対差の絶対値, |𝜆𝑀𝑎𝑡ℎ𝑒𝑚𝑎𝑡𝑖𝑐𝑎/𝜆𝐹𝑜𝑟𝑡𝑟𝑎𝑛 − 1| ×  100 をエネルギーの関数として⽰す。これより，30 eV- 2 keV の領域で，𝛾 ≤ 0.1 eV のときは両者の差は 0.5 %以下であることがわかる。とくに，Cu, Au では 10 eV 以上の全エネルギー範囲で 0.5% 以下の差となっている。また，𝛾 =  1 eV では 30 eV - 2 keVの領域では両者の差は 5%程度であった。xv  図(A-F)から，Mathematica による IMFP 値は𝛾が⼩さくなるに従って，Fortran の計算結果に近づいており，𝛾 →  0 では両者は計算誤差範囲で⼀致するといってよいだろう。  ５．おわりに  FPA による IMFP の計算⽅法を Mathematica と Fortran を例に取り，実際にプログラミングを⾏うことを意識して説明した。特に damping 項を含む Landhard のモデル誘電関数を使⽤すると Mathematica では，⼤変に簡便にプログラミングをすることができ，その計算結果も Fortran による結果とよく⼀致する。⼀⽅，Fortran など⼿続き型プログラミングで⻑波⻑極限における Lindhard 誘電関数を⽤いると，計算は⼤変複雑となる。今回は両者を初めて⽐較した報告としても貴重であると著者らは考えている。このような⽐較は，通常の論⽂には現れることはないが，複雑な数値計算では２つの異なるプログラム⾔語で，かつ計算⽅法を変えて計算した値を⽐較することは，正確さ（誤差範囲）を確認する⼿段として⼤変重要であると考える。読者の感想はいかがであろうか？  今回は，Mathematica の使⽤を前提として話を進めたが，この⼀番の意図はその簡便さ xv 実⽤的には，計算時間を考慮すると，Mathematica では適した計算条件と考える。   17 と，柔軟さにある。例えば，ω, ω𝑝, q 積分の順番の変更（Fortran では⼤変に複雑な問題となる）や電⼦交換効果を Born-Ochur近似として組み込むなどは⼤変簡単に実⾏できる。ただし，計算の⾼速化や繰り返し計算の柔軟性は Fortran が勝るように思う。               18 参考⽂献   1 D. R. Penn, Phys. Rev. 8 35,482 (1987). 計算⽅法に限って⾔えば，その概要は IMFP Iの論⽂に詳しい。https://doi.org/10.48505/nims.4168 2 Francesc Salvat, Laia Barjuan, Pedro Andreo, Phys. Rev. A 105, 042813 (2022). https://doi.org/10.1103/PhysRevA.105.042813 3 ⽥沼繁夫, J. Surf. Anal., 30, 2 (2023). https://mdr.nims.go.jp/concern/publications/5t34sr49n?locale=en 4 H. Shinotsuka, S. Tanuma, C. J. Powell, D. R. Penn, Surf. Interface Anal. 47, 871 (2015). ibid, 47, 1132 (2015). 5 ⽥沼繁夫, Journal of Surface Analysis, 29, 72 (2022).  https://doi.org/10.1384/jsa.29.72 6  C. J. Powell, Surf. Interface Anal. 7, 263 (1985). 7 (a) A. Howie and R. M. Stern, Z. Naturforsch 27a. 382 (1972); (b) C. J. Powell, Surf. Science 44, 29 (1974). 8 R. H. Ritchie, Phys. Rev. 114, 644 (1959). 9 H. Shinotsuka, S. Tanuma, C.J. Powell, D.R. Penn, Nucl. Instrum. Meth. Phys. Res. B, 270, 75 (2012).  19 付録 A. Mathematica による ELF および計算に必要なパラメータ (例) (1) ELFファイル １⾏⽬ コメント, ２⾏⽬から energy (eV)  elf{= Im[−1/ε(𝑞 = 0, ω)]} エネルギーは昇順に並べてある。  例：Al_ELF13.el metallic aluminium from the book of Palik HOC I< 10keV  + Henke  0.04 1.530E-05 0.045 1.800E-05 0.05 2.080E-05 ……… 9320751 1.76145E-18 9760024 1.53456E-18  (2) パラメータファイル：params.dat 1          " spline の次数 " 4           "⽬標精度 デフォルト 4 " 500000    " MaxErrorIncreases  デフォルト 2000" 11.2        " Fermi energy (eV) " 0.1        " Damping coefficient (eV) "        20   Figure 1. Plots of 𝑀tot# for atoms and elemental solids as a function of atomic number. Blue solid circles show the calculated results with equation (1) from ELF. The red solid squares show the calculated results for atom by Salvat et al.[2] with the relativistic plane-wave Born approximation using self-consistent Dirac-Hartree-Fock-Slater potential.    0510152025300 20 40 60 80 100Elemental solidsAtom (by Salvat)Mtot2 Atomic numberLiNaRbSrKCaCsBaFrRaMg 21   Figure 2. (a) Plots of 𝑀$%$#  for 41 elemental solids obtained from equation (4) versus 𝑆$%$(−1) calculated from optical ELF with equation (1). (b) The ratios for 𝑀$%$# /𝑆$%$(−1) as a function of 𝑆$%$(−1).    0.90.9511.051.12 4 6 8 10Mtot2/Stot(-1)Stot(-1)(b)ScCuGe02468100 2 4 6 8 10(a) Elemental solidsMtot2 from Bethe fitStot(-1) from ELF 22                       Figure 3. Regions of the x-z plain. According to equation (14b), different 𝜀! equation is used; I: 𝜀2𝐿 = 𝑥8𝑘𝐹 𝑧3   II : 𝜀2𝐿 = 1−(𝑧−(𝑥/4𝑧))28𝑘𝐹 𝑧3   III: 𝜀2𝐿 = 0 .  Black line indicate the volume plasmon dispersion.      02468100 0.5 1 1.5 2 2.5 3x=4z(1-z)x=4z(z+1)x=4z(z-1)x=ω/EFz = q/2kFIIIIIIIII 23      Figure 4  Energy loss function and IMFPs for Al calculated with Full Penn algorithm using γ = 0.2 eV in the energy range between 3 eV and 1000 eV. (a) Energy loss function for Al. (b) IMFPs for Al calculated from equation (21) using γ = 0. 2 eV for "Normal" in caption and IMFPs calculated from equations (21) and (28) for "MC" in caption.    10-510-410-310-210-110010110210-1 100 101 102 103AlEnergy loss functionEnergy loss (eV)(a)10-1100101102100 101 102 103NormalMCInelastic mean free path (nm)Electron energy above Fermi level (eV)(b) 24  Figure.5 Comparison of IMFPs for Al, Cu, and Au in the energy range from 3 eV (or 10 eV) to 2000 eV, calculated from optical ELFs using Mathematica and Fortran.