【ディラック方程式】水素原子の解の可視化
0.概要
以前書いた記事で,シュレディンガー方程式(非相対論)に従う水素原子の電子波動関数を描画した(参考[1]).今回は主に下記2点を実施する.
- ディラック方程式(相対論)に従う水素原子の電子確率密度を可視化
- 非相対論と相対論の解の違いについて議論
水素原子のような原子核電荷の小さな原子の場合,非相対論的エネルギーと相対論的エネルギーの差異は0.001%以下であるが,電子確率密度の空間分布には大きな差異が生まれることが知られている(参考[2,3]).
ディラック方程式に従う水素原子のエネルギー固有関数の中には,upスピンとdownスピンが混じりあった状態のものが存在する(図1).今回行った計算で,参考[2,3]で言及されているように電子確率密度の差異を確認できた.その結果を図2に示す.シュレディンガー方程式に基づく計算において存在していた電子確率密度の"節"が,ディラック方程式に基づく計算では一部,消失していることが分かる.
図1. upスピンとdownスピンの重ね合わせ状態(球面スピノルのLarge成分)
図2. 非相対論と相対論の差異:電子確率密度
図2に示したような差異は電子スピンを考慮した結果生じたものである.よって,シュレディンガー方程式に相対論補正項を追加し,2成分波動関数としてエネルギー固有関数を計算するとシュレディンガー方程式とディラック方程式の計算結果の差異はほぼ無くなると推定される.ただ,今回はそこまでの確認は行わなかった.
また,非相対論(シュレディンガー方程式から規定される状態)と相対論(ディラック方程式から規定される状態)の比較は表1のような規則で行った.表1のような対応は(参考[3])の文献で行われている比較方法を参考にした.この記事では,相対論において全角運動量量子数
表1. 量子状態の対応表:非相対論と相対論
以下で,ディラック方程式に基づく水素原子の電子確率密度の計算方法とその結果について説明する.また,この記事の数式はCGSガウス単位系での表記になっているが,計算結果は全てハートリー原子単位系に変換して出力している.
1.ディラック方程式
この節では下記を確認する.
- ディラック方程式に従う場合の水素原子の電子波動関数の解析解
原子核と電子の相互作用の影響をクーロンポテンシャルで置き換えると,電子波動関数は式(1)のようなディラック方程式に従う.
水素原子の電子波動関数には解析解が存在し,そのエネルギー固有関数
2.計算結果:束縛エネルギー
この節では下記2点を確認する.
- 相対論的束縛エネルギーと非相対論的束縛エネルギーの違い
- 実験で確認されるエネルギー準位のずれを相対論を取り入れると,一部説明可能になる
ディラック方程式から導くことができる電子の相対論的束縛エネルギーは式(19)のようになる.式(19)中の
式(20)からわかるように非相対論において,束縛エネルギーは主量子数
分光記号 |
|
|
||||
---|---|---|---|---|---|---|
1 | 0 | 1/2 | ±1/2 | -5.0000665659e-01 | 0.00000000e+00 | |
2 | 0 | 1/2 | ±1/2 | -1.2500208018e-01 | 0.00000000e+00 | |
2 | 1 | 1/2 | ±1/2 | -1.2500208018e-01 | 0.00000000e+00 | |
2 | 1 | 3/2 | ±1/2, ±3/2 | -1.2500041602e-01 | 1.66416066e-06 | |
3 | 0 | 1/2 | ±1/2 | -5.5556295174e-02 | 0.00000000e+00 | |
3 | 1 | 1/2 | ±1/2 | -5.5556295174e-02 | 0.00000000e+00 | |
3 | 1 | 3/2 | ±1/2, ±3/2 | -5.5555802089e-02 | 4.93084372e-07 | |
3 | 2 | 3/2 | ±1/2, ±3/2 | -5.5555802089e-02 | 4.93084372e-07 | |
3 | 2 | 5/2 | ±1/2, ±3/2, ±5/2 | -5.5555637733e-02 | 6.57440978e-07 | |
4 | 0 | 1/2 | ±1/2 | -3.1250338026e-02 | 0.00000000e+00 | |
4 | 1 | 1/2 | ±1/2 | -3.1250338026e-02 | 0.00000000e+00 | |
4 | 1 | 3/2 | ±1/2, ±3/2 | -3.1250130010e-02 | 2.08015990e-07 | |
4 | 2 | 3/2 | ±1/2, ±3/2 | -3.1250130010e-02 | 2.08015990e-07 | |
4 | 2 | 5/2 | ±1/2, ±3/2, ±5/2 | -3.1250060670e-02 | 2.77355866e-07 | |
4 | 3 | 5/2 | ±1/2, ±3/2, ±5/2 | -3.1250060670e-02 | 2.77355866e-07 | |
4 | 3 | 7/2 | ±1/2, ±3/2, ±5/2, ±7/2 | -3.1250026000e-02 | 3.12025804e-07 |
図3. 相対論により算出されるエネルギー準位のずれ(イメージ図)
3.計算結果:動径方向
この節ではディラック方程式の動径方向の解について,下記3点を実施する.
- ディラック方程式の動径方向の解の"形"をグラフ化して確認
- Large成分の絶対値はSmall成分の絶対値よりも圧倒的に大きな値をとる
- ディラック方程式の動径方向の解のLarge成分は,シュレディンガー方程式の動径方向の解とほぼ一致する.
ディラック方程式では,電子波動関数の動径成分を式(9),(10)のようにあらわすことができる.
動径方向波動関数を計算するプログラム(python)
import numpy as np
from scipy.special import hyp1f1
from scipy.special import gamma as Gamma
from scipy.special import factorial
import matplotlib.pyplot as plt
def calc_g_f(r,n,j,kappa_j,c,alpha,Z):
E_nj = c**2 *(1+ (Z*alpha)**2/(n-j-1/2+((j+1/2)**2 - (Z*alpha)**2)**(1/2))**2 )**(-1/2)
lambda_nj = (c**2 - (E_nj/c)**2)**(1/2)
gamma_j = np.sqrt((j+1/2)**2 -(Z*alpha)**2)
nd_nj = n - (j+1/2)
A_nj = (nd_nj + gamma_j)*c**2/E_nj
#Large成分の動径波動関数
g_coff1 = (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
g_coff2 = (((c**2+E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
g1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
g2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g = g_coff1 * g_coff2 * g1 * (g2 - g3)
#Small成分の動径波動関数
f_coff1 = - (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
f_coff2 = (((c**2-E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
f1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
f2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f = f_coff1 * f_coff2 * f1 * (f2 + f3)
return g, f
#物理定数:ハートリー原子単位系
c = 137.03599160 #光速
alpha = 1/c #微細構造定数
#量子数,物性値(Z:原子核電荷)
Z = 1
n = 4
l = 3
#j = l + 1/2
#kappa_j = -(j + 1/2)
j = l - 1/2
kappa_j = j + 1/2
#計算点
N = 201
rmin, rmax = 0.0001, 50
r = np.linspace(rmin,rmax,N)
#相対論の動径波動関数:Large成分とSmall成分
g_nlj, f_nlj = calc_g_f(r, n, j, kappa_j, c, alpha, Z)
#グラフ描画
plt.rcParams["font.size"] = 16
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(1,1,1)
fig.subplots_adjust(left=0.2)
#ax.set_aspect('equal')
ax.plot(r, g_nlj, label='$g_{nlj}$')
ax.plot(r, f_nlj/(Z*alpha), label='$f_{nlj}/(Z \\alpha)$')
ax.grid()
ax.legend()
ax.set_xlabel('r [bohr]')
ax.set_ylabel('Radial wave function [-]')
plt.show()
#画像の保存
name = "動径波動関数_" + str(n)
name = name + str(l)
name = name + "_" + str(j)
name = name + ".png"
fig.savefig(name)
図4. 動径波動関数:
図5. 動径波動関数:
図6. 動径波動関数:
図7. 動径波動関数:
図8. 動径波動関数:
図9. 動径波動関数:
図10. 動径波動関数:
図11. 動径波動関数:
図12. 動径波動関数:
図13. 動径波動関数:
図14. 動径波動関数:
図15. 動径波動関数:
図16. 動径波動関数:
図18. 動径波動関数:
図19. 動径波動関数:
図20. 動径波動関数:
動径確率密度分布について,相対論(ディラック方程式に従う電子)と非相対論(シュレディンガー方程式に従う電子)の比較をした.動径確率密度分布は相対論では式(23)のように定義され,非相対論では式(24)のように定義される.式(24)中の
動径確率密度分布を非相対論と相対論で比較した結果を図21~26に示す.非相対論と相対論は表2のような対応付けをし,比較した.相対論では
表2. 非相対論と相対論の対応表:動径方向
動径確率密度分布の比較(pythonプログラム)
import numpy as np
from scipy.special import hyp1f1
from scipy.special import gamma as Gamma
from scipy.special import factorial
from scipy.special import assoc_laguerre
import matplotlib.pyplot as plt
def calc_g_f(r,n,j,kappa_j,c,alpha,Z):
E_nj = c**2 *(1+ (Z*alpha)**2/(n-j-1/2+((j+1/2)**2 - (Z*alpha)**2)**(1/2))**2 )**(-1/2)
lambda_nj = (c**2 - (E_nj/c)**2)**(1/2)
gamma_j = np.sqrt((j+1/2)**2 -(Z*alpha)**2)
nd_nj = n - (j+1/2)
A_nj = (nd_nj + gamma_j)*c**2/E_nj
#Large成分の動径波動関数
g_coff1 = (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
g_coff2 = (((c**2+E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
g1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
g2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g = g_coff1 * g_coff2 * g1 * (g2 - g3)
#Small成分の動径波動関数
f_coff1 = - (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
f_coff2 = (((c**2-E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
f1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
f2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f = f_coff1 * f_coff2 * f1 * (f2 + f3)
return g, f
#動径波動関数
def calc_R(r, n, l, Z):
a0 = 1.0
a = a0/Z
zeta = 2.0/n/a * r
#修正係数
n_d = n-l-1
L_coff = (-1)**(2*l+1) * factorial(n+l)
f_coff = - ((2.0/n/a)**3 * factorial(n-l-1)/2.0/n/factorial(n+l)**3)**(1.0/2.0)
f = f_coff * np.exp(-zeta/2.0) * zeta**l * L_coff * assoc_laguerre(zeta, n_d, 2*l+1)
return f
#物理定数:ハートリー原子単位系
c = 137.03599160 #光速
alpha = 1/c #微細構造定数
#量子数,物性値(Z:原子核電荷)
Z = 1
n = 3
l = 2
j = l + 1/2
kappa_j = -(j + 1/2)
E_nj = c**2 *(1+ (Z*alpha)**2/(n-j-1/2+((j+1/2)**2 - (Z*alpha)**2)**(1/2))**2 )**(-1/2)
lambda_nj = (c**2 - (E_nj/c)**2)**(1/2)
gamma_j = np.sqrt((j+1/2)**2 -(Z*alpha)**2)
nd_nj = n - (j+1/2)
A_nj = (nd_nj + gamma_j)*c**2/E_nj
#計算点
N = 401
rmin, rmax = 0.0001, 40
r = np.linspace(rmin,rmax,N)
# --- 相対論 ---
#動径波動関数:Large成分とSmall成分
g_nlj, f_nlj = calc_g_f(r, n, j, kappa_j, c, alpha, Z)
#動径確率密度分布
P_nlj_rel = r**2 * (g_nlj**2 + f_nlj**2)
# --- 非相対論 ---
#動径波動関数
R_nl = calc_R(r, n, l, Z)
#動径確率密度分布
P_nl = (r * R_nl)**2
#グラフ描画
plt.rcParams["font.size"] = 16
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(1,1,1)
fig.subplots_adjust(left=0.2)
ax.plot(r, P_nlj_rel, label='relativistic')
ax.plot(r, P_nl, label='non-relativistic')
ax.grid()
ax.legend()
ax.set_xlabel('r [bohr]')
ax.set_ylabel('Radial probability density distribution [-]')
fig.show()
name = "相対論と非相対論の比較_"+ "R2_" + str(n) + str(l)+ "_" + str(j)+ "_" + ".png"
fig.savefig(name)
図21. 動径確率密度分布の比較,相対論:
図22. 動径確率密度分布の比較,相対論:
図23. 動径確率密度分布の比較,相対論:
図24. 動径確率密度分布の比較,相対論:
図25. 動径確率密度分布の比較,相対論:
図26. 動径確率密度分布の比較,相対論:
4.計算結果:球面スピノル
この節ではディラック方程式の角度方向の解について,下記4点を実施する.
- ディラック方程式の角度方向の解(球面スピノル)を可視化
- ディラック方程式とシュレディンガー方程式の角度方向の解を比較
- ディラック方程式とシュレディンガー方程式の電子確率密度の違いは,角度方向の解の違いに起因していることを確認
- ディラック方程式とシュレディンガー方程式の角度方向の解が一致する場合と一致しない場合について考察
4成分波動関数
球面スピノルの描画に使用したプログラム(python)
import sys
import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
# 球面調和関数のノーテーションや文字を物理の教科書用にカスタマイズ
# l < m のときnanを返さないように場合分け
def my_sph_harm(l, m, theta, phi):
if l < abs(m) :
f = 0.0
else:
f = sph_harm(m, l, phi, theta)
return f
# (r,Θ,φ)空間から(x,y,z)空間への射影
def polar(r, theta, phi):
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
return x, y, z
flag = 1 # j = l + 1/2のときは1, j = l - 1/2のときはそれ以外
l = 2 # 方位量子数, 0以上の整数値
if flag == 1:
j = l + 1/2
else:
j = l - 1/2
m = 1 + 1/2 # 磁気量子数, -j<= m <=j を満たす半整数
N = 101 # theta, phiの分割数, 0以上の整数値
if abs(m) > j:
print("mの絶対値はj以下に設定する必要があります.")
print("mの値を見直して下さい.")
print("プログラムを強制終了します.")
sys.exit()
# theta, phi: 定義域設定
theta_min = 0.0
theta_max = np.pi
phi_min = 0.0
phi_max = 2.0 * np.pi
# theta, phi: メッシュ作成
theta_a = np.linspace(theta_min, theta_max, N)
phi_a = np.linspace(phi_min, phi_max, N)
theta, phi = np.meshgrid(theta_a, phi_a)
# j = l+1/2 or l-1/2 かの場合分け
if flag == 1:
# j = l+1/2 のとき
omega_1 = np.sqrt((j+m)/(2*j)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j-m)/(2*j)) * my_sph_harm(l, m+1/2, theta, phi)
omega_3 = - np.sqrt((j-m+1)/(2*j+2)) * my_sph_harm(l+1, m-1/2, theta, phi)
omega_4 = np.sqrt((j+m+1)/(2*j+2)) * my_sph_harm(l+1, m+1/2, theta, phi)
else:
# j = l-1/2 のとき
omega_1 = - np.sqrt((j-m+1)/(2*j+2)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j+m+1)/(2*j+2)) * my_sph_harm(l, m+1/2, theta, phi)
omega_3 = np.sqrt((j+m)/(2*j)) * my_sph_harm(l-1, m-1/2, theta, phi)
omega_4 = np.sqrt((j-m)/(2*j)) * my_sph_harm(l-1, m+1/2, theta, phi)
omega_L_u = np.sqrt(np.abs(omega_1)**2)
omega_L_d = np.sqrt(np.abs(omega_2)**2)
omega_S_u = np.sqrt(np.abs(omega_3)**2)
omega_S_d = np.sqrt(np.abs(omega_4)**2)
omega_L = np.sqrt(np.abs(omega_1)**2 + np.abs(omega_2)**2)
omega_S = np.sqrt(np.abs(omega_3)**2 + np.abs(omega_4)**2)
#グラフ共通設定
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["font.size"] = 20
fig = plt.figure(figsize=(21,7))
ax_lim = 0.5
alpha = 0.5
rst = 3
cst = 3
#Large, up
r = omega_L_u
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 0), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Large}_{\\rm up}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#Large, down
r = omega_L_d
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 1), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Large}_{\\rm down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#Large, up+down
r = omega_L
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 2), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Large}_{\\rm up+down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#画像の保存:Large成分
name = "球面スピノル_" + str(l)
name = name + "_" + str(j)
name = name + "_" + str(m)
name = name + "_Large.png"
fig.savefig(name)
#Small, up
r = omega_S_u
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 0), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Small}_{\\rm up}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#Small, down
r = omega_S_d
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 1), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Small}_{\\rm down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#Small, up+down
r = omega_S
x, y, z = polar(r, theta, phi)
ax = plt.subplot2grid((1, 3), (0, 2), projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Smal}_{\\rm up+down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#画像の保存:Small成分
name = "球面スピノル_" + str(l)
name = name + "_" + str(j)
name = name + "_" + str(m)
name = name + "_Small.png"
fig.savefig(name)
図27-1. Large成分の球面スピノル
図27-2. Small成分の球面スピノル
図28-1. Large成分の球面スピノル
図28-2. Small成分の球面スピノル
図29-1. Large成分の球面スピノル
図29-2. Small成分の球面スピノル
図30-1. Large成分の球面スピノル
図30-2. Small成分の球面スピノル
図31-1. Large成分の球面スピノル
図31-2. Small成分の球面スピノル
図32-1. Large成分の球面スピノル
図32-2. Small成分の球面スピノル
図33-1. Large成分の球面スピノル
図33-2. Small成分の球面スピノル
図34-1. Large成分の球面スピノル
図34-2. Small成分の球面スピノル
図35-1. Large成分の球面スピノル
図35-2. Small成分の球面スピノル
ここで,非相対論(シュレディンガー方程式に従う電子)と相対論(ディラック方程式に従う電子)における電子確率密度がどのように記述されるかを整理し,非相対論と相対論の角度方向の解を比較する際に具体的に何を比較すれば意味のある議論ができるかを考える.
シュレディンガー方程式に従う電子波動関数を
相対論(ディラック方程式)から導かれれる電子確率密度は
非相対論と相対論において電子確率密度にどのような違いが出るかを議論する際,式(35),(36)から電子確率密度の違いは主にその角度成分
非相対論と相対論は表3のような対応付けをし,比較した.相対論では
以前書いた記事では非相対論(シュレディンガー方程式)の角度方向の解を図示する際に,球面調和関数
表3. 非相対論と相対論の対応表:角度方向
角度方向の解の比較(python)
import sys
import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
# 球面調和関数のノーテーションや文字を物理の教科書用にカスタマイズ
# l < m のときnanを返さないように場合分け
def my_sph_harm(l, m, theta, phi):
if l < abs(m) :
f = 0.0
else:
f = sph_harm(m, l, phi, theta)
return f
# (r,Θ,φ)空間から(x,y,z)空間への射影
def polar(r, theta, phi):
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
return x, y, z
# 相対論の量子数指定
flag = 1 # j = l + 1/2のときは1, j = l - 1/2のときはそれ以外
l = 2 # 方位量子数, 0以上の整数値
if flag == 1:
j = l + 1/2
else:
j = l - 1/2
m = 2 + 1/2 # 磁気量子数, -j<= m <=j を満たす半整数
# 非相対論の量子数指定
ms = 2
N = 101 # theta, phiの分割数, 0以上の整数値
if abs(m) > j:
print("mの絶対値はj以下に設定する必要があります.")
print("mの値を見直して下さい.")
print("プログラムを強制終了します.")
sys.exit()
# theta, phi: 定義域設定
theta_min = 0.0
theta_max = np.pi
phi_min = 0.0
phi_max = 2.0 * np.pi
# theta, phi: メッシュ作成
theta_a = np.linspace(theta_min, theta_max, N)
phi_a = np.linspace(phi_min, phi_max, N)
theta, phi = np.meshgrid(theta_a, phi_a)
# j = l+1/2 or l-1/2 かの場合分け
if flag == 1:
# j = l+1/2 のとき
omega_1 = np.sqrt((j+m)/(2*j)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j-m)/(2*j)) * my_sph_harm(l, m+1/2, theta, phi)
else:
# j = l-1/2 のとき
omega_1 = - np.sqrt((j-m+1)/(2*j+2)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j+m+1)/(2*j+2)) * my_sph_harm(l, m+1/2, theta, phi)
omega_L = np.sqrt(np.abs(omega_1)**2 + np.abs(omega_2)**2)
Y_lm = my_sph_harm(l, ms, theta, phi)
Y_lm_abs = np.abs(Y_lm)
#グラフ共通設定
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["font.size"] = 20
fig = plt.figure(figsize=(14,14))
ax_lim = 0.5
alpha = 0.5
rst = 3
cst = 3
#non-relativistic
r = Y_lm_abs
x, y, z = polar(r, theta, phi)
ax = fig.add_subplot(2,2,1,projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = |Y_l^m (\\theta, \phi)|$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
ax = fig.add_subplot(2,2,3,projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = |Y_l^m (\\theta, \phi)|$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.view_init(elev=0)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#Large, up+down:relativistic
r = omega_L
x, y, z = polar(r, theta, phi)
ax = fig.add_subplot(2,2,2,projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Large}_{\\rm up+down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
ax = fig.add_subplot(2,2,4,projection='3d')
ax.set_box_aspect((1, 1, 1))
ax.set_title('$r(\\theta, \phi) = r^{\\rm Large}_{\\rm up+down}(\\theta, \phi)$')
ax.set_xlim(-ax_lim, ax_lim)
ax.set_ylim(-ax_lim, ax_lim)
ax.set_zlim(-ax_lim, ax_lim)
ax.view_init(elev=0)
ax.plot_surface(x, y, z, color="red" , alpha = alpha, linewidth=0, rstride=rst, cstride=cst)
#画像の保存
name = "相対論と非相対論の比較_角度方向_" + str(l)
name = name + "_" + str(j)
name = name + "_" + str(m)
name = name + ".png"
fig.savefig(name)
図36~41から,相対論における量子数
逆に,相対論における量子数
図36. 角度方向の解の比較,非相対論:
図37. 角度方向の解の比較,非相対論:
図38. 角度方向の解の比較,非相対論:
図39. 角度方向の解の比較,非相対論:
図40. 角度方向の解の比較,非相対論:
図41. 角度方向の解の比較,非相対論:
5.計算結果:電子確率密度
この節では下記2点を実施する.
- ディラック方程式に従う場合の電子確率密度を可視化
- ディラック方程式とシュレディンガー方程式の電子確率密度の差異を確認
シュレディンガー方程式に従う電子の電子確率密度は式(35)のようにあらわせる.式(35)は下に再掲.ディラック方程式に従う電子の電子確率密度は式(39)のようにあらわせる.式(39)は式(36)に
3節と4節で行った議論から下記2点が予想される.表4中で「〇」をつけているのは相対論と非相対論の電子確率密度がほぼ一致すると予想されもので,「×」をつけているのは電子確率密度が角度方向の解の分だけ違いが出ると予想されるものである.
- ディラック方程式(相対論)における,量子数
がm となるとき,相対論とそれに対応する非相対論(シュレディンガー方程式)の電子確率密度はほぼ一致する.m=j - ディラック方程式(相対論)における,量子数
がm となるとき,相対論とそれに対応する非相対論(シュレディンガー方程式)の電子確率密度は角度方向の解の分だけ違いがでる.m\ne j
表4. 非相対論と相対論の比較表:電子確率密度
電子確率密度を図42~55に示す.上述したような規則で,ディラック方程式(相対論)とシュレディンガー方程式(非相対論)の電子確率密度の差異を確認できる.可視化に使用したプログラムは下記に折り畳んで置いた.
電子確率密度を可視化するためのプログラム(python):非相対論
import plotly.graph_objects as go
import numpy as np
from scipy.special import sph_harm
from scipy.special import factorial
from scipy.special import assoc_laguerre
import plotly.io as pio
#pio.renderers.default='svg'
pio.renderers.default='browser'
#動径波動関数
def R_nl(r, n, l):
a0 = 1.0
zeta = 2.0/n/a0 * r
#修正係数
n_d = n-l-1
L_coff = (-1)**(2*l+1) * factorial(n+l)
f_coff = - ((2.0/n/a0)**3 * factorial(n-l-1)/2.0/n/factorial(n+l)**3)**(1.0/2.0)
f = f_coff * np.exp(-zeta/2.0) * zeta**l * L_coff * assoc_laguerre(zeta, n_d, 2*l+1)
return f
# 球面調和関数のノーテーションや文字を物理の教科書用にカスタマイズ
def my_sph_harm(l, m, theta, phi):
return sph_harm(m, l, phi, theta)
# 量子数の指定
n = 2 # 主量子数: 1 <= n
l = 1 # 方位量子数 0 <= l <= n-1
m = 0 # 磁気量子数 -l <= m <= l
# 計算領域設定
N = 50 # 分割数
area = 10
xmin, ymin, zmin = -area, -area, -area
xmax, ymax, zmax = area, area, area
# メッシュ作成
x_a = np.linspace(xmin,xmax,N)
y_a = np.linspace(ymin,ymax,N)
z_a = np.linspace(zmin,zmax,N)
x, y, z = np.meshgrid(x_a, y_a, z_a)
# 極座標化
r = np.sqrt(x**2+y**2+z**2)
theta = np.arccos(z/r)
phi = np.sign(y) * np.arccos(x/np.sqrt(x**2 + y**2))
# 電子波動関数の計算
R_wave = R_nl(r, n, l) # 動径波動関数
Y_wave = my_sph_harm(l, m, theta, phi) # 複素数の球面調和関数
psi = R_wave * Y_wave # 電子波動関数
psi2 = np.abs(psi)**2
# 配列の1次元化
X = x.flatten()
Y = y.flatten()
Z = z.flatten()
Psi = psi2.flatten()
# カラーバー範囲設定
Max = np.fabs(np.amax(psi2))
Max = 1/2 * Max
# グラフの描画
fig = go.Figure(data=go.Volume(
x=X,
y=Y,
z=Z,
value=Psi,
isomax=Max,
opacity=0.1, # 透明度:0~1
surface_count=20, # 描画する等値面の数
opacityscale="uniform", # 透明度のカスタマイズ
))
# グラフレイアウトの調整
fig.update_layout(
scene_camera= dict(eye=dict(x=1.5, y=1.5, z=0.5)) # カメラアングル
)
# 図の表示
fig.show()
# グラフの保存
name = "psi2_" + str(n)
name = name + str(l)
name = name + str(m)
name = name + ".png"
fig.write_image(name)
電子確率密度を可視化するためのプログラム(python):相対論
import sys
import numpy as np
from scipy.special import hyp1f1
from scipy.special import gamma as Gamma
from scipy.special import factorial
from scipy.special import sph_harm
import plotly.graph_objects as go
import plotly.io as pio
#pio.renderers.default='svg'
pio.renderers.default='browser'
def calc_g_f(r,n,j,kappa_j,c,alpha,Z):
E_nj = c**2 *(1+ (Z*alpha)**2/(n-j-1/2+((j+1/2)**2 - (Z*alpha)**2)**(1/2))**2 )**(-1/2)
lambda_nj = (c**2 - (E_nj/c)**2)**(1/2)
gamma_j = np.sqrt((j+1/2)**2 -(Z*alpha)**2)
nd_nj = n - (j+1/2)
A_nj = (nd_nj + gamma_j)*c**2/E_nj
#Large成分の動径波動関数
g_coff1 = (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
g_coff2 = (((c**2+E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
g1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
g2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
g = g_coff1 * g_coff2 * g1 * (g2 - g3)
#Small成分の動径波動関数
f_coff1 = - (2*lambda_nj)**(3/2)/Gamma(2*gamma_j+1)
f_coff2 = (((c**2-E_nj)*Gamma(2*gamma_j+nd_nj+1))/(4*c**2*A_nj*(A_nj-kappa_j)*factorial(nd_nj)))**(1/2)
f1 = (2*lambda_nj*r)**(gamma_j-1) * np.exp(-lambda_nj*r)
f2 = (A_nj-kappa_j)*hyp1f1(-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f3 = nd_nj*hyp1f1(1-nd_nj, 2*gamma_j+1, 2*lambda_nj*r)
f = f_coff1 * f_coff2 * f1 * (f2 + f3)
return g, f
def my_sph_harm(l, m, theta, phi):
if l < abs(m) :
f = 0.0
else:
f = sph_harm(m, l, phi, theta)
return f
#物理定数:ハートリー原子単位系
c = 137.03599160 #光速
alpha = 1/c #微細構造定数
#量子数,物性値(Z:原子核電荷)
Z = 1
n = 2
l = 1 # 方位量子数, 0以上の整数値
flag = 1 # j = l + 1/2のときは1, j = l - 1/2のときはそれ以外
if flag == 1:
j = l + 1/2
kappa_j = -(j + 1/2)
else:
j = l - 1/2
kappa_j = j + 1/2
m = 1 + 1/2 # 磁気量子数, -j<= m <=j を満たす半整数
#m = -(1 + 1/2) # 磁気量子数, -j<= m <=j を満たす半整数
if abs(l) > n :
print("l ≦ n-1 のように l は設定する必要があります.")
print("プログラムを強制終了します.")
sys.exit()
if abs(m) > j:
print("mの絶対値はj以下に設定する必要があります.")
print("mの値を見直して下さい.")
print("プログラムを強制終了します.")
sys.exit()
# 計算領域設定
N = 50 # 分割数
area = 10
xmin, ymin, zmin = -area, -area, -area
xmax, ymax, zmax = area, area, area
# メッシュ作成
x_a = np.linspace(xmin,xmax,N)
y_a = np.linspace(ymin,ymax,N)
z_a = np.linspace(zmin,zmax,N)
x, y, z = np.meshgrid(x_a, y_a, z_a)
# 極座標化
r = np.sqrt(x**2+y**2+z**2)
theta = np.arccos(z/r)
phi = np.sign(y) * np.arccos(x/np.sqrt(x**2 + y**2))
# 相対論の動径波動関数:Large成分とSmall成分
g_nj, f_nj = calc_g_f(r, n, j, kappa_j, c, alpha, Z)
# 球面スピノルの計算
# j = l+1/2 or l-1/2 かの場合分け
if flag == 1:
# j = l+1/2 のとき
omega_1 = np.sqrt((j+m)/(2*j)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j-m)/(2*j)) * my_sph_harm(l, m+1/2, theta, phi)
omega_3 = - np.sqrt((j-m+1)/(2*j+2)) * my_sph_harm(l+1, m-1/2, theta, phi)
omega_4 = np.sqrt((j+m+1)/(2*j+2)) * my_sph_harm(l+1, m+1/2, theta, phi)
else:
# j = l-1/2 のとき
omega_1 = - np.sqrt((j-m+1)/(2*j+2)) * my_sph_harm(l, m-1/2, theta, phi)
omega_2 = np.sqrt((j+m+1)/(2*j+2)) * my_sph_harm(l, m+1/2, theta, phi)
omega_3 = np.sqrt((j+m)/(2*j)) * my_sph_harm(l-1, m-1/2, theta, phi)
omega_4 = np.sqrt((j-m)/(2*j)) * my_sph_harm(l-1, m+1/2, theta, phi)
# 4成分波動関数の計算
psi_1 = 1.j * g_nj * omega_1
psi_2 = 1.j * g_nj * omega_2
psi_3 = - f_nj * omega_3
psi_4 = - f_nj * omega_4
# 電子確率密度の計算
Psi2 = np.abs(psi_1)**2 + np.abs(psi_2)**2 + np.abs(psi_3)**2 + np.abs(psi_4)**2
# 配列の1次元化
XX = x.flatten()
YY = y.flatten()
ZZ = z.flatten()
PPsi2 = Psi2.flatten()
# カラーバー範囲設定
Max = np.fabs(np.amax(Psi2))
Max = 1/2 * Max
# グラフの描画
fig = go.Figure(data=go.Volume(
x=XX,
y=YY,
z=ZZ,
value=PPsi2,
isomax=Max,
opacity=0.1, # 透明度:0~1
surface_count=20, # 描画する等値面の数
opacityscale="uniform", # 透明度のカスタマイズ
))
# グラフレイアウトの調整
fig.update_layout(
scene_camera= dict(eye=dict(x=1.5, y=1.5, z=0.5)) # カメラアングル
)
# 図の表示
fig.show()
# グラフの保存
name = "相対論_"+ "Psi2_" + str(n) + str(l)+ "_" + str(j)+ "_" + str(m) + ".png"
fig.write_image(name)
図42. 電子確率密度,非相対論:
図43. 電子確率密度,非相対論:
図44. 電子確率密度,相対論:
図45. 電子確率密度,非相対論:
図46. 電子確率密度,非相対論:
図47. 電子確率密度,非相対論:
図48. 電子確率密度,相対論:
図49. 電子確率密度,非相対論:
図50. 電子確率密度,非相対論:
図51. 電子確率密度,相対論:
図52. 電子確率密度,相対論:
図53. 電子確率密度,非相対論:
図54. 電子確率密度,非相対論:
図55. 電子確率密度,非相対論:
6.結論
- ディラック方程式のエネルギー固有関数の中にはupスピンとdownスピンの重ね合わせ状態をとるものが存在する.
- その結果,ディラック方程式とシュレディンガー方程式から導かれる電子確率密度に,明確な差異が生じる場合がある.
- その差異は主に角度方向の解の違いに起因している.
環境
・python 3.9.13
・IPython 8.4.0
・scipy 1.7.3
・numpy 1.21.5
・plotly 5.9.0
・matplotlib 3.5.2
・Anaconda 4.14.0
・Spyder 5.9.2
・Windows 10 Home 64bit
参考
[1]. 【Plotly】水素原子の図鑑:電子波動関数の可視化(参照2022-10-23)
[2]. H.E.White. "Pictorial Representations of the Dirac Electron Cloud for Hydrogen-Like Atoms". Phys. Rev., 38, 513(1931). (https://doi.org/10.1103/PhysRev.38.513)
[3]. 秦野, 山本, 舘脇. "Dirac方程式の厳密解の可視化",Journal of Computer Chemistry. Japan, 2016, 15, 105-117. (https://doi.org/10.2477/jccj.2016-0014)
[4]. W.Greiner. "Relativistic Quantum Mechanics. Wave Equation 3rd Edition". Springer, 2000, p.210-233.
[5]. 猪木, 川合. 量子力学Ⅰ. 講談社サイエンティフィク, 2010, p.151.
[6]. 【Python】水素原子の可視化:ラゲール陪多項式計算モジュール(参照2022-10-23)
[7]. 【Python】水素原子の可視化:球面調和関数(参照2022-10-23)
Discussion