【Python】ステップ関数型ポテンシャル障壁に向けて照射したガウス波束の時間発展（量子力学）
先に得られた平面波は波数$${k}$$が異なっても元のシュレディンガー方程式を満たすため、ステップ関数型ポテンシャル障壁が存在する場合でも、平面波の重ね合わせで波束の運動を調べることができるよ。異なる波数（角振動数）をもつ平面波を次のとおりに表すことができるね。
$$
\psi_{\rm I}(x,t)= \frac{1}{\sqrt{L}} \sum\limits_{n}c_n \left[e^{ik_{{\rm I}n}x}+r(k_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[t(k_n)e^{ik_{{\rm II}n}x}\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}
$$
ガウス波束アニメーション
上図は照射したガウス波束が障壁に衝突したときのスナップショットだよ。３つの線は波動関数の実部（青）と虚部（橙）と絶対値（緑）を表しているよ。次に示す動画のとおり、ガウス波束は広がりながら進んでいき、障壁に衝突した後に反射波と透過波が生じているね。
入射波の中心エネルギーとポテンシャルエネルギーが等しい場合（$${E=V_0}$$）、透過後の波束の運動が面白いね。ピークの位置は少しずつ移動するけれども、波束の左端は境界面に固定されているね。
プログラムソース：Python
上記アニメーションを生成するPythonプログラムソースを公開するよ。もしよければ、試してみてくださーい。
#################################################################
## ステップ関数型ポテンシャル障壁へ照射したガウス波束の時間発展（E = 30.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
#################################################################
## 物理系の設定
#################################################################
#波束の中心エネルギー
E0 = 30.0 * eV
#重ね合わせの数
NK = 2000
#ガウス分布の幅
sigma = 2 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
#空間分割サイズ
dx = 1.0 * 10**-9
#空間分割数
NX = 500
#描画範囲
x_min = -10.0 * dx
x_max = 10.0 * dx
#壁の高さ
V = 30.0 * eV
#ガウス波束の初期位置
x0 = -7.5 * dx
#計算時間の幅
ts = 0
te = 500
#時間間隔
dt = 0.2 * 10**-16
#ガウス波束の値を計算する関数
def Psi( x, t, region ):
#波動関数値の初期化
psi = x * (0.0 + 0.0j)
#各波数ごとの寄与を足し合わせる
for kn in range(NK):
#各波数を取得
k1 = k0 + dk * (kn - NK/2)
#波数から各振動数を取得
omega = hbar / (2.0 * me) * k1**2
#平面波のエネルギー
E = (hbar * k1)**2 / (2.0 * me)
#領域IIの波数
k2 = np.sqrt(np.complex128( 2.0 * me * (E - V) / hbar**2 ))
#反射係数
rc = ( k1 - k2 ) / ( k1 + k2 )
#透過係数
tc = 2.0 * k1 / ( k1 + k2 )
#入射波と反射波と透過波
phi_I = np.exp( I * k1 * x - I * omega * t )
phi_R = rc * np.exp( - I * k1 * x - I * omega * t )
phi_T = tc * np.exp( I * k2 * x - I * omega * t )
#ガウス分布
Ck = np.exp( - I * k1 * x0 ) * np.exp( -( (k1 - k0) / (2.0 * sigma) )**2 )
#平面波を足し合わせる
if(region == 1):
psi += (phi_I + phi_R) * Ck
elif(region == 2):
psi += phi_T * Ck
return psi
#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0, 1 ))
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, x_max, NX)
#アニメーション作成用
imgs = []
plt.text(-9.9, 1.28, r"$\psi_{\rm I}(x,t)= \frac{1}{\sqrt{L}} \sum_{n}c_n \left[e^{ik_{{\rm I}n}x}+r(k_n) e^{-ik_{{\rm I}n}x}\right] e^{-i\omega_n t} $" , ha='left', va='center', fontsize=25)
plt.text(1.5, -0.85, r"$\psi_{\rm II}(x,t)= \frac{1}{\sqrt{L}} \sum_{n}c_n \left[t(k_n)e^{ik_{{\rm II}n}x}\right] e^{-i\omega_n t} $" , ha='left', va='center', fontsize=25)
### 波束の空間分布の計算
for tn in range(ts, te + 1):
print(tn)
#実時間の取得
t = dt * tn
#時刻tの波動関数を取得
psi1 = Psi( x1, t, 1 )
psi2 = Psi( x2, t, 2 )
#波動関数の規格化
psi1 /= psi_abs
psi2 /= psi_abs
#各コマを描画
img = plt.plot(x1/dx, psi1.real, colors[0], linestyle='solid', linewidth = 3.0)
img += plt.plot(x1/dx, psi1.imag, colors[1], linestyle='solid', linewidth = 3.0)
img += plt.plot(x1/dx, abs(psi1), colors[2], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, psi2.real, colors[3], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, psi2.imag, colors[4], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, abs(psi2), colors[5], linestyle='solid', linewidth = 3.0)
#壁の描画
img += plt.plot([-10, 0, 0, 10], [0, 0, 1.0, 1.0], color ="black", linestyle='solid', linewidth = 5.0)
#time = plt.text( 8.5, 1.57, "t = " + str(round(t/10**-15, 1)) + r"$[{\rm fs}]$" , ha='left', va='center', fontsize=18)
#テキストをグラフに追加
#img.append(time)
#アニメーションに追加
imgs.append(img)
#グラフの描画（波動関数）
plt.title( u"ステップ関数型ポテンシャル障壁に向けて照射したガウス波束の波動関数（" + r"$ E_0 = 30.0 [{\rm eV}] , V = 30.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.085, right = 0.97, bottom=0.15, 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/dx, x_max/dx])
plt.ylim([-1.1, 1.5])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, imgs, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.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
#波数座標点配列の生成
k = np.linspace(k_min, k_max, NK)
#正規分布
c_n = np.exp( -((k - k0) / (2.0 * sigma))**2 )
#グラフの描画（固有関数）
plt.title( u"ガウス分布（正規分布）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ (k_n - \bar{k})/\sigma $", 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([(k_min- k0)/sigma, (k_max- k0)/sigma])
plt.ylim([0, 1.05])
#波数分布グラフの描画
plt.plot((k-k0)/sigma, c_n, linestyle='solid', linewidth = 5)
#グラフの表示
#plt.show()