2018-10-12

3価のTiイオンのLII, LIII端における原子多重項遷移

量子力学

Ti ${}^{3+}$ の自由イオン状態の $L_{2,3}$ edgeは三本のピークが立つ。

Ti ${}^{3+}$ は基底状態で3dが空っぽなので、基底状態の多重項は ${}^1S_0$ となる。
（多重項は ${}^{2S+1}L_J$ で表される。ただし、 $L=0, 1, 2, 3, ...$ を $S, P, D, F, ...$ で表す。）

$L_{2,3}$ edgeでは、2pにホールが開いて、3dに電子が一個足されるので、

電子配置：2p ${}^5$ 3d ${}^1$

出現する多重項は、ホールと電子で同じ（エネルギー準位の順番は逆になる）
なので、2p ${}^1$ 3d ${}^1$ として捉えると、各軌道に電子が一個だけ入っているので、

スピン角運動量 $S_i$ ： $S_{2p} = 1/2$ ； $S_{3d} = 1/2$
軌道角運動量 $L_i$ ： $L_{2p} = 1$ ； $L_{3d} = 2$

イオン全体の角運動量はスピン、軌道それぞれで合成すれば良いので、

全スピン角運動量 $| S_{2p} - S_{3d} | \le S \le S_{2p} + S_{3d}$ ： $S = 0, 1$
全軌道角運動量 $| L_{2p} - L_{3d} | \le L \le L_{2p} + L_{3d}$ ： $L = 1, 2, 3$

この時点で、許される励起後の多重項は、
$\displaystyle {}^1P, {}^1D, {}^1F, {}^3P, {}^3D, {}^3F$

更に、スピン軌道相互作用が入って $S$ と $L$ が混ざるので、全角運動量 $|L-S| \le J \le L+S$ を合わせた多重項は、
$\displaystyle {}^1P_1, {}^1D_2, {}^1F_3, {}^3P_0, {}^3P_1, {}^3P_2, {}^3D_1, {}^3D_2, {}^3D_3, {}^3F_2, {}^3F_3, {}^3F_4$

よって、全角運動量 $J$ に対するに双極子遷移の選択則は $\Delta J=\pm1, 0$ （ただし、 $J=0 \rightarrow 0$ は禁制）なので、許容される遷移は
$\displaystyle {}^1S_0 \rightarrow {}^1P_1, {}^3P_1, {}^3D_1$
の三つとなる。

2018-10-12

Green関数の固有関数展開

数学量子力学

前回、フーリエ変換の視点で自由Green関数を弄った。
koideforest.hatenadiary.com

今回は、もっと一般的（？）にハミルトニアンの固有関数で展開することを考える。
$\displaystyle \left( E - H \right) G = 1 \\ \displaystyle \left( E - \epsilon \right) < \epsilon | G | \epsilon' > = < \epsilon | \epsilon' > \\ \displaystyle \left( E - \epsilon \right) G( \epsilon, \epsilon' ) = \delta_{ \epsilon, \epsilon' } \\ \displaystyle \therefore G^{\pm}( \epsilon, \epsilon' ) = \frac{ \delta_{ \epsilon, \epsilon' } }{ E - \epsilon \pm i \eta } = G^{\pm}( \epsilon ) \delta_{ \epsilon, \epsilon' } \\ \displaystyle ( \eta \rightarrow 0+ )$

したがって、Green関数の位置表示は、
$\displaystyle G^{\pm}( {\bf r}, {\bf r}' ) = < {\bf r} | G | {\bf r}' > \\ \displaystyle = \sum_{\epsilon} \sum_{\epsilon}' < {\bf r} | \epsilon > < \epsilon | G^{\pm} | \epsilon' > < \epsilon' | {\bf r}' > \\ \displaystyle = \sum_{\epsilon} \sum_{\epsilon}' \psi_\epsilon( {\bf r } ) G^{\pm}( \epsilon, \epsilon' ) \psi^*_{\epsilon'}( {\bf r }' ) \\ \displaystyle = \sum_{\epsilon} \sum_{\epsilon}' \psi_\epsilon( {\bf r } ) G^{\pm}( \epsilon ) \delta_{\epsilon, \epsilon' } \psi^*_{\epsilon'}( {\bf r }' ) \\ \displaystyle = \sum_{\epsilon} \psi_\epsilon( {\bf r } ) G^{\pm}( \epsilon ) \psi^*_{\epsilon}( {\bf r }' ) \\ \displaystyle = \sum_{\epsilon} \frac{ \psi_\epsilon( {\bf r } ) \psi^*_{\epsilon}( {\bf r }' ) }{ E - \epsilon \pm i \eta }$

自由電子の時には、運動エネルギーの固有状態表示で展開すれば良く、その時には、 $\psi_{\bf k}( {\bf r } ) \psi^*_{\bf k}( {\bf r }' ) \propto e^{ i {\bf k} \cdot ( {\bf r} -{ \bf r}' ) }$ のように、相対位置にのみ依存するようになる。

2018-10-12

自由電子Green関数のFourier変換

数学量子力学

パッと探したところ $G_0( {\bf r}, {\bf r}' ) = G_0( {\bf r} - {\bf r}' )$ と出来る理由が見当たらなかったので、真面目に計算してみた。
（もちろん、ハミルトニアンの並進対称性から明らかなのはわかっているが、真面目な方法でやって同じ答えになるはずである）

Diracのbraket表示を使って、演算子の方程式から始める。
$\displaystyle \left( E - H_0 \right) G_0 = 1 \\ < {\bf r} | \left( E - H_0 \right) G_0 | {\bf r}' > = < {\bf r} | {\bf r}' > \\ \left( E - H_0( {\bf r} ) \right) < {\bf r} | G_0 | {\bf r}' > = < {\bf r} | {\bf r}' > \\ \left( E - H_0( {\bf r} ) \right) G_0( {\bf r}, {\bf r}' ) = \delta( {\bf r} - {\bf r}' )$

Hartree原子単位系での（位置表示での）ハミルトニアン。
$\displaystyle H_0 = - \frac{ \nabla^2 }{ 2 }$

一歩前に戻って、運動量基底の完備関係式を挟む。

$\displaystyle \int d{\bf k} d{\bf k}' \, \left( E - H_0( {\bf r} ) \right) < {\bf r} | {\bf k} > < {\bf k} | G_0 | {\bf k}' > < {\bf k}' | {\bf r}' > \\ \qquad = \int d{\bf k} < {\bf r} | {\bf k} > < {\bf k} | {\bf r}' > \\ \displaystyle \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, \left( E - H_0( {\bf r} ) \right) e^{ i {\bf k} \cdot {\bf r} } G_0( {\bf k}, {\bf k}' ) e^{ - i {\bf k}' \cdot {\bf r}' } = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} \, e^{ i {\bf k} \cdot {\bf r} } e^{ - i {\bf k} \cdot {\bf r}' } \\ \displaystyle \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, \left( E - \frac{k^2}{2} \right) e^{ i {\bf k} \cdot {\bf r} } G_0( {\bf k}, {\bf k}' ) e^{ - i {\bf k}' \cdot {\bf r}' } = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} \, e^{ i {\bf k} \cdot {\bf r} } e^{ - i {\bf k} \cdot {\bf r}' } \\ \displaystyle \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, \left( E - \frac{k^2}{2} \right) e^{ i {\bf k} \cdot ( {\bf r} - {\bf r}' ) } G_0( {\bf k}, {\bf k}' ) e^{ i ( {\bf k} - {\bf k}' ) \cdot {\bf r}' } = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} \, e^{ i {\bf k} \cdot ( {\bf r} - {\bf r}' ) } \\ \displaystyle \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, \left( \left( E - \frac{k^2}{2} \right) G_0( {\bf k}, {\bf k}' ) e^{ i ( {\bf k} - {\bf k}' ) \cdot {\bf r}' } \right) e^{ i {\bf k} \cdot ( {\bf r} - {\bf r}' ) } \\ \displaystyle \qquad = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \delta( {\bf k}- {\bf k}' ) \, e^{ i {\bf k} \cdot ( {\bf r} - {\bf r}' ) }$

したがって、
$\displaystyle \left( E - \frac{k^2}{2} \right) G_0( {\bf k}, {\bf k}' ) e^{ i ( {\bf k} - {\bf k}' ) \cdot {\bf r}' } = \delta( {\bf k}- {\bf k}' )$

解析性を指定する為に、極を僅かに外せば、
$\displaystyle G^{\pm}_0( {\bf k}, {\bf k}' ) = \frac{ \delta( {\bf k} - {\bf k}' ) e^{ i ( {\bf k} - {\bf k}') \cdot {\bf r}' } }{ E - k^2 / 2 \pm i \eta } = \frac{ \delta( {\bf k} - {\bf k}' ) }{ E - k^2 / 2 \pm i \eta } = G^{\pm}_0( {\bf k} ) \delta( {\bf k} - {\bf k}' ) \\ \displaystyle ( \eta \rightarrow 0+ )$

このデルタ関数が、 $G_0$ のフーリエ成分が ${\bf k}$ に対して対角的であることを示している。
これを使って、逆フーリエ変換で戻せば、 $G_0$ が相対位置にしか依存しないことがわかる。
$\displaystyle G^{\pm}_0( {\bf r}, {\bf r}' ) = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, G^{\pm}_0( {\bf k}, {\bf k}' ) e^{ i {\bf k} \cdot {\bf r} } e^{ - i {\bf k}' \cdot {\bf r}' } \\ \displaystyle = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} d{\bf k}' \, G^{\pm}_0( {\bf k} ) \delta( {\bf k} - {\bf k}' ) e^{ i {\bf k} \cdot {\bf r} } e^{ - i {\bf k}' \cdot {\bf r}' } \\ \displaystyle = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} \, G^{\pm}_0( {\bf k} ) e^{ i {\bf k} \cdot {\bf r} } e^{ - i {\bf k} \cdot {\bf r}' } \\ \displaystyle = \frac{ 1 }{ ( 2 \pi )^3 } \int d{\bf k} \, G^{\pm}_0( {\bf k} ) e^{ i {\bf k} \cdot ( {\bf r} - {\bf r}' ) } \\ \displaystyle = G^{\pm}_0( {\bf r} - {\bf r}' )$

2018-10-10

XMCDのSum ruleをpythonで処理

プログラミング固体物理

とあるスクールで、講師が事前に用意したOriginのマクロを使ってXMCDのSum ruleから磁気モーメントを求めるという講習があったが、Originを使いたいと全く思わなかったため、一人でpython スクリプトを黙々と作っていた。

ローレンツ関数やら誤差関数やらで適当に作った架空のXASとXMCD。
f:id:koideforest:20181010044729p:plain
これを処理してSum ruleを適用する。

手順としては

pre edgeを差し引く。
引いたものをpost edgeで割る。
XMCDを（再度）作る。
ピークの位置を求める。
step関数を作る。
step関数をXASから引いて積分する。
$L_3$ と $L_2$ の境目を決めてXMCDを積分する。

という感じだろう。

pre edge、post edgeを以下の関数を使って、求める。

import numpy as np

def pre_edge( x, y, start, end ):
    pre_x = []
    pre_y = []
    for ix, xx in enumerate( x ):
        if start < xx:
            pre_x.append( xx )
            pre_y.append( y[ ix ] )
        if end < xx:
            break
    pre_x = np.array( pre_x )
    pre_y = np.array( pre_y )
    A = np.c_[ pre_x, np.ones( len( pre_x ) ) ]
    a, b = np.linalg.lstsq( A, pre_y, rcond=None )[0]
    return a * x + b

def post_edge( x, y, start, end ):
    post_x = []
    post_y = []
    for ix, xx in enumerate( x ):
        if start < xx:
            post_x.append( xx )
            post_y.append( y[ ix ] )
        if end < xx:
            break
    post_x = np.array( post_x )
    post_y = np.array( post_y )
    A = np.c_[ post_x, np.ones( len( post_x ) ) ]
    a, b = np.linalg.lstsq( A, post_y, rcond=None )[0]
    return a * x + b

要は、各エネルギー領域を指定してあげて、その中のデータ点を用いて最小二乗法で直線を定めるというもの。

今の場合だと、

pre edge: 680 ~ 700 eV
pre edge: 760 ~ 800 eV

くらいで設定した。
とりあえずプラスhelicityのXASにのみ適用してみる。

# energy: array for energy points
# spec_p: array for XAS spectrum with plus helicity

pre   = pre_edge( energy, spec_p, 680, 700 )
norm0 = spec_p - pre

post = post_edge( energy, norm0, 760, 800 )
norm_p = norm0 / post

fig = plt.figure( figsize=( 8, 2.5 ) )
ax1 = fig.add_subplot( 1, 3, 1 )
ax2 = fig.add_subplot( 1, 3, 2 )
ax3 = fig.add_subplot( 1, 3, 3 )

ax1.plot( energy, spec_p, "r", label="spec_p" )
ax1.plot( energy, pre, "green" )
ax1.set_xlabel( "Energy [eV]")

ax2.plot( energy, norm0, "r", label="norm0" )
ax2.plot( energy, post, "purple" )
ax2.set_xlabel( "Energy [eV]")

ax3.plot( energy, norm_p, "r", label="norm_p" )
ax3.set_xlabel( "Energy [eV]")

plt.tight_layout()
plt.show()

f:id:koideforest:20181010051755p:plain
うまく引けている感じがする。

規格化後のXASとXMCD。

xas = ( norm_p + norm_m ) / 2.
xmcd = norm_p - norm_m

f:id:koideforest:20181010053707p:plain

ピークの位置を、次の関数で求める。
最初のピークは710 ~ 730 eV、次のピークは730 ~ 750 eVの範囲にあるので、そのように設定。

def white_line_position( x, y, x0, x1 ):
    p = max([ yy for xx, yy in zip( x, y ) if x0 < xx < x1 ])
    for ix, xx in enumerate( x ):
        if x0 < xx < x1:
            if abs( y[ ix ] - p ) < 1e-7:
                break
    return ix, xx

_, e1 = white_line_position( energy, xas, 710, 730 )
_, e2 = white_line_position( energy, xas, 730, 750 )
print( e1, e2 )
# 720.36..., 740.54...

d状態以外の寄与を差し引くためのstep関数を次のように用意。

from scipy.special import erf

def step_func( x, x0, d ):
    return ( erf( ( x - x0 ) / d ) + 1. ) / 2.

def subtract_func( x, x01, x02, d, c ):
    return c * step_func( x, x01, d ) + ( c / 2. ) * step_func( x, x02, d)

元々の誤差関数は $\pm1$ の値を持つため、step_funcでは $[0,1$ ]になるように調整した。
強度比は軌道モーメントがなければ $L_3/L_2 = 2$ であるため、大きさを決める係数 $c$ は一つで良い。（この関数は、（遷移金属原子における）d状態以外の寄与のためのモデルであり、それらは一般に軌道モーメントは小さい。）

求めたピークの位置を使って、step関数を作る。
step関数のfitting領域は760 eV以降のpost edgeに設定。

def index_energy( energy, e0 ):
    for i, e in enumerate( energy ):
        if e >= e0:
            if abs( energy[ i ] - e0 ) <= abs( energy[ i - 1 ] - e0 ):
                return i
            else:
                return i - 1

i_post = index_energy( energy, 760 )

from scipy.optimize import curve_fit

c, _ = curve_fit(
    lambda e, c: subtract_func( e, e1, e2, de, c ),
    energy[ i_post : ],
    xas[ i_post : ]
    )
subtract = subtract_func( energy, e1, e2, de, c )

plt.plot( energy, subtract )
plt.plot( energy, xas )
plt.xlabel( "Energy [eV]" )
plt.ylabel( "Intensity [arb.unit]" )
plt.show()

f:id:koideforest:20181021055450p:plain
それらしく階段関数が置けている。
curve_fitの使い方は以下のサイトを参照。
pythonでフィッティングをする - おっぱいそん！

階段関数を差し引いたXAS、そしてXMCDを次のように積分する。
ここではSimpson法を用いる。

from scipy.integrate import simps

def integrate( x, y ):
    result = [ 0. ]
    for ix in range( len( x ) )[ 1: ]:
        result.append( simps( y[ :ix ], x[ :ix ] ) )
    return np.array( result )

int_xas = integrate( energy, xas - subtract )
int_xmcd = integrate( energy, xmcd )
plt.plot( energy, int_xas )
plt.plot( energy, int_xmcd )
plt.xlabel( "Energy [eV]" )
plt.show()

print( int_xas[-1] )
# 18.914...

f:id:koideforest:20181010060456p:plain
ちょっと後ろの方でフニャッとしているが、こんなもんだろう。

$L_3$ と $L_2$ の境目は、XASの一回微分がゼロになるところを持って来た。

def get_index( xs, x0 ):
    for i, x in enumerate( xs ):
        if x >= x0:
            if abs( xs[ i ] - x0 ) <= abs( xs[ i - 1 ] - x0 ):
                return i
            else:
                return i - 1

def index_zero( x, y, x_start, x_end ):
    iy_start = get_index( x, x_start )
    iy_end   = get_index( x, x_end )
    y_zero       = min( abs( y[ iy_start : iy_end ] ) )
    for iy, yy in enumerate( y ):
        if not iy_start <= iy <= iy_end: continue
        if abs( yy ) - y_zero <= 1e-7:
            return iy

i_zero = index_zero( energy[ 1 : ], np.diff( xas, n=1 ), 730, 735 )
plt.plot( energy, xas )
plt.plot( energy[ im ], xas[ im ], "o" )
plt.show()

print( int_xmcd[ i_zero ] )
# -4.5744...

f:id:koideforest:20181021064816p:plain
差分による近似的に微分した配列が得られるnp.diffの使い方は、以下のサイトを参照。
要素の差分、足し合わせを計算するNumPyのdiff関数とcumsum関数の使い方 - DeepAge

ひとまずはこれでSum ruleに必要なXAS及びXMCDの積分値を求めることが出来た。

ちなみに、適当に作ったスペクトルは以下のスクリプトを使用。

def lorentzian_norm( x, x0, gamma ):
    return gamma**2 / ( ( x - x0 )**2 + gamma**2 ) 

energy = np.linspace( 680, 800, 1000 )
de = energy[1] - energy[0]

peak1 = 720
peak2 = 740

back  = ( 0.01 * ( energy - 900 ) )**2
base  = subtract_func( energy, peak1, peak2, 10 * de, 1 )

peak_p  = 2 * lorentzian_norm( energy, peak1, 2 )
peak_p += 1 * lorentzian_norm( energy, peak2, 2 )
spec_p = base + peak_p + back

peak_m  = 3.0 * lorentzian_norm( energy, peak1, 2 )
peak_m += 0.5 * lorentzian_norm( energy, peak2, 2 )
spec_m = base + peak_m + back

ローレンツ関数の大きさ等は以下を参照。
ローレンツ関数とは？

2018-10-06

自由電子Green関数の部分波展開

数学量子力学

※訂正: マイナスの境界条件で内向き波になるように訂正(2019/03/08)

Hartree原子単位系

$\displaystyle k = \sqrt{ 2 E } \\ \displaystyle J_L( k, {\bf r } ) = i^l \, j_l( k r ) \, Y_L( \hat{\bf r }) \\ \displaystyle H^{(\pm)}_L( k, {\bf r } ) = i^l \, h^{(\pm)}_l( k r ) \, Y_L( \hat{\bf r }) \\ \displaystyle G^{\pm}( k, {\bf r } - {\bf r }' ) = \mp 2 i k \, \sum_L \, J^*_L( k, {\bf r }_< ) \, H^{(\pm)}_L( k, {\bf r }_> )$

以下、これを証明する。

$\displaystyle H_0 = - \frac{ \nabla^2 }{ 2 } \\ \displaystyle \left( \frac{ k^2 }{ 2 } - H_0 \right) G( k, {\bf r } - {\bf r }' ) = \delta( {\bf r } - {\bf r }' )$

フーリエ変換を使って、微分方程式を解く。
$\displaystyle G( k, {\bf r } - {\bf r }' ) = \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } \, G( k, {\bf q } ) \, e^{ i { \bf q } \cdot ( {\bf r } - {\bf r }' ) } \\ \displaystyle \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } \left( \frac{ k^2 }{ 2 } - \frac{ q^2 }{ 2 } \right) G( k, {\bf q } ) e^{ i { \bf q } \cdot ( {\bf r } - {\bf r }' ) } = \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } e^{ i { \bf q } \cdot ( {\bf r } - {\bf r }' ) } \\ \displaystyle \left( \frac{ k^2 }{ 2 } - \frac{ q^2 }{ 2 } \right) G( k, {\bf q } ) = 1 \\ \displaystyle G^{ \pm }( k, {\bf q } ) = 2 \, \frac{ 1 }{ k^2 - q^2 \pm i \eta } \quad ( \eta \rightarrow 0^+ ) \\ \displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = 2 \, \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } \, \frac{ e^{ i { \bf q } \cdot ( {\bf r } - {\bf r }' ) } }{ k^2 - q^2 \pm i \eta } \, = 2 \, \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } \, \frac{ e^{ i { \bf q } \cdot {\bf r } } e^{ - i { \bf q } \cdot {\bf r }' } }{ k^2 - q^2 \pm i \eta } \,$

平面波の角運動量展開
$\displaystyle e^{ i {\bf k} \cdot {\bf r} } = 4 \pi \sum_L J_L( k, {\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 J_L( k, {\bf r} ) Y^*_L( \hat{\bf k} ) \right)$

$\displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = 2 \, \int \frac{ d {\bf q} }{ ( 2 \pi )^3 } \, \frac{ ( 4 \pi )^2 }{ k^2 - q^2 \pm i \eta } \, \sum_L \sum_{L'} J_L( q, {\bf r} ) Y^*_L( \hat{\bf q} ) J^*_{L'}( q, {\bf r}' ) Y_{L'}( \hat{\bf q} ) \\ \displaystyle = 2 \, \int^{\infty}_0 \frac{ q^2 \, dq }{ ( 2 \pi )^3 } \, \frac{ ( 4 \pi )^2 }{ k^2 - q^2 \pm i \eta } \, \sum_L \sum_{L'} J_L( q, {\bf r} ) J^*_{L'}( q, {\bf r}' ) \int d \hat{\bf q} Y^*_L( \hat{\bf q} ) Y_{L'}( \hat{\bf q} ) \\ \displaystyle = 2 \, \int^{\infty}_0 \frac{ q^2 \, dq }{ ( 2 \pi )^3 } \, \frac{ ( 4 \pi )^2 }{ k^2 - q^2 \pm i \eta } \, \sum_L J_L( q, {\bf r} ) J^*_{L}( q, {\bf r}' )$

パリティ対称性 $j_l( - \rho ) = (-1)^l j_l( \rho )$ より、被積分関数は $q$ に対して偶関数。よって、
$\displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = \int^{\infty}_{-\infty} \frac{ q^2 \, dq }{ ( 2 \pi )^3 } \, \frac{ ( 4 \pi )^2 }{ k^2 - q^2 \pm i \eta } \, \sum_L J_L( q, {\bf r} ) J^*_{L}( q, {\bf r}' )$

球ベッセル関数は特異点を持たないので、分母にのみ特異性がある。
これを利用して、留数定理を用いて積分を計算する（そのために積分区間を広げた）。
分母の $q$ に対する極は、
$\displaystyle k^2 - q^2 \pm i \eta = - ( q^2 - ( k^2 \pm i \eta ) )= - \left( q^2 - \sqrt{ k^4 + \eta^2 }e^{ \pm i \theta } \right) = - \left( q^2 - k^2 e^{ \pm i \theta } \right) \\ \displaystyle = - \left( q + k e^{ \pm i \theta / 2 } \right) \left( q - k e^{ \pm i \theta / 2 } \right) \\ \displaystyle = - \left( q + ( k \pm i \eta ) \right) \left( q - ( k \pm i \eta ) \right) \\ \displaystyle = - \left( q - ( - k \mp i \eta ) \right) \left( q - ( k \pm i \eta ) \right) \\ \displaystyle \left( \theta = \tan^{-1}( \eta / k^2 ) > 0 \right)$

ややこしいが、結局は以下の通りである。

$G^{+}$ : $q^{\uparrow} = k$ ; $q^{\downarrow} = -k$
$G^{-}$ : $q^{\uparrow} = -k$ ; $q^{\downarrow} = k$

Jordanの補助定理を使うことを見越して、無限遠 $\rho \rightarrow \infty$ における各関数の漸近系をまとめておくと、
$\displaystyle j_l( \rho ) = \frac{ 1 }{ 2 } \left( h^{(1)}_l( \rho ) + h^{(2)}_l( \rho ) \right) \\ \displaystyle h^{(\pm)}_l( \rho ) \rightarrow \pm \frac{ 1 }{ i \rho } e^{ \pm i ( \rho - l \pi /2 ) }$

以下、多少恣意的にバラす。
$\displaystyle j_l( q r ) j_l( q r' ) = \frac{ 1 }{ 2 } \, j_l( q r ) \left( h^{(+)}_l( q r' ) + h^{(-)}_l( q r' ) \right) \\ \rightarrow \frac{ 1 }{ 4 i q^2 r r' } \left( \left( e^{ i q ( r + r' ) } e^{ - i l \pi } - e^{ i q ( - r + r' ) } \right) + \left( - e^{ i q ( r - r' ) } + e^{ i q ( - r - r' ) } e^{ - i l \pi } \right) \right)$

指数関数の肩に着目し、 $r < r'$ を仮定すると、

$j_l( q r ) h^{(+)}_l( q r' )$ ：上半面で指数関数的に減少。
$j_l( q r ) h^{(-)}_l( q r' )$ ：下半面で指数関数的に減少。

したがって、上半面と下半面に分けて積分を行う。
下半面で積分路を取る時に、左回りにするために $\infty \rightarrow -\infty$ と変更する必要があることから、マイナスが出てくることに注意。

$\displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = \frac{1}{2} \frac{ ( 4\pi )^2 }{ ( 2 \pi )^3 } \sum_L Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf r}' ) \, \int^{\infty}_{-\infty} \, \frac{ dq \, q^2 \, j_l( q r ) \left( h^{(+)}_l( q r' ) + h^{(-)}_l( q r' ) \right) }{ - ( q - ( - k \mp i \eta ) ) ( q - ( k \pm i \eta ) ) } \\ \displaystyle = \frac{1}{2} \frac{ ( 4\pi )^2 }{ ( 2 \pi )^3 } \sum_L Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf r}' ) \, \left( \oint_{\uparrow} \, \frac{ dq \, q^2 \, j_l( q r ) h^{(+)}_l( q r' ) }{ - ( q - ( - k \mp i \eta ) ) ( q - ( k \pm i \eta ) ) } - \oint_{\downarrow} \, \frac{ dq \, q^2 \, j_l( q r ) h^{(-)}_l( q r' ) }{ - ( q - ( - k \mp i \eta ) ) ( q - ( k \pm i \eta ) ) } \right) \\ \displaystyle = \frac{1}{2} \frac{ ( 4\pi )^2 }{ ( 2 \pi )^3 } \sum_L Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf r}' ) \, 2 \pi i \left( \frac{ (\pm k)^2 \, j_l( \pm k r ) h^{(+)}_l( \pm k r' ) }{ \mp 2 k } - \frac{ ( \mp k )^2 \, j_l( \mp k r ) h^{(-)}_l( \mp k r' ) }{ \pm 2 k } \right)$

$h^{(-)}_l( - \rho ) = (-1)^l h^{(+)}_l( \rho )$ より、

$\displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = \frac{1}{2} \frac{ ( 4\pi )^2 }{ ( 2 \pi )^3 } \sum_L Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf r}' ) \, 2 \pi i \left( \frac{ k^2 \, j_l( k r ) h^{(\pm)}_l( k r' ) }{ \mp 2 k } - \frac{ k^2 \, j_l( k r ) h^{(\pm)}_l( k r' ) }{ \pm 2 k } \right) \\ \displaystyle = \frac{1}{2} \frac{ ( 4\pi )^2 }{ ( 2 \pi )^3 } \sum_L Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf r}' ) \, 2 \pi i \left( \mp k \, j_l( k r ) h^{(\pm)}_l( k r' ) \right) \\ \displaystyle = \mp 2 i k \sum_L j_l( k r ) Y_{L}( \hat{\bf r} ) h^{(\pm)}_l( k r' ) Y^*_{L}( \hat{\bf r}' ) \\ \displaystyle = \mp 2 i k \sum_L j_l( k r ) Y^*_{\bar{L}}( \hat{\bf r} ) h^{(\pm)}_l( k r' ) Y_{ \bar{L}}( \hat{\bf r}' ) \\ \displaystyle = \mp 2 i k \sum_{\bar{L}} j_l( k r ) Y^*_{\bar{L}}( \hat{\bf r} ) h^{(\pm)}_l( k r' ) Y_{ \bar{L}}( \hat{\bf r}' ) \\ \displaystyle = \mp 2 i k \sum_{L} j_l( k r ) Y^*_{L}( \hat{\bf r} ) h^{(\pm)}_l( k r' ) Y_{ L }( \hat{\bf r}' ) \\ \displaystyle = \mp 2 i k \sum_L J^*_L( k, {\bf r} ) H^{(\pm)}_L( k, {\bf r}' )$
$L$ の和を取っていることと、 $m$ に依存するのが $Y_L$ しかいないため、 $L \rightarrow \bar{L}$ の置き換えが可能。

$r > r'$ の場合は、 $j_l( q r ) j_l( q r' )$ より、明らかに $r \leftrightarrow r'$ の交換が可能であるから、 $r < r'$ の場合と同じ結果が得られる。

したがって、

$\displaystyle G^{ \pm }( k, {\bf r } - {\bf r }' ) = \mp 2 i k \sum_L J^*_L( k, {\bf r}_< ) H^{(\pm)}_L( k, {\bf r}_> )$

2018-10-04

MacOSでpythonをインストールし直す

プログラミング

python2を主に使っていたが、仕事上ちゃんとpython3の環境を整えなければならなくなり、その際の記録。

いろいろお試し感覚でごちゃごちゃ入れていて、

System由来のpython2
手動経由のpython2
MacPorts経由のpython2
pyenv経由のmicroconda (python3)

が混在していた。
メインでは手動で入れたものを使用していた。それもあり、一度整理しておきたいと思っていた。

手動で入れたpython2は /Library/Python/... に入っていたので、Pythonのフォルダごと消去した。

sudo rm -r /Library/Python

.../bin/の中のpip等も手動で消した。

MacPortsでは以下のコマンドでアンインストールが行える。

port installed
sudo port uninstall python

しかし、python_select などの細かいものは一気に消してはくれない。残しておくのはそれはそれで気持ち悪いので、一個一個消した。
その他、MacPortsが主に散らかす /opt/local/ の各フォルダの中にも（残しておいても問題ないが）ゴミが残るので、完全に綺麗に消すのは難しいと思う。

microcondaも、pipとcondaは混ぜない方が良いという記事を発見したので、pyenvのコマンドでアンインストールした。
condaとpip：混ぜるな危険 - onoz000’s blog

何も考えずにpyenvを使うことは推奨されていないが、簡単にバージョンごとアンインストール等が出来るのは、リソースの限られたノートパソコン上では悪くないと感じた。
pyenvが必要かどうかフローチャート
なので、pyenvはこのまま採用。

ただし、何故だか知らないが、（少なくともMacOSでは）pyenvにはupdateのコマンドが内蔵されていない。
幸い、git経由でpyenvをインストールしていたため、git pull で更新することが出来た。

cd ~/.pyenv
git pull

これによって、python3.7.0がインストールリストに表示され、インストールすることが出来た。
ただし、MacOSでは、matplotlib等を使う場合、pythonはFramework上でインストールしなければエラーが出るため、（仕組みは知らないが）以下のコマンドでインストールをした。

pyenv install -l
env PYTHON_CONFIGURE_OPTS="--enable-framework" pyenv install 3.7.0
pyenv local 3.7.0
pyenv global 3.7.0
pyenv rehash

Home · pyenv/pyenv Wiki · GitHub

pyenvでpython3.7.0をインストールした後は、pipでモジュールをインストールすることになるが、（pyenvのフォルダがhome直下にあるので必要ないが）sudoを付けるとpipがいないと怒られた。
これは環境変数がsudoだとリセットされることで起こるらしいので、自分はリセットしないように変更した。

sudo visudo

sudo でコマンド打ったら「コマンドが見つかりません」と言われたときの気持ち - 彼女からは、おいちゃんと呼ばれています

jupiter notebookは以下のサイトを参考にインストール
Jupyter Notebookのインストール
他のpythonのバージョンを加える場合には以下のサイトが参考になった。
Jupyter Notebookでpython3 Kernelを追加するのにはまったメモ

ちなみに、Juliaをjupiter notebookのカーネルに追加するには、以下のサイトを参照。
Julia v0.6.1のインストールとJupyter Notebookで使うまで

2018-09-21

横軸と縦軸のデータをすぐに保存&読込 (python)

プログラミング

pythonで作業していて、パッとデータを保存する、もしくは読み込む方法。
結局は、横軸と縦軸の各一次元データ（横ベクトル）を結合して縦長のデータにする部分が面倒くさい。
しかし、それは numpy.c_ で解決することがわかった。
Pythonで配列や行列の結合

import numpy as np

# save data
datas1 = np.c_[ x1, y1 ]  # same as np.vstack([ x1, y1 ]).T
np.savetxt( "datas1.dat", datas )

# load data
datas2 = np.loadtxt( "datas2.dat" )
x2 = datas2[ :, 0 ]  # Caution: not same as datas2[ : ][ 0 ]
y2 = datas2[ :, 1 ]

nano_exit

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

3価のTiイオンのLII, LIII端における原子多重項遷移

Green関数の固有関数展開

自由電子Green関数のFourier変換

XMCDのSum ruleをpythonで処理

自由電子Green関数の部分波展開

MacOSでpythonをインストールし直す

横軸と縦軸のデータをすぐに保存&読込 (python)