【Python】凸型ポテンシャル障壁へ照射した平面波の時間発展（量子力学）

2025年1月20日 13:49

今回は上図のように有限の高さと幅をもつポテンシャル障壁が存在する「箱型ポテンシャル障壁」に量子粒子を照射するよ。このような箱型ポテンシャル障壁が存在する系でもそれぞれの領域で存在しうる平面波の重ね合わせで波束の運動をシミュレーションすることができるよ。箱型ポテンシャル障壁を次のように定義するよ。

$$
V(x) = \left\{ \begin{matrix} V_{\rm I} =0 & (x<0) \cr V_{\rm II}=V & (0 \le x < d) \\ V_{\rm III}=0 & (d\le x) \end{matrix} \right.
$$

平面波の解析解

上図は３つの領域における単一平面波の模式図だよ。入射波、反射波、透過波の波動関数をそれぞれ$${\varphi_I(x,t)}$$、$${\varphi_R(x,t)}$$ 、$${\varphi_T(x,t)}$$と表し、$${0\le x < d}$$の前進波と後進波の波動関数をそれぞれ$${\varphi_+(x,t)}$$、$${\varphi_-(x,t)}$$と表した場合、各領域の波動関数は次の通りになるよ。

$$
\psi(x,t) =\left\{ \begin{matrix} \psi_{\rm I}(x,t) = \psi_I(x,t) +\psi_R(x,t) \ \ \, & (x<0) \cr \psi_{\rm II}(x,t) = \psi^{(+)}(x,t) +\psi^{(-)} (x,t)& (0\leq x\leq d) \cr \psi_{\rm III}(x,t) = \psi_T(x,t) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \, & (x>d) \end{matrix} \right.
$$

上記の具体的な表式は次のとおりだよ。

$$
\psi_{\rm I}(x,t) = \left[ e^{ik_{\rm I}x} +re^{-ik_{\rm I}x}\right] e^{-i\omega t}
$$

$$
\psi_{\rm II}(x,t) = \left[Ae^{ik_{\rm II}x} +Be^{-ik_{\rm II}x}\right] e^{-i\omega t}
$$

$$
\psi_{\rm III}(x,t) = \left[te^{ik_{\rm III}(x-d)} \right] e^{-i\omega t}
$$

これらの未知の係数$${r, t, A, B}$$は各領域間でなめらかに接続する条件（連続条件と傾きが連続する条件）

$$
\psi_{\rm I}(0,t) = \psi_{\rm II}(0,t)
$$

$$
\psi_{\rm II}(d,t) = \psi_{\rm III}(d,t)
$$

$$
\left.\frac{ \partial \psi_{\rm I}(x,t)}{\partial x}\right|_{x=0} = \left.\frac{ \partial \psi_{\rm II}(x,t)}{\partial x}\right|_{x=0}
$$

$$
\left.\frac{ \partial \psi_{\rm II}(x,t)}{\partial x}\right|_{x=d} = \left.\frac{ \partial \psi_{\rm III}(x,t)}{\partial x}\right|_{x=d}
$$

から決定することができるね。これらから得られる方程式は次のとおりだよ。

$$
1+r= A+B
$$

$$
Ae^{ik_{\rm II}d} +Be^{-ik_{\rm II}d} =t
$$

$$
k_{\rm I}(1-r)= k_{\rm II}(A-B)
$$

$$
k_{\rm II} \left[ Ae^{ik_{\rm II}d} -Be^{-ik_{\rm II}d} \right] =k_{\rm III} t
$$

この連立方程式を解くことで、$${A, B, r, t}$$が以下の通り得られるね。

$$
A =\frac{ 2k_{\rm I}(k_{\rm II} +k_{\rm III}) }{(k_{\rm I} -k_{\rm II})(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d} +(k_{\rm I} +k_{\rm II})(k_{\rm II} +k_{\rm III})}
$$

$$
B =\frac{ 2k_{\rm I}(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d}}{(k_{\rm I} -k_{\rm II})(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d} +(k_{\rm I} +k_{\rm II})(k_{\rm II} +k_{\rm III})}
$$

$$
r =\frac{ (k_{\rm I} +k_{\rm II})(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d} + (k_{\rm I} -k_{\rm II})(k_{\rm II} +k_{\rm III}) }{(k_{\rm I} -k_{\rm II})(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d} +(k_{\rm I} +k_{\rm II})(k_{\rm II} +k_{\rm III})}
$$

$$
t =\frac{ 4k_{\rm I}k_{\rm II} e^{ik_{\rm II}d}}{(k_{\rm I} -k_{\rm II})(k_{\rm II} -k_{\rm III}) e^{2ik_{\rm II}d} +(k_{\rm I} +k_{\rm II})(k_{\rm II} +k_{\rm III})}
$$

ただし、各領域の波数は

$$
k_{\rm I} =\frac{\sqrt{2m(E-V_{\rm I})}}{\hbar} \ , \ k_{\rm II} =\frac{\sqrt{2m(E-V_{\rm II})}}{\hbar}\ , \ k_{\rm III} =\frac{\sqrt{2m(E-V_{\rm III})}}{\hbar}
$$

で与えられるよ。上記の係数を用いて各領域の波動関数を得ることもできるのだけれども、任意の層状構造に適用できる一般的な方法である「転送行列法」もあとで紹介するよ。

平面波アニメーション

上図は照射した際の領域Ⅰ、領域Ⅱ、領域Ⅲの実部（青）、虚部（橙）と絶対値（緑）をそれぞれ描画したスナップショットだよ。入射波のエネルギーを$${E=1.0[{\rm eV}]}$$とし、箱型ポテンシャルエネルギーの大きさを$${V_{\rm II}=0.8[{\rm eV}]}$$と与えているよ。この場合、$${E>V_{\rm II}}$$のため領域Ⅱの波数は実数となるのだけれども、もし、$${E<V_{\rm II}}$$となるポテンシャルエネルギーを与えると、領域Ⅱの波数は虚数となるよ。

【量子力学】凸型ポテンシャル障壁へ照射した平面波の時間発展 #量子力学 #python #散乱問題 pic.twitter.com/btwHQPAEtW
— シミュ君 (@simu_kun) January 20, 2025

プログラムソース：Python

上記アニメーションを生成するPythonプログラムソースを公開するよ。もしよければ、試してみてくださーい。

#################################################################
## 箱型ポテンシャル障壁へ照射した平面波の時間発展（E = 1.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#虚数単位
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19

######################################
#　物理系の設定
######################################
#電子のエネルギー
E = 1.0 * eV
#角振動数
omega = E/hbar
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 100

#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 2.0 * nm
#壁の高さ
V1 = 0 * eV
V2 = 0.8 * eV
V3 = 0.0 * eV

#各領域の波数ベクトル
k1 = np.sqrt(np.complex128( 2.0 * me * (E - V1) / hbar**2 ))
k2 = np.sqrt(np.complex128( 2.0 * me * (E - V2) / hbar**2 ))
k3 = np.sqrt(np.complex128( 2.0 * me * (E - V3) / hbar**2 ))

#描画範囲
x_min = -4.0 * nm
x_max =  6.0 * nm
#描画区間数
NX = 500
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, d,     NX)
x3 = np.linspace(d, x_max, NX)

######################################
# 各種係数を計算する関数定義
######################################
def calculateCoefficient( k1, k2, k3 ):
	#分母
	C = (k1 - k2) * (k2 - k3) * np.exp( 2.0 * I * k2 * d ) + (k1 + k2) * (k2 + k3)
	#係数
	A = 2.0 * k1 * (k2 + k3) / C
	B = 2.0 * k1 * (k2 - k3) * np.exp( 2.0 * I * k2 * d ) / C
	rc = ( (k1 + k2) * (k2 - k3) * np.exp( 2.0 * I * k2 * d ) + (k1 - k2) * (k2 + k3) )/ C
	tc = 4.0 * k1 * k2 * np.exp( I * k2 * d ) / C
	return A, B, rc, tc

#係数の取得
A, B, rc, tc = calculateCoefficient( k1, k2, k3 )
#################################################################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画
#################################################################

#アニメーション作成用
imgs=[]
	
#各時刻における計算
for tn in range(NT):
	t = dt * tn
	#波動関数の計算
	#領域１
	phi_I =      np.exp(   I * k1 * x1 - I * omega * t )
	phi_R = rc * np.exp( - I * k1 * x1 - I * omega * t )
	phi1 = phi_I + phi_R
	#領域２
	phi2_p = A * np.exp(   I * k2 * x2 - I * omega * t )
	phi2_m = B * np.exp( - I * k2 * x2 - I * omega * t )
	phi2 = phi2_p + phi2_m 
	#領域３
	phi3 = tc * np.exp( I * k3 * (x3 - d) - I * omega * t )

	#各コマを描画
	img  = plt.plot(x1/nm, phi1.real, colors[0], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x1/nm, phi1.imag, colors[1], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x1/nm, abs(phi1), colors[2], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x2/nm, phi2.real, colors[0], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x2/nm, phi2.imag, colors[1], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x2/nm, abs(phi2), colors[2], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x3/nm, phi3.real, colors[0], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x3/nm, phi3.imag, colors[1], linestyle='solid', linewidth = 5.0)
	img += plt.plot(x3/nm, abs(phi3), colors[2], linestyle='solid', linewidth = 5.0)
	#壁の描画
	img += plt.plot([-4, 0, 0, 2, 2, 6], [0, 0, 2*0.8, 2*0.8, 0, 0], color ="black", linestyle='solid', linewidth = 5.0)
	#アニメーションに追加
	imgs.append( img )

#グラフの描画（波動関数）
plt.title( u"箱型ポテンシャル障壁へ照射した平面波の時間発展（" + r"$ E = 1.0 [{\rm eV}], d = 2.0[{\rm nm}], V = 0.8[{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ \psi(x, t) $", fontsize=30)

#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.12, top = 0.95)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)


#描画範囲を設定
plt.xlim([x_min/nm, x_max/nm])
plt.ylim([-2.0, 2.5])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, imgs, interval=30)

#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()

【Python】凸型ポテンシャル障壁へ照射した平面波の時間発展（量子力学）

平面波の解析解

平面波アニメーション

プログラムソース：Python

いいなと思ったら応援しよう！