はじめに
XRR（X-Ray Reflectivity）は、薄膜試料の表面すれすれにX線を入射し、反射強度の角度依存性から膜厚・密度・表面/界面粗さを非破壊で推定する手法です。数nm〜数百nmの薄膜評価に広く使われ、ALDやスパッタ成膜の膜質管理には欠かせません。
解析にはrefl1dやGenXといった優れたOSSがありますが、Windows + Anaconda環境では依存関係のトラブルに悩まされがちです。また既製ツールはモデル定義がコード内に埋もれやすく、「成膜条件を変えた大量の試料を同じモデルで一括解析したい」という現場のニーズと噛み合わないことがあります。
そこで、以下の方針でXRR解析アプリを自作しました。
- 物理エンジン（Parratt漸化式）をNumPyだけで自前実装し、重い依存を排除
- 膜構造モデルはCSVファイルで定義（1行=1層、探索範囲とフィットフラグ付き）
- 最適化は差分進化法（scipy.optimize.differential_evolution）
- UIはStreamlit
コードはGitHubで公開しています。
https://github.com/yharada520/XRR-structure-analysis-app
XRRの原理（最小限）
物質のX線に対する屈折率は1よりわずかに小さく、
n = 1 − δ + iβ （δ ~ 10⁻⁵〜10⁻⁶）
δは電子密度（ほぼ質量密度）に比例します。臨界角 θc = √(2δ) 以下で全反射が起き、臨界角の位置から密度がわかります。膜があると膜上面と下面からの反射波が干渉し、Kiessigフリンジと呼ばれる振動が現れます。その周期から膜厚が、高角側の減衰の速さから粗さがわかります。
多層膜の反射率は、Parratt（1954）の漸化式で厳密に計算できます。界面粗さはNevot-Croce因子としてFresnel係数に乗じるのが標準的です。
設計1: 膜構造をCSVで定義する
このアプリの中心的な設計です。1行=1層、上から表面層→基板の順に書きます。
layer,formula,density_init,density_min,density_max,thickness_init,thickness_min,thickness_max,roughness_init,roughness_min,roughness_max,fit_density,fit_thickness,fit_roughness
Surface,C2H4O,1.2,0.8,1.8,1.0,0.3,3.0,0.3,0.1,1.0,1,1,1
Al2O3,Al2O3,3.0,2.4,3.5,30.0,20.0,40.0,0.5,0.1,2.0,1,1,1
SiOx,SiO2,2.2,2.0,2.4,1.5,0.5,3.0,0.3,0.1,1.0,0,1,1
Si_sub,Si,2.33,,,inf,,,0.2,0.1,1.0,0,0,1
ポイントは3つあります。
-
fit_*フラグで固定/可変を層×パラメータ単位で制御できる。「界面SiOxの密度は文献値固定、膜厚だけ振る」といった実務的なモデリングがCSV編集だけで済む - min/maxがそのまま差分進化の探索範囲になる。事前知識（成膜レートから膜厚は±30%以内、など）を自然に組み込める
- 光学定数δ・βは
formula（化学式）と密度からperiodictableパッケージで自動計算。ユーザーが散乱因子表を引く必要がない
設計2: Parratt漸化式のNumPy実装
エンジン本体は実質40行程度です。θ軸方向を完全ベクトル化し、ループは層方向のみにします。
def reflectivity_from_arrays(theta_deg, delta, beta, thickness_nm, roughness_nm, wavelength):
theta = np.radians(theta_deg)
n_media = len(delta) + 1 # 真空を含む媒質数
n = np.ones(n_media, dtype=np.complex128)
n[1:] = 1.0 - delta + 1j * beta
d_ang = thickness_nm * 10.0
sigma_ang = roughness_nm * 10.0
# kz[j] 形状: (n_media, n_theta)
cos2 = np.cos(theta)[np.newaxis, :] ** 2
kz = (2.0 * np.pi / wavelength) * np.sqrt(n[:, np.newaxis] ** 2 - cos2)
r_amp = np.zeros_like(kz[0])
for j in range(n_media - 2, -1, -1): # 基板側から表面へ
kj, kj1 = kz[j], kz[j + 1]
fresnel = (kj - kj1) / (kj + kj1)
fresnel = fresnel * np.exp(-2.0 * kj * kj1 * sigma_ang[j] ** 2) # Nevot-Croce
if j == n_media - 2:
r_amp = fresnel
else:
phase = np.exp(2j * kj1 * d_ang[j])
r_amp = (fresnel + r_amp * phase) / (1.0 + fresnel * r_amp * phase)
return np.abs(r_amp) ** 2
580点×4層で1回の評価が1ミリ秒以下。差分進化の数万回評価に余裕で耐えます。
設計3: δ・βの密度線形性で高速化
最適化ループ内で毎回periodictableを呼ぶと遅い、という問題があります。ここでδ・βは密度に厳密に線形（δ = 定数 × ρ）であることを使います。フィット開始前に単位密度あたりのδ・βを1回だけ計算しておき、ループ内は掛け算だけにします。
# 事前計算（層ごとに1回だけ）
delta_per_rho, beta_per_rho = delta_beta(formula, density=1.0, wavelength=1.5418)
# 最適化ループ内
delta = delta_per_rho * rho # periodictable呼び出しなし
これで6パラメータのフィットが約2〜4秒（評価回数 約1〜2万回）に収まります。
設計4: 目的関数はlog10スケール
XRRの反射率は6桁以上のダイナミックレンジを持つため、線形残差では臨界角近傍だけが支配的になり高角側のフリンジ（膜厚情報の本体）が無視されます。目的関数はlog10残差の二乗和にします。
chi2 = np.sum((np.log10(I_exp) - np.log10(I_calc))**2)
フィット対象は「CSVでフラグを立てたパラメータ + scale + background」を1本のベクトルにpackし、scipy.optimize.differential_evolution（init="sobol", polish=True）に渡します。DEは大域探索性が高く、XRRのような多峰性の目的関数（フリンジの「1本ずれ」に対応する局所解が無数にある）に適しています。
装置効果も忘れずに入れます: footprint補正（低角でビームが試料からはみ出す効果）、2θ分解能のGaussian畳み込み、スケール因子、バックグラウンドです。
検証: 物理テストと合成データ復元
数値エンジンのテストは「物理の解析解と照合する」のが鉄則です。pytestで以下を検証しました。
| テスト | 判定基準 | 結果 |
|---|---|---|
| Si臨界角（CuKα） | 0.223° ± 0.005° | 0.2233° |
| 単一界面Fresnel解析解 | 相対誤差 | < 1e-10 |
| Kiessigフリンジ周期（SiO₂ 20nm） | 屈折補正後 2π/d ± 1% | 一致 |
| Q⁻⁴減衰、粗さによる単調減衰 | 定性 | 一致 |
仕上げは合成データ復元テストです。真値既知の構造から3%ノイズ付きデータを生成し、初期値をずらして6パラメータ同時フィット:
| パラメータ | 真値 | 復元値 | 誤差 |
|---|---|---|---|
| SiO₂膜厚 | 21.30 nm | 21.299 nm | −0.00% |
| SiO₂密度 | 2.150 g/cm³ | 2.1497 | −0.01% |
| SiO₂粗さ | 0.40 nm | 0.400 nm | +0.03% |
ハマりどころ
(1) フリンジ間隔は低角側で2π/dより最大8%狭い
Kiessigフリンジ周期のテストが最初は通りませんでした。原因はバグではなく屈折効果です。フリンジ極小は膜内部の波数で等間隔に並ぶため、真空側のQで見ると低角ほど間隔が縮みます。屈折補正 Q′ = √(Q² − Qc²) を適用すると厳密に等間隔になります。実データでフリンジ間隔から膜厚を概算するときも、この補正を忘れると薄い膜ほど過大評価します。
(2) 極小位置のグリッド量子化
角度ステップが有限なので、検出した極小位置は±1グリッドの誤差を持ちます。個々のフリンジ間隔で判定するとばらつきで落ちるため、フリンジ次数に対する回帰直線の傾きで判定するようにしました。
(3) Rigaku .rasの3列目
.rasのデータ行は「2θ 強度 アッテネータ係数」の3列のことがあります。実強度に戻すには 強度×係数 が必要です。2列固定で読むと減衰器が入った角度域で強度が桁で狂います。
Streamlit UI
データアップロード → 構造の表エディタ編集 → フィット（収束カーブをライブ表示）→ 結果（実測vs計算、残差、密度深さプロファイル、パラメータ表）→ Excel/CSVエクスポート、という流れをワンページにまとめています。
実用上効いているのは**「探索範囲の端に張り付いたパラメータの警告」**です。フィット値が範囲端にあるのはモデルか範囲設定に無理があるサインなので、自動検出して表示します（background=0のような下限0への正常収束は除外）。
今後
- BICによるモデル自動選択（表面層・界面層の有無を4モデル比較）
- pygadによる本格的なGA選択肢
- ブートストラップ/MCMCによる不確かさ定量化
まとめ
- Parratt漸化式はNumPyで40行、自前実装は十分現実的
- δ・βの密度線形性を使えばDEの数万回評価が数秒で回る
- 目的関数はlog10残差一択、装置効果（footprint・分解能）も必須
- テストは物理の解析解と合成データ復元で固める
参考文献
- 表 和彦, X線薄膜測定法 第5回 X線反射率測定法, リガクジャーナル 40(2), 1–9 (2009)
- L. G. Parratt, "Surface Studies of Solids by Total Reflection of X-Rays", Phys. Rev. 95, 359 (1954)
- L. Névot, P. Croce, Rev. Phys. Appl. 15, 761 (1980)