fortran66のブログ

fortran について書きます。

ベッセル関数 Jn が 陪ルジャンドル関数 Plm の極限である件

 \lim_{l\rightarrow\infty} l^{-m} P^m_l(\cos(x/l))=J_m(x)

[追記:Mathematica の Legendre 関数と位相因子の定義の違いがあったので直しました。また Fortran/Julia プログラムも修正しました。用語等も修正しました。]

前置き: ベッセル関数J_0 にむかつくw

その昔、物理数学で特殊関数を色々学んでいく過程で、第一種(円筒)ベッセル関数 J_n が、J_0 だけは \cos(x) っぽく 1 から始まるのに、残りの J_n (n\geq1)\sin(x) っぽく 0 から始まっているのにイライラしたものです。(下図参照)

J_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')

その後、球座標ばかり扱うようになって、「大きい J_n(x) お前はどうでもいい!小さい j_n(x) かわいいよ!」という気分になりました。(j_0(x) も一人だけ身勝手なんですが、\sin(x)/x なのでなんとなく許せます。)

しかし、ある日布団に寝そべって天井の模様を眺めながら、空想の中でいとしい球座標中の球体を上からペチャンコに潰してせんべい状にすると、球表面の球面調和関数の模様が、ベロンとせんべいの上に平面的に広がって、円筒座標上の模様となりました。

無ポテンシャル自由空間で解けば、円筒座標の動径方向解は第一種ベッセル関数 J_n(r) で、球座標での解は球ベッセル関数 j_n(r) と球面調和関数の積で、いま円筒座標の動径方向の小さい所を考えると、球の動径方向の依存性は無視できて(すなわち小さい j_n(r) ちゃんの出番なし)、 J_n(r)P_l^m(\theta) の何らかの極限になっていることが予想されます。

なお z 軸周りには、円筒座標・球座標とも \exp(im\phi) の依存性を持ち、共通なので考えなくて大丈夫です。

結論図

Wolfram Development Platform
軌道角運動量 l=100、磁気量子数 m=25 での Bessel 関数と associated Legendre 関数の比較。( J_{25}(x) ,  (-1)^{25} 100^{-25} P_{100}^{25}(\cos(x/100)) )

二つはおおむね重なり合います。 |m|\leq l/5 でおおむねピタリと重なります。ここでは、わざと違いがはっきりするようにしました。)

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}
と書かれます。これを微分方程式として解くと、球面調和関数が出てきます。つまり、陪ルジャンドル関数が出てきます。

ここで \hbar の値の大小は問題ではなく、単に座標系のスケールの取り方の問題に過ぎません。また、虚数単位も、演算子をエルミートにとるか反エルミートにとるか程度の差で、さほど本質的ではないと思います。物理では物理量演算子としてエルミートに取りますが、数学的には反エルミートに取った方がリー代数としてつじつまが合わせやすくなっています。

次に交換関係の物理的意味を考えてみます。物理的には角運動量は各軸周りの微小回転演算子になっているので、[l_x,l_y]=l_xl_y-l_yl_x=i\hbar\,l_z の関係式は、ある点をy軸まわりに少し回転させた後に新しいx軸まわりに少し回転させたものと、x軸まわりに少し回転させた後に新しいy軸まわりに少し回転させたものは、同じ位置に来ずに緯度は等しいが経度のズレた点に行くので、その差は z軸まわりの回転で埋め合わせられるという意味になっています。

これは角運動量の交換関係が、球体を暗示し球体の曲がり具合に関係していることを意味しているとみなせます。もしすべての交換関係が可換であれば、空間は各軸ごとに直交してユークリッド空間様になっているとみなせます。結局、リー代数の構造定数は、空間の曲がり具合の関係としても解釈できることを意味していると思われます。

さてここで、球座標中の球を z 軸方向に潰してペチャンコにする場合は、x, y 平面は平らになって並進対称性をもつようになり、z 軸方向の回転対称性は残ると考えられます。すると l_x,l_y は極限を取ることで l_x\rightarrow -p_y, l_y\rightarrow -p_x と回転軸と 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}
という形になると考えられます。

また、これは \hbar\rightarrow0 の古典論への極限の問題と考えると、l\rightarrow \inftyl_x,l_y だけが大きくなって古典論に近づく状況ともみなせると思います。

なお Wigner が 1950 年代に、より一般的に contraction という名で交換関係の半古典論への極限を考えています。また W. Miller がリー群・リー代数の観点から扱っています。

推測

以上を踏まえると、角運動量 l\rightarrow \infty の極限で、陪ルジャンドル関数が円筒ベッセル関数に近づくと考えられ、球面上の北極点から経線にそって降りてゆく円周上の距離が、球面がつぶれた時の平面上の距離になるはずだから、平面での距離 x は、球面では角運動量に規格化される形( \cos(x/l))で陪ルジャンドル関数に関係すると予想されます。

また、原点付近でのベッセル関数のべき依存性は x^m なのに対し、陪ルジャンドル関数は m が小さくなるにつれ全領域で全体が小さくなってゆくので l^m の係数がかかるとと予想されます。(まぁ先に図を描かせて合うまで適当にいじるんですけど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}

\theta=x/l とおいて、関数 F(x) に関する微分方程式とみなすと、

\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}
です。

ここで、l\rightarrow\infty\sin(x/l)\rightarrow x/l, \cos(x/l)\rightarrow1, l(l+1)\rightarrow l^2 の極限をとると、
\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}
となり、辺々に x^2/l^2 を掛けると、
\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}
となります。これは、丁度円筒ベッセル関数の微分方程式となっています。よって、F(x)J_m(x) であることが分かります。

変数を整えてやると、再び、
\begin{align} \lim_{l\rightarrow\infty} l^{-m} P^m_l(\cos(x/l))&=J_m(x)
\end{align}
を得ます。係数と位相因子は、まぁ適当にw

結論:ベッセル関数J_0(x) 許すw

こうして、第一種(円筒)ベッセル関数が陪ルジャンドル関数の極限である事が分かると、球対称下の波動関数の球面調和関数部分を思い浮かべることで、J_0(x) は磁気量子数 m=0 の成分で、\cos^l(x) 型の主要成分をもち、それ以外の磁気量子数の場合 \cos^{(l-m)}(x)\sin^m(x) 型の主要成分を持つ点に鑑みて、 J_0(x) だけが \cos(x) っぽくなっても仕方ないね!という気分になってきます。

また、下図に示すように、磁気量子数 m が軌道角運動量 l に対して大きくなると急激に一致が悪くなってゆくのは、m=l は、コマでいえば回転軸を垂直に立てて回っている状態で、lz が大きく lx, ly は小さい状態で、今回の近似の前提を完全に崩しており、m=0 は、コマが倒れかけて軸を横に寝かせてクルクルしている状況で、lz は小さく lx, ly が大きいと思えば、近似の前提に適っていて当然の気がします。

l=100, m=30

なお、この極限関係は、Mehler–Heine formula と呼ばれているらしく、幕末 1860 年代に知られていたようです。
Mehler–Heine formula - Wikipedia

Fortran 計算

Maple/Mathematica の類で計算すればいいのですが、それでは立つ瀬がないので Fortran でも計算しておきます。Bessel 関数は組み込み関数としてありますが、ルジャンドル関数は自分で用意する必要があります。本棚から埃のかぶった物理関数の本を取り出して、適当に漸化式で計算することにします。

ルジャンドル関数 P_l(x)

\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}

ルジャンドル関数 P_l(x)

\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}
角運動量l=150 位を超えるあたりから倍精度ではたりなくなるので4倍精度にする必要があります。

ルジャンドル関数の定義域は有界なので、ベッセル関の定義域は 0~{\pi\over2}l となります。l 無限大の極限で全域を表せるようになります。

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")