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 的计算了。