nano_exit

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

ポテンシャルのルジャンドル関数展開

量子力学 プログラミング

「ポテンシャルの角運動量展開」の二次元版だと思って頂いて差し支えない。

\displaystyle V( r, \theta ) = \sum_l V_l( r ) P_l( \cos\theta ) \\ \displaystyle V_l( r ) = \frac{ 2 l + 1 }{ 2 } \int^\pi_0 d\theta \, \sin\theta \, P_l( \cos\theta ) V( r, \theta )

係数 (2l+1)/2 はルジャンドル関数の直交関係から出て来る。
koideforest.hatenadiary.com

上記の式を用いて、 l_{max}=10まで展開してみる。

import numpy as np
from math import radians
from scipy.special import legendre
from scipy.integrate import simps
from scipy.interpolate import interp1d, interp2d
from matplotlib import pyplot as plt

def Plx( l, x ):
    # legendre( l ) >>> poly1d([ ... ])
    # poly1d([ ... ])( x ) >>> float
    return legendre( l )( x )

def pw_expand( func, l, r ):
    N = 300
    theta = np.linspace( 0., np.pi, N )
    coeff = ( 2 * l + 1 ) / 2.
    pw = []
    for rr in r:
        integrand = []
        for t in theta:
            x = rr * np.cos( t )
            y = rr * np.sin( t )
            integrand.append( coeff * np.sin( t ) * func( x, y ) * Plx( l, np.cos( t ) ) ) 
        print( np.array( integrand ) )
        pw.append( simps( np.array( integrand ), theta ) )
    return np.array( pw )

# For an interpolated function
def pw_expand_interp( func, l, r ):
    N = 300
    theta = np.linspace( 0., np.pi, N )
    coeff = ( 2 * l + 1 ) / 2.
    pw = []
    for rr in r:
        integrand = []
        for t in theta:
            x = rr * np.cos( t )
            y = rr * np.sin( t )
            integrand.append( coeff * np.sin( t ) * func( x, y )[0] * Plx( l, np.cos( t ) ) ) 
        pw.append( simps( np.array( integrand ), theta ) )
    return np.array( pw )

最初に、真四角の定数ポテンシャルを展開してみる。

N = 500
x = np.linspace( -2, 2, N )
y = np.linspace( -2, 2, N )
X, Y = np.meshgrid( x, y )

mat_v_const = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_const[ iy ][ ix ] = potential_const( xx, yy )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_v_const, cmap = 'Greys' )
fig.colorbar( pcm )
plt.savefig( "v_const.png" )
plt.show()

f:id:koideforest:20181217080205p:plain
v_const

角運動量は、球対称の時に保存量となるため、球対称に近い時に展開の収束が良いと言える。
これは、フーリエ変換の時に正弦波に近いほど収束が良いのと同じ話である。
そのため、不連続な境界を円周上に持つ真四角ポテンシャルの収束は、最悪だと予想出来る。
実際、全然収束しない。

r = np.linspace( 0, 2 * np.sqrt(2), 200 )
pws = []
for l in range( 10 + 1 ):
    pws.append( pw_expand( potential_const, l, r ) )

r = np.linspace( 0, 2 * np.sqrt(2), 200 )
f_const = [ interp1d( r, pws[ l ], kind = "cubic" ) for l in range( 10 + 1 ) ]
mat_pw_const = np.zeros( ( N, N ) )
l_max = 10
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        rr = np.linalg.norm([ xx, yy ])
        tt = np.arctan2( yy, xx )
        for l in range( l_max + 1 ):
            if l % 2 == 0:
                mat_pw_const[ iy ][ ix ] += f_const[l]( rr ) * Plx( l, np.cos( tt ) )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = ax.pcolormesh( X, Y, mat_pw_const, cmap = "bwr", vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( "pw_const_plot.png" )
plt.show()

f:id:koideforest:20181217080304p:plain
pw_const_plot

元のポテンシャルとの差を取ると、展開が上手く行っていない部分がより分かり易くなる。

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_pw_const - mat_v_const, cmap = "bwr", vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( "pw_const_plot_diff.png" )
plt.show()

f:id:koideforest:20181217080328p:plain
pw_const_diff_plot

各 lの成分に分解した V_l( r )を見てみると、 lが大きくなっても寄与が消える気配がない。

for ipw, pw in enumerate( pws ):
    plt.plot( r, pw, label = "{:d}".format( ipw ) )
plt.legend()
plt.savefig( "pw_const.png" )
plt.show()

f:id:koideforest:20181217080356p:plain
pw_const
真四角ポテンシャルが隅のパリティ（空間反転に対して対称）を持つため、 lが偶数の時だけ寄与する（ lが奇数の時、パリティは奇）。

次に、境界は球対称だが、中身はランダムポテンシャルで対称性もクソもない場合。（後でポテンシャルを再利用するため、cubic splineで補間して関数として取っておく。）

import random

def potential_random( x, y ):
    if x**2 + y**2 <=1:
        return random.random()
    else:
        return 0.

mat_v_rand = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_rand[ iy ][ ix ] = potential_random( xx, yy )

f_mat_v = interp2d( x, y, mat_v_rand, kind = "cubic" )
mat_v_interp = np.zeros( ( N, N ) )
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        mat_v_interp[ iy ][ ix ] = f_mat_v( xx, yy )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_v_interp, cmap = 'Greys' )
fig.colorbar( pcm )
plt.savefig( "v_rand.png" )
plt.show()

f:id:koideforest:20181217080418p:plain
v_rand

境界が丸いため、その部分は問題が無いが、やはり中身は苦労する。

pws_rand = []
for l in range( 10 + 1 ):
    pws_rand.append( pw_expand_interp( f_mat_v, l, r ) )

f_rand = [ interp1d( r, pws_rand[ l ], kind = "cubic" ) for l in range( 10 + 1 ) ]
mat_rand_pw = np.zeros( ( N, N ) )
l_max = 10
for ix, xx in enumerate( x ):
    for iy, yy in enumerate( y ):
        rr = np.linalg.norm([ xx, yy ])
        tt = np.arctan2( yy, xx )
        for l in range( l_max + 1 ):
            mat_rand_pw[ iy ][ ix ] += f_rand[l]( rr ) * Plx( l, np.cos( tt ) )

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_rand_pw, cmap = 'Greys', vmin = 0, vmax = 1 )
fig.colorbar( pcm )
plt.savefig( )
plt.show()

f:id:koideforest:20181217080443p:plain
pw_rand_plot

全然似ても似つかない。
差を取ると、何を再現出来ているのか不明なくらい差が大きい。

fig = plt.figure( figsize = ( 5, 4 ) )
ax = fig.add_subplot( 1, 1, 1 )
pcm = plt.pcolormesh( X, Y, mat_rand_pw - mat_v_interp, cmap = 'bwr', vmin = -1, vmax = 1 )
fig.colorbar( pcm )
plt.show()

f:id:koideforest:20181217080512p:plain
pw_rand_diff_plot

V_l(r)もほぼ全ての lが寄与を持つ。

for ipw, pw in enumerate( pws_rand ):
        plt.plot( r, pw, label = "{:d}".format( ipw ) )
plt.legend( loc = 1 )
plt.savefig( "pw_rand.png" )
plt.show()

f:id:koideforest:20181217080537p:plain
pw_rand

これらは、かなり非現実的なポテンシャルであるため、次回はもう少しマシなポテンシャルを展開してみたいと思う。

ボルン・オッペンハイマー近似

量子力学

Born-Oppenheimer近似の説明で、電子波動関数が原子核の位置に依存する部分の導入に違和感を感じていたので、自分なりにまとめる。

ハットは演算子が残っていることを表す。
例えば、電子の添え字を eとし、電子の位置を \vec{ r } = ( \vec{ r }_1, \vec{ r }_2, \cdots )とすると、
\displaystyle < \vec{r} | \hat{ r } | \Phi_e > = \vec{r} < \vec{r} | \Phi_e > = \vec{r} \Phi_e( \vec{ r } ) \\ \displaystyle < \vec{r} | \hat{T}_e | \Phi_e > = T_e ( \vec{r} ) \Phi_e( \vec{ r } ) = - \sum_i \frac{\hbar^2}{2m} \nabla_{r_i}^2 \Phi_e( \vec{ r } ) \\ \displaystyle < \vec{r} | \hat{V}_{ee} | \Phi_e > = V_{ee} ( \vec{r} ) \Phi_e( \vec{ r } ) = \sum_{i>j} \frac{ e^2}{| \vec{r}_i - \vec{r}_j |} \Phi_e( \vec{ r } )
のように、位置表示（位置を固定）等の表示を指定すると、具体的な形が設定されるとする。
例えば、運動量表示にすると、運動量エネルギー演算子は微分を含まず \hbar ^2k^2 / 2mというように表され、逆に位置演算子が運動量の微分を表すようになる。ここではそんなことはしないが、一応、断っておく。

原子核（添字： n）も合わせた、全ハミルトニアンを次のように表す。
\displaystyle \hat{H} = \hat{H}_{e} + \hat{H}_{n} + \hat{V}_{en} \\ \displaystyle \hat{H}_e = \hat{T}_e + \hat{V}_{ee} \\ \displaystyle \hat{H}_n = \hat{T}_n + \hat{V}_{nn}

したがって、シュレーディンガー方程式は、
\displaystyle \hat{H} | \Psi > = E | \Psi > \\ \left( < \vec{r} | \otimes < \vec{R} | \right) \hat{H} | \Psi > = H( \vec{ r }, \vec{ R } ) \Psi( \vec{r}, \vec{ R } ) = \left( H_{e}( \vec{r} ) + H_{n}( \vec{ R } ) + V_{en}( \vec{r}, \vec{ R } ) \right) \Psi( \vec{r}, \vec{ R } )
\vec{R} = ( \vec{R}_1, \vec{R}_2, \cdots ) は原子核の位置ベクトルである。

ここで、電子-原子核相互作用の原子核の位置ベクトル \vec{R}をパラメータにした \hat{V}^{R}_{en}を定義する。
\displaystyle < \vec{r} | \hat{V}^{R}_{en} | \Phi^{R}_{e} > = V^{R}_{en}( \vec{r} ) \, \Phi^{R}_{e}( \vec{r} ) = V_{en}( \vec{r}, \vec{R}) \, \Phi^{R}_{e}( \vec{r} )

個人的に重要なので、以下を強調する。
\displaystyle | \Phi^{R}_{e} > \neq | \Phi_{e} > \, \left( \hat{H}_e | \Phi_{e} > = E_e | \Phi_{e} > \right) \\ \displaystyle \hat{V}^{R}_{en} \neq \hat{V}_{en} \\ \displaystyle V^{R}_{en}( \vec{r} ) = V_{en}( \vec{r}, \vec{R})
つまり、演算子そのものは \hat{V}_{en}とは異なるが、パラメータを上手く原子位置に合わせると具体的な関数形が一致するような演算子 \hat{V}^{R}_{en}を作るのである。

こうすることで、 V_{en}の効果を取り込みながらも、基底空間として電子座標のみを考えたシュレーディンガー方程式が作れる。
\displaystyle < \vec{r} | \hat{H}^R_e | \Phi^R_{e} > \equiv < \vec{r} | \left( \hat{H}_e + \hat{V}^R_{en} \right) | \Phi^R_{e} > \\ \displaystyle \quad = \left( \hat{H}_e( \vec{r} ) + \hat{V}^R_{en}( \vec{r} ) \right) \Phi^{R}_{e}( \vec{r} ) = E^R_{e} \Phi^R_{e}( \vec{r} )

これにより、任意の関数 \Psi( \vec{r}, \vec{R} )はパラメータ R_0に依存した \Phi^{R_0}_{e}( \vec{r} )の線形結合で表せる。
\displaystyle \Psi( \vec{r}, \vec{R} ) = \sum_{k} C^{R_0}_k( \vec{R} ) \Phi^{R_0}_{e, k}( \vec{r} ) \\ \displaystyle H^{R_0}_e( \vec{r} ) \Phi^{R_0}_{e,k}( \vec{r} ) = E^{R_0}_{e,k} \Phi^{R_0}_{e,k}( \vec{r} )

追記：
この様に、電子系のハミルトニアンに対して対角的な電子波動関数による基底は「断熱的(adiabatic)電子基底」と呼ばれる。
当然、線形結合の係数を弄って（ユニタリー変換）、対角的でない電子波動関数の組で展開することも出来る。これは「透熱的(diabatic)電子基底」と呼ばれる。
透熱的な場合、上手く基底を選ぶことで、後述するBorn-Oppenheimer近似や断熱近似の影響を小さくすることが出来る。この辺りは次回以降に記述予定。

訂正：
電子位置 \vec{r}を受け止める状態空間は \Phi^{R_0}_{e, k}で描き切るので、残りの原子核位置 \vec{R}は展開係数が受けることになる。
\displaystyle | \Psi > = \sum_{k} | C^{R_0}_k > \otimes | \Phi^{R_0}_{e, k} > \\ \displaystyle \left( < \vec{r} | \otimes < \vec{R} | \right) | \Psi > = \sum_k C^{R_0}_k( \vec{R} ) \, \Phi^{R_0}_{e, k} ( \vec{r} )

これらの結果を、パラメータ R_0が元々の変数 \vec{R}に一致する様に調節して、全系のシュレーディンガー方程式に代入すれば、
\displaystyle H( \vec{ r }, \vec{ R } ) \Psi( \vec{r}, \vec{ R } ) = \left( H^{R}_{e}( \vec{r} ) + H_{n}( \vec{ R } ) \right) \Psi( \vec{r}, \vec{ R } ) \\ \displaystyle \quad = \left( H^{R}_{e}( \vec{r} ) + H_{n}( \vec{ R } ) \right) \sum_k C^R_k( \vec{R} ) \Phi^{R}_{e, k} ( \vec{r} ) = \sum_k \left( E^R_{e,k} + H_{n}( \vec{ R } ) \right) C^R_k( \vec{R} ) \Phi^{R}_{e, k} ( \vec{r} ) \\ \displaystyle \quad = E \sum_k C^R_k( \vec{R} ) \Phi^{R}_{e, k} ( \vec{r} )

ここで注意するべきは、 H_n( \vec{ R } ) の中の T_n( \vec{ R } )が微分演算を持っていて、 \Phi^{R}_{e, k} ( \vec{r} )にも作用する点である。
左から \Phi^{R*}_{e, k'} ( \vec{r} )を掛けて、 \vec{r}で積分すると、
\displaystyle \sum_k \left( T_{n,k'k}( \vec{ R } ) + \left( E^R_{e,k} + V_{nn}( \vec{ R } ) \right) \delta_{k'k} \right) C^R_k( \vec{R} ) = E C^R_{k'}( \vec{R} ) \\ \displaystyle T_{n,k'k}( \vec{ R } ) = \int d\vec{r} \, \left( \Phi^{R}_{e, k'}( \vec{r} ) \right)^* T_{n}( \vec{ R } ) \, \Phi^{R}_{e, k}( \vec{r} )

ここまで来れば、パラメータ Rを完全に変数と同一と扱っても問題無いだろう。
そうすれば、これまで単純に展開係数として扱っていた C^R_k( \vec{R} )が、正に原子核の波動関数 \Phi^{(k)}_{n}( \vec{R} )であることに気がつく。

\displaystyle \sum_k \left( T_{n,k'k}( \vec{ R } ) + V^{(k)}_{n}( \vec{ R } ) \delta_{k'k} \right) \Phi^{(k)}_n ( \vec{R} ) = E \Phi^{(k')}_n ( \vec{R} ) \\ \displaystyle V^{(k)}_{n}( \vec{R} ) = E_{e,k}( \vec{R} ) + V^{(k)}_{nn}( \vec{ R } )

T_{n,k'k}( \vec{ R } )の中身を詳しくみると、
\displaystyle T_{n,k'k}( \vec{ R } ) = - \sum_{i} \frac{\hbar^2}{ 2 M_i } \int d\vec{r} \, \left( \Phi_{e, k'}( \vec{r}, \vec{R} ) \right)^* \nabla_{R_i}^2 \, \Phi_{e, k}( \vec{r}, \vec{R} ) \\ \displaystyle \quad = - \sum_{i} \frac{\hbar^2}{ 2 M_i } \left( \nabla_{R_i}^2 \delta_{k'k} + T^{(1)}_{n,k'k}\nabla_{R_i} + T^{(2)}_{n,k'k} \right) \\ \displaystyle T^{(1)}_{n,k'k} = 2 \int d\vec{r} \, \left( \Phi_{e, k'}( \vec{r}, \vec{R} ) \right)^* \left( \nabla_{R_i} \, \Phi_{e, k}( \vec{r}, \vec{R} ) \right) \\ \displaystyle T^{(2)}_{n,k'k} = \int d\vec{r} \, \left( \Phi_{e, k'}( \vec{r}, \vec{R} ) \right)^* \left( \nabla_{R_i}^2 \, \Phi_{e, k}( \vec{r}, \vec{R} ) \right)

ここまでは、厳密である。
ここで、 T^{(1)} = T^{(2)} = 0と置くのがBorn-Oppenheimer近似である。
\displaystyle \left( - \sum_i \frac{\hbar^2}{ 2 M_i } \nabla_{ R_i }^2 + V^{(k)}_{n}( \vec{ R } ) \right) \Phi^{(k)}_n ( \vec{R} ) = E \Phi^{(k)}_n ( \vec{R} )
繰り返しておくと、 kは（多体）電子系の固有状態のラベルであり、例えば k=0を基底状態と思えば、電子系を基底状態に固定した時の原子核の動きをこの方程式から求めることが出来る。
基底状態に「固定する」とは言っても、 Rが変わるごとに基底状態の電子波動関数を解き直すので、電子波動関数が原子の振動等で変わらないという意味ではない。

wikipediaによれば、 T^{(i)}_{k'k} = T^{(i)}_{kk} \delta_{k'k}と置くのが断熱近似で、厳密にはBorn-Oppenheimer近似とは異なる。
断熱近似 - Wikipedia
この場合にも、異なる電子状態は混ざらない（電子状態のラベル kに対して対角的）から、原理的には電子状態を固定して原子核の運動を求めることが出来る。

上記でわかるとおり、Born-Oppenheimer近似やら断熱近似やらは、原子核の波動関数を求める時の話であって、原子核の位置を固定して電子波動関数を求めること自体は別段何も近似ではない。
電子波動関数の場合の「凍結フォノン近似」と呼ばれているものは、原子核が平衡位置にある時の電子波動関数だけを計算して満足するものであり、結果的に原子波動関数は求めないで済ます。

次回以降、Crude-Born-Oppenheimer or Crude adiabaticについても触れてみたい。

参考文献（ネットに落ちていたpdf）
http://vergil.chemistry.gatech.edu/notes/bo/bo.pdf

ポテンシャルの角運動量展開

量子力学

非球対称ポテンシャルでは、波動関数を球面調和函数で展開すると角運動量を添字とする行列になることを示した。
koideforest.hatenadiary.com
\displaystyle V_{LL'}( r ) = \int d\hat{r} Y^*_L( \hat{r} ) V( \vec{r} ) Y_L( \hat{r} )
角度積分するのに掛かる時間を t_I、必要な最大の軌道角運動量と動径メッシュの数をそれぞれ l_{max}、 N_rとすると、全ての行列要素を求める時間は t_I ( l_{max} + 1 )^4 N_rと表せる。
( l_{max} + 1 )^2は、例えば、 l_{max}=1まで必要とすると、 L = ( l, m ) = (0,0), (1,-1), (1,0), (1,1)で4つ必要であり、それを表せている。

次に、直接求めるのではなく、ポテンシャルを角運動量（球面調和関数）で展開することを考える。
\displaystyle V( \vec{r} ) = \sum_{L} V_L( r ) Y_L( \hat{r} ) \\ \displaystyle V_L( r ) = \int d\hat{r} \, Y^*_L( \hat{r} ) V( \vec{r} )
この展開に掛かる時間は t_I (l_{max}+1)^2 N_rと表せる。

これを使って、波動関数を角運動量展開した時の行列要素は、以下のように表せる。
\displaystyle V_{LL'}( r ) = \int d\hat{r} Y^*_L( \hat{r} ) V( \vec{r} ) Y_L( \hat{r} ) = \sum_{L''} V_{L''}( r ) \int d\hat{r} Y^*_L( \hat{r} ) Y_{L''}( \hat{r} ) Y_{L'}( \hat{r} ) \\ \displaystyle \quad = \sum_{L''} V_{L''}( r ) \, G( L'' L' | L )

G( L'' L' | L )はGaunt係数である。 l''は、最大で l,l'の2倍が必要であるため、Gaunt係数を求める時間を t_Gとすると、行列要素全てを求める時間は t_G (l_{max} + 1)^4 ( 2 l_{max} + 1 )^2 となる。
ポテンシャルを展開するところと合わせれば、 t_I (2l_{max}+1)^2 N_r + t_G (l_{max} + 1)^4 ( 2 l_{max} + 1 )^2と表せる。

したがって、
\displaystyle t_I (2l_{max}+1)^2 N_r + t_G (l_{max} + 1)^4 ( 2 l_{max} + 1 )^2 < t_I (l_{max}+1)^4 N_r \\ \displaystyle t_G < t_I N_r \frac{ (l_{max}+1)^4 - (2l_{max}+1)^2 }{ (l_{max} + 1)^4 ( 2 l_{max} + 1 )^2 } \\ \displaystyle t_G < t_I N_r \frac{ l_{max}^2 ( ( l_{max} + 2)^2 - 2 ) }{ (l_{max} + 1)^4 ( 2 l_{max} + 1 )^2 }
の時、Gaunt積分の方法の方が速くなる。

ただ、直感的には、どっちもあまり変わらなそうな気もする。

非球対称ポテンシャルにおけるT行列

量子力学

前回、非球対称ポテンシャルの時には、位相シフトがT行列を表すのにあまり役に立たないことを示した。
koideforest.hatenadiary.com

今回は、どうやってT行列（の行列成分）を求めるかを考える。
波動関数はT行列を使って次のように書けることを前回示した。
\displaystyle \psi_{\vec{k} }(\vec{r}) = \sqrt{\frac{2}{\pi}} \sum_{L'} \sum_{L} i^l \left( j_l(kr) \, \delta_{L'L} - ik \, h^+_{l'}(kr) \, \tilde{t}_{L'L} \right) Y_{L'}( \hat{r} ) Y^*_L( \hat{k} )

これは、ポテンシャルの外側の解である。
そこで、内側を含めた一般の動径波動関数を R_{L'L}( r ) \equiv R_{L'L}( k, r ) と定義すると、
\displaystyle \psi_{\vec{k} }(\vec{r}) = \sqrt{\frac{2}{\pi}} \sum_{L'} \sum_{L} i^l R_{L'L}( r ) Y_{L'}( \hat{r} ) Y^*_L( \hat{k} )

これをシュレーディンガー方程式に入れると、
\displaystyle \left( -\frac{1}{2}\frac{d^2}{d r^2} + \frac{ \hat{L}^2 }{ r^2 } + V( \vec{r} ) - \epsilon \right) \sqrt{\frac{2}{\pi}} \sum_{L'} \sum_{L} i^l r R_{L'L}( r ) Y_{L'}( \hat{r} ) Y^*_L( \hat{k} ) = 0 \\ \displaystyle \sum_{L'} \left( -\frac{1}{2}\frac{d^2}{d r^2} + \frac{ l' ( l' + 1 ) }{ r^2 } + V( \vec{r} ) - \epsilon \right) \sum_{L} i^l r R_{L'L}( r ) Y_{L'}( \hat{r} ) Y^*_L( \hat{k} ) = 0 \\ \displaystyle \sum_{L'} \left( \left( -\frac{1}{2}\frac{d^2}{d r^2} + \frac{ l' ( l' + 1 ) }{ r^2 } - \epsilon \right) \delta_{L_0 L'} + V_{L_0 L'}( r ) \right) \sum_{L} i^l r R_{L'L}( r ) Y^*_L( \hat{k} ) = 0 \\ \displaystyle V_{L_0 L'}( r ) \equiv \int d\hat{r} \, Y^*_{L_0}( \hat{r} ) V( \vec{r} ) Y_{L'}( \hat{r} )

ここで、 Y_{L_1}( \hat{k} )を掛けて、 \hat{k}で角度積分すると、
\displaystyle \sum_{L'} \left( \left( -\frac{1}{2}\frac{d^2}{d r^2} + \frac{ l' ( l' + 1 ) }{ r^2 } - \epsilon \right) \delta_{L_0 L'} + V_{L_0 L'}( r ) \right) i^{l_1} r R_{L' L_1}( r ) = 0

このシュレーディンガー方程式だけを眺める限りでは、「 L_1のラベルって意味あるの？」という感じしかしないが、外側の解を再度眺めると、 L' = L_1の時に裸の球ベッセル関数が入ることが要請されている。
原点近傍では、遠心力ポテンシャルの寄与が支配的であるため、自由解と同じ振る舞いをすると期待出来る。
そのため、微分方程式の初期条件として、 R_{L' L_1}( r \rightarrow 0) \rightarrow j_{l_1}(kr) \, \delta_{L' L_1}が課せられる。
つまり、 L_1は初期条件のラベルとして働いている。

したがって、行列表示で簡潔に表せば、
\displaystyle (H - E) \left( r \left( c_{L_1} \vec{R}_{L_1} \right) \right) = 0
L_1を固定した時に、定数倍を含めた波動関数ベクトル c_{L_1} \vec{R}_{L_1} が求まり、それを全ての L_1に対して実行すれば良いということになる。
数値計算上、地味に気を使わないといけないのは、定数倍の c_{L_1}の部分である。
微分方程式では定数倍しても方程式が変わらないため、例えば初期条件の値がずれると、得られた波動関数が解析解に比べて定数倍だけズレて出て来たりする。
これは、次にT行列を求めるところで、少し注意が要る部分だと思う。

波動関数 \tilde{R}_{L'L} \equiv c_L R_{L'L} が上記の方程式から得られたとする。
R_{L' L}が外側の解の表記と一致することから、
\displaystyle \tilde{R}_{L'L}( r ) = c_L \left( j_l(kr) \, \delta_{L'L} - ik \, h^+_{l'}(kr) \, \tilde{t}_{L'L} \right)
問題はここからどうやって \tilde{t}を抜き出すかということである。

まず、 c_Lを求めることを考える。
これにはWronskianが便利である。
詳細は省くが、球ベッセル関数と球ハンケル関数のWronskian W( j_l, h^+_l )は、解析的な関数 wで表すことが出来る。
\displaystyle W( j_l, h^+_l )( k, r ) = j_l \frac{ d h^+_l }{ d r } - \frac{ d j_l }{ d r } h^+_l = w( k, r )

したがって、
\displaystyle W( \tilde{R}_{LL}, h^+_l )( k, r ) = c_L W( R_{LL}, h^+_l )( k, r ) = c_L W( j_l, h^+_l )( k, r ) = c_L w( k, r ) \\ \displaystyle \therefore c_L = \frac{ W( \tilde{R}_{LL}, h^+_l )( k, r ) }{ w( k, r ) }

あとは、簡単な式変形から、
\displaystyle \tilde{t}_{L'L} = \frac{ \tilde{R}_{L'L}( r ) / c_L - j_l(kr) \, \delta_{L'L} }{ ik \, h^+_{l'}(kr) }

行列の微分方程式を解いたり、Wronskianを求めたり、面倒臭いことが多いが、逆に球対称な時がどれだけ簡単かを教えてくれている気がする。

散乱理論：位相シフトとT行列：非球対称ポテンシャル

量子力学

以前、球対称ポテンシャルの時の位相シフトとT行列についてまとめた。
koideforest.hatenadiary.com

今回は、非球対称の時に両者がどのように結ばれるかを調べる。

前回と同様、外側で値を持たないポテンシャルに対する外側の波動関数は、一般に次のように書ける。
\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_L i^l ( A_l \, j_l( kr ) + B_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad = \sqrt{ \frac{2}{\pi} } \sum_L i^l A_l( j_l( kr ) - \tan\tilde{\delta}_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \left( \frac{ 1 }{ (2\pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r} } = \sqrt{ \frac{2}{\pi} } \sum_L i^l j_l( kr )Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \right)

前回と同様、Lippman-Scwinger方程式を用いて解を表す（Hartree原子単位）。
\displaystyle \psi_{\bf k}( {\bf r} ) = \phi_{\bf k}( {\bf r} ) + \int d{\bf r}' \, G_0( {\bf r} - {\bf r}' ) V( {\bf r}' ) \psi_{\bf k}( {\bf r}' ) \\ \displaystyle \quad = \phi_{\bf k}( {\bf r} ) + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) \phi_{\bf k}( {\bf r}'' ) \\ \displaystyle \quad = \frac{ 1 }{(2 \pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r} } + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) \frac{ 1 }{(2 \pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r}'' } \\ \displaystyle \quad = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( {\bf r} ) \\ \displaystyle \qquad \qquad \qquad + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( \hat{\bf r} ) \\ \displaystyle \qquad \qquad - 2 i k \sum_{L'} \int d{\bf r}' d{\bf r}'' \, j_{l'}( kr_< ) Y^*_{L'}( \hat{\bf r}_< ) h^{(1)}_{l'}( kr_> ) Y_{L'}( \hat{\bf r}_> ) \\ \displaystyle \qquad \qquad \qquad \qquad \qquad \qquad \times T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} )

次に {\bf r}がポテンシャルの範囲の外側の時を考える。
ここまでの流れは、前回と完全に同じだが、ここから非球対称ポテンシャル、すなわち T行列が角運動量に対して非対角成分を持つ場合を考える。
\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( \hat{\bf r} ) \\ \displaystyle \qquad \qquad - 2 i k \sum_{L'} h^{(1)}_{l'}( kr ) Y_{L'}( \hat{\bf r} ) \int d{\bf r}' d{\bf r}'' \, j_{l'}( kr' ) Y^*_{L'}( \hat{\bf r}' ) \\ \displaystyle \qquad \qquad \qquad \qquad \qquad \qquad \times T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad \equiv \sqrt{ \frac{ 2 }{ \pi } } \sum_{L'}\sum_{L} i^l \left( j_l( kr ) \delta_{L'L} - i k h^{(1)}_{l'}( kr ) \, \tilde{ t }_{L'L}( k ) \right) Y_{L'}( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \tilde{ t }_{L'L}( k ) = 2 \int d{\bf r}' d{\bf r}'' \, j_{l'}( kr' ) Y^*_{L'}( \hat{\bf r}' ) T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'')

これを参考にして、外側の解の一般解を拡張する。
\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_{L'} \sum_L i^l A_{L'L}( j_l( kr ) - \tan\tilde{\delta}_{L'L} \, n_{l'}( kr ) ) Y_{L'}( \hat{\bf r} ) Y^*_L( \hat{\bf k} )

球対称の時、 A_{L'L} = A_l \, \delta_{L'L}、 \tan\tilde{\delta}_{L'L} = \tan\tilde{\delta}_l \, \delta_{L'L}となる。

したがって、 \tilde{\delta}_{L'L}と \tilde{t}_{L'L}( k )の関係は、
\displaystyle A_{L'L} = \delta_{L'L} - i k \tilde{t}_{L'L}( k ) \\ \displaystyle - A_{L'L} \tan \tilde{\delta}_{L'L} = k \tilde{t}_{L'L}( k ) \\ \displaystyle \therefore \tan \tilde{\delta}_{L'L} = - \frac{ k \tilde{t}_{L'L}( k ) }{ A_{L'L} } = i \frac{ i k \tilde{t}_{L'L}( k ) }{ \delta_{L'L} - i k \tilde{t}_{L'L}( k ) }

一見、球対称の時の自然な拡張っぽくなっていて、上手く行きそうだが、 L' \neq Lの時、
\displaystyle \tan \tilde{\delta}_{L'L} = i \frac{ i k \tilde{t}_{L'L}( k ) }{ - i k \tilde{t}_{L'L}( k ) } = -i
となり、非対角項を上手く記述することが出来ない。
これはある意味当たり前で、「元々の自由解は純粋な球ベッセル関数だけどポテンシャルの影響でどれくらい球ノイマンが混ざったか」を見るのが位相シフトであるにも関わらず、自由解では存在しない非対角項を無理矢理に球ベッセル関数と比較しようとしたところで、球ハンケル関数の性質がただ反映されるのがオチである。

したがって、T行列は対角項のみが位相シフトと関係を持ち、非対角項を位相シフトで表すことは出来ない。

永年方程式は無限空間には使えない。

量子力学

第一原理計算でよくある平面波展開は、よく考えると無限空間には使えない。
というのも、平面波の規格化が箱の規格化ではなく、デルタ関数規格化だからである。

\displaystyle \left( -\frac{ \nabla^2 }{ 2 } + V( \vec{r} ) - \epsilon \right) \psi( \vec{r} ) = 0

\psiを無限空間の平面波で展開すると、
\displaystyle \psi( \vec{r} ) = \frac{ 1 }{( 2 \pi )^{3/2}} \int d\vec{p} \, \tilde{ \psi }( \vec{p} ) e^{ i \vec{p} \cdot \vec{ r } } \\ \displaystyle \int d\vec{p} \, \left( \frac{ p^2 }{ 2 } + V( \vec{r} ) - \epsilon \right) \tilde{\psi}( \vec{p} ) e^{ i \vec{p} \cdot \vec{ r } } = 0

これの左から任意の平面波をかけて、位置で積分することで直交性から行列が上手いこと得られるというものであるが、
\displaystyle \frac{ 1 }{( 2 \pi )^3} \int d\vec{r} \int d\vec{p} \, \left( \frac{ p^2 }{ 2 } + V( \vec{r} ) - \epsilon \right) \tilde{\psi}( \vec{p} ) e^{ i \left( \vec{p} - \vec{p}_0 \right) \cdot \vec{ r } } = 0 \\ \displaystyle \int d\vec{p} \, \left( \left( \frac{ p^2 }{ 2 } - \epsilon \right) \delta( \vec{p} - \vec{p}_0 ) + V_{\vec{p}_0, \vec{p}} \right) \tilde{\psi}( \vec{p} ) = 0 \\ \displaystyle V_{\vec{p}_0, \vec{p}} = \frac{ 1 }{( 2 \pi )^3} \int d\vec{r} \, V(\vec{r}) e^{ i \left( \vec{p} - \vec{p}_0 \right) \cdot \vec{ r } }

周期系の結晶運動量の平面波の場合と比較すると、
\displaystyle \sum_{\vec{k}} \left( \left( \frac{ k^2 }{ 2 } - \epsilon \right) \delta_{\vec{k}_0, \vec{k} } + V_{\vec{k}_0, \vec{k}} \right) \tilde{\psi}( \vec{k} ) = \sum_{\vec{k}} M_{\vec{k}_0, \vec{k} } \tilde{\psi}_{\vec{k}} = 0 \\ \displaystyle V_{\vec{k}_0, \vec{k}} = \frac{ 1 }{L^3} \int d\vec{r} \, V(\vec{r}) e^{ i \left( \vec{k} - \vec{k}_0 \right) \cdot \vec{ r } }

周期系の場合には、 \delta_{\vec{k}_0, \vec{k} } が使われていて、値が発散することはないため、行列にして固有値を求めることが安全に出来るが、無限系の場合、 \delta( \vec{p} - \vec{p}_0 ) であるから、仮に行列にしたとしても対角成分が無限大に発散して、処理出来ない。

単振動を（完全）Green関数を使って解く。

数学 力学

これまでは以下の非摂動Green関数 G_0を使って来た。
koideforest.hatenadiary.com
\displaystyle \frac{ d^2 }{ d t^2 } x( t ) = f( t ) x( t ) \\ \displaystyle \frac{ d^2 }{ d t^2 } x_0( t ) = 0 \\ \displaystyle \frac{ d^2 }{ d t^2 } G_0( t, t' ) = \delta( t - t' ) \\ \displaystyle x( t ) = x_0( t ) + \int dt\, G_0( t, t' ) f( t' ) x( t' ) \\ \displaystyle \because \frac{ d^2 }{ d t^2 } x( t ) = \int dt' \, \delta( t - t' ) f( t' ) x( t' ) = f( t ) x( t )

また、T行列を使うと、積分の中を x_0を使って表すことが出来る。
koideforest.hatenadiary.com
\displaystyle T( t, t' ) = f( t ) \delta( t - t' ) + f( t ) G_0( t - t' ) f( t' ) \\ \displaystyle \qquad \qquad + \int dt'' \, f( t ) G_0( t - t'' ) f( t'' ) G_0( t'' - t' ) f( t' ) \cdots \\ \displaystyle f( t ) x( t ) = \int dt' \, T( t, t' ) x_0( t' ) \\ \displaystyle x( t ) = x_0( t ) + \int dt' dt''\, G_0( t, t' ) T( t', t'') x_0( t'' )

今回は、 f(t)を取り込んだ（完全） Green関数 Gを使う。
\displaystyle \left( \frac{ d^2 }{ d t^2 } - f( x ) \right) x( t ) = 0 \\ \displaystyle \left( \frac{ d^2 }{ d t^2 } - f( x ) \right) G( t, t' ) = \delta( t - t' ) \\ \displaystyle x( t ) = x_0( t ) + \int dt\, G( t, t' ) f( t' ) x_0( t' ) \\ \displaystyle \because \left( \frac{ d^2 }{ d t^2 } - f( x ) \right) x( t ) = -f( x ) x_0( t ) + \int dt\, \delta( t - t' ) f( t' ) x_0( t' ) \\ \displaystyle \quad = -f( x ) x_0( t ) + f( t ) x_0( t ) = 0 \\ \displaystyle \therefore G( t, t' ) = G_0( t, t' ) + \int dt_1 dt_2 \, G_0( t, t_1 ) T( t_1, t_2 ) G_0( t_2, t' )

単振動 f(t) = - \omega^2 のGreen関数は以前に既に求めてあるが、境界条件を定めていなかったため、それを含めて決定する。
境界条件は x(0) = l, \dot{x}( t ) = 0であるため、 G( 0, t' ) = \dot{G}( 0, t' ) = 0が課せられる。
koideforest.hatenadiary.com
\displaystyle G( t, t' ) = {\rm sgn}( t - t' ) \frac{ \sin( \omega( t - t' ) ) }{ 2 \omega } + a \cos( \omega t ) + b \sin( \omega t ) \\ \displaystyle G( 0 , t' ) = - \frac{ \sin( \omega( - t' ) ) }{ 2 \omega } + a = 0 \rightarrow a = - \frac{ \sin( \omega t' ) }{ 2 \omega } \\ \displaystyle \left. \frac{ d }{ d t } G( t , t' ) \right|_{ t = 0 } = - \frac{ \cos( \omega( - t' ) ) }{ 2 } + \omega b = 0 \rightarrow b = \frac{ \cos( \omega t' ) }{ 2 \omega } \\ \displaystyle \therefore G( t, t' ) = {\rm sgn}( t - t' ) \frac{ \sin( \omega( t - t' ) ) }{ 2 \omega } - \frac{ \cos( \omega t ) \sin( \omega t' ) }{ 2 \omega } + \frac{ \sin( \omega t ) \cos( \omega t' ) }{ 2 \omega } \\ \displaystyle \qquad = {\rm sgn}( t - t' ) \frac{ \sin( \omega( t - t' ) ) }{ 2 \omega } + \frac{ \sin( \omega ( t - t' ) ) }{ 2 \omega } \\ \displaystyle \qquad = \frac{ \sin( \omega( t - t' ) ) }{ \omega } \theta( t - t' )
とりあえず、遅延条件としてちゃんとGreen関数が求まった。
斉次解が今回求めたい解そのものなので、わざわざ Gを x(t)を解くのに必要なのかは疑問ではあるが、今は教育上必要ということで目を瞑る。

これを、積分の中に放り込めば、 f( t ) = - \omega^2より
\displaystyle x( t ) = l + \int^t_0 dt'\, \frac{ \sin( \omega( t - t' ) ) }{ \omega } ( - \omega^2 ) l = l - l \omega \int^t_0 dt'\, \sin( \omega( t - t' ) ) \\ \displaystyle \qquad = l + l \omega \int^{-t}_{0} dt''\, \sin( \omega( t + t'' ) ) \\ \displaystyle \qquad = l + l \omega \frac{1}{\omega} \left( - \cos( 0 ) + \cos( \omega t ) \right) \\ \displaystyle \qquad = l \cos( \omega t )

よって、ちゃんと求めることが出来た。