【Python】凸型ポテンシャル障壁に向けて照射したガウス波束の時間発展（量子力学）
先に得られた平面波は波数$${k}$$が異なっても元のシュレディンガー方程式を満たすので、箱型ポテンシャル障壁が存在する場合でも、平面波の重ね合わせで波束の運動を調べることができるね。異なる波数（角振動数）をもつ平面波を次のとおりに表わすことができるよ。
$$
\psi_{\rm I}(x,t) = \frac{1}{\sqrt{L}} \sum\limits_{n}c_n \left[e^{ik_{{\rm I}n}x} +r_n e^{-ik_{{\rm I}n}x}\right] e^{-i\omega_n t}
$$
$$
\psi_{\rm II}(x,t) = \frac{1}{\sqrt{L}} \sum\limits_{n}c_n \left[A_ne^{ik_{{\rm II}n}x} +B_ne^{-ik_{{\rm II}n}x}\right] e^{-i\omega_n t}
$$
$$
\psi_{\rm III}(x,t) = \frac{1}{\sqrt{L}} \sum\limits_{n}c_n \left[t_ne^{ik_{{\rm III}n}(x-d)} \right] e^{-i\omega_n t}
$$
ただし、エネルギーと角振動数、波数の関係性（分散関係）は以下のとおりだよ。
$$
E_n=\hbar\omega_n=\frac{\hbar^2k_{n}^2}{2m}
$$
$$
k_{{\rm I}n} =\frac{\sqrt{2m(E_n-V_{{\rm I}})}}{\hbar} \ ,\ \ k_{{\rm II}n} =\frac{\sqrt{2m(E_n-V_{{\rm II}})}}{\hbar}\ ,\ \ k_{{\rm III}n} =\frac{\sqrt{2m(E_n-V_{\rm III})}}{\hbar}
$$
ガウス波束は波数の重ね合わせの重みをガウス分布（正規分布）とするよ。
$$
c_n = c_0\exp\left[ -\left(\frac{ k_n - \bar{k} }{2\sigma} \right)^2 \right]
$$
ガウス波束アニメーション
下図は中心エネルギー$${E=20[{\rm eV}]}$$のガウス波束を凸型ポテンシャル障壁（$${V=10[{\rm eV}]}$$）に照射したときのスナップショットだよ。３つの線は各領域の波動関数の実部（青）と虚部（橙）と絶対値（緑）を表しています。ガウス波束は広がりながら進んでいき、障壁に衝突した後に反射波と透過波が生じるね。障壁の中に閉じ込められている様子もわかるよ。
【量子力学】凸型ポテンシャル障壁へ照射したガウス波束の時間発展（V=20[eV]） pic.twitter.com/oOwyTxkB8e
— シミュ君 (@simu_kun) January 29, 2025
プログラムソース：Python
上記アニメーションを生成するPythonプログラムソースを公開するよ。もしよければ、試してみてくださーい。
#################################################################
## 凸型ポテンシャル障壁へ照射したガウス波束の時間発展
#################################################################
import os
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
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20250129/"
if( OUTPUT_FLAG and os.path.isdir(OUTPUT_PASS) == False) : os.mkdir(OUTPUT_PASS)
######################################
# 物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19
######################################
# 物理系の設定
######################################
#波束の中心エネルギー
E0 = 20.0 * eV
#重ね合わせの数
NK = 2000
#ガウス分布の幅
sigma = 2 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
#計算時間の幅
ts = 0
te = 250
#時間間隔
dt = 0.4 * 10**-16
#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 2.0 * nm
#壁の高さ
V1 = 0 * eV
V2 = 20.0 * eV
V3 = 0.0 * eV
#ガウス波束の初期位置
x0 = -5.0 * nm
#描画範囲
x_min = -5.0 * nm
x_max = 5.0 * nm
#描画区間数
NX = 100
#座標点配列の生成
x1 = np.linspace(x_min, 0, int(abs(x_min/nm)*NX))
x2 = np.linspace(0, d, int(abs(d/nm)*NX))
x3 = np.linspace(d, x_max, int(abs((x_max-d)/nm)*NX))
######################################
# 各種係数を計算する関数定義
######################################
def calculateCoefficient( k1, k2, k3 ):
#分母
C = (k1 - k2) * (k2 - k3) * np.exp( 2.0 * I * k2 * d ) + (k1 + k2) * (k2 + k3)
if( C != 0 ):
#係数
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
else:
return 0, 0, 1, 0
#ガウス波束の値を計算する関数
def Psi( x, t, region ):
#波動関数値の初期化
psi = x * (0.0 + 0.0j)
#各波数ごとの寄与を足し合わせる
for kn in range(NK):
#各波数を取得
k = k0 + dk * (kn - NK/2)
E = (hbar * k)**2 / (2.0 * me)
#波数から各振動数を取得
omega = hbar / (2.0 * me) * k**2
#各領域の波数
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 ))
#係数の取得
A, B, rc, tc = calculateCoefficient( k1, k2, k3 )
#ガウス分布
Ck = np.exp( - I * k1 * x0 ) * np.exp( -( (k - k0) / (2.0 * sigma) )**2 )
#平面波を足し合わせる
if(region == 1):
#入射波と反射波と透過波
psi_I = np.exp( I * k1 * x - I * omega * t )
psi_R = rc * np.exp( - I * k1 * x - I * omega * t )
psi += (psi_I + psi_R) * Ck
elif(region == 2):
psi2_p = A * np.exp( I * k2 * x - I * omega * t )
psi2_m = B * np.exp( - I * k2 * x - I * omega * t )
psi += (psi2_p + psi2_m) * Ck
elif(region == 3):
psi_T = tc * np.exp( I * k3 * ( x - d ) - I * omega * t )
psi += psi_T * Ck
return psi
#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0, 1 ))
######################################
# 波動関数の計算
######################################
#アニメーション作成用
imgs=[]
#各時刻における計算
for tn in range(ts, te + 1):
print(tn)
t = dt * tn
#波動関数の計算
psi1 = Psi( x1, t, 1 ) / psi_abs
psi2 = Psi( x2, t, 2 ) / psi_abs
psi3 = Psi( x3, t, 3 ) / psi_abs
#各コマを描画
img = plt.plot(x1/nm, psi1.real, colors[0], linestyle='solid', linewidth = 5.0)
img += plt.plot(x1/nm, psi1.imag, colors[1], linestyle='solid', linewidth = 5.0)
img += plt.plot(x1/nm, abs(psi1), colors[2], linestyle='solid', linewidth = 5.0)
img += plt.plot(x2/nm, psi2.real, colors[0], linestyle='solid', linewidth = 5.0)
img += plt.plot(x2/nm, psi2.imag, colors[1], linestyle='solid', linewidth = 5.0)
img += plt.plot(x2/nm, abs(psi2), colors[2], linestyle='solid', linewidth = 5.0)
img += plt.plot(x3/nm, psi3.real, colors[0], linestyle='solid', linewidth = 5.0)
img += plt.plot(x3/nm, psi3.imag, colors[1], linestyle='solid', linewidth = 5.0)
img += plt.plot(x3/nm, abs(psi3), colors[2], linestyle='solid', linewidth = 5.0)
#壁の描画
img += plt.plot([-5, 0, 0, 2, 2, 6], [0, 0, 1, 1, 0, 0], color ="black", linestyle='solid', linewidth = 5.0)
imgs.append( img )
if( OUTPUT_FLAG ):
path = OUTPUT_PASS + str( tn ) + '.txt'
f = open(path, 'w')
_x1 = x1/nm
_x2 = x2/nm
_x3 = x3/nm
for i in range(len(_x1)):
x = _x1[ i ]
psi = psi1[ i ]
f.write( str( round( x, 3)) + "\t")
f.write( str( round( psi.real, 3)) + "\t")
f.write( str( round( psi.imag, 3)) + "\t")
f.write("\n")
for i in range(len(_x2)):
x = _x2[ i ]
psi = psi2[ i ]
f.write( str( round( x, 3)) + "\t")
f.write( str( round( psi.real, 3)) + "\t")
f.write( str( round( psi.imag, 3)) + "\t")
f.write("\n")
for i in range(len(_x3)):
x = _x3[ i ]
psi = psi3[ i ]
f.write( str( round( x, 3)) + "\t")
f.write( str( round( psi.real, 3)) + "\t")
f.write( str( round( psi.imag, 3)) + "\t")
f.write("\n")
f.close()
#グラフの描画（波動関数）
plt.title( u"凸型ポテンシャル障壁へ照射したガウス波束の時間発展（" + r"$ E_0 = 20.0 [{\rm eV}], d = 2.0[{\rm nm}], V = 20.0[{\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([-1.2, 1.2])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, imgs, interval=30)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output2.mp4", writer="ffmpeg", dpi=300)
#################################################################
## ガウス分布の描画
#################################################################
fig_gaussian = plt.figure(figsize=(15, 8))
#波数の計算範囲
k_min = k0 - dk * NK / 2.0
k_max = k0 + dk * NK / 2.0
#波数座標点配列の生成
kl = np.linspace(k_min, k_max, NK)
El = (hbar * kl)**2 / (2.0 * me) / eV
#正規分布
c_n = np.exp( -((kl - k0) / (2.0 * sigma))**2 )
#グラフの描画（固有関数）
plt.title( u"ガウス波束のエネルギー分布", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ E[{\rm eV}] $", fontsize=30)
plt.ylabel(r"$ c_n/c_0 $", fontsize=30)
#罫線の描画
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([8.5, 36])
plt.ylim([0, 1.05])
#波数分布グラフの描画
plt.plot(El, c_n, linestyle='solid', linewidth = 5)
#グラフの表示
plt.show()