前文我们已经知道，当 \(H=H_{0}+H_{I}\) 且通常 \([H_{0},H_{I}]\neq0\) 时，我们同样的将虚时分为 \(L_{\tau}\) 份： \(\beta=L_{\tau} \Delta \tau\) ，且做Trotter分解：
\[ \begin{aligned} Z &=\operatorname{Tr}\left[e^{-\beta \hat{H}}\right] \\ &=\operatorname{Tr}\left[\left(e^{-\Delta \tau \hat{H}}\right)^{L_{\tau}}\right] \end{aligned} \]，其中 \(e^{-\Delta \tau \hat{H}} \approx e^{-\Delta \tau \hat{H}_{I}} e^{-\Delta \tau \hat{H}_{0}}\) ，
即：
\[ \begin{aligned} 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{aligned} \]而在前文我们已经知道：对于Hubbard U项，我们有：
\[ \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} \]其中 \(\gamma=\frac{1}{2} e^{-\Delta \tau U / 4}, \quad \cosh (\alpha)=e^{\Delta \tau U / 2}\) 。 \(l\) 是虚时方向的label。 此时\(N\) 为格点数 \(L\times L\) 。
为了形式简便，我们有如下记号： \(\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}+h . c .\) ，其中 \(T\) 为只在最近邻处有矩阵元 \(-t\) 的矩阵。而实际上更简便的写法是只写出自旋up或down的部分。之后的符号 \(T\) 本文中只表示大矩阵中的分块的小矩阵，which is：
\(\boldsymbol{c}_{\uparrow}^{\dagger} T \boldsymbol{c}_{\uparrow}=-t \sum_{\langle i j\rangle } \hat{c}_{i \uparrow}^{\dagger} \hat{c}_{j \uparrow}+h . c .\) ， \(\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 .\)
此时配分函数写作：
\[ \begin{aligned} 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, l} \hat{n}_{i \uparrow}}\right) \left( e^{-\Delta \tau \boldsymbol{c}_{\uparrow}^{\dagger} T \boldsymbol{c}_{\uparrow}}\right) \left(\prod_{i} e^{-\alpha s_{i, l} \hat{n}_{i \downarrow}}\right) \left( e^{-\Delta \tau \boldsymbol{c}_{\downarrow}^{\dagger} T \boldsymbol{c}_{\downarrow}}\right)\right] \right\} \end{aligned} \]显然可以看到，e指数上的形式关于自旋上下是分块对角的，所以可以写成更为简便的形式，而根据 \(\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)\) ，我们知道。对费米子部分求Trace得到的结果，等效为相应矩阵的exp的行列式。此时有：
\(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]\)
其中：
\(\begin{array}{l} \mathbf{B}^{\uparrow}(l_2 \Delta \tau, l_1 \Delta \tau)=\prod_{l=l_{1}+1}^{l_{2}}e^{\alpha \operatorname{Diag}\left(\vec{S}_{l}\right)} e^{-\Delta \tau T} \\ \mathbf{B}^{\downarrow}(l_2 \Delta \tau, l_1 \Delta \tau)=\prod_{l=l_{1}+1}^{l_{2}}e^{-\alpha \operatorname{Diag}\left(\vec{S}_{l}\right)} e^{-\Delta \tau T} \end{array}\) ，
其中 \(\operatorname{Diag}\left(\vec{S}_{l}\right)\) 表示对角元分别为 \(s_{i,l}\) 的矩阵。
现在我们看到，我们将费米子求trace求掉了，配分函数的求和部分，只是关于一个 \(L \times L \times L_{\tau}\) 的伊辛辅助场的求和。我们暂且不管辅助场的事情，就和伊辛模型一样，在伊辛模型的时候，我们给定某一个伊辛场的构型，得到该构型所对应的权重（概率）。此时同样的，给定 \(L \times L \times L_{\tau}\) 的伊辛场的取值，我们有相应的权重：
\(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}\)
其中加 \(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的时候，真正影响计算的是权重的比值，所以相同的常数因子我们可以略去，重新写权重为：
\(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]\)
现在我们看到，给定了某个具体的辅助场：2+1的伊辛场的构型，我们都能算出来该构型的权重（虽然很麻烦）。我们需要算每个 \(B\) ，此时需要算一个大矩阵的exp。矩阵的exp的计算通常就要对角化来算了，利用线代就学的 \(e^{-A}=e^{-U^\dagger \mathrm{diag} (\lambda) U}=U^\dagger \mathrm{Diag} (e^{-\lambda_i}) U\) 。
此时已经是很费力的操作了，然而我们还要将得到的 \(2 L_\tau\) 矩阵个矩阵乘起来（在乘的过程中可能会导致数值的误差）。乘完之后，还要加上一个单位阵，再算这个巨大的 \(N \times N\) 矩阵的行列式。这又是一项费力的操作。而再进行了这么多操作之后，现在只计算出来了一个数：该构型所对应的权重，而当你尝试更新到新的构型的时候，想算新构型的权重，trivial的想，要重新分别计算各个 \(B\) 再乘起来再算行列式…这显然是很可怕的计算量。而 \(2 L_\tau\) 个奇异性可能很大的矩阵直接相乘显然很容易造成数值上的误差。
这些问题将在之后逐一进行解决，但至少你现在可以无比暴力的来进行Hubbard Model的DQMC的计算了（应该没有人这么暴力的来做233毕竟大家的机时费都不是大风刮来的，时间上也等不起）。
去年买的辣鸡暗影精灵5，在家修不好跑去郑州修，赶在一年保修期结束前，修了俩星期终于搞定了…命途多舛的砖头本哭哭。