소개

정상 상태의 슈레딩거 방정식은 차분법에 의해 행렬의 (일반) 고유치 문제에 귀착시킬 수 있다.
본 논문에서는 이 방법을 이용하여 3차원 등방조화 진동자 전위에서 전자의 고유 에너지 및 고유 함수를 결정하는 것을 목적으로 한다.

행렬법의 개요에 대해서는

[Python에 의한 과학·기술 계산] 행렬 형식에 의한 상미분 방정식의 경계값 문제의 해법, 수치 계산 [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 원자 단위을 사용합니다.
이것은

전자 질량 $m=1/2$
디락 상수 $\hbar =1$
길이를 보어 $a_{B}=(0.529177 Å)$ 단위,
에너지를 $1Ry=13.6058 eV

하는 것이다.

"""
行列法による境界値問題
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의 함수형을 아래 그림에 나타낸다.



L=1
eigen values_3points= [ 2.4997526 4.49899205 6.49760609]

엄밀치 2.5, 4.5, 6.5와는 수%의 오차내에서 일치하고 있다.

$u_n(r)$n=1, 3, 5의 함수형을 아래 그림에 나타낸다.

부록

잠재력의 함수 유형이 변경되면 v [i]의 공식을 변경하면됩니다.

참고문헌

[1] [Python에 의한 과학·기술 계산] 행렬 형식에 의한 상미분 방정식의 경계값 문제의 해법, 수치 계산

[2] 고토 겐이치 외, 『상세해리론 응용 양자역학연습』, 공립 출판, 1982.