【ディラック方程式】水素原子の解の可視化

この記事では，電子スピンを考慮しない1成分シュレディンガー方程式に基づく議論のことを「非相対論」，ディラック方程式に基づく議論のことを「相対論」と便宜上呼んでいます．

0.概要

以前書いた記事で，シュレディンガー方程式（非相対論）に従う水素原子の電子波動関数を描画した(参考[1])．今回は主に下記2点を実施する．

• ディラック方程式（相対論）に従う水素原子の電子確率密度を可視化
• 非相対論と相対論の解の違いについて議論

水素原子のような原子核電荷の小さな原子の場合，非相対論的エネルギーと相対論的エネルギーの差異は0.001%以下であるが，電子確率密度の空間分布には大きな差異が生まれることが知られている(参考[2,3])．

ディラック方程式に従う水素原子のエネルギー固有関数の中には，upスピンとdownスピンが混じりあった状態のものが存在する（図1）．今回行った計算で，参考[2,3]で言及されているように電子確率密度の差異を確認できた．その結果を図2に示す．シュレディンガー方程式に基づく計算において存在していた電子確率密度の"節"が，ディラック方程式に基づく計算では一部，消失していることが分かる．

図1. upスピンとdownスピンの重ね合わせ状態(球面スピノルのLarge成分)

図2. 非相対論と相対論の差異：電子確率密度

図2に示したような差異は電子スピンを考慮した結果生じたものである．よって，シュレディンガー方程式に相対論補正項を追加し，2成分波動関数としてエネルギー固有関数を計算するとシュレディンガー方程式とディラック方程式の計算結果の差異はほぼ無くなると推定される．ただ，今回はそこまでの確認は行わなかった．

また，非相対論（シュレディンガー方程式から規定される状態）と相対論（ディラック方程式から規定される状態）の比較は表1のような規則で行った．表1のような対応は(参考[3])の文献で行われている比較方法を参考にした．この記事では，相対論において全角運動量量子数$j$が$j=l-1/2$となるような量子状態をとるとき，それに対応する非相対論の量子状態は無いものとしている．これは，電子スピンを考慮した量子状態と，電子スピンを考慮していない量子状態を1対1で対応付けることが困難なためである．

表1. 量子状態の対応表：非相対論と相対論

以下で，ディラック方程式に基づく水素原子の電子確率密度の計算方法とその結果について説明する．また，この記事の数式はCGSガウス単位系での表記になっているが，計算結果は全てハートリー原子単位系に変換して出力している．

1.ディラック方程式

この節では下記を確認する．

• ディラック方程式に従う場合の水素原子の電子波動関数の解析解

原子核と電子の相互作用の影響をクーロンポテンシャルで置き換えると，電子波動関数は式(1)のようなディラック方程式に従う．$c$，$m_0$，$e$，$Z$，$E_{nj}$，$\bm{\Psi}_{nljm}(x,y,z)$はそれぞれ，光速，電子の静止質量，電気素量，原子核電荷，エネルギー固有値，電子波動関数のエネルギー固有関数である（水素原子の場合は$Z=1$）．非相対論との比較が行いやすいディラック表示では，$\bm{\hat{\alpha}}$，$\hat{\beta}$を式(3)，(4)のような表現行列で書くことができる．この記事では電子波動関数をディラック表示で記述する．

$$\begin{gather} \left(c \bm{\hat{\alpha}} \cdot \bm{\hat{p}} + \hat{\beta} m_0 c^2 - \frac{Ze^2}{r} \right) \bm{\Psi}_{nljm}(x,y,z) = E_{nj} \bm{\Psi}_{nljm}(x,y,z)\\ \hat{p}_1 = - i \hbar \frac{\partial}{\partial x}, \ \ \hat{p}_2 = - i \hbar \frac{\partial}{\partial y}, \ \ \hat{p}_3 = - i \hbar \frac{\partial}{\partial z}\\ \hat{\alpha}_1 = \begin{pmatrix} 0 &0 &0 &1\\ 0 &0 &1 &0\\ 0 &1 &0 &0\\ 1 &0 &0 &0 \end{pmatrix} , \ \ \hat{\alpha}_2 = \begin{pmatrix} 0 &0 &0 &-i\\ 0 &0 &i &0\\ 0 &-i &0 &0\\ i &0 &0 &0 \end{pmatrix} , \ \ \hat{\alpha}_3 = \begin{pmatrix} 0 &0 &1 &0\\ 0 &0 &0 &-1\\ 1 &0 &0 &0\\ 0 &-1 &0 &0 \end{pmatrix}\\ \hat{\beta} = \begin{pmatrix} 1 &0 &0 &0\\ 0 &1 &0 &0\\ 0 &0 &-1 &0\\ 0 &0 &0 &-1 \end{pmatrix}\\ \end{gather}$$

水素原子の電子波動関数には解析解が存在し，そのエネルギー固有関数$\bm{\Psi}_{nljm} (r, \theta, \phi)$は極座標表示で，式(5)~(18)のように書くことができる(例えば(参考[4]))．ただし，式(7)，(8)中の$Y_l^{m} (\theta, \phi)$は球面調和関数で，式(9)，(10)中の$\Gamma(x)$，$F(a,\ b;\ x)$はそれぞれガンマ関数，合流型超幾何関数である．$n$，$l$は式(18)を満たす整数で，$j$，$m$は式(18)を満たす半整数である．

$$\begin{gather} \bm{\Psi}_{nljm} (r, \theta, \phi) = \begin{pmatrix} i g_{nlj}(r) \bm{\Omega}_{jlm} (\theta, \phi) \\ - f_{nlj}(r) \bm{\Omega}_{jl'm} (\theta, \phi)\\ \end{pmatrix}\\ l' = \begin{cases} l+1 &(j = l + \frac{1}{2}) \\ l-1 &(j = l - \frac{1}{2})\\ \end{cases}\\ \begin{cases} \bm{\Omega}_{jlm} (\theta, \phi) = \begin{pmatrix} \sqrt{\frac{j+m}{2j}} Y_l^{m-\frac{1}{2}} (\theta, \phi) \\ \sqrt{\frac{j-m}{2j}} Y_l^{m+\frac{1}{2}} (\theta, \phi) \\ \end{pmatrix} &(j = l + \frac{1}{2}) \\ \bm{\Omega}_{jl'm} (\theta, \phi) = \begin{pmatrix} -\sqrt{\frac{j-m+1}{2j+2}} Y_{l+1}^{m-\frac{1}{2}} (\theta, \phi) \\ \sqrt{\frac{j+m+1}{2j+2}} Y_{l+1}^{m+\frac{1}{2}} (\theta, \phi) \\ \end{pmatrix} &(j = l + \frac{1}{2}) \\ \end{cases}\\ \begin{cases} \bm{\Omega}_{jlm} (\theta, \phi) = \begin{pmatrix} -\sqrt{\frac{j-m+1}{2j+2}} Y_l^{m-\frac{1}{2}} (\theta, \phi) \\ \sqrt{\frac{j+m+1}{2j+2}} Y_l^{m+\frac{1}{2}} (\theta, \phi) \\ \end{pmatrix} &(j = l - \frac{1}{2}) \\ \bm{\Omega}_{jl'm} (\theta, \phi) = \begin{pmatrix} \sqrt{\frac{j+m}{2j}} Y_{l-1}^{m-\frac{1}{2}} (\theta, \phi) \\ \sqrt{\frac{j-m}{2j}} Y_{l-1}^{m+\frac{1}{2}} (\theta, \phi) \\ \end{pmatrix} &(j = l - \frac{1}{2}) \\ \end{cases}\\ \begin{split} g_{nlj} (r) = \frac{\left(2 \lambda_{nj}\right)^{\frac{3}{2}}}{\Gamma (2 \gamma_j +1)} \times \left( \frac{(m_0 c^2 + E_{nj}) \Gamma(2\gamma_{j}+n'_{nj})}{4 m_0 c^2 A_{nj} (A_{nj} -\kappa_j) n'_{nj}!} \right)^{\frac{1}{2}} \times \left(2 \lambda_{nj} r \right)^{\gamma_j -1} \times e^{- \lambda_{nj} r} \\ \times \left\{ (A_{nj} - \kappa_j) F(-n'_{nj},\ 2\gamma_j+1;\ 2\lambda_{nj}r) -n'_{nj} F(1-n'_{nj},\ 2\gamma_j+1;\ 2\lambda_{nj}r) \right\} \end{split}\\ \begin{split} f_{nlj} (r) = \frac{-\left(2 \lambda_{nj}\right)^{\frac{3}{2}}}{\Gamma (2 \gamma_j +1)} \times \left( \frac{(m_0 c^2 - E_{nj}) \Gamma(2\gamma_{j}+n'_{nj})}{4m_0 c^2 A_{nj} (A_{nj} -\kappa_j) n'_{nj}!} \right)^{\frac{1}{2}} \times \left(2 \lambda_{nj} r \right)^{\gamma_j -1} \times e^{- \lambda_{nj} r} \\ \times \left\{ (A_{nj} - \kappa_j) F(-n'_{nj},\ 2\gamma_j+1;\ 2\lambda_{nj}r) +n'_{nj} F(1-n'_{nj},\ 2\gamma_j+1;\ 2\lambda_{nj}r) \right\} \end{split}\\ \alpha = \frac{e^2}{\hbar c}\\ E_{nj} = m_0 c^2 \left( 1+ \frac{(Z\alpha)^2}{\left( n-j-\frac{1}{2} + \left( \left(j + \frac{1}{2}\right)^2 - \left(Z \alpha \right)^2 \right)^{\frac{1}{2}} \right)^2} \right)^{-\frac{1}{2}}\\ \lambda_{nj} = \frac{\left(m_0^2 c^4 - E_{nj}^2 \right)^{\frac{1}{2}}}{\hbar c}\\ \gamma_j = \sqrt{\left(j+\frac{1}{2}\right)^2 - \left(Z \alpha \right)^2}\\ n'_{nj} = n -\left(j+\frac{1}{2}\right)\\ A_{nj} = \frac{(n'+\gamma_j)m_0 c^2}{E_{nj}}\\ \kappa_j = \begin{cases} -(j + \frac{1}{2}) &(j = l + \frac{1}{2}) \\ j + \frac{1}{2} &(j = l - \frac{1}{2})\\ \end{cases}\\ \begin{cases} n \geq 1 \\ l \geq n-1 \\ j = l + \frac{1}{2} \ \ {\rm{or}} \ \ l -\frac{1}{2} \ \ (j \geq \frac{1}{2}) \\ -j \leq m \leq j \end{cases}\\ \end{gather}$$

2.計算結果：束縛エネルギー

この節では下記2点を確認する．

• 相対論的束縛エネルギーと非相対論的束縛エネルギーの違い
• 実験で確認されるエネルギー準位のずれを相対論を取り入れると，一部説明可能になる

ディラック方程式から導くことができる電子の相対論的束縛エネルギーは式(19)のようになる．式(19)中の$E_{nj}$は式(12)から算出することができる．シュレディンガー方程式から導くことができる電子の非相対論的束縛エネルギーは式(20)のようになる(例えば(参考[5]))．$a_0$はボーア半径である．

$$\begin{gather} E^{\rm (bind:rela)}_{nj} = E_{nj} - m_0 c^2 \\ E^{\rm (bind:non-rela)}_{n} = - \frac{e^2}{2 a_0 n^2} \\ \end{gather}$$

式(20)からわかるように非相対論において，束縛エネルギーは主量子数$n$のみに依存する．それに対し，相対論的束縛エネルギーは主量子数$n$と全角運動量量子数$j$に依存する．相対論では，全角運動量量子数$j$が異なる状態でのエネルギー準位の分裂を説明することができる．下表に相対論で計算される束縛エネルギー$E^{\rm (bind:rela)}_{nj}$と，式(21)で定義されるエネルギー準位のずれ$\Delta E$を示す．また，図3にエネルギー準位のずれのイメージ図を示す．

$$\begin{gather} \Delta E = E^{\rm (bind:rela)}_{nj} - E^{\rm (bind:rela)}_{nj=1/2}\\ \end{gather}$$

分光記号	$\bm n$	$\bm l$	$\bm j$	$\bm m$	${\bm E^{\rm rela}_{\rm bind}}$ [hartree]	${\bm \Delta \bm E}$ [hartree]
$1S_{1/2}$	1	0	1/2	±1/2	-5.0000665659e-01	0.00000000e+00
$2S_{1/2}$	2	0	1/2	±1/2	-1.2500208018e-01	0.00000000e+00
$2P_{1/2}$	2	1	1/2	±1/2	-1.2500208018e-01	0.00000000e+00
$2P_{3/2}$	2	1	3/2	±1/2, ±3/2	-1.2500041602e-01	1.66416066e-06
$3S_{1/2}$	3	0	1/2	±1/2	-5.5556295174e-02	0.00000000e+00
$3P_{1/2}$	3	1	1/2	±1/2	-5.5556295174e-02	0.00000000e+00
$3P_{3/2}$	3	1	3/2	±1/2, ±3/2	-5.5555802089e-02	4.93084372e-07
$3D_{3/2}$	3	2	3/2	±1/2, ±3/2	-5.5555802089e-02	4.93084372e-07
$3D_{5/2}$	3	2	5/2	±1/2, ±3/2, ±5/2	-5.5555637733e-02	6.57440978e-07
$4S_{1/2}$	4	0	1/2	±1/2	-3.1250338026e-02	0.00000000e+00
$4P_{1/2}$	4	1	1/2	±1/2	-3.1250338026e-02	0.00000000e+00
$4P_{3/2}$	4	1	3/2	±1/2, ±3/2	-3.1250130010e-02	2.08015990e-07
$4D_{3/2}$	4	2	3/2	±1/2, ±3/2	-3.1250130010e-02	2.08015990e-07
$4D_{5/2}$	4	2	5/2	±1/2, ±3/2, ±5/2	-3.1250060670e-02	2.77355866e-07
$4F_{5/2}$	4	3	5/2	±1/2, ±3/2, ±5/2	-3.1250060670e-02	2.77355866e-07
$4F_{7/2}$	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)のようにあらわすことができる．$g_{nlj}(r)$は4成分波動関数のLarge成分(上側2成分)の動径方向の解に対応し，$f_{nlj}(r)$はSmall成分(下側2成分)の動径方向の解に対応する．$g_{nlj}(r)$と$f_{nlj}(r)$の大きさの比は，式(9)，(10)，(12)からある程度見積もることができ，その比は式(22)のようになる．水素原子の場合$Z=1$で，$\alpha$は微細構造定数で$\alpha\fallingdotseq 1/137$であるため，$g_{nlj}(r)$は$f_{nlj}(r)$よりも100倍以上大きな値を持つことが分かる($n$は1以上の整数)．

$$\begin{gather} \frac{g_{nlj}(r)}{f_{nlj}(r)} \sim \left(\frac{m_0 c^2 + E_{nj}}{m_0 c^2 - E_{nj}}\right)^{\frac{1}{2}} \sim \frac{2 n^2}{Z \alpha} \end{gather}$$

$n=1 \sim 4$の場合について，$g_{nlj}(r)$と$f_{nlj}(r)$を図4~20に図示する．ただし，$f_{nlj}(r)$は$1/Z\alpha$倍したものを表示している．また，グラフ作成に使用したプログラムは下記に折り畳んでいる．図4~20からも，$f_{nlj}(r)$は$1/Z\alpha$倍して$g_{nlj}(r)$と同程度のオーダーを持っていることが分かる．このことからも，4成分波動関数のうち，$g_{nlj}(r)$が寄与する上側2成分であるLarge成分が電子確率密度に大きく寄与することが確認できた．

動径方向波動関数を計算するプログラム(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. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$1S_{1/2}: n=1, l=0, j=1/2$）

図5. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$2S_{1/2}:n=2, l=0, j=1/2$）

図6. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$2P_{1/2}:n=2, l=1,j=1/2$）

図7. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$2P_{3/2}:n=2, l=1, j=3/2$）

図8. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$3S_{1/2}:n=3, l=0, j=1/2$）

図9. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$3P_{1/2}:n=3, l=1, j=1/2$）

図10. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$3P_{3/2}:n=3, l=1, j=3/2$）

図11. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$3D_{3/2}:n=3, l=2, j=3/2$）

図12. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$3D_{5/2}:n=3, l=2, j=5/2$）

図13. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4S_{1/2}:n=4, l=0, j=1/2$）

図14. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4P_{1/2}:n=4, l=1, j=1/2$）

図15. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4P_{3/2}:n=4, l=1, j=3/2$）

図16. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4D_{3/2}:n=4, l=2, j=3/2$）

図18. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4D_{5/2}:n=4, l=2, j=5/2$）

図19. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4F_{5/2}:n=4, l=3, j=5/2$）

図20. 動径波動関数：$g_{nlj}(r)$，$f_{nlj}(r)$ （$4F_{7/2}:n=4, l=3, j=7/2$）

動径確率密度分布について，相対論（ディラック方程式に従う電子）と非相対論（シュレディンガー方程式に従う電子）の比較をした．動径確率密度分布は相対論では式(23)のように定義され，非相対論では式(24)のように定義される．式(24)中の$R_{nl} (r)$は非相対論における動径波動関数である．非相対論における動径波動関数$R_{nl} (r)$は式(25)式のようになる．式(25)中の$L_{n+l}^{2l+1} (\zeta)$はラゲール陪多項式である．ラゲール陪多項式の定義は資料によって"揺れ"が存在するので((参考[6]))，特殊関数の計算ルーチンを使用する際は注意されたい．

$$\begin{align} P^{rela}_{nlj} (r) &= r^2 \left[\left\{g_{nlj}(r)\right\}^2 + \left\{f_{nlj}(r)\right\}^2 \right]\\ P^{non-rela}_{nl} (r) &= r^2 \left[R_{nl} (r)\right]^2 \\ R_{nl} (r) &= -\left\{ \left(\frac{2}{na_0} \right)^3 \frac{(n-l-1)!}{2n \left[ \left(n+l\right)!\right]^3} \right\}^{1/2} \cdot {\rm exp}\left(\frac{1}{2}\zeta\right) \cdot \zeta^l \cdot L_{n+l}^{2l+1} (\zeta) \end{align}$$

動径確率密度分布を非相対論と相対論で比較した結果を図21~26に示す．非相対論と相対論は表2のような対応付けをし，比較した．相対論では$j=l+1/2$と$j=l-1/2$の2種類の状態が存在するが，$j=l+1/2$となるときのみ非相対論との対応付けをしている．図21~26から非相対論と相対論の結果はほぼ一致していることが分かる．また，水素原子の場合，$\left[g_{nlj}(r)\right]^2 \gg \left[f_{nlj}(r)\right]^2$であるから，$\left[g_{nlj}(r)\right]^2 \fallingdotseq \left[R_{nl} (r)\right]^2$であることも分かる．また，計算に使用したプログラムは下記に置いている．

表2. 非相対論と相対論の対応表：動径方向

動径確率密度分布の比較(pythonプログラム)

動径確率密度分布の比較(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. 動径確率密度分布の比較，相対論：$n=1, l=0, j=1/2$ $(1S_{1/2})$，非相対論：$n=1, l=0$ $(1s)$

図22. 動径確率密度分布の比較，相対論：$n=2, l=0, j=1/2$ $(2S_{1/2})$，非相対論：$n=2, l=0$ $(2s)$

図23. 動径確率密度分布の比較，相対論：$n=2, l=1, j=3/2$ $(2P_{3/2})$，非相対論：$n=2, l=1$ $(2p)$

図24. 動径確率密度分布の比較，相対論：$n=3, l=0, j=1/2$ $(3S_{1/2})$，非相対論：$n=3, l=0$ $(3s)$

図25. 動径確率密度分布の比較，相対論：$n=3, l=1, j=3/2$ $(3P_{3/2})$，非相対論：$n=3, l=1$ $(3p)$

図26. 動径確率密度分布の比較，相対論：$n=3, l=2, j=5/2$ $(3D_{5/2})$，非相対論：$n=3, l=2$ $(3d)$

4.計算結果：球面スピノル

この節ではディラック方程式の角度方向の解について，下記4点を実施する．

• ディラック方程式の角度方向の解（球面スピノル）を可視化
• ディラック方程式とシュレディンガー方程式の角度方向の解を比較
• ディラック方程式とシュレディンガー方程式の電子確率密度の違いは，角度方向の解の違いに起因していることを確認
• ディラック方程式とシュレディンガー方程式の角度方向の解が一致する場合と一致しない場合について考察

4成分波動関数$\bm{\Psi}_{nljm} (r, \theta, \phi)$の角度成分$\bm{\tilde{\Omega}}_{jlm} (\theta, \phi)$を式(26)のように定義し，$\bm{\tilde{\Omega}}_{jlm} (\theta, \phi)$を式(27)で定義される曲面を描画することで球面スピノルを可視化した．Large成分，Small成分，upスピン，downスピンのそれぞれの寄与を見るために，式(27)中で使用される$r(\theta, \phi)$には式(28)~(34)で定義されるものを代入し，ある$j$，$l$，$m$に対し，合計6パターンの図を描画した．図27~35に$l=0,1,2$の値を取る場合について図示した．

$r(\theta, \phi)$の値として，式(32)~(33)のものを採用した場合，upスピンとdownスピンの位相が混ざるが，その位相が混ざった状態を上手く描画する方法が思いつかなかった．今回は位相の情報を描画することを諦めた．また，$a>0$としたとき，$m=a$と$m=-a$は位相の違いを無視すれば同一のグラフになるため，$m$の値については$m>0$のもののみ図示することにした．グラフの描画に使用したプログラムは下記に折り畳んで置いている．

$$\begin{gather} \bm{\tilde{\Omega}}_{jlm} (\theta, \phi) = \def\arraystretch{1.8} \begin{pmatrix} \omega^{(1)}_{jlm} (\theta, \phi)\\ \omega^{(2)}_{jlm} (\theta, \phi)\\ \omega^{(3)}_{jl'm} (\theta, \phi)\\ \omega^{(4)}_{jl'm} (\theta, \phi) \end{pmatrix} = \begin{pmatrix} \bm{\Omega}_{jlm} (\theta, \phi)\\ \bm{\Omega}_{jl'm} (\theta, \phi)\\ \end{pmatrix}\\ \begin{cases} x = r(\theta, \phi) {\rm sin}\theta {\rm cos}\phi \\ y = r(\theta, \phi) {\rm sin}\theta {\rm sin}\phi \\ z = r(\theta, \phi) {\rm cos} \phi \end{cases}\\ r^{\rm Large}_{\rm up} (\theta, \phi) = \sqrt{\left|\omega^{(1)}_{jlm} (\theta, \phi)\right|^2} \\ r^{\rm Large}_{\rm down} (\theta, \phi) = \sqrt{\left|\omega^{(2)}_{jlm} (\theta, \phi)\right|^2} \\ r^{\rm Small}_{\rm up} (\theta, \phi) = \sqrt{\left|\omega^{(3)}_{jl'm} (\theta, \phi)\right|^2} \\ r^{\rm Small}_{\rm down} (\theta, \phi) = \sqrt{\left|\omega^{(4)}_{jl'm} (\theta, \phi)\right|^2} \\ r^{\rm Large}_{\rm up+down} (\theta, \phi) = \sqrt{ \sum_{i = 1}^2{\left|\omega^{(i)}_{jlm} (\theta, \phi)\right|^2} } = \sqrt{\bm{\Omega}_{jlm}^\dagger (\theta, \phi) \bm{\Omega}_{jlm} (\theta, \phi) } \\ r^{\rm Small}_{\rm up+down} (\theta, \phi) = \sqrt{ \sum_{i = 3}^4{\left|\omega^{(i)}_{jl'm} (\theta, \phi)\right|^2} } = \sqrt{\bm{\Omega}_{jl'm}^\dagger (\theta, \phi) \bm{\Omega}_{jl'm} (\theta, \phi) } \\ \end{gather}$$

球面スピノルの描画に使用したプログラム(python)

球面スピノルの描画に使用したプログラム(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成分の球面スピノル$(S_{1/2})$：$l=0, j=1/2, m=1/2\ \ (j=l+1/2)$

図27-2. Small成分の球面スピノル$(S_{1/2})$：$l=0, j=1/2, m=1/2\ \ (j=l+1/2)$

図28-1. Large成分の球面スピノル$(P_{1/2})$：$l=1, j=1/2, m=1/2\ \ (j=l-1/2)$

図28-2. Small成分の球面スピノル$(P_{1/2})$：$l=1, j=1/2, m=1/2\ \ (j=l-1/2)$

図29-1. Large成分の球面スピノル$(P_{3/2})$：$l=1, j=3/2, m=1/2\ \ (j=l+1/2)$

図29-2. Small成分の球面スピノル$(P_{3/2})$：$l=1, j=3/2, m=1/2\ \ (j=l+1/2)$

図30-1. Large成分の球面スピノル$(P_{3/2})$：$l=1, j=3/2, m=3/2\ \ (j=l+1/2)$

図30-2. Small成分の球面スピノル$(P_{3/2})$：$l=1, j=3/2, m=3/2\ \ (j=l+1/2)$

図31-1. Large成分の球面スピノル$(D_{3/2})$：$l=2, j=3/2, m=1/2\ \ (j=l-1/2)$

図31-2. Small成分の球面スピノル$(D_{3/2})$：$l=2, j=3/2, m=1/2\ \ (j=l-1/2)$

図32-1. Large成分の球面スピノル$(D_{3/2})$：$l=2, j=3/2, m=3/2\ \ (j=l-1/2)$

図32-2. Small成分の球面スピノル$(D_{3/2})$：$l=2, j=3/2, m=3/2\ \ (j=l-1/2)$

図33-1. Large成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=1/2\ \ (j=l+1/2)$

図33-2. Small成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=1/2\ \ (j=l+1/2)$

図34-1. Large成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=3/2\ \ (j=l+1/2)$

図34-2. Small成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=3/2\ \ (j=l+1/2)$

図35-1. Large成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=5/2\ \ (j=l+1/2)$

図35-2. Small成分の球面スピノル$(D_{5/2})$：$l=2, j=5/2, m=5/2\ \ (j=l+1/2)$

ここで，非相対論(シュレディンガー方程式に従う電子)と相対論(ディラック方程式に従う電子)における電子確率密度がどのように記述されるかを整理し，非相対論と相対論の角度方向の解を比較する際に具体的に何を比較すれば意味のある議論ができるかを考える．

シュレディンガー方程式に従う電子波動関数を$\psi_{nlm}(r,\theta, \phi)$とすると，$\psi_{nlm}(r,\theta, \phi)$は式(34)のように動径成分$R_{nl}(r)$と角度成分$Y_l^m (\theta, \phi)$の積であらわせる．非相対論(シュレディンガー方程式)から導かれる電子確率密度$\rho^{\rm non-rela}_{nlm} (r,\theta, \phi)$は式(35)のようにあらわせる．

$$\begin{gather} \psi_{nlm}(r,\theta, \phi) = R_{nl}(r) Y_l^m (\theta, \phi) \\ \begin{split} \rho^{\rm non-rela}_{nlm} (r,\theta, \phi) &= \psi_{nlm}^\dagger(r,\theta, \phi) \psi_{nlm}(r,\theta, \phi) \\ &= \left[ R_{nl}(r) \right]^2 \left|Y_l^m (\theta, \phi) \right|^2 \end{split} \end{gather}$$

相対論(ディラック方程式)から導かれれる電子確率密度は$\rho^{\rm rela}_{nljm} (r,\theta, \phi)$は式(36)のようにあらわせる．ここで，式(36)中の式変形には$\left[g_{nlj}(r)\right]^2 \gg \left[f_{nlj}(r)\right]^2$であることと，$\left[g_{nlj}(r)\right]^2 \fallingdotseq \left[R_{nl} (r)\right]^2$を使用している．この記事では$n=3$までの場合についてしか，$\left[g_{nlj}(r)\right]^2 \fallingdotseq \left[R_{nl} (r)\right]^2$であることを具体的に示していないが(3節での議論)，水素原子のような原子核電荷の小さな原子系でかつ(水素原子の場合$Z=1$)，$j=l+1/2$となるような$j$を考えた場合，$\left[g_{nlj}(r)\right]^2 \fallingdotseq \left[R_{nl} (r)\right]^2$は一般に成り立つ．

$$\begin{gather} \begin{split} \rho^{\rm rela}_{nljm} (r,\theta, \phi) &= \bm{\Psi}^\dagger_{nljm}(r, \theta, \phi) \bm{\Psi}_{nljm}(r, \theta, \phi) \\ &= \left[g_{nlj}(r)\right]^2 \left(\sum_{i = 1}^2{\left|\omega^{(i)}_{jlm} (\theta, \phi)\right|^2}\right) + \left[f_{nlj}(r)\right]^2 \left(\sum_{i = 3}^4{\left|\omega^{(i)}_{jl'm} (\theta, \phi)\right|^2}\right) \\ &= \left[g_{nlj}(r)\right]^2 \left[r^{\rm Large}_{\rm up+down} (\theta, \phi)\right]^2 + \left[f_{nlj}(r)\right]^2 \left[r^{\rm Small}_{\rm up+down} (\theta, \phi)\right]^2 \\ &\fallingdotseq \left[g_{nlj}(r)\right]^2 \left[r^{\rm Large}_{\rm up+down} (\theta, \phi)\right]^2 \ \ \left(\because \left[g_{nlj}(r)\right]^2 \gg \left[f_{nlj}(r)\right]^2 \right)\\ &\fallingdotseq \left[R_{nl} (r)\right]^2 \left[r^{\rm Large}_{\rm up+down} (\theta, \phi)\right]^2 \ \ \left(\because \left[g_{nlj}(r)\right]^2 \fallingdotseq \left[R_{nl} (r)\right]^2 \right)\\ \end{split} \end{gather}$$

非相対論と相対論において電子確率密度にどのような違いが出るかを議論する際，式(35)，(36)から電子確率密度の違いは主にその角度成分$|Y_l^m (\theta, \phi)|^2$と$[r^{\rm Large}_{\rm up+down} (\theta, \phi)]^2$から生じていることが分かる．そこで，非相対論と相対論における角度成分の解の違いを比較する際は，非相対論と相対論とでそれぞれ式(37)，(38)のような$r(\theta, \phi)$を採用し，そのときの曲面の違いを議論する．

$$\begin{align} r(\theta, \phi) &= |Y_l^m (\theta, \phi)|\\ r(\theta, \phi) &= r^{\rm Large}_{\rm up+down} (\theta, \phi) \\ \end{align}$$

非相対論と相対論は表3のような対応付けをし，比較した．相対論では$j=l+1/2$と$j=l-1/2$の2種類の状態が存在するが，$j=l+1/2$となるときのみ非相対論との対応付けをしている．波動関数の角度方向の解を非相対論と相対論で比較した結果を図36~41に示す．

以前書いた記事では非相対論（シュレディンガー方程式）の角度方向の解を図示する際に，球面調和関数$Y_{l}^m (\theta, \phi)$を実数化した基底関数系を可視化したが(参考[7])，今回は相対論との比較をするために$Y_{l}^m (\theta, \phi)$を実数化せずにそのまま図示している．また，非相対論と相対論の曲面の違いが見にくいものもあったため，図36~41には斜め上の視点から見た場合と水平方向から見た場合の2パターンの視点で図示している．図36~41の図を作成するのに使用したプログラムは下記に折り畳んで置いている．

表3. 非相対論と相対論の対応表：角度方向

角度方向の解の比較(python)

角度方向の解の比較(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から，相対論における量子数$j$と$m$が$m=j$の関係を満たすとき，角度方向の解をプロットすると，非相対論と相対論のプロット結果は同一のものになっていることが分かる．これは，$m=j$となるときLarge成分の球面スピノル$\bm{\Omega}_{jlm} (\theta, \phi)$のdownスピン成分がゼロになり，upスピン成分が球面調和関数$Y_l^l(\theta, \phi)$そのものになり，$|Y_l^m (\theta, \phi)| = r^{\rm Large}_{\rm up+down} (\theta, \phi)$となるためである．このことは式(7)から分かる．式(7)を下に再掲する．

逆に，相対論における量子数$j$と$m$が$m \ne j$のとき，$|Y_l^m (\theta, \phi)| \ne r^{\rm Large}_{\rm up+down} (\theta, \phi)$となるため，それに対応する図36~41のプロット結果も異なったものになっていることが分かる．

$$\bm{\Omega}_{jlm} (\theta, \phi) = \begin{pmatrix} \sqrt{\frac{j+m}{2j}} Y_l^{m-\frac{1}{2}} (\theta, \phi) \\ \sqrt{\frac{j-m}{2j}} Y_l^{m+\frac{1}{2}} (\theta, \phi) \\ \end{pmatrix} \ \ \ (j = l + \frac{1}{2}) \tag{7}$$

図36. 角度方向の解の比較，非相対論：$(l=0,m=0)$，相対論：$(l=0,j=1/2,m=1/2)$

図37. 角度方向の解の比較，非相対論：$(l=1,m=0)$，相対論：$(l=1,j=3/2,m=1/2)$

図38. 角度方向の解の比較，非相対論：$(l=1,m=1)$，相対論：$(l=1,j=3/2,m=3/2)$

図39. 角度方向の解の比較，非相対論：$(l=2,m=0)$，相対論：$(l=2,j=5/2,m=1/2)$

図40. 角度方向の解の比較，非相対論：$(l=2,m=1)$，相対論：$(l=2,j=5/2,m=3/2)$

図41. 角度方向の解の比較，非相対論：$(l=2,m=2)$，相対論：$(l=2,j=5/2,m=5/2)$

5.計算結果：電子確率密度

この節では下記2点を実施する．

• ディラック方程式に従う場合の電子確率密度を可視化
• ディラック方程式とシュレディンガー方程式の電子確率密度の差異を確認

シュレディンガー方程式に従う電子の電子確率密度は式(35)のようにあらわせる．式(35)は下に再掲．ディラック方程式に従う電子の電子確率密度は式(39)のようにあらわせる．式(39)は式(36)に$[r^{\rm Large}_{\rm up+down} (\theta, \phi)]^2 = |\bm{\Omega}_{jlm} (\theta, \phi)|^2$と$[r^{\rm Small}_{\rm up+down} (\theta, \phi)]^2 = |\bm{\Omega}_{jl'm} (\theta, \phi)|^2$を代入して求めている．$\rho^{\rm non-rela}_{nlm} (r,\theta, \phi)$と$\rho^{\rm rela}_{nljm} (r,\theta, \phi)$は半透明の等値面を複数枚描画し，可視化した．詳細な可視化方法は参考[1]の記事を参照されたい．

$$\rho^{\rm non-rela}_{nlm} (r,\theta, \phi) = \left[ R_{nl}(r) \right]^2 \left|Y_l^m (\theta, \phi) \right|^2 \tag{35}$$

$$\begin{align} \rho^{\rm rela}_{nljm} (r,\theta, \phi) = \left[g_{nlj}(r)\right]^2 \left|\bm{\Omega}_{jlm} (\theta, \phi) \right|^2 + \left[f_{nlj}(r)\right]^2 \left|\bm{\Omega}_{jl'm} (\theta, \phi) \right|^2 \\ \end{align}$$

3節と4節で行った議論から下記2点が予想される．表4中で「〇」をつけているのは相対論と非相対論の電子確率密度がほぼ一致すると予想されもので，「×」をつけているのは電子確率密度が角度方向の解の分だけ違いが出ると予想されるものである．

• ディラック方程式(相対論)における，量子数$m$が$m=j$となるとき，相対論とそれに対応する非相対論(シュレディンガー方程式)の電子確率密度はほぼ一致する．
• ディラック方程式(相対論)における，量子数$m$が$m\ne j$となるとき，相対論とそれに対応する非相対論(シュレディンガー方程式)の電子確率密度は角度方向の解の分だけ違いがでる．

表4. 非相対論と相対論の比較表：電子確率密度

電子確率密度を図42~55に示す．上述したような規則で，ディラック方程式（相対論）とシュレディンガー方程式（非相対論）の電子確率密度の差異を確認できる．可視化に使用したプログラムは下記に折り畳んで置いた．

電子確率密度を可視化するためのプログラム(python)：非相対論

電子確率密度を可視化するためのプログラム(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)：相対論

電子確率密度を可視化するためのプログラム(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. 電子確率密度，非相対論：$(n=1,l=0,m=0)$，相対論：$(n=1,l=0,j=1/2,m=1/2)$

図43. 電子確率密度，非相対論：$(n=2,l=0,m=0)$，相対論：$(n=2,l=0,j=1/2,m=1/2)$

図44. 電子確率密度，相対論：$(n=2,l=1,j=1/2,m=1/2)$

図45. 電子確率密度，非相対論：$(n=2,l=1,m=0)$，相対論：$(n=2,l=1,j=3/2,m=1/2)$

図46. 電子確率密度，非相対論：$(n=2,l=1,m=1)$，相対論：$(n=2,l=1,j=3/2,m=3/2)$

図47. 電子確率密度，非相対論：$(n=3,l=0,m=0)$，相対論：$(n=3,l=0,j=1/2,m=1/2)$

図48. 電子確率密度，相対論：$(n=3,l=1,j=1/2,m=1/2)$

図49. 電子確率密度，非相対論：$(n=3,l=1,m=0)$，相対論：$(n=3,l=1,j=3/2,m=1/2)$

図50. 電子確率密度，非相対論：$(n=3,l=1,m=1)$，相対論：$(n=3,l=1,j=3/2,m=3/2)$

図51. 電子確率密度，相対論：$(n=3,l=2,j=3/2,m=1/2)$

図52. 電子確率密度，相対論：$(n=3,l=2,j=3/2,m=3/2)$

図53. 電子確率密度，非相対論：$(n=3,l=2,m=0)$，相対論：$(n=3,l=2,j=5/2,m=1/2)$

図54. 電子確率密度，非相対論：$(n=3,l=2,m=1)$，相対論：$(n=3,l=2,j=5/2,m=3/2)$

図55. 電子確率密度，非相対論：$(n=3,l=2,m=2)$，相対論：$(n=3,l=2,j=5/2,m=5/2)$

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)