0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

EDSスペクトル・シミュレータを作った ― 計数統計と薄膜マトリクス効果をブラウザで可視化する(Streamlit + xraylib)

0
Posted at

🇯🇵 日本語版(この下) / 🇬🇧 English version below

  • SEM-EDS(エネルギー分散型X線分析)の現場でありがちな2つの誤解を、定量的に体験できる Web アプリを作りました。
    • 「ピークが見えない=その元素は無い」→ 実は √N の統計ノイズに埋もれているだけかもしれない
    • 「加速電圧は高いほど良い」→ 薄膜では情報体積の大半が基板。電圧を下げると薄膜の S/N が相対的に上がる
  • スタックは Python / Streamlit / Plotly / xraylib / SciPy / Docker
  • 物理は Kramers 制動X線、ポアソンノイズ、Fano 分解能、Packwood-Brown の φ(ρz)、Kanaya-Okayama 飛程、Bragg 加算則、電子衝突電離(Bethe形)。
  • 実装中に踏んだ2つの落とし穴(XRF断面積を使いかけた話 / scipy.quad の無限区間で積分が消える話)も共有します。
  • 触ってみる(インストール不要): https://edssimulator-xjd8gyxkgdwezujsfevvha.streamlit.app/
  • リポジトリ: https://github.com/yharada520/EDS_simulator (Apache-2.0)

低カウントと十分なカウントの比較

同じ試料(Si 99% + 微量 Ti 1%)でも、測定条件(線量)を変えるだけで微量ピークが統計ノイズに埋もれたり出現したりする。左は Max 約2,000カウント、右は約90,000カウント。


なぜ作ったか

材料分析の現場では、EDS の計数統計や電子線の相互作用領域についての誤解が意外と根強く残っています。典型例を2つ挙げます。

  1. 微量元素の「不在」の誤認。1% 含有している元素が、Max 2,000カウント程度の不十分な条件では、バックグラウンドのポアソンゆらぎ(√N)に埋もれて見えない。これを「ピークが無いから存在しない」と結論してしまう。
  2. 加速電圧の取り過ぎ。ALD などで成膜した極薄膜を分析するとき、情報体積の大部分が基板由来(マトリクス効果)であることを理解せず、加速電圧を高く設定してしまう。

どちらも「言葉で説明されても腹落ちしにくい」タイプの話です。そこで、スライダーを動かすと理由が目に見える教育・R&D 用ツールを作りました。

何を作ったか

Streamlit の Web アプリで、2つのモードを持ちます。

  • 均質バルク(統計)モード: 加速電圧・プローブ電流・積算時間・組成を振ると、理論スペクトルにポアソンノイズが重畳され、微量ピークが埋もれる/出現する様子が見えます。
  • 薄膜/基板(多層)モード: 表面から基板まで最大5層を積み、各層を化学式(TiN, TiO2, Al2O3 など)で指定できます。加速電圧を下げると基板シグナルが減り、埋もれた密着層が相対的に見えやすくなる挙動を再現します。

アプリの画面

左サイドバーで条件を、中央にスペクトル、右に Max カウントや各ピークの簡易 S/N を表示します。

技術スタック

領域 使用技術
UI / 描画 Streamlit, Plotly
数値・物理 NumPy, SciPy
原子データ xraylib(ライン energy、吸収端、蛍光収率、質量減衰係数)
環境 Docker / conda / pip、GitHub Actions(CI)

xraylib は conda-forge が定番ですが、実は PyPI にも各 OS・Python 3.10〜3.13 の wheel が揃っており、pip install xraylib でそのまま入ります。未導入環境でもフェーズ1が動くよう、遅延インポート+内蔵フォールバックにしてあります。

物理モデルの実装

1. 計数統計(ポアソンノイズ)

制動X線(バックグラウンド)は Kramers の法則の近似です。

$$I(E) \propto \bar{Z},\frac{E_0 - E}{E}\quad(E < E_0)$$

特性X線はエネルギー依存の分解能を持つ Gaussian。EDS 分解能は Fano 統計の標準式で、Mn Kα (5.9 keV) の 130 eV を基準に較正します。

$$\mathrm{FWHM}(E)^2 = (2.3548)^2,\varepsilon,F,E + N^2$$

($\varepsilon=3.85$ eV は Si の電子正孔対生成エネルギー、$F=0.12$ は Fano 係数)。理論スペクトルを組み立てたら、最後に numpy.random.poisson でショットノイズを重畳します。

rng = np.random.default_rng(seed)
observed = rng.poisson(np.clip(theoretical, 0.0, None)).astype(float)

Max カウントは intensity_scale × (プローブ電流 × 積算時間) に線形連動させ、最強ピークをこの値に正規化しています。これで「Max カウントを増やすと √N に対してピークが立ってくる」挙動が直感的に効きます。

2. 特性X線の相対強度 ― ここで最初の落とし穴

多殻元素(Au の M 線と L 線など)の相対強度を正しく出したい。最初、xraylib の CS_FluorLine_Kissel(Z, line, E) を使おうとしました。ところが検証すると、Ti Kα が 5 kV(吸収端直上)で最大という挙動をする。

これは CS_FluorLine_Kissel光子励起(XRF) の蛍光生成断面積だからです。光電吸収は吸収端直上で最大、その上で E⁻³ 的に減衰します。しかし SEM-EDS は電子線励起。電子衝突電離は閾値(過電圧 $U=E_0/E_c=1$)でゼロ、$U\approx e$ 付近で最大という逆の挙動です。そのまま使えば非物理になるところでした。

正しくは Bethe 形の電子衝突電離断面積 × 蛍光収率 × 遷移確率です。

def electron_impact_weight(xl, z, series, edge_kev, rad_rate, e0_kev):
    u0 = e0_kev / edge_kev
    if u0 <= 1.0:
        return 0.0  # 過電圧不足で励起されない
    n_i = _SHELL_ELECTRONS[series]           # 副殻の電子数
    ionization = n_i * math.log(u0) / (u0 * edge_kev**2)  # ∝ 電子衝突電離
    omega = _fluor_yield(xl, z, series)      # 蛍光収率(xraylib.FluorYield)
    return ionization * omega * rad_rate     # × 遷移確率(xraylib.RadRate)

結果、Au は 15 kV で M/L 比 ≈ 11.6、10 kV 以下では L 殻(端 ~11.9 keV)が励起されず M 線のみ ―― 実機どおりの挙動になりました。**「便利な関数があるからといって、対象の物理(光子励起 vs 電子励起)を取り違えない」**という良い教訓でした。

3. 薄膜/基板の φ(ρz) ― ここで2つ目の落とし穴

X線発生の深さ分布 φ(ρz) は Packwood-Brown 型の表面中心 Gaussian を採用しました。

$$\phi(\rho z) = \gamma,e^{-\alpha^2(\rho z)^2}\left[1 - \frac{\gamma-\phi_0}{\gamma}e^{-\beta\rho z}\right]$$

深さスケール(α)は材料屋に馴染みの深い Kanaya-Okayama 電子飛程に固定し、加速電圧依存性 $R \propto E_0^{1.67}$ を担保。薄膜層・基板層の放出強度を scipy.integrate.quad で深さ積分します。

ここで盛大にハマりました。基板積分の上限を np.inf にしていたところ、加速電圧を下げて φ の Gaussian 幅(1/α)が狭くなると、積分値が桁で消失したのです(15 kV で本来 5.5e-4 になるはずが 8.5e-13)。

原因は quad の適応サンプリング。無限区間 [a, ∞) に対して鋭いピークが下端付近に立つと、変数変換後のサンプル点がピークを取りこぼし、被積分関数がアンダーフローして 0 と評価されていました。

対策は単純で、上限を φ のガウス裾に応じた有限値($10/\alpha$、$e^{-100}\approx0$)にすること。

def _gaussian_depth_cap(p):     # exp(-alpha^2 rhoz^2) が消える深さ
    return 10.0 / max(p.alpha, 1e-9)

# quad の上限は必ず有限に
val, _ = quad(integrand, lo, min(hi, cap), limit=200)

副次効果として「膜厚が電子飛程を超えると基板シグナルが 0 になる」という物理も正しく表現できるようになりました。無限区間 × 鋭いピークの数値積分は、有限化するかピーク位置を points で与える、という定番の注意点です。

4. 化合物層 ― Bragg 加算則

各層を化学式でパースし、質量分率に応じて構成元素が発光、吸収は Bragg の加算則で計算します。

$$\left(\frac{\mu}{\rho}\right)_\text{compound} = \sum_i w_i \left(\frac{\mu}{\rho}\right)_i$$

各層の放出X線は、上側の全層の(質量厚 × その層の μ/ρ)で減衰させます ―― 埋もれた層ほど見えにくく、深すぎれば信号ゼロ。密度はプリセット表→単一元素の元素密度→仮の既定値、の順で決めます。

なお密度でも小さくハマりました。プリセット外化合物で体積加算則 $1/\rho=\sum w_i/\rho_i$ を使ったら CuO で 0.01 g/cm³ という異常値。xraylib が O や N を気体密度(O = 0.00133)で返すためでした。ここは素直に、プリセットに無ければ仮の既定値+「密度は手入力推奨」の警告にしています。

アーキテクチャ

物理演算は eds_sim/ パッケージに機能分割し、app.py は UI と描画に専念させています。

eds_sim/
├── config.py         # 定数・データクラス設定
├── continuum.py      # Kramers 制動X線
├── characteristic.py # 特性X線ライン → Gaussian(バルク/層構造)
├── composition.py    # 化学式パース・化合物密度/質量吸収(Bragg 加算則)
├── detector.py       # 分解能・窓吸収・効率
├── depth.py          # φ(ρz) 深さ分布・薄膜/基板の深さ積分
├── elements.py       # xraylib ヘルパ(遅延import・フォールバック)
└── spectrum.py       # 合成 + ポアソンノイズ

xraylib を遅延インポートし、未導入なら内蔵テーブルにフォールバックする設計にしたことで、素の pip 環境や CI でもフェーズ1が完全動作します。テストは xraylib の有無で skipif を切り替えています。

テスト・CI・実行

スモークテストを pytest で用意し、GitHub Actions で Python 3.10 / 3.11 / 3.12 のマトリクス実行しています(requirements.txt に xraylib を含めたので、CI 上でも電子衝突電離モデルの回帰テストが走ります)。

起動は Docker が手軽です。

docker compose up --build   # → http://localhost:8501

pip でも動きます(xraylib 込みでフル精度)。

pip install -r requirements.txt
streamlit run app.py

モデルの限界(正直な注記)

  • P/B 比・絶対カウントは可視化向けの経験較正であり、第一原理の定量値ではありません。
  • 電子散乱マトリクスは基板組成で近似(薄膜 ρz ≪ 電子飛程を想定)。厚い重元素最上層では電子減速を過小評価しえます。
  • 電子衝突電離は Bethe 形の近似。より高精度には Casnati / Bote-Salvat や Coster-Kronig の導入余地があります。
  • 化合物式は単純式のみ(括弧・水和物は今後対応)。

まとめ

「なんとなく」で流されがちな計数統計とマトリクス効果を、スライダー操作で腹落ちさせるツールを作りました。物理の実装以上に、光子励起 vs 電子励起の取り違え無限区間の数値積分といった実装上の落とし穴が学びになりました。OSS(Apache-2.0)で公開しているので、教育や社内勉強会に使えそうであれば覗いてみてください。



A browser-based EDS spectrum simulator ― visualizing counting statistics and thin-film matrix effects (Streamlit + xraylib)

🇬🇧 English version / 🇯🇵 日本語版は上部にあります。

TL;DR

  • A web app that lets you feel, quantitatively, two common misconceptions in SEM-EDS (energy-dispersive X-ray spectroscopy):
    • "No peak, so the element is absent." A 1 % element can simply be buried in the √N counting noise.
    • "Higher kV is always better." For a thin film, most of the interaction volume is the substrate. Lowering the accelerating voltage improves the film's S/N.
  • Stack: Python / Streamlit / Plotly / xraylib / SciPy / Docker.
  • Physics: Kramers bremsstrahlung, Poisson noise, Fano resolution, Packwood-Brown φ(ρz), Kanaya-Okayama range, Bragg additivity, and electron-impact ionization (Bethe form).
  • I also share two implementation gotchas: almost using an XRF cross-section for an electron beam, and a scipy.quad integral that silently collapsed over an infinite interval.
  • Live demo (no install): https://edssimulator-xjd8gyxkgdwezujsfevvha.streamlit.app/
  • Repository: https://github.com/yharada520/EDS_simulator (Apache-2.0)

Under- vs well-counted spectra

The same sample (Si 99 % + a 1 % Ti trace). Only the measurement condition (dose) differs — left ≈ 2,000 max counts, right ≈ 90,000. The trace peak hides, then emerges.

Why I built it

Two failure modes show up again and again in the analysis lab:

  1. Mistaking "invisible" for "absent." At a few thousand counts, a 1 % element sits below the Poisson fluctuation (√N) of the background, and people conclude "there is no peak, so it isn't there."
  2. Cranking up the accelerating voltage. When analyzing an ALD-thin film, much of the signal actually comes from the substrate (a matrix effect), yet the voltage is set too high.

Both are hard to internalize from words alone. So I built an educational / R&D tool where moving a slider makes the reason visible.

What it is

A Streamlit web app with two modes:

  • Homogeneous bulk (statistics): vary accelerating voltage, probe current, live time and composition; Poisson noise is applied to the theoretical spectrum, so you watch a trace peak drown or surface.
  • Thin film / substrate (multilayer): stack up to 5 layers (surface → substrate), each entered as a chemical formula (TiN, TiO2, Al2O3, …). Lowering the voltage suppresses the substrate signal and relatively enhances a buried adhesion layer.

The app UI

Controls on the left sidebar, spectrum in the center, and max counts / per-peak S/N on the right.

Tech stack

Area Tools
UI / plotting Streamlit, Plotly
Numerics / physics NumPy, SciPy
Atomic data xraylib (line energies, edges, fluorescence yields, mass attenuation)
Env Docker / conda / pip, GitHub Actions (CI)

xraylib is usually installed from conda-forge, but it also ships PyPI wheels for every OS and Python 3.10–3.13, so pip install xraylib just works. It is lazily imported with a built-in fallback table, so Phase 1 runs even without it.

Implementing the physics

1. Counting statistics (Poisson noise)

The continuum (bremsstrahlung) uses Kramers' law:

$$I(E) \propto \bar{Z},\frac{E_0 - E}{E}\quad(E < E_0)$$

Characteristic lines are Gaussians with an energy-dependent width. The EDS resolution follows the standard Fano-statistics formula, calibrated to 130 eV at Mn Kα (5.9 keV):

$$\mathrm{FWHM}(E)^2 = (2.3548)^2,\varepsilon,F,E + N^2$$

(ε = 3.85 eV is the Si electron–hole pair energy, F = 0.12 the Fano factor). Finally, shot noise is applied with numpy.random.poisson:

rng = np.random.default_rng(seed)
observed = rng.poisson(np.clip(theoretical, 0.0, None)).astype(float)

Max counts scale linearly with intensity_scale × (probe current × live time), and the strongest peak is normalized to that value — so "more counts → the peak rises above √N" works intuitively.

2. Relative line intensities ― the first gotcha

I wanted correct shell-to-shell intensities for multi-shell elements (e.g. Au M vs L). My first attempt used xraylib's CS_FluorLine_Kissel(Z, line, E). But on inspection, Ti Kα peaked at 5 kV (just above the edge) — backwards.

The reason: CS_FluorLine_Kissel is a photon-excited (XRF) fluorescence cross-section. Photoionization is maximal just above the edge and falls off as ~E⁻³. But SEM-EDS is electron-excited, and electron-impact ionization is zero at threshold (overvoltage U = E₀/E_c = 1) and peaks around U ≈ e — the opposite trend. Using it as-is would have been unphysical.

The correct quantity is a Bethe-form electron-impact cross-section × fluorescence yield × radiative rate:

def electron_impact_weight(xl, z, series, edge_kev, rad_rate, e0_kev):
    u0 = e0_kev / edge_kev
    if u0 <= 1.0:
        return 0.0  # below threshold: not excited
    n_i = _SHELL_ELECTRONS[series]           # electrons in the subshell
    ionization = n_i * math.log(u0) / (u0 * edge_kev**2)  # ∝ electron impact
    omega = _fluor_yield(xl, z, series)      # fluorescence yield (xraylib.FluorYield)
    return ionization * omega * rad_rate     # × radiative rate (xraylib.RadRate)

Now Au has an M/L ratio ≈ 11.6 at 15 kV, and below 10 kV the L shell (edge ~11.9 keV) isn't excited so only M lines remain — exactly what a real instrument shows. Lesson: a convenient function is no substitute for matching the actual physics (photon vs electron excitation).

3. The thin-film φ(ρz) ― the second gotcha

The depth distribution of X-ray generation, φ(ρz), uses a Packwood-Brown surface-centered Gaussian:

$$\phi(\rho z) = \gamma,e^{-\alpha^2(\rho z)^2}\left[1 - \frac{\gamma-\phi_0}{\gamma}e^{-\beta\rho z}\right]$$

The depth scale (α) is anchored to the Kanaya-Okayama electron range — familiar to materials people — giving the voltage dependence $R \propto E_0^{1.67}$. Emission from each layer is integrated over depth with scipy.integrate.quad.

This is where I got burned. With the substrate integral's upper limit set to np.inf, lowering the voltage narrows φ's Gaussian (1/α), and the integral collapsed by orders of magnitude (8.5e-13 where it should be ~5.5e-4 at 15 kV).

The cause is quad's adaptive sampling: over a semi-infinite interval [a, ∞) with a sharp peak near the lower bound, the transformed sample points miss the peak and the integrand underflows to 0.

The fix is trivial — cap the upper limit at a finite value set by the Gaussian tail ($10/\alpha$, $e^{-100}\approx 0$):

def _gaussian_depth_cap(p):     # depth where exp(-alpha^2 rhoz^2) vanishes
    return 10.0 / max(p.alpha, 1e-9)

# always give quad a finite upper limit
val, _ = quad(integrand, lo, min(hi, cap), limit=200)

A nice side effect: "a film thicker than the electron range yields zero substrate signal" now falls out naturally. The takeaway: for a sharp peak on an infinite interval, make it finite or pass the peak location via points.

4. Compound layers ― Bragg additivity

Each layer is parsed from a chemical formula; its elements emit in proportion to their mass fractions, and absorption uses Bragg's additivity rule:

$$\left(\frac{\mu}{\rho}\right)_\text{compound} = \sum_i w_i \left(\frac{\mu}{\rho}\right)_i$$

Emission from each layer is attenuated by all overlying layers (mass thickness × that layer's μ/ρ) — buried layers are harder to see, and vanish if too deep. Density is resolved as preset table → single-element density → nominal default.

Density had its own trap: for compounds outside the preset, the volume-additivity rule $1/\rho=\sum w_i/\rho_i$ gave 0.01 g/cm³ for CuO — because xraylib returns O and N at gas density (O = 0.00133). So I fall back to a nominal value plus a "please enter the density" warning.

Architecture

Physics lives in the eds_sim/ package; app.py only does UI and plotting.

eds_sim/
├── config.py         # constants / dataclasses
├── continuum.py      # Kramers bremsstrahlung
├── characteristic.py # characteristic lines → Gaussians (bulk / layered)
├── composition.py    # formula parsing, compound density / MAC (Bragg)
├── detector.py       # resolution, window absorption, efficiency
├── depth.py          # phi(rho z) and thin-film / substrate depth integrals
├── elements.py       # xraylib helpers (lazy import, fallback)
└── spectrum.py       # assembly + Poisson noise

Lazy-importing xraylib with a fallback table means Phase 1 fully runs on plain pip or in CI; tests toggle with skipif on xraylib availability.

Tests, CI, running it

Smoke tests run under pytest on GitHub Actions across Python 3.10 / 3.11 / 3.12 (xraylib is in requirements.txt, so the electron-impact regression test runs in CI too).

Docker is the easiest way to run it:

docker compose up --build   # → http://localhost:8501

pip works too (full accuracy with xraylib):

pip install -r requirements.txt
streamlit run app.py

Honest limitations

  • P/B ratio and absolute counts are an empirical calibration for visualization, not first-principles quantities.
  • The electron-scattering matrix is approximated by the substrate composition (thin films, ρz ≪ range); thick heavy top layers underestimate electron slowing-down.
  • Electron-impact ionization is a Bethe-form approximation; Casnati / Bote-Salvat and Coster-Kronig would refine it.
  • Compound formulas are simple only (parentheses / hydrates are future work).

Wrap-up

I built a tool that turns "vaguely known" counting statistics and matrix effects into something you can grasp by moving a slider. Beyond the physics, the real lessons were the implementation traps: photon vs electron excitation and numerical integration over an infinite interval. It's open source (Apache-2.0) — take a look if it's useful for teaching or a study group.


Tags: Python, Streamlit, MaterialsScience, Simulation, xraylib

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?