#*追記
1 May. 2018: トンネル効果(内容(3))ならびに調和振動子ポテンシャル下での運動 (内容(5)&(6))を追加しました。
#はじめに
電子などのミクロな粒子の運動は近似的にシュレディンガー方程式に従うことが知られています。
相互作用ポテンシャル$V(x)$を受ける粒子の波動関数$\psi_(x,t)$に対する時間に依存する一次元シュレディンガー方程式は，
$$i\hbar\frac{\partial}{\partial t} \psi(x,t) = \hat{H} \psi(x,t) $$
$$ \hat{H} = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} +V(x){\tag 1}$$
と表されます。$\psi_(x,t)$は一般に複素数です。時間に依存するシュレディンガー方程式は複素数値関数に対する拡散方程式の形をしていることがわかります。
粒子の時間発展は波動関数を得ることによって知ることができます。時間に依存しない場合に比べてネットでも文献が少ないですが，粒子のダイナミクスや散乱現象を波動関数の時間発展を通じて直接調べることは量子力学の直感的な理解を深めることができると思います。
本稿では，拡散方程式を数値的に解く技術を活用して時間依存シュレディンガー方程式を解くPythonコードを実装します。それを用いて粒子散乱現象を波動関数(複素数)の時間発展から調べます。
本コードは1粒子系に限られていますが，量子力学のエッセンスを学ぶには十分だと思います。
#どんな人向け?
● 量子力学を勉強中の学生さま
● 波動関数の時間発展(量子ダイナミクス)に興味のある方
● 散乱現象に興味がある方
● 波動関数が動く様子を観たい方
#内容
(1)コードを実装し，自由粒子($V(x)=0$)に対して適用します。
(2)次いで強い散乱体を設定し散乱現象を調べます。
(3)トンネル現象をシミュレートします。
(4)弱い散乱体を設定し散乱現象を調べます。
(5)調和振動子ポテンシャル下における定常状態(基底状態)の運動を調べます。
(6)(5)のポテンシャルを変化させたときの基底状態の変化を調べます。これは量子力学における時間に依存する摂動論の理解に役立ちます。
初期条件としてガウス関数(波束)を考えます。
十分遠方での波動関数の値(境界条件)を$\psi =0$とします。
空間・時間解像度はほどほどの条件にしています。興味がある方はより解像度を上げたシミュレーションをしてみてください。
#方法
空間のグリッド幅を$\Delta x$として，波動関数\psi(x,t) を等分点での値($n=0,1,2,...$)$$\psi_n(t)=\psi(n\Delta x, t){\tag 2}$$
で置き換えます。
このときシュレディンガー方程式は以下の差分方程式にしたがいます。
$$i\hbar \frac{\partial}{\partial t}\psi_n(t) = \sum_m H_{nm}\psi_m(t) $$
$$H_{nm} = -\frac{\hbar^2}{2m (\Delta x)^2}(\delta_{n,m-1}-2\delta_{n,m}+\delta_{n,m+1}+V(n\Delta_x)\delta_{n,m}$$
$\delta_{n,m}$はクロネッカーのデルタを表します。
系の時間発展にはクランク-ニコルソンの陰的スキーム(補遺)を，空間に対しては中心差分法を用います。
すると解くべき行列方程式として
$$(I+\frac{i\Delta t}{2\hbar} H), \psi(t+\Delta t)=(I-\frac{i\Delta t}{2\hbar} H), \psi(t) $$
が得られます。ここで$I$は単位行列，$H$はハミルトニアン行列， $\psi$はベクトル ($\psi_0, \psi_1, ...$)を表します。
与えられた初期条件$\psi(0)$から出発しこの連立方程式を繰り返し解くことで波動関数の時間発展を計算できます。
$Aを=(I+\frac{i\Delta t}{2\hbar} H)$として，波動関数の時間発展を
$$\psi(t+\Delta t) = A^{-1}A^\dagger\psi(t)$$
として求めます。なお，行列$A$の性質を考えるともう少しエレガントに計算できますが，わかりやすさを重視してこの形で計算しています。
#(1) 自由粒子の伝播
#計算コード
注意: jupyter-notebookを利用しています。また，高速計算を主眼としていませんので，ハイパフォーマンスコンピューティングの観点からは非効率的な記述箇所が見受けられると思います。
"""
時間に依存するシュレディンガー方程式
自由粒子: 波束
1 May. 2018
"""
%matplotlib notebook
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
mass=1/2 # mass of electron
hbar=1
Nx =800 # x方向のグリッド点数
Nt =400 # t 方向のグリッド点数
Lt =100
x0, x1 = -200, 200 #境界のx座標
delta_x=(x1-x0)/Nx
delta_t=Lt/Nt
print("delta_x=",delta_x)
print("delta_t=",delta_t)
r=delta_t/(delta_x**2)
print("r=",r)
#
psi = np.zeros([Nx,Nt],dtype=np.complex) # 波動関数。複素数。
# 初期条件
width=5
xc= x0+delta_x*(0.1*Nx)
p0=2 # 初期運動量
for n in range(Nx):
xx = x0+n*delta_x
psi [n,0] = ((np.pi**(-1/4))/np.sqrt(width))*np.exp(-((xx-xc)**2)/(2*width**2))*np.exp(1j*p0*(xx-xc)) # 初期条件
#関数定義
def krondel(nn,mm):
if nn==mm:
return 1
else:
return 0
#
I_mat=np.identity(Nx, dtype=float) #単位行列
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
for n in range(Nx):
V[n] = 0
## ハミルトニアン行列の作成
Hmat=np.zeros([Nx,Nx])
for n in range(Nx):
for m in range(Nx):
Hmat[n,m] = (-hbar**2/(2*mass*delta_x**2))*(krondel(n,m-1)-2*krondel(n,m)+krondel(n,m+1))+V[n]*krondel(n,m)
## 係数行列の作成
Amat = I_mat+(1j*delta_t/(2*hbar))*Hmat
#AAmat = np.conjugate(Amat.T) # Aのエルミート共役: A^dagger
AAmat = I_mat-(1j*delta_t/(2*hbar))*Hmat
A_invmat=np.linalg.inv(Amat) #Aの逆行列
# メイン
psi_vec_t =np.zeros([Nx],dtype=np.complex)
dum_psi =np.zeros([Nx],dtype=np.complex)
# 境界条件 (境界でψ=0)
psi [0,0] = 0
psi [-1,0] = 0
for jj in range(1,Nt):
psi_vec_t = psi[:,jj-1]
dumA=np.dot(A_invmat,AAmat)
psi[:,jj] = np.dot(dumA, psi_vec_t)
psi [0,jj] = 0 # 境界条件 (境界でψ=0)
psi [-1,jj] = 0
psi_Real = np.real(psi) # 波動関数の実部
psi_Imag= np.imag(psi) # 波動関数の虚部
# 確率密度の計算
elec_density =np.absolute(psi)**2
##結果(1) 自由粒子の伝播
波動関数の実部・虚部，そしてそれらの積である確率密度の時間発展を可視化します。jupyter-notebook環境上で実行しました。
%matplotlib nbagg
from matplotlib.animation import ArtistAnimation # アニメーション作成のためのメソッドをインポート
fig = plt.figure()
anim = [] #アニメーション用に描くパラパラ図のデータを格納するためのリスト
for i in range(Nt):
wavefunc_Real=list(psi_Real[:,i])
wavefunc_Imag=list(psi_Imag[:,i])
density=list(elec_density[:,i])
x = np.arange(x0, x1, delta_x)
if i % int(Nt*0.04) ==0:
im=plt.plot(x, wavefunc_Real, '-', color='blue', linewidth =0.1, aa=True)
im2=plt.plot(x, wavefunc_Imag, '-', color='red',markersize=10, linewidth =0.1, aa=True)
im3=plt.plot(x,density, '-', color='black',markersize=10, linewidth =0.4, aa=True)
anim.append(im+im2+im3)
anim = ArtistAnimation(fig, anim) # アニメーション作成
plt.xlabel('x')
plt.ylabel('')
plt.xlim(-200, 200)
plt.ylim(-0.3, 0.3)
plt.plot(x,0.02*V, '--', color='black',markersize=10, linewidth =1, aa=True)
plt.legend(('Psi_Real', 'Psi_Imag', 'Density (x4)'), loc='upper left')
fig.show()
#anim.save("t.gif", writer='imagemagick') #アニメーションをt.gifという名前で保存し，gifアニメーションファイルを作成する。
自由波束の伝播の様子。青色が波動関数の実部。赤色が虚部。黒線が確率密度， 黒点線がポテンシャル障壁(見やすくするために0.02倍しています)を表しています。時間が経つにつれ少しづつ波形が崩れていきますが，波束の中心は等速運動をしています。
これは解析解から得られる挙動と一致しています。
#(2) 散乱現象 強い散乱 (入射粒子のエネルギー << 障壁高さ)
次に，波束と散乱体との相互作用を考えます。ポテンシャル障壁を設定し，波束がどのように散乱されるのか調べます。
それには(1)の計算コード中で相互作用ポテンシャルVの定義を
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
i0=int(0.8*Nx)
i1=int( 0.9*Nx)
for n in range(i0,i1): #散乱体を置く
V[n] =50
などと変更します。
この例では，空間に$0.1Nxdelta_x$の幅の散乱体(高さ50)を設定しています。これは粒子のエネルギーよりも大きいです。
##結果(2) 散乱現象1 強い散乱 (入射粒子のエネルギー << 障壁高さ)
図. 強い散乱の様子。青色が波動関数の実部。赤色が虚部。黒線が確率密度， 黒点線がポテンシャル障壁(見やすくするために0.02倍しています)を表します。
散乱体により粒子が完全に散乱されている様子がわかります。
#(3) トンネル現象 (入射粒子のエネルギー < 障壁高さ)
入射粒子のエネルギーよりも障壁の高さが"ほどほどに"大きい場合は，
エネルギーが障壁高さ以下でも障壁を越えることがあります。これは量子力学固有の現象で，「トンネル効果」と呼ばれていますｌ。この現象をシミュレートしてみます。
障壁の高さを入射粒子のエネルギー ($p^2=16$)よりもわずかに大きな値(17)にします。また，壁の幅を狭くします。
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
i0=int(0.6*Nx)
i1=int( 0.6*Nx+5)
for n in range(i0,i1): #散乱体を置く
V[n] =17
##結果 (3) トンネル現象 (入射粒子のエネルギー < 障壁高さ)
図. トンネル効果.青色が波動関数の実部。赤色が虚部。黒線が確率密度， 黒点線がポテンシャル障壁(見やすくするために0.02倍しています)を表します。
入射粒子のエネルギーが障壁高さよりも小さく，古典力学では粒子は絶対に壁を貫通できませんが，量子力学的粒子ではあり得ます。これがトンネル効果です。
(4) 散乱現象 弱い散乱
次に，散乱障壁を低くしてみます。
(1)のコード中で相互作用ポテンシャルの記述部分を以下のように書き換えます。(2)の場合よりもずっと小さく，入射粒子のエネルギーよりも小さくなっています。
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
i0=int(0.6*Nx)
i1=int( 0.7*Nx)
for n in range(i0,i1): #散乱体を置く
V[n] = 2
結果 (4) 散乱現象 弱い散乱
図. 弱い散乱 (入射粒子のエネルギー > ポテンシャル障壁の高さ )
青色が波動関数の実部。赤色が虚部。黒線が確率密度， 黒点線がポテンシャル障壁(見やすくするために0.02倍しています)を表します。
強い散乱の場合と異なり，ポテンシャル障壁を透過する波束があることがわかります。
なお，古典力学的には反射が生じない状況ですが，量子力学的粒子の持つ波動性に起因して反射(波)が生じます。
そのため，散乱体の左側でも粒子が観測され得ます。量子理学特有の性質です。
この状況がシミュレーション結果でも再現されています。
#(5) 調和振動子ポテンシャル: 定常状態を初期値に選んだ場合
調和振動子ポテンシャル$$V(x)=\frac{m \omega^2x^2}{2}$$を考えます。
初期条件とポテンシャルを以下のようにセットします。
初期条件は定常状態の厳密解です。確率密度が時間に依存しないことを確かめます。
# 初期条件
width=1
xc= 0
p0=0 # 初期運動量
for n in range(Nx):
xx = x0+n*delta_x
psi [n,0] = ((np.pi**(-1/4))/np.sqrt(width))*np.exp(-((xx-xc)**2)/(2*width**2))*np.exp(1j*p0*(xx-xc)) # 初期条件: 波束
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
xc=0
omega=2
for n in range(Nx): #散乱体を置く
xx = x0+delta_x*n
V[n] =mass*omega**2*(xx-xc)**2/2
#結果 (5) 調和振動子ポテンシャル: 定常状態を初期値に選んだ場合
波動関数の実部・虚部は変動しますが，確率密度は変化しません。これは初期状態として定常状態のシュレディンガー方程式の解を選んだためです。
#(6) 調和振動子ポテンシャルの変調に伴う状態変化
(5)で設定した調和振動子ポテンシャル$V(x)=\frac{m \omega^2x^2}{2}$に4次の項を入れ少し変調した$$V'(x)=\frac{m \omega^2x^2}{2}+ x^4$$のもとでの運動を考えます。
初期状態は(5)で用いたものと同じにします。このとき，ポテンシャル$V(x)$から$V'(x)$への変化に伴い状態がどのように変わるのかを調べます。これは量子力学における時間に依存する摂動論の理解に役立ちます。
#相互作用ポテンシャルの定義
V = np.zeros([Nx]) #相互作用ポテンシャル
xc=0
omega=2
for n in range(Nx): #散乱体を置く
xx = x0+delta_x*n
V[n] =mass*omega**2*(xx-xc)**2/2+(xx-xc)**4
##(6) 結果 調和振動子ポテンシャルの変調に伴う状態変化
図 変調ポテンシャル下での波動関数と確率密度(Density)。青色が波動関数の実部，赤色が虚部，黒線が確率密度， そして黒点線がポテンシャル(見やすくするために0.02倍しています)を表します。
(5)調和振動子ポテンシャル下の定常状態が変調され，時間に伴い変化しています。これは新しいポテンシャル下ではもはや定常状態ではなく，他の状態が混ざることを意味しています。他の状態がどの程度混ざるのかはポテンシャルの性質に依ります。ポテンシャルの変化が小さい場合は摂動計算([5])が有効です。
#参考文献
[1] [Pythonによる科学・技術計算] クランク-ニコルソン法(陰解法)とFTCS法(陽解法)による1次元非定常熱伝導方程式の数値解法，放物型偏微分方程式
[2] 逆行列の計算法です [Pythonによる科学・技術計算] 連立一次方程式の解法, 数値計算, numpy
[3] 可視化です [Pythonによる科学・技術計算] 放物運動のアニメーションを軌跡(locus)付きで描画, matplotlib，
こちらも参考になりましたmatplotlibのArtistAnimationで二つ以上のアニメーションを描く方法
シュレディンガー方程式に関する初等的な内容のおさらいには以下の本が手頃です。
[4] 鈴木克彦，シュレディンガー方程式, 共立出版, 2013.
量子散乱現象に関しては以下の本が詳しいです。初等的な量子力学を読んだあとなら読み進められると思います。
[5]砂川重信，散乱の量子論，岩波オンデマンドブックス，2015
補遺
時間発展
ハミルトニアンが時間に依存しない場合，方程式の形式解は演算子の指数関数表示(時間発展演算子)，
$$\hat{U}(t-t_0) = \exp\left(\frac{t-t_0}{i\hbar}\hat{H}\right){\tag 2}$$
を用いて，
$$ \psi(x,t)= \hat{U}(t-t_0)\psi(x,t_0) {\tag 3}$$
と表せます。
ユニタリー性。確率保存。
粒子はどこかにいます。量子力学ではそれが以下のように表されます。
$$\int_{-\infty}^{\infty} |\psi(x,t)|^2dx=1$$
これに対応して時間発展演算子$U$が
$$UU^\dagger =1 $$
という性質を持っていることが分かります。この性質を持つ演算子をユニタリー演算子と言います。
時間に対する単純な差分スキームの問題点
時間微分は前進差分から
$$\psi(x,t+\Delta t)\approx \frac{1}{i\hbar}\hat{H} \psi(x,t) \Delta t+\psi(x,t)=(1+\frac{1}{i\hbar}\hat{H}\Delta t)\psi(x,t)$$
となります。これは時間発展演算子(2)を$\Delta t= t-t0$として$\Delta t$の1次まで展開したことに対応します。
このとき
$$U^\dagger U= 1+\frac{1}{\hbar ^2}\hat{H^2} (\Delta t)^2 \gt 1$$
であり，ユニタリー性(確率保存性)は満たされません。
$$\int_{-\infty}^{\infty} |\psi(x,t)|^2dx \neq 1$$
また，$U^\dagger U=1$から$U^{-1}=U^\dagger$となることを利用すると
$$U^\dagger \psi(x,t+\Delta t) = \psi(x,t)$$
が得られます。この式を$\Delta t$の1次まで解くと，形式的な解として
$$ \psi(x,t+\Delta t) \approx \frac{1}{1+\frac{1}{i\hbar}\hat{H}\Delta t}\psi(x,t)$$
が得られます。これは陰的差分となっています。
この場合も$U^\dagger U \lt 1$となりユニタリー性が崩れ粒子数が保存されません。
#ケーリー(Cayley)形式による時間発展
クランク-ニコルソン法の差分スキームを用いて波動関数の時間発展を計算します。時刻$t+\Delta t/2)$における関数評価を用いて
$$\psi(x,t+\Delta t/2) \approx \frac{\psi(x,t+\Delta t)+\psi(x,t)}{2}$$
とし，上述の陽的・陰的差分スキームによる解を用いて，
$$\psi(x,t+\Delta t) \approx \frac{1-\frac{1}{i\hbar}\hat{H}\frac{1}{2}\Delta t}{1+\frac{1}{i\hbar}\hat{H} \frac{1}{2}\Delta t}\psi(x,t) $$
となります。これをケーリー(Cayley)形式といいます。時間差分について二次精度になっています。
注目すべきは，右辺の$\frac{1-\frac{1}{i\hbar}\hat{H}\frac{1}{2}\Delta t}{1+\frac{1}{i\hbar}\hat{H} \frac{1}{2}\Delta t}$が近似的に時間発展演算子とみなせて，この演算子はユニタリー性を満たしています。したがって，この形式を使うことで，ユニタリー性が厳密に保存され，粒子数が保存されます。また，エネルギー等の物理量も厳密に保存されます。これらの理由から時間に依存するシュレディンガー方程式の時間発展の差分スキームとしてケーリー形式が適していることが分かります。
(参考)ハミルトニアンが時間に依存する場合(たとえば，系が外場である電磁場と相互作用するとき)の時間発展演算子
\begin{align}
\hat{U}(t-t_0) &= I + \frac{1}{i\hbar} \int_{t_0}^{t} \hat{H}(t_1) dt_1
+ \cdots
\\&+ \left(\frac{1}{i\hbar}\right)^n \int_{t_0}^{t}\cdots \int_{t_0}^{t_{n-1}}
\hat{H}(t_1) \cdots \hat{H}(t_n) dt_1 \cdots dt_n
+ \cdots\,.
\end{align}
となります。