量子力学では、境界条件を課すことで有限の体積における物理現象とそのエネルギー固有値を決めることができる。

例えば、長さ $L$ の一次元無限井戸型ポテンシャル中の平面波解は一般に
$\psi(x) = A e^{iEx}$

と表せる。
Dirichlet境界条件とは物理系の境界（$x=L$）で波動関数がゼロになるためにエネルギーが満たすべき条件のことである。
この条件を波動関数に課すと、

$Ae^{iEL}=0 \to E_n=\frac{2n\pi}{L}$
とエネルギー固有値を求めることができる。

無限井戸型ポテンシャルが半径 $R$ の球型になっている場合、極座標でシュレディンガー方程式を解かなければならない。今回はシュレディンガー方程式の解き方が分かっていることを前提に話を進める。

極座標シュレディンガー方程式の解は
$l$次の球ベッセル関数 $j_l(pr)$ 、球面調和関数 $Y^m_l (\theta,\phi)$を用いて

$\psi(x)=Nj_l(pr)Y^m_l(\theta,\phi)$
と表せる。
ここで $l$ は角運動量演算子 $L$ の固有値を構成する方位量子数であり、$m$ は磁気量子数、$p$は運動量の動径成分である。

Dirichlet境界条件を課すことで、運動量の動径成分も量子化することができる。そのためには

$j_l(pR)=0$

これを解析的に解くことはできないため、数値計算で具体的な値を求める必要がある。
猪木・河合著の「量子力学Ⅰ」によると、

$l=0$ のとき、$pR=3.14,6.28,9.42,…$
$l=1$ のとき、$pR=4.49,7.73,…$
$l=2$ のとき、$pR=5.76,9.10,…$
であった。
今回Pythonでの数値計算を通してこれらを確認し、そしてその他の状態における運動量の動径成分についても確認してみる。

今回用いるライブラリは以下の通りである。

import numpy as np 
import matplotlib.pyplot as plt
from scipy.special import *
from scipy.optimize import root_scalar

まずは、球ベッセル関数の様子を $l=0,1,2$ の範囲でプロットして確認する。

#x軸の描画
x = np.linspace(0,20,100)
y = np.zeros(100)
plt.plot(x,y,c='black')

#球ベッセル関数の描画
r = np.linspace(0,20,100)
for l in range(3):
    bessel = spherical_jn(l,r)
    plt.plot(r,bessel,label=f'l={l}')

plt.legend()
plt.show()

実行結果は以下の通りである。

これを見てみると、二つの仮説が立てられる。
1.$~~~j_l(pR)=0$ を満たす $pR$ の値は $1$ 以上離れている
2.$~~~j_l(pR)$は少なくとも $pR=l$ 以下で $0$ にならない。

したがって以下のような方針で解を求める。

解を格納するための空の配列momを用意する。また、point = 0 を定義する。
区間$[i,i+1]$ で $j_l(i)j_l(i+1)<0$ を満たすとき、root_scalar を用いて解を求め、 mom に格納する。また、ここで point に1加える。
$~j_l(i)j_l(i+1)<0$を満たさないとき、何もしない。
2,3が終わったあとで $i\to i+1$ とし、同じことをする。
point = 10 となった時点で解は10個集まったのでここで終了する。

とりあえず $l=10$ までの範囲で解を集める。コードは以下の通りとなる。

def f(l,x):
    return spherical_jn(l,x)

lm=10
R=1
SphBesselJZero_n = [] #各 l ごとに解をまとめて格納する。

for l in range(lm):
    #print(f'l:{l}')
    SphBesselJZero_ln = [] #同じ l に対して複数の解を格納する。
   
    point = 0
    i = l
    while(point<10):
        if(f(i)*f(i+1)<0):
            ans = root_scalar(f,bracket = [i,i+1])
            SphBesselJZero_ln.append(ans.root)
            point +=1
            
        i += 1
    SphBesselJZero_n.append(SphBesselJZero_ln)

mom = []
for l in range(lm):
    mom.append(SphBesselJZero_n[l])
mom

実行結果は以下の通りである。

array([[ 3.14159265,  6.28318531,  9.42477796, 12.56637061, 15.70796327,
        18.84955592, 21.99114858, 25.13274123, 28.27433388, 31.41592654],
       [ 4.49340946,  7.72525184, 10.90412166, 14.06619391, 17.22075527,
        20.37130296, 23.5194525 , 26.66605426, 29.81159879, 32.95638904],
       [ 5.7634592 ,  9.09501133, 12.32294097, 15.51460301, 18.68903636,
        21.85387422, 25.0128032 , 28.16782971, 31.32014171, 34.47048833],
       [ 6.987932  , 10.41711855, 13.69802315, 16.92362129, 20.12180617,
        23.30424699, 26.47676366, 29.64260454, 32.80373239, 35.9614058 ],
       [ 8.18256145, 11.70490715, 15.03966471, 18.30125596, 21.52541773,
        24.72756555, 27.9155762 , 31.09393321, 34.26539009, 37.43173677],
       [ 9.35581211, 12.96653017, 16.35470964, 19.6531521 , 22.90455065,
        26.12775014, 29.33256258, 32.52466129, 35.70757695, 38.88363096],
       [10.51283541, 14.20739246, 17.64797487, 20.98346307, 24.26276804,
        27.50786836, 30.73038073, 33.9371083 , 37.13233172, 40.31889251],
       [11.65703219, 15.43128921, 18.9229992 , 22.29534802, 25.60285595,
        28.87037335, 32.11119624, 35.33319418, 38.54136485, 41.73905287],
       [12.79078171, 16.64100288, 20.18247076, 23.59127482, 26.92704078,
        30.21726271, 33.47680082, 36.71452913, 39.93612781, 43.14542502],
       [13.91582261, 17.8386432 , 21.42848697, 24.87321392, 28.23713436,
        31.55018838, 34.82869654, 38.08247909, 41.31786469, 44.53914463]])

猪木・河合著の「量子力学Ⅰ」
$l=0$ のとき、$pR=3.14,6.28,9.42,…$
$l=1$ のとき、$pR=4.49,7.73,…$
$l=2$ のとき、$pR=5.76,9.10,…$
と比較しても良い結果になっていることが分かる。

今回は無限井戸型ポテンシャル中の境界条件を用いたいために関数を

$f(l,x)=j_l(x)$

としたが、用途によってはどのように変えても良い。

例えば MIT Bag modelではクォークは半径 $R$ の球の表面でカレントがゼロになるという条件から $j_l(pR)-j_{l+1}(pR)=0$ の解を求めることがある。
そのときには

$f(l,x)=j_l(x)-j_{l+1}(pR)$

とすればよい。