DQMC Note 2:从Hubbard Model开始的DQMC

为什么叫DQMC?

DQMC主要用来计算有限温的费米子系统。在实际运作的过程中,其中只有二费米子的部分。那么对于Hubbard Model这种含有四费米相互作用的哈密顿量,我们通常采取Hubbard-Stratonovich变换的方式,引入辅助场进行decouple。然后就变成了二费米子的形式,多了一个辅助场。(后面我们会知道不同的decouple的方式可能会导致符号问题的有无)

DQMC之所以叫行列式蒙特卡洛,就在于他将求Trace转化为了求一个矩阵的行列式。如我们要求形如以下形式的Trace:

\[\begin{equation}\operatorname{Tr}\left[e^{-\sum_{i, j} \hat{c}_{i}^{\dagger} A_{i, j} \hat{c}_{j}}\right]=\operatorname{Det}\left[\mathbf{1}+e^{-\mathbf{A}}\right]\end{equation}\]

其中 $A_{i, j}$ 为矩阵 $A$ 的元素。我们可以看到求Trace的操作转化为了求一个单位阵  $\mathbf{1}$  加上矩阵   $e^{-\mathbf{A}}$  的行列式。

证明的方式多种多样,比如我们可以把e指数上的部分进行对角化,Trace即为对于  $n_{i}=0,1$  的所有可能性求和:

\[\begin{equation}\begin{aligned}&\operatorname{Tr}\left[e^{-\sum_{i, j} \hat{c}_{i}^{\dagger} A_{i, j} \hat{c}_{j}}\right]\\ =&\operatorname{Tr}\left[e^{-\sum_{k} \epsilon(k) \hat{c}_{k}^{\dagger} \hat{c}_{k}}\right]\\ =&\Pi_{k} \operatorname{Tr}\left[e^{-\epsilon(k) \hat{c}_{k}^{\dagger} \hat{c}_{k}}\right]\\ =&\Pi_{k}\left[1+e^{-\epsilon(k)}\right]\end{aligned}\end{equation}\]

行列式就是本征值的乘积嘛,上述等式显然成立。

或者如果觉得e指数上带着算符不太直接,可以进行展开:

\[\begin{equation}\exp \left(-c_{k}^{\dagger} \epsilon(k) c_{k}\right)=1-c_{k}^{\dagger} \epsilon(k) c_{k}+\frac{1}{2 !} c_{k}^{\dagger} \epsilon(k) c_{k} c_{k}^{\dagger} \epsilon(k) c_{k}-\frac{1}{3 !} \cdots\end{equation}\]

再利用  $c_{k}^{\dagger} c_{k}+c_{k} c_{k}^{\dagger}=1$ 和 $c_{k} c_{k}|\psi\rangle=0$  ,有:

\[\begin{equation}\begin{aligned}
\exp \left(-c_{k}^{\dagger} \epsilon(k) c_{k}\right) &=1-c_{k}^{\dagger} \epsilon(k) c_{k}+\frac{1}{2 !} c_{k}^{\dagger} \epsilon(k)^{2} c_{k}-\frac{1}{3 !} \cdots \\
&=1+\left[e^{-\epsilon(k)}-1\right] c_{k}^{\dagger} c_{k}
\end{aligned}\end{equation}\]

(聪明的同学显然早就知道费米子的话 $\hat{n}_{k}^{m}=\hat{n}_{k}$ 了...)

而对DQMC重要的是,上式的推广形式。对于新的形式我们依旧能进行对角化 PhysRevB.31.4403

\[\begin{equation}e^{-\sum_{i, j} c_{i}^{\dagger} A_{i, j} c_{j}} e^{-\sum_{i, j} c_{i}^{\dagger} B_{i, j} c_{j}}=e^{-\sum_{\nu} c_{\nu}^{\dagger} l_{\nu} c_{\nu}}\end{equation}\]

从单粒子态开始: $e^{-A} e^{-B}$ 作用在 $|\psi\rangle=c_{\nu}^{\dagger}|0\rangle$  :

\[\begin{equation}e^{-\sum_{i, j} c_{i}^{\dagger} A_{i, j} c_{j}} e^{-\sum_{i, j} c_{i}^{\dagger} B_{i, j} c_{j}}|\psi\rangle=\left(e^{-A} e^{-B}\right)_{\nu \nu} c_{\nu}^{\dagger}|0\rangle=e^{-l_{\nu}} c_{\nu}^{\dagger}|0\rangle\end{equation}\]

考虑一个多粒子态比如二粒子态 $|\phi\rangle=c_{\mu_{1}}^{\dagger} c_{\mu_{2}}^{\dagger}|0\rangle$ ,有:

\[\begin{equation}\begin{aligned}
e^{-c_{i}^{\dagger} B_{i j} c_{j}}|\phi\rangle &=\prod_{\mu}\left[1+\left(e^{-B_{\mu}}-1\right) c_{\mu}^{\dagger} c_{\mu}\right] c_{\mu_{1}}^{\dagger} c_{\mu_{2}}^{\dagger}|0\rangle \\
&=e^{-B_{\mu_{1}}} e^{-B_{\mu_{2}}} c_{\mu_{1}}^{\dagger} c_{\mu_{2}}^{\dagger}|0\rangle
\end{aligned}\end{equation}\]

因此我们有:

\[\begin{equation}\operatorname{Tr}\left[e^{-\sum_{i, j} c_{i}^{\dagger} A_{i, j} c_{j}} e^{-\sum_{i, j} c_{i}^{\dagger} B_{i, j} c_{j}}\right]=\operatorname{Det}\left(1+e^{-\mathbf{A}} e^{-\mathbf{B}}\right)\end{equation}\]

对于更一般的情况,形式不言自明。

Hubbard Model 和 Hubbard-Stratonovitch 变换

我们知道对于自由费米子系统,可以解析的直接求出松原格林函数等信息:

\[\begin{equation}G\left(i \omega_{n}\right)=\frac{1}{i \omega_{n}-\left(\epsilon_{k}-\mu\right)}\end{equation}\]

而对于像Hubbard Model这样的含有四费米子相互作用的模型,我们需要离散的Hubbard-Stratonovitch分解,通过引入辅助场来进行蒙卡的求解。

历史上,此类DQMC最早是用来解决玻色费米耦合的模型, 又叫Blankenbecler-Scalapino-Sugar(BSS)算法, 而之后1983年Hirsch提出了针对onsite库伦相互作用的Hubbard模型的HS变换,即通过引入辅助场,将相互作用decouple为无相互作用费米子和辅助场的奉合, 才逐渐的能 够简单的对相互作用费米子的系统进行DQMC的计算。

Hubbard Model的哈密顿量写作:

\[\begin{equation}
\hat{H}=-t \sum_{\langle i j\rangle \sigma} \hat{c}_{i \sigma}^{\dagger} \hat{c}_{j \sigma}+h . c .+U \sum_{i}\left(\hat{n}_{i \uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i \downarrow}-\frac{1}{2}\right)
\end{equation}\]

我们通常会把虚时分为 $M$  份:  $\beta=M \Delta \tau$  。然后再做Trotter分解。这样做的原因是因为$\hat{H}_{0}=-t \sum_{\langle i j\rangle \sigma} \hat{c}_{i \sigma}^{\dagger} \hat{c}_{j \sigma}+h.c.$ 和 $\hat{H}_{I}=U \sum_{i}\left(\hat{n}_{i \uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i \downarrow}-\frac{1}{2}\right)$ 不对易。

我们想对  $e^{-\Delta \tau \hat{H}_{I}}$  做HS变换, 于是就要把他分离出来:

\[\begin{equation} e^{-\Delta \tau \hat{H}}=e^{-\Delta \tau\left(\hat{H}_{0}+\hat{H}_{I}\right)}=e^{-\Delta \tau \hat{H}_{0}} e^{-\Delta \tau \hat{H}_{I}}+\mathcal{O}\left[(\Delta \tau)^{2}\right]\end{equation}\]

至于为什么是 $\mathcal{O}\left[(\Delta \tau)^{2}\right]$  的误差,有多种理解方式比如从BCH公式出发, 这里不再多说了,留作QM习题

那么到底啥是HS变换呢? 某种意义上说,就是逆的高斯积分。

在大一的热学课上我们应该就学过高斯积分的公式:

\[\begin{equation} \quad e^{A^{2} / 2}=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{+\infty} \mathrm{d} \phi e^{-\frac{\phi^{2}}{2}-\phi A}\end{equation}\]

而此时我们知道, Hubbard U项可以写成平方的形式:

\[\begin{equation} \quad H_{U}=-\frac{U}{2}\left(n_{\uparrow}-n_{\downarrow}\right)^{2}+\frac{U}{4}\end{equation}\]

那么, 此时把  $A^{2}=\Delta \tau U\left(n_{\uparrow}-n_{\downarrow}\right)^{2}$  代入前面的高斯积分的公式即可。

这样我们发现通过引 入  $\phi,$ 我们把 $A^{2}$  的形式变为了  $A \phi$  的形式, 是  $A$  的一次项。

有趣的是,对于Hubbard来说,我们有一种更巧妙的decouple方式:

\[\begin{equation} e^{-\Delta \tau H_{U}}=\gamma \sum_{s=\pm 1} e^{\alpha s\left(n_{\uparrow}-n_{\downarrow}\right)}\end{equation}\]

其中  $\gamma=\frac{1}{2} e^{-\Delta \tau U / 4}, \; \cosh (\alpha)=e^{\Delta \tau U / 2}$ 

系数可以通过考虑  $e^{-\Delta \tau H_{U}}$  在单个格点的希尔伯特空间得到:

\[\begin{equation}
\begin{aligned}
e^{-\Delta \tau U / 4}|0\rangle &=2 \gamma|0\rangle \\
e^{-\Delta \tau U / 4}|\uparrow \downarrow\rangle &=2 \gamma|\uparrow \downarrow\rangle \\
e^{\Delta \tau U / 4}|\uparrow\rangle &=2 \gamma \cosh (\alpha)|\uparrow\rangle \\
e^{\Delta \tau U / 4}|\downarrow\rangle &=2 \gamma \cosh (\alpha)|\downarrow\rangle
\end{aligned}
\end{equation}\]

这种方式破坏了自旋的  $S U(2)$  对称性, 我们可以选择

\[\begin{equation} e^{-\Delta \tau H_{U}}=\tilde{\gamma} \sum_{s=\pm 1} e^{i \tilde{\alpha} s\left(n_{\uparrow}+n_{\downarrow}-1\right)}\end{equation}\]

其中  $\cos (\tilde{\alpha})=e^{-\Delta \tau U / 2}, \quad \tilde{\gamma}=\frac{1}{2} e^{\Delta \tau U / 4},$  代价就是引入了复数。

此处Hubbard Model的情况比较特殊才可以通过引入  $\pm 1$  的辅助场来进行decouple,更一般的 四费米相互作用(要求厄米),我们有更加一般的HS变换的方式, 如引入  $\pm 1$,$\pm 2$ 的辅助场:

\[\begin{equation}
e^{\Delta \tau \lambda A^{2}}=\sum_{l=\pm 1,\pm 2} \frac{\gamma(l)}{4} e^{\sqrt{\Delta \tau \lambda} \eta(l) O} +\mathcal{O}\left(\Delta \tau^{4}\right)\end{equation}\]
, 其中
\[\begin{equation}\begin{aligned}
\gamma(\pm 1)=1+\sqrt{6} / 3, & \qquad \gamma(\pm 2)=1-\sqrt{6} / 3 \\
\eta(\pm 1)=\pm \sqrt{2(3-\sqrt{6})}, &\qquad \eta(\pm 2)=\pm \sqrt{2(3+\sqrt{6})}
\end{aligned}\end{equation}\]

这种方式不是严格的,有和  $\Delta \tau$  有关的误差, 对于观测量来说, 带来的误差是正比于  $\Delta \tau^{3}$  的。
而现在的这些系数是怎么求解得到的呢? 感觉到处都直接么结果很少有人动手说说盾来的,这里我们不妨设:

\[\begin{equation}\gamma(1)=\gamma(-1)=a, \quad \gamma(2)=\gamma(-2)=b, \quad \eta(1)=\sqrt{c}=-\eta(1), \quad \eta(2)=\sqrt{d}=-\eta(2)\end{equation}\]

然后将两边都泰勒展开到  $O\left(p^{4}\right)$ ,有:

\[\begin{equation}\begin{array}{c}
1=\frac{1}{2}(a+b) \\
1=\frac{1}{4}(a c+b d) \\
\frac{1}{2}=\frac{1}{48}\left(a c^{2}+b d^{2}\right) \\
\frac{1}{6}=\frac{1}{1440}\left(a c^{3}+b d^{3}\right)
\end{array}\end{equation}\]

解得:

\[\begin{equation}\begin{array}{l}
a=1+\sqrt{6} / 3, \quad b=1-\sqrt{6} / 3 \\
c=2(3-\sqrt{6}), \quad d=2(3+\sqrt{6})
\end{array}\end{equation}\]

(我们后面会知道,不同的decouple方式,可能会导致符号问题的有无。所谓的符号问题此处就是对于辅助场的某个构型,对应的权重为复数或负数。

下一节将会推导一下引入辅助场之后,对于辅助场的具体构型,怎么计算相应的Weight。

Hubbard Model的权重

前文我们已经知道,  $\quad$ 当 $H=H_{0}+H_{I}$  且通常  $\left[H_{0}, H_{I}\right] \neq 0$  时,我们同样的将虚时分为  $L_{\tau}$  份:  $\quad \beta=L_{\tau} \Delta \tau,$  且做Trotter分解:

\[\begin{equation}
\begin{aligned}
Z &=\operatorname{Tr}\left[e^{-\beta \hat{H}}\right] \\
&=\operatorname{Tr}\left[\left(e^{-\Delta \tau \hat{H}}\right)^{L_{\tau}}\right] \cdot \text { 其中 } e^{-\Delta \tau \hat{H}} \approx e^{-\Delta \tau \hat{H}_{I}} e^{-\Delta \tau \hat{H}_{0}}
\end{aligned}
\end{equation}\]

即:

\[\begin{equation}Z=\operatorname{Tr}\left[e^{-\Delta \tau \hat{H}_{I}} e^{-\Delta \tau \hat{H}_{0}} e^{-\Delta \tau \hat{H}_{I}} e^{-\Delta \tau \hat{H}_{0}} \cdots e^{-\Delta \tau \hat{H}_{I}} e^{-\Delta \tau \hat{H}_{0}}\right]
\end{equation}\]

而在前文我们已经知道:对于Hubbard U项, 我们有:

\[\begin{equation}
\begin{aligned}
e^{-\Delta \tau \hat{H}_{I}} &=\prod_{i} e^{-\Delta \tau U\left(\hat{n}_{i \uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i \downarrow}-\frac{1}{2}\right)} \\
&=\prod_{i} \gamma \sum_{s_{i, l}=\pm 1} e^{\alpha s_{i, l}\left(\hat{n}_{i \uparrow}-\hat{n}_{i \downarrow}\right)} \\
&=\gamma^{N} \sum_{s_{i, l}=\pm 1}\left(\prod_{i} e^{\alpha s_{i, l} \hat{n}_{i \uparrow}} \prod_{i} e^{-\alpha s_{i, l} \hat{n}_{i \downarrow}}\right)
\end{aligned}
\end{equation}\]

其中  $\gamma=\frac{1}{2} e^{-\Delta \tau U / 4}, \quad \cosh (\alpha)=e^{\Delta \tau U / 2}$  。  $l$  是虚时方向的label。 此时  $N$  为格点数 $L \times L$

为了形式简便, 我们有如下记号:  $\quad \hat{H}_{0}=\boldsymbol{c}^{\dagger} T \boldsymbol{c}=-t \sum_{\langle i j\rangle \sigma} \hat{c}_{i \sigma}^{\dagger} \hat{c}_{j \sigma}+\boldsymbol{h} . \boldsymbol{c} .,$  其中  $T$  为只在
最近邻处有矩阵元  $-t$  的矩阵。而实际上更简便的写法是只写出自旋up或down的部分。之后的符号 $T$ 本文中只表示大矩阵中的分块的小矩阵, which is:

\[\begin{equation}
\boldsymbol{c}_{\uparrow}^{\dagger} T \boldsymbol{c}_{\uparrow}=-t \sum_{\langle i j\rangle} \hat{c}_{i \uparrow}^{\dagger} \hat{c}_{j \uparrow}+\boldsymbol{h} . c . \quad \boldsymbol{c}_{\downarrow}^{\dagger} T \boldsymbol{c}_{\downarrow}=-t \sum_{\langle i j\rangle} \hat{c}_{i \downarrow}^{\dagger} \hat{c}_{j \downarrow}+h . c .
\end{equation}\]

此时配分函数写作:

\[\begin{equation}
Z =\sum_{s_{i, l}=\pm 1} \gamma^{N L_{\tau}} \operatorname{Tr}_{F}\left\{\prod_{l=1}^{L_{\tau}}\left[\left(\prod_{i} e^{\alpha s_{i,}, \hat{n}_{i \uparrow}}\right)\left(e^{-\Delta \tau_{\uparrow}^{\dagger} T c_{\uparrow}}\right)\left(\prod_{i} e^{-\alpha s_{i, l} \hat{n}_{i l}}\right)\left(e^{-\Delta \tau c^{\dagger} T c_{\downarrow}}\right)\right]\right\}
\end{equation}\]

显然可以看到,e指数上的形式关于自旋上下是分块的,所以可以写成更为简便的形式, 而根据

\[\begin{equation} \operatorname{Tr}\left[e^{-\sum_{i, j} c_{i}^{\dagger} A_{i, j} c_{j}} e^{-\sum_{i, j} c_{i}^{\dagger} B_{i, j} c_{j}}\right]=\operatorname{Det}\left(1+e^{-\mathbf{A}} e^{-\mathbf{B}}\right)\end{equation}\]

我们知道。对费米子部分求Trace得到的结果,等效为相应矩阵的exp的行列式。此时有:

\[\begin{equation}Z=\gamma^{N L_{\tau}} \sum_{\left\{s_{i, l}\right\}} \prod_{\sigma=\uparrow \downarrow} \operatorname{Det}\left[\mathbf{I}+\mathbf{B}^{\sigma}\left(L_{\tau} \Delta \tau,\left(L_{\tau}-1\right) \Delta \tau\right) \cdots \mathbf{B}^{\sigma}(\Delta \tau, 0)\right]\end{equation}\]

其中:

\[\begin{equation}\begin{array}{l}
\mathbf{B}^{\uparrow}\left(l_{2} \Delta \tau, l_{1} \Delta \tau\right)=\prod_{l=l_{1}+1}^{l_{2}} e^{\alpha \operatorname{Diag}\left(\vec{S}_{l}\right)} e^{-\Delta \tau T} \\
\mathbf{B}^{\downarrow}\left(l_{2} \Delta \tau, l_{1} \Delta \tau\right)=\prod_{l=l_{1}+1}^{l_{2}} e^{-\alpha \operatorname{Diag}\left(\vec{S}_{l}\right)} e^{-\Delta \tau T}
\end{array}\end{equation}\]

其中  $\operatorname{Diag}\left(\vec{S}_{l}\right)$  表示对角元分别为  $s_{i, l}$  的矩阵。

现在我们看到, 我们将费米子求trace求掉了,配分函数的求和部分,只是关于一个 $L \times L \times L_{\tau}$ 的伊辛辅助场的求和。我们暂目不管辅助场的事情, 就和伊辛模型一样,在伊辛 模型的时候,我们给定某一个伊辛场的构型,得到该构型所对应的权重(概率)。此时同样的,给 定  $L \times L \times L_{\tau}$  的伊辛场的取值,我们有相应的权重:

\[\begin{equation}W_{C}=\frac{\gamma^{N L_{\tau}} \prod_{\sigma=\uparrow \downarrow} \operatorname{Det}\left[\mathbf{I}+\mathbf{B}_{C}^{\sigma}\left(L_{\tau} \Delta \tau,\left(L_{\tau}-1\right) \Delta \tau\right) \cdots \mathbf{B}_{C}^{\sigma}(\Delta \tau, 0)\right]}{Z}\end{equation}\]

其中加 C 指的是将具体的Configuration的值代入得到相应的矩阵再进行操作,比如我们看  $B_{C}^{\uparrow}(\Delta \tau, 0)=e^{\alpha \operatorname{Diag}\left(\vec{S}_{1}\right)} e^{-\Delta \tau T}$  这一部分,我们发现后面的一半和伊辛场无关, 前面的一半和  $s_{i, 1}$  的取值是 +1 还是 -1 有关。

我们知道,在进行MCMC的时候,真正影响计算的是权重的比值,所以相同的常数因子我们可以 略去,重新写权重为:

\[\begin{equation}W_{C}=\prod_{\sigma=\uparrow \downarrow} \operatorname{Det}\left[\mathbf{I}+\mathbf{B}_{C}^{\sigma}\left(L_{\tau} \Delta \tau,\left(L_{\tau}-1\right) \Delta \tau\right) \cdots \mathbf{B}_{C}^{\sigma}(\Delta \tau, 0)\right]=\prod_{\sigma=\uparrow \downarrow} \operatorname{Det}\left[\mathbf{I}+ \mathbf{B}_{C}^{\sigma}(\beta, 0)\right]\end{equation}\]

现在我们看到,给定了某个具体的辅助场:2+1的伊辛场的构型,我们都能算出来该构型的权重 (周然很麻烦)。我们需要算每个 B, 此时需要算一个大矩阵的exp。矩阵的exp的计算通常就 要对角化来算了,利用线代就学的  $e^{-A}=e^{-U^{\dagger} \operatorname{diag}(\lambda) U}=U^{\dagger} \operatorname{Diag}\left(e^{-\lambda_{i}}\right) U$ 

此时已经是很费力的操作了,然而我们还要将得到的  $2 L_{\tau}$  矩阵个矩阵乘起来(在乘的过程中可能 会导致数值的误差)。乘完之后,还要加上一个单位阵,再算这个巨大的  $N \times N$  矩阵的行列 式。这又是一项费力的操作。而再进行了这么多操作之后,现在只计算出来了一个数:该构型所对 应的权重,而当你尝试更新到新的构型的时候,想算新构型的权重,trivial的想,要電新分别计算 各个 $B$ 再乘起来再算行列式...这显然是很可怕的计算量。而  $2 L_{\tau}$  个奇异性可能很大的矩阵直接相乘显然很容易造成数值上的误差。

这些问题将在之后逐一进行解决, 但至少你现在可以无比暴力的来进行Hubbard Model的DQMC 的计算了。