線形応答理論でバンドから物性値を計算(理論パート)

エレクトロニクス、材料科学や物性物理学では電気伝導度(電気抵抗の逆数)や誘電率、帯磁率といった応答係数を計算する方法がいくつかあり、よく知られた手法の一つが線形応答理論です。これを使うと、金属における物性値をバンド図から計算することができるようになります。材料の電気特性を調べるのによく使う基本的な方法なのですが、僕は理解するまでに苦労したのでここで紹介させていただきます。
記事が長くなってしまうので、ここではバンドとは何か分かる程度の量子力学の知識を前提とさせていただきます(おそらく物理学専攻で大学3年を終わった時点くらいの知識量です)。ここではコードに落としやすいように、少々複雑ですが一般の場合を考えたいと思います。まずはどういう計算をするのかを説明します。

実際のシミュレーションは以下のページでやります：
https://qiita.com/hiro949/items/ee7ae41d1aca8903f4ba

1. タイトバインディングモデル

一般に金属を表す結晶格子のタイトバインディングモデルのハミルトニアンは一般に以下のように書けます：

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

ここで各変数を説明すると、$\mathbf{R},\mathbf{R'}$は格子上にある各原子の座標で、$\alpha, \beta$はそれ以外のパラメーターを指します：例えば、電子スピンの$\uparrow, \downarrow$や$s,p,d,f,...$軌道のうちどの原子軌道にいる電子なのか、ブラべ格子に対する副格子自由度などです。
$c^\dagger,c$は各状態${\mathbf{r},\alpha}$の電子の生成消滅演算子です。
$t_{\mathbf{R},\alpha \leftarrow \mathbf{R'}, \beta}$は重なり積分値で、状態$\mathbf{R'}, \beta$から$\mathbf{R},\alpha$へのホッピングパラメーターになっています。金属結晶では並進対称性があるので、$t$は$\Delta \mathbf{R} = \mathbf{R} - \mathbf{R'}$についてのみ座標依存性をもちます。
$<\mathbf{R},\mathbf{R'}>$は原子が結合している場合に和をとることを表しています。結合していないときに$t=0$とおくことで代用可能です。
原子の位置は基本並進ベクトル$\mathbf{a}_i$と副格子自由度の分の座標移動$ (\mathbf{r}_a)$ $(a=1\cdots m)$を用いて、

\mathbf{R} = \mathbf{R}_a( \mathbf{n}) = \sum_{i=1,2,3}(  n_i \mathbf{a}_i + \mathbf{r}_a )

と表せます。
以上のことを考慮して、周期的境界条件$n_i \mod L$のもとで逆格子ベクトル$ \mathbf{b}_i　$を用いたフーリエ変換

\mathbf{k} = \sum_{i} q_i \mathbf{b}_i \\
c_{\mathbf{r}_a, \alpha}(\mathbf{k}) = \frac{1}{\sqrt{L}}\sum_{\mathbf{n}} c_{\mathbf{R}_a( \mathbf{n}), \alpha} e^{-i \mathbf{k}\cdot\mathbf{R}_a( \mathbf{n})}

を行うと、

H = \sum_{\mathbf{k}}\sum_{<\mathbf{R}_a(\Delta \mathbf{n}), \mathbf{r}_b>}\sum_{\alpha,\beta} t_{\mathbf{R}_a(\Delta \mathbf{n}),\alpha \leftarrow 0, \mathbf{r}_b \beta} c^\dagger_{\mathbf{r}_a,\alpha}(\mathbf{k}) c_{\mathbf{r}_b, \beta}(\mathbf{k}) e^{i \mathbf{k}\cdot(\sum_i \Delta n_i \mathbf{a}_i + \mathbf{r}_a - \mathbf{r}_b)}

ここで、副格子自由度は高々可算個なので、その他のパラメーター$\alpha$に含めてしまいます：$ (\mathbf{r}_a,\alpha) \rightarrow \alpha $。

また記号が複雑になってきたので $ t_{\alpha,\beta}(\Delta \mathbf{n}) = t_{\mathbf{R}_a(\Delta \mathbf{n}),\alpha \leftarrow 0, \mathbf{r}_b \beta} $ とおき直します。するとハミルトニアンは

H = \sum_{\mathbf{q}}\sum_{\alpha,\beta}\sum_{<\mathbf{R}_\alpha(\Delta \mathbf{n}), \mathbf{r}_\beta>} t_{\alpha,\beta}(\Delta \mathbf{n}) c^\dagger_\alpha(\mathbf{q}) c_\beta(\mathbf{q})  e^{i \mathbf{q} \cdot \Delta \mathbf{n} } e^{i \sum_{j} q_j\mathbf{b}_j \cdot (\mathbf{r}_\alpha - \mathbf{r}_\beta) }

となり、

H_{\alpha\beta}(\mathbf{q}) = \sum_{<\mathbf{R}_\alpha(\Delta \mathbf{n}), \mathbf{r}_\beta>} t_{\alpha,\beta}(\Delta \mathbf{n})  e^{i \mathbf{q} \cdot \Delta \mathbf{n} } e^{i \sum_{j} q_j\mathbf{b}_j \cdot (\mathbf{r}_\alpha - \mathbf{r}_\beta) }

とおけば、

H = \sum_{\mathbf{q}}\sum_{\alpha,\beta} H_{\alpha\beta}(\mathbf{q}) c^\dagger_\alpha(\mathbf{q}) c_\beta(\mathbf{q})

です。
ちなみにユニタリ行列

U_{phase} = Diag[(e^{i \sum_{j} q_j\mathbf{b}_j \cdot \mathbf{r}_\alpha })_\alpha] \\
= \begin{pmatrix}
e^{i \sum_{j} q_j\mathbf{b}_j \cdot \mathbf{r}_1} &  &  \\
 & \ddots & \\
 &  & e^{i \sum_{j} q_j\mathbf{b}_j \cdot \mathbf{r}_m}
\end{pmatrix}

を使って、

U_{phase}^\dagger H(\mathbf{q}) U_{phase} =  \sum_{<\mathbf{R}_\alpha(\Delta \mathbf{n}), \mathbf{r}_\beta>} t_{\alpha,\beta}(\Delta \mathbf{n})  e^{i \mathbf{q} \cdot \Delta \mathbf{n} }

とユニタリ変換したものを代わりに考えてもOKです。

2. 久保公式

次に応答係数を計算します。詳細は省きます。
まずハミルトニアン$H$で表されるシステムにおける物理量$X$の熱力学的平均は

<X> = <X>(\beta) = \frac{Tr(Xe^{-\beta H})}{Tr(e^{-\beta H})}

と書けます。ここで$\beta=1/T$は逆温度です($k_B=1$としています)。
また、物理量は第二量子化して

X = \sum_{\mathbf{r}}\sum_{\alpha,\beta} X_{\alpha\beta}c^\dagger_\alpha(\mathbf{r}) c_\beta(\mathbf{r}) \\
= \sum_{\mathbf{k}}\sum_{\alpha,\beta} X_{\alpha\beta}c^\dagger_\alpha(\mathbf{k}) c_\beta(\mathbf{k})

と書けます。

時間に依存する外場$\mathbf{E}(t)$を加えることで、ハミルトニアンに$\Delta H(t) = \sum_\mathbf{r} \mathbf{A}(\mathbf{r}) \cdot \mathbf{E}(t)$だけ変化が生じたとします。このときの物理量$\mathbf{B(\mathbf{r},t)}$の熱力学的平均は$E$の1次のオーダーで

<B_i>(\mathbf{r},t) = <B_i(\mathbf{r})>_0 -i \sum_j\sum_\mathbf{r'} \int_{-\infty}^t dt' <[B_i(\mathbf{r}),A_j^\dagger(\mathbf{r'})](t)>E_j(t)

と表せます($i,j=x,y,z$)。$ \langle B_i(\mathbf{r}) \rangle_0 $は０次の量で、$E=0$とした場合の熱力学平均です。よって$\mathbf{E}$による$\mathbf{B}$の期待値の変化分は

\Delta <\mathbf{B}>(\mathbf{r},t) = <\mathbf{B}>(\mathbf{r},t) - <\mathbf{B}(\mathbf{r})>_0 =  \chi_{BA}(\mathbf{r},t)\mathbf{E}(t)

とすれば、応答係数

\chi_{BA}(\mathbf{r},t) = -i \sum_\mathbf{r'}\int_{-\infty}^t dt' <[B_i(\mathbf{r-r'}),A_j^\dagger(0)](t)>

を得られます(久保公式)。
なお、$O(t)$は演算子$ O=[B,A] $のハイゼンベルグ表示です。
さらに周期性よりフーリエ変換して

\chi_{BA}(\mathbf{q},t) = -i \int_{-\infty}^t dt' <[B_i(\mathbf{q}),A_j^\dagger(\mathbf{-q})](t)>

となります。

3.　応答係数のレーマン表示

では応答係数の具体的な形を計算していきましょう。

温度グリーン関数

第二量子化してあるので基本的に$ < c^\dagger c > $のような量を計算していくことになります。そこで温度グリーン関数 $G$を導入します。

c_a(\tau) = e^{\tau H} c_a e^{-\tau H} \\
G_{ab}(\tau) = -i <T_\tau[c_a(\tau) c^\dagger_b] > \\
= -i<c_a(\tau) c^\dagger_b>\theta(\tau) -i<c^\dagger_b c_a(\tau)>\theta(-\tau) \\
= \begin{cases}
-i<c_a(\tau) c^\dagger_b> & (\tau>0) \\
-i<c^\dagger_b c_a(\tau)> & (\tau<0)
\end{cases}

ここで$\theta$は階段関数です。
実はこの温度グリーン関数には周期性があって

G_{ab}(\tau + \beta)= \begin{cases}
G_{ab}(\tau) & (ボゾン) \\
-G_{ab}(\tau) & (フェルミオン)
\end{cases}

となります。だからフーリエ級数展開して

G_{ab}(\tau) = \frac{1}{\beta}\sum_{n}G_{ab}(i\omega_n)

と書けます。ここで$\omega_n$ $(n=\cdots,-2,-1,0,1,2,\cdots)$は松原周波数といって

\omega_n = \begin{cases}
\frac{2n\pi}{\beta} & (ボゾン) \\
\frac{(2n+1)\pi}{\beta} & (フェルミオン)
\end{cases}

また、$H = \sum_{\mathbf{q}}\sum_{\alpha,\beta} H_{\alpha\beta}(\mathbf{q}) c^\dagger_\alpha(\mathbf{q}) c_\beta(\mathbf{q})$と書けるとき、$G(i\omega_n)$は

G_{ab}(\mathbf{q};i\omega_n) = G_{(a,\mathbf{q})(b,\mathbf{q})}(i\omega_n) = [i\omega_n - H(\mathbf{q})]^{-1}|_{ab}

と逆行列計算で求められます($\mathbf{q}$についてはすでにブロック対角化されています)。

よって$H(\mathbf{q})$の固有値が$E_n(\mathbf{q})$、固有ベクトルを並べたユニタリ行列が$U(\mathbf{q}) = [\mathbf{u}_1(\mathbf{q}), \mathbf{u}_2(\mathbf{q}), \cdots]$とすると、

G_{ab}(\mathbf{q};i\omega_n) = \sum_{n}U^*_{an}(\mathbf{q}) U_{bn}(\mathbf{q}) \frac{1}{i\omega_n - E_n(\mathbf{q})}

と書けます。

レーマン表示

突然ですが、

\chi_{BA}(\mathbf{q},\tau) =  \int_0^\beta d\tau <[B_i(\mathbf{q},\tau),A_j^\dagger(-\mathbf{q},0)]> 

を計算します。逆温度についてフーリエ級数展開しますが、$<[B,A]>$のフーリエ展開には粒子の種類に関係なくボゾンの松原周波数$\omega_l$を使います。これはウィック分解や留数積分を使って

\chi_{BA}(\mathbf{q},i\omega_l) = \frac{1}{\beta L^3} \sum_{i_1\cdots i_4}\sum_{\mathbf{k},\omega_n} B_{i_1i_2} A_{i_3i_4}G_{i_3i_1}(\mathbf{k},i\omega_n)G_{i_2i_4}(\mathbf{k+q},i\omega_n+i\omega_l)\\
= \frac{1}{L^3}\sum_{i_1\cdots i_4}\sum_\mathbf{k}\sum_{j_1,j_2} B_{i_1i_2} A_{i_3i_4} U^*_{i_1j_1}(\mathbf{k})U_{i_3j_1}(\mathbf{k})U^*_{i_4j_2}(\mathbf{k+q})U_{i_2j_2}(\mathbf{k+q}) \frac{ f[E_{j_2}(\mathbf{k+q})] - f[E_{j_1}(\mathbf{k})] }{ E_{j_1}(\mathbf{k}) - E_{j_2}(\mathbf{k+q})  + i\omega_l } \\
=  \sum_{i_1\cdots i_4} B_{i_1i_2} A_{i_3i_4} \Phi_{i_1 \cdots i_4}(\mathbf{q},i\omega_l)

となります。$f(E)$はフェルミ/ボーズ分布です。
これは解析接続して$i\omega_n \rightarrow \omega + i\delta$とすると$\chi(\mathbf{q},t)$のフーリエ変換に一致して

\chi_{BA}(\mathbf{q},\omega) = \lim_{\delta \rightarrow +0}\chi_{BA}(\mathbf{q},\omega + i\delta)

です。直流成分の場合は $\omega \rightarrow 0$とすればよいです。

#電気伝導度について
電流演算子は

J = e\frac{\partial \hat{r}}{\partial t}

と時間微分が入っています。このような場合は時間についてフーリエ変換する際に部分積分が必要となり最終的な結果が少し変わってきます。

電流演算子

位置演算子は

\hat{r} = i\frac{\hbar}{m}\frac{\partial}{\partial \mathbf{k}}

なのでハイゼンベルグ方程式より

J(\mathbf{k}) = e\frac{\partial \hat{r}}{\partial t} = \frac{e}{m}[\frac{\partial}{\partial \mathbf{k}},H(\mathbf{k})] = \frac{e}{m} \frac{\partial H(\mathbf{k})}{\partial \mathbf{k}} \\

J_{ab}(\mathbf{k}) =  \frac{e}{m} \sum_{<\mathbf{R}_\alpha(\Delta \mathbf{n}), \mathbf{r}_\beta>} t_{\alpha,\beta}(\Delta \mathbf{n}) [R_\alpha(\Delta \mathbf{n})-\mathbf{r}_\beta] e^{i \mathbf{q} \cdot \Delta \mathbf{n} } e^{i \sum_{j} q_j\mathbf{b}_j \cdot (\mathbf{r}_\alpha - \mathbf{r}_\beta) }

となります（$U_{phase}(\mathbf{k})$の$\mathbf{k}$微分も含まれているので、この部分も考える）。

レーマン表示

電気伝導度は

\sigma(\mathbf{q},\omega + i\delta) = -\frac{1}{i\omega}[ \chi_{JJ}(\mathbf{q},\omega+i\delta) - \chi_{JJ}(\mathbf{q},i\delta)]

となり、直流の場合は

\sigma(\mathbf{q},\omega + i\delta) = i\frac{ \partial \chi_{JJ}(\mathbf{q},\omega+i\delta) }{\partial \omega}

です。

#　注意点

• $\delta \rightarrow +0$、$\omega \rightarrow 0$など極限操作がありますが、この極限の順番によって結果が異なるので気を付けてください。
• 伝導度などはブリリュアンゾーンを$-\mathbf{b}_i/2$ ~ $\mathbf{b}_i/2$でとるか、$0$ ~ $\mathbf{b}_i$でとるかで結果が変わってきます。後者でやります。

参考資料

• 伏屋雄紀、福山秀敏「久保公式とグリーン関数法の実践的基礎」（固体物理誌上セミナー）

#関連
https://qiita.com/hiro949/items/ee7ae41d1aca8903f4ba