[Pythonによる科学・技術計算]行列法による3次元等方調和振動子ポテンシャル中の定常状態のシュレディンガー方程式の解法，境界値問題，量子力学

#はじめに

**定常状態のシュレディンガー方程式は，差分法により行列の(一般)固有値問題に帰着させることが可能である。
本稿ではその方法を用いて3次元等方調和振動子ポテンシャル中の電子の固有エネルギーおよび固有関数を決定することを目的とする。**


行列法のあらましについては

[[Pythonによる科学・技術計算]行列形式による常微分方程式の境界値問題の解法，数値計算](http://qiita.com/sci_Haru/items/2d238fcadbe7990239df)[1]

を参考にしていただければ幸いです。


***
#内容
3次元等方調和振動子ポテンシャル
$V(r) = \frac{m_e \omega^2r^2}{2} {\tag {1}}$
に対する動径方向の波動関数$R(r)$に対する定常状態のシュレディンガー方程式は，

$u(r) = r R(r) {\tag{2}}$

として，

$-u''(r)+\frac{2 m_e}{\hbar^2}(V(r)+\frac{L(L+1)\hbar^2}{2m_e r^2}) = E \ u(r) {\tag {3}}$
であたえられる。

**この方程式は合流型超幾何微分方程式へと変換でき，原点で正則な解のエネルギー固有値$E$の厳密解は**

$E_n = \frac{\hbar \omega}{2} (n+\frac{3}{2}),\  ( \  n = 4N+2L+3 {\tag{4}})$

**であることが知られている[2]。**ここでLは軌道量子数， Nは自然数 $N=0,1,2,...$である。

したがってエネルギースペクトルは
**$E = \hbar \omega [1.5, 2.5, 3.5, ...] {\tag{5}})$**

となる。そのほんとんどが縮退している。



本稿では方程式(3)を行列の固有値問題として解きエネルギー固有値$E_n$および固有関数$u_n(r)$を計算する。

***
#コード

コードはすべて[Rydberg原子単位](https://ja.wikipedia.org/wiki/%E5%8E%9F%E5%AD%90%E5%8D%98%E4%BD%8D%E7%B3%BB)を用いる。
これは， 

電子の質量$m=1/2$
ディラック定数$\hbar =1$
長さをボーア$a_{B}=(0.529177 Å)$単位，
エネルギーを$1Ry=13.6058 eV 

とするものである。

```py3

"""
行列法による境界値問題
3次元等方調和振動子ポテンシャル
"""
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt

delta_x = 0.05
x0, x1 = 0.001, 10
N=int((x1-x0)/delta_x)
print("N=",N)

L=0 # 軌道量子数。適宜変更する。

hbar=1
m_elec=1/2
omega=1


y = np.zeros([N-1,N+1])

y[:,0] = 0
y[:,-1] = 0

A=np.zeros([N-1,N-1])
B=np.identity(N-1)

v = np.zeros([N-1])
vcent = np.zeros([N-1])
veff = np.zeros([N-1])





for i in range(N-1): #調和振動子ポンシャル
    x = x0 + i*delta_x
    vcent[i] = L*(L+1)/x**2 # 遠心力ポテンシャル
    v[i] = m_elec*omega**2*(x**2)/2 
    
    veff[i] = v[i] +vcent[i]

for i in range(N-1):  # 3重対角行列なのでindexの位置を気をつけること

    if i == 0:
        A[i,i] = 2/(delta_x**2) + veff[i]
        A[i,i+1] = -1/(delta_x**2)
    elif i == N-2:
        A[i,i-1] = -1/(delta_x**2)
        A[i,i] = 2/(delta_x**2) + veff[i]
    else:
        A[i,i-1] = -1/(delta_x**2)
        A[i,i] = 2/(delta_x**2) + veff[i]
        A[i,i+1] = -1/(delta_x**2)
        
eigen, vec=  scipy.linalg.eigh(A,B)



print("eigen values_3points=",eigen[0:4])

for j in range(N-1):
    for i in range(1,N):
        y[j, i] = vec[i-1,j]

#
# for plot
X= np.linspace(x0,x1, N+1)
plt.plot(X, y[0,:],'-',markersize=5,label='y1')
plt.plot(X, y[1,:],'-',markersize=5,label='y2')
plt.plot(X, y[2,:],'-',markersize=5,label='y3')
plt.legend(loc='upper right')
plt.xlabel('X') # ｘ軸のラベル
plt.ylabel('Y') # y軸のラベル

plt.show()
```
***

#結果
L=0
eigen values_3points= [ 1.46118173  3.44087896  5.42483146]

**厳密値 1.5, 3.5, 5.5とは数%の誤差内で一致している。**


$u_n(r)$ n=0, 2, 4の関数形を下図に示す。

![t.png](https://qiita-image-store.s3.amazonaws.com/0/192457/75c5f19d-2cea-410b-ac88-b19e2b397fbf.png)



L=1
eigen values_3points= [ 2.4997526   4.49899205  6.49760609]

**厳密値 2.5, 4.5, 6.5とは数%の誤差内で一致している。**


$u_n(r)$ n=1, 3, 5の関数形を下図に示す。
![t.png](https://qiita-image-store.s3.amazonaws.com/0/192457/16605c70-40a7-65f7-5cb1-54f172e89d39.png)


***

#補遺
**ポテンシャルの関数形が変わる場合，単にv[i]の計算式を変えるだけでよい。**

***

#参考文献

[1] [[Pythonによる科学・技術計算]行列形式による常微分方程式の境界値問題の解法，数値計算](http://qiita.com/sci_Haru/items/2d238fcadbe7990239df)

[2]後藤 憲一 ほか， [『詳解理論応用量子力学演習』](https://www.amazon.co.jp/%E8%A9%B3%E8%A7%A3%E7%90%86%E8%AB%96%E5%BF%9C%E7%94%A8%E9%87%8F%E5%AD%90%E5%8A%9B%E5%AD%A6%E6%BC%94%E7%BF%92-%E5%BE%8C%E8%97%A4-%E6%86%B2%E4%B8%80/dp/4320031717), 共立出版, 1982.