以前、線形応答でバンドから応答係数を算出するということをやってみました。
- 理論パート：https://qiita.com/hiro949/items/bdb5516fbd9a3b8806b6
- 計算パート：https://qiita.com/hiro949/items/ee7ae41d1aca8903f4ba
前回は相互作用を考えていませんでしたが、今回はそこに相互作用の効果も入れてみたいと思います。
今回も物理学専攻で学部３年くらいまでの量子力学の知識を前提としてお話ししたいと思います。
記号などは線形応答計算の理論パートに従うことにします。また、煩雑なので計算の詳細は省きます。
イジング模型を学ばれた方は低温に行くとシステムが転移温度$T_c$を境に相転移して磁石になるという藻をやったことがあると思います。そのとき、$T_c$近傍で帯磁率$\chi$が発散する現象があったかと思います。今回の計算でも同様に相互作用を加味して帯磁率、誘電率などを計算すると、対応する量に疎転移があると転移店で発散するという現象を見ることができます。
ハバード模型
タイトバインディングモデルは相互作用を考えていませんでしたが、今回は同じ原子にいる電子同士にはクーロン力が働くとします。この分の効果を$H_{int}$とします。
\begin{align}
H &= H_0 + H_{int} \\
H_0 &= \sum_{<\mathbf{R},\mathbf{R'}>}\sum_{\alpha,\beta} t_{\mathbf{R},\alpha \leftarrow \mathbf{R'}, \beta} c^\dagger_{\mathbf{R},\alpha} c_{\mathbf{R'}, \beta} \\
H_{int} &= \frac{!}{2} \sum_\mathbf{R} \sum_{\alpha_1\cdots\alpha_4} U_{\alpha_1\cdots\alpha_4} c^\dagger_{\mathbf{R},\alpha_1} c^\dagger_{\mathbf{R},\alpha_2} c_{\mathbf{R},\alpha_3} c_{\mathbf{R},\alpha_4}
&= \frac{1}{2}\sum_\mathbf{R} \sum_{\alpha_1\cdots\alpha_4} U_{\alpha_1\cdots\alpha_4} c^\dagger_{\mathbf{R},\alpha_1} c^\dagger_{\mathbf{R},\alpha_2} c_{\mathbf{R},\alpha_3} c_{\mathbf{R},\alpha_4} \\
& = \frac{1}{2L^2} \sum_\mathbf{k_1,\cdots,k_4} \sum_{\alpha_1\cdots\alpha_4} U_{\alpha_1\cdots\alpha_4} c^\dagger_{\mathbf{k}_1,\alpha_1} c^\dagger_{\mathbf{k}_2,\alpha_2} c_{\mathbf{k}_3,\alpha_3} c_{\mathbf{k_4},\alpha_4}
e^{-i\mathbf{R} \cdot(\mathbf{k_1 + k_2-k_3-k_4})} e^{-i\mathbf{r_1} \cdot \mathbf{k_1} -i\mathbf{r_2} \cdot \mathbf{k_2} + i\mathbf{r_2} \cdot \mathbf{k_3} + i\mathbf{r_1} \cdot \mathbf{k_4} } \\
& = \frac{1}{2L} \sum_\mathbf{k_1,k_2,k_3} \sum_{\alpha_1\cdots\alpha_4} U_{\alpha_1\cdots\alpha_4} c^\dagger_{\mathbf{k}_1,\alpha_1} c^\dagger_{\mathbf{k}_2,\alpha_2} c_{\mathbf{k}_3,\alpha_3} c_{\mathbf{k_1 + k_2-k_3},\alpha_4}
e^{i(\mathbf{r_1} - \mathbf{r_2} ) \cdot (\mathbf{k_2} - \mathbf{k_3}) } \\
& = \frac{1}{2L} \sum_\mathbf{k,k'} \sum_{\alpha_1\cdots\alpha_4} \sum_\mathbf{q} U_{\alpha_1\cdots\alpha_4} e^{i(\mathbf{r_1} - \mathbf{r_2} ) \cdot \mathbf{q} } c^\dagger_{\mathbf{k},\alpha_1} c^\dagger_{\mathbf{k'},\alpha_2} c_{\mathbf{k'-q},\alpha_3} c_{\mathbf{k + q},\alpha_4}
\end{align}
もっと一般的には$U$にも$\mathbf{q}$依存性があって
H_{int} = \frac{1}{2L} \sum_\mathbf{k,k'} \sum_{\alpha_1\cdots\alpha_4} \sum_\mathbf{q} U_{\alpha_1\cdots\alpha_4}(\mathbf{q}) c^\dagger_{\mathbf{k},\alpha_1} c^\dagger_{\mathbf{k'},\alpha_2} c_{\mathbf{k'-q},\alpha_3} c_{\mathbf{k + q},\alpha_4}
と書けます。
各原子の波動関数を$\psi_\alpha$として、相互作用の大きさは
U_{\alpha_1\cdots\alpha_4} = \int\int d\mathbf{r} d\mathbf{r'} \psi^*_{\alpha_1}(\mathbf{r+R}) \psi^*_{\alpha_2}(\mathbf{r'+R})\frac{e^2}{|\mathbf{r-r'}|} \psi^*_{\alpha_3}(\mathbf{r'+R}) \psi^*_{\alpha_4}(\mathbf{r+R})
のように計算されますが、s,p,d,f,...軌道で波動関数の形状は違い、パウリの排他率で電子スピン依存性もあります。したがって相互作用の大きさは軌道やスピンのパラメーターに依存するテンソルになります。
簡単な場合
一番簡単なのはモデルで考える軌道がs軌道のみの場合です。このときはスピンの自由度$\sigma$だけが関係していて、定数$U$を用いて
H_{int} = \sum_\mathbf{R} \sum_{\sigma\sigma'} \frac{U}{2} c^\dagger_{\mathbf{R},\sigma} c^\dagger_{\mathbf{R},\sigma'} c_{\mathbf{R},\sigma'} c_{\mathbf{R},\sigma}
となります。
摂動計算
応答係数
\chi_{BA}(\mathbf{q},i\omega_l) = \int_0^\beta d\tau <[B_i(\mathbf{q},\tau),A_j^\dagger(-\mathbf{q},0)]>
を相互作用表示$O(\tau)_I = e^{\tau H_0 } O e^{-\tau H_0 }$で表すと
\begin{align}
\chi_{BA}(\mathbf{q},i\omega_l) &= \sum_{i_1\cdots i_4} B_{i_1i_2}A_{i_3i_4}\frac{1}{L^3}\sum_\mathbf{k} \int_0^\beta d\tau e^{i\omega_l \tau} \sum_{n=0}^\infty \frac{(-1)^n}{n!}
\int_0^\beta d\tau_1 \cdots \int_0^\beta d\tau_n <T H_{int}(\tau_1)_I \cdots H_{int}(\tau_n)_I c^\dagger_{\alpha_1}(\mathbf{k};\tau)_I c_{\alpha_2}(\mathbf{k+q};\tau)_I c^\dagger_{\alpha_4}(\mathbf{k+q}) c_{\alpha_3}(\mathbf{k+q})> \\
&=\sum_{i_1\cdots i_4} B_{i_1i_2}A_{i_4i_3} \sum_{n=0}^\infty \Phi^{(n)}_{i_1 \cdots i_4}(\mathbf{q},i\omega_l)
\end{align}
となります。$n=0$の項は線形応答理論のときと同じです。次に$n=1$の項を計算します。
ウィック展開や松原周波数でのフーリエ級数展開を用いて計算を進めると
\begin{align}
\Phi^{(1)}_{i_1 \cdots i_4}(\mathbf{q},i\omega_l) &= -\frac{1}{\beta^2 L^6} \sum_{n,m}\sum_\mathbf{k,k'}\sum_{j_1\cdots j_4}(U_{i_1 i_2 i_3 i_4} + U_{i_2 i_1 i_4 i_3}) [
G_{j_4 i_1}(\mathbf{k},i\omega_n)G_{i_2 j_1}(\mathbf{k+q},i\omega_n+i\omega_l)G_{i_3 j_2}(\mathbf{k'},i\omega_m)G_{j_3 i_4}(\mathbf{k'+q},i\omega_m+i\omega_l) - G_{j_3 i_1}(\mathbf{k},i\omega_n)G_{i_2 j_1}(\mathbf{k+q},i\omega_n+i\omega_l)G_{i_3 j_2}(\mathbf{k'},i\omega_m)G_{j_4 i_4}(\mathbf{k'+q},i\omega_m+i\omega_l)
] \\
&= -\sum_{j_1\cdots j_4}( U_{i_1 i_2 i_3 i_4} + U_{i_2 i_1 i_4 i_3} - U_{i_1 i_2 i_4 i_3} - U_{i_2 i_1 i_3 i_4}) \Phi^{(0)}_{i_1i_2 j_4 j_1}(\mathbf{q},i\omega_l)\Phi^{(0)}_{j_2j_4 i_3 i_4}(\mathbf{q},i\omega_l)
\end{align}
ここで、
V_{j_1j_2j_3j_4} = U_{j_2j_3j_4j_1} + U_{j_3j_2j_1j_4} - U_{j_2j_3j_1j_4} - U_{j_3j_2j_4j_1}
とすると
\Phi^{(1)}_{i_1 \cdots i_4}(\mathbf{q},i\omega_l) = -\sum_{j_1\cdots j_4} \Phi^{(0)}_{i_1i_2 j_1 j_2}(\mathbf{q},i\omega_l) V_{j_1j_2j_3j_4} \Phi^{(0)}_{j_3j_4 i_3 i_4}(\mathbf{q},i\omega_l)
となります。$\Phi^{(1)}$のうち、$V$の最初の２つを含む項がハートリー項、後ろの２項がフォック項と呼ばれています。
$(j_1,j_2)$、$(j_3,j_4)$をそれぞれ組にして辞書順に並べて番号付けすると、$ V_{j_1j_2j_3j_4}$は行列表示できます。これは$\Phi^{(0)}$も同様です。行列として書くと、
\Phi^{(1)}(\mathbf{q},i\omega_l) = -\Phi^{(0)}(\mathbf{q},i\omega_l)V \Phi^{(0)}(\mathbf{q},i\omega_l)
と書けます。
乱雑位相近似(RPA)
まとめると、１次元のオーダーで
\Phi(\mathbf{q},i\omega_l) = \Phi^{(0)}(\mathbf{q},i\omega_l) - \Phi^{(0)}(\mathbf{q},i\omega_l)V \Phi^{(0)}(\mathbf{q},i\omega_l)
となります。これを
\Phi^{RPA}(\mathbf{q},i\omega_l) = \Phi^{(0)}(\mathbf{q},i\omega_l) - \Phi^{(0)}(\mathbf{q},i\omega_l)V \Phi^{RPA}(\mathbf{q},i\omega_l)
と書き換えれば(これは式変形ではない)、
\begin{align}
\Phi^{RPA}(\mathbf{q},i\omega_l) &= [ \mathbf{1} + \Phi^{(0)}(\mathbf{q},i\omega_l)V ]^{-1} \Phi^{(0)}(\mathbf{q},i\omega_l) \\
\chi_{BA}(q,\omega+i\delta) &= \sum_{i_1\cdots i_4}B_{i_1i_2}\Phi^{RPA}_{i_1\cdots i_4}(\mathbf{q},\omega+i\delta)A^*_{i_3i_4}
\end{align}
のようになり、これを応答関数の乱雑位相近似(RPA)といいます(ここまで長かった...)。
完全ではありませんがRPAには高次のオーダーの量も含まれていて、磁気的な相転移に対しては割とうまくいくことが多いです。
正方格子模型の例
https://qiita.com/hiro949/items/ee7ae41d1aca8903f4ba
でやった線形応答計算をRPAにしてやってみましょう。
$(\uparrow\uparrow),(\uparrow\downarrow),(\downarrow\uparrow),(\downarrow\downarrow)$の順に並べると
\begin{align}
U &= \frac{1}{2}\begin{pmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & U & 0 \\
0 & U & 0 & 0 \\
0 & 0 & 0 & 0
\end{pmatrix} \\
V &= \begin{pmatrix}
0 & 0 & 0 & U \\
0 & -U & 0 & 0 \\
0 & 0 & -U & 0 \\
U & 0 & 0 & 0
\end{pmatrix}
\end{align}
となります。
また０次の項はユニリタ行列が$U(\mathbf{k})=I$だから対角成分しかなく、また$E_\uparrow(k)=E_\downarrow(k)$なので
\Phi^{(0)} = \Phi_0\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
です。よってRPA相関関数は
\begin{align}
\Phi &= (1 + \Phi^{(0)}V)^{-1} \Phi_0 \\
&= \begin{pmatrix}
1 & 0 & 0 & U\Phi_0 \\
0 & 1-U\Phi_0 & 0 & 0 \\
0 & 0 & 1-U\Phi_0 & 0 \\
U\Phi_0 & 0 & 0 & 1
\end{pmatrix}^{-1} \Phi_0 \\
&= \begin{pmatrix}
[1-U^2\Phi_0^2]^{-1} & 0 & 0 & -U\Phi_0[1-U^2\Phi_0^2]^{-1} \\
0 & [1-U\Phi_0]^{-1} & 0 & 0 \\
0 & 0 & [1-U\Phi_0]^{-1} & 0 \\
-U\Phi_0[1-U\Phi_0^2]^{-1} & 0 & 0 & [1-U^2\Phi_0^2]^{-1}
\end{pmatrix} \Phi_0
\end{align}
です。このとき物理量$A,B$はベクトルとして表示されて、
\begin{align}
n &= \begin{pmatrix} 1 & 0 & 0 & 1 \end{pmatrix} \\
\sigma_x &= \begin{pmatrix} 0 & 1 & 1 & 0 \end{pmatrix} \\
\sigma_y &= \begin{pmatrix} 0 & -i & i & 0 \end{pmatrix} \\
\sigma_z &= \begin{pmatrix} 1 & 0 & 0 & -1 \end{pmatrix}
\end{align}
だから(電荷密度演算子は単位行列、磁化演算子はパウリ行列を使う)、
RPA感受率を計算すると
誘電率が
\chi^c = \frac{2\chi_0}{1+U\chi_0}
帯磁率がx,y,z成分すべて
\chi^s = \frac{2\chi_0}{1-U\chi_0}
となります（$\chi_0=\Phi_0$）。
物理的考察
ちなみに、これら４つの物理量は$\Phi$の固有ベクトルになっていて、$\chi^s,\chi^c$はこれらの量が属する固有値になっています($\chi^s$については三重縮退しています。)。また、$\chi^s(\pi,\pi) > \chi^c(\pi,\pi)$であることから
\begin{align}
S^x(n,m) &= \sigma_xe^{i(n+m)\pi} = (-1)^{n+m} \sigma_x \\
S^y(n,m) &= \sigma_ye^{i(n+m)\pi} = (-1)^{n+m} \sigma_y \\
S^z(n,m) &= \sigma_ze^{i(n+m)\pi} = (-1)^{n+m} \sigma_z
\end{align}
の秩序が発達しやすいと分かります。ここで$(n,m)$はサイトの位置の番号です。つまり、反強磁性が発達すると予測されます(実際このハバード模型では反強磁性相転移が起きます)。
実装
以上の計算からRPA感受率のコードは線形応答とほとんど同じで、最後に$\chi_0/(1 \pm U\chi_0)$の処理を追加するだけです。
ちなみに、ここではこの置き換えで横着をしましたが、ちゃんと$V$をつくって今の手計算と同様の計算過程を踏むコードの方が汎用性があります。
でそのようなコードを実装し、もっと一般的な場合を取り扱いました。
rpa.py
import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
t = 1.
U = 4.
system_size = 32
dimension = 2
volume = system_size**dimension
temp_limit = 4.
temp_iteration = 128
temp_widtg = temp_limit/temp_iteration
temp_offset = 1e-4
def calc_boltzmann_factor(energy, kbt):
inv_temp = 1/(kbt)
a = -inv_temp*(energy)
if a > 709: #prevent overflow
a = 709
elif a < -745: #prevent underflow
a = -745
return np.exp( a )
def fermi_distribution(kbt,energy,chemical_potential):
return 1/(calc_boltzmann_factor(chemical_potential-energy, kbt)+1)
def extend_k_mesh(k_mesh_list, ki_mesh_list):
return [ k + ki for k in k_mesh_list for ki in ki_mesh_list ]
def make_k_mesh_list():
global system_size, dimension
k0 = 2*np.pi/system_size
k_mesh_list = [ [ n*k0 ] for n in range(system_size) ]
ki_mesh_list = copy.deepcopy( k_mesh_list )
for i in range( dimension-1 ):
k_mesh_list = extend_k_mesh(k_mesh_list, ki_mesh_list)
return k_mesh_list
def make_temperature_list():
global temp_limit, temp_iteration, temp_width, temp_offset
return [ n*temp_width + temp_offset for n in range(temp_iteration+1) ]
#(1)
def culc_green_function_core( E_kq, E_k, omega, kbt, chemical_potential):
global volume
delta = 0.1/volume
numerator = fermi_distribution(kbt,E_kq,chemical_potential) - fermi_distribution(kbt,E_k,chemical_potential)
denominator = ( E_kq - E_k + omega - 1j*delta )/ ( ( E_kq - E_k + omega )**2 + delta**2 )
return -numerator*denominator
# eigenenergy
def hop(k):
global t
return -2*t*( np.cos(k[0]) + np.cos(k[1]) )
#(2)
def susceptibility_loc_k(kbt,omega,q,k):
E_k= hop(k)
E_kq = hop(k+q)
return culc_green_function_core( E_kq, E_k, omega, kbt, 0)
#(3)
def susceptibility( susceptibility_index, q, omega, kbt, chemical_potential, k ):
chi0 = (1/volume)*np.sum( np.array( [ test_susceptibility_loc_k(kbt,omega,q,k) for k in make_k_mesh_list() ] ) )
#### 変更した箇所 ####
return (2*chi0)/( 1 - U*chi0 )
def calc_kbt_dependence_of_susceptibilities():
q = np.array( [np.pi, np.pi] )
omega = 0
return [ susceptibility(kbt,omega,q) for kbt in make_temperature_list() ]
temperature = make_temperature_list()
chi1 = calc_kbt_dependence_of_susceptibilities()
plt.style.context('classic')
plt.plot(temperature, chi1)
plt.xlim([0,temp_limit])
plt.xlabel("temperature")
plt.ylabel("susceptibility")
plt.ylim([-0.001,2.5])
plt.show()
結果は次のようになります。
図から確かに帯磁率の発散がみられます。
注意
- RPAは相転移する前の高温無秩序状態でしか意味のない量なので、転移点より低温は無視して大丈夫です。
- 電気伝導度などには使えません。$J$が$\mathbf{k}$に依存するため、うまく式を分解できないからです。