[直観的に量子力学入門、電子の存在確率を可視化] 時間依存1次元シュレーディンガー方程式をルンゲクッタ法で解いて坂道ポテンシャル中での電子波束の挙動をシミュレーション。

#0. 初めに#
2021/9/18 初投稿
電子などのミクロな世界の粒子は波としての性質を持っておりこのことは二重スリット実験によって確かめられている。量子力学ではすべての状態それぞれある割合で重なり合っていると解釈する。それゆえに粒子のいる位置が確率的にしか定まらないし、その割合が確率に相当する。またこの波の振幅の絶対値¹の2乗²値がその場所に電子が存在する確率密度³に相当する。この波は波動関数と呼ばれており、シュレーディンガー方程式を解くことによって求まる。
この記事では波動関数を平面波の有限和で近似しその展開係数の時間変化を計算することで時間依存1次元シュレーディンガー方程式を数値的に解き、存在確率分布の時間変化を可視化する。

シュレーディンガー方程式(1次元時間依存)

i\hbar\frac{\partial\psi(x,t)}{\partial t} = \hat{H}\psi(x,t) 

###リンク###
1.物理的イメージと問題設定
2.解くべき方程式
3.方針
4.計算方法とアルゴリズム
5.計算コード
6.波束アニメーション

#1.物理的イメージと問題設定##
##波動関数の定性的なイメージ##

赤線:\psi(x) \space\space\space
青線: |\psi(x)|^2

波の振幅が大きいところほど発見されやすい

####例えば粒子の位置を10回位置を測定し、3か所で観測されたとき⁴。####

この例では$x=x_1$で2回、$x=x_2$で5回、$x=x_3$で3回測定された。
上の図の場合波動関数の値は

|\psi(x_1)|^2 = 0.2 \space\space\space\space |\psi(x_2)|^2 = 0.5 \space\space\space\space |\psi(x_2)|^2 = 0.3

位置	存在確率
$x_1$	0.2 (20%)
$x_2$	0.3 (30%)
$x_3$	0.5 (50%)

この時「$x=x_1$にいる状態」、「$x=x_2$にいる状態」、「$x=x_3$にいる状態」がそれぞれ20%、50%、30%の割合で重なり合っている。割合が実現する確率(存在確率)に相当。

##問題設定##
質量$m$の電子が坂道ポテンシャル中$V(x)=Ax$を1次元運動するとき、波動関数の時間変化を計算する。
このようなポテンシャルは例として一様な電場から受ける力$-qE$によるポテンシャルが挙げられる。

####初期状態t=0(ガウス波束)####

\psi(x,0) = (\pi \sigma^2)^{-\frac{1}{4}} \exp({{-\frac{(x-x_c)^2}{2\sigma}+ik_0(x-x_c)}})\\
x_c: 中心位置\space\space, E_0:中心エネルギー\space\space\\
\sigma: 波束広がりパラメーター, k_0: k_0=\sqrt{\frac{2mE_0}{\hbar^2}}

存在確率が$x_c$を中心に正規分布になっており、$x_c$から離れると存在確率は急速に0に収束する。存在確率が$x_c$付近に集中しており古典的な粒子の描像に近い状態である。⁵

#2.解くべき方程式#

i:虚数単位 \space\space\space\space
\hbar: ディラック定数 \space \space\space\space
m: 粒子の質量 \space \space\space\space
\psi: 波動関数

1次元シュレーディンガー方程式

i\hbar\frac{\partial\psi(x,t)}{\partial t} = \hat{H}\psi(x,t) 
\space \space \space \space
\hat{H}:=-\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x)

#3.方針#
1.波動関数を平面波で展開する(フーリエ級数展開)
2.展開した波動関数をシュレーディンガー方程式に代入する
3.展開係数に関する連立常微分方程式をルンゲクッタ法で解く

#4.計算方法とアルゴリズム#
###定義###

基底関数(平面波):\phi_n(x):=\sqrt{\frac{1}{L}} \exp({ik_nx})\space\space\space波数: k_n:=\frac{2\pi n}{L}\\
ハミルトニアン演算子: \hat{H}:=-\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x)

###波動関数を平面波で展開###

\psi(x,t) = \sum_{i=1}^{n} c_i(t) \phi_i(x)

###シュレーディンガー方程式に代入し、変形###

 i\hbar\sum_{i=1}^{n} \frac{d c_i(t)}{dt} \phi_i(x) = \sum_{i=1}^{n} c_i(t) \hat{H} \phi_i(x)\\

両辺に$\phi_j^*(x)$をかけて$x=-L/2,L/2$で積分する

 i\hbar\sum_{i=1}^{n} \frac{d c_i(t)}{dt} \int_{-\frac{L}{2}}^{\frac{L}{2}} \phi_j^*(x) \phi_i(x) dx = \sum_{i=1}^{n} c_i(t) \int_{-\frac{L}{2}}^{\frac{L}{2}} \phi_j^*(x)  \hat{H} \phi_i(x) dx

###基底関数の規格直交性###

 \int_{-\frac{L}{2}}^{\frac{L}{2}} \phi_i^*(x) \phi_j(x)dx = \delta_{ij}

クロネッカーのデルタ記号

\delta_{ij} := \left\{
\begin{array}{ll}
1 & (i = j) \\
0 & (i \neq j)
\end{array}
\right.

###ハミルトニアンの行列要素###

 h_{ji} := \int_{-\frac{L}{2}}^{\frac{L}{2}} \phi_j^*(x)  \hat{H} \phi_i(x) dx \\

規格直交性を用いると

 h_{ji} = \frac{\hbar^2 k_i^2}{2m}\delta_{ji}+\frac{1}{L} {I_V}_{ji}\\
I_{Vji} := \int_{-\frac{L}{2}}^{\frac{L}{2}} \exp(i(k_i-k_j)x) V(x) dx\\

###ハミルトニアンの行列要素を具体的に計算する###
今回はポテンシャルとして$V(x) = Ax$を用いた。上の式を用いて計算する

I_{Vji} = A\int_{-\frac{L}{2}}^{\frac{L}{2}} x \exp(i(k_i-k_j)x) dx

部分積分

A\int_{-\frac{L}{2}}^{\frac{L}{2}} \frac{d}{dx}( \frac{1}{i(k_i-k_j)} \exp(i(k_i-k_j)x)) x dx\\
=\frac{A}{i(k_i-k_j)}[\exp(i(k_i-k_j)x]_{-\frac{L}{2}}^{\frac{L}{2}}\\
-\frac{A}{i(k_i-k_j)}\int_{-\frac{L}{2}}^{\frac{L}{2}} \exp(i(k_i-k_j)x))dx\\
= \frac{A}{i(k_i-k_j)}\frac{L}{2}(\exp(i(k_i-k_j)L/2) + \exp(-i(k_i-k_j)L/2))\\
+\frac{A}{(k_i-k_j)^2}\frac{L}{2}(\exp(i(k_i-k_j)L/2) - \exp(-i(k_i-k_j)L/2))

Eulerの公式を用いると

\exp(ix) = \cos(x) + i\sin(x)\\
\exp(-ix) = \cos(x) - i\sin(x)

結果

i \neq j: I_{Vji} 
= \frac{AL}{k_i-k_j} \sin{\frac{(k_i-k_j)L}{2}} \\
+ \frac{2A}{(k_i-k_j)^2}\cos{\frac{(k_i-k_j)L}{2}} \\
- \frac{2A}{(k_i-k_j)^2}
\\
i=j: I_{Vji} = 0 \\
h_{ji} = \frac{\hbar^2 k_j^2}{2m}\delta_{ij} + \frac{1}{L} I_{Vji}

###連立常微分方程式###

\frac{d c_j(t)}{dt} =  -\frac{i}{\hbar} \sum_{i=1}^{n} h_{ji} c_i(t) \\

これは行列形式で書くと

\frac{d}{dt}

\begin{pmatrix}
c_{1}(t)\\
c_{2}(t)\\
\vdots
\end{pmatrix}

= -\frac{i}{\hbar} \begin{pmatrix}
h_{1,1} & h_{1,2} & \cdots \\
h_{2,1} & h_{2,2} & \cdots\\
\vdots & \vdots & \ddots
\end{pmatrix}
 
\begin{pmatrix}
c_{1}(t)\\
c_{2}(t)\\
\vdots
\end{pmatrix} 

ルンゲ=クッタ法

$c_i(t)$から$c_i(t+\Delta t)$を求める

\Delta_1:=0, \space\space \Delta_2:=\frac{\Delta t }{2}, \space\space \Delta_3:=\frac{\Delta t }{2},\space\space \Delta_4:=\Delta t\\
\delta_1=\delta_4=1, \space\space \delta_2=\delta_3=2

漸化式

{k_0}_{i} = 0\\
{k_m}_{j} = -\frac{i}{\hbar} \sum_{i=1}^{n} h_{ji} (c_i(t) + \Delta_i {k_i}_{m-1})\\
c_j(t+\Delta t) = c_j(t) + \frac{\Delta t}{6}\sum^{4}_{m=1} \delta_m k_m

#5.計算コード#
###環境###
OS: windows10
Julia: Julia 1.62

####パラメーターなど####

SchrödingerEq.jl

# 2021 8/26
# ディレクトリ移動
cd("Documents\\programming\\Julia\\Quantum-time-evolution\\V=Ax")
# ライブラリ
using LinearAlgebra
using Plots

# スケール
L = 500
small_l = 30
x_lim_min = (-1*small_l/2)
x_lim_max = (small_l/2)
dx = 0.06
# 時間変数
tmin = 0
tmax = 12
dt = 0.001
Nt = Int(round((tmax - tmin)/dt))
# ポテンシャルplotサイズ
pote_small_big = 0.5
#
# 物理定数
A_pa = 1.5
h_bar = 1
m_pa = 1
E0 = 18
k0 = ((2*m_pa*E0)/(h_bar^2))^0.5
xc = -5
σ = 1
nmin = -1000
nmax = 1000
#
# 波数
n_num = nmin:1:nmax
k_i = ((2*π).*n_num)./L 
Nk = size(k_i)[1]

####ハミルトニアンの行列要素####

i \neq j : I_{Vji} 
= \frac{AL}{k_i-k_j} \sin{\frac{(k_i-k_j)L}{2}} \\
+ \frac{2A}{(k_i-k_j)^2}\cos{\frac{(k_i-k_j)L}{2}} \\
- \frac{2A}{(k_i-k_j)^2}
\\
i=j: I_{Vji} = 0 \\
h_{ji} = \frac{\hbar^2 k_j^2}{2m}\delta_{ij} + \frac{1}{L} I_{Vji}

SchrödingerEq.jl

#ハミルトニアン行列
Iv_mat = complex(zeros(Nk, Nk))
AL_i = (A_pa*L)/1im
_2iA = 2*1im*A_pa
_0_5L = (0.5*L)
Iv_jj = 0.0
for j in 1:Nk
    # 波数 
    ki_kj = (k_i .- k_i[j])
    # 非対角成分
    Iv_mat[j, :] = (
        ((AL_i ./ ki_kj) .* cos.(_0_5L .* ki_kj))
        .+ ((_2iA ./ ki_kj.^2) .* sin.(_0_5L .* ki_kj))
    )
    # 対角成分
    Iv_mat[j, j] = Iv_jj
end 
H_mat = (
    ((h_bar^2)/(2*m_pa)) .* diagm(0 => k_i.^2)
    .+ (1/L).*(Iv_mat)
)

####数値計算####

SchrödingerEq.jl

#----------------------配列等定義------------------
# 展開係数
c_lis = complex(zeros(Nk, 2))
_4πσ2_L2_025 = ((4*π*σ^2)/(L^2))^0.25
m_ik0xc = -1im.*k_i.*xc
m_05_σ2 = -0.5*(σ^2)
kn_k02 = (k_i .- k0).^2
#
# 初期条件
c_lis[:, 1] = _4πσ2_L2_025*(
    exp.(m_05_σ2.*kn_k02)
    .* 
    exp.(m_ik0xc)
)
#
# サブデータ(数値解格納用 アニメーションに使う)
sb_dt = 0.05
sb_Nt = Int(round(Nt * dt/sb_dt))
cat_val = Int(round(sb_dt/dt))
sb_c = complex(zeros(Nk, sb_Nt))
sb_i = 2
sb_c[:, 1] = c_lis[:, 1]
#
#
print(" ", 0, "/", Nt)
print(" |c| = ", norm(c_lis[:, 1]), "\n\n")
#
# ルンゲクッタ変数用
m_i_h = (-1im/h_bar)
k_Rk = complex(zeros(5, Nk))
#
#
#------------------------ルンゲクッタ法-------------------
#
@time for i in 2:Nt
    sum_arr = complex(zeros(Nk))
    for j in 2:5
        if j == 2
            delta = 0
        elseif j == 5
            delta = dt
        else
            delta = dt/2
        end
        #
        k_Rk[j, :] = m_i_h .* (H_mat * (c_lis[:, 1] .+ (delta.*k_Rk[j-1, :])))
        #
        if j == 2 || j == 5
            sum_arr .+= k_Rk[j, :]
        else
            sum_arr .+= (2).*k_Rk[j, :]
        end
    end
    #
    c_lis[:, 2] = c_lis[:, 1] .+ ((dt/6).*sum_arr)
    #
    # 初期化
    c_lis[:, 1] = c_lis[:, 2]
    #
    #print(" ", i, "/", Nt)
    #print(" |c| = ", norm(c_lis[:, 2]), "\n\n")
    # データ格納
    if i%cat_val == 0 && sb_i < sb_Nt
        print(" ", i, "/", Nt)
        print(" |c| = ", norm(c_lis[:, 2]), "\n\n")
        sb_c[:, sb_i] = c_lis[:, 2]
        sb_i += 1
    end 
end

####波動関数の空間分布計算####
座標を離散化し基底関数の値を行列にまとめると次のように計算できる

x_i := x_0 + i \Delta x \\

\begin{pmatrix}
\psi(x_1, t)  \\
\psi(x_2, t)  \\
\vdots
\end{pmatrix}
= 
(c_1(t), c_2(t), \cdots)
\begin{pmatrix}
\phi_1(x_1) & \phi_1(x_2) & \cdots \\
\phi_2(x_1) & \phi_2(x_2) & \cdots\\
\vdots & \vdots & \ddots
\end{pmatrix}

SchrödingerEq.jl

print(" bacefunc: φ = e^(ikx)\n")
# 基底関数
x = x_lim_min:dx:x_lim_max
Nx = size(x)[1]
φ = complex(zeros(Nk, Nx))
for i in 1:Nk
   ik = 1im*k_i[i]
   φ[i, :] = (L^-0.5).*exp.(ik.*x)
end
#
print(" wavefunc: ψ = Σ c_k * e^(ikx)\n")
# 波動関数
ψ = complex(zeros(Nx, sb_Nt))
for i in 1:sb_Nt
    ψ[:, i] = transpose(sb_c[:, i]) * φ
end

####可視化 、gif保存####

SchrödingerEq.jl

# ポテンシャル
psi2_min = findmin((abs.(ψ)).^2)[1]
psi2_max = findmax((abs.(ψ)).^2)[1]
V_x = A_pa.*x
V_x = (pote_small_big * (psi2_max/findmax(abs.(V_x))[1])) .* V_x
# 波動関数をポテンシャルの上に乗せる
plot_wave = zeros(Nx, sb_Nt)
for i in 1:sb_Nt
   plot_wave[:, i] = (abs.(ψ[:, i])).^2 .+ V_x
end
wave_min = findmin(plot_wave)[1]
wave_max = findmax(plot_wave)[1]
#
#
#
# リム用変数
cmin = findmin((abs.(sb_c)).^2)[1]
cmax = findmax((abs.(sb_c)).^2)[1]
V_x_min = findmin(V_x)[1]
V_x_max = findmax(V_x)[1]
y_lim_min = findmin([wave_min, V_x_min])[1]
y_lim_max = findmax([wave_max, V_x_max])[1]
#
# 波動関数プロット
anim = @animate for i in 1:sb_Nt
    plot_data = zeros(Nx, 2)
    plot_data[:, 1] = plot_wave[:, i]
    plot_data[:, 2] = V_x
    plot(
        x, 
        plot_data,
        xlims = (x_lim_min, x_lim_max),
        ylims = (y_lim_min, y_lim_max),
        xlabel = "x",   #X軸のラベル
        ylabel = "|ψ|^2",   #Y軸のラベル
        title  = "wavefunc",     #タイトル
    )
    print("wave plot: ", i, "/", sb_Nt, "\n\n")
end
gif(anim, "wavefunc2.gif", fps=10)

#6.波束アニメーション#
####アニメーション####

####アニメーション2####

####物理的な解釈####
ポテンシャルから力を受けて加速させられていることが分かります。
最初は電子の存在領域が狭い範囲に集中していたけど、だんだんと波束が広がっていき位置の量子的な不確定性が目立ってきているのが分かる。
また、波束の形が崩れずに形を保っているの波束の中心位置(位置の期待値)は近似的に古典力学に従っていると思われます。
アニメーション2では最初、初期条件として電子に運動エネルギー(正確には期待値で波束の中心エネルギーに相当)を与えたがポテンシャルから受ける力に負けて逆方向に移動していくことが見て取れる。また、1枚目と比べてポテンシャルに逆らっている時の方がより波束が大きく広がり位置の不確定性が目立っている。

1. 波動関数は一般的に複素数の値を取る ↩

2. 波の振幅の二乗値が波のエネルギーに比例するから ↩

3. 区間x=[a,b]で確率密度を積分した値がその区間で粒子が見いだされる確率そのものになる。これを全空間で積分したら1になるように規格化する ↩

4. 例として簡単のため三か所で飛び飛びに発見されたとしているが、普通は連続的な確率分布になる。 ↩

5. あくまでも古典的粒子みたいに完全に位置が決まっているわけではない。 ↩