Pythonで学ぶ金属元素蛍光X線のラインプロファイル : 特に Fe Kα・Kβ を"ちゃんと"理解する

蛍光X線プロファイル

蛍光X線の中心エネルギーや強度比、線の「タイプ」や自然幅は、最近では xraylib のようなライブラリを使えば簡単に取得できるようになりました。

しかし、実際にスペクトルをフィットして装置分解能を評価しようとすると、

• 「Kα1, Kα2 を 2 本 Voigt で足しただけ」では足りない
• 微細構造や衛星線（shake-up / shake-off など）を含めて、7〜8 本の Voigt を足さないと形が合わない

という世界になります。

この記事では、

1. まず Fe の Kα・Kβ を題材にして、物理的に何が起きているのかを丁寧に説明し、
2. そのうえで、Python で複数 Voigt + テイルを足したラインプロファイルを生成するコードを紹介します。

実測スペクトルを実際にフィットしたい方は、姉妹記事
👉 pythonで複数のvoigt関数にレスポンスを入れてフィットする方法
もあわせて参照してください。

最近は，蛍光X線の中心エネルギー，強度比，タイプ，自然幅は，xraylibを使えば，簡単に取得できる．しかし，実際上は，量子多体系の効果でラインプロファイルが広がる効果を加味する必要があり，例えば，MnのKalpha線はKalpha1,2の強度比が1:2で終わりのはずが，7つや8つのVoigt関数をいれてフィットしないと，正しく装置の分解能を評価できません。以下では、いくつか代表的なラインプロファイルのパラメータと簡単な生成方法をpythonでまとめておきました。

その前に、ラインの呼び方には２通りあることを知っておいた方がよいです。 Siegbahn と IUPAC という読み方がある、Siegbahn notation など。簡単に考えると、IUPAC は、KL1,KL2,KL3のように数字が大きくなるほど、エネルギーが高くなる。例えば、L1, L2, and L3 が、2S1/2, 2P1/2, and 2P3/2 に対応しているためである。2S1/2の読み方は、蛍光X線スペクトルの読み方について など参照ください。主量子数nを冒頭の数字、軌道角運動量をアルファベット、全角運動量(spin+orbit)を３つめの数字で表記したものでです。

1. Kα・Kβ ってそもそも何？

1.1 シェルと軌道の復習

原子の束縛電子は、量子数でラベリングされます：

• 主量子数： $n = 1, 2, 3, \dots$

• 軌道角運動量： $l = 0, 1, 2, 3, \dots$

  • $l = 0 \to s$軌道, $l = 1 \to p$軌道, $l = 2 \to d$軌道, ...

• 全角運動量（スピンを含む）：$j$

物理的には、

• $n=1$ のシェル（K殻）
• $n=2$ のシェル（L殻）
• $n=3$ のシェル（M殻）

という言い方をします。

たとえば Fe (Z=26) だと、基底状態では

• 1s 軌道が 2 個の電子で満たされている
• 2p, 3p などにも電子が詰まっている

などと思っておけばとりあえずは十分です。

1.2 Kα と Kβ の定義

何かしらの高エネルギー光子や電子が当たって、K殻の 1s 電子が叩き出されて空孔（ホール）ができると、その原子は不安定になります。

このホールを埋めるために、上の軌道から電子が落ちてきて、その差分エネルギーを X線として放出します。これが蛍光X線です。

• Kα線

  • K 殻（1s）にできたホールを、L 殻（n=2）からの電子が埋める遷移
  • 大雑把にいうと「2p → 1s」の遷移

• Kβ線

  • K 殻（1s）にできたホールを、M 殻（n=3）からの電子が埋める遷移
  • 大雑把にいうと「3p → 1s」や「3d → 1s」などの遷移

Kβ線は、Kβ1とKβ3のことを暗黙に意味してることが多いです。Kα線:Kβ線=10:1の強度比であることは覚えておきましょう。これも厳密には、Kα線:Kβ1,3線:Kβ5線=10:1:0.1の強度比となることも覚えておくと良いでしょう。

K 殻への遷移なので、どちらも「K 何か」という名前がつきますが、

• 「どの殻から落ちてきたか」によって α, β, γ... などに分かれます。

2. Siegbahn 表記と IUPAC 表記

蛍光X線の名前の付け方には、大きく分けて 2 通りあります。

• Siegbahn 表記

  • Fe Kα$_1$, Kα$_2$, Kβ$_1$ ... のような馴染みのある記号

• IUPAC 表記

  • K-L$_3$, K-L$_2$, K-M$_3$ ... のように書く

2.1 IUPAC 表記の基本

IUPAC では、

• K, L, M ... が 最終状態（落ちていった先）のシェル

• L$_1$, L$_2$, L$_3$ がそれぞれ

  • L1 → 2s_{1/2}
  • L2 → 2p_{1/2}
  • L3 → 2p_{3/2}

に対応します。

たとえば Fe の代表的な遷移は：

• Fe Kα(*1) ≒ K-L3（2p$_{3/2} \to$ 1s）
• Fe Kα(*2) ≒ K-L2（2p$_{1/2} \to$ 1s）
• Fe Kβ(*1) ≒ K-M3（3p$_{3/2} \to$ 1s）

という対応になります。

IUPAC 表記まできちんと押さえておくと、

• xraylib などのデータベース
• 原子分光の成書・論文

と話が通じやすくなるのでおすすめです。

量子数の読み方$2p(_{3/2}$ など）は
「蛍光X線スペクトルの読み方について」などを参照ください。

詳細やxraylibの使い方は下記を参考にしてください。

3. Fe Kα・Kβ はなぜ 2本ではなく「束」なのか？

教科書的には

• Fe Kα：6.40 keV 付近
• Fe Kβ：7.06 keV 付近

と 1 本ずつ書かれていることが多いですが、高分解能スペクトルで見ると、Fe Kα だけで何本もピークが並んでいるように見えます。

理由は大きく分けて：

1. スピン軌道相互作用による細かい準位分裂

   • 2p_{3/2}, 2p_{1/2} など

2. 多電子効果（マルチプレット構造）

   • ホールを 2 個以上含む状態、shake-up / shake-off など

3. 外殻電子配置の違い

   • Fe の価数や化学状態によって、わずかにエネルギーがシフト

など、「電子が 1 個と空孔 1 個だけ」では済まない量子多体効果が効いてきます。

3.1 このコードで使っている Fe Kα の成分

記事中の Python コードでは、Fe Kα を次のような配列で与えています。

name="Fe Kalpha"
energy = np.array((6404.148, 6403.295, 6400.653, 6402.077, 6391.190, 6389.106, 6390.275))
lgamma = np.array((1.613, 1.965, 4.833, 2.803, 2.487, 2.339, 4.433))
amp = np.array([697, 376, 88, 136, 339, 60, 102])

ここで：

• energy：各成分の中心エネルギー [eV]
• lgamma：各成分の自然幅（ローレンツ幅）の HWHM [eV]
• amp：各成分の相対強度（後で規格化）

となっています。

物理的には、

• 上のほうの数本が 主線（Kα$_1$, Kα$_2)$の微細構造
• 下のほうの数本が 衛星線（多電子効果を反映した遷移）

に相当すると考えればイメージしやすいと思います。

重要なのは：

「Fe Kα」と言った瞬間に、実際には 1 本ではなく複数本の Voigt の足し合わせが必要になる

という点です。
装置分解能を真面目に評価したいときは、ここをサボると 装置側の分解能と線プロファイル由来の広がりを取り違えてしまう危険があります。

4. Voigt プロファイルと装置応答

4.1 基本：自然幅＋ガウシアン応答

1 本の遷移線を考えると、

• 自然幅：原子核・電子状態の有限寿命 → ローレンツ型の広がり
• 装置分解能：検出器のエネルギー分解能 → ほぼガウシアンとみなせる

ことが多いので、両者の畳み込みである Voigt 関数で表現するのが標準的です。

• ローレンツ成分：線の「裾」が長く伸びる
• ガウシアン成分：装置分解能が悪いと線全体がボケる

このコードでは、

• lgamma がローレンツ成分（自然幅）の HWHM
• gw がガウシアン成分（装置分解能）の $\sigma$

として Voigt を構成しています。

4.2 低エネルギーテイル

実際の検出器（特に半導体検出器や一部のマイクロカロリメータ）では、

• 電荷収集の不完全さ
• L 電子 / M 電子のオージェ放出
• クロストークやデッドレイヤーなど装置由来

の影響で、低エネルギー側に引きずられたテイルが出ることがあります。

このコードでは、smear() 関数の中で、

• テイルの割合：P_tailfrac
• テイルの「長さ」：P_tailtau

というパラメータで、指数関数的なテイルを FFT で畳み込む形でモデル化しています。

5. 全体コード：複数 Voigt + テイルによるライン形状生成

ここからが実際の Python コードです。
Fe の Kα・Kβ だけでなく、Ti, V, Cr, Mn, Co, Ni, Cu, As なども一気に描画します。

Google Colab からも動作チェックできます。

5.1 コードの概要

Lineprofile クラスが，ラインのエネルギー，強度，自然幅，装置の分解能を加味して，Voigt関数を生成する．mymodel クラスの中では，それをrawfuncという関数にして，必要なら低エネルギー側のテイルの計算も行えるようにしている．

5.2 Fe Ka と Fe Kb をプロットするコードとスペクトル

import matplotlib.pyplot as plt
import numpy as np
import scipy.special
plt.rcParams['font.family'] = 'serif'

# ---- ここは元コードと同じ ----
def voigt(xval,params):
    center,lw,gw = params
    z = (xval - center + 1j*lw)/(gw * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    model_y = (w.real)/(gw * np.sqrt(2.0*np.pi))
    return model_y

def smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly = False):
    if P_tailfrac <= 1e-5:
        return rawfunc(x)

    dx = x[1] - x[0]
    freq = np.fft.rfftfreq(len(x), d=dx)
    rawspectrum = rawfunc(x)
    ft = np.fft.rfft(rawspectrum)
    if tailonly:
        ft *= P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 0)
    else:
        ft += ft * P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 1)

    smoothspectrum = np.fft.irfft(ft, n=len(x))
    if not tailonly:
        smoothspectrum[smoothspectrum < 0] = 0
    return smoothspectrum

def Lineprofile(xval, params, energy, lgamma, amp, consts=None):
    norm,gw,gain,bkg1,bkg2 = params

    if consts is None:
        consts = np.ones(len(energy))

    prob = (amp * lgamma) / np.sum(amp * lgamma)

    model_y = 0.0
    for ene, lg, pr, con in zip(energy, lgamma, prob, consts):
        voi = voigt(xval, [ene*gain, lg*0.5, gw])
        model_y += norm * con * pr * voi

    background = bkg1 * np.ones_like(xval) + (xval - np.mean(xval)) * bkg2
    model_y += background
    return model_y

def mymodel(x, params, energy, lgamma, amp,
            consts=None, tailonly=False):
    norm,gw,gain,P_tailfrac,P_tailtau,bkg1,bkg2 = params

    def rawfunc(xx):
        return Lineprofile(xx, [norm,gw,gain,bkg1,bkg2],
                           energy, lgamma, amp, consts=consts)

    return smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly=tailonly)

# ---- ここから Fe 専用設定 ----

# 共通パラメータ（ガウシアン幅など）
gfwhm = 2.0              # ガウシアンのFWHM [eV]
gw    = gfwhm / 2.35
norm  = 5e4
gain  = 1.0
bkg1  = 0.0
bkg2  = 0.0
P_tailfrac = 0.0         # とりあえずテイルなし
P_tailtau  = 10.0
init_params = [norm, gw, gain, P_tailfrac, P_tailtau, bkg1, bkg2]

nbin   = 2000

# --- Fe Kα ---
energy_ka = np.array((6404.148, 6403.295, 6400.653, 6402.077,
                      6391.190, 6389.106, 6390.275))
lgamma_ka = np.array((1.613, 1.965, 4.833, 2.803,
                      2.487, 2.339, 4.433))
amp_ka    = np.array([697, 376, 88, 136, 339, 60, 102])

emin_ka = 6380.0
emax_ka = 6420.0
x_ka = np.linspace(emin_ka, emax_ka, nbin)
y_ka = mymodel(x_ka, init_params, energy_ka, lgamma_ka, amp_ka)

# --- Fe Kβ ---
energy_kb = np.array((7046.90, 7057.21, 7058.36, 7054.75))
lgamma_kb = np.array((14.17, 3.12, 1.97, 6.38))
amp_kb    = np.array([107, 448, 615, 141])

emin_kb = 7030.0
emax_kb = 7075.0
x_kb = np.linspace(emin_kb, emax_kb, nbin)
y_kb = mymodel(x_kb, init_params, energy_kb, lgamma_kb, amp_kb)

# --- プロット ---
fig, axes = plt.subplots(2, 1, figsize=(8, 8), sharex=False)

# Fe Kα
axes[0].plot(x_ka, y_ka, '-')
axes[0].set_title("Fe Kα complex (zoomed)")
axes[0].set_xlabel("Energy [eV]")
axes[0].set_ylabel("Arbitrary units")
axes[0].grid(alpha=0.3)

# Fe Kβ
axes[1].plot(x_kb, y_kb, '-')
axes[1].set_title("Fe Kβ complex (zoomed)")
axes[1].set_xlabel("Energy [eV]")
axes[1].set_ylabel("Arbitrary units")
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig("FaKa_Kb.png")
plt.show()

5.3 この図を見るときのポイント

• Kα パネル

  • 6.404 keV 付近に一番高い山 → Kα$_1$ 主成分
  • そのすぐ左（6.40 keV 付近）に少し低い山 → Kα$_2$や微細構造
  • さらに左、6.39 keV 付近に小さな山 → サテライト成分
    → もし装置分解能が良ければ、「主線＋肩」構造として実際のスペクトルにも現れる

• Kβ パネル

  • 7.057–7.058 keV あたりがメインの塊
  • 7.046 keV 付近や 7.054 keV 付近に、幅の広い成分が混ざっていて
    → 全体として少し非対称な形になる
    → Kβ 側のプロファイルも「1本線」ではなく、複合線になっている

5.4 実測スペクトルと比べるときの指針

Fe の Kα/Kβ を Resolve や TESなどで測るとき、このモデルと見比べて注目すると良いポイントは：

1. ピークの中心エネルギー

   • エネルギー校正のチェック
   • 化学シフトがあれば、Kα/Kβ の位置が少しずつずれる

2. Kα 複合体の「非対称性」

   • 左側（低エネルギー側）のサテライトやテイルの入り方
   • 検出器由来のテイルと、線リスト由来のサテライトをどう切り分けるか

3. Kα vs Kβ の強度比

   • Kβ/Kα 比が異常に高い／低い →
     実験条件や自己吸収、モデルの取り込み漏れなどを疑うポイント (宇宙物理特有のズレかもしれない。)

4. 装置分解能の評価

   • gfwhm を変えながらこのモデルをフィットして、
     **「この gfwhm だと実測と一番よく合う」**という値を探すと、
     分解能の実効的な評価になる。

5. 電離状態によって Kα/Kβ のエネルギーが変わる

（宇宙 X 線天文学で特にここが重要）

宇宙空間では、鉄は固体ではなく、高温の気体（プラズマ）状態が多く、そこでは電子が複数個失われた状態（電離状態） になりえます。(宇宙ダストの中の鉄は固体でXAFSなどの化学的な効果が期待できます。)

電離が進むと：

• 内殻電子の束縛エネルギーが変わる
• 結果として Kα / Kβ の遷移エネルギーが少しずつ変化する
• → 実測した Kα/Kβ の中心エネルギーから、イオン化度を推定できる

たとえば Fe II ～ Fe IX のように電子が数個失われただけの場合でも
Kα/Kβ の中心位置は ~1 eV 程度の系統的変化を示します。

XSTAR でも、電離パラメータ ξ の違いにより線位置が変化する様子が示されている：

Modeling of iron K lines: Radiative and Auger decay data for Formula –Formula
P. Palmeri, C. Mendoza, T. R. Kallman, M. A. Bautista and M. Meléndez
A&A, 410 1 (2003) 359-364
DOI: https://doi.org/10.1051/0004-6361:20031262

このPalmeri et al. (2003) 論文においても、電離段階に応じて Kα と Kβ の中心波長（∼エネルギー）が変化することが示されており、モデル化に極めて重要とされています：

📌 Fe II～IX の Kα/Kβ 中心波長の電離依存：

• Palmeri et al. (2003) 論文のFig.1 (p.361) に N（束縛電子数）とともに Kα/Kβ 計算値が変化
• 電子数が減る（電離が進む）と Kα がわずかにブルーシフトする傾向
• Kβ は逆にわずかにレッドシフト

📌 Kβ/Kα 強度比も電離度に依存

• Palmeri et al. (2003) 論文のFig.2 (p.361) のように、Kβ/Kα 比は電子数（電離度）によって連続的に変化

つまり、観測された Fe Kα/Kβ の中心位置や強度比は、鉄の電離段階（イオン化度）の診断指標になります。

🔍 なぜ電子が取れると蛍光 X線のエネルギーが変わるのか？（物理の本質）

◆ ① 内殻電子は「外殻電子の遮蔽（screening）」を受けている

鉄の Kα 遷移は

• 1s → 2p の電子遷移

鉄の Kβ は

• 1s → 3p の遷移

本質的には 内殻電子が陽子 26 個から受ける引力 − 外殻電子の遮蔽効果 で決まる。

外殻電子が存在することで、内殻電子が感じる実効核電荷 $Z_{\mathrm{eff}}$ が弱くなる。

◆ ② 電離が進むと遮蔽が弱くなる

鉄が電離して

• 3d, 4s, 3p などの外殻電子が失われる
  ↓
  遮蔽が弱くなる
  ↓
  内殻電子が感じる核電場が強くなる
  ↓
  内殻電子の束縛エネルギーが増す → 遷移エネルギー（蛍光エネルギー）が上昇（ブルーシフト）

◆ ③ Kα と Kβ は異なる軌道が関与

Kα：2p → 1s
Kβ：3p → 1s

電離が進むと

• 外殻の 3p, 3d 電子数の変動が
• Kβ → Kα より強く影響する場合もある

そのため、

• Kα：主として 2p 電子の遮蔽変化
• Kβ：より外側の 3p 電子の遮蔽変化

によって、電離度ごとに Kα と Kβ の変化率が異なる。

論文（Palmeri et al. 2003）はまさにこれを HFR 計算で示し、
N(束縛電子数) の関数として Kα/Kβ の中心波長や強度比が系統的に変化することを報告。

◆ ④ 量子力学的には、全電子系の固有状態の再構成

電子が減ると

• 多電子軌道が再配置され
• 有効ポテンシャルが変わり
• スピン軌道相互作用も変化し
• 遷移エネルギーが変動

つまり「線の位置が変わる」のは単なる遮蔽の変化以上に、全電子系の再計算値として当然の結果と言える。

◆ ⑤ 宇宙での応用例

観測された Kα の中心値が

• 6.400 keV → 中性近い (Fe I–IX)
• 6.45 keV → 部分電離 (Fe XVI–XXIII)
• 6.7 keV → He-like (Fe XXV)
• 6.97 keV → H-like (Fe XXVI)

というように、線中心位置だけで Fe の電離度が推定できる。

◆ ⑥ 結論

電離によって Kα/Kβ のエネルギーが変わる理由は：

• 外殻電子の喪失 → 遮蔽低下 → 実効核電荷上昇
• 束縛エネルギー増大 → 蛍光遷移エネルギー上昇
• 多電子効果・スピン軌道相互作用の変化
• イオン化度ごとに異なる安定構成・軌道構造

という量子力学的な必然である。

そのため Kα/Kβ は電離診断線として最も強力なツールになる。

全ラインを一気にプロット

plot_manylines_wletail_qiita.py

#!/usr/bin/env python
__version__= '1.0'

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import numpy as np
import scipy.special

def mymodel(x,params, consts=[], tailonly = False):
    norm,gw,gain,P_tailfrac,P_tailtau,bkg1,bkg2 = params    
    # norm : nomarlizaion
    # gw : sigma of gaussian
    # gain : gain of the spectrum 
    # P_tailfrac : fraction of tail 
    # P_tailtau : width of the low energy tail
    # bkg1 : constant of background
    # bkg2 : linearity of background    
    initparams = [norm,gw,gain,bkg1,bkg2]
    def rawfunc(x): # local function, updated when mymodel is called 
        return Lineprofile(x,initparams,consts=consts)               
    model_y = smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly=tailonly)
    return model_y

def Lineprofile(xval,params,consts=[]):
    norm,gw,gain,bkg1,bkg2 = params
    # norm : normalization 
    # gw : sigma of the gaussian 
    # gain : if gain changes
    # consttant facter if needed 
    prob = (amp * lgamma) / np.sum(amp * lgamma) # probabilites for each lines. 
    model_y = 0 
    if len(consts) == 0:
        consts = np.ones(len(energy))
    else:
        consts = consts

    for i, (ene,lg,pr,con) in enumerate(zip(energy,lgamma,prob,consts)):
        voi = voigt(xval,[ene*gain,lg*0.5,gw])
        model_y += norm * con * pr * voi

    background = bkg1 * np.ones(len(xval)) + (xval - np.mean(xval)) * bkg2
    model_y = model_y + background
    # print "bkg1,bkg2 = ", bkg1,bkg2, background
    return model_y

def voigt(xval,params):
    center,lw,gw = params
    # center : center of Lorentzian line
    # lw : HWFM of Lorentzian (half-width at half-maximum (HWHM))
    # gw : sigma of the gaussian 
    z = (xval - center + 1j*lw)/(gw * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    model_y = (w.real)/(gw * np.sqrt(2.0*np.pi))
    return model_y

def smear(rawfunc, x, P_tailfrac, P_tailtau, tailonly = False):
    if P_tailfrac <= 1e-5:
        return rawfunc(x)

    dx = x[1] - x[0]
    freq = np.fft.rfftfreq(len(x), d=dx)
    rawspectrum = rawfunc(x)
    ft = np.fft.rfft(rawspectrum)
    if tailonly:
        ft *= P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 0)
    else:
        ft += ft * P_tailfrac * (1.0 / (1 - 2j * np.pi * freq * P_tailtau) - 1)

    smoothspectrum = np.fft.irfft(ft, n=len(x))
    if tailonly:
        pass
    else:
        smoothspectrum[smoothspectrum < 0] = 0
    return smoothspectrum    

class aline:
    def __init__(self,x,y,name):
        self.x = x
        self.y = y
        self.name = name

# global variables
gfwhm = 2
gw = gfwhm / 2.35
norm = 500000.0
gain = 1.0
bkg1 = 1.0
bkg2 = 0.0
P_tailfrac = 1e-6 # no tail 
P_tailtau = 10
nbin=1000
ewidth=500
init_params=[norm,gw,gain,P_tailfrac,P_tailtau,bkg1,bkg2]

linelist = []


#################################################################
name="Ti Kalpha"
energy = np.array((4510.918, 4509.954, 4507.763, 4514.002, 4504.910, 4503.088))
lgamma = np.array((1.37, 2.22, 3.75, 1.70, 1.88, 4.49))
amp = np.array((4549, 626, 236, 143, 2034, 54))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Ti KBeta"
energy = np.array((25.37, 30.096, 31.967, 35.59)) + 4900
lgamma = np.array((16.3, 4.25, 0.42, 0.47))
amp = np.array((199, 455, 326, 19.2))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="V Kalpha"
energy = np.array((4952.237, 4950.656, 4948.266, 4955.269, 4944.672, 4943.014))
lgamma = np.array((1.45, 2.00, 1.81, 1.76, 2.94, 3.09))
amp = np.array((25832, 5410, 1536, 956, 12971, 603))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="V KBeta"
energy = np.array((18.19, 24.50, 26.992)) + 5400
lgamma = np.array((18.86, 5.48, 2.499))
amp = np.array((258, 236, 507))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Cr Kalpha"
energy = 5400 + np.array([14.874, 14.099, 12.745, 10.583, 18.304, 5.551, 3.986])
lgamma = np.array([1.457, 1.760, 3.138, 5.149, 1.988, 2.224, 4.4740])
amp = np.array([882, 237, 85, 45, 15, 386, 36])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Cr KBeta"
energy = 5900 + np.array((47.00, 35.31, 46.24, 42.04, 44.93))
lgamma = np.array([1.70, 15.98, 1.90, 6.69, 3.37])
amp = np.array([670, 55, 337, 82, 151])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Mn Kalpha"
energy = 5800 + np.array((98.853, 97.867, 94.829, 96.532, 99.417, 102.712, 87.743, 86.495))
lgamma = np.array([1.715, 2.043, 4.499, 2.663, 0.969, 1.553, 2.361, 4.216])
amp = np.array([790, 264, 68, 96, 71, 10, 372, 100])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Mn KBeta"
energy = 6400 + np.array((90.89, 86.31, 77.73, 90.06, 88.83))
lgamma = np.array((1.83, 9.40, 13.22, 1.81, 2.81))
amp = np.array([608, 109, 77, 397, 176])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Fe Kalpha"
energy = np.array((6404.148, 6403.295, 6400.653, 6402.077, 6391.190, 6389.106, 6390.275))
lgamma = np.array((1.613, 1.965, 4.833, 2.803, 2.487, 2.339, 4.433))
amp = np.array([697, 376, 88, 136, 339, 60, 102])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Fe Kbeta"
energy = np.array((7046.90, 7057.21, 7058.36, 7054.75))
lgamma = np.array((14.17, 3.12, 1.97, 6.38))
amp = np.array([107, 448, 615, 141])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Co Kalpha"
energy = np.array((6930.425, 6929.388, 6927.676, 6930.941, 6915.713, 6914.659, 6913.078))
lgamma = np.array((1.795, 2.695, 4.555, 0.808, 2.406, 2.773, 4.463))
amp = np.array((809, 205, 107, 41, 314, 131, 43))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Co Kbeta"
energy = np.array((7649.60, 7647.83, 7639.87, 7645.49, 7636.21, 7654.13))
lgamma = np.array((3.05, 3.58, 9.78, 4.89, 13.59, 3.79))
amp = np.array((798, 286, 85, 114, 33, 35))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Ni Kalpha"
energy = np.array((7478.281, 7476.529, 7461.131, 7459.874, 7458.029))
lgamma = np.array((2.013, 4.711, 2.674, 3.039, 4.476))
amp = np.array((909, 136, 351, 79, 24))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Ni Kbeta"
energy = np.array((8265.01, 8263.01, 8256.67, 8268.70))
lgamma = np.array((3.76, 4.34, 13.70, 5.18))
amp = np.array((722, 358, 89, 104))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="Cu Kalpha"
energy = np.array([8047.8372, 8045.3672, 8027.9935, 8026.5041])
lgamma = np.array([2.285, 3.358, 2.667, 3.571])
amp = np.array([957, 90, 334, 111])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="Cu KBeta"
energy = np.array([8905.532, 8903.109, 8908.462, 8897.387, 8911.39])
lgamma = np.array([3.52, 3.52, 3.55, 8.08, 5.31])
amp = np.array([757, 388, 171, 68, 55])
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


#################################################################
name="As Kalpha"
energy = np.array((10543.2674,10507.50))
lgamma = np.array((3.08, 3.17))
amp = np.array((1.00, 0.51))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))

#################################################################
name="As KBeta"
energy = np.array((11725.73,11719.86))
lgamma = np.array((2.09+2.25, 2.09+2.25))
amp = np.array((0.13, 0.06))
xmin=np.mean(energy) - ewidth
xmax=np.mean(energy) + ewidth
x = np.linspace(xmin,xmax,nbin)
model_y = mymodel(x,init_params)
linelist.append(aline(x,model_y,name))


plt.figure(figsize=(10,8))
plt.title("Line Profiles")

for oneline in linelist:

    plt.xlabel("Energy (eV)")
    plt.plot(oneline.x, oneline.y, '-', label = oneline.name)
    plt.legend(numpoints=1, frameon=False, loc="best")
    plt.grid(linestyle='dotted',alpha=0.5)

plt.savefig("linelist.png")
plt.show()
 

使い方

単純に実行すると，このような図が生成される．

縦軸をログの方が見やすいかもしれない。

plt.figure(figsize=(10,8))
plt.title("Line Profiles")
plt.yscale("log")
for oneline in linelist:

    plt.xlabel("Energy (eV)")
    plt.plot(oneline.x, oneline.y, '-', label = oneline.name)
    plt.legend(numpoints=1, frameon=False, loc="best")
    plt.grid(linestyle='dotted',alpha=0.5)

plt.savefig("linelist_log.png")
plt.show()

もし，装置の低エネルギー側のテイルを入れたい場合は，P_tailfrac = 1e-6 # no tail の部分を大きな数字にすればよい．

energy，lgamma，amp，という global 変数を用意して，それを書き換えるという邪道な書き方をしているが，正しくは，データベースとなる情報はそもそもコードと一体管理しないほうが良いし，ラインごとにちゃんとしたオブジェクトを生成した方がよい．

まとめ

ここまで読むと、「論文や xraylib や NISTなどの数値をそのまま使うだけではなくて、そもそも、どうやってその数値が計算されているのかを知りたい」という欲が出てくると思います。

最後に、「Fe Kα/Kβ まわりの理論計算を自分でできるようになりたい人」向けに少しだけ紹介です。

多電子原子・原子分光をもう少し真面目にやりたくなった人向け

「教科書にいきなり行くのはちょっと重い…」という場合、まずは日本語の総説から雰囲気をつかむのがおすすめです。

たとえば原子衝突学会誌のシリーズ

「多電子原子の構造とダイナミックス ― 独立粒子モデルの来し方行く末 ―」（小池 文博 先生）

は、“独立粒子モデルってそもそも何なのか？” から、“実際に電子状態計算をするときのコード選び” までを丁寧につないでくれるシリーズです。

ここでは、そのうち 2 本を入口として紹介しておきます。

(1) 第1回「擬制としての 1 電子軌道概念」([J-STAGE][1])

• 多電子原子の歴史的な理解の流れ

  • 水素原子 → 多電子原子 → 自動電離やオージェ過程など「多電子ならではの現象」へ、という発展の流れを俯瞰。

• 平均場近似と 1 電子軌道の考え方

  • Hartree–Fock / Dirac–Fock 型の「平均場＋1電子軌道」の枠組みが、
    原子・分子・固体の “殻模型／軌道” の直感を与えてきたことを整理。

• 1 電子軌道はあくまで「擬制（フィクション）」である、という視点

  • 独立粒子モデルは、厳密解ではなく「展開基底のひとつ」にすぎない。
  • 強く相関した系に対して、どこまで 1 電子軌道で押し通せるか、
    どこから破綻して別の記述（超球座標モデルなど）が必要になるか、という問題意識が提示される。

• メッセージとしては：

  「便利だからみんな使っている 1 電子軌道だけれど、何がうまく説明できて、どこから先は危ないかを意識しよう」

Fe Kα/Kβ の議論でも、「なぜ ‘1 本の線’ ではなく複雑な多電子構造になるのか」 を考えるうえで、この「独立粒子モデル vs 多電子相関」の視点がそのまま効いてきます。

(2) 第6回「原子の電子状態と遷移の計算のために」

こちらは、「実際に電子状態や遷移を計算したい人向けの実務ガイド」 になっています。

主なポイントは：

• 多電子原子の計算で必ず意識すべき 2 本柱

  1. 相対論効果（relativistic effect）

     • Z が大きくなると 1s 電子の速度が光速に近づき、
       Dirac–Coulomb ハミルトニアンで扱わないと 1％レベルでも怪しくなる。
     • 内殻励起（まさに Fe K など）では、上準位・下準位で相対論補正の入り方が違うため、
       「単なるオフセット」として無視できない。

  2. 電子相関（electron correlation）

     • HF/DF はあくまで独立粒子モデルであり、
       真の値との差が「電子相関エネルギー」として定義される。
     • 遷移エネルギー（上下準位の差）を問題にするときには、
       高 Z の多価イオンでも電子相関を無視できないことが強調される。

• 電子相関を扱う代表的な枠組み

  • 多配置 Hartree–Fock (MCHF)、配置間相互作用 (CI)
    → 電子配置をたくさん混ぜて相関を表現する方法。
  • 多体摂動法 (MBPT)
    → HF 解を基準に摂動次数を上げていくアプローチ。
  • Hylleraas 型波動関数や超球座標など、
    電子間距離 $r_{ij}$ を直接波動関数に書き込む相関の扱いにも触れている。

• 具体的な計算コードの紹介と考え方

  • 相対論をどこまで入れるか（Schrödinger vs Dirac）
  • 相関をどう入れるか（MCHF/CI, MBPT など）
  • その組み合わせで、原子構造コードを分類・評価する視点が整理されている。

この論文を読むと、

「Fe K の 1 eV レベルの議論をするときに、どの程度の相対論・相関を入れた計算が必要か？」

といった判断がなぜ必要で、なぜ難しいのかが理解しやすくなるでしょう。将来的に Cowan コードや FAC などを触ってみたい人にとって、“なぜそのコードがその近似を採用しているか” を日本語で把握できる入口としても有用かもしれません。