第5章 · 量子蒙特卡洛入门#

本章我们将扩展上一章讨论的思想和技术，介绍几种用于量子多体系统蒙特卡洛模拟的基本方法。出于教学目的，我们仅考虑非零温度下的简单量子自旋-1/2模型。我们讨论的一些技术，例如带有局部更新和虚时离散化的路径积分蒙特卡洛方法，只是为介绍关键思想和概念，为后续讨论最先进的连续时间团簇算法和虫洞算法（worm algorithm）做准备。在后续章节中，我们将讨论更高级的量子自旋算法以及费米子和玻色子系统的零温和有限温算法。

5.1 经典表示#

目前，所有蒙特卡洛模拟都是在经典计算机上进行的。使用这些计算机，我们只能直接操作经典数，其乘积彼此对易。然而，量子力学处理的是算符，大多数算符的乘积彼此不对易。因此，我们需要为量子系统找到一种经典表示，将量子问题重新表述为经典数的形式。我们从一个简单直观的方法开始。

我们都熟悉量子力学的矩阵表述，特别是将哈密顿算符表示为一个厄米矩阵。这种表示是量子问题最基本的一种"经典化"。虽然一般来说它们的乘积不对易，但矩阵只是经典数的有序排列，经典计算机可以轻松处理其代数运算。特别是，经典计算机可以对这样的矩阵进行对角化，并给出所有特征值和特征向量。这些信息就"解决"了量子问题。

在感兴趣的量子多体问题中，相空间与经典情况一样，是许多局部相空间的直积。如果局部空间的大小为n，且我们有N个格点，则相空间随系统物理尺寸呈指数增长，为$n^{N}$。哈密顿矩阵的阶数等于相空间的大小，因此也呈现指数增长。由于计算机内存限制，我们只能获得、事实上也只能期望获得小尺寸格点的所有特征值和特征向量（第10.1节）。然而，小尺寸格点的解通常会受到边界效应的影响，往往不能反映热力学极限下的物理性质。因此，我们必须利用蒙特卡洛方法对非常大的相空间中的重要部分进行采样的能力。

还有其他两种经典化的方法：高温级数展开（Handscomb, 1962a）和费曼路径积分表述（Suzuki, 1976a,b）。费曼路径积分将d维量子问题转化为(d+1)维经典问题。d维中的原始哈密顿量作用于沿着新维度堆叠的层上，并有一组新的相互作用将各层耦合起来。在量子统计力学中，这个额外维度是虚时轴。

路径积分算法有离散时间和连续时间两种版本，这些算法通常用于模拟量子多体系统的有限温平衡性质。第一代路径积分算法使用离散化的虚时。为了消除离散化误差（该误差在虚时步长上呈多项式增长），我们必须外推到连续虚时极限。通常，我们可以在算法层面做到这一点，因此现代算法的模拟结果没有离散化误差。¹

高温展开导致一种经典表示，其在形式上类似于路径积分方法的离散虚时表示。一个重要的区别是，高温级数方法的离散化误差是指数小的，在大多数实际情况下可以忽略。

在第6章中，我们将详细讨论量子高温级数展开和路径积分表述。在作为量子蒙特卡洛算法初步发展的本章中，我们将简要介绍路径积分方法。

5.2 量子自旋#

脚注

¹ 然而，在处理自旋-1/2系统时，通过外推到零虚时步长来完全消除离散化误差可能很棘手，尤其是对于较大自旋S的情形。在S=1/2的特殊情况下，使用环回算法可以完全消除离散化误差；而对于较大的S，即使使用环回算法也可能存在残留的离散化误差。

现在我们着手建立第一个量子蒙特卡洛算法。我们从自旋-1/2量子伊辛模型和XY模型的离散虚时路径积分表示开始。通过取连续时间极限，我们最终用环回/团簇算法和虫洞算法取代了传统的离散时间(世界线)蒙特卡洛方法。通过这样做，我们实现了类似于我们在第4章对经典多体问题所展示的算法演变。我们的目标是揭示量子蒙特卡洛算法中使用的几个基本特征和技术。虽然我们的讨论是在特定的简单有限温度实例背景下进行的，但其中一些特征和技术对于零温问题也很有用。

5.2.1 纵向场伊辛模型#

经典一维伊辛模型 (1.3) 的量子版本由哈密顿量给出

$$ H = - \sum_{i = 1}^{N} \left( J^{z} S_{i}^{z} S_{i + 1}^{z} + H^{z} S_{i}^{z} \right) .\tag{5.1} $$

这里，$S_{i}^{z}$ 是格点 $i$ 处 $S = \frac{1}{2}$ 角动量算符的 z 分量，因此其特征值为 $\frac{1}{2}$ 或 $-\frac{1}{2}$。如果 $H^{z} = 0$，该哈密顿量在每一个格点 i 处执行变换 $S_{i}^{z} \rightarrow - S_{i}^{z}$ 时保持不变。如果 $H^{z} \neq 0$，只要 $H^{z}$ 的符号反转，它在该变换下也保持不变。因此，这个量子模型与经典模型具有相同的对称性。如下文将看到的，这个“量子模型”实际上等效于经典伊辛模型。

在量子统计力学中，配分函数为

$$ Z = \mathrm{Tr} \left[ \mathrm{e}^{- \beta H} \right] ,\tag{5.2} $$

其中 $\beta = 1 / k T$ 表示逆温度。(下文中，我们设玻尔兹曼常数 $k = 1$。) 迹可以用一组完备的正交基态来计算，我们将其记为 $\{| C \rangle \}$，

$$ \langle C | C^{\prime} \rangle = \delta ( C , C^{\prime} ) , \quad \sum_{C} | C \rangle \langle C | = I . $$

利用这样的基，(5.2) 变为

$$ Z = \sum_{C} \langle C | e^{- \beta H} | C \rangle . $$

如果 $| C \rangle$ 同时也是 $H$ 的本征态，即 $H | C \rangle = E ( C ) | C \rangle$，那么可得 $\langle C | e^{- \beta H} | C \rangle = e^{- \beta E ( C ) }$，因此

$$ Z = \sum_{c} e^{- \beta E ( C )} .\tag{5.3} $$

类似地，某个物理可观测量 X 的热期望值变为

$$ \langle X \rangle = {\frac{1} {Z}} \mathrm{Tr} \Big [ X e^{- \beta H} \Big ] = {\frac{1} {Z}} \sum_{C} \sum_{C^{\prime}} \langle C | X | C^{\prime} \rangle \langle C^{\prime} | e^{- \beta H} | C \rangle .\tag{5.4} $$

特别地，当 $\{| C \rangle\}$ 是哈密顿量的正交本征态时，我们得到

$$ \langle X \rangle = {\frac{1} {Z}} \sum_{C} X ( C ) e^{- \beta E ( C )} , $$

其中 $X ( C ) = \langle C | X | C \rangle$，这个表达式让人想起经典统计力学中的期望值。

因为知道 H 的本征能量和本征态等同于知道问题的解，所以将 (5.2) 简化为 (5.3) 通常只能是形式上的。然而，对于模型 (5.1)，本征态是

$$ | C \rangle = | s_{1} \rangle \otimes | s_{2} \rangle \otimes \dots \otimes | s_{N} \rangle = | s_{1} s_{2} \cdot \cdot \cdot s_{N} \rangle ,\tag{5.5} $$

其中 $s_{i} = \pm \frac{1}{2}$ 且 $| s_{i} \rangle$ 是 $S_{i}^{z}$ 的特征对。那么，特征值 $E ( C )$ 就是经典模型的能量，只是经典自旋 $s_{i}$ 的取值是 $\pm \frac{1}{2}$ 而不是 $\pm 1$，² 而对配分函数 (5.3) 进行采样的问题与经典伊辛问题相同。

5.2.2 横向场伊辛模型#

现在我们考虑一个不那么平凡的量子自旋模型——处于横向磁场中的自旋-1/2链，并通过离散化虚时区间显式推导费曼路径积分。横向场伊辛模型的哈密顿量为

$$ H = - \sum_{i = 1}^{N} \left( J^{z} S_{i}^{z} S_{i + 1}^{z} + H^{x} S_{i}^{x} \right) ,\tag{5.6} $$

我们将其分解为两个非对易部分

$$ H = H_{0} + H_{1}\tag{5.7} $$

定义为

$$ H_{0} = - J^{z} \sum_{i = 1}^{N} S_{i}^{z} S_{i + 1}^{z} , \quad H_{1} = - H^{x} \sum_{i = 1}^{N} S_{i}^{x} .\tag{5.8} $$

² 在经典模型(1.3)中，自旋变量是 $\sigma_i$，取其值为 $\pm 1$。

构建有限温度量子蒙特卡洛算法的难点在于对哈密顿量取指数。由于当前模型存在精确解 (Sachdev, 1999)，原则上我们可以精确地进行取指数运算，因此不需要蒙特卡洛模拟。然而，我们的目的是教学，所以暂时假设我们无法进行该取指数运算。

我们首先将配分函数 $Z = \mathrm{Tr} \ e^{- \beta \cal H}$ 改写为具有经典自由度的形式。注意到任何算符都与自身对易，我们写出

$$ e^{- \beta H} = \underbrace{e^{- \Delta \tau H} e^{- \Delta \tau H} \cdot \cdot \cdot e^{- \Delta \tau H}}_{M \mathrm{个因子}} ,\tag{5.9} $$

其中 $\Delta \tau = \beta / M$，因此

$$ {\cal Z} = \sum_{C_{1} , C_{2} , \ldots , C_{M}} \langle C_{1} | e^{- \Delta \tau H} | C_{M} \rangle \langle C_{M} | e^{- \Delta \tau H} | C_{M - 1} \rangle \cdot \cdot \cdot \langle C_{2} | e^{- \Delta \tau H} | C_{1} \rangle .\tag{5.10} $$

该方程是量子力学配分函数的费曼路径积分（Feynman path-integral）表达式。哈密顿量指数的矩阵乘积遍历一个有序的基矢序列。每个基态对应相空间中的一个构型，因此我们可以将 (5.10) 中的求和理解为对虚时（imaginary-time）有序构型序列的求和。构型序列 $C_{1} \ \ C_{2} \ \ \cdot \cdot \cdot C_{M} \ \ C_{1}$ 定义了相空间中的一条路径。因此，配分函数是对所有可能（闭合）路径的求和。在采样构型时，蒙特卡洛算法会选择出这些路径中最重要的一些。

目前来看，(5.10) 是一个平凡的结果。如果我们知道 $\exp ( - \Delta \tau H )$ 在所选基矢下的矩阵元，就可以轻松计算它。如果知道这些矩阵元，路径积分表述反而是一种退步，因为我们必须计算所有矩阵乘积，这比直接求迹 $\exp ( - \beta H )$ 的矩阵元工作量更大。表述 (5.10) 适用于我们无法精确表达 $\exp ( - \Delta \tau H )$ 矩阵元的情况。通过将 $\exp ( - \beta H )$ 写成多个 $\exp ( - \Delta \tau H )$ 因子的乘积，我们将一个小参数 $\Delta \tau$ 引入了问题，这使得我们可以为 $\exp ( - \Delta \tau H )$ 的每个矩阵元获得一个好的近似。

对 $H$ 取指数的困难源于 $S_{i}^{x}$ 和 $S_{j}^{z}$ 在同格点 $( i = j )$ 上不具有对易性（noncommutativity）。因此，$H$ 的两个部分不能同时对角化。相应地，(5.5) 中定义的态不是哈密顿量 (5.6) 的本征态；它们仅是 $H_{0}$ 的本征态。然而，在这个基矢下（或者在 $H_{1}$ 的本征基下），$\langle C^{\prime} | e^{- \Delta \tau H_{0}} | C \rangle$ 和 $\langle C^{\prime} | e^{- \Delta \tau H_{1}} | C \rangle$ 都是已知的。不幸的是，如果 $[ H_{0} , H_{1} ] ~ \neq ~ 0$，则 $e^{- \Delta \tau ( H_{0} + H_{1} )} \ \ne$ $e^{- \Delta \tau H_{0}} e^{- \Delta \tau H_{1}}$ 。为了继续，我们需要制定一种策略来处理 $[ A , B ] \neq 0$ 情况下形如 $e^{- \Delta \tau ( A + B )}$ 的指数。

让我们考虑如果我们采用近似

$$ e^{- \Delta \tau ( A + B )} \approx e^{- \Delta \tau A} e^{- \Delta \tau B} \approx e^{- \Delta \tau B} e^{- \Delta \tau A} .\tag{5.11} $$

所引入的误差。对于算符（和矩阵），$\exp ( A ) = I + A + {\textstyle{\frac{1} {2 !}}} A^{2} + {\textstyle{\frac{1} {3 !}}} A^{3} + \cdot \cdot \cdot$ 由其幂级数定义。因此，

$$ e^{- \Delta \tau ( A + B )} = e^{- \Delta \tau A} e^{- \Delta \tau B} + \mathcal{O} ( ( \Delta \tau )^{2} ) ,\tag{5.12} $$

这意味着误差量级为 $( \Delta \tau )^{2}$ 。一个更高阶（对称）的近似是

$$ e^{- \Delta \tau ( A + B )} = e^{- \Delta \tau B / 2} e^{- \Delta \tau A} e^{- \Delta \tau B / 2} + \mathcal{O} ( ( \Delta \tau )^{3} ) .\tag{5.13} $$

像(5.12)和(5.13)这样的近似称为铃木-Trotter近似（Suzuki-Trotter approximations, Hatano and Suzuki, 2005），或简称为Trotter近似。

自然地，我们选择一组基，在其中可以方便地求出$H_{0}$和$H_{1}$指数算符的矩阵元。基矢(5.5)使得对$H_{0}$求指数算符变得简单，同时也便于计算$\exp ( - \Delta \tau H_{1} )$的矩阵元。利用低阶Trotter近似(5.12)，我们写出$e^{- \Delta \tau H} \approx e^{- \Delta \tau H_{0}} e^{- \Delta \tau H_{1}}$，然后有

$$ \langle{\cal C}_{k + 1} | e^{- \Delta \tau H} | C_{k} \rangle \approx e^{- \Delta \tau E_{0} \left( C_{k + 1} \right)} \langle C_{k + 1} | e^{- \Delta \tau H_{1}} | C_{k} \rangle , $$

其中$E_{0} ( C_{k} )$是零磁场下对应态$| C_{k} \rangle$的经典伊辛能量（同样，自旋$s_{i}$取值为$\pm \frac{1} {2}$而非±1）。

我们必须为每个时间步$\tau_{k} = k \Delta \tau$计算能量$E_{0} ( C_{k} )$，并通过为自旋本征值附加额外下标来记录每个自旋本征值对应哪个时间步，从而写出

$$ E_{0} ( C_{k} ) = - J^{z} \sum_{i} s_{i , k} s_{i + 1 , k} . $$

因此，

$$ e^{- \Delta \tau \sum_{k} E_{0} ( C_{k} )} = e^{\Delta \tau J^{z} \sum_{i k} s_{i , k} s_{i + 1 , k}} , $$

这对应于M条解耦的经典伊辛类链的玻尔兹曼权重。

剩下的任务是确定哈密顿量中横场部分的指数算符的矩阵元。首先，我们注意到

$$ \left. C | e^{- \Delta \tau H_{1}} | C^{\prime} \right. = \prod_{i = 1}^{N} \left. s_{i} | e^{\Delta \tau H^{x} S_{i}^{x}} | s_{i}^{\prime} \right. .\tag{5.14} $$

接下来，

$$ \begin{array} {r l} {e^{\Delta \tau H^{x} S_{i}^{x}}} & {= I + \Delta \tau S_{i}^{x} H^{x} + \frac{1} {2 !} \big ( \Delta \tau H^{x} S_{i}^{x} \big )^{2} + \frac{1} {3 !} \big ( \Delta \tau H^{x} S_{i}^{x} \big )^{3} + \cdots} \\ & {= I + ( \Delta \tau H^{x} ) S_{i}^{x} + \frac{1} {4 \cdot 2 !} ( \Delta \tau H^{x} )^{2} I + \frac{1} {4 \cdot 3 !} ( \Delta \tau H^{x} )^{3} S_{i}^{x} + \cdots ,} \end{array} $$

这里我们利用了$\begin{array} {r} {( S^{x} )^{2} = \frac{1} {4}} \end{array}$这一事实。累加幂级数得到

$$ \begin{array} {r} {\langle s | e^{\Delta \tau S^{x} H^{x}} | s^{\prime} \rangle = \langle s | s^{\prime} \rangle \cosh{\frac{1} {2}} \Delta \tau H^{x} + 2 \langle s | S^{x} | s^{\prime} \rangle \sinh{\frac{1} {2}} \Delta \tau H^{x} .} \end{array} $$

注意到第二项中唯一的非零矩阵元为

$$ \begin{array} {r} {\langle \downarrow | S^{x} | \uparrow \rangle = \frac{1} {2} , \quad \langle \uparrow | S^{x} | \downarrow \rangle = \frac{1} {2} ,} \end{array} $$

我们将指数算符的矩阵元重新表示为

$$ \langle s \vert e^{\Delta \tau S^{x} H^{x}} \vert s^{\prime} \rangle = Q e^{K^{\tau} s s^{\prime}} ,\tag{5.15} $$

其中$s , s^{\prime} = \pm \frac{1} {2}$，且

$$ \begin{array} {r} {Q = \sqrt{\sinh \frac{1} {2} \Delta \tau H^{x} \cosh \frac{1} {2} \Delta \tau H^{x}} , \quad K^{\tau} = 2 \ln \coth \frac{1} {2} \Delta \tau H^{x} .} \end{array}\tag{5.16} $$

最终，

$$ Z \propto \sum_{C_{1} , C_{2} , \ldots , C_{M}} \exp \left( \Delta \tau J^{z} \sum_{i , k} s_{i , k} s_{i + 1 , k} + K^{\tau} \sum_{i , k} s_{i , k} s_{i , k + 1} \right) .\tag{5.17} $$

这正是一个零场、各向异性、二维经典伊辛类模型的配分函数（除了因子$\frac{1} {2}$在场中的差异）。如果我们将这个经典模型的逆温度$\beta_{\mathrm{classical}}$设为1，则x方向上的耦合为$K^{x} = \Delta \tau J^{z}$，τ方向上的场依赖交换为$K^{\tau}$：

$$ H_{\mathrm{classical}} = - K^{x} \sum_{i , k} s_{i , k} s_{i + 1 , k} - K^{\tau} \sum_{i , k} s_{i , k} s_{i , k + 1} .\tag{5.18} $$

该格点在空间维度上的大小为$L$，在时间维度上的大小为$M = \beta / \Delta \tau$。在$\Delta \tau$的主阶上，空间和时间方向的耦合常数为

$$ \begin{array} {r} {K^{x} = \Delta \tau J^{z} , \qquad K^{\tau} = - 2 \ln \frac{1} {2} \Delta \tau H^{x} .} \end{array}\tag{5.19} $$

见（5.16）。

5.2 量子自旋#

让我们总结一下刚刚发生的事情。通过将 $\exp ( - \beta H )$ 的矩阵元表示为 $\exp ( - \Delta \tau H )$ 的矩阵乘积，利用 Trotter 近似，并以 $H_{0}$ 的本征态为基，我们得到了一个由伊辛类变量 $s_{i , k}$ 构成的二维格点。表示 $\exp ( - \Delta \tau H_{0} )$ 的矩阵产生了无外场、独立伊辛类链的类玻尔兹曼因子。这些链中的自旋变量仅在 x 方向相互作用。表示 $\exp ( - \Delta \tau H_{1} )$ 的矩阵则在 τ 方向产生了这些链之间的相互作用。由路径积分表示可知，这个 τ 方向实际上是一个虚时方向。将 $\exp ( - \beta H )$ 表示为因子 $\exp ( - \Delta \tau H )$ 的乘积，得到了路径积分的离散虚时表示。根据费曼的说法（Feynman and Hibbs, 1965），将场论转化为量子统计力学是通过对路径积分的动力学因子 $\exp ( - i t H )$ 进行解析延拓，令 $i t \to \tau$，并要求在虚时方向上具有周期性来实现的。迹表达了这种周期性。在将 d 维量子系统映射到 (d+1) 维经典模型时，哈密顿量并非简单地提升到额外维度，而是沿着这个额外维度引入了新的相互作用，以描述虚时动力学。

在我们基于 Trotter 分解的推导中，策略是将哈密顿量分解成若干部分，使其指数形式能在某个基下（通常是哈密顿量某一部分的本征基）简便地求得。最直接的后果是引入了离散化误差，然而，如下文所述，我们可以通过重新表述算法来完全避免这一误差。另一个更棘手的问题（即使在连续时间表述中依然存在）是可能出现符号问题。由于哈密顿量是厄密的，其指数总是一个正定算子，这意味着它的所有特征值都大于零。算子的这一性质与我们所选择的用以表示其矩阵元的基无关。例如，在纵场伊辛模型中，表示哈密顿量指数的矩阵就明显地体现了这种正定性。它在所选的基下是对角化的，并且因为其对角元是某个实数的指数，所以它们都是正的。

正定矩阵不同于正矩阵（即所有元素都大于零的矩阵），也不同于非负矩阵（即没有负元素的矩阵）。正定性和非负性取决于我们选择用来表示算子的基。对于蒙特卡洛模拟，我们希望得到正矩阵或非负矩阵。然而，量子虚时传播子仅能保证是正定的。只有在对角化表示中（而我们几乎从不使用这种表示），才能确保表示 $\exp ( - \beta H )$ 的矩阵是非负的。如果我们无法找到一种可行的表示来消除负的矩阵元，就可能出现符号问题。这个问题会影响我们验证结果的能力，并且可能严重到使量子蒙特卡洛算法无法为计算平均值相关的统计误差提供可靠的估计。我们将在后续章节中进一步讨论这个问题。

对于横场模型，可能出现了符号问题，但先前未加评论地跳过了。在 (5.15) 中，如果 $H^{x}$ 是负的，非对角矩阵元就是负的。如果是这样，某些构型会对配分函数做出负贡献。针对当前情况，我们可以通过要求（正如我们隐含假设的那样）$H^{x} > 0$ 来避免符号问题。对于 $H^{x} < 0$ 的情况，我们可以通过对哈密顿量应用一个幺正变换 $( S_{i}^{x} , S_{i}^{y} , S_{i}^{z} ) \to ( - S_{i}^{x} , S_{i}^{y} , - S_{i}^{z} )$ 来恢复物理，这会导致磁场 $H^{x} \to - H^{x}$ 的反转。这种不变性体现在：对于每个自旋，在虚时方向上由 $S_{i}^{x}$ 引起的自旋翻转次数总是偶数，这使我们能够简单地忽略与单个自旋翻转事件相关的符号。

5.2.3 连续时间极限#

将一维量子横向场伊辛模型转化为（1+1）维经典伊辛模型，意味着我们可以将任何用于经典伊辛模型的蒙特卡洛算法作为量子蒙特卡洛算法，用于横向场模型。然而，除非我们取极限 $M \to \infty ( \Delta \tau \to 0 )$ ，否则 (5.17) 并不精确，因此我们的模拟结果会因时间离散化而存在系统误差。我们需要消除这种误差。一个直接的解决方案是使用不同的 $M$ 值进行多次蒙特卡洛模拟，然后外推至 $M \to \infty$ 。尽管此方法行之有效，并在实践中得到应用，但对于大多数问题来说，它已不再是必需的。在本节和第5.2.5节中，我们将第4.4节讨论的团簇算法的基本思想扩展到两个简单的量子自旋模型。然后，在第6章中，我们将讨论其在各式其他系统中的应用。

现在我们解释如何对横向场伊辛模型 (5.6) 取连续时间极限。前一小节的讨论表明，该模型的离散虚时表示不过是具有各向异性相互作用的二维经典伊辛模型，相互作用强度为 $\begin{array} {r l} {K^{x}} & {{} =} \end{array}$ $\Delta \tau J^{x}$ 和 $K^{\tau} = - 2 \ln ( \Delta \tau H^{x} / 2 )$ ，且“温度”为 $1 / \beta_{\mathrm{classical}} ~ = ~ 1$ 。我们观察到，随着 $\Delta \tau$ 的减小，时间方向的有效耦合常数会增加。这种增加是量子到经典映射的一个普遍特征，它确保了构型（在此情形下，是自旋向上和向下时间间隔的集合）具有明确定义的极限 $\Delta \tau \to 0$ 。我们可以通过注意到离散时间表述中的自旋构型仅仅是“真实”连续费曼路径的“低分辨率”图像来理解这一点。假设在连续极限下，自旋翻转（扭结）之间的平均距离为 $\tau_{\mathrm{kink}}$ 。那么，在离散时间中，我们平均每 $\tau_{\mathrm{kink}} / \Delta \tau$ 个时间步长就有一个扭结。为了在时间方向获得如此长的特征长度（以 $\Delta \tau$ 为单位），该方向的有效耦合必须很大。

我们已经引入了Swendsen-Wang算法作为模拟伊辛模型的高效算法（第4.4节）。为了将此方法应用于横向场伊辛模型的离散时间表示 (5.17)，我们切换到伊辛变量 $\tilde{s}_{i} = \pm 1$ ，并采用耦合常数 $\begin{array} {r} {\tilde{K}^{x} = \frac{1} {4} K^{x} = \frac{1} {4} \Delta \tau J^{z}} \end{array}$ 和 $\begin{array} {r} {\tilde{K}^{\tau} = \frac{1} {4} K^{\tau} = - \frac{1} {2} \ln ( \frac{1} {2} \Delta \tau H^{x} )} \end{array}$ 。在离散时间团簇模拟中，我们随后

以概率 $\begin{array} {r} {1 - e^{- 2 \tilde{K}^{x}} \approx \frac{1} {2} \Delta \tau J^{z}} \end{array}$ 放置一个水平键，该键绑定空间方向上的两个相邻平行自旋。
以概率 $\begin{array} {r} {1 - e^{- 2 \tilde{K}^{\tau}} \approx 1 - \frac{1} {2} \Delta \tau H^{x}} \end{array}$ 连接时间方向上的相邻平行自旋。

我们可以用不涉及 $\Delta \tau$ 但在极限 $\Delta \tau \to 0$ 下统计等价的不同过程来替换这些过程。

为此，我们改变构型的“编码系统”。在离散化虚时表示中，我们在每个时空格点上定义了一个自旋。当虚时连续时，这显然不可能做到，因此，我们不再指定每个时空点上的自旋值，而是为每个格点只指定扭结（局部自旋值发生变化的点）的虚时位置以及这些扭结之后的自旋值。因为在给定离散时间点出现扭结的概率与 $1 / M$ 成正比，而时间点的数量是 M，所以扭结的总数不应强烈依赖于 M，并且应该在极限 $M \infty$ 下收敛到某个有限值。因此，我们可以用新的马尔可夫过程中的一个构型 C 表示为：

$$ \begin{array} {c} {C = ( ( \tau_{1} ( 1 ) , \tilde{s}_{1} ( 1 ) ) , ( \tau_{1} ( 2 ) , \tilde{s}_{1} ( 2 ) ) , \cdots , ( \tau_{1} ( n_{1} ) , \tilde{s}_{1} ( n_{1} ) ) ,} \\ {( \tau_{2} ( 1 ) , \tilde{s}_{2} ( 1 ) ) , \cdots ( \tau_{N} ( n_{N} ) , \tilde{s}_{N} ( n_{N} ) ) ) ,} \end{array} $$

其中 $\tau_{i} ( k )$ 是格点 i 上第 k 个扭结（即第 k 次自旋翻转）的虚时值，$\tilde{s}_{i} ( k )$ 是该扭结之后伊辛自旋的值，$n_{i}$ 是格点 i 上的扭结数。（若格点 i 上没有扭结，则只存储自旋态。）

Swendsen-Wang 演算法的第二步相当于在沿时间方向上的相邻对之间以概率 $\begin{array} {r l r} {\bar{e}^{- 2 \tilde{K}^{\tau}} =} \end{array}$ $\frac{1} {2} \Delta \tau H^{x}$ 放置切割（cut）。在 $\Delta \tau \to 0$ 极限下，这相当于通过按均匀概率密度 $\scriptstyle{\frac{1} {2}} H^{x}$ 在虚时轴上随机分布的切割，将恒定自旋的片段切成更短的片段（图 5.1）。

图5.1 横向伊辛模型连续时间团簇演算法示意图。竖直线代表四个格点的虚时间隔。自旋向上片段用粗线标记。四个步骤为：(i) 在虚时间隔上插入切割，(ii) 插入最近邻键（只有当两侧自旋态相同时才可能插入键，用实水平线表示；被拒绝的键用虚线标记），(iii) 随机翻转由最近邻键连接的团簇，以及 (iv) 移除切割并定义新片段。在第二和第三面板中，我们用灰色阴影突出显示其中一个团簇。

另一方面，在空间方向上，离散时间 Swendsen-Wang 演算法的第一步是以概率 $\begin{array} {r} {1 - e^{- 2 \tilde{K}^{x}} \sim \frac{\bar{1}} {2} \Delta \tau J^{z}} \end{array}$ 在两个相邻平行自旋之间放置键。在连续时间极限下，这些键以密度 $\scriptstyle{\frac{1} {2}} J^{z}$ 绑定两个重叠的相邻片段——前提是它们具有相同的自旋值。实际操作中，我们以均匀概率密度 $\scriptstyle{\frac{1} {2}} J^{z}$ 提出潜在的键插入时间（不论自旋态），然后检查由这些键连接的两个片段上的自旋态是否相同。仅当键所在位置的自旋态相同时，才插入该键（图 5.1 中的实水平线）；否则，拒绝提出的键（图 5.1 中的虚线）。

Swendsen-Wang 演算法的下一步是识别并随机翻转由空间键连接在一起的片段团簇。为此，我们使用并查集演算法，例如第 13.4 节讨论的 Hoshen-Kopelman 演算法（Hoshen and

演算法 11 横向伊辛模型的连续虚时团簇演算法。 输入: 由自旋向上和自旋向下片段组成的构型。 对每个格点 do 移除所有自旋态不改变的已有切割；通过均匀随机地从区间 $[0, \beta)$ 中以密度 $\scriptstyle{\frac{1} {2}} H^{x}$ 选取切割时间 $\tau$ 生成新切割； 结束对 对每对相互作用的格点 do 以密度 $\scriptstyle{\frac{1} {2}} J^{z}$ 均匀随机地从区间 $[0, \beta)$ 选取潜在的键时间 $\tau$； 对每个选中的 τ do 如果两个格点在 τ 时刻具有相同的自旋取向，则在 τ 时刻插入一条连接两个片段的水平键； 结束如果 结束对 结束对 识别连接的片段团簇； 对每个团簇 do 以概率 ${\frac{1} {2}}$ 翻转团簇中所有片段的自旋； 结束对 返回更新后的构型。

Kopelman，1976）。图 5.1 中灰色阴影标出了一个团簇实例。最后，我们通过移除过时的切割来定义新片段。演算法 11 给出了这些步骤的伪代码（Rieger and Kawashima，1999）。

由于我们只记录片段端点的位置以及与每个片段相关联的自旋状态，该算法（在低温且 $H^{x} \neq 0$ 条件下）所需的内存与切割数量成正比，即 $\mathcal{O} ( \beta H^{x} N )$。因此，与离散时间表示相比，其内存需求至少减少了 $M / ( \beta H^{x} )$ 倍。尽管只要该因子大于 1，离散时间模拟就能产生定性地正确的结果，但我们通常需要使其远大于 1，以减少系统误差并获得定量的正确结果。

现在讨论如何在算法 11 中生成均匀且随机分布的时间。需注意这些时间及其数量 n 是随机变量，并且这些时间最终要按升序排列，即 $0 < \tau_{1} < \tau_{2} < \cdots < \tau_{n} < \beta$。一种生成具有均匀概率密度 $\lambda$ 的离散时间有序事件的随机过程称为齐次泊松过程（Karlin and Taylor，1975）。毫不意外地，我们将需要用到泊松分布 (2.11)，

$$ P ( N ( \beta ) = n ) = \frac{( \lambda \beta )^{n}} {n !} e^{- \lambda \beta} ,\tag{5.20} $$

其中 $N ( \beta )$ 是在区间 $[ 0 , \beta )$ 内发生的事件数量。

我们可以通过以下方式说明该分布的适用性：如第 2.3.1 节所述，区间 $[ 0 , \beta )$ 内发生的平均事件数为 $\lambda \beta$，因此它们以均匀密度 $\lambda$ 发生。从其函数形式很容易看出，泊松分布具有以下短时行为：

当区间长度 $\Delta \tau$ 很小时，在时间段 $[ \tau , \tau + \Delta \tau )$ 内恰好发生一个事件的概率为 $\lambda \Delta \tau$，即它与区间长度成正比。
在此区间内没有事件发生的概率为 $1 - \lambda \Delta \tau$。
两个事件发生的概率为零，直到 $\mathcal{O} ( ( \Delta \tau )^{2} )$ 量级。

重要的观察结果是，这种短时行为与我们基于极限情况的概率一致。关于在时间 $\Delta \tau$ 内可发生事件数量的条件与我们隐含假设的相同。因此，为了开发对时间进行采样的算法，我们无需从短时性质逐步构建到有限区间，而可以直接在目标有限时间区间 $[ 0 , \beta )$ 上使用泊松分布。

泊松分布的函数形式引出了许多有趣的性质，为采样算法提供了严格基础。根据定义可直接得到 $P ( N ( \beta ) = 1 ) = \lambda \beta \exp ( - \lambda \beta )$（指数分布）和 $P ( N ( \beta ) = 0 ) = \exp ( - \lambda \beta )$。以下重要关系式（不加证明给出，Karlin and Taylor，1975）：

$$ P ( N ( \tau_{1} + \tau_{2} ) - N ( \tau_{2} ) = n ) = P ( N ( \tau_{1} ) = n ) = \frac{( \lambda \tau_{1} )^{n}} {n !} e^{- \lambda \tau_{1}} , \quad \tau_{1} , \tau_{2} > 0 . $$

这一关系表明，我们可以将一个泊松过程分割成多个区间，并且每个区间内的分布也是一个泊松分布。该结果隐含的事实是，不重叠的区间在统计上是独立的。

现在我们利用上述性质来生成算法 11 所需的时间序列。我们一次生成一个事件。在任何有限时间区间内发生一个事件的概率是一个指数概率，其速率常数为 λ。在第 2.3.2 节中，我们展示了如何通过生成一个均匀随机数 $\zeta$ 并计算下式来从该连续概率密度中进行采样：

$$ \tau_{1} = - \ln \zeta / \lambda $$

算法12 以概率密度 $\lambda$ 顺序生成事件 $\tau_{k} < \tau_{k + 1} \in [ 0 , \beta )$

输入: $\lambda , \beta .$ $\tau = 0$ $k = 0$ ; 重复生成一个均匀随机数 $\zeta \in ( 0 , 1 ]$ $\tau \leftarrow \tau - \ln \zeta / \lambda$ ; $k \gets k + 1$ ; $\tau_{k} \gets \tau :$ 直到 $\tau \geq \beta .$ $n = k - 1$ ; 返回 $n$ 和 $\{\tau_{1} , \tau_{2} , \ldots , \tau_{n} \}$

以获得第一个事件时间 $\tau_{1}$ 。利用泊松时间区间的统计独立性，我们简单地重复该过程以获得 $\tau_{2} - \tau_{1}$ ，即从当前事件时间 $\tau_{1}$ 到下一个事件的时间。然后，通过相同的过程，我们得到 $\tau_{3} - \tau_{2}$ ，依此类推，直到 $\sum_{i} \tau_{i}$ 大于 $\beta$ 。此策略总结在算法12中。

虽然该算法是有效的，但存在一种通常更高效的替代方案，因为它省去了对数的计算。该替代方案也是一个简单的算法，但其合理性并非显而易见。为了阐述这个想法，我们首先考虑区间 $[ 0 , \beta )$ 中仅发生一个事件的情形，并探讨其位置是如何分布的。我们分两步回答这个问题。首先，考虑一个速率为 $\lambda$ 的泊松过程，并以其 $N ( \beta ) = 1$ 为条件，即 $[ 0 , \beta )$ 中仅发生一个事件。然后，考察到达时间 $\tau_{1}$ 相对于区间 $[ 0 , \beta )$ 内均匀分布的任意其他随机时间 $\xi$ 的关系。根据条件概率的定义以及非重叠区间内事件的独立性：

$$ \begin{array} {r l} & {P \left( \tau_{1} \le \xi \left| N ( \beta ) = 1 \right. \right) = \displaystyle \frac{P ( \tau_{1} \le \xi , N ( \beta ) = 1 )} {P ( N ( \beta ) = 1 )}} \\ & {\quad \quad \quad = \displaystyle \frac{P ( N ( \xi ) = 1 , N ( \beta ) - N ( \xi ) = 0 )} {P ( N ( \beta ) = 1 )}} \\ & {\quad \quad \quad = \displaystyle \frac{\lambda \xi e^{- \lambda \xi} e^{- \lambda ( \beta - \xi )}} {\lambda \beta e^{- \lambda \beta}}} \\ & {\quad \quad \quad = \displaystyle \frac{\xi} {\beta} .} \end{array} $$

该结果表明，区间 $[ 0 , \beta )$ 内单一事件的到达时间在该区间上均匀分布。泊松点位于 $\xi$

如果我们在 $[ 0 , \beta )$ 上均匀生成一组 $n$ 个随机数 $\{\xi \}$ ，我们将生成 $n$ 个时间 $\{t_{\xi_{n}} \}$ ，每个时间都均匀分布在区间 $[ 0 , \beta )$ 上。在这些时间的 $n!$ 种排列中，只有一种是升序排列。通过排序，这组时间总可以排列成此顺序。有序序列 $0 < t_{1} < t_{2} < \cdot \cdot \cdot < t_{n} < \beta$ 的联合分布（概率密度）为：

$$ P ( \tau_{1} = t_{1} , \dots , \tau = t_{n} | N ( \beta ) = n ) = \frac{n !} {\beta^{n}} . $$

如果我们按照上述方式生成所需的序列，这个分布是否与来自泊松分布的 $n$ 个事件的有序样本一致？

我们可以再次利用条件概率和泊松分布的统计独立性来证明这一点：

$$ \begin{array} {r l} & {P ( \tau_{1} = t_{1} , \dots , \tau_{n} = t_{n} | N ( \beta ) = n ) = \frac{P ( t_{1} , t_{2} , \dots , t_{n} , N ( \beta ) = n )} {P ( N ( \beta ) = n )}} \\ & {\qquad = \frac{P ( t_{1} , t_{2} - t_{1} , \dots , t_{n} - t_{n - 1} , t_{n + 1} > \beta - t_{n} )} {P ( N ( \beta ) = n )}} \\ & {\qquad = \frac{\lambda^{n} e^{- \lambda \beta}} {P ( N ( \beta ) = n )}} \\ & {\qquad = \frac{n !} {\beta^{n}} .} \end{array} $$

从第一式到第二式，我们利用了时间和时间间隔的泊松性质。从第二式到第三式，我们注意到在每个间隔内分布是泊松的，并且除了最后一个间隔外，每个间隔内只发生一个事件。通过泊松分布的定义，我们从第三式得到第四式。算法 13 表述了这一结果。给定 $\beta$ 和 $\lambda$，第一部分计算 $N(\beta)=n$，第二部分生成无序集合 $\{\tau_i\}$，最后一步通过对该集合中的值排序产生所需结果。

5.2.5 零场 XY 模型#

在本节及随后的三节中，我们以零场一维自旋 $-\frac12$ 量子 XY 模型为例讨论各种基本概念：

$$ H = - \sum_{i = 1}^{N} J^{x y} \big ( S_{i}^{x} S_{i + 1}^{x} + S_{i}^{y} S_{i + 1}^{y} \big )\tag{5.21} $$

我们针对该模型讨论的算法，即环回算法（Section 5.2.5）和蠕虫算法（Section 5.2.6），代表了具有广泛应用的量子蒙特卡洛算法的重要类别，如我们将在第 6 章所见。

算法 13 生成密度为 $\lambda$ 的均匀分布随机事件 $\tau_{k} \in [ 0 , \beta )$ 的更快方法输入：$\lambda , \beta$ $n = 0$; $d = e^{- \beta \lambda}$ $p = d$: 生成均匀随机数 $\zeta \in [ 0 , 1 ]$ 当 $\zeta > p$ 时执行 $n \gets n + 1$; $d \gets d \beta \lambda / n$; $p \gets p + d$ 结束循环对于 $k = 1$ 到 n 生成均匀随机数 $\tau_{k} \in [ 0 , \beta )$ 结束循环将集合 $\{\tau_{k} \}$ 按升序排序; 返回 $\{\tau_{k} \}$

再次，我们首先将配分函数用经典自由度表示。虽然该模型中没有任何项在 $S^{z}$ 基（5.5）中对角化，但我们仍因其方便性而继续使用该基。注意到不同格点上的自旋算符对易，因此所有起始于奇数（偶数）格点的自旋算符对都对易，我们将 H 分成两部分：

$$ H = H_{\mathrm{odd}} + H_{\mathrm{even}} , $$

其中

$$ \begin{array} {l} {{\displaystyle{\cal H}_{\mathrm{odd}} = - \sum_{i = 1}^{N / 2} J^{x y} \big ( S_{2 i - 1}^{x} S_{2 i}^{x} + S_{2 i - 1}^{y} S_{2 i}^{y} \big )} ,} \\ {{\displaystyle ~}} \\ {{\cal H}_{\mathrm{even}} = - \sum_{i = 1}^{N / 2} J^{x y} \big ( S_{2 i}^{x} S_{2 i + 1}^{x} + S_{2 i}^{y} S_{2 i + 1}^{y} \big ) .} \end{array} $$

对配分函数应用 Trotter 近似，我们得到

$$ Z = \mathrm{Tr} \Bigl [ \left( e^{- \Delta \tau H_{\mathrm{odd}}} e^{- \Delta \tau H_{\mathrm{even}}} \right)^{M} \Bigr ] , $$

然后在每个因子之间插入完备基，我们发现需要计算的基本矩阵元为

$$ \begin{array} {r} {w \left( \begin{array} {c c} {s_{i , k + 1}} & {s_{i + 1 , k + 1}} \\ {s_{i , k}} & {s_{i + 1 , k}} \end{array} \right) \equiv \big \langle s_{i , k + 1} s_{i + 1 , k + 1} \big \vert e^{\Delta \tau J^{\mathrm{xy}} \big ( S_{i}^{x} S_{i + 1}^{x} + S_{i}^{y} S_{i + 1}^{y} \big )} \big \vert s_{i , k} s_{i + 1 , k} \big \rangle .} \end{array} $$

计算它相对直接。首先，对于 $i \neq j$，且利用 $S^{+} = ( S_{x} + i S_{y} )$ 和 $S^{-} = ( S_{x} - i S_{y} )$，

$$ \begin{array} {r} {S_{i}^{x} S_{j}^{x} + S_{i}^{y} S_{j}^{y} = \frac{1} {2} \left( S_{i}^{+} S_{j}^{-} + S_{i}^{-} S_{j}^{+} \right) .} \end{array} $$

将其代入指数函数并将指数展开为幂级数后，我们逐项显式分析可能的矩阵元，得到

$$ \begin{array} {r l} {\left( \begin{array} {c c c c c} {w \left( \frac{\uparrow \uparrow} {\uparrow \uparrow} \right)} & {w \left( \frac{\uparrow \downarrow} {\uparrow \uparrow} \right)} & {w \left( \frac{\downarrow \uparrow} {\uparrow \uparrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\uparrow \uparrow} \right)} \\ {w \left( \frac{\uparrow \uparrow} {\uparrow \downarrow} \right)} & {w \left( \frac{\uparrow \downarrow} {\uparrow \downarrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\uparrow \downarrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\uparrow \downarrow} \right)} \\ {w \left( \frac{\uparrow \uparrow} {\uparrow \uparrow} \right)} & {w \left( \frac{\uparrow \downarrow} {\uparrow \downarrow} \right)} & {w \left( \frac{\downarrow \uparrow} {\downarrow \uparrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\downarrow \downarrow} \right)} \\ {w \left( \frac{\uparrow \uparrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\uparrow \downarrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\downarrow \uparrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\downarrow \downarrow} \right)} \\ {w \left( \frac{\uparrow \uparrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\uparrow \downarrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\downarrow \downarrow} \right)} & {w \left( \frac{\downarrow \downarrow} {\downarrow \downarrow} \right)} \end{array} \right)} & {{}} \\ {= \left( \begin{array} {c c c c} {1} & {0} & {0} & {0} \\ {0} & {\cosh \frac{1} {2} \Delta \tau J^{\mathrm{op}}} & {\sinh \frac{1} {2} \Delta \tau J^{\mathrm{op}}} & {0} \\ {0} & {\sinh \frac{1} {2} \Delta \tau J^{\mathrm{op}}} & {\cosh \frac{1} {2} \Delta \tau J^{\mathrm{op}}} & {0} \\ {0} & {0} & {0} & {1} \end{array} \right) .} \end{array}\tag{5.22} $$

其中发生了若干变化。哈密顿量的奇偶位点分解与Trotter分解在虚时方向上创建了一个由$2M$层组成的结构。从偶位点出发的自旋算子对作用于$k=1,3,5,\dots$层，而从奇位点出发的自旋算子对作用于$k=2,4,6,\ldots$层。两个相邻层对应于完整哈密顿量推进一个时间步$\Delta \tau$的传播。

局部转移矩阵仅依赖于一层上两个相邻自旋变量以及下一层上相同位点对的值。这四个自旋位于$(i,k)$、$(i+1,k)$、$(i+1,k+1)$和$(i,k+1)$位置，定义了一个小方块（plaquette）。我们注意到这些小方块自然地将时空分解为棋盘结构（图5.2），并且我们仅需对亮色或暗色小方块进行采样。

每个小方块关联16种可能自旋构型。然而，(5.22)表明其中仅有6种允许存在。非零矩阵元表达了一个守恒律：

$$ s_{i , k} + s_{i + 1 , k} = s_{i , k + 1} + s_{i + 1 , k + 1} ,\tag{5.23} $$

对于每个小方块上的总角动量z分量。因此，总角动量z分量 $S^{z} \equiv \textstyle \sum_{i} S_{i}^{z}$ 在每个虚时均守恒。换言之，空间平面上自旋向上的数目与虚时无关。这是对称性导致相空间缩减的一个实例。量子伊辛模型在自旋反转下不变，而XY模型在所有自旋绕z轴任意旋转下不变。这种不变性使 $S^{z}$ 守恒。该守恒律要求：若我们在蒙特卡洛模拟过程中改变一个空间平面上的自旋向上数目，则必须在所有其他平面上以相同方式改变。不变性与守恒律的一个独特推论是：模拟初始化时产生的 $S^{z}$ 值将全程保持不变，除非我们引入某种“全局翻转”程序。

图5.2 一维空间和虚时中的棋盘结构（左）。棋盘上的世界线（右）。

需要强调的是，(5.23)不仅意味着总 $S^{z}$ 守恒，还表达了一个局域守恒律，仅涉及某小方块中的四个自旋变量。这种相空间的限制在定义采样策略时带来了挑战。逐个格点并借助Metropolis算法提议局部自旋反转的做法与守恒律不一致。遍历性采样必须同时涉及多个自旋。

若我们用点标记有自旋向上的阴影小方块的角点，并遵循如下规则将这些点连成线：在全向上方块中，线沿左边缘和右边缘竖直延伸（因此线不相交），则得到一组在虚时方向上回环的连续路径（见图5.2右图）。每条这样的路径称为一条世界线（world line）。于是蒙特卡洛采样就变成了形变世界线的问题。有一种特别便利的方法可同时保持自旋守恒：我们依次访问每个未加阴影的小方块，若有一条线段沿某条边延伸而非另一条边，则尝试将其从一条边移动到另一条边。我们在图5.3中示意了这些移动，其中假设未加阴影的小方块位于

图5.3 基本的世界线更新。世界线可从实线移动到虚线，反之亦然。

$$ \left( \begin{array} {c c} {{s_{i , k + 1}}} & {{s_{i + 1 , k + 1}}} \\ {{s_{i , k}}} & {{s_{i + 1 , k}}} \end{array} \right) . $$

随后基于下式通过Metropolis判定接受或拒绝该移动：

$$ P_{\mathrm{accept}} \equiv \mathrm{min} [ 1 , R_{\mathrm{north}} R_{\mathrm{south}} R_{\mathrm{east}} R_{\mathrm{west}} ] , $$

其中

$$ \begin{array} {r l} {R_{\mathrm{north}} = \frac{w \left( \begin{array} {c c} {s_{i , k + 2}} & {s_{i + 1 , k + 2}} \\ {- s_{i , k + 1}} & {- s_{i + 1 , k + 1}} \end{array} \right)} {w \left( \begin{array} {c c} {s_{i , k + 2}} & {s_{i + 1 , k + 2}} \\ {s_{i , k + 1}} & {s_{i + 1 , k + 2}} \end{array} \right)} ,} & {R_{\mathrm{south}} = \frac{w \left( \begin{array} {c c} {- s_{i , k}} & {- s_{i + 1 , k}} \\ {s_{i , k - 1}} & {s_{i + 1 , k - 1}} \end{array} \right)} {w \left( \begin{array} {c c} {s_{i , k}} & {s_{i + 1 , k}} \\ {s_{i , k - 1}} & {s_{i + 1 , k - 1}} \end{array} \right)} ,} \\ {R_{\mathrm{east}} = \frac{w \left( \begin{array} {c c} {- s_{i + 1 , k + 1}} & {s_{i + 2 , k + 1}} \\ {- s_{i + 1 , k}} & {s_{i + 2 , k}} \end{array} \right)} {w \left( \begin{array} {c c} {s_{i + 1 , k + 1}} & {s_{i + 2 , k + 1}} \\ {s_{i + 1 , k + 1}} & {s_{i + 2 , k + 1}} \end{array} \right)} ,} & {R_{\mathrm{west}} = \frac{w \left( \begin{array} {c c} {s_{i - 1 , k + 1}} & {- s_{i , k + 1}} \\ {s_{i - 1 , k}} & {- s_{i , k}} \end{array} \right)} {w \left( \begin{array} {c c} {s_{i - 1 , k + 1}} & {s_{i , k + 1}} \\ {s_{i - 1 , k}} & {s_{i , k + 1}} \\ {s_{i , k - 1}} & {s_{i , k}} \end{array} \right)}} \end{array}\tag{5.24} $$

分别是未着色方形以北、南、东、西方向的小方格的局部权重在世界线构型中提议变化前后的比值。该提议是翻转未着色方形边界上的自旋，前提是两侧的自旋态不同（算法 14）。

我们从 (5.22) 中注意到，当 $\ J^{x y} \ < \ 0$ 时存在符号问题。如果我们的晶格是二分的，即格点分为 A 和 B 两个子晶格，使得只有 B 格点是 A 格点的最近邻，反之亦然，那么我们可以通过执行一个规范变换（canonical transformation），反转一个子晶格上的自旋来消除该问题（Marshall, 1955）。方形晶格是二分的，而三角晶格则不是。非二分晶格常常导致局部自旋构型受挫。对于阻挫系统，符号问题通常持续存在。

算法 14 一维自旋-$\frac{1}{2}$ XY 模型局部更新世界线蒙特卡洛方法。 (5.22) 定义了 (5.24) 中的权重 w。输入：一个自旋构型 $\left\{{{s}_{i , k}} \right\}$ 对于每个未着色方形 $( s_{i , k} , s_{i + 1 , k} , s_{i + 1 , k + 1} , s_{i , k + 1} )$: 如果 $s_{i , k} = s_{i , k + 1}$ 且 $s_{i + 1 , k} = s_{i + 1 , k + 1}$ 且 $s_{i , k} = - s_{i + 1 , k}$，则：通过 (5.24) 计算 $R_{\mathrm{north}} , R_{\mathrm{south}} , R_{\mathrm{east}} , R_{\mathrm{west}}$；生成一个均匀随机数 $\zeta \in [ 0 , 1 ]$；如果 $\zeta < R_{\mathrm{north}} R_{\mathrm{south}} R_{\mathrm{east}} R_{\mathrm{west}}$，则： $s_{i , k} \gets - s_{i , k}$ $s_{i + 1 , k} \gets - s_{i + 1 , k}$ $s_{i + 1 , k + 1} \gets - s_{i + 1 , k + 1}$ $s_{i , k + 1} \gets - s_{i , k + 1}$ 需要考虑边界条件。结束如果结束如果结束循环返回更新后的自旋构型。

与横场伊辛模型离散时间算法的情况类似，我们必须通过外推到 $M \to \infty$ 的极限来消除由 Trotter 近似引起的系统误差。幸运的是，我们可以再次将离散时间算法替换为连续时间算法。由于我们对 XY 模型配分函数的重构并未产生类似于经典模型的形式，因此无法直接应用 Swendsen-Wang 算法并由此继续。然而，如下一节所述，我们可以采用类似的方法。

5.2.5 环回模拟#

我们现在将团簇方案（最初在 4.4 节引入，并在 5.2.3 节针对横场伊辛模型重新表述）推广到量子系统配分函数的路径积分表示中。这里我们以 XY 模型为例。5.2.3 节讨论的横场伊辛模型的推广，通过将离散虚时的路径积分映射到 $( d + 1 )$ 维经典伊辛模型，是相当直接的。而本节对 XY 模型的推广则涉及更多内容，即根据相互作用项的性质引入各种类型的图元。本节我们将侧重于“如何操作”而非“为何如此”，并将所提议团簇类型算法的详细推导和证明推迟到第 6 章。

图 5.4 自旋-$\frac{1}{2}$ XY 模型环回算法的图元，以及一个蒙特卡罗步骤。

与横场伊辛模型的情况类似，我们只需通过重新解释随机分配，以特定密度将割点和键分布到虚时间间隔上，而非分配到单个方格上，即可得到该算法的连续时间版本。虽然我们侧重于研究一维自旋-$\frac{1}{2}$ XY 模型，但在许多情况下，将团簇更新方案推广到更高维度和其它类型的哈密顿量是直接的（第 6 章）。

第一步是用克罗内克δ函数（Kronecker deltas）表示局部虚时演化算符。我们从虚时演化算符的矩阵元开始，这些矩阵元本质上与 (5.22) 相同，但仅修正到 $\Delta \tau$ 的一阶³。为方便起见，我们向成对哈密顿量中添加一个常数 $- {\frac{1} {4}} J^{x y}$，然后将矩阵元表示为

$$ \begin{array} {r l} & {\langle s_{i}^{\prime} s_{j}^{\prime} | 1 - \Delta \tau H_{i j} | s_{i} s_{j} \rangle} \\ & {\qquad = \langle s_{i}^{\prime} s_{j}^{\prime} | s_{i} s_{j} \rangle + \frac{1} {2} \Delta \tau J^{x y} \langle s_{i}^{\prime} s_{j}^{\prime} | ( S_{i}^{+} S_{j}^{-} + S_{i}^{-} S_{j}^{+} ) + \frac{1} {2} | s_{i} s_{j} \rangle} \\ & {\qquad = \delta_{s_{i}^{\prime} , s_{i}} \delta_{s_{j}^{\prime} , s_{j}} + \frac{1} {4} \Delta \tau J^{x y} ( \delta_{s_{i}^{\prime} , s_{j}} \delta_{s_{j}^{\prime} , s_{i}} + \delta_{s_{i} , - s_{j}} \delta_{s_{i}^{\prime} , - s_{j}^{\prime}} ) .} \end{array}\tag{5.25} $$

这三个项对应三种不同图元，如图5.4中圆圈所示。为确认这种图元分配是合理的，让我们检查不同小片构型的贡献。第二个 $\delta_{s_{i}^{\prime} , s_{j}} \delta_{s_{j}^{\prime} , s_{i}}$ 项在小片构型

$$ {\left( \begin{array} {l l} {\uparrow} & {\uparrow} \\ {\uparrow} & {\uparrow} \end{array} \right)} , \quad{\left( \begin{array} {l l} {\downarrow} & {\downarrow} \\ {\downarrow} & {\downarrow} \end{array} \right)} , \quad{\left( \begin{array} {l l} {\uparrow} & {\downarrow} \\ {\downarrow} & {\uparrow} \end{array} \right)} , \quad{\left( \begin{array} {l l} {\downarrow} & {\uparrow} \\ {\uparrow} & {\downarrow} \end{array} \right)} , $$

中贡献 $\scriptstyle{\frac{1} {4}} \Delta \tau J^{x y}$，而第三个 $\delta_{s_{i} , - s_{j}} \delta_{s_{i}^{\prime} , - s_{j}^{\prime}}$ 项在构型

$$ \left( \begin{array} {c c} {{\uparrow}} & {{\downarrow}} \\ {{\uparrow}} & {{\downarrow}} \end{array} \right) , \quad \left( \begin{array} {c c} {{\downarrow}} & {{\uparrow}} \\ {{\downarrow}} & {{\uparrow}} \end{array} \right) , \quad \left( \begin{array} {c c} {{\uparrow}} & {{\downarrow}} \\ {{\downarrow}} & {{\uparrow}} \end{array} \right) , \quad \left( \begin{array} {c c} {{\downarrow}} & {{\uparrow}} \\ {{\uparrow}} & {{\downarrow}} \end{array} \right) . $$

中贡献 $\frac{1} {4} \Delta \tau J^{x y}$。

在贡献第二项的所有构型中，对角线上的自旋变量是平行的，而贡献第三项的构型中，同一行的自旋变量是反平行的。

表达式 (5.25) 引导我们得出所讨论的XY模型的连续时间团簇算法。具体来说，我们将第一项识别为“竖直”图元，第二项识别为“对角”图元，第三项识别为“水平”图元。我们首先考虑世界线没有扭结的位置。竖直图元对应于恒等算符。对于对角图元，因为当 $s_{i} = s_{j}$ 时其权重为 $\frac{1} {4} J^{x y} \Delta \tau$，如果两个自旋相互平行，我们在虚时间间隔上以密度 ${\scriptstyle{\frac{1} {4}}} J^{x y}$ 生成该图元。该图元对局部自旋构型施加限制，即位于对角顶点的两个自旋必须相互平行，因此我们像图5.4那样，用一对（不相连的）对角线来表示它。类似地，(5.25) 中的第三项对 $s_{i} = - s_{j}$ 给出权重 $\frac{1} {4} J^{x y} \Delta \tau$。因此，我们将相应的图元以密度 ${\frac{1} {4}} J^{x y}$ 插入到具有反平行自旋的时间间隔中。代表该项的克罗内克 delta 函数表明，我们在水平方向连接反平行自旋，同样如图5.4所示。在世界线有扭结的位置，只有对角或水平图元可以匹配。由于这两个图元的权重相等，我们在每个扭结处以相等的概率选择其中一个图元。

由于所有三个图元都将方格上的四个点成对连接，因此生成的团簇是环，如图5.4中间面板所示（因此得名环回算法）。团簇算法是否会变成环回算法取决于哈密顿量的性质。

与横向场伊辛模型的情况类似，取连续时间极限 $\Delta \tau \to 0$ 只是将点表示切换为线段表示，并以适当密度放置图元（扭结位置除外）的问题。算法15总结了连续时间版本的环回算法。

算法15 自旋-$\cdot \frac{1} {2}$ XY 模型的环回算法。输入：一个世界线构型；对于所有扭结 do 以概率 $\frac{1} {2}$ 选择一个对角或水平图元，并将其分配给该扭结；结束对于所有相邻位点对 do 以密度 ${\frac{1} {4}} J^{x y}$ 生成图元位置 $\tau \in [ 0 , \beta )$ 如果自旋平行（反平行），则放置对角（水平）图元；结束识别环；对于所有环 do 以概率 $\frac{1} {2}$ 翻转该环；结束返回更新后的世界线构型。

5.2.6 使用虫洞进行模拟#

我们不必将整个晶格分解为环，然后通过随机翻转这些环来计算新的世界线构型，而是可以通过创建虫洞并让其移动来更新世界线构型，就像第4.5节中讨论的伊辛模型的虫洞更新一样。这里的基本思想相同：通过引入不连续点来扩展构型空间，相应地推广权重，并构造满足与广义权重的细致平衡的马尔可夫转移矩阵。在一些简单的例子中，包括本节展示的 $S ~ = ~ {\textstyle{\frac{1} {2}}} ~ X Y$ 模型，得到的算法与前一小节讨论的环回算法非常接近，主要区别仅在于虫洞是一个接一个地顺序创建和翻转的，而不像环回算法那样同时构建所有环。因此，该算法也被称为有向环回算法。⁴ 然而，一般而言，虫洞算法和环回算法之间没有直接的对应关系。我们将在第6章详细研究虫洞和有向环回算法。下面的简单例子是一个“预览”，没有严格推导该过程。

该算法的第一步是定义潜在扭结位置的位置。但不是在特定位置上插入特定的图元，而是插入“顶点”。顶点类似于环回算法中的图元，但顶点上的图元类型是预先不确定的。它更像是一个散射中心，虫洞头部会在该处以非确定性的方式改变其轨迹。

具体而言，在移除上一周期的顶点后，我们在每个已有的扭结上放置一个顶点，并以密度 ${\frac{1} {4}} J^{x y}$ 在每对相邻格点之间放置额外的顶点。然后，我们将虫洞的头部和尾部放置在时空中的一个随机选择点，并随机选择头部运动的初始方向（时间上的向前或向后）。假设选择的方向是“向前”。在这种情况下，头部在时间上向前移动，直到碰到下一个顶点。⁵ 如果该顶点位于扭结上，则我们以概率 $\frac{1} {2}$ 为该顶点分配一个图元（对角图元或水平图元）。如果顶点不在扭结上，则图元由世界线构型决定：如果顶点位于两个平行自旋之间，则图元为对角图元；否则，为水平图元。无论哪种情况，虫洞头部都会在顶点处发生散射。在对角图元的情况下，散射意味着头部在相邻格点上继续向前移动；而水平图元则意味着头部反转方向，在相邻格点上向后移动。遵循这些规则，虫洞头部最终会回到其原始位置，并与尾部湮灭。此时，世界线构型通过一个环得到了更新，然后我们随机放置另一对虫洞头部和尾部，并重复该过程。这个循环如图5.5所示。

一个值得注意的微妙之处在于，分配给顶点的图元可能在虫洞穿越时空的过程中发生变化，并且虫洞头部可能以不同方式在同一顶点处散射多次。例如，在图5.5中，只有一个顶点被头部实际访问，并且它被访问了两次。第一次相遇时，头部被对角散射；而第二次，头部则被水平散射。分配给这个顶点的图元不是固定的，而是在每次头部访问该顶点时随机确定的。从这个意义上说，该算法并不完全等同于环回算法的顺序执行。虽然我们在前面的讨论中使用了“向顶点分配图元”这一表述，以强调与环回算法的联系，但“头部在顶点处的散射”是对该过程更恰当的描述。

经过 $N_{\mathrm{cycle}}$ 个虫洞更新周期后，我们移除这些顶点并重新开始，即放置一组新的顶点，对这组顶点进行 $N_{\mathrm{cycle}}$ 次虫洞更新，依此类推。该算法的伪代码在算法16中给出。

在原始的虫洞算法中，如 Prokov’ev 等人（1998）所提出的，扭结位置是在虫洞运动过程中即时插入的，使用了特殊的更新来插入或移除扭结。将顶点生成与虫洞运动分开与否，通常在效率上没有太大区别。然而，虫洞算法的并行化实现（Masaki 等人，2013）是基于预先确定的顶点位置的。另一方面，当系统非常稀疏时，大部分顶点从未被头部访问过，因为头部的运动总是伴随着粒子的产生/湮灭，这会严重降低效率。在这种情况下，即时生成顶点可能非常有利，因为无人区域的顶点甚至不会被创建。在第6.4.6节中，我们描述了一种本质上与原始虫洞算法相同，但更类似于图5.5所示算法的方法。

图5.5 自旋-¹ XY模型中使用虫洞更新的有向环回算法示意图。

最后需要指出，尽管我们以相同的自旋-¹ XY模型为例来使环回算法和虫洞算法之间的对应关系更清晰，但对于该模型，引入虫洞实际上并非必要，因为环回算法在此情形下工作良好。环回更新克服了（计算上的）临界慢化问题，而离散时间世界线算法中的局部更新则受此困扰。此外，通过简单地将整个时空域划分为小块，将环回更新的计算任务分布到多个处理器上也相对容易。⁶ 然而，对于其他模型，虫洞算法可能高效得多。当哈密顿量包含相互竞争的项时，环回更新的性能通常非常差。一个典型的例子是存在均匀外磁场的反铁磁海森堡模型，其中交换耦合倾向于相邻自旋反平行排列，而外磁场则倾向相反排列。如果温度较低，环回算法会遇到一个瓶颈：反铁磁排列与铁磁排列之间的路径被一个高自由能势垒阻断，即使两者对热力学分布都有不可忽略的贡献。因此，尽管此例中不存在臭名昭著的负符号问题，但环回更新在低温下存在平衡化问题。可以通过利用上述“微妙之处”来解决，即虫洞更新不像环回算法那样涉及“固定”的图元。也就是说，在虫洞算法中，所有相互竞争的矩阵元都可以在虫洞的散射概率中被考虑（从而消除了阻挫）。在将相互作用分解为图元方面，虫洞算法比环回算法更灵活。

算法16 自旋- $- \frac 12$ XY模型中使用虫洞更新的有向环回算法。输入：世界线构型。在所有结点处放置顶点；对于所有相邻位点对：生成新的顶点位置 $\tau \in [ 0 , \beta )$ ，密度为 ${\scriptstyle{\frac{1} {4}}} J^{x y} ,$ 并在这些位置放置顶点；结束循环；对于 $k = 1$ 到 $N_{\mathrm{cycle}}$：在时空域中均匀随机选取一个点；在该点放置头部和尾部；选择头部初始方向，概率为 $\frac{1} {2}$；循环：将头部移动到下一个顶点；如果头部碰到尾部：让头部和尾部湮灭并退出循环；否则如果头部碰到结点上的顶点：以概率 $\frac{1} {2}$ 将头部对角或水平散射；否则：根据世界线构型散射头部；只有一种方式满足局部守恒律。结束判断；结束循环；结束循环；返回更新后的世界线构型。

5.2.7 遍历性与绕数#

严格来说，使用局部更新的世界线方法不是遍历的。它守恒绕数。如我们之前所讨论的，环回或虫洞算法用全局采样取代了局部面采样，这种变化使得绕数可以波动，并使不等时关联函数的测量以及巨正则系综中的模拟更易进行。虽然算法的效率可能不强烈依赖于系综之间的差异，但通过使用巨正则系综，我们可以直接访问一些有趣的物理量。

如果使用空间周期性边界条件（虚时方向总是周期性的），时空晶格位于一个 $( d + 1 )$ 维环面上。绕数是路径在回到其起点之前绕环面的圈数。空间方向上的绕数与玻色系统中的超流密度和自旋系统中的螺旋模量（helicity modulus）密切相关（Pollock 和 Ceperley, 1987）。具体来说，我们通过 $\mu \left( = x , y , \ldots \right)$ 方向的绕数 $W_{\mu}$ 得到超流密度 $\rho_{s}$ 为：

$$ \rho_{s} = \frac{L_{\mu}^{2}} {\Omega} \langle W_{\mu}^{2} \rangle ,\tag{5.26} $$

其中 $\Omega \equiv L_{\tau} L_{x} L_{y} \cdot \cdot \cdot$，在玻色子模型中跳跃常数 $t$ 时，$L_{\tau} \equiv 2 t \beta$。在 XY 模型中，$L_{\tau} = J^{x y} {\beta}$ 是时间方向上的无量纲长度。类似地，我们得到无量纲压缩率 $\tilde{\kappa} \equiv 2 t \partial \langle n \rangle / \partial \mu$ 为

$$ \tilde{\kappa} = \frac{L_{\tau}^{2}} {\Omega} \langle W_{\tau}^{2} \rangle ,\tag{5.27} $$

其中 $W_{\tau}$ 是时间绕数。这些量在正常态（$\rho_{s}$）和不可压缩态（$\tilde{\kappa}$）中呈指数级小。因此，对于大晶格，在这些态中非零绕数构型的概率极低，所以看似灾难性的局部更新非遍历性实际上并不构成严重问题。即使在超流态或可压缩态中，热力学量通常对正则系综和巨正则系综之间的差异也不敏感。然而，使用巨正则系综通常是有利的，因为 $\rho_{s}$ 和 $\kappa$ 通常是表征量子系统态的相关响应函数。

5.3 玻色子和费米子#

量子蒙特卡罗方法不局限于有限温度下量子自旋系统的模拟。玻色子和费米子系统，以及量子多体系统的基态性质，也具有相当大的研究意义。并非所有有限温度算法都能通过环回/团簇或虫算法方便地以连续时间形式表述。即使可以，这些算法的效率也无法保证。大多数量子蒙特卡洛算法，无论是针对量子自旋、玻色子还是费米子，都可能面临符号问题。

费米子系统的算法具有不同于量子自旋和玻色子算法的特征。然而，有些玻色子系统是用与量子自旋算法非常相似的算法来模拟的。具有不同量子统计的一维系统可以相互映射。在某些情况下，这种等价性可能允许直接采用现有算法，而非重新创建算法。在以下小节中，我们将简要讨论这些主题。

5.3.1 玻色子#

如今，环回和虫算法不仅是模拟一维，也是模拟二维和三维量子自旋模型的标准方法。它们与这些模型典型的短程交换相互作用很好地契合。环回和虫算法同样适用于硬核玻色子模型的模拟。⁷ 针对这种量子统计的 world-line 算法的推导过程与量子自旋情形类似，基态是玻色子数算符的本征态。

在一维情况下存在一条捷径。我们无需重新推导玻色子 world-line 算法的转移概率，而是使用松原-松田变换（Matsubara-Matsuda 变换，1956）将硬核玻色子模型映射到自旋- ¹模型，然后采用自旋算法。

具体来说，松原-松田变换

$$ S_{i}^{+} a_{i}^{\dagger} , S_{i}^{-} a_{i} , S_{i}^{z} + \textstyle{\frac{1} {2}} a_{i}^{\dagger} a_{i} = n_{i} $$

将硬核玻色子映射到自旋 $-\frac{1}{2}$，反之亦然。例如，如果我们将此变换应用于纵向场 XXZ 模型

$$ H = - \sum_{i = 1}^{N} \left[ J^{x y} \left( S_{i}^{x} S_{i + 1}^{x} + S_{i}^{y} S_{i + 1}^{y} \right) + J^{z} S_{i}^{z} S_{i + 1}^{z} + H^{z} S_{i}^{z} \right] $$

我们得到

$$ H = - t \sum_{i = 1}^{N} \left( a_{i}^{\dagger} a_{i + 1} + a_{i + 1}^{\dagger} a_{i} \right) - \mu \sum_{i = 1}^{N} n_{i} + V \sum_{i = 1}^{N} n_{i} n_{i + 1} + \mathrm{c o n s t .} , $$

其中 ${t} = {\textstyle{\frac{1} {2}}} J^{x y} , \ V = J^{z}$，且 $\mu = H^{z} - J^{z}$。由此得到的模型称为扩展硬核玻色-哈伯德模型（extended hard-core Bose-Hubbard model）。“扩展”一词源于模型中存在最近邻耦合，而“硬核”则表示在位排斥无穷大。在第6.6节中，我们将讨论标准的玻色-哈伯德模型，其中仅包含有限大小的在位排斥。

由此可见，这种映射非常直接。毫不意外，用于消除量子自旋算法中潜在符号问题的技术在此同样适用。对于这两类量子统计，哈密顿量矩阵元的符号通常由模型参数的符号决定。这使得自旋和玻色子算法的构造与费米子算法所面临的挑战有所不同。

5.3.2 费米子#

费米子世界线算法的推导过程与量子自旋类似。在对哈密顿量各部分进行指数化时，为了体现不同的量子统计，会出现细微差别。基态通常是费米子数算符的本征态。问题不在于无法构造算法，而在于很难构造出能高效工作的算法。

为了更好地理解困难的根源，我们考虑一个特例，其中约旦-维格纳变换（Jordan-Wigner transformation）（Jordan and Wigner, 1928）将一维费米子模型映射到量子自旋-$\frac{1}{2}$模型。我们应用该变换，定义如下：

$$ c_{i} = K_{i} S_{i}^{+} , \quad c_{i}^{\dagger} = S_{i}^{-} K_{i}^{\dagger} , \quad K_{i} \equiv e^{i \pi \sum_{j = 1}^{i - 1} S_{j}^{+} S_{j}^{-}} , $$

将其应用于一维无自旋无相互作用费米子模型

$$ H = - t \sum_{i = 1}^{N} \Big ( c_{i}^{\dagger} c_{i + 1} + c_{i + 1}^{\dagger} c_{i} \Big ) . $$

算符 $K_{i}$ 称为纽结算符（kink operator），它统计了位置 $i$ 左侧自旋从下到上翻转的次数。这自动计入了由两个粒子交换产生的费米子符号。纽结算符是酉且自共轭的，

$$ K_{i}^{\dagger} = K_{i} , \quad K_{i}^{2} = I , $$

并且与 $j \geq i$ 时的 $S_{j}^{\pm}$ 对易。利用

$$ c_{i}^{\dagger} c_{i + 1} = S_{i}^{+} S_{i + 1}^{-} , \quad c_{i + 1}^{\dagger} c_{i} = S_{i + 1}^{+} S_{i}^{-} , \quad c_{i}^{\dagger} c_{i} = S_{i}^{+} S_{i}^{-} ,\tag{5.28} $$

我们得到

$$ H = - 2 t \sum_{i = 1}^{N} \big ( S_{i}^{x} S_{i + 1}^{x} + S_{i}^{y} S_{i + 1}^{y} \big ) , $$

即 $J^{x y} = 2 t$ 的 XY 模型。因此，我们可以利用针对该自旋模型讨论的算法来模拟无自旋费米子。⁸ 当我们在问题中加入费米子自旋时，就必须为自旋向上和自旋向下的粒子分别引入世界线。然后，模拟还需要额外的蒙特卡罗步骤，根据自旋向上和自旋向下之间的相互作用来接受或拒绝世界线的移动（Kawashima et al., 1994）。

让我们强调一下自旋升降算符与产生湮灭算符之间映射中的纽结算符。在松原-松田变换中并没有类似的算符。回顾一下，我们的基态是直接乘积态，以特定顺序列出格点。只要格点按顺序编号，约旦-维格纳变换就适用于任意维度。之所以难以将上述一维解推广到更高维度，是因为不可能以这样一种方式对所有格点进行编号，使得所有相互作用的格点对在列表中彼此相邻。(如果可能，则意味着该晶格是一维的。) 更具体地说，(5.28) 的一般形式是

$$ c_{i}^{\dagger} c_{j} = e^{i \pi \sum_{k = i + 1}^{j - 1} S_{k}^{+} S_{k}^{-}} S_{i}^{+} S_{j}^{-} , $$

⁸ 但需注意，边界条件可能会改变：周期边界的费米子模型会映射到反铁磁的自旋模型，反之亦然。

带着一个丑陋的相位因子，当 $| j - i | > 1$ 时该因子不为零，并使问题在一般情况下不可解。这个符号问题实际上限制了世界线方法只能用于一维系统。

简言之，模拟费米子问题的困难源于费米子间排列产生的符号，该符号由相位因子 $e^{i \pi \sum_{k = i + 1}^{j - 1} n_{k}} = \pm 1$ 表示。如果我们暂时忽略这个相位因子，费米子哈密顿量与硬核玻色子或 $\begin{array} {r} {S = {\frac{1} {2}}} \end{array}$ 自旋模型完全相同。因此，我们可以将针对自旋或玻色子的环回/团簇（或任何其他）算法也应用于费米子。然而，要以一种不会导致不同蒙特卡洛构型之间严重抵消的方式来考虑相位因子

$$ S_{i}^{+} = c_{i} K_{i} , \quad S_{i}^{-} = c_{i}^{\dagger} K_{i} , $$

是非常困难的，这被称为"费米子符号问题"或"负符号问题"。迄今为止，这个问题没有通用的解决方案。

5.4 负符号问题#

负符号问题是量子模型数值模拟中最大的瓶颈。对于我们这里考虑的算法，它源于哈密顿量的负矩阵元。当局部哈密顿量 H 的某些非对角矩阵元为负时，对于某些构型 C，权重 $w ( C )$ 通常为负。在这种情况下，我们执行蒙特卡洛模拟，其目标分布是 $| w ( C ) |$，而不是 $w ( C )$。然后，我们利用恒等式计算可观测量 Q 的期望值

$$ \begin{array} {r l} & {\langle Q \rangle = \frac{\sum_{C} \left| w ( C ) \right| \mathrm{sign} ( C ) Q ( C )} {\sum_{C} \left| w ( C ) \right| \mathrm{sign} ( C )}} \\ & {\quad =frac{\sum_{C} \left| w ( C ) \right| \mathrm{sign} ( C ) Q ( C ) / \sum_{C} \left| w ( C ) \right|} {\sum_{C} \left| w ( C ) \right| \mathrm{sign} ( C ) / \sum_{C} \left| w ( C ) \right|}} \\ & {\quad =frac{\langle \mathrm{sign} ( C ) Q ( C ) \rangle_{\mathrm{MC}}} {\langle \mathrm{sign} ( C ) \rangle_{\mathrm{MC}}} ,} \end{array}\tag{5.29} $$

其中 sign $( C ) = \pm 1$ 表示 $w ( C )$ 的符号，$\langle \cdots \rangle_{\mathrm{MC}}$ 是以 $| w ( C ) |$ 为权重的蒙特卡洛平均。

配分函数的负贡献 $Z_{-}$ 抵消了正贡献 $Z_{+}$ 的一部分，而总和 ${Z = Z_{+} - Z_{-}}$ 必须始终为正。事实上，在许多情况下，负贡献几乎完全抵消了正贡献。我们可以通过考虑一个虚构的哈密顿量 $H^{\prime}$ 来理解这一点（Loh 等人，1990；Hatano 和 Suzuki，1992），该哈密顿量的矩阵元是原始哈密顿量 H 相应矩阵元的绝对值：$\langle C^{\prime} | H^{\prime} | C \rangle \equiv$ $| \langle C^{\prime} | H | C \rangle |$。当两个哈密顿量之间的差异是广延量时，每个格点的自由能差 $\Delta f \equiv ( F^{\prime} - F ) / N$ 为 $\mathcal{O} ( 1 )$。由于 $\beta{\cal F} = \ln Z$ 且 $\beta{\cal F}^{\prime} = \ln Z^{\prime}$，其中 $Z^{\prime}$ 是虚构系统的配分函数，因此（5.29）中的分母可以表示为

$$ \frac{Z_{+} - Z_{-}} {Z_{+} + Z_{-}} = \frac{Z} {Z^{\prime}} = e^{- \beta N \Delta f} . $$

从这个表达式可以清楚地看出，在低温和系统尺寸较大的情况下，$Z$ 的正部和负部几乎完全抵消，使得从（5.29）计算期望值 $\langle Q \rangle$ 实际上变得不可能，因为分子和分母中的统计误差都将远大于平均值。

5.5 动力学#

平衡量子蒙特卡洛结果拥有物理上有意义的动力学信息，而经典蒙特卡洛结果则没有。能够获取真实的动力学信息是对我们在模拟中必须处理额外虚时维度的回报。要获取这些信息，必须进行特定类型的测量。根据定义，量子力学算符 X 的期望值为

$$ \langle X \rangle = {\frac{\operatorname{Tr} e^{- \beta H} X} {\operatorname{Tr} e^{- \beta H}}} . $$

虚时 $\tau$ 的算符为 $X ( \tau ) = e^{\tau H} X e^{- \tau H}$，我们有

$$ \langle X ( \tau ) \rangle = \frac{\mathrm{Tr} e^{- ( \beta - \tau ) H} X e^{- \tau H}} {\mathrm{Tr} e^{- \beta H}} . $$

由于平衡态的时间平移不变性，有 $\begin{array} {r l} {\langle X \rangle} & {{} =} \end{array}$ $\begin{array} {r} {{\frac{1} {\beta}} \int_{0}^{\beta} d \tau \langle X ( \tau ) \rangle} \end{array}$ 。这类期望值称为等时期望值（equal-time expectation values）。

动力学探针涉及两个（或多个）不同时刻算符的乘积

$$ \left. A ( \tau ) B ( 0 ) \right. = \frac{\mathrm{Tr} e^{- ( \beta - \tau ) H} A e^{- \tau H} B} {\mathrm{Tr} e^{- \beta H}} . $$

这些期望值同样具有时间平移不变性：$\left. A ( \tau ) B ( 0 ) \right. =$ $\langle A ( \tau + \tau^{\prime} ) B ( \tau^{\prime} ) \rangle$

非等时测量的重要性在于它们与谱密度（spectral density）的关系，

$$ \langle A ( \tau ) B ( 0 ) \rangle = \int_{- \infty}^{\infty} d \omega{\frac{e^{- \tau \omega} D ( \omega )} {1 \pm e^{- \beta \omega}}} , $$

其中 $\pm$ 取决于 $A$ 和 $B$ 是反对易还是对易。谱密度 $D ( \omega )$ 将虚时动力学与真实动力学信息联系起来。更具体地说，谱密度包含了关于系统元激发（elementary excitations）的可测量信息。例如，测量 $\langle S_{i}^{z} S_{j}^{z} \rangle$ 的空间傅里叶变换的虚时依赖性，并对积分方程求逆，得到 $D ( k , \omega )$ （Kawashima et al., 1996），这通常通过非弹性中子散射实验测量。该函数中的峰是 $z$ 分量自旋激发的标志，峰的位置作为 $k$ 的函数给出激发的色散关系 $\omega ( k )$ 。即使没有蒙特卡洛估计伴随的统计误差，对积分方程求逆也是一项困难的任务。我们将在第12章讨论这种求逆技术。更易于处理的是静态磁化率（static susceptibilities），通过对虚时积分得到。例如，静态纵向和横向自旋磁化率为

$$ \begin{array} {l} {{\displaystyle \chi_{i j}^{\parallel} = \int_{0}^{\beta} d \tau \left. S_{i}^{z} ( \tau ) S_{j}^{z} ( 0 ) \right.} ,} \\ {{\displaystyle \chi_{i j}^{\perp} = \int_{0}^{\beta} d \tau \left. S_{i}^{x} ( \tau ) S_{j}^{y} ( 0 ) + S_{i}^{y} ( \tau ) S_{j}^{x} ( 0 ) \right.} .} \end{array} $$

对于对角关联函数的测量，我们使用恒等式

$$ \langle S_{i}^{z} ( \tau ) S_{j}^{z} ( 0 ) \rangle \equiv \langle s_{i} ( \tau ) s_{j} ( 0 ) \rangle_{\mathrm{MC}} , $$

其中 $\langle \cdots \rangle_{\mathrm{MC}}$ 表示对蒙特卡洛模拟生成的环状位形（world-line configurations）的简单平均。虽然非对角关联函数的测量不像这样直接，但当我们使用环/团簇更新（loop/cluster update）或蠕虫更新（worm update）时，公式几乎同样简单。例如，我们可以使用以下方法估计 $\operatorname{spin} {-} \frac{1} {2} X Y$ 模型中 $S^{x}$ 算符之间的关联函数

$$ \langle S_{i}^{x} ( \tau ) S_{j}^{x} ( 0 ) \rangle \equiv \frac{1} {4} \langle \theta ( ( i , \tau ) , ( j , 0 ) ) \rangle_{\mathrm{{M C}}} , $$

其中 $\theta ( X , Y )$ 在时空点 $X$ 和 $Y$ 位于同一环上时为 1，否则为 0。通过对 $X$ 和 $Y$ 积分，我们得到均匀、静态的横向磁化率，

$$ \begin{array} {l} {{\displaystyle \chi^{\perp} = \frac{1} {N \beta} \sum_{i , j} \int_{0}^{\beta} d \tau_{1} d \tau_{2} \langle T S_{i}^{x} ( \tau_{2} ) S_{i}^{x} ( \tau_{1} ) \rangle}} \\ {{\displaystyle \quad = \frac{1} {4} \frac{1} {N \beta} \sum_{i , j} \int_{0}^{\beta} d \tau_{1} d \tau_{2} \langle \theta ( ( i , \tau_{2} ) , ( j , \tau_{1} ) ) \rangle_{\mathrm{MC}}}} \\ {{\displaystyle \quad = \frac{1} {4} \frac{1} {N \beta} \left. \sum_{i} \Lambda_{l}^{2} \right._{\mathrm{MC}} ,}} \end{array} $$

其中 $\tau$ 表示编时乘积，$\Lambda_{l}$ 是由环索引 l 指定的环的长度。类似地，在蠕虫算法中，我们通过简单测量头部相对于尾部访问特定位置的频率来估计关联函数。对角关联函数的简单高效测量是环/团簇算法和蠕虫算法的重要优势。我们将在第6章更详细地讨论这些要点。

建议阅读#

M. Creutz and B. Freedman, “A statistical approach to quantum mechanics,” Ann. Phys. 132, 427 (1981).

R. T. Scalettar, “World-line quantum Monte Carlo,” in Quantum Monte Carlo Methods in Physics and Chemistry, series C, volume 525, NATO Science Series, ed. M. P. Nightingale and C. J. Umrigar (Dordrecht: Kluwer Academic, 1999).

H. G. Evertz, “The loop algorithm,” Adv. Phys. 52, 1 (2003).

N. Prokof’ev and B. Svistunov, “Worm algorithm for problems in quantum and classical systems,” in Understanding Phase Transitions, ed. L. D. Carr (Boca Raton, FL: Taylor and Francis, 2010).

A. W. Sandvik, “Computational studies of quantum spin systems,” in Lectures on the Physics of Strongly Correlated Systems XIV, ed. A. Avella and F. Mancini, AIP Conference Proceedings, 1297, 135 (2010).

习题#

5.1 守恒法则。证明：如果量子自旋哈密顿量在 U(1) 旋转下不变，

$$ S_{i}^{+} \to e^{i \phi} S_{i}^{+} , S_{i}^{-} \to e^{- i \phi} S_{i}^{-} ( \mathrm{f o r ~ a l l} i \in V ) , $$

则该哈密顿量在 $\textstyle \sum_{i \in V} S_{i}^{z}$ 的不同特征值特征态之间的矩阵元为零，其中 V 是哈密顿量定义的区域。

5.2 泊松过程最熟悉的例子是具有衰变率 λ 的放射性原子的衰变。对于这种衰变：

证明原子存活时间超过某个 τ 的概率为

$$ P_{\mathrm{survive}} ( \tau ) = ( 1 - \lambda \Delta \tau )^{\tau / \Delta \tau} e^{- \lambda \tau} \quad \text{当} \quad \Delta \tau \to 0 . $$

证明第一次衰变事件发生在短时间区间 $[ \tau , \tau + \Delta \tau )$ 内的概率为

$$ P_{\mathrm{survive}} ( \tau ) \lambda \Delta \tau = \lambda \Delta \tau e^{- \lambda \tau} . $$

通过使用两个泊松分布的乘积得到相同结果：第一个计算在区间 $[ 0 , \tau )$ 内无事件的概率，第二个计算在短时间区间 $[ \tau , \tau + \Delta \tau )$ 内恰好发生一个事件的概率。

5.3 假设有两个独立的泊松过程，时间区间为 τ，其特征分别为事件数 $(N_1, \lambda_1)$ 和事件密度 $(N_2, \lambda_2)$。令 $N = N_1 + N_2$ 且 $\lambda = \lambda_1 + \lambda_2$。

证明 $P(N = 0) = e^{-\lambda \tau}$，$P(N = 1) = \lambda \tau e^{-\lambda \tau}$。
将上述结果推广，证明两个独立的泊松过程具有密度为 λ 的联合分布。
证明在给定 N 的条件下，$N_1$ 的分布是指数为 N、参数为 $p = \lambda_1 / (\lambda_1 + \lambda_2)$ 的二项分布（公式 2.34）。

5.4 假设 x 服从参数为密度 λ 的泊松分布，且给定 x 时，y 服从指数为 k、参数为 p 的二项分布（习题 2.2）。证明 y 的无条件分布是密度为 $p \lambda = E[E[y | x]]$ 的泊松分布。

5.5 验证公式 (5.12) 和 (5.13)。

5.6 证明

$$ e^{A} C e^{- A} = C + [ A , C ] + \frac{1} {2 !} \left[ A , [ A , C ] \right] + \frac{1} {3 !} \left[ A \left[ A \left[ A , C \right] \right] \right] + \cdots \tag{5.30} $$

5.7 将算法15中的步骤与图5.4中的面板联系起来，注意有些图的分配是概率性的。在图5.4的中间面板中，识别出所有的环。共有四个。哪些环被翻转得到右图？

5.8 表示 $\exp(-\tau S^z)$ 和 $\exp(-\tau S^x)$ 的矩阵是什么？

5.9 为了满足正则性，松原-松田（Matsubara-Matsuda）变换和约旦-维格纳（Jordan-Wigner）变换必须保持算符的对易关系。验证这两种变换及其逆变换是否满足这一条件。利用这两种变换，建立直接在费米子和玻色子之间进行变换的关系。

第二部分有限温度