不純物アンダーソン模型を量子モンテカルロ法で解いてみよう（理論パート）

Last updated at 2025-02-28Posted at 2025-02-25

\newcommand{\tr}{\mathrm{Tr}}
\newcommand{\texp}{\mathrm{Texp}}
\newcommand{\ks}{{k\sigma}}
\newcommand{\Dtau}{{\Delta \tau}}
\newcommand{\flip}{\mathrm{flip}}

本記事では量子モンテカルロ法を紹介したいと思います。量子モンテカルロ法は温度方向に離散的なメッシュを切ってその方向の格子に対してイジング模型などのモンテカルロ法と同様に状態を更新していくやり方です。
今回はその一番簡単な例として不純物アンダーソン模型の不純物サイトのシステムを考えます。（１サイトしかないので、量子モンテカルロ法でも１次元で済むからです。）

不純物アンダーソン模型

バス(伝導電子$c$)と不純物中の電子$f$が相互作用する以下のモデルです：

\begin{align}
    H_{IAM} &= H_C + H_I + H_{hyb} + H_{U} \\
    H_C &= \sum_{\ks}  \epsilon_{\ks}c^\dagger_\ks c_\ks \\
    H_I &= \sum_{\sigma} \left( \epsilon_{f\sigma} - \mu +\frac{U}{2} \right) f^\dagger_\sigma f_\sigma \\
    &= \sum_\sigma \left( \epsilon_{f\sigma} - \mu +\frac{U}{2} \right) n_{f\sigma} \\
    H_{hyb} &=  \sum_{\ks} \left( V^*_k f^\dagger_\sigma c_\ks + V_kc^\dagger_\ks f_\sigma \right) \\
    H_U &=\frac{U}{2} \sum_{\sigma} f^\dagger_\sigma f_\sigma f^\dagger_{-\sigma} f_{-\sigma} - \frac{U}{2} \sum_\sigma  n_{f\sigma} \\
    &= \frac{U}{2} \sum_{\sigma} n_{f\sigma} n_{f -\sigma} - \frac{U}{2} \sum_\sigma  n_{f\sigma}
\end{align}

ここで、後々のために$U/2$だけポテンシャルをシフトした項を$H_U$に付けました。
また、外場のないハーフフィリングの場合は化学ポテンシャル$\mu$が以下のようになります。

\begin{align}
\epsilon_{f\sigma} = const. &= \epsilon \\
\epsilon_{f\sigma} - \mu +\frac{U}{2} &= 0
\end{align}

この分配関数は以下の補遺のように

伝導電子について部分トレースをとると以下のように書けます。

\begin{align}
Z_{eff} &= \tr \exp \left[ - \int_ 0^\beta d\tau \mathcal{A}_0(\tau) - \int_0^\beta d\tau H_U(\tau) \right] \\
\mathcal{A}_0(\tau) &= \sum_\sigma \int_ 0^\beta d\tau' \mathcal{K}_\sigma(\tau,\tau') \\
\mathcal{K}_\sigma(\tau,\tau') &= -[\mathcal{G}_\sigma(\tau - \tau')]^{-1} f^\dagger_\sigma(\tau) f_\sigma(\tau') \\
\mathcal{G}_\sigma(\tau - \tau') &= \left(\sum_n e^{-i\omega_n (\tau-\tau')} (i\omega_n - \epsilon_{f\sigma} - \Delta_\sigma(i\omega_n) ) \right)^{-1}
\end{align}

を考えることになります。

Hirsch-Fyeアルゴリズム

量子モンテカルロ法の一番簡単な方法です。

鈴木・トロッター分解

トロッターの公式 ($A,B$は非可換演算子)

e^{t(A+B)} = \lim_{n\rightarrow \infty} \left( e^{\frac{t}{n}A} e^{\frac{t}{n}B} \right)^n

を緩めて ($M \gg 1$)に対して

e^{t(A+B)} \simeq \left( e^{\frac{t}{M}A} e^{\frac{t}{M}B} \right)^M = e^{\frac{t}{M}A} e^{\frac{t}{M}B} e^{\frac{t}{M}A} e^{\frac{t}{M}B} \cdots e^{\frac{t}{M}A} e^{\frac{t}{M}B}

という近似式を使って分配関数を分解して計算します。

まず、$\Lambda \gg 1$をとって

\begin{align}
\beta &= \Lambda \Dtau \\
\int_0^\beta d\tau &\rightarrow \sum_{l=0}^{\Lambda} \Dtau
\end{align}

と差分近似すると、

\begin{align}
Z_{eff} &\simeq \tr \exp \left[ -\sum_{l} (\Dtau)^2 \mathcal{A}_0(l\Dtau) - \sum_{l} \Dtau H_U(l\Dtau) \right] \\
\mathcal{A}_0(l\Dtau) &\simeq \sum_{\sigma,l'} (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \\
H_U(\tau) &\simeq \frac{U}{2} \sum_{\sigma} \left[ n_{f\sigma}(l\Dtau) n_{f-\sigma}(l\Dtau) - n_{f\sigma}(l\Dtau) \right]
\end{align}

と書けます。よって、鈴木・トロッター分解して

Z_{eff} \simeq \tr \prod_l \exp \left[ -(\Dtau)^2 \mathcal{A}_0(l\Dtau) \right]\exp\left[ - \Dtau H_U(l\Dtau) \right]

が得られます。

ハバード・ストラトのビッチ変換

今$n_\sigma = 0$ or $1$であることを利用すると

\begin{align}
\Dtau \frac{U}{2} (n_\uparrow - n_\downarrow)^2 &= \Dtau \frac{U}{2}(\sum_\sigma n^2_\sigma - 2n_\uparrow n_\downarrow) \\
&= \Dtau \frac{U}{2} \sum_\sigma n_\sigma - \Dtau U n_\uparrow n_\downarrow \\
&= - \Dtau H_U
\end{align}

です。
ここで、

e^{\frac{\Dtau U}{2} (n_\uparrow - n_\downarrow)^2} = \sum_{s = \pm 1} w(s) e^{-\lambda s (n_\uparrow  - n_\downarrow)}

という形になるように$\lambda,w(s)$を決めます。

$(n_\uparrow, n_\downarrow) = (0,0),(0,1),(1,0),(1,1)$しかパターンはないので、

\begin{align}
1 &= w(1) + w(-1) \\
e^{\frac{\Dtau U}{2}} &= w(1) e^{\lambda} + w(-1) e^{-\lambda} \\
e^{\frac{\Dtau U}{2}} &= w(1) e^{-\lambda} + w(-1) e^{\lambda} \\
1 &= w(1) + w(-1)
\end{align}

が条件となり、これを解くと、

\begin{align}
w(1) &= w(-1) = \frac{1}{2} \\
\lambda &= \mathrm{arccosh}\left( e^{\frac{\Dtau U}{2}} \right) \sim \omicron(\sqrt{\Dtau})
\end{align}

です。以上から

\begin{align}
e^{- \Dtau H_U} &= e^{\frac{\Dtau U}{2} (n_\uparrow - n_\downarrow)^2} \\
&= \frac{1}{2} \sum_{s = \pm 1}  e^{-\lambda s (n_\uparrow  - n_\downarrow)} \\
&= \frac{1}{2} \sum_{s = \pm 1}  e^{-\lambda s \sum_\sigma \sigma n_\sigma}
\end{align}

を得て、

\begin{align}
Z_{eff} &\simeq \tr \prod_l \exp \left[ -\mathcal{A}_0(l\Dtau) \right]\exp\left[ - \Dtau H_U(l\Dtau)\right] \\
&= \tr \prod_l \exp \left[ -\mathcal{A}_0(l\Dtau) \right]\left[ \frac{1}{2} \sum_{s_l = \pm 1}  \exp \left(\lambda s_l \sum_\sigma \sigma n_\sigma(l\Dtau) \right) \right] \\
&= \tr \prod_l \left[ \frac{1}{2} \sum_{s_l = \pm 1} \exp \left[ -\mathcal{A}_0(l\Dtau) \right] \exp \left(\lambda s_l \sum_\sigma \sigma n_\sigma(l\Dtau) \right) \right] \\
&= \left(\frac{1}{2}\right)^{\Lambda+1} \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr \prod_l \exp \left[ -\sum_{\sigma,l'} (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \right] \exp \left(\lambda s_l \sum_\sigma \sigma n_\sigma(l\Dtau) \right) \\
&= \left(\frac{1}{2}\right)^{\Lambda+1} \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr \prod_{\sigma,l} \exp \left[ -\sum_{l'} (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \right] \exp \left(\lambda s_l \sigma \sum_{l'} \delta_{ll'}  f^\dagger_\sigma(l\Dtau) f_\sigma(l'\Dtau) \right) \\
&= \left(\frac{1}{2}\right)^{\Lambda+1} \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr \prod_{\sigma,l,l'} \exp \left[ - (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \right] \exp \left[\lambda s_l \sigma \delta_{ll'}  f^\dagger_\sigma(l\Dtau) f_\sigma(l'\Dtau) \right]
\end{align}

です。

ここでトロッター分解の逆を使うと、

\begin{align}
Z_{eff} &\simeq \left(\frac{1}{2}\right)^{\Lambda+1} \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr \exp \left[ -\sum_{\sigma,l.l'} (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \right] \exp \left[ -\sum_{\sigma,l'} (\Dtau)^2 \mathcal{K}_\sigma(l\Dtau,l'\Dtau) \right] e^{n_\sigma(l\Dtau)}
\end{align}

ベーカー・キャンベル・ハウスドルフの公式を使うと、

M =

トロッター分解の逆を使って

\begin{align}
Z_{eff} &\simeq \left(\frac{1}{2}\right)^{\Lambda+1} \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr prod_\sigma \exp \left[ \sum_{l,l'} \left\{ (\Dtau)^2 [\mathcal{G}_\sigma((l - l')\Dtau)]^{-1} + \delta_{l,l'}\lambda s_l \sigma \right\}f^\dagger_\sigma(l\Dtau) f_\sigma(l'\Dtau) \right] \\
 &\propto \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \tr \prod_\sigma \exp \left[ \sum_{l,l'} \ M_{l,l'}(\sigma,s_l) f^\dagger_\sigma(l\Dtau) f_\sigma(l'\Dtau) \right] \\
 M_{l,l'}(\sigma,s_l) &= (\Dtau)^2 [\mathcal{G}_\sigma((l - l')\Dtau)]^{-1} + \delta_{l,l'}\lambda s_l \sigma
\end{align}

を得られます。

より厳密な式

の式(218)によるとより厳密に計算すると

M_{l,l'}(\sigma,s_l) = (\Dtau)^2 [\mathcal{G}_\sigma((l - l')\Dtau)]^{-1} + \delta_{l,l'}(1-e^{\lambda \sigma s_l})

と書けるらしいですが、ここまで示すことはできませんでした。
導出過程が理解できたら更新したいと思います。

ここで再びガウス積分を使って、

\begin{align}
Z_{eff} &\simeq \sum_{\{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \prod_\sigma \int \mathcal{D} f_\sigma \mathcal{D} f^*_\sigma \exp \left[ \sum_{l,l'} \ M_{l,l'}(\sigma,s_l) f^*_\sigma(l\Dtau) f_\sigma(l'\Dtau) \right] \\
&= \sum_{\mathbf{s} = \{s_l\} \in (1,-1)^{\otimes (\Lambda + 1)}} \prod_\sigma \det{M(\sigma,\mathbf{s})}
\end{align}

となります。

更新方法

メトロポリスアルゴリズムを使います。
統計力学によるとシステムが取りうる状態$\omega \in \Omega$に対して、状態$\omega$になる確率$P(\omega)$と分配関数には

Z \propto \sum_{\omega \in \Omega} P(\omega)

という関係があります。
ここから今回はの場合、状態は補助場のアンザッツ：

\Omega = \{\pm 1\}^\Lambda

で各状態をとる確率が

P(\mathbf{s}) = \frac{\prod_\sigma \det{M(\sigma,\mathbf{s})}}{Z_{eff}}

と分かります。(実際の不純物系のスピン$\sigma$は１サイトしかないので特に乱択は使わず計算します。)

補助場$s_m$を反転させる操作

\flip: (s_0,\cdots,s_m,\cdots) \rightarrow (s_0,\cdots,-s_m,\cdots)

を考えると、この反転による確率重みは

\begin{align}
W\left(\flip(\mathbf{s}) \mid \mathbf{s}\right) &= \frac{P(\flip(\mathbf{s})) }{P(\mathbf{s}) } \\
&= \prod_\sigma \frac{\det{M(\sigma,\flip(\mathbf{s}))} }{ \det{M(\sigma,\mathbf{s})} } \\
&= \prod_\sigma \det{M(\sigma,\flip(\mathbf{s}))} \det{M^{-1}(\sigma,\mathbf{s})} \\
&= \prod_\sigma \det{\left[ M(\sigma,\flip(\mathbf{s})) M^{-1}(\sigma,\mathbf{s}) \right]} \\
&= \frac{\det{M(\uparrow,\flip(\mathbf{s}))} \det{M(\downarrow,\flip(\mathbf{s}))} }{ \det{M(\uparrow,\mathbf{s})} \det{M(\downarrow,\mathbf{s})} }
\end{align}

です。ここで、

\begin{align}
M_{l,l'}(\sigma,\flip(\mathbf{s})) &= M_{l,l'}(\sigma,\mathbf{s}) - 2\delta_{l,l'}\delta_{l,m}\lambda s_l \sigma \\
&= M_{l,l'}(\sigma,\mathbf{s}) + \Delta^m_{l,l'}(\sigma,\mathbf{s}) \\
\det{\left[ M(\sigma,\flip(\mathbf{s})) M^{-1}(\sigma,\mathbf{s}) \right]} &= \det{\left[ (M(\sigma,\mathbf{s}) + \Delta^m(\sigma,\mathbf{s})  ) M^{-1}(\sigma,\mathbf{s}) \right]} \\
&= \det{\left[ 1 + \Delta^m(\sigma,\mathbf{s})  M^{-1}(\sigma,\mathbf{s}) \right]}\\
&= 1 + \Delta^m_{mm}(\sigma,\mathbf{s})  \left( M^{-1}(\sigma,\mathbf{s}) \right)_{mm}
\end{align}

です。(最後は$m$行、$m$列目についてブロック行列の性質を使いました。$m$列目は$(m,m)$の値以外は0であることに注意してください。)

ちなみに厳密な方でやると

\begin{align}
M_{l,l'}(\sigma,\flip(\mathbf{s})) &= M_{l,l'}(\sigma,\mathbf{s}) + 2\delta_{l,l'}\delta_{l,m}\sinh{(\lambda \sigma s_l)} \
\end{align}

Woodburyの公式

更新確率の重みを計算するには毎ステップ$M^{-1}(\sigma,\mathbf{s})$の計算が必要ですが、実は$\mathcal{O}(\Lambda^3)$のフルの逆行列計算は最初の１回だけで良くて２回目以降はWoodburyの公式を使って$\mathcal{O}(\Lambda^2)$で計算できます。

\begin{align}
M^{-1}(\sigma,\flip(\mathbf{s}))
&= \left[ M(\sigma,\mathbf{s}) + \Delta^m(\sigma,\mathbf{s}) \right]^{-1}\\
&= \left[ M(\sigma,\mathbf{s}) + \Delta^m_{mm}(\sigma,\mathbf{s})\mathbf{e}^T_m \otimes \mathbf{e}_m \right]^{-1}\\
&= M^{-1}(\sigma,\mathbf{s}) - \frac{\Delta^m_{mm}(\sigma,\mathbf{s})}{1 + \Delta^m_{mm}(\sigma,\mathbf{s}) \mathbf{e}^T_m M^{-1}(\sigma,\mathbf{s}) \mathbf{e}_m } M^{-1}(\sigma,\mathbf{s}) (\mathbf{e}^T_m \otimes \mathbf{e}_m) M^{-1}(\sigma,\mathbf{s}) \\
&= M^{-1}(\sigma,\mathbf{s}) - \frac{\Delta^m_{mm}(\sigma,\mathbf{s})}{1 + \Delta^m_{mm}(\sigma,\mathbf{s}) [ M^{-1}(\sigma,\mathbf{s}) ]_{mm} } (M^{-1}(\sigma,\mathbf{s})\mathbf{e}_m)^T \otimes (M^{-1}(\sigma,\mathbf{s})\mathbf{e}_m)
\end{align}

ただし$\mathbf{e}_m$は単位ベクトル

(e_m)_l = \delta_{lm}

です。

Reference

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up