[追記:Mathematica の Legendre 関数と位相因子の定義の違いがあったので直しました。また Fortran/Julia プログラムも修正しました。用語等も修正しました。]
前置き: ベッセル関数 にむかつくw
その昔、物理数学で特殊関数を色々学んでいく過程で、第一種(円筒)ベッセル関数 が、 だけは っぽく 1 から始まるのに、残りの は っぽく 0 から始まっているのにイライラしたものです。(下図参照)
「 何故おまえは、好き好んで集団の和を乱してwひとりヨガリするのか!」
program Bessel implicit none integer, parameter :: n = 100 real , parameter :: pi = 4 * atan(1.0) real :: x(0:n), y(0:n, 0:5) integer :: i x = [(0.1 * i , i = 0, n)] do i = 0, 5 y(:, i) = Bessel_jn(i, x) end do do i = 0, n write(9, '(*(g0, 1x))') x(i), y(i, :) end do stop 'end' end program Bessel
#%fig: import matplotlib.pyplot as plt import numpy as np x, j0, j1, j2, j3, j4, j5 = np.loadtxt('fort.9', unpack = True) plt.title(r"Bessel function of the first kind $J_n$", fontsize=18) plt.xlabel(r"x", fontsize=15) plt.ylabel(r"$J_n(x)$", fontsize=15) plt.plot(x, j0, label=r"$J_0$") plt.plot(x, j1, label=r"$J_1$") plt.plot(x, j2, label=r"$J_2$") plt.plot(x, j3, label=r"$J_3$") plt.plot(x, j4, label=r"$J_4$") plt.plot(x, j5, label=r"$J_5$") plt.legend(fontsize=11) plt.grid() plt.savefig('Bessel')
その後、球座標ばかり扱うようになって、「大きい お前はどうでもいい!小さい かわいいよ!」という気分になりました。( も一人だけ身勝手なんですが、 なのでなんとなく許せます。)
しかし、ある日布団に寝そべって天井の模様を眺めながら、空想の中でいとしい球座標中の球体を上からペチャンコに潰してせんべい状にすると、球表面の球面調和関数の模様が、ベロンとせんべいの上に平面的に広がって、円筒座標上の模様となりました。
無ポテンシャル自由空間で解けば、円筒座標の動径方向解は第一種ベッセル関数 で、球座標での解は球ベッセル関数 と球面調和関数の積で、いま円筒座標の動径方向の小さい所を考えると、球の動径方向の依存性は無視できて(すなわち小さい ちゃんの出番なし)、 が の何らかの極限になっていることが予想されます。
なお z 軸周りには、円筒座標・球座標とも の依存性を持ち、共通なので考えなくて大丈夫です。
結論図
Wolfram Development Platform
軌道角運動量 l=100、磁気量子数 m=25 での Bessel 関数と associated Legendre 関数の比較。( , )
二つはおおむね重なり合います。 でおおむねピタリと重なります。ここでは、わざと違いがはっきりするようにしました。)
l=100 m=25 Plot[{BesselJ[m,x], (-1)^m*l^(-m)*LegendreP[l,m,Cos[x/l]]}, {x, 0, l*Pi/2}]
メモ帳
m= 10 Manipulate[Plot[{BesselJ[m,x], (-1)^m*l^(-m)*LegendreP[l,m,Cos[x/l]]}, {x, 0, l*Pi/2}], {l, m, 100}]
角運動量の交換関係の極限
さて少し迂遠ながら、角運動量の交換関係から考え始めたいと思います。量子力学で出てくる三次元空間での角運動量演算子の交換則は普通、
\begin{align} [l_x, l_y]&=i\hbar\, l_z, \\
[l_y, l_z]&=i\hbar\, l_x,\\
[l_z, l_x]&=i\hbar\, l_y
\end{align}
と書かれます。これを微分方程式として解くと、球面調和関数が出てきます。つまり、陪ルジャンドル関数が出てきます。
ここで の値の大小は問題ではなく、単に座標系のスケールの取り方の問題に過ぎません。また、虚数単位も、演算子をエルミートにとるか反エルミートにとるか程度の差で、さほど本質的ではないと思います。物理では物理量演算子としてエルミートに取りますが、数学的には反エルミートに取った方がリー代数としてつじつまが合わせやすくなっています。
次に交換関係の物理的意味を考えてみます。物理的には角運動量は各軸周りの微小回転演算子になっているので、 の関係式は、ある点をy軸まわりに少し回転させた後に新しいx軸まわりに少し回転させたものと、x軸まわりに少し回転させた後に新しいy軸まわりに少し回転させたものは、同じ位置に来ずに緯度は等しいが経度のズレた点に行くので、その差は z軸まわりの回転で埋め合わせられるという意味になっています。
これは角運動量の交換関係が、球体を暗示し球体の曲がり具合に関係していることを意味しているとみなせます。もしすべての交換関係が可換であれば、空間は各軸ごとに直交してユークリッド空間様になっているとみなせます。結局、リー代数の構造定数は、空間の曲がり具合の関係としても解釈できることを意味していると思われます。
さてここで、球座標中の球を z 軸方向に潰してペチャンコにする場合は、x, y 平面は平らになって並進対称性をもつようになり、z 軸方向の回転対称性は残ると考えられます。すると は極限を取ることで と回転軸と 90 度ずれた並進運動になるとして、交換関係は、
\begin{align} [p_y, p_x]&=0, \\
[p_x, l_z]&= i\hbar\,p_y,\\
[l_z, p_y]&= i\hbar\,p_x
\end{align}
という形になると考えられます。
また、これは の古典論への極限の問題と考えると、 で だけが大きくなって古典論に近づく状況ともみなせると思います。
なお Wigner が 1950 年代に、より一般的に contraction という名で交換関係の半古典論への極限を考えています。また W. Miller がリー群・リー代数の観点から扱っています。
推測
以上を踏まえると、角運動量 の極限で、陪ルジャンドル関数が円筒ベッセル関数に近づくと考えられ、球面上の北極点から経線にそって降りてゆく円周上の距離が、球面がつぶれた時の平面上の距離になるはずだから、平面での距離 x は、球面では角運動量に規格化される形( )で陪ルジャンドル関数に関係すると予想されます。
また、原点付近でのベッセル関数のべき依存性は なのに対し、陪ルジャンドル関数は が小さくなるにつれ全領域で全体が小さくなってゆくので の係数がかかるとと予想されます。(まぁ先に図を描かせて合うまで適当にいじるんですけどw)
これによって、数式ソフトに図を描かせると
\begin{align} \lim_{l\rightarrow\infty} (-1)^m l^{-m} P^m_l(\cos(x/l))&=J_m(x)
\end{align}
が推測されます。(上には Mathematica の例を示しました。)
微分方程式の極限
またもう少しまっとうに、微分方程式の極限から導くこともできます。
\begin{align}
\left[{1\over\sin\theta}{d\over d\theta}\left(\sin\theta {d\over d\theta}\right)-{m^2\over\sin^2(\theta)}+l(l+1)\right]
P^m_l(\cos\theta)&=0\\
\end{align}
とおいて、関数 に関する微分方程式とみなすと、
\begin{align}
\left[
{l\over\sin{x\over l}}
\left(
\cos{x\over l}{d\over dx}+ l \sin{x\over l}{d^2\over dx^2}
\right) -{m^2\over\sin^2{x\over l}}+l(l+1)
\right]
F(x)&=0\\
\end{align}
です。
ここで、 で の極限をとると、
\begin{align}
\left[
\left(
{l^2\over x}{d\over dx}+ l^2 {d^2\over dx^2}
\right)
+
\left( -{l^2m^2\over x^2}+l^2
\right)
\right]
F(x)&=0,\\
\end{align}
となり、辺々に を掛けると、
\begin{align}
\left[
x^2{d^2\over dx^2}+x{d\over dx}+
\left(x^2-m^2\right)
\right]
F(x)&=0\\
\end{align}
となります。これは、丁度円筒ベッセル関数の微分方程式となっています。よって、 は であることが分かります。
変数を整えてやると、再び、
\begin{align} \lim_{l\rightarrow\infty} l^{-m} P^m_l(\cos(x/l))&=J_m(x)
\end{align}
を得ます。係数と位相因子は、まぁ適当にw
結論:ベッセル関数 許すw
こうして、第一種(円筒)ベッセル関数が陪ルジャンドル関数の極限である事が分かると、球対称下の波動関数の球面調和関数部分を思い浮かべることで、 は磁気量子数 の成分で、 型の主要成分をもち、それ以外の磁気量子数の場合 型の主要成分を持つ点に鑑みて、 だけが っぽくなっても仕方ないね!という気分になってきます。
また、下図に示すように、磁気量子数 が軌道角運動量 に対して大きくなると急激に一致が悪くなってゆくのは、 は、コマでいえば回転軸を垂直に立てて回っている状態で、 が大きく は小さい状態で、今回の近似の前提を完全に崩しており、 は、コマが倒れかけて軸を横に寝かせてクルクルしている状況で、 は小さく が大きいと思えば、近似の前提に適っていて当然の気がします。
l=100, m=30
なお、この極限関係は、Mehler–Heine formula と呼ばれているらしく、幕末 1860 年代に知られていたようです。
Mehler–Heine formula - Wikipedia
Fortran 計算
Maple/Mathematica の類で計算すればいいのですが、それでは立つ瀬がないので Fortran でも計算しておきます。Bessel 関数は組み込み関数としてありますが、ルジャンドル関数は自分で用意する必要があります。本棚から埃のかぶった物理関数の本を取り出して、適当に漸化式で計算することにします。
ルジャンドル関数
\begin{align}
P_0(x)&=1,\\
P_1(x)&=x,\\
(l+1)P_{l+1}(x)&=(2l+1)xP_l(x)-nP_{n-1}(x),\,\,\,n\geq1.\\
\end{align}
陪ルジャンドル関数
\begin{align}
P_l^m(1)&=0 \,\,\,\, m\neq0,\\
P_l^0(1)&=1,\\
P_l^l(x)&=(2l-1)!!(1-x^2)^{l/2},\\
(l(l+1)-m(m -1) )P_l^{m -1}(x)&={2mx\over(1-x^2)^{1/2}}P_l^m(x)-P_l^{m+1}(x).
\end{align}
角運動量が 位を超えるあたりから倍精度ではたりなくなるので4倍精度にする必要があります。
ルジャンドル関数の定義域は有界なので、ベッセル関の定義域は となります。 無限大の極限で全域を表せるようになります。
l=100, m=20
!%fcflags: -O2 module m_legendre use, intrinsic :: iso_fortran_env implicit none integer, parameter :: kd = real64 contains pure elemental real(kd) function pn(n, x) integer , intent(in) :: n real(kd), intent(in) :: x integer :: i real(kd) :: p0, p1 p0 = 1.0_kd p1 = x select case(n) case(0) pn = p0 case(1) pn = p1 case(2:) do i = 1, n - 1 pn = ( (2 * i + 1) * x * p1 - i * p0 ) / (i + 1) p0 = p1 p1 = pn end do case default pn = 0 !0.0 / 0.0 ! NaN end select end function pn pure elemental real(kd) function dfact(n) integer, intent(in) :: n integer :: i dfact = 1.0_kd do i = n, 1, -2 dfact = dfact * i end do end function dfact pure elemental real(kd) function pn_assoc(n, m, x) integer , intent(in) :: n, m real(kd), intent(in) :: x if (abs(x) == 1.0_kd) then ! exceptional case |x| = 1 if (m == 0) then pn_assoc = 1.0_kd else pn_assoc = 0.0_kd end if else pn_assoc = dfact(2 * n - 1) * pn_a(n, m, x) end if end function pn_assoc pure elemental real(kd) function pn_a(n, m, xx) ! Plm integer , intent(in) :: n, m real(kd), intent(in) :: xx integer :: i real(kd) :: x, pm0, pm1 ! x = acos(abs(xx)) ! pm0 = sin(x)**n pm1 = 0.0_kd pm0 = sqrt(1.0_kd - xx*xx)**n if (m == n) then pn_a = pm0 else if (abs(m) <= n ) then do i = n, m + 1, -1 pn_a = ( 2 * i * abs(xx) / sqrt(1.0_kd - xx*xx) * pm0 - pm1 ) / ( n * (n + 1) - i * (i - 1) ) ! pn_a = ( 2 * i / tan(x) * pm0 - pm1 ) / ( n * (n + 1) - i * (i - 1) ) pm1 = pm0 pm0 = pn_a end do else pn_a = 0 !0.0 / 0.0 ! NaN end if pn_a = pn_a * sign(1.0_kd, xx)**(n-m) end function pn_a end module m_legendre program Legendre use m_legendre implicit none real, parameter :: pi = 4 * atan(1.0_kd) real(kd) :: x, y, t integer :: i, n, l, m l = 100 m = 20 ! n = 10 * l do i = 0, n x = i * 1.0d0 * l * pi /2 / n write(9, '(3(g0, 1x))') x, bessel_jn(m, x), real(l, kd)**(-m) * pn_assoc(l, m, cos(x / l)) ! check ! t = x * x ! y = pn(1, x) !2 print *, y - (3 * t - 1) / 2 !3 print *, y - x * (5 * t - 3) / 2 !4 print *, y - (35 * t**2 - 30 * t + 3) / 8 !5 print *, y - x * (63 * t**2 - 70 * t + 15) / 8 !6 print *, y - (231 * t**3 - 315 * t**2 + 105 * t - 5) / 16 !7 print *, x, y - x * (429 * t**3 - 693 * t**2 + 315 * t - 35) / 16 !8 print *, x, y - (6435 * t**4 - 12012 * t**3 + 6930 * t**2 - 1260 * t + 35) / 128 end do stop 'normal end' end program Legendre
注:mが大きい負の時有効桁が足りなくなるのでうまく処理できていませんw 正のmを反転させる処理が必要です。
#%fig: import matplotlib.pyplot as plt import numpy as np x, j0, j1 = np.loadtxt('fort.9', unpack = True) plt.title(r"$J_m(x)=\lim_{l\rightarrow\infty}l^{-m}P_l^m(\cos(x/l))$", fontsize=15) plt.xlabel(r"x", fontsize=15) plt.ylabel(r"y", fontsize=15) plt.plot(x, j0, label=r"$J_m$") plt.plot(x, j1, label=r"$P_l^m$") plt.legend(fontsize=11) plt.grid() plt.savefig("fort.png")
Julia
クソあばずれ julia は型変換ルールが謎過ぎて、高々10行程度の Fortran サブルーチンを移植するのに2時間くらいかかってしまいましたw 多倍長精度をやりたかったので、実数の型宣言をしないようにしたのですが、泥沼にはまりました。
using Plots using SpecialFunctions
function pn(n, x) #Legendre function $P_n(x)$ p0 = 1 p1 = x p2 = 0 if n == 0 return p0 elseif n == 1 return p1 elseif n > 1 for i = n-1:-1:1 p2 = ((2 * i + 1) * x * p1 - i * p0 ) / (i + 1) p0 = p1 p1 = p2 end return p2 else return 0 # out of domain end end #function pn function dfact(n::Int) # n!! double factorial s = 1.0 for i = n:-2:1 s *= i end return s end function pn_assoc(n::Int, m::Int, x) # associated Legendre function $P_n^m(x)$ factor = x >= 0 ? 1 : (-1)^(n+m) if (abs(x) == 1.0 && m == 0) return 1 elseif (abs(x)==1) return 0.0 else return dfact(2*n-1) * pn_a(n, m, x) * factor end end function pn_a(n::Int, m::Int, xx) # 0<=x<1, associated Legendre function $P_n^m(x)$ x = acos(xx) pm1 = 0.0 pm0 = sin(x)^n pm2 = 0.0 if m==n return pm0 elseif abs(m) <= n for i=n:-1:m+1 pm2 = ( 2 * i / tan(x) * pm0 - pm1 ) / ( n*(n+1) - i*(i-1) ) pm1 = pm0 pm0 = pm2 end return pm2 else return 0 # out of domain end end
l = 100 m = 10 x = 0:0.1:pi/2*l y1 = besselj.(m, x) y2 = y1*0 for i = 1:length(x) y2[i] = pn_assoc(l, m, cos(x[i]/l)) end plot(x, y1) plot!(x, y2) savefig("julia.png")