nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

EXAFS振動のフーリエ変換について

「EXAFS振動をフーリエ変換する」というフレーズはX線吸収分光をやっていると腐る程見かけるが、実際にはArtemisやらLarchやらのソフトで流し込んで、実際に気にするところはパラメータのフィッティングだけということがほとんどな気がする。
ここでは「そもそもフーリエ変換するとはどういうことか?」について触れて行きたいと思う。

EXAFSについては、田渕先生の以下のスライドが非常に分かりやすい。
http://titan.nusr.nagoya-u.ac.jp/Tabuchi/BL5S1/lib/exe/fetch.php?media=nssr-long-201406.pdf

EXAFS振動 \chi( k )を超簡略化して、一つのsin関数(かつ余計な項(位相シフト)を含まない)で与えられるとする。これは最も理想的なEXAFS振動である。
「理想的な」という意味は、フーリエ変換(波数 k \rightarrow 距離 R)すると、原子間距離(求めたい量) R_0のところで(原理的に)厳密に鋭いピークを持つ。(位相シフトなどゴチャゴチャ入ると、 R_0からピーク位置が少しズレる)

 \displaystyle
\chi( k ) = \sin{( 2 k R_0 )}

("2"が付いているのは距離 R_0 を電子が「行って返って来る」ことに依る)

これをいきなりフーリエ変換する前に、通常のやり方(信号解析(時間 t \rightarrow 周波数 f))を復習して、比較出来るようにする。
周波数 f_0の信号 s(t)がsin関数で次のように与えられているとする。

 \displaystyle
s( t ) = \sin{( 2 \pi f_0 t ) } = \sin{( \omega_0 t ) }

これをpythonを使ってフーリエ変換してみる。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

N = 600
delta_t = 1.0 / 600.0

min_t = 0.0
max_t = N * delta_t
t = np.linspace( min_t, max_t, N )

f0 = 10
s = np.sin( 2. * np.pi * f0 * t )

plt.plot( t, s, "-o" )
plt.savefig( "s.png" )
plt.close()

tilde_s = fft( s )
tilde_s = fftshift( tilde_s )
tilde_t = fftfreq( N, delta_t )
tilde_t = fftshift( tilde_t )

plt.plot( tilde_t, 1.0/N * np.abs( tilde_s ), "-o" )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

plt.plot( tilde_t, 2.0 * 1.0/N * np.abs( tilde_s ), "-o" )
plt.xlim( 0, 50 )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

f:id:koideforest:20180716015720p:plain
f:id:koideforest:20180716020903p:plain
f:id:koideforest:20180716020921p:plain

  •  f_0 = 10に対応したピークがキチンと得られている。
  • sin関数のフーリエ変換なので、負の振動数成分と対で現れる( cos関数も同様)。
  • 規格化は 1/N Nはメッシュ数(今の場合、 N = 600)。
  • そのため、正の振動成分のみをプロットする場合、二倍しないと一見規格化が保たれていないように見える。
  • フーリエ変換後の \tilde{t}の範囲は [ - 1 / 2 \Delta_t, \,  1 / 2 \Delta_t ](今の場合、 \Delta_t = 1/600、[ -300, 300 ])。
  • これに対し、 \Delta_{ \tilde{t} } = 1 / \Delta_t N (今の場合、 \Delta_{ \tilde{t} } = 1 )。
  • したがって、幾ら細かく( \Delta_t \ll 1 )データを取っても、元データのデータ範囲  N \Delta_tが変わらないとフーリエ変換後のメッシュは細かくならない。

単純にフーリエ変換すると、固有振動数 f_0の位置にピークが得られたわけだが、これはつまり横軸が振動数( \tilde{t} = f)に変換されており、角振動数 \omegaではない。
 \omegaに変換するには、以下の2パターンが考えられる。

1.  \omega_0 / 2 \pi の位置にピークが出て来ると解釈する。

 \displaystyle
s( t ) = \sin( \omega_0 t ) = \sin \left( 2 \pi \frac{\omega_0}{2\pi} t \right )
より f_0 = \omega_0 / 2 \pi と思える。

2. 横軸 \tilde{t} 2\piを掛けて、 \omegaに変えてしまう。( \omega = 2 \pi f


さて、ここでEXAFSの話に戻る。 \chi( k ) = \sin( 2 k R_0 ) = \sin( R_0 ( 2 k ) ) だったので、この形は \sin( \omega_0 t )に近い。
やりたいことは、「 \chi( k )フーリエ変換して、 R_0にピークが出る動径分布関数を得る」ことである。なので、

  1.  kではなく、 2kフーリエ変換する。
  2.  R_0 f_0ではなく \omega_0に対応しているので、横軸 \tilde{2k} 2\piを掛ける。

これを踏まえて、 \chi(k) から動径分布関数を得る。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

min_k = 2
max_k = 12
delta_k = 0.1
k = np.arange( min_k, max_k, delta_k )
N = len( k )

R0 = 3
k2 = 2. * k
delta_k2 = 2. * delta_k
chi = np.sin( R0 * k2 )
plt.plot( k, chi, "-o")
plt.savefig( "chi.png" )
plt.close()

tilde_chi = fft( chi )
tilde_chi = fftshift( tilde_chi )
tilde_k2 = fftfreq( N, delta_k2 )
tilde_k2 = fftshift( tilde_k2 )
R = 2. * np.pi * tilde_k2
plt.plot( R, 2.0 / N * np.abs( tilde_chi ), "-o" )
plt.xlim( 0, 10 )
plt.ylim( 0, 1 )
plt.grid()
plt.savefig( "tilde_chi.png" )
plt.close()

f:id:koideforest:20180716031948p:plain
f:id:koideforest:20180716032005p:plain

 R_0=3の付近にピークが得られているのがわかる。しかし、綺麗なピークとは言い難い。
これは kの範囲を広げると改善される。
 k = [ 2, 24 ]にすれば、
f:id:koideforest:20180716033511p:plain
f:id:koideforest:20180716033527p:plain

しかし、現実的には k>12の実験データを得ることは難しい。
そこで、「窓関数が掛かっている」と見做して、ゼロ値のデータを足したものをフーリエ変換する。
f:id:koideforest:20180716034843p:plain
f:id:koideforest:20180716034922p:plain

ピーク強度が小さくなり、周りに小さなピーク構造が現れたが、その代わりにデータ点が増えて滑らかな関数を表すようになった。
これが、Artemis等のソフトで行われている有限フーリエ変換である。窓関数を掛けるのはデータの一部「のみ」でフーリエ変換するのではなく、追加でゼロデータを足しまくって無理矢理滑らかなピークにすることの宣言でもある。
これによって、一見、十分なデータ数があるように見えてしまい、「fittingパラメータを増やすためには広い kの範囲と十分大きな Nが必要」と言われても全然しっくり来なくなってしまっているのである。

実態は、汚いピーク(少ないデータ点)を与えたフーリエ変換であり、それに対してfittingをすると思うと、「もっとたくさん(フーリエ変換後のピーク近傍の)データ点が欲しい!」と思えるのではないだろうか?

メタン(CH4)の分子軌道とLiebフラットバンド

四つ足で御馴染みのメタンの分子軌道が、固体におけるLiebフラットバンドが現れる構造と同じということで、テンションが上がった。

以下の簡単な飛び移りだけを考えたハミルトニアン行列の固有値固有ベクトルを求める。
対角成分(各原子上でのエネルギー)は簡単のためゼロにしてある。

import numpy as np
import sympy as sm

t = -1.0
h = np.zeros( ( 5, 5 ) )
for i in range( 5 ):
  for j in range( 5 ):
    if i == 0 and j != 0:
      h[ i ][ j ] = t
    elif i !=0 and j == 0:
      h[ i ][ j ] = t

# [ [ 0, t, t, t, t ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ], ]

H = sm.Matrix( h )
H.eigenvects()

#1. -2.0, [  2.0,  1.0, 1.0, 1.0, 1.0 ]
#2.  0.0, [  0, -1.0, 1.0, 0, 0 ]
#3.  0.0, [  0, -1.0, 0, 1.0, 0 ]
#4.  0.0, [  0, -1.0, 0, 0, 1.0 ]
#5.  2.0, [ -2.0,  1.0, 1.0, 1.0, 1.0 ]

今、遷移行列 tが負なので、1.の固有状態において節が無く広がっている(飛び移る)と安定になり、5.の固有状態において節が出来ると不安定になる様になっている。

Liebフラットバンドにおいて重要なのは、ゼロ固有値を持つのにヌルでない2., 3., 4.の状態である。
メタン的に言えば、単に「炭素上に節を持った状態」もしくは「非結合軌道」で終わってしまうところだが、もうちょっと言えば「水素上に波動関数が局在した状態」であり、「飛び移れるのに飛び移らない状態」が生じている。
今の計算は完全に有限系であるが、酸化物のReO3は正にLiebフラットバンドを持つ理想的な三次元構造を持っていて、バンド構造にもめちゃくちゃ平らなバンドが現れている。
でもパッとググったところ、誰も何も言っていないっぽい。まぁ普通にフェルミエネルギーがフラットバンドよりも上で金属だから、フラットバンド由来の物性が多分無いからだと思うが。。。

前に「飛び移りがあるのに飛び移らない状態があるのは普通なのか?」と悩んだことがあったので、自分の疑問が正当化された感じがした次第である。

立体角積分 - Lebedev quadrature

数値計算で立体角積分をしたいとき、安直に二重積分するよりも、Lebedev求積法を使った方が効率が良い。
Lebedev quadrature - Wikipedia
http://people.maths.ox.ac.uk/beentjes/Essays/QuadratureSphere.pdf

効率が良いという意味は、メッシュ数が少なくても精度が出るという意味で、言い換えれば同じ精度を出すのに二重積分よりもずっと速く計算が出来るということである。

求積法は、つまるところ、積分を「上手く精度が出るようにメッシュ間隔等を調整しかつ各点に重み付けをした和」に変換したものである。
Lebedev法の場合には、

 \displaystyle
\int sin \theta \, d\theta \, \int d\phi \, f(\theta,\phi) \approx 4 \pi \sum^N_{i=1} \omega_i \, f(\theta_i, \phi_i )

というように 4 \pi は後で掛けることを想定している。

Lebedev求積法のFortranのフリーのソースコードとしては、以下のものが知られている。
http://ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

subroutineの構成としては、

get_oh( ... )
LD0006( X, Y, Z, W, N )
LD0014( X, Y, Z, W, N )
...
LD5810( X, Y, Z, W, N )

という感じで、使うのはLDXXXX(...)の方である。get_ohは外から使うことは無い。

使い方のポイントを箇条書きすると、

  • 引数の"X, Y, Z, W, N"は全てoutputであるため、こちらから何か値を用意する必要は無い。
  • ただし、"X, Y, Z, W"は配列サイズを事前に指定する必要がある
  • 例えばLD0006(...)の"0006"はメッシュ数を表していて、これを使いたければ少なくとも6より大きい配列サイズの"X, Y, Z, W"を用意しておく。
  • "N"はメッシュ数であるためLD0006(...)ではN=6が返って来る。
  • "X Y Z"がLebedevメッシュであり、 ( \theta, \phi )に直したければ、"acos(z)"や"atan2(y,x)"で自分で変換する。
  • "W"が重み付けで、"sum( W )"はほぼ1になるようになっている。
  • どのLDXXXX(...)を使うべきかは、計算時間と計算精度との相談になる。

gfortranでコンパイルしてWの和が1になるかチェックする簡単なテストプログラム。

program check_Lebedev
implicit none
double precision, allocatable :: x(:), y(:), z(:), w(:)
! double precision, allocatable :: theta(:), phi(:)
integer :: n, n_lbdv

write( *, * ) 'number of mesh:'
read( *, * ) n
allocate( x(n), y(n), z(n), w(n) )
select case( n )
    case(    6 )
        call LD0006( x, y, z, w, n_lbdv )
    ...
    case( 5810 )
        call LD5810( x, y, z, w, n_lbdv )
    case default
        write( *, * ) ' n /= n_lbdv'
        stop
end select

write( *, * ) sum( w )

! allocate( theta( n_lbdv ), phi( n_lbdv ) )
! theta = acos( z )
! phi = atan2( y, x )

end program check_Lebedev

include "Lebedev-Laikov_CommentOutBy!.F"

リンク先のLebedev-Laikov.Fはコメントアウトが行頭の"c"で行われているため、sedか何かを使って"!"でコメントアウトするようにする必要がある(適当なコンパイルオプションを使えば、そんなことしなくても良いのかもしれないが)。

これで球面調和関数とかを掛けた計算とかが比較的簡単に出来るようになると思われる。

Møller wave operator における 加藤の連鎖律(?)

N粒子系における二体ポテンシャル Vを、粒子ペアに番号を付けて、その間の相互作用の和( M=N(N-1)/2)として定義する。

 \displaystyle
\begin{align}
V &\equiv \sum^{M}_{k=1} v_k = V_i + \bar{ V }_i = V_M = \bar{ V_0 }\\
V_i &\equiv \sum^{ i }_{k=1} v_k \\
\bar{ V }_i &\equiv \sum^{M}_{k= i + 1 } v_k
\end{align}

二体ポテンシャルとしたが、結局は相互作用に何かしら番号が付けられて、全体がその和として与えらえるのであれば、いろんな応用が考えられる(と思う)。
このような相互作用に対するN粒子系のGreen関数を以下のように定義。

 \displaystyle
\begin{align}
G &\equiv ( z - H_0 - V )^{-1} = G_M \\
G_i &\equiv ( z - H_0 - V_i )^{-1}
\end{align}

ここで次の量を考える。

 \displaystyle
\begin{align}
1 + G \bar{ V }_i &= G G^{-1} \left( 1 + G \bar{ V }_i \right)
    = G \left( G^{-1} + \bar{ V }_i \right) = G \left(  z - H_0 - V + \bar{ V }_i \right) \\
    &= G \left(  z - H_0 - V_i \right) = G \left( G_i \right)^{-1} \\
1 + G_i v_i &= G_i \left( G_i \right)^{-1} \left( 1 + G_i v_i \right)
    = G_i \left( \left( G_i \right)^{-1} + v_i \right) = G_i \left(  z - H_0 - V_i + v_i \right) \\
    &= G_i \left(  z - H_0 - V_{i-1} \right) = G_i \left( G_{i-1} \right)^{-1}
\end{align}

したがって、 G = G_Mに注意すると、

 \displaystyle
\begin{align}
G_M \left( G_i \right)^{-1} 
    &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \left( G_i \right)^{-1} \\
    &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \cdots \left( G_{i+1} \right)^{-1} G_{i+1} \left( G_i \right)^{-1}
\end{align}

すなわち、

 \displaystyle
\begin{align}
1 + G_M \bar{ V }_i
    &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_{i+1} v_{i+1} \right)
\end{align}

これを i = 0の場合に適用すれば、

 \displaystyle
\begin{align}
1 + G V &= 1 + G \sum^{M}_{k=1} v_k \\
    &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_1 v_1 \right)
\end{align}

という様に、和を積で表現し直すことが出来る。

 1 + GVはMøller wave operator  \Omegaとして知られており、 \Omega_i = 1 + G_i v_i と定義すれば、

 \displaystyle
\begin{align}
\Omega = \prod^{M}_{k=1} \Omega_k
\end{align}

と表せる。これを加藤の連鎖律と呼ぶらしい。

相互作用のペアを大きめに取っておいて、上手い具合に分割したら上手いこと相互作用を繰り込んだりして遊べそうな気がする。

球対称関数を別の位置で球平均する

 rを原点から測った距離とし、原点から見て球対称な関数を F( r )とする。
原点から見て位置 {\bf R}_iにあるサイト iがあるとする。
サイト iを原点に取り直した任意の位置ベクトルを {\bf r}_iと定義する。
やりたいことは、 F( r )をサイト iから見ると球対称ではないので、サイト iから見たときの球対称成分 f( r_i )を作りたいということである。
それはすなわちサイト i上で球平均することに対応する。
 {\bf r}_i = {\bf r} - {\bf R}_i、もしくは {\bf r} = {\bf r}_i + {\bf R}_i(こっちの方が直観的かも)であるから、

 \displaystyle
\begin{align}
f( r_i ) = \left( \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | ) \right) / \left( \int d \hat{\bf r}_i \right)
    = \frac{ 1 }{ 4 \pi } \int d \hat{\bf r}_i \, F( | {\bf r}_i + {\bf R}_i | )
\end{align}

球対称であるから、 {\bf R}_i z軸上に取っても結果は変わらない。
この時、 F( | {\bf r}_i + {\bf R}_i | )は方位角に依存しないから、

 \displaystyle
\begin{align}
f( r_i ) = \frac{ 1 }{ 2 } \int^{\pi}_{0} {\rm sin} \theta \, d \theta \, F( | {\bf r}_i + {\bf R}_i | )
\end{align}

ここで、 \theta {\bf r}_i - {\bf R}_iの成す角である。
マイナスが気になるかも知れないが、最終的にはその動径成分の R_iのみを使うので、符号はあまり重要でない。
 | {\bf r}_i + {\bf R}_i | = \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta }であることに注意すると、

 \displaystyle
\begin{align}
&\frac{ d | {\bf r}_i + {\bf R}_i | }{ d \theta }
    = \frac{ r_i R_i }{ \sqrt{ r^2_i + R^2_i - 2 r_i R_i {\rm cos} \theta } } {\rm sin} \theta
    = \frac{ r_i R_i }{ | {\bf r}_i + {\bf R}_i | } {\rm sin} \theta
\\
&\therefore \; {\rm sin} \theta \, d \theta
    = \frac{ | {\bf r}_i + {\bf R}_i | }{ r_i R_i  } \, d | {\bf r}_i + {\bf R}_i | 
\end{align}

よって、

 \displaystyle
\begin{align}
f( r_i ) &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, d | {\bf r}_i + {\bf R}_i | \, | {\bf r}_i + {\bf R}_i | \, F( | {\bf r}_i + {\bf R}_i | ) \\
    &= \frac{ 1 }{ 2 r_i R_i } \int^{ R_i + r_i }_{ R_i - r_i } \, dr \, r \, F( r )
\end{align}

 r_iの関数として f( r_i )を求めたいので、ここでの扱いは言わばパラメータのように外から与えられるものであり、積分変数ではないため積分の外に出せる。

例として、クーロンポテンシャル F( r ) = Z / r \equiv V( r )を考えると、

 \displaystyle
\begin{align}
V( r_i ) &= \frac{ Z }{ 2 r_i R_i } \left( R_i + r_i - ( R_i - r_i ) \right) = \frac{ Z }{ R_i } \equiv \bar{ V }
\end{align}

となり、違うサイトから見ると定数ポテンシャルとして与えられることになる。
したがって、点電荷が一次元に正負(絶対値は同じ大きさ Z)で交互に等間隔 Rで無限に並んでいて、ある点電荷に働くトータルの静電ポテンシャルは、球平均すれば

 \displaystyle
\begin{align}
V = - \frac{ 2 Z }{ R } \left( 1 - \frac{ 1 }{ 2 } + \frac{ 1 }{ 3 } - \cdots \right) = - \frac{ 2 Z }{ R } {\rm ln}2
\end{align}

と求めることが出来る。
ただ、これは ln( 1 + x )マクローリン展開、つまり x \approx 0に対して速く収束する展開に、 x = 1を代入して得られたものを使用しているわけだから、本当に点電荷がほぼ無限に並んでいないと(ほぼ無限に展開しないと)値の一致は悪いと思われる。クーロンポテンシャルの収束が遅いのがよく分かる。

デルタ関数のエンタングルメントを完備関係式で断つ?

例えば球面調和関数

 \displaystyle
\begin{align}
\delta( \hat{\bf r} - \hat{\bf r}' )
    &= \delta( \hat{\bf r}' - \hat{\bf r} ) \\
    &= \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) 
    = \sum_L Y^*_{L} ( \hat{\bf r} ) Y_L ( \hat{\bf r}' ) \\
    &= \sum_L Y_{L} ( \hat{\bf r}' ) Y^*_L ( \hat{\bf r} ) 
    = \sum_L Y^*_{L} ( \hat{\bf r}' ) Y_L ( \hat{\bf r} ) 
\end{align}

これはある意味 \hat{\bf r},\, \hat{\bf r} エンタングルメントを断ったと言えるか?
添字が一個なので、ベクトルに拡張すれば、

 \displaystyle
\begin{align}
\delta( \hat{\bf r} - \hat{\bf r}' )
 = \left(
        Y_{L_1} ( \hat{\bf r} ), Y_{L_2} ( \hat{\bf r} ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y^*_{L_1}( \hat{\bf r}' ) \\ Y^*_{L_2}( \hat{\bf r}' ) \\ \vdots
        \end{array}
    \right)
 = \left(
        Y^*_{L_1} ( \hat{\bf r} ), Y^*_{L_2} ( \hat{\bf r} ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r}' ) \\ Y_{L_2}( \hat{\bf r}' ) \\ \vdots
        \end{array}
    \right) \\
 = \left(
        Y_{L_1} ( \hat{\bf r}' ), Y_{L_2} ( \hat{\bf r}' ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y^*_{L_1}( \hat{\bf r} ) \\ Y^*_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
 = \left(
        Y^*_{L_1} ( \hat{\bf r}' ), Y^*_{L_2} ( \hat{\bf r}' ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
\end{align}

よって、ある種の行列化によって変数分離が可能である。ポイントは、もともと行列ではなくただの数だったので、「左に横ベクトル・右に縦ベクトル」の順番は変更出来ないという点である。
あまりエンタングルメントについて詳しくない(というか無知に等しい)が、エンタングルメントを断つ(解す)とは多分こういうことなのかと。
上では非局所相関を局所相関の積で表すというようなことをしたわけだが、この結果を用いると、いわゆる多重極展開と呼ばれる、動径成分と立体角成分に分ける展開が可能になる。
(これもある種のエンタングルメント解しか?)

 \displaystyle
\begin{align}
F( r, \hat{\bf r} )
    &= \int F( r, \hat{\bf r}' ) \delta( \hat{\bf r} - \hat{\bf r}' ) d \hat{\bf r}'
    = \int F( r, \hat{\bf r}' ) \sum_L Y_{L} ( \hat{\bf r} ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \\
    &= \sum_L \left( \int F( r, \hat{\bf r}' ) Y^*_L ( \hat{\bf r}' ) d \hat{\bf r}' \right) Y_{L}( \hat{\bf r} ) 
    \equiv \sum_L f_L( r ) Y_{L}( \hat{\bf r} ) \\
    &= \left(
        f_{L_1}( r ), f_{L_2}( r ), \cdots
    \right)
    \left(
        \begin{array}{c}
            Y_{L_1}( \hat{\bf r} ) \\ Y_{L_2}( \hat{\bf r} ) \\ \vdots
        \end{array}
    \right)
\end{align}

これは「連続量の相関を離散量の話に置き換える」と言えるだろうか?
ちなみに逆もおそらく可能(離散量の相関を連続量の話に置き換える)。

 \displaystyle
\begin{align}
\delta_{ L, L' }
    &= \int d \hat{\bf r} Y^*_L( \hat{\bf r} ) Y_{L'}( \hat{\bf r} ) \\
g_{N L}
    &= \sum_{L'} g_{NL'} \delta_{L' L}
    = \sum_{L'} g_{N L'} \int d \hat{\bf r} Y^*_{L'}( \hat{\bf r} ) Y_{L}( \hat{\bf r} ) \\
    &= \int d \hat{\bf r} \left( \sum_{L'} g_{N L'} Y^*_{L'}( \hat{\bf r} ) \right) Y_{L}( \hat{\bf r} )
    \equiv \int d \hat{\bf r} G_{N}( \hat{\bf r} ) Y_{L}( \hat{\bf r} )
\end{align}

連続と離散の絡み合いは楽しいので、もっと深めて行きたい。

Gaunt積分のまとめ

Gaunt積分の自分的メモ。
Clebsch–Gordan coefficients - Wikipedia
3-j symbol - Wikipedia

Gaunt積分(多分一般的な呼び名ではない。が、自分の分野ではよく出てくる。)
 \displaystyle
\begin{align}
 & \int d\hat{\bf r} \, Y^*_{ l_3 m_3 }(\hat{\bf r}) Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) \\
 & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) }{ 4 \pi ( 2 l_3 + 1 ) } }
        < l_1,0;l_2,0|l_3,0 > < l_1,m_1;l_2,m_2|l_3,m_3 >
\\
 & \quad = (-1)^{ m_3 } \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right)
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right)
\end{align}

Gaunt係数(Wikipediaとかに載っているのはこっち。)
 \displaystyle
\begin{align}
 & \int d\hat{\bf r} \, Y_{ l_1 m_1 }(\hat{\bf r}) Y_{ l_2 m_2 }(\hat{\bf r}) Y_{ l_3 m_3 }(\hat{\bf r}) \\
 & \quad = \sqrt{ \frac{ (2 l_1 + 1 ) ( 2 l_2 + 1 ) ( 2 l_3 + 1 ) }{ 4 \pi } }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{array} \right)
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right)
\end{align}

Clebsh-Gordan係数と3j記号の関係
 \displaystyle
\begin{align}
 < l_1,m_1;l_2,m_2|l_3,m_3 > =
        (-1)^{ l_1 - l_2 } (-1)^{ m_3 } \sqrt{ 2 l_3 + 1 }
        \left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & - m_3 \end{array} \right)
\end{align}

3j記号の列の交換
 \displaystyle
\begin{align}
\left( \begin{array}{ccc}  l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{array} \right)
    &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc}  l_2 & l_1 & l_3 \\ m_2 & m_1 & m_3 \end{array} \right)
\\
    &= (-1)^{ l_1 + l_2 + l_3 } \left( \begin{array}{ccc}  l_1 & l_3 & l_2 \\ m_1 & m_3 & m_2 \end{array} \right)
\end{align}