はじまるよ！
どうもこんにちは。前回(かなり前ですが)に引き続き、今回もPythonで結晶構造解析です。コードはGitHub¹

次の記事: Pythonで電子顕微鏡シミュレーション：マルチスライス法(2)

注意点：　趣味でシミュレーション(趣味レーション)していると、結果が本当に正しいのか判別がつきにくく、ここに載っている結果も物理的に正しくない可能性があります。間違いがあったらご指摘をいただけると幸いです。

#電子顕微鏡
　電子顕微鏡は、電子線を対象にぶつけることで、微細な物を見たり分析したりする装置です。光学顕微鏡や前回の記事で取り上げたX線回折など、微細な物を観察する手法は他にもありますが、電子線は可視光線やX線と比較して波長が短く、より高い分解能を発揮します。その分解能を生かして、物性物理における構造解析や生化学における生体組織観察など、様々な場面で使われています。
　電子顕微鏡には大きく分けて走査型電子顕微鏡(SEM)、透過型電子顕微鏡(TEM)、走査透過型電子顕微鏡(STEM)の三種類があります。

#透過型電子顕微鏡(Transmission Electron Microscope, TEM)
　透過型電子顕微鏡(TEM)は、対象物に電子を照射し、透過した電子で像を作る顕微鏡です。透過電子線(透過波)の強度から内部構造を調べたり、電子線の回折(回折波)から結晶構造を特定したりします。
　さらにその発展として、透過波と回折波を干渉させて像を作る、高分解能電子顕微鏡法(High-Resolution electron microscopy: HREM)という手法があります。HREMを用いることで、個々の原子配列をそのまま映し出すほどの高い分解能が得られます。HREMによって作られる像を結晶構造像と言います。今回は、結晶構造像をシミュレーションで作ることを目的とします。

#回折波と透過波の解説
　電子線は波の性質を持っていますので、障害物を廻り込んで伝播する、回折現象を起こします。結晶中の障害物とは、原子の周りにある電子雲(ポテンシャル)です。結晶では原子は規則的に並んでいるので、ポテンシャルも周期的に分布しており、それによる回折は
$$2dsin\theta = n \lambda$$
で示されるブラッグ条件を満たしたときのみ回折波として出てきます。このような反射をブラッグ反射と言います。dは結晶面間隔、θは結晶面と入射電子線のなす角度、λは電子線の波長です。

　ところで、TEMで対象とする試料はとても薄く、nmオーダーの厚みしかありません。そのような試料に電子線を照射すると、反射や吸収、散乱を起こさずに透過するものも出てきます。このような電子線を透過波と言います。
　回折波と透過波では波の位相が違うので、干渉させることで位相の差を像に反映させることができます(実はそのまま干渉させるだけでは位相差は反映されないのですが、後述する球面収差を用いることで可能です)。

#マルチスライス法
ここからは、シミュレーションの具体的な手法についての説明です。
電子線が試料に入射し、透過または回折して試料下面で像を作るとき、電子線は試料内のポテンシャルから相互作用を受けて出力されます。シミュレーションでそれを再現するには、試料内部のポテンシャルを再現する必要がありますが、数十、数百 nmの厚さのある試料のポテンシャルを完全に再現するには、途方もない数の原子の、それに伴う電子雲からの相互作用を計算しなければなりません。
そこで、試料を電子線の入射方向に対して垂直にスライスし、薄いスライスが積み重なったものとして計算する、マルチスライス法が考えられました。電子線はスライスを通過することで散乱され位相変化を受け、次のスライスまで空間を伝播し、次のスライスでまた位相変化を受け……と繰り返すことで試料下面では結晶構造像が作られます。下はちょっとした模式図です。

　簡単のために二枚のスライスを通過して下面に像を作ることを考えてみましょう。
　入射波$\psi_0$が一枚目のスライスに入射すると、スライス内のポテンシャル$V_p$から位相変化のみを受けるとする、位相物体近似を適用します。
投影ポテンシャルに比例した位相変化を受けるとき、スライス下面での波は
$$ \psi_{1out}(x, y) = e^{i\sigma V_{p1}(x, y)} $$
となります。$\sigma = \pi/\lambda E$($\lambda$: 波長, E: 加速電圧)です。さらに、$V_p(x, y)$から受ける変化は小さいとする弱位相物体近似を適用すると、
$$ \psi_{1out} (x, y) = \{1 + i\sigma V_{p1}(x, y) \} $$
$$ q_1(x, y) = 1 + i\sigma V_{p1}(x, y) $$
$$ \psi_{1out}(x, y) = q_1(x, y)\psi_0(x, y) $$
となります。$q(x, y)$を透過関数と言います。
　一枚目のスライスを通過した$\psi_1$が二枚目のスライスまでの空間を伝播します。波の伝播を表す伝播関数$p$は
$$ p(x, y) = -\frac{1}{i\lambda \Delta z}e^{-ik\frac{(x^2 + y^2)}{2\Delta z}} $$
です。これはフラウンホーファー回折の式でもあります。ここで、$p(x, y)$をフーリエ変換しておきます。後々効いてきます。
$$ P(h, k) = -\frac{1}{i\lambda \Delta z}\int e^{-ik\frac{(x^2 + y^2)}{2\Delta z}} e^{-2\pi i(hx + ky)} dxdy = e^{2\pi i \zeta(h, k) \Delta z } $$
ここで$\zeta(h, k)$は
$$ \zeta(h, k) = \frac{\lambda}{2}(h^2 + k^2) $$
で、励起誤差と呼ばれています。以上から、二枚目のスライスに到達した波は
$$ \psi_{2in}(x, y) = \psi_{1out}(x, y) \otimes p(x, y) $$
となります。$\otimes$は畳み込みを表します。
　二枚目のスライスを通過した波は、一枚目と同様に
$$ \psi_{2out}(x, y) = q_2(x, y)\psi_{2in}(x, y) $$
です。このように、スライスの数だけ計算を繰り返します。
　ところが、先ほどの計算では、スライスの数だけ畳み込み演算をしなければなりません。一般に畳み込みは計算コストが高く、直接計算するのは避けるべきです。そこで、以下の式を使い、畳み込みをフーリエ変換に置き換えます。
$$ F[f(x)\otimes g(x)] = F[f(x)]F[g(x)] $$
すると、今までの式は
$$ \Psi_{out}(h, k) = F[q_2(x, y)F^{-1}[F[q_1(x, y)\psi_0(x, y)]P(h, k)]] $$
となり、フーリエ変換と逆フーリエ変換を繰り返すことで計算することができます。

#プログラム
　マルチスライス法での計算では、まず投影ポテンシャルを用意しなければなりません。投影ポテンシャルは結晶構造因子を逆フーリエ変換して得られます(そもそも、結晶構造因子は内部ポテンシャルをフーリエ変換したものです)。よってまず、結晶構造因子を計算して、逆フーリエ変換します。結晶構造因子の説明や計算方法については前回の記事にあります。
　今回は試料をα-Fe(bcc)とし、入射方向はz軸方向である[001]とします。スライスの間隔$\Delta z$は格子定数/2とします。原子の座標はpos_1.csvとpos_2.csvに予め入力してあります。
　逆格子はnp.meshgrid()で生成し、np.vectorize()で適用します。(前回はnp.vectorize()について知らなかったので、わざわざ列挙してました。無駄でしたね……)

fa = 0
hkl = [h, k, 0]
thkl = np.transpose(hkl)
dk = 1/((np.matmul(np.matmul(invG, thkl), hkl))**(1/2))
if(not np.isinf(dk)):
 sinthetalamb = lamb/(2*dk)
 for i in range(4):
  fa = fa + param_fe[i][0]*np.exp((-param_fe[i][1]*(sinthetalamb**2)))
            
 F = 0
 for i in r:
  F = F + fa*np.exp(1.0j*2*np.pi*(h*i[0] + k*i[1]))

結晶構造因子$F$をnp.fft.ifft2()で投影ポテンシャル$V_p$とし、透過関数$q$を計算します。

Vp = np.fft.ifft2(F)
q = 1 + 1.0j*sigma*Vp

伝播関数Pは最初からフーリエ変換したかたちで計算します。

zeta = lamb/2
zeta = zeta*(h**2 + k**2)
        
p = np.exp(1.0j*2*np.pi*zeta*deltaZ)

マルチスライス計算を行います。入射電子線$\psi_1$は1とします。

psi = np.ones((self.N, self.N))
psi = np.fft.fft2(q*psi)*p
psi = np.fft.ifft2(psi)

これをスライスごとに繰り返します。マルチスライス計算から出力される$\psi$は複素数です。実際に像として現れる強度は波の振幅$I$ですので、

I = np.abs(psi)**2

とすると、結晶構造像が得られる……わけではありません。

#(2)へつづく
　すいません、ここまで長々と説明をしてきましたが、今回の計算では結晶構造像が(正確には、ちゃんとした結晶構造像が)得られません。実はまだ考慮すべき要因があります。長くなってしまいましたので、次の記事に続きます。

#参考文献
　記事を書くにあたり参考にした記事や書籍を勝手ながら紹介させていただきます。

日本電子株式会社　用語集
電子顕微鏡メーカーといえば日本電子さんですが、ホームページの用語集がとても充実しています。電子顕微鏡について勉強している人にとってはかなり役に立つんじゃないでしょうか。
レンズの点像分布関数とMTF曲線の数値計算(回折計算編)
回折理論についてわかりやすい説明をされています。この記事でnp.vectorize()の存在を知ったのですが、すごい便利ですね。

#コード

https://github.com/PicricAcid/multislice_method

https://github.com/PicricAcid/multislice_method
に載せています。 ↩