第5章 · 量子蒙特卡洛模拟与随机级数展开方法#

费曼的量子统计力学路径积分表述[175]在QMC方法的发展中发挥了重要作用。对于自旋系统及相关晶格模型，基于虚时路径积分的方法通常被称为世界线方法（world line methods）[176]。这些技术最初基于虚时的近似离散化——玻尔兹曼算符$\exp ( - \beta H )$的Suzuki–Trotter分解[177, 178]。后来，直接在虚时连续统中运行的精确算法也被开发出来[179, 32]。然而，第一个实际有用的QMC方法并未使用路径积分。早在20世纪60年代初，Handscomb基于$\exp ( - \beta H )$的幂级数展开以及置换算符（海森堡交换可以用其表示）乘积迹的可精确计算性，为海森堡铁磁体开发了一种无近似的方法[180, 181]。尽管这一方案也被推广到海森堡反铁磁体[182, 183, 184]和其他一些系统[185]，但世界线方法通常更高效，并在很长一段时间内主导了该领域。随着更普遍适用的精确随机级数展开（SSE）公式[186, 187]的引入（其中也对迹进行采样，从而绕过了先前对置换算符代数的依赖），QMC计算的幂级数方法得以复兴。Handscomb方法的这种推广除了提高了计算实用性和效率外，还清晰地展示了离散幂级数与连续虚时路径积分之间的密切关联。基于这两种公式的现代QMC算法最终也颇为相似[33, 31]。特别是，对于某些模型（如海森堡模型），自旋构型都可以使用高效的环簇更新（loop-cluster updates）[188, 189, 179, 190]进行采样，这些是经典Swendsen-Wang[118]簇算法对量子系统的推广[191]。环概念进一步推广到“虫”（worms）[32]和“有向环”（directed loops）[190, 33]，具有更广泛的适用性（例如，在存在外场的情况下），并已经实现了对广泛量子自旋和玻色子模型的大规模研究。其中一些模型现在可以达到接近经典系统最先进水平的细节程度进行研究。

对于阻挫自旋系统，QMC方法通常会受到“符号问题”（sign problems）的困扰，即存在非正定的路径积分或级数展开[192]。虽然原则上仍可进行计算，并且在高温下控制符号涨落方面取得了一些进展[193]，但在低温下，由混合符号引起的统计误差变得不可控制。除了讨论阻挫模型符号问题的起源，我们这里只考虑无阻挫（二分）系统。尽管这仅构成了我们原则上感兴趣系统的一个子集，但这类模型仍然展现出丰富的有趣物理，并继续在前沿研究中提供重要的见解。

在本笔记中，我们将仅讨论SSE方法的实现和应用，对于自旋系统，它通常比连续时间路径积分更高效且技术上更简单。然而，理解这两种方案之间的关系仍然很有用。因此，在5.1节中，我们首先回顾量子统计力学中的路径积分和世界线QMC方法，包括离散和连续虚时两种情况。在5.1.3节中，我们介绍量子统计力学的一般级数展开公式，并研究其与路径积分的数学关系。在5.2节中，我们描述了针对$S = 1 / 2$海森堡模型的SSE方法的高效实现。在5.3节中，我们讨论了几个一维和二维系统的说明性计算结果；包括5.3.1节中的单链、5.3.2节中的梯子系统、5.3.3节中的均匀二维晶格、5.3.4节中的二聚化系统（重点讨论奈尔态与非磁态之间的量子相变），以及5.3.5节中的J-Q模型（包含连续和一级Néel-VBS相变的例子）。

5.1 量子统计力学的路径积分与级数展开公式#

量子统计力学的主要技术问题是如何处理玻尔兹曼算符 $\exp ( - \beta H )$。正如我们在上一节中看到的精确对角化方法，对于超过几十个自旋的系统，直接构造对应的矩阵是不可行的。路径积分方法提供了一种方法，将此算符的迹（配分函数）转化为可以使用蒙特卡洛方法采样的形式。另一种方法是从指数函数的幂级数展开出发。这里我们介绍这两种方法，作为在第5.2节中实现海森堡模型SSE方法的基础。

5.1.1 虚时路径积分#

路径积分表述的出发点是，将逆温度 $\beta$ 下的指数算符写作指数中包含 $\Delta_{\tau} = \beta / L$ 的 L 个算符的乘积；

这个迹可以表示为任意基底下对角矩阵元的和。我们还可以在每个指数算符之间插入一组完备态。于是，配分函数取 L 维矩阵元乘积和的形式；

形式上，指数算符等价于虚（欧几里得）时间 $t = - i \Delta_{\tau}$ 下的薛定谔时间演化算符 $\exp ( - i H t )$（取 $\hbar = 1$）。因此，可以将 (223) 视为一个 $\bar{d}$ 维量子系统到 $d + 1$ 维等效系统的映射，其中新增的维度是虚时间。$^{6}$ 状态指标 $l$ 对应一组离散的虚时间点 $\tau_{l} = l \Delta_{\tau}$，满足 $0 \le \tau_{l} \le \beta$ 且状态具有周期性时间边界条件 $( \alpha_{L} = \alpha_{0} )$。离散时间点通常被称为"时间片"。

从技术角度看，将 $Z$ 写成 (223) 形式的目的在于，虽然 $\exp ( - \beta H )$ 的矩阵元难以计算，但当 $\Delta_{\tau}$ 很小（即时间片数量很大）时，我们可以使用某种近似来评估 $\exp ( - \Delta_{\tau} H )$ 的矩阵元。然后，我们可以计算系统在所选基下可能演化的不同时间周期"路径" $\alpha_{0} \to \alpha_{1} \dotsb \to \alpha_{L - 1} \to \alpha_{0}$ 的权重。在QMC计算中，根据 (223) 中的权重对这些路径进行重要性采样。

为了讨论可能的时空路径及其权重，我们先考虑一种看似极端的近似指数算符的方法，即仅将其泰勒展开至 $\Delta_{\tau} H$ 的线性阶：

其中 $\{\alpha \}$ 统指所有态 $| \alpha_{0} \rangle , . . . , | \alpha_{L - 1} \rangle$。由于每个因子现在具有 $\Delta_{\tau}^{2}$ 量级的误差，并且有 $L = \beta / \Delta_{\tau}$ 个因子，因此在固定 $\beta$ 下 $Z$ 的相对误差为 $\Delta_{\tau}$ 量级。每个展开因子中被忽略的首项也含有 $H^{2}$，因此我们可能推测离散化误差也应与 $N^{2}$ 成比例。然而，要准确找出误差如何随 $\beta$ 和 $N$ 标度，需要比单纯数被忽略项在 $\langle H^{2} \rangle \propto N^{2}$ 和 $L = \beta / \Delta_{\imath}$ 层面上的数目更为细致的分析，因为对误差有贡献的路径与对近似 $Z$（写成(224)式）有贡献的路径并不相同。基于后面讨论的论证，为了获得与 $\beta$ 和 $N$ 无关的误差，在线性近似(224)中我们预计至少需要使用 $N \beta$ 量级的时间切片数，即 $\Delta_{\tau} \propto 1 / N$。

更好的时间切片算符近似可用于减少切片数量（使其与 $N$ 无关且与 $\beta$ 成正比，误差为 $\Delta_{\tau}^{2}$ 或更小）。然而，为了说明路径积分及基于它的QMC方法的许多基本方面，首先使用线性近似是方便的。事实上，正如我们稍后将看到的，甚至可以在该近似中取极限 $\Delta_{\tau} \to 0$，并基于精确的连续时间路径积分 [179, 32, 33] 来制定QMC算法。因此，我们对线性版本进行更详细的讨论，尽管在实践中不推荐使用 $\Delta_{\tau} > 0$ 的方案来实现实际程序。在理解了线性近似中路径的性质之后，将很容易理解如何取极限 $\Delta_{\tau} \to 0$ 或使用高阶离散近似，例如第5.1.2节讨论的Suzuki-Trotter分解。

玻色子路径积分和世界线。为了对路径及其如何采样有一个具体的例子，此时最好以特定的哈密顿量来继续讨论。首先考虑仅具有动能项的玻色子系统。我们在格点上工作$^{7}$，纯动能哈密顿量可以写为（将无关紧要的前因子设为1）；

其中 $a_{i}^{+}$ 和 $a_{i}$ 分别是格点 $i = 1 , \ldots , N$ 上的玻色子产生和湮灭算符（boson creation and destruction operators），而 $\langle i , j \rangle$ 是任意晶格（任意维度，但为便于可视化，我们将明确考虑一维链）上的一对最近邻格点。为简单起见，我们将考虑硬核玻色子（hard-core bosons），其格点占据数限制为 $n_{i} = 0 , 1$（可以认为这源于极强的格点内排斥）。硬核玻色子也与 $S = 1 / 2$ 量子自旋有简单的关系。占据和空置格点分别直接对应于自旋向上和向下，并且在无相互作用的情况下，玻色子哈密顿量等价于 XY 模型，其中 $H_{i j} = - ( S_{i}^{+} S_{j}^{-} + S_{i}^{-} S_{j}^{+} )$ 。(225) 中没有化学势——我们使用具有固定玻色子数 $N_{B}$ 的正则系综（密度 $\rho = N_{B} / N$，或在自旋语言中，磁化强度 $m = \rho - 1 / 2$）。稍后，我们将看到加入相互作用也很容易。

图50. 离散路径积分中各项的图形表示，针对一个具有14个格点的一维链上的三个玻色子系统。共有20个时间切片，标记为 $l = 0 , \ldots , 19$，对应于方程(224)中的态。玻色子用圆圈表示。将它们连接起来可以更清晰地显示时空晶格上世界线的形成。注意时间（垂直）方向上的周期性边界条件，即第一个和最后一个 $( l = 0 )$ 态是相同的。还需注意，使用线性时间切片算符 $\left( 1 - \Delta H \right)$ 时，在一个时间切片内不能有超过一个粒子发生跳跃（对角线线段）。在空间周期性边界条件下，全同粒子可以进行循环置换（缠绕），如右侧构型所示。

我们在格点占据数 $n_{i} , i = 1 , \ldots , N$ 的基矢下工作。考虑出现在方程(224)中的矩阵元 $\langle \alpha_{l + 1} | 1 + \Delta_{\tau} K_{i j} | \alpha_{l} \rangle$，此时我们也将其 H 的每次出现写成对所有单个跳跃项 $K_{i j}$ 的和（并注意负号相消）。那么，对于非零矩阵元，态 $| \alpha_{l} \rangle$ 和 $| \alpha_{l + 1} \rangle$ 之间有两种可能的关系：要么 $| \alpha_{l + 1} \rangle = | \alpha_{l + 1} \rangle$，要么 $| \alpha_{l + 1} \rangle = K_{i j} | \alpha_{l} \rangle$。在后一种情况下，最初必须在格点 i 处有一个粒子，该粒子被 $K_{i j}$ 移动到先前空置的格点 j，反之亦然。由于这些条件必须对所有连续的矩阵元 $l = 0 , . . . , L - 1 , 0$ 成立，每个粒子都必须遵循一条“世界线”（借用了相对论中类似时空路径的概念），在虚时传播的每一步，它要么停留在相同的空间位置上，要么移动（跳跃）一个晶格常数。图50展示了一维系统的两种这样的世界线构型。

再次注意虚时方向上的周期性边界条件，这源于(224)中对 $\alpha_{0}$ 的迹。这些边界条件适用于玻色子态（占据数），但并非严格适用于世界线。由于玻色子是全同粒子，它们在从 $\tau = 0$ 到 $\tau = \beta$ 的演化过程中可以被置换，并且仍然满足所需的时间周期性。在$D=2$的量子自旋系统中，唯一可能的置换是实现一个围绕系统周期性

最近邻跳跃的一维硬核系统是周期链上涉及所有粒子的循环置换，导致环上产生净粒子流。图50右侧展示了一个这种“缠绕”（winding）构型的例子。软核玻色子（即对占据数$n_{i}$没有限制）以及二维或三维中的硬核玻色子（甚至包括包含次近邻跳跃的一维情况）可以发生其他置换，但只有缠绕才会产生净电流。世界线缠绕系统的净次数（总电流除以系统长度）称为缠绕数（winding number）。它是一个拓扑量，表征构型的非局域性质。在更高维度中，每个空间方向都可能发生缠绕，对应有不同的缠绕数。玻色子系统缠绕的一个有趣方面是，它（经过适当归一化的缠绕数方差）对应着超流性[175, 196]。在自旋系统中，缠绕与自旋刚度（spin stiffness）或螺旋模量（helicity modulus）有关，我们将在后面进一步讨论。

现在我们将恰当世界线构型（即满足上述约束条件的构型）在配分函数中的权重记为$W ( \{\alpha \} )$；

在本文考虑的情况下，权重的形式非常简单：每个动力学操作$K_{i j}$都乘以(224)中的时间步长$\Delta_{\tau}$，并且在硬核情况下玻色子算符的矩阵元始终为1。因此，路径权重为

其中 $n_{K}$ 是动力学跳跃的总次数，对应于图50中的对角世界线段。该权重可以简单推广到软核玻色子情形，其中产生和湮灭算符与因子 $\sqrt{n_{i}}$ 和 $\sqrt{n_{i} - 1}$ 相关联。下面我们还将考虑玻色子之间的相互作用，这会导致动力学贡献(227)乘以更复杂的因子。

在QMC模拟中，人们对世界线进行修改，使得细致平衡成立，且构型根据 $W ( \{\alpha \} )$ 分布。世界线的局部更新如图51所示。注意，虽然这些局部移动在固定缠绕数的扇区内是遍历的，但它们无法改变拓扑缠绕数。在更高维度中，单根世界线的局部移动也无法导致任何其他需要包含的置换。涉及两根世界线的局部更新可用于对置换进行采样，但不足以改变与周期性边界相关的缠绕数。基于简单局部更新的算法在历史上曾成功研究了一些系统，但后来发展出了更复杂、更强大的环路更新 [31] 以及扩展构型空间中的相关更新方法 [32, 33]。这些方法在遍历地演化世界线方面效率更高（而不是“卡住”或在相空间中某些受限部分花费很长时间），并且也能改变缠绕数。

在本节中，我们主要关注路径积分表示，暂不讨论世界线采样的实际实现（我们稍后在本节中会对世界线采样进行一些更定量的讨论，但对于序列展开形式化中出现的相当类似的构型，我们只会详细实现其采样）。接下来，我们考虑我们可能想要计算的期望值的估计量。这也会让我们对世界线本身的性质有更深入的理解。

图51. 世界线的局部更新，它将动力学跳跃（对角世界线段）的数量改变 $\pm 2$。在(a)中，对角线段直接相连，而在(b)中，它们被一个时间片隔开。可以考虑跳跃间隔任意的移动。

期望值及其世界线估计量。考虑任意算符 $O$ 及其热力学期望值：

像我们在式(223)中对 $Z$ 所做的那样处理分子，期望值的精确时间切片形式可以写为

我们希望以适合于蒙特卡洛重要性采样的形式来表达这个期望值：

然而，这并不总是可行的，因为（229）中分子贡献的路径可能与 $Z$ 贡献的路径不同。如果对于贡献期望值的构型，$W(\{\alpha\})=\bar{0}$ 成立，那么就无法在（230）中定义估计量（estimator）$O(\{\alpha\})$。此时必须采用不同的方式继续计算。我们首先讨论那些适用于（230）的期望值类型。

对角算符（Diagonal operators）。 对于在占据数上对角的量，例如密度关联函数 $\langle n_i n_j \rangle$，形式（230）显然成立，因为在（229）中 $O|\alpha_0\rangle = O(\alpha_0)|\alpha_0\rangle$，其中 $O(\alpha_0)$ 表示 $O$ 的一个本征值。此时可以将估计量取为 $O(\{\alpha\})=O(\alpha_0)$。无论时间片演化算符 $\mathrm{e}^{-\Delta_\tau H}$ 如何近似，对角算符估计量的这种形式始终有效。由于迹的循环性质，算符 $o$ 可以被插入（229）中时间片算符乘积的任意位置，因此我们还可以对所有时间片进行平均，并使用

这通常会改善模拟的统计量。在实际应用中，为了节省时间而不显著损失统计量，可以在（231）中进行部分求和，例如每第N个时间片求和一次（通常 $L > N$），因为邻近状态的本征值 $O(\alpha_l)$ 和 $O(\alpha_{l+1})$ 差异很小（最多对应于一条世界线的跳跃）。而间隔 $\propto N$ 个时间片的状态通常有显著差异，至少在一定程度上贡献独立的统计量（具体程度可以通过虚时关联函数 $\langle {\cal O}(\tau) {\cal O}(0) \rangle$ 进行统计量化）。

动能（The kinetic energy）。 非对角算符的期望值通常更为复杂，但动能 $\langle K \rangle$ 是一个例外。它不仅易于计算，也是理解路径积分性质特别重要的量。在对 $Z$ 中的 $\mathrm{e}^{-\Delta \tau H}$ 进行线性近似的层面上，我们可以将（229）中 $O=K$ 的 $\mathrm{e}^{-\Delta_\tau H}$ 仅仅近似为 $K$，因为这仅在一个时间片上完成，并导致 $\Delta_\tau$ 阶的相对误差。该误差与（224）和（229）中所有其他时间片上的线性近似造成的总误差是同阶的。给定一个对 $Z$ 有贡献的构型，特定动能算符 $K_{ij}$ 的估计量为

这里，只有当第一个时间片上的格点 $i$ 和 $j$ 之间存在世界线跳跃时，分子才非零。在这种情况下，$K_{ij}(\{\alpha\})=1/\Delta_\tau = L/\beta$。在所有其他情况下，该估计量为零。我们还可以对（229）中算符 $K_{ij}$ 的所有时间片位置进行平均，得到

其中 $n_{ij}$ 表示在格点 $i$ 与 $j$ 之间的世界线构型中动能跳跃的次数。因此，总动能由总动能跳跃次数 $n_{K}$ 的平均值给出：$\langle K \rangle = - \langle n_{K} \rangle / \beta$。

我们本可以通过更简单的方式推导出动能的表达式，即使用内能（这里只有动能）的热力学公式 $\mathbf{\bar{\Delta}} E = \bar{\partial} \ln ( Z ) / \partial \beta$，其中 $Z$ 由式（226）和（227）给出。更复杂的推导过程明确展示了形式（230）如何涉及分子和分母构型空间的匹配。对于一般的非对角算符，这种匹配是不可能的；只有那些是哈密顿量一部分的算符才可能。在考虑其他情况之前，让我们讨论一下动能估计量的另一个重要方面。

表达式（233）的价值不仅在于它能让我们计算动能。它还传递了关于路径积分和世界线本质的基本信息。一个自然的问题是，随着时间片数量的增加，典型的世界线将如何演化？具体来说，在典型的世界线构型中，我们可以预期有多少个动能跳跃？

在低温下，动能应该几乎与温度无关（当 $\beta \to \infty$ 时趋近于常数）。它应与晶格尺寸 $N$ 成正比。利用式（233），我们可以因此推导出动能跳跃次数的预期大小和温度标度关系：

图52. 连续时间世界线构型。这里，动力学事件（跳跃）发生在任意虚时 $0 \leq \tau < \beta$ 处。构型 (a) 对配分函数有贡献，而构型 (b) 包含一对间隔超过一个格点间距的产生（实心圆）和湮灭（空心圆）算符，因此对 $Z$ 或对角期望值没有贡献。这种构型反而对非对角期望值 $\langle a_{i}^{\bar{+}} a_{j} \rangle$ 有贡献。

这告诉我们，世界线的粗糙度存在一个典型尺度；对于大的 $L$，跳跃之间的典型时间间隔为 $\sim 1 / N$。如本节开头所述，在对时间切片算符采用线性近似的离散路径积分中，我们应该使用约 $N \beta$ 量级的时间切片，以避免误差随 $N$ 和 $\beta$ 增长。我们现在对此有了解释，因为从式 (234) 可以清楚看出，更少的时间切片无法容纳正确描述动能（进而间接影响任何其他量）行为所必需的动能跳跃次数。

从另一个角度来看，$n_K$ 与 $\langle K \rangle$ 之间的关系也很重要；它表明，即使在连续极限下，如果我们让时间切片数量 $L \to \infty$，对有贡献构型中的动能跳跃（“事件”）数仍然有限，那么直接在此极限下构建模拟算法是可行的。连续时间世界线因此可以仅以一个状态 $\left| \alpha_{0} \right.$ 和 $n_{K}$ 个事件的形式存储：即事件发生的时间 $\tau_{i} (i=1,\dots,n_K)$ 以及每次跳跃的方向。由于 $n_K$ 的标度行为为 $N \beta$，因此应该有可能构建出时间和内存需求同样按系统大小和温度标度的蒙特卡洛采样算法。

图52展示了连续时间中的世界线构型示例。左侧的构型对配分函数有贡献，而右侧的构型则没有，因为它不满足所有时间周期性约束条件。该构型反而对算符 $a_{i}^{+} a_{j}$ 的期望值有贡献，其傅里叶变换是动量分布函数 $\rho ( k ) = \langle a_{k}^{+} a_{k} \rangle$。我们接下来简要讨论这类非对角的量。

非对角算符与断裂的世界线。如果我们考虑一个对 $Z$ 有贡献的路径 $\{\alpha \}$，并按照处理式 (232) 中动能项相同的方式来处理一般非对角算符 $a_{i}^{+} a_{j}$ 的期望值，那么很明显，除非 $i, j$ 是最近邻（此时该算符是动能的一部分），否则结果始终为零。另一方面，我们可以构建其他路径，这些路径对 $\langle a_{i}^{+} a_{j} \rangle$ 有贡献但对 $Z$ 无贡献，如图52(b)所示。如果格点 $i, j$ 间隔超过一个晶格间距，期望值 $\langle a_{i}^{+} a_{j} \rangle$ 就不能用标准的蒙特卡洛采样形式 (230) 来写，因为对于所有 $O ( \{\alpha \} ) \neq 0$ 的情况，都有 $W ( \{\alpha \} ) = 0$。然而，通过在对 $Z$ 有贡献的周期性世界线构型空间与那些包含与产生和湮灭算符相对应的两个缺陷的构型所组成的联合空间中进行操作，仍然可以计算 $\langle a_{k}^{+} a_{k} \rangle$ 以及相应的虚时依赖函数 $\langle a_{k}^{+} ( \tau ) a_{k} ( 0 ) \rangle$（单粒子格林函数）[32]。

包含相互作用的路径积分。将上述玻色子路径积分推广到除动能 $K$（式 (225)）之外还包含势能项的哈密顿量并不困难。用 $V$ 表示在占据数基下是对角的任何相互作用，例如形式为 $v_{ij} n_i n_j$ 的吸引 ($v_{ij} < 0$) 或排斥 ($v_{ij} > 0$) 项。令 $H = K + V$，我们可以将 $\mathrm{e}^{-\Delta_\tau H}$ 分解为：

该误差源于对易子 $[K, V] \neq 0$。由于 $V$ 是对角的，对于 (223) 中的每个时间片，我们可以写出

其中 $V_{l}$ 表示在第 $l$ 个时间片上计算的势能。因此，世界线构型的权重 (227) 在存在相互作用时被修改为

权重的相互作用部分可以很容易地在如图51所示的基本世界线移动中考虑，也可以在更复杂的连续虚时处理方法中考虑 [32]。

连续极限。我们已经讨论过，路径积分的连续极限可以直接用于 QMC 算法。然而，检查构型权重——对于纯动能哈密顿量是 (227)，存在相互作用时是 (237)——我们面临一个明显的问题：当 $\Delta_{\tau} \rightarrow 0$ 时，权重趋于零。这显然必须以某种方式通过构型数量的增加来补偿，使得配分函数保持有限。在蒙特卡洛计算中，我们不直接处理配分函数，只需要构型权重的比值来计算世界线更新的接受概率。考虑最简单类型更新的此类比值，即图51中所示的插入和移除相反动能跳跃的操作，它们在连续极限下也有问题：形式上，插入和移除的概率分别为零和无穷大。这不仅是连续极限下的问题，在 $\Delta_{\tau}$ 很小的情况下也同样存在，此时接受插入的概率虽然非零但非常小。

在基于离散路径积分的模拟首次引入 [177, 178] 二十年后，才在“虫”算法 [32] 和圈算法 [179] 的背景下引入了处理连续极限的正确方法。这里我们仅考虑一个玩具模型的局域世界线移动，以此说明实际上如何在连续极限下相当容易地构建蒙特卡洛抽样。

图51展示了引入两个事件、将世界线移动到相邻格点或从其移开的两种不同方式。考虑一段在时间 $\tau_{1}$ 和 $\tau_{2}$ 之间初始为直线的时间片，其长度为 $m$ 个时间片，即 $\tau_{2} = \tau_{1} + m \Delta_{\tau}$。假设在 $[ \tau_{1} , \tau_{2} ]$ 范围内相邻格点上没有世界线，那么引入这两个事件时的权重变化为 $\Delta_{\tau}^{2}$（对于纯动能哈密顿量），无论这些事件在此范围内的何处发生。插入事件的方式总共有 $m ( m - 1 ) / 2$ 种不同方式，这为此类修改构型子空间提供了总相对权重 $\Delta_{\tau}^{2} m ( m - 1 ) / 2$，而相对于少了两个事件的原始构型，其权重为 1。在连续极限下，所有此类可能更新构型的权重变为 $( \tau_{2} - \tau_{1} )^{2} / 2$。正如我们上面所讨论的，当 $\Delta_{\tau} \rightarrow 0$ 时，动能事件的平均总数保持有限。那么直线世界线段的标准长度 $\tau_{2} - \tau_{1}$ 也保持有限，因此我们应该能够构建一个具有有限插入概率的方案。关键在于，我们不仅考虑单个特定的更新构型，而是考虑一个连续的、可能的更新范围。

图53. 连续时间下的局域世界线移动。在 (a) 中，两个相反的动能跳跃被插入到现有跳跃 $\tau_{1}$ 和 $\tau_{2}$ 之间随机选定的时间 $\tau_{a}$ 和 $\tau_{b}$ 处。在 (b) 中，两个相反的跳跃要么被移除（），要么两个事件的时间被随机改变为 $\tau_{1}$ 和 $\tau_{2}$ 之间新的任意时间 $\tau_{a}^{\prime}$ 和 $\tau_{b}^{\prime}$（）。接受概率由式 (238) 给出。

我们还需要考虑移除两个事件。为简化讨论，这里仅考虑两个格点上的单个玻色子。这样就避免了还需要考虑相邻格点上世界线所施加约束的复杂性。我们将为其构造满足细致平衡概率的世界线移动如图53所示。前面讨论的事件插入如图(a)所示，需要注意的是，时间 $\tau_{1}$ 和 $\tau_{2}$ 对应两个连续的现有时事件（如果没有事件，我们取 $\tau_{1} = 0$ 和 $\tau_{2} = \beta$）。我们将两个新事件的时间称为 $\tau_{a}$ 和 $\tau_{b}$，它们在 $\tau_{1}$ 和 $\tau_{2}$ 之间的任意位置随机选取。在移除两个事件的相反更新中，我们考虑两个现有时事件 $\tau_{a}$ 和 $\tau_{b}$ 在 $\tau_{1}$ 和 $\tau_{2}$ 处的两个周围事件之间的总权重。然后我们可以再次计算现有时事件的总权重，不仅是在固定的当前时间 $\tau_{a}$ 和 $\tau_{b}$，而是在 $\bar{[} \tau_{1} , \tau_{2} ]$ 内的任意时间 $\tau_{a}^{\prime}$ 和 $\tau_{b}^{\prime}$。这个相对权重仍然是 $( \tau_{2} - \tau_{1} )^{2} / 2$，相比于缺失这两个事件的单一构型的权重为1。

因此，一旦通过检查当前构型（可以以事件列表的形式存储）识别出 $\tau_{1}$ 和 $\tau_{2}$，双事件插入和移除的接受概率为：

对于接受的插入，我们随机生成 $\tau_{a}$ 和 $\tau_{b}$，而在被拒绝的移除中，我们也为现有时事件生成新的时间 $\tau_{a}^{\prime}$ 和 $\tau_{b}^{\prime}$。

在更大的系统中（具有多个玻色子），上述方案必须加以修改，以考虑在移动选定世界线时其他世界线的约束。时间 $\tau_{1}$ 和 $\tau_{2}$ 应反映这些约束。存在相互作用的情况下，接受概率需要修改以考虑式(237)中的势能因子。然而，在实践中不会使用这些局部更新，因为圈算法和虫算法[32, 179, 33]的实现并不复杂多少，但效率却高得多。此处的讨论旨在说明，原则上可以如何克服连续极限的明显问题。

5.1.2 铃木-特洛特分解#

虽然现代晶格模型的量子蒙特卡洛世界线算法通常基于连续极限下的路径积分，但从历史上看，基于铃木-特洛特近似的时间切片算子 $\mathrm{e} ( - \Delta_{\tau} H )$ 的时间离散化变体是最早发展起来的。离散方法有时仍在使用，并且在某些情况下是最实用的选择[197]。这里我们讨论基于铃木-特洛特方法的主要特征，以及它们如何应用于玻色子和自旋系统。

铃木-特洛特（或分裂算子）近似[177, 198]的一个例子已在式(235)中给出。这只是两个（或多个）非对易算子指数分解为指数乘积的众多可能分解之一[199]。最常用的近似公式是，对于任意一对算子A和B以及一个小因子 $\Delta$：

这里的误差同样正比于对易子 $[ A , B ]$ 。通过使用更大数量、巧妙选取的关于 A 和 B 函数的指数，剩余误差可以进一步减小，原则上可以达到 ∆ 的任意高阶幂次 [200]。然而，高阶近似式在实践中往往过于复杂而难以操作。

如果我们像在路径积分中那样，关心的是指数乘积的迹，那么不难看出两种低阶近似式 (239) 实际上是等价的，因为

这源于迹的循环性质。因此，虽然世界线方法通常基于式 (239) 的第一行进行讨论，其中每个因子的误差 $\propto \Delta_{\tau}^{2}$，但实际误差更小，为 $\propto \Delta_{\tau}^{3}$。然而，由于路径积分中指数因子的个数为 $L = \beta / \Delta_{\tau}$，对于固定的 N 和 $\beta$，形如 (240) 的乘积中的总误差 $\propto \Delta_{\tau}^{2}$。此外，可以证明，对于足够大的 $L$，在固定 $\Delta_{\tau}$ 的情况下，误差也不会随 N 增长。因此，在大多数情况下，我们可以使 $\Delta_{\tau}$ 与 N 和 $\beta$ 无关（尽管对于某些类型的期望值，必须小心以避免当 $\beta \to \infty$ 时误差发散 [201, 202]）。

在路径积分的语境下，铃木-特洛特近似的实用性在于：尽管配分函数 (223) 中的 $\exp ( - \Delta_{\tau} H )$ 不容易计算，但如果我们能找到 H 的一种分解，将其拆为两项（或少量几项），即 $H = H_{A} + H_{B}$，那么就能精确计算 $\exp ( - \Delta_{\tau} H_{A} )$ 和 $\exp ( - \Delta_{\tau} H_{B} )$。典型的例子是一维系统，其只包含最近邻的动能和相互作用，对于一对格点 $i, i+1$，我们将它们统称为 $H_{i , i + 1}$。然后，我们按照下式将 $H$ 分解为 $H = H_{A} + H_{B}$：

通过这样的安排，每个求和 $H_{A}$ 和 $H_{B}$ 中的所有项都是相互对易的；对于任意 i 和 m，有 $[ H_{i , i + 1} , H_{i + 2 m , i + 2 m + 1} ] = 0$。因此，我们可以将 $\exp ( - \Delta_{\tau} H_{A} )$ 和 $\exp ( - \Delta_{\tau} H_{B} )$ 拆分成两体算符而不引入任何额外误差，例如：

在时间切片配分函数 (223) 中使用这样的乘积，我们可以在所有指数之间插入完备集。由于每个算符现在只涉及两个格点，矩阵元可以简化为以下形式：

对于硬核玻色子，这些元素构成一个 $4 \times 4$ 矩阵。由于键算符 $H_{i,j}$ 守恒所涉及两个格点上的玻色子数，该矩阵是块对角化的，包含两个单元素块（对应 $n_{\alpha i} = n_{\alpha j} = n_{\beta i} = n_{\beta j}$）和一个 $2 \times 2$ 矩阵，其对角和非对角元素均连接键上具有单个玻色子的态。因此，即使包含相互作用（只要不超过最近邻；否则需要采用更复杂的 $H$ 分解），所有矩阵元也很容易计算。

从图形上看，每个矩阵元 (243) 对应一个四格点面块（plaquette），根据占据数，该面块中有零条、一条或两条世界线片段穿过。考虑所有时间切片以及 $\exp( -\Delta_{\tau} H_{A} )$ 和 $\exp( -\Delta_{\tau} H_{B} )$ 中的所有格点对 $i, i+1$ ，这些面块形成一种棋盘图案，其中相邻两行对应一个时间切片，如图 54 所示。显然，哈密顿量分解 (241) 也常被称为棋盘分解。

对于守恒总粒子数的哈密顿量，世界线构型的约束条件与时间切片算符线性近似中的约束类似，但有一些重要修改。动力学跳跃只允许出现在图 54 的阴影面块上，但与图 50 所示的线性时间切片近似不同（每个时间切片中最多只能有一个跳跃），这里对单个时间切片内的跳跃次数没有进一步限制。世界线构型看起来仍然非常相似。当 $\Delta_{\tau} \to 0$ 时，它们在物理相关的演化意义上变得等价，即从某个时刻 $\tau$ 演化到更晚的 $\tau'$ ，且 $\tau' - \tau \gg \Delta_{\tau}$。

在蒙特卡洛模拟中，可以使用图 51 所示类型的简单更新来更新世界线，但需满足一个约束：对角环段（动力学跳跃）仅允许出现在图 54 的阴影面块上。更复杂的“环”和“有向环”更新（可在一次操作中同时移动多条世界线的大段片段）用于现代算法中 [31, 191, 33]（我们将在后面讨论随机级数展开方法时详细论述）。

图 54. 基于棋盘分解和Suzuki-Trotter近似的一维世界线构型。玻色子的动力学跳跃（或一对自旋的翻转）仅允许穿过阴影方块（面块）。宽度为 $\Delta_{\tau}$ 的一个时间切片由连续两行面块组成。右侧显示的六个孤立面块对应非零矩阵元，对于具有海森堡相互作用的自旋模型（世界线和空位分别对应 $\uparrow$ 和 $\downarrow$ 自旋），这些矩阵元由式 (244) 给出。

应用于海森堡模型。考察Suzuki-Trotter方法中路径权重的具体例子是有用的。让我们计算反铁磁海森堡相互作用的面积块矩阵元： $H_{i,i+1} = \mathbf{S}_i \cdot \mathbf{S}_{i+1}$ 。在这种情况下，式 (243) 中的玻色子占据数被替换为自旋态 $\uparrow$ 和 $\downarrow$ 。我们可以考虑在自旋之间形成的世界线（并且注意到，在图 54 这样的图片中，我们也可以为 $\downarrow$ 自旋绘制世界线；它们占据世界线未覆盖的所有格点，并在每个对角段处与这些世界线相交）。计算仅涉及直接的代数运算，我们仅列出六个允许（非零）的矩阵元结果：

世界线构型的权重是这些矩阵元的乘积，所有这些矩阵元都在图 54 的右侧用图形表示。注意非对角矩阵元前面的负号。对于允许的世界线构型，由于世界线的周期性约束，所有符号都会相互抵消。这也适用于应用于更高维二分晶格的世界线方法，但对于阻挫系统，由于同时存在负权重和正权重，会出现"符号问题"（我们将在第 5.1.3 节进一步讨论）。因此，在实践中，世界线方法和类似的 QMC 方法主要用于研究二分自旋系统和玻色子模型。对于费米子系统，世界线的置换也会导致符号问题，但一维情况除外，在一维中只存在全局循环置换（缠绕），并且可以通过选择周期或反周期边界条件来避免相关符号 [176, 187]。

在更高维中，可以使用 Suzuki-Trotter 方法，并进行类似于上述讨论的各种分解。我们为一维系统讨论的局域更新也可以适用于更高维。值得阅读文献 [79, 203]，这是两项开创性的世界线研究，它们牢固地确立了二维 $S = \bar{1 / 2}$ 海森堡反铁磁体的基本性质。

虽然在 Suzuki-Trotter 方法中使用的时间片数量不必依赖于系统尺寸（因为在一个时间片内可能发生 ∝ $N$ 个事件），但对世界线构型进行蒙特卡洛采样的计算复杂度与连续时间公式中的标度相同，即 $N \beta$ 。在精确的连续统公式中，通过仅存储和操作事件的发生时间来实现这种标度，而在 Suzuki-Trotter 方法中，通常使用完整的时空晶格（尽管原则上也可以只存储其中的事件）。因此，连续时间方法实际上通常在计算机上运行得更快，此外还不受离散化误差的影响。使用离散路径积分的唯一理由是在那些实际上难以取连续极限的问题中，例如包含耗散的有效模型，这会导致世界线之间的含时相互作用 [197]。

5.1.3 级数展开表示#

作为离散时间切片方法或路径积分连续极限的替代方案，还可以通过使用玻尔兹曼算符的泰勒展开来构建一个适合于蒙特卡洛采样的构型空间；

该方法由Handscomb首创[180, 181]，他开发了一种研究$S = 1 / 2$铁磁体的方法。后来，泰勒展开被更广泛地视为众多模型精确QMC算法的起点[182, 185, 186, 187]。对有限晶格而言，在有限$\beta$下，指数算符的幂级数展开是收敛的（这从技术上源于哈密顿量H的谱有界）。实际上，如下文所见，级数展开使得虚时连续统能够以离散形式表示，从而避免了近似，同时完全保留了有限可枚举构型空间的优势。

选取一个基组后，配分函数可写为

其中$\{\alpha \}_{n}$的下标表示需要对n个态求和。此表达式可与线性阶离散路径积分(224)进行对比。取H为(225)中的玻色子动能，我们显然可以绘制出与图50非常相似的世界线图来表示贡献项。主要区别在于“切片”数量，即展开幂次$n$是可变的，并且对于给定的n，每一“切片”上都会发生粒子跳跃。在这种表示中引入的新“传播”维度不同于虚时，但显然，这些态的标签$p = 0 , \ldots , n - 1$与路径积分中的虚时密切相关。我们将在下文进一步讨论这种确切的关系。

在只有动能的情况下，对于允许构型，(246)中的所有矩阵元都等于1，负号相互抵消，权重为

在更一般的情况下，该权重将根据(246)中矩阵元的乘积进行修正。与路径积分一样，对于非对角项中无阻挫的玻色子系统和量子自旋系统，这些权重是正定的（或者可以通过一些简单技巧使其变为正定）。我们在此假设权重正定，并将在下文进一步讨论其确切条件。

推导总内能$U = \langle H \rangle$的一般表达式是很有用的。对于任意哈密顿量，可将其写为

注意，每一项中都包含一个额外的H矩阵元；对于给定的n，求和涉及$n + 1$个态，而配分函数(246)中只有n个态。在(248)中有贡献的所有构型也对配分函数有贡献，但权重不同。我们可以通过以下书写形式明确匹配这些构型：

我们可以将关于n的求和扩展到包含$n = 0$项，因为该项为零。这样，(246)和(249)中的项就完全匹配，并且可以将$n / \beta$视为能量估计量。因此，根据配分函数中各构型的权重进行采样后，能量简单地由下式给出：

将H写为求和形式$\begin{array} {r} {H = - \sum_{i} H_{i}} \end{array}$ ，并将其应用于(246)和(248)中所有H的出现，我们可以推导出单个项$H_{i}$的期望值的类似表达式。该估计量的结果是该算符在配分函数展开的算符串中出现的平均次数；

对于动能项，这与路径积分表达式(233)相同。这表明，如果我们所谓路径只是指动能事件发生的顺序，而不涉及离散或连续路径积分中事件的时间点$\tau_{l}$，那么两种表述中的路径存在一一对应关系。这种对应关系非常普遍。区别在于，在级数展开中，势能项的处理方式与动能算符相同（即，不像路径积分中那样被重新指数化）。它们的存在对应于"无事件"，不影响世界线，但会关联到修正路径权重的对角矩阵元。而在路径积分中，则存在方程(237)中的势能因子，它依赖于动能事件的实际时间点。

方程(250)表明，尽管(246)中对展开阶数n的求和延伸到$n = \infty$，实际上可以截断到某个阶数$n_{\mathrm{max}} \propto N \beta$。为了进一步说明为何这只会导致指数级小且完全可忽略的误差，可以利用比热的估计量，该估计量通过取(250)对温度的导数得到；

当温度$T \to 0$时，$C$应趋于零，而(252)表明n分布的方差等于$\langle n \rangle$。因此显然，在某个幂次$n \propto N \beta$之外，分布呈指数级衰减。实际上，基于这种表述的量子蒙特卡洛（随机级数展开，SSE）算法会根据不同n区间的相对权重自动对n分布进行采样。

通过明确截断泰勒展开，可以使级数表示与路径积分表示之间的关系更加清晰。在$n = L$处截断，我们可以通过用$L - n$个单位算符I扩充所有幂次$H^{n}$（其中$n < L$）来形式化地构建一个固定大小的采样空间。考虑单位算符在$L$个算符乘积中的所有可能位置，我们可以定义一个修正算符串${\cal S} = {\cal S}_{1} , {\cal S}_{2} , \ldots , {\cal S}_{L}$，其中$S_{i} \in \{H , I \}$，并对所有这些序列进行求和。然后，我们需要通过等效项的数量$\binom{L} {n}$来补偿权重，得到

其中 $n$ 现在指的是算符串中 $S_{i} = H$ 的元素个数（我们不再需要对 $n$ 进行显式求和）。如果现在我们考虑玻色子动能式(225)并取极限 $L \to \infty$，权重简化为 $( \beta / L )^{n}$，这与 $L$ 个时间切片的路径积分权重式(227)相同。在此极限下，显然指标 $p = 0 , \ldots , L - 1$ 与虚时 $\tau = p \beta / L$ 是相关的。区别在于，采用完整权重 $\beta^{n} ( L - n ) ! / L !$ 的级数展开在实际中对于 $L \propto \beta N$ 是精确的，而在路径积分中我们必须取连续极限以避免离散化误差。对于有限的 $L$，级数指标 $p$ 并不精确对应虚时，而是代表一个虚时的分布[187]。由于传播指标 $p$ 与虚时之间的紧密关系，在级数展开中也将 $p$ 空间称为"时间"维度是恰当的。

级数展开的实际优势在于它提供了虚时连续统的一种精确但离散的表示。基于它的算法通常比连续时间世界线更容易实现，并且在计算上可能更高效。另一方面，对于势能项较大的哈密顿量，展开中算符的数量 $\langle n \rangle$ 可能远大于路径积分中动能事件的数量 $\left. n_{K} \right.$。后一种方法应该更高效，尽管精确处理需要操作浮点时间变量。对于量子自旋系统，动能（非对角）和势能（对角）能量之间通常存在有利的平衡，因此级数展开通常是更优的选择。接下来，我们将基于它开发一种QMC算法。

5.2 $S = 1 / 2$ 海森堡模型的SSE方法#

现在详细讨论如何针对 $S = 1 / 2$ 海森堡模型实现SSE方法。最初我们不需要指定晶格，而是考虑将哈密顿量写为键算符的和，

如式(204)所示，晶格编码为由键连接的位点列表 $[ i ( b ) , j ( b ) ]$，$b = 1 , \ldots , N_{b}$。图46给出了二维中的一个标记方案示例。

对于任何二分晶格，即所有位点 $i ( b )$ 和 $j ( b )$ 分别属于子晶格 A 和 B（在任何维度下），都可以构造正定的SSE。这构成了一大类有趣且重要的系统。最初我们只考虑单一的反铁磁耦合常数 $J_{b} = J > 0$，但对非均匀系统所需的修改非常简单。

5.2.1 构型空间#

将海森堡相互作用式(254)分解为它在对角 $z$ 自旋分量标准基中的对角部分和非对角部分是很有用的。然后我们定义具有两个指标的算符 $H_{a , b}$，其中 $a = 1 , 2$ 分别指对角和非对角，而 $b = 1 , \ldots , N_{b}$ 与之前一样是键指标；

这里我们在对角算符中引入了一个负号和一个常数，使得完整的哈密顿量可以写为

在 $H_{1 , b}$ 中包含这一常数的原因是为了使级数展开正定，我们很快就会看到。常数 $J N_{b} / 4$ 在算法中无关紧要，但我们在计算能量时会将其包含在内。

配分函数的级数展开。 SSE 方法中配分函数的一般形式可以写为 (246)，但这里我们最初不写出插入算符之间的所有完备态集。我们转而关注算符，将 (257) 中 $H$ 的所有实例展开为对所有键算符的求和。因此，海森堡模型 SSE 算法的起点是

这里 $J$ 已被吸收进 $\beta = J / T$，$S_{n}$ 指来自 $H^{n}$ 的海森堡键算符 (255) 和 (256) 的乘积（字符串）

而 $n_{2}$ 是非对角算符的总数，即字符串中 $b ( i ) = 2$ 的元素个数。注意，此处用于指代字符串中算符位置的标签 $p$ 取值 $0 , \ldots , n - 1$；在程序实现中，这种标记比 $p = 1 , \ldots , n$ 更方便。

当算符字符串作用于态 $| \alpha \rangle = | S_{i}^{z} , \ldots , S_{N}^{z} \rangle$ 时，我们会得到一系列其他基态，且不会分支为多个态的叠加。我们将这些称为传播态 $| \alpha ( p ) \rangle$

我们没有明确写出这些态的归一化，但将所有 $| \alpha ( p ) \rangle$ 视为已正确归一化。这些态当然对应于 (246) 中对完备集求和中的那些态，但在实践中，SSE 方法是围绕 (258) 中对算符字符串的求和来构建的，我们并不总是需要显式写出这些态。

由于 (255) 中的常数 $1 / 4$，所有对平行自旋的操作都会消灭该态；

因此，对 $Z$ 有贡献的算符-态构型 $( \alpha , S_{n} )$ 只能涉及对反平行自旋的操作，并且传播态 (260) 是在此假设下定义的。相应的矩阵元为

所有这些矩阵元都相等这一事实在取样算法中非常有用，这也是在对角算符中包含常数 $1 / 4$ 的原因。

除了仅作用于反平行自旋的局部约束外，对于 (258) 中完整算符乘积的矩阵元非零，传播还必须满足周期性条件 $\left| \alpha ( n ) \right. = \left| \alpha ( 0 ) \right.$ ，其中 $| \alpha ( 0 ) \rangle = | \alpha \rangle$。

虽然可以在长度波动的算符字符串空间 [180, 187, 204] 中制定蒙特卡罗取样程序，但在实践中通常更容易在上一节讨论的固定大小截断空间中进行操作（从数学上讲，在这样的空间中也更容易证明细致平衡条件）。$^{8}$ 在实际模拟中，截止值 $L$ 将由程序自动确定，如下文所述，使得 $L$ 安全地超过曾经取样的最大 n 值。再次强调，这种截断不会引起任何可检测的误差，不应被视为近似。将用于扩充长度为 $n < L$ 的字符串的幺正算符定义为 $H_{0 , 0} = I$，并将 $[ a ( p ) , b ( p ) ] = [ 0 , 0 ]$ 作为索引列表 (259) 中允许的元素，则配分函数可以写为：

其中 n 指的是固定长度算符字符串 $S_{L}$ 中非 $[0,0]$ 元素的数量（并且 n 通过对 $S_{L}$ 的求和隐式求和）。由于所有矩阵元都等于 $1 / 2$，因此允许的构型的权重由下式给出：

这里我们省略了符号因子 $(-1)^{n_{2}}$，如下文将进一步讨论的，对于二分（无阻挫）系统，该因子总是为正。还要注意 $L!$ 是一个不重要的归一化因子，也可以省略。

自旋 $\uparrow$ 和自旋 $\downarrow$ 分别由硬核玻色子的存在和缺失表示，我们可以绘制类似图50的图形来展示贡献SSE构型的规则。然而，现在我们还有不改变世界线但修改路径权重的对角算符。在路径积分方法中，对角算符全部被归并到一个公共的指数因子中，即式（237），但在级数展开中没有简单的方法做到这一点。如上文5.1节所述，当对角能量占主导时不应使用SSE方法，但在我们考虑的情况下，有 $\langle H_{1,b} \rangle \approx \langle H_{2,b} \rangle$，对角算符并不占主导。如下文所见，数目可观的对角算符的存在实际上在SSE算法中得到了利用。在讨论SSE方法时，绘制明确显示对角算符和非对角算符的图形非常有用。图55给出了一个示例，其中还包含了对以后在计算机实现该方法时会用到的数据结构的图示。状态将由整数 $\sigma(i) = \pm 1$ 表示，对应于 $S_i^z = \pm 1/2$。为了紧凑地存储算符字符串，将索引 $[a(p), b(p)]$ 组合成一个整数列表 $s(p) = 2b(p) + a(p) - 1$ 也是有用的。所有传播后的状态无需同时存储。它们可以根据需要从单个存储的状态 $|\alpha(0)\rangle$ 和算符字符串生成。稍后，我们还会引入一种不同的紧凑存储方式，其中也包含传播状态的一些自旋。

图55. 一个8自旋链的SSE构型，展示了所有传播后的状态。空心和实心条分别表示对角算符 $H_{1,i}$ 和非对角算符 $H_{2,i}$，而状态之间没有条则表示"...填充..."单位算符 $H_{0,0}$。状态 $|\alpha\rangle$ 的自旋 $\uparrow$ 和 $\downarrow$ 存储为 $\sigma(i) = \pm 1$，算符字符串 $S_L$ 使用偶数和奇数整数分别编码对角和非对角算符，根据 $s(p) = 2b(p) + a(p) - 1$。

图56. 一个由三个非对角操作（以条表示）将三角形上所有自旋带回其原始状态的示例。每个自旋翻转关联一个负号，导致路径权重为负，并在此阻挫系统及其他阻挫系统的QMC模拟中产生"符号问题"（由于不同符号的构型相互抵消）。

阻挫相互作用与"符号问题"。乍一看，由于式(258)中的因子 $(-1)^{n_2}$，我们似乎面临一个符号问题——非正定展开。实际上，对于二分晶格，所有项都是正的。这是因为在每一个允许的构型中，为了满足"时间"周期性 $|\alpha(L)\rangle = |\alpha(0)\rangle$，需要偶数个 $n_2$ 非对角算符。我们已经在世界线方法的上下文中讨论过这一点，其中式(244)中的非对角矩阵元为负，但世界线的连续性要求这些矩阵元有偶数个。这是两种方法之间密切关系的又一个例证。

对于阻挫系统，级数展开不是正定的（路径积分也完全出于同样的原因不是正定的）。这一点可以很容易地用三角形上的三个自旋系统来证明。如图56所示，在这种情况下，一个允许的构型可以包含三个非对角算符，导致整体负号。对于任何可以在反铁磁相互作用自旋之间形成具有奇数个位点环路的系统，这一点都成立——这可以作为阻挫的定义。

对于二分系统，正定性还可以通过另一种方式证明，即对其中一个子晶格（设为 $B$）上的自旋算符进行幺正变换，使得对于 $j \in B$，有 $S_{j}^{+} {-} \bar{S}_{j}^{+}$ 和 $S_{j}^{-} - S_{j}^{-}$（而对角算符 $S_{j}^{z}$ 不变）。这不会影响模型的能谱（因为所有自旋算符之间的对易关系保持不变），但哈密顿量（257）中对角项前的符号变为正。于是（258）式中的因子 ${\bar{( - 1 )}}^{n_{2}}$ 就不存在了。对于阻挫系统，不存在这样的变换可以移除所有负号。

注意，只有相互作用中的非对角部分会导致符号问题。我们可以研究对角部分存在阻挫的系统，但此时相互作用不再是自旋各向同性的。简单地忽略阻挫海森堡系统中的符号，对应于反铁磁的对角耦合和铁磁的非对角耦合。这里我们只考虑各向同性相互作用，并将讨论限制在无阻挫系统。

对于具有周期边界的二分系统，只有当所有方向 $\gamma = 1 , \ldots , d$ 上的晶格尺寸 $L_{\gamma}$ 都是偶数时，才严格保证不存在符号问题。绕线构型，如图50(b)所示的路径积分情形，如果绕行奇数个位点奇数次，则贡献为负。显然，在这种情况下，边界也会引入阻挫。因此，我们将只考虑（在所有周期晶格方向上）具有偶数尺寸的系统。

构型的链表存储。我们已讨论过将SSE构型存储为一个状态 $| \alpha ( 0 ) \rangle$ 和一个算符串 $S_{L}$。通过顺序作用算符，可以轻松生成如图55所示的所有传播状态 $\left| \alpha ( p ) \right.$。然而，在蒙特卡洛抽样过程中，我们还需要以非顺序的方式访问算符和传播状态的一些属性——给定一个算符及其作用的位点，我们需要知道下一个作用于这些位点的是哪些算符（无论是在算符序列中向上还是向下移动）。仅仅为获取这些信息而前后传播单个状态将极其耗时，并且存储所有传播状态也不现实。因此，我们还使用另一种数据结构，其中算符的“连通性”是显式的，并以紧凑的网络形式表示。这种链表顶点结构如图57所示，使用了与图55相同的SSE构型。在图形表示中，算符之间每个自旋的恒定状态被替换为直线，在计算机程序中这对应于链接（指针）。算符仅展示其作用前后的两个自旋状态。我们将这些自旋-算符对象称为顶点，并将这四个自旋称为顶点的腿。蒙特卡洛抽样过程的一个阶段将涉及对顶点进行修改（而它们的位点则在另一阶段改变）。链接将使我们能够在顶点之间快速移动，并进行一系列修改以维持所有约束条件。

允许的顶点由哈密顿量决定。以这里考虑的各向同性海森堡系统为例，共有四种顶点，如图58所示。它们当然对应于非零矩阵元（262），再次强调，对角算符中的常数 $\frac{1} {4}$ 是为什么没有四个自旋都相同的顶点的原因。各向异性相互作用或外磁场需要包含额外的顶点[190, 33]（并且此处讨论的算法也必须相应修改）。虽然四个腿的自旋状态可以唯一标识顶点，但为了更清晰，我们在图中仍继续使用空心和实心条，分别表示对角和非对角顶点。

图57. 图55中构型的链接顶点存储。在左侧的图形表示中，算符之间不变的自旋态已被替换为连接算符作用前后自旋的连线（链接）。这些链接可以存储在一个列表 $X ( \nu )$ 中，其中四个元素 $\nu = 4 p + l , l = 0 , 1 , 2 , 3$ ，对应于序列 $S_{L}$ 中位置 $p$ 处的顶点腿（编号约定如图58所示）。对于两个相连的腿 $\nu$ 和 $\nu^{\prime}$，有 $X ( \nu ) = \nu^{\prime}$ 且 $X ( \nu^{\prime} ) = \nu$。

图58. 各向同性 S = 1/2 海森堡模型允许的顶点。顶点腿的编号 l = 0,1,2,3 对应于图57所示链式存储中的位置 $\nu = 4 p + l$。

对于算符序列 $S_{L}$ 中的给定位置 $p$，对应的列表元素 $s ( p )$ 告诉我们算符类型（对角或非对角）以及其作用的键 $b$（如图55所示）。下面将看到，除了这些信息外，我们只需存储顶点的连接性，而不需要存储它们的自旋态。允许我们在相连的顶点腿之间跳转的链接存储为列表 $X ( \nu )$，如图57所示。为了图示清晰，这里的一维列表被排列成四列，元素标记为 $\nu = 4 p + l$，对应于每种类型的腿 $l = 0 , 1 , 2 , 3$，其标记方式在图58中指定。我们将在后面描述一种在给定算符序列的情况下构造此链表的高效方法。目前，只需知道对于序列中位置 $p$ 处的给定算符，其第 $l$ 条腿在链接顶点列表中的位置为 $\nu = 4 p + l$。该腿与另一个列表地址为 $\nu^{\prime} = X ( \nu )$ 的顶点腿相连。这种结构构成了一个双向链表，满足 $X [ X ( \nu ) ] = \nu$，我们可以在其中"向上"和"向下"移动。从列表中的位置 $\nu$，我们可以提取 $S_{L}$ 中对应的算符位置 $p = \nu / 4$（取整数部分）和腿索引 $l = \mathrm{mod} [ \nu , 4 ]$。我们还可以通过更改腿的标签 $01$ 或 $23$ 来"横向"移动到同一顶点的另一条腿，这在两种情况下都可以表示为列表地址 $\nu$ 的奇偶变化的简单规则。这些移动将使我们能够构建由改变的自旋和算符组成的闭合回路。

5.2.2 蒙特卡洛抽样过程#

当蒙特卡洛抽样配分函数 (258) 时，我们必须同时改变算符序列 $S_{L}$ 和存储的状态。然而，从图 55 和图 57 可以清楚地看出，这两个数据结构中的更新并非必然独立。位置 $p \in \{p_{1} , \ldots , p_{2} \}$ 处的算符所允许的更改（即保持所有构型约束的更改）会导致传播态在范围 $| \alpha ( p_{1} + 1 ) \rangle , . . . , | \alpha ( p_{2} ) \rangle$ 内发生相应变化。图 59 展示了这样一个两个算符的更新实例，其中仅算符的类型发生了改变，即从对角变为非对角（这导致的状态变化类似于图 51 中的世界线更新）。由于 $p_{1}$ 处改变的算符作用于 $| \alpha ( p_{1} ) \rangle$，因此状态 $\vert \alpha ( p_{1} + 1 ) \rangle$ 以及后续状态也随之改变（在本例中为两个自旋翻转），但当最后一个在 $p_{2}$ 处更新的算符作用后，这些修改被“修复”，最终的状态 $| \alpha ( p_{2} + 1 ) \rangle$ 与更新前相同。由于传播态构成一个循环置换的周期结构，存储的状态 $\mathinner{| {\alpha} \rangle} = \mathinner{| {\alpha ( 0 )} \rangle}$ 并无特殊之处，它也可能根据更新算符的位置而改变。

由于 $| \alpha ( p ) \rangle$ 具有时间（p）周期性，算符更新强制改变状态的方式并非唯一。在上述假设 $p_{2} > p_{1}$ 的例子中，我们也可以先用 $p_{2}$ 处的新算符作用于 $| \alpha ( p_{2} ) \rangle$。这时，状态 $| \alpha ( p_{2} + 1 ) \rangle , \dots , | \alpha ( L - 1 ) \rangle$ 以及 $| \alpha ( 0 ) \rangle , . . . , | \alpha ( p_{1} ) \rangle$ 都将被影响。总是有两种以这种方式改变状态的途径，但由于约束的存在，在任何特定情况下可能只有一种被允许，如图 59 所示。在大多数情况下，在随机位置进行的此类更新根本不能满足约束，因此必须专门寻找两个（或更多）可以更新的算符 [205]。

图 60. 对角更新。算符要么被插入，$[ 0 , 0 ]_{p} [ 1 , b ]_{p} ,$ 于随机选的键上，要么被移除，$[ 1 , b ]_{p} \to [ 0 , 0 ]_{p} .$。在前一种情况下，如果所选键上的自旋是平行的，则更新被取消，因此必须知道 $| \alpha ( p - 1 ) \rangle$ 的状态。接受概率取决于尝试更新前序列中的算符数量 n，如公式 (265) 所示。

如图 55 和图 57 所示，状态 $| \alpha \rangle$ 有时可以在不对算符序列做任何更改的情况下更新——因为没有任何算符作用于自旋 $\sigma ( 1 )$，其状态是任意的。在低温下，当每个自旋上的平均算符数量很大时，这种自由自旋非常罕见。

关注算符序列，我们将对非单位算符 $[ a , b ]_{p} \neq [ \bar{0 , 0} ]$ 的数量 n，它们在 $S_{L}$ 中的位置 $p \in \{0 , \ldots , L - 1 \}$，键索引 $b \in \{1 , \ldots , \bar{N}_{b} \}$ 以及类型指数 $a \in \{1 , 2 \}$ 进行抽样。一般策略是只允许通过填充单位算符与对角算符之间的交换来改变展开阶数 n；$[ 0 , 0 ]_{p} [ 1 , b ]_{p}$（其中 $n n \pm 1$）。然后，保持所有算符的格点位置 b 不变，对于某组算符（如上面两个算符的例子，但我们会用更高效的方式处理任意数量的算符）改变其类型指数，$[ 1 , b ]_{p} [ 2 , b ]_{p}$。如果这些非对角替换执行得当，那么结合对角 $n n \pm 1$ 更新，抽样过程是遍历的。接下来我们考虑这些更新过程的细节。我们首先假设展开截断 $L$ 已经被正确确定，然后讨论如何在实践中保证这一点。

对角更新。单个对角算符的更新如图60所示，可以在$S_{L}$中的位置$p = 0 , \ldots , L - 1$处顺序进行。在$[ 1 , b ]_{p} \to [ 0 , 0 ]_{p}$（删除一个对角算符）方向上进行此类更新时没有任何约束，而插入一个对角算符$[ 0 , 0 ]_{p} \to [ 1 , b ]_{p}$仅在键b上的自旋在其作用的传播态$| \alpha ( p ) \rangle$中反平行时才允许，即$\sigma ( i ( b ) ) \neq \sigma ( j ( b ) )$。因此，我们必须生成这些态（只存储当前所需的一个），这只需每次遇到非对角算符$[ 2 , b ]_{p}$时翻转两个自旋$\sigma ( i ( b ) )$和$\sigma ( j ( b ) )$即可（在这种情况下，不能在$p$处执行单算符更新）。

用于Metropolis接受概率的权重比很容易从(264)得到。但请注意，我们还需要校正这些更新尝试固有的不平衡——在$[ 1 , b ]_{p} \to [ 0 , 0 ]_{p}$方向上进行更新只有一种唯一方式，而当将$[ 0 , 0 ]_{p}$更改为$[ 1 , b ]_{p}$时，键b应从$N_{b}$个可能位置中随机生成[并且，如果对于选定的$b$，自旋$\sigma ( i ( b ) ) = \sigma ( j ( b ) )$，则立即拒绝该更新，我们继续处理下一个$p$]。因此，遵循第3.2节讨论的细致平衡准则，特别是将转移概率(43)写为选择概率和接受概率的乘积，我们必须包含涉及特定键b的更新的选择概率之比$N_{b}$。然后我们得到以下接受概率：

其中 $n$ 是更新前的算符数目（且在 $n \to n + 1$ 或 $n \to n - 1$ 之后）。

要证明对角更新的细致平衡，似乎有必要在随机位置 $p$ 而非顺序执行这些更新。否则，在某个位置 $p$ 执行更新后，下一步执行逆更新的概率为零。然而，如果考虑从 $p = 0$ 到 $P = L - 1$ 的整个更新序列（扫描），然后再执行从 $P = L - 1$ 到 $p = 0$ 的逆更新序列，则扫描满足细致平衡。在实际操作中，不必在算符序列的两个更新方向之间切换。

非对角更新。涉及非对角算符的更新显然必须至少涉及两个算符，以维持传播状态上的周期性约束。最简单的一类成对更新已在图59中讨论并说明。对于一维系统，此类更新在开边界条件下是遍历的，但在周期边界下，局域更新无法改变拓扑缠绕数（从图50可以清楚看到，局域形变不会影响循环置换）。原则上，这可能仍然不是一个非常严重的问题，因为对于任何 $T > 0$，当系统尺寸 $N \to \infty$ 时，只有零缠绕数扇区有贡献（实际上，对于较大的 $N$，温度越低 $N$ 需越大）。如果我们关注基态性质，在固定 $N$ 且 $T \to 0$ 时，也能得到正确结果，因为在此极限下，在传播状态 $\left| \alpha ( p ) \right.$ 的有限范围 $p , . . . , p + D_{p}$ 内定义的等效缠绕数（$D_{p} \gg N$，对应有限虚时范围）仍可自由涨落，不受全局零缠绕数约束的影响。物理量因此对缠绕数不敏感[206]。然而，为了收敛到基态，必须使用更低的温度。

除了在一维中无法改变缠绕数外，在高维中，图59所示的那类算符替换也无法导致应考虑的相同自旋的任何置换（玻色子语言中的粒子置换）（而在我们考虑的一维系统中，所有自旋的循环置换是唯一可能的置换）。为了解决此问题，可以构造其他类型的局域更新，例如，涉及作用于二维正方格子上 $2 \times 2$ 位点小方块角上的两个算符[205]，这种更新可以在给定缠绕数扇区内的所有此类置换之间采样。改变拓扑缠绕数的更新必须涉及环绕边界的一系列位点，但对于长度大于 $\approx 20$ 的格点，这些针对性的更新变得非常不可能（即，当同时改变大量算符时，满足约束变得过于困难）[186, 205]。

幸运的是，我们不必处理专门针对置换和缠绕数扇区的更新，有一种非常高效且简单的环更新（loop update）方法可以自动完成所有这些任务。这类更新最初是作为伊辛模型（Ising model）团簇算法（cluster algorithm）的推广而引入的，适用于被翻转的团簇呈环状形式的模型：经典六顶点模型（six-vertex model）[191]。利用离散铃木-特罗特分解（Suzuki-Trotter decomposition）构建的 $S = 1 / 2$ 海森堡模型的有效世界线系统，恰好等价于各向异性六顶点模型，因此针对它的环更新同时也将经典团簇更新推广到了量子力学系统。这些思想后来也被应用于连续时间世界线[179]以及 SSE 方法中的非对角更新[190]。相对于局部更新，性能提升是巨大的（如在经典情形中，这会大幅降低动力学指数），并将量子自旋系统的模拟带入了一个全新的水平。与经典团簇算法类似，环更新在实际中仅限于特定的几类模型，各向同性海森堡系统便是其中之一。将环的概念推广到虫（worms）[32]和有向环（directed loops）[33]（两者都可视为在构造过程中允许自交的环，而原始的环更新不允许自交），可以应用于更广泛的系统。

图61. 一个链接顶点（linked-vertex）SSE构型，其中画出了一个环及其两种“取向”，以及相应的算符索引序列。环所覆盖的所有自旋都被翻转，并且每次环经过时算符都在对角与非对角之间变化（如果一个算符被两次经过，则净变化为零）。每个顶点腿（自旋）唯一地属于一个环，而不被任何算符作用的自旋（此处为 $i = 1$ 处的自旋）也可视为形成了自己的环。

对于此处考虑的 $S = 1/2$ 模型，甚至没有必要详细讨论局部的非对角更新，我们将只专注于如何实现强大得多的环更新。在 SSE 方法中，算符串仍然是主要关注点，而环更新对应于构建一个由链接列表表示（linked-list representation）中的链接连接起来的算符（顶点）环。

算符环更新。 图61展示了一个算符环及其翻转方式的例子。这里的“翻转”指的是沿环的自旋（明确来说是顶点腿上的自旋，隐含地则是环所覆盖的所有传播态中的自旋），以及当我们将这些变化映射回序列 $S_{L}$ 时算符本身的变化；一个两条腿连接到环上的对角顶点会变成非对角顶点，反之亦然。当一个顶点的所有四条腿都属于同一个环时（图中便有一个这样的例子），算符类型不变。

环更新的关键在于，权重 (264) 只依赖于算符的数量，而算符环翻转时这个数量不变。因此，根据细致平衡规则，这样的翻转总是可以接受的。同样清楚的是，环可以很大，因此环翻转可能导致全局自旋置换的变化。前面讨论的成对替换（pair substitutions）对应于翻转同一键上两个算符的环。任何专门为采样置换而构造的算符替换，也可以表述为算符环，包括那些改变缠绕数的替换。

请注意，一旦算符位置确定，环路就是完全确定性的，并且全时空构型中的每个自旋都唯一地属于一个环路。对角更新是实际改变环路的机制，而环路更新的目的仅仅是识别并翻转其中一些环路。因此，最好的做法是先构建所有环路，然后以概率 $1 / 2$ 翻转每一个（这也保持了细致平衡），而不是随机构建环路并总是翻转它们（在那种情况下，某些环路会被构建多次并被不必要地翻转）。自由自旋（即那些没有顶点作用的自旋，如图61中的第一个自旋）也可以被视为环路，因为翻转这样一个自旋隐式地翻转了时间周期性传播态全构型中的一整条自旋线。

5.2.3 计算机实现#

在SSE方法中，蒙特卡洛扫描的一个方便定义是：先进行完整的对角更新序列，然后构建所有环路并（以概率 $1 / 2$）翻转。从图61的环路示意图中可以清楚看出，SSE构型的链接顶点存储 $X ( )$ 为何是有用的；构建环路相当于在这个列表中根据非常简单的规则移动（沿着链接移动，或横向移动到同一顶点上的相邻自旋）。另一方面，在对角更新中，使用单个自旋状态 $\sigma ( )$ 和算符串 $s ( )$ 的原始存储（如图55所示）更为方便。构型将始终以这种方式存储，而链接顶点将在每组环路更新之前构建。环路更新通过在链接列表 $X ( )$ 中移动来执行，同时更新 $s ( )$ 中相应的算符（如图57所示）。存储的自旋 $\sigma ( )$ 也可能受到环路更新的影响，这可以在所有环路翻转完成后处理。现在，我们讨论蒙特卡洛扫描中所有主要步骤的伪代码实现。

对角更新。我们假设当前存储了一个有效构型（模拟开始时，可以是随机自旋状态和仅包含 $[ 0 , 0 ]$ 元素的算符序列），使用 $\sigma ( i ) = \pm 1$ 表示单个自旋，并将 $[ a , b ]_{p}$ 紧凑存储为单个整数 $s ( p )$ ，其中偶数和奇数 2b 和 $2 b + 1$ 分别对应对角算符和非对角算符（如图55所示）。然后，可以通过以下方式实现一次对角更新扫描：

do $p = 0$ to $L - 1$
if (s(p) = 0) then
b = random[1, . . ., N_b]; if (i(b)) = ( j(b)) skip to next p
if (random[0 1] < P_insert(n)) then s(p) = 2b; n = n + 1 endif
elseif (mod[s(p), 2] = 0) then
if (random[0 - 1] < P_remove(n)) then s(p) = 0; n = n - 1 endif
else
b = s(p)/2; (i(b)) = - (i(b)); ( j(b)) = - ( j(b))
endif
enddo

{27}

这里 $P_{\mathrm{insert}} ( n )$ 和 $P_{\mathrm{remove}} ( n )$ 分别是当前展开阶为 n 时，插入和移除一个对角算符的接受概率（265）和（266）。这些概率可以预先计算或在运行时计算。也可以仅预先计算与 n 无关的常数 $\beta N_{b} / 2$ ，如果 random $[ 0 - 1 ] \times ( L - n ) < \beta N_{b} / 2$ 则接受移除，如果 random $[ 0 - 1 ] \times \mathsf{\bar{\beta}} N_{b} / 2 < L - n + 1$ 则接受插入。这避免了使用完整接受概率（265）和（266）所需的更耗时的除法，或存储大量依赖于 n 的概率。

属于键 b 的格点可以存储在列表 $[ i ( b ) , j ( b ) ]$ 中。除了这个列表，SSE采样算法是与晶格无关的。也很容易修改上述代码以适应非均匀的（依赖于 b 的）海森堡耦合强度。对角更新是唯一显式出现相互作用的阶段，在环路更新中完全不需要任何修改。

构建链接顶点列表。我们将从$p = 0$开始，遍历算符$s ( p )$来构建该链表。回想一下，对于$p$处算符的顶点支腿编号$l = 0 , 1 , 2 , 3$，相应的链接存储为$X ( \nu )$，其中$\nu = 4 p + l$，如图57和图58所示。为了构建该列表，我们将利用另外两个数据结构。自旋$i$的第一个顶点支腿在$X ( )$中的位置$v$将存储为$V_{\mathrm{first}} ( i )$，而$V_{\mathrm{last}} ( \ r_{i} )$将是迄今为止遇到的自旋$i$上最后一个支腿的位置。注意，给定自旋上的第一个支腿总是在算符下方（在算符作用之前，$l = 0 , 1 )$，而最后一个支腿则在算符上方$( l = 2 , 3 )$。初始时，所有$V_{\mathrm{first}} ( i )$和$V_{\mathrm{last}} ( i )$都设置为-1。之后，当遇到作用于自旋$i$的算符$s ( p )$且$V_{\mathrm{last}} ( i ) \geq 0$时，$V_{\mathrm{last}} ( i )$就是对应自旋$i$上一个操作的列表元素。然后我们使用新算符的相应下部支腿（$l = 0$或1，取决于所作用两个自旋中的哪一个等于$i$），并设置$\nu = 4 p + l$与$V_{\mathrm{last}} ( i )$之间的链接：$X ( \nu ) = \bar{V_{\mathrm{last}} ( i )}$ 和 $X ( V_{\mathrm{last}} ( i ) ) = \nu$。另一方面，如果$V_{\mathrm{last}} ( i ) = - 1$，那么我们找到了自旋$i$的第一个操作，并设置$V_{\mathrm{first}} ( i ) = \nu$。在这两种情况下，我们都将自旋$i$上最后一个操作的元素设置为$V_{\mathrm{last}} ( i ) = \nu + 2$，这里加2对应于操作后自旋$i$上的顶点支腿（支腿索引$l = 2 \mathrm{or} 3 )$。对于每个算符，我们必须检查它所作用的两个自旋$i_{1}$和$i_{2}$。在下面的伪代码中，对于$p$处的算符，要填充的顶点列表位置是$\nu_{0} + l$，其中$\nu_{0} = 4 p$，而支腿对$l = 0 , 2$和$l = 1 , 3$分别对应于自旋$i_{1}$和$i_{2}$。我们用$\nu_{1}$和$\nu_{2}$来表示两个自旋上最后一次（前一个）操作的列表元素。除了跨越边界$p = L - 1 , 0$的链接外，所有链接都可以按照下式构建：

为了连接跨越时间边界的链接，我们可以使用包含第一个和最后一个列表位置的数组：

由于 $X ( )$ 中的某些位置也不包含任何链接，对应于填充算符 $s ( p ) = 0$，我们需要某种方式来区分它们。我们可以在构建链接之前将所有 $X ( \nu )$ 设为负数，例如 $X ( \nu ) = - 1$。这样，任何满足 $X ( \nu ) \geq 0$ 的 $\nu$ 都可以用作追踪环路的起点。

算符环路更新的实现。我们要追踪所有环路，并以 $1 / 2$ 的概率翻转每个环路。为了确保每个环路不被重复考虑，我们应该总是从 $X ( )$ 中一个之前未被访问过的位置开始新环路。然后我们需要某种标志来指示某个位置是否已被访问。与其为这些标志分配单独的存储空间，我们实际上可以复用列表 $X ( )$ 本身——其中存储的链接只使用一次，一旦 v 处的链接被使用，我们就可以将 $X ( \nu )$ 设为负数，以指示该位置不应再次使用。为了开始一个新环路，我们可以寻找第一个满足 $X ( \nu_{0} ) \geq 0$ 的位置 $\nu_{0}$。在遍历环路之前，我们随机决定是否翻转它。如果不翻转，则仍需执行环路追踪，以将所有相应的 $X ( \nu )$ 标记为已访问。对于被翻转的环路，我们还应该更改 $s ( )$ 中的算符，以反映每个被访问顶点处对角到非对角的改变。之后我们再处理存储态 $| \alpha \rangle$ 中自由自旋的可能翻转。用于追踪和翻转环路的代码结构如下：

do $\nu_{0} = 0$ to $4 L - 1$ step 2

if $( X ( \nu_{0} ) < 0 )$ 跳到下一个 v_0

{30}

if (random[0 - 1] < 1/2) then

遍历该环路；对环路中所有 v，设 $X ( \nu ) = - 1$

else

遍历该环路；对环路中所有 v，设 $X ( \nu ) = - 2$

• 翻转该环路（在遍历环路时更改算符类型 $s ( p = \nu / 4 )$）

endif enddo

这里 $\nu_{0}$ 是一个尝试性的新环路起点，只有当 $X ( \nu_{0} )$ 是非负数时，我们才真正开始一个新环路。注意，由于顶点列表和环路的结构，我们只需要考虑偶数起点（$\{30 \}$ 第一行中的步长 2），即 $X ( \nu )$ 中两个连续的（偶数，奇数）元素总是对应于同一个环路。为了不重复访问同一个环路，我们将仅被访问的环路设为 $X ( \nu ) = - 1$，将被翻转的环路设为 $X ( \nu ) = - 2$。我们稍后需要区分被翻转的环和仅访问的环路，以便更新存储态来反映被翻转的环路。

为了从某个 v_0 开始追踪一个环路，我们可以使用支腿索引 $l_{0} = \mathbf{mod} [ \nu_{0} , 4 ]$，并移动到同一个顶点上的相邻支腿 $l_{0}^{\prime}$。根据图 58 中的支腿标记约定，与支腿 l 相邻的支腿的规则可以概括为 $l = ( \bar{0} , 1 , 2 , 3 ) l^{\prime} = ( \bar{1} , 0 , 3 , 2 )$ 该规则可以通过对 l 进行位操作非常高效地实现，因为 $l^{\prime}$ 对应于翻转 $( 01 )$ l 的最低有效位。这里我们假设所用编程语言中具备此功能 [例如，在 Fortran 90 中，可以对整数 l 使用异或操作；$l^{\prime} = \mathbf{ieor} ( l , 1 ) ]$，并设 $l_{0}^{\prime} = \mathbf{fipbit} ( l_{0} , 0 )$，其中 0 是被翻转的位。我们实际上不需要提取支腿索引本身，而可以通过操纵其位置 $\nu_{0} = 4 p_{0} + l_{0}$ $X ( )$ 来进行到相邻支腿的相应移动。我们可以使用相同的位翻转方法来找到它；$\nu_{0}^{\prime} = \mathbf{flipbit} ( \nu_{0} , 0 )$。为了移动到下一个顶点，我们使用链接；$\nu_{1} = X ( \nu_{0}^{\prime} )$。这完成了环路追踪过程的一个步骤。接下来我们以同样的方式找到相邻支腿的位置 $\nu_{1}^{\prime}$，并从那里移动到与之相连的支腿 $\nu_{2}$。如此继续，直到在某一步 $k$ 我们找到 $\nu_{k} = \nu_{0}$。至此，我们完成了一个完整的环路。以伪代码形式，环路的遍历及其翻转可以通过以下方式完成（注意被访问的列表位置 v 不需要下标或上标）：

这里，对算符类型（对角或非对角）的翻转也通过 flipbit 过程完成，这恰好对应于将算符代码 $s ( p )$ 从 2b（对角算符）改为 2b + 1（非对角算符），反之亦然。在代码 30 中不翻转环路的情况下，唯一的区别是 31 中涉及 p 和 $s ( p )$ 的两条语句被省略，并且我们将所有被访问的位置标记为 $X ( \nu ) - 1$ 而不是 2（这些信息将用于更新状态 $| \alpha \rangle$）。

在所有环路追踪完成后，我们还需要更新存储状态 $| \alpha \rangle$ 中的自旋。对于每个有翻转环路通过的格点 i，存储的自旋 $\sigma ( i )$ 应当被翻转。我们还应以 $1 / 2$ 的概率翻转任何没有环路连接到其上的自旋（即没有算符作用的自旋）。这两个任务都可以在构建链接顶点列表时利用存储在 $V_{\mathrm{first}} ( i )$ 中的信息来处理。如果 $\nu = V_{\mathrm{first}} ( i ) = - 1$，这表示没有算符作用于自旋 i。否则，我们可以利用现在存储在 $X ( \nu )$ 中的标志 1 或 2 来确定通过 i 的环路是否被翻转 [对于翻转的环路，$X ( \nu ) = - 2$]。伪代码形式如下：

do i = 1 to N v = V_first(i) if (v = 1) then if (random[0-1]< 1/2) (i) = (i) else if (X (v) = 2) (i) = (i) endif enddo

{32}

这完成了环路更新，我们可以重复以下周期：(i) 执行一次对角更新的扫描，(ii) 构建链接顶点列表，(iii) 追踪/翻转所有环路并相应更新。这些步骤构成了 SSE 方法中的一次蒙特卡洛扫描。

平衡化与展开截断。 我们现在讨论展开截断长度 L（当使用固定长度方案时）。我们应该确保 L 超过将被采样的最大展开阶数 n，这样实际上的截断就不是一种近似。由于新的对角算符通过与填充的单位算符交换而插入，我们也希望有相当数量的这些 $s ( p ) = 0$ 元素存在。但由于内存和 CPU 时间与 $L$ 成比例，我们不应使 L 过大。原则上可以定义一个最优的 $L$，使得对角更新的接受率最大化。然而，实际上 L 是否真正最优并不重要，因为在决定自相关时间上，非对角更新更为重要。因此，我们只需确保 $s ( p ) = 0$ 的元素占据一个显著的比例，使得采样的 n 永远不会达到 L，并且在对角更新中可以尝试大量的算符插入。

我们可以设定 $L \approx ( 1 + a ) \langle n \rangle$，例如取 $a = 1 / 3$。由于 n 的相对涨落较小，约为 $\sim \sqrt{n}$，我们可以通过以下方式实现：令 $L \to a \times n$，并在每次蒙特卡洛扫描后，当 $a \times n > L$ 时，向 $S_{L}$ 中添加额外的 $s ( p ) = 0$ 元素 [对于 $L_{\mathrm{old}} + 1 < p < L_{\mathrm{new}}$ 设置 $s ( p ) = 0$，或者随机分布这些零]。图 6.2 展示了一个在逆温度 $\beta = 16$ 下 $16 \times 16$ 系统平衡化过程中该过程实际运行的例子。初始时 L 设为 $N / 2 = 128$。L 的调整迅速收敛到一个可接受的值，并且如插图所示，在后续的长时间模拟中，n 从未接近 L。

5.2.4 物理观测量的估计量#

在采用环更新的SSE模拟中，自相关时间通常非常短；通常仅需几个蒙特卡洛扫描步或更少[33, 31]。因此，物理观测量（“测量”）的计算通常应在每次（或每隔几次）蒙特卡洛扫描后执行。如第3.2节所述，自相关时间并不决定测量频率，仅决定独立数据的生成频率。因此，只要总模拟时间远大于自相关时间，自相关会影响误差棒，但不会影响计算期望值的正确性（即测量相关数据并非错误，但若测量耗时较长，可能造成浪费）。为可靠计算误差棒，应采用数据分组方法，这在第3.2节也有讨论。本节我们将讨论SSE计算中几种重要的观测量类型。更详细的内容和表达式推导可参考文献[186, 187, 207]。

图62给出了一个$16 \times 16$海森堡系统在逆温度$\beta = 16$的平衡过程中，展开截断L初期的演化。图中同时显示了每次蒙特卡洛扫描后弦中算符数n的变化，以及迄今为止达到的最大n值。经过5000次扫描后的最终截断值为$L = 6764$。插图显示了后续$10^{7}$次蒙特卡洛扫描过程中n的分布。

能量与比热。我们在第5.1.3节已讨论了一些观测量：内能(250)、哈密顿量中单个算符的相关期望值(251)，以及比热(252)。已知$\langle H \rangle = - \langle n \rangle / \beta$，对比热$C = ( \langle H^{2} \rangle - \langle H \rangle^{2} ) / T^{2}$并不仅仅由n的涨落给出，还包含一项n这一结果可能令人惊讶。这是因为$\langle H^{2} \rangle = \langle n ( n - 1 ) \rangle / \beta^{2}$，这一点可用第5.1.3节推导$\langle H \rangle = - \langle n \rangle / \beta$ 的相同方法轻松证明。注意该表达式包含了从每个键算符中减去的常数$1 / 4$（见(256)式），但该常数在比热表达式(252)中抵消了。在实际计算中，低温下比热（此时数值很小）难以可靠计算（即统计误差很大），因为它是由两个大数（$\sim N^{2} \beta^{2}$）之差得到的。

对角算符。在z分量基矢下对角算符的期望值也易于计算，只需对传播态求平均即可：

无需对所有$p$求此平均，因为相邻的传播态最多只相差两个翻转的自旋。为了节省时间，可以对$p$进行部分求和，例如$O_{z} ( 0 ) + O_{z} ( N ) + O_{z} ( 2 N ) + \ldots$。然而，在某些情况下，也可以直接从$p = 0$到$p = L - 1$传播$O_{z} ( p )$：先计算出$O_{z} ( 0 )$，然后根据传播$| \alpha ( p ) \rangle$时发生的自旋翻转来更新$O_{z} ( p )$的值。例如，从一个SSE构型计算交错磁化强度平方的期望值的全平均贡献，可按如下方式实现：

这里 $\phi_{i} = \pm 1$ 是格点 i 的交错相位因子，m 包含交错磁化强度，它随着遍历算符字符串而变化。当翻转两个反平行自旋时，m 的变化可以表示为 $2 \phi_{i} \sigma ( i ( b ) )$，因为 $\phi_{j} \sigma ( j ( b ) ) =$ $\phi_{i} \sigma ( i ( b ) )$。m 平方的平均值累积为 $m_{s 2}$，在此示例中仅使用由原始索引序列生成的 n 个传播态，而不包含填补单位算符 [即，跳过算符 $s ( p ) = 0 ]$ 的步骤]。也可以对所有 L 个状态实例求和，此时结果将在最后一行除以 L 而不是 n。如上所述，特殊情况 $n = 0$ 没有得到正确处理，但代码很容易修改，只需在这种情况下使用 $\dot{m_{s 2}} = m^{2} / N^{2}$ 即可。

33 中测量的计算量缩放为 $\beta N_{;}$，即与 SSE 采样算法相同。若要计算任意动量下自旋关联（静态结构因子）的傅里叶变换 [注意 $m_{s}^{2}$ 对应 $S ( \pi )$，其中 $\pi$ 表示交错波矢，例如二维中的 $\pi = ( \pi , \pi )$]，例如用于关联长度定义 (70) 或 (73)，可以对 $m ( \mathbf{q} )$ 的实部和虚部使用类似过程

其中 33 中的 $\phi_{i}$ 被相应的正弦和余弦因子取代 [并明确考虑来自 $i ( b )$ 和 $j ( b )$ 处翻转自旋两者的变化]。然后根据 $S ( \bar{\bf q} ) = S ( \bar{\bf q} ) + \mathrm{Re}^{2} \{m ( {\bf q} ) \} + \mathrm{Im}^{2} \{m ( {\bf q} ) \}$ 累积结构因子

为了计算完整关联函数，

对所有 $i , j$ 求平均值 [或对所有 q 的傅里叶变换，可写为 $\langle m ( {\bf q} ) m ( - {\bf q} ) \rangle ]$，像 33 中那样在每次自旋翻转后更新所有这些 $\propto N$ 个不同函数是不切实际的。更好的方法是，例如在每个第 N 个传播态从头计算所有关联。这仍然导致测量所有关联的计算量缩放为 ∝ $\beta N^{2}$——比采样差一个因子 N。那么应该判断是只计算关联的一个子集 [例如，沿 $( x , y )$ 平面中的某些线] 更好，还是仅仅减少测量频率，以便测量不会完全主导计算。

（磁化率）。另一类重要的观测量是广义磁化率，即线性响应函数，形式为

其中 b 是添加到哈密顿量中的场项 bB 的前因子，A 是我们想要计算其对这种微扰响应的算符。例如，我们可能对磁场仅作用在格点 i 时格点 j 处的响应 $\chi_{i j}$ 感兴趣，此时 $A = S_{j}^{z}$ 且 $B = S_{i}^{z}$。这种磁化率由库伯公式（Kubo formula）[186] 给出

其中 $A ( \tau ) = \mathsf{e}^{- \tau H} A \mathsf{e}^{\tau H}$。如果 A 和 B 都是对角的，则此库伯积分可以在 SSE 模拟中使用通用公式计算

对 $A(p)$ 和 $B(p)$ 的求和可使用与代码33类似的代码计算。

我们通常关注 $A = B$ 的磁化率，例如波矢q处的磁响应 $\chi(\mathbf{q})$，此时 $A = B = m(\mathbf{q})$，其中 $m(\mathbf{q})$ 为自旋构型的傅里叶变换，见(268)。由于 $\chi(\mathbf{q})$ 的精确值必为实数，只需计算(272)的实部。若模拟中计算了所有局域响应函数 $\chi_{i j} = \chi(\mathbf{r}_{i j})$，则可利用它们随后计算傅里叶变换 $\chi(\mathbf{q})$。这些响应函数在核磁共振实验中可直接获取；示例见文献[208]。

注意，若A和B与哈密顿量对易，则 $A(\tau) = A(0) = A$ [且 $A(p) = A$ 与 $B(p) = B$ 在(272)中与p无关]，磁化率退化为经典表达式 $\chi_{AB} = \beta \big( \langle AB \rangle - \langle A \rangle \langle B \rangle \big)$。实际中，对于海森堡模型，唯一适用该形式的情形是均匀磁化率：

此处(271)中的 $\beta \langle M_z \rangle^2 $ 项消失，因为哈密顿量H不含磁场。

若(271)中的算符A和B非对角，磁化率通常更为复杂。此处仅讨论一个重要特例，其中估计量实际上非常简单。若涉及的算符为哈密顿量的两项（如(255)和(256)所定义，此处简称为 $H_A$ 和 $H_B$，对应任意两个此类算符——任意两个对角或非对角键算符），则磁化率测量仅涉及对采样的SSE算符序列中这些算符的数量 $N(A)$ 和 $N(B)$ 进行计数：

此类最重要的例子是电流磁化率 $\chi_{I_x I_x}$，其中自旋电流算符 $I_x$（此处为明确起见，指晶格x方向）定义为：

为简化起见，假设哈密顿量仅含最近邻相互作用。对于长程相互作用，哈密顿量中相同位点对之间将存在对应的电流项。虽然(275)并非SSE采样中使用的非对角算符(256)的直接求和，但这些算符可写为两部分之和：

其中每个情况下的位点对 $[i(b), j(b)]$ 的排序使得+部分和-部分分别沿键方向 $\mathbf{r}_j - \mathbf{r}_i$ 的正向和负向传输一个单位的自旋。这里的关键点在于，尽管SSE构型包含完整的非对角算符 $H_{2,b}$，但在每个实例中仅有(276)中的+或-部分有贡献（另一部分会破坏其所作用的态）。因此可利用(274)评估如下形式的局域电流-电流响应函数：

其中 $I_b$ 表示键b处的电流算符。利用(274)可得：

其中 $N^{+} ( b )$ 和 $N^{-} ( b )$ 分别表示 SSE 算符序列中通过键 $b$ 沿正方向和负方向传输自旋的算符数目。当查看 SSE 构型的图形表示（例如图 55）时，这一含义就变得清晰了。在该图中，只需对自旋向右移动（正方向）的非对角事件数减去向左移动的事件数进行求和。除非由键 $b$ 定义的晶格方向上的缠绕数不为零，否则 $N^{+} ( b ) - N^{-} ( b )$ 项为零，这可以从图 50 中缠绕数的图示看出。值得注意的是，电流响应的形式 (278) 在 SSE 和路径积分（世界线）方法中是相同的。

自旋刚度。 我们在第 3.5 节讨论经典 XY 模型时曾涉及自旋刚度。该量在 $T = 0$ 时的基本定义（式 (83)）以及在 $T > 0$ 时的定义（式 (86)）对于量子自旋系统（XY 模型或海森堡模型）以及更广泛的玻色子系统（其中类似量为超流相位刚度，与超流密度成正比）是相同的。为 SSE 或世界线方法推导蒙特卡洛估计量时发现，式 (97) 中的项 $\beta \langle I_{x}^{2} \rangle$（一个经典响应函数）应替换为相应的量子力学久保公式，从而自旋刚度表示为

这里假定相位扭曲沿晶格 $x$ 方向施加。其中的久保积分由形式为 (277) 的键电流响应函数求和组成，利用结果 (278) 可以发现能量项 $\left. H_{x} \right.$ 恰好被来自 (278) 的 $\delta$ 函数项抵消。这导致了用 SSE（或世界线）构型表示的自旋刚度的非常简单的表达式：

这里 $N_{x}^{+}$ 和 $N_{x}^{-}$ 分别表示沿 $x$ 正方向和负方向传输自旋的算符总数。对于二维或三维的空间各向同性系统，这当然可以在所有维度上取平均值，而对于各向异性系统，刚度则依赖于方向。

自旋刚度通常表示为缠绕数的涨落（此处为 $x$ 方向，其他方向表达式类似），

The winding number is an integer (positive, negative, or zero). The spin stiffness is therefore often written as

即，尺寸归一化电流 $W_{x}=0,\pm1,\pm2,...$。对于一个 $d$ 维系统，若 $N = L^{d}$，则任意方向 $a$ 的刚度即为 $\rho_{s} = \langle W_{a}^{2} \rangle / \beta$，但对于其他形状（例如 $N = L_{x} \times L_{y}$，且 $L_{x} \neq L_{y}$），长度也会进入公式，此时 (280) 式更为简洁。

在推导自旋刚度时，通常为简单起见，假设自旋在（自旋空间中的）XY 平面内有序；然而在有限系统的 SSE 模拟中，对于海森堡模型而言，序参量在自旋空间中的方向并不受限于该平面（对称性未被破坏）。一个简单的旋转平均论证表明，为了得到对称性破缺态（即在考虑热力学极限时）的正确刚度，应将表达式 (280) 乘以 $3/2$。表达式 (280) 也适用于量子 XY 模型[209]，在该模型中结果不应乘以 $3/2$。

一个与 (280) 完全类似的缠绕数估计量最初是在连续空间中玻色子的超流相刚度的背景下推导出来的 [210]。全局循环置换粒子在超流中的作用可追溯到费曼在路径积分方面的开创性工作 [175]。

图 63. 由穿过传播态 $| \alpha ( p ) \rangle$ 的时空环形成的团簇示例（对于任意固定的 $p$）。这里有六个这样的团簇，标记为 $1, \ldots , 6$。空心圆和实心圆分别对应于 $| \alpha ( p ) \rangle$ 中的自旋。当一个环被翻转时，对应团簇中的所有自旋也会被翻转，如这里所示的两个不同构型，对应团簇 1 的两种状态（由大圆圈包围的格点）。注意，每个团簇内的自旋总是处于两种交错构型之一。

5.2.5 改进估计量#

SSE 方法中的算符环更新（以及更广泛的环更新 [31, 191]）是团簇更新的一个例子。这种非局域更新最初是为伊辛模型的经典蒙特卡洛模拟开发的 [118]。团簇方法的一个方面在于，可以对物理量的估计量在所有团簇取向取平均值，因为团簇翻转后构型权重不变。这在 $S = 1/2$ 海森堡模型的 SSE 模拟中非常明显，因为权重 (264) 仅取决于序列中的算符数量 $n$，而环翻转时该数量不变。如果团簇（此处为算符环）的数量为 $m$，那么等权构型的总数为 $2^{m}$，对这些所有构型取平均值所能得到的估计量，其噪声远小于仅依赖单个构型的估计量。这里的关键在于，对于许多重要的物理量，这个平均值可以解析计算，并且得到的改进估计量形式简单，可以在模拟中快速评估。这里我们只讨论静态（等时）结构因子和均匀磁化率这两个相对简单的情况。关于其他物理量的改进估计量，参见 Evertz 的综述文章 [31]。

考虑一个传播态 $| \alpha ( p ) \rangle$，例如存储的 $| \alpha ( 0 ) \rangle$。在图 61 所示的 SSE 构型的链接顶点表示中，有一个环穿过该态中的每个自旋（没有算符作用的自旋也视为单独的环）。同一个环可以穿过 $| \alpha ( 0 ) \rangle$ 中的许多自旋，并且所有属于同一个环的自旋形成一个团簇，其意义在于：如果该环被翻转，则团簇中所有自旋同时翻转。注意，环是时空中的对象，而这里讨论的团簇是在固定时间（此处为传播指数 $p$）的截面上定义的。一个团簇可以由在空间上看似不连续的几个部分组成，因为这些部分可以在环存在的更大时空体积内相连接。图 63 展示了一个二维晶格上团簇的例子。

由于我们处理的是二分格子，并且环的构造方式使得在沿着环移动时，每次改变方向，垂直环段上的自旋（参考图61等图示）都会发生变化，因此在给定态 $| \alpha ( p ) \rangle$ 下形成的团簇内的自旋总是具有交错结构。标记为 $j$ 的团簇的交错磁化强度 $m_{s} ( j )$（其中 $j = 1 , \ldots , C$，C 是团簇总数）为 $m_{s} ( j ) = \pm n_{j}$，其中 $n_{j}$ 是团簇 $j$ 中的自旋数。对于给定的构型，总交错磁化强度为 $\begin{array} {r} {M_{s} = \sum_{i = 1}^{C} m_{s} ( j )} \end{array}$。当对所有不同的团簇取向实现平均这个和的平方时，交叉项 $\langle m_{s} ( i ) m_{s} ( j ) \rangle = 0 ( \mathrm{for} i \neq j )$。因此只剩下 $i = j$ 的贡献，交错结构因子简单给出为

其他波矢 q 处的结构因子仅略微复杂，需要将团簇大小 $n_{j}$ 替换为对每个团簇的相位 $\phi_{\mathbf{r}} \mathrm{exp} ( i \mathbf{q} \cdot \mathbf{r} )$ 求和，其中 r 指给定团簇上的格点，$\phi_{\mathbf{r}} = \pm 1$ 是交错相位因子，负责处理团簇内部的自旋交错结构（分母 4 用于修正自旋值为 $\pm 1 / 2$ 的事实）。这里也可以利用真实结构因子对于任何 q 都必须是实数的性质。

原则上，如式 (267) 所示，等时关联函数（例如结构因子）也可以（完全或部分地）对传播指数 $p$ 取平均。然而，与代码33中简单估计量的情形相比，这需要为改进估计量做更多工作，因为在环更新过程中需要记录某些信息，以便为几个固定的 $p$ 构建团簇，而且这样做可能并不总是值得。然而，如果没有这种平均，一个简单的 p 平均估计量（例如代码33中实现的）实际上可能在低温下给出更好的结果，因为平均带来的增益可能非常显著。在这方面，$\mathbf q = 0$ 的情况是特殊的，因为这对应于总磁化强度平方，它是一个守恒量（即独立于 SSE 传播指数 $p$），因此无法通过对 $p$ 进一步平均来改进统计。因此，均匀磁化率 (273) 的最优估计量为

其他波矢的磁化率涉及完整的时空环结构，而不仅仅是固定状态下形成的团簇（环的截面）。例如，交错磁化率由所有环尺寸的平方和给出 [31]。

5.2.6 程序验证#

QMC 程序应始终通过将小系统的结果与精确对角化（exact diagonalization）数据进行比较来验证。当正确实现时，SSE 方法应为精确的，这意味着计算量与其真实值的偏差应纯粹是统计性的，源于有限数量的采样构型。我们在第 3.2 节讨论了如何根据“误差棒”量化统计涨落。超出误差棒的偏差要么源于编程错误，要么源于所使用的随机数生成器的缺陷。虽然大多数编程错误会导致明显错误的结果，但也可能存在细微错误，这些小系统上这些错误仅会导致与精确结果有微小偏差。此类系统误差可能随系统尺寸增长而产生严重后果。随机数生成器也是如此；不完美随机数导致的偏差影响在小晶格上可能非常微小，但在较大系统中可能变得更加显著。因此，基于非常长的 SSE 运行与精确对角化结果进行比较以检测可能存在的微小偏差至关重要。虽然不可能严格证明程序在所有方面都是正确的，但与精确数据在非常小的误差棒内一致使这一点变得非常可能。

图64. 一个 $4 \times 4$ 海森堡系统的均匀磁化率随温度的变化。曲线是完整对角化的精确结果。点是 SSE 结果，基于每个温度 $10^{10}$ 次更新扫描。在此尺度下无法分辨误差棒。插图显示了 $\chi$ 与精确结果的偏差（带误差棒），包括简单估计量和改进（环）估计量。为了使两个估计量的误差棒不重叠，数据点相对于其实际 T 值略有偏移。

测试结果。我们现在讨论 1D 和 2D 海森堡模型的一些测试结果。这些计算（以及本讲义中讨论的大多数其他 SSE 计算）中使用的随机数生成器是一种简单的 64 位线性同余生成器，乘数为 2862933555777941757 [211]，加法为 1013904243。

图 64 显示了一个 $4 \times 4$ 系统的每自旋磁化率。对于每个温度，执行了 $10^{10}$ 次更新扫描和测量，并将数据细分为 100 个组（bin）以计算统计误差。对于 $T = J / 10$，在一台中档 PC 工作站上这大约需要 100 CPU 小时。当使用 $\chi$ 的标准估计量时，误差棒（如第 3.2 节所讨论，定义为估计平均值波动的标准差）通常约为 $\approx 10^{- 6}$，而对于改进估计量 (283) 则更小。对于 $T / J \approx 1$，相对误差（误差棒除以 $\chi$）约为 $5 \times 10^{- 6}$。平均值的偏差与误差棒的大小完全一致。回想一下，对于正确计算的统计误差，大约应有 $2/3$ 的数据点将其真值包含在一个误差棒内。对于更大的系统，使用改进估计量获得的精度提升可能比图 64 中大得多。

虽然在小格点上确认 SSE 计算的无偏性至关重要，但对于更大系统，人们仍可能担心随机数生成器的潜在问题。在一维情况下，我们还可以将 SSE 计算结果与基于精确 Bethe ansatz（贝特假说）的非常长链基态结果进行对比。为了在 SSE 计算中逼近基态，方便的做法是使用 $\beta = 2^{m}$ 形式的逆温度，并增大 $m$ 直到计算量变得与 $m$ 无关。图 65 展示了这一方法，其中绘出了长度为 $N = 1024$ 和 $4096$ 的链，使用简单展开阶估计量 (250) 计算得到的内能 $E = \langle H \rangle$。这些链长的数值 Bethe ansatz 结果列于文献 [212] 中。两者在误差棒内完美符合，且对于每个 $N$，最后三到四个点未检测到温度依赖性。这种良好符合，在相对统计误差低至 $5 \times \mathrm{10^{- 7}}$ 的范围内，相当令人信服地表明，该计算实际上完全无偏，并且可以研究甚至相当大的系统的基态。

图 65. Heisenberg 链在长度 $N = 1024$ 和 $4096$ 时，多个逆温度 $\beta = 2^{m}$ 下每个自旋的内能的 SSE 结果，与相应的精确（Bethe ansatz）基态能量的比较。SSE 结果通过每个 $m$ 进行 $1 - 2 \times 10^{6}$ 次更新扫描获得，对于最大的 $N , \beta$ 需要数百 CPU 小时。

5.3 SSE 在一维和二维系统中的应用#

接下来我们讨论使用 SSE 方法得到的一维和二维 Heisenberg 模型的几个具有说明性的结果。在 5.3.1 节中，给出了清晰显示一维链临界行为对数修正的结果。在 5.3.2 节中，讨论了具有偶数和奇数腿数的梯子系统在基态上的定性差异。在 5.3.3 节中，针对标准二维 Heisenberg 模型计算了子晶格磁化强度以及其他一些与奈尔态相关的低能参数，并在 5.3.4 节中引入二聚化以驱动该模型进入一个量子相变，转变为普通的非磁态。分析了各种量的量子临界有限尺寸标度行为。在 5.3.5 节中，使用两种分别表现为连续相变和一级相变的 J-Q 模型，研究了更复杂的 Néel-VBS 相变情况。

图 66. 链长 $N = 1024$ 和 $4096$ 时，距离 $r = 2^{m}$ 处的自旋关联函数。误差棒小于符号大小。实线为预期形式 $A r^{- 1} \ln ( r / r_{0} )^{1 / 2}$，其中 $A = 0.21$，$r_{0} = 0.08$。虚线显示了形式 $\propto r^{- 1}$ 作为比较。这些结果分别使用逆温度 $\beta = 2^{13}$（$N = 1024$）和 $\beta = 2^{14}$（$N = 4096$）获得，这对于 $T \to 0$ 的收敛是足够的。

5.3.1 海森堡链#

$T = 0$ 时的自旋关联。在第 4.3.1 节中，我们讨论了 Heisenberg 链自旋关联函数的 Lanczos 结果，并看到了一些预期中的对 $\sim 1 / r$ 临界行为的对数修正的迹象（图 34）。然而，Lanczos 方法所能达到的系统尺寸不足以定量研究这些标度修正。正如我们上面所见，使用 SSE 方法可以对包含数千个自旋的链（通过仔细检查向 $T \to 0$ 极限的收敛性）进行无偏的基态研究。图 66 展示了 $N = 1024$ 和 $4096$ 时，在距离 $r = 2^{m}$ 处的自旋关联，以对数-对数坐标绘制。为节省时间，仅计算了这些距离处的关联 [执行空间平均测量所需的时间按 $N \log ( N )$ 标度]。两种系统尺寸的结果在 $r$ 达到 $2^{7}$ 时非常吻合，表明对于 $N = 1024$，结果收敛到无限尺寸值直到该距离（因此对于 $N = 4096$，由于收敛行为应大致与 N 成线性比例，结果收敛直至 $r \approx 2^{9}$）。

预期的形式 $\vert C ( r ) \vert = A \ln^{1 / 2} ( r / r_{0} ) r^{- 1} [ 57 , 58 , 59 ]$ 在 $N = 4096$ 时直到 $r = 2^{9}$ 都得到了非常好的复现。通过拟合得到的参数 A 和 $r_{0}$ 列在图注中。如果将对数项的指数 $\sigma = 1 / 2$ 作为自由参数，基于数据通过拟合确定，该指数确实接近 0.5，但误差条相当大，大约为 0.1。为了真正仔细研究指数，应进一步增加链长（这在原则上是可行的）。

如果事先不知道对数修正的存在，并试图仅基于数值计算来提取自旋关联的形式，根据图 66 中的数据，最初可能会得出衰变 $\propto 1 / r^{\alpha}$，其中 $\alpha \approx 0.85$ 的结论。然而，存在虽小但显著的偏离纯幂律的情况，只有当相对统计误差足够小时才能检测到。在图中显示的数据中，误差条太小以致无法分辨。虽然绝对误差很小，对于 $N = 4096$ 通常约为 $\approx 2 \times 10^{- 6}$（基于大约 $10^{6}$ 次蒙特卡洛扫描，在 $\dot{\beta} = 2^{14}$ 并使用改进估计量），但相对误差要大得多，对于 $r = 2^{7}$ 约为 0.002，对于最长距离 $r = 2^{11}$ 则增长到 0.02。

图 67. 圆圈显示了海森堡链在温度 T 为 $2^{- m} , m = 0 , 1 , 2 , \ldots$ 形式下的均匀磁化率的 SSE 结果，对于不同 2 的幂次链长直至 $\bar{N} = 4096$。误差条远小于圆圈。精确的 $N = 16$ 结果以实线曲线显示，虚线曲线显示了低温形式 (284)。渐近 $T \to 0$ 值 $J \chi = 1 / \pi^{2}$ 由图形左侧的短水平线段指示。

低温磁化率。对数修正同样重要，例如在 $T > 0$ 的海森堡链中。最突出的例子是均匀磁化率，对其低能场论的重整化群研究预测了如下形式 [7, 213]

其中 $T \to 0$ 的值 $\chi ( 0 ) = 1 / 2 \pi c$，且 $c = J \pi / 2$ 为自旋子速度，与海森堡链基态的贝特拟解结果一致。贝特拟解也可推广到 $T > 0$ 的情况。在低温 $( T / J < 0.02 )$ 下，渐近表达式（284）与参数 $T_{0} = 7.7$（通过拟合贝特拟解结果确定）吻合良好 [7]。需要注意的是，（284）中的对数修正意味着收敛到 $T = 0$ 极限非常缓慢，因此若不了解这类修正的存在，将难以外推数值结果。（284）的高阶对数修正亦已知晓 [213]，但此处我们仅考虑主阶修正。

图67展示了使用 $N = 2^{n}$ 形式的链长和温度 $T / J = 2^{- \bar{m}}$ 通过SSE方法得到的磁化率结果。由于有限长链中存在激发隙（在4.3.1节讨论），在温度 $T \sim 1 / N$ 以下磁化率指数衰减至零。图中使用的对数-线性尺度使这些有限尺寸效应非常清晰，同时也展示了当有限尺寸能隙变得小于 $T$ 时，结果如何迅速收敛到热力学极限形式。使用文献[7]中确定的 $T_{0} = 7.7$ 值，（284）形式在 $T / J$ 低于0.05时吻合良好（实际上，出于某些未知原因，吻合程度似乎甚至优于文献[7]）。

除证实低能场论的正确性外，（284）形式的主要用途在于它不仅适用于此处考虑的简单海森堡链，还适用于任何与海森堡链处于相同相的1D自旋系统（在足够低温度下）。这包括4.3.2节讨论的阻挫 $\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$ 链，其中 $J_{2} / J_{1}$ 小于二聚化点 $( J_{2} / J_{1} )_{c} \approx 0.241$。恰好在二聚化点，主导对数修正应消失 [164]（但仍存在其他高阶修正），而在相变点之上，自旋能隙打开，即使 $N = \infty$ 时 $\chi \to 0$。渐近形式（284）应在包含长程相互作用的链的QLRO( )相 [4.3.3节讨论的模型] 以及许多其他系统中有效。参数 $c$ 和 $T_{0}$ 依赖于模型参数，拟合数值数据提供了一种提取速度 $c$ 的方法。如我们将在下一节所见，由奇数条链组成的梯子系统也与单链处于相同的临界相，故（284）同样适用。

5.3.2 梯子系统#

梯形晶格由固定数量的 $L_{y}$ 条耦合1D链（常称为梯子的腿）组成，链长 $L_{x}$ 趋于无穷（或实际中，$L_{x} / L_{y}$ 足够大使得结果收敛到此极限）。这种几何结构中的强关联量子系统 [9, 214, 215] 作为介于1维与2维之间的插值手段扮演着重要角色。许多材料呈现出 $L_{y}$ 较小的弱耦合梯子结构，这促使了对 $L_{y} = 2, 3$ 等情况的研究。在更基础的层面上，观察1D系统的特殊性质如何随 $L_{y}$ 增加而演化，并趋近2D极限 [216]，是极具意义的。此处我们将研究 $L_{\nu}$ 最高为6的海森堡梯子的一些最基本性质。

$S = 1/2$海森堡梯子物理最关键的方面在于，对于偶数 $L_{y}$（梯子每级梯级上的格点数）和奇数 $L_{y}$，其低能性质完全不同。当 $L_{y}$ 为偶数时，总存在一个自旋能隙（该能隙在 $L_{y} \to \infty$ 时消失），而 $L_{y}$ 为奇数的系统则是无能隙的，其性质类似于单条一维链（在某个随 $L_{y} \to \infty$ 而消失的能量或温度标度以下）。这种差异可以根据价键的简单图像 [217] 粗略理解，如图 68 所示。如果梯级之间没有耦合（$J_{x} = 0$, $J_{y} > 0$），双腿梯子的基态只是梯级单重态的乘积——这是一个唯一的（非简并）态，与第一激发态（其中一个单重态被激发为三重态）之间存在能隙 $J_{y}$。另一方面，对于三腿梯子，单个梯级的基态是二重简并的，$S_{z} = \pm 1/2$。当梯级之间耦合时（$J_{x} > 0$），这种简并被解除，但对于任意比值 $J_{x}/J_{y}$，简并的残余仍然存在，这可以理解为来自 $L_{x}$ 个相互作用的 $S = 1/2$ 自由度。这导致其低能性质类似于单条 $S = 1/2$ 链。在双腿梯子中，对于任意 $J_{x}/J_{y}$，单重态-三重态激发能隙保持非零，该能隙之上的低能激发形成了一个传播的梯级三重态能带 [218]。自旋关联随距离呈指数衰减，关联长度在 $L_{y} \to \infty$ 时发散。

图 68. 在梯级间耦合 $J_{x} = 0$ 时，双腿和三腿海森堡梯子基态的示意图。竖线表示价键（单重态），圆圈表示未配对的自由自旋。双腿梯子的基态是唯一的单重态乘积态，而在三腿梯子中，每个梯级态是二重简并的，$S_{\mathrm{rung}}^{z} = \pm 1/2$。这些态中的每一个都是未配对自旋位于上链和下链两种态的对称组合。这里为每个梯级展示了未配对自旋位置和 $S^{z}$ 的一种随机构型。

在价键基矢（在第 2.2 节讨论过）中，$J_{x} > 0$ 的双腿梯子仍然以短键为主（可以说它是一个 RVB 自旋液体，但由于梯子几何的限制，称其为 VBS 也是合适的），而三腿梯子的临界态则需要键概率随键长呈幂律衰减 [217]（这种态可以称为临界 RVB 态）。对于更宽梯子的基态，同样的图像仍然有效，但与梯子行为相关的能量标度在 $L_{y} \to \infty$ 时消失，并在更高能量处发生向二维行为的交叉。相关的能量标度对于偶数 $L_{y}$ 是自旋能隙，对于奇数 $L_{y}$ 是自旋刚度——正如我们在第 3.5 节讨论 KT 相变时所述，幂律关联足以维持非零的自旋刚度，在单链和奇数腿梯子中出现的 $\sim 1/r$ 关联正是这种情况。

此处我们将只讨论空间各向同性耦合的结果：$J_{x} = J_{y}$。对于偶数 $L_{y}$，在 $y$ 方向采用周期边界条件和开放边界条件的行为类似（对于 $L_{y} = 2$，周期边界条件仅相当于加倍 $J_{y}$），而对于奇数 $L_{y}$，系统对所有 $J_{x}/J_{y}$ 都表现为无能隙性质仅适用于开放的 $y$ 边界。在"管状"几何（$y$ 方向周期边界条件）中，奇数 $L_{y}$ 会导致阻挫，此时行为更复杂，可能出现多个基态相 [219]。在下文讨论的 SSE 计算中，所有情况都使用了开放的 $y$ 边界，而 $x$ 边界是周期的。

自旋关联。 在前一节中，我们研究了单链的自旋关联，并确认了对 $\propto 1/r$ 临界行为存在对数修正。现在我们检验奇 $L_{y} > 1$ 的情况，同时也考察偶 $L_{y}$ 情况下指数衰减的关联。

在y方向具有开放边界条件的情况下，关联函数 $\langle \mathbf{S}_{i} \cdot \mathbf{S}_{j} \rangle$ 不仅依赖于两个自旋之间的间距 $\mathbf{r}_{i j}$，还取决于两个自旋的y坐标。这里我们使用了两个自旋之间最大反射对称（在y方向上）的距离，即对于给定的位于 $( x , y )$ 的自旋，其中 $x = 1 , \ldots , L_{x}$，$y = 1 , \ldots , L_{y}$，关联函数是通过与位于 $( x + L_{x} / 2 , L_{y} - y + 1 )$ 的自旋计算得到的，同时在x方向考虑了周期边界条件。然后对所有 $( x , y )$ 进行平均。这只是一种便捷的选择，具体如何处理y方向上的间隔并不重要，因为当x间隔很大时，关联函数对y方向的依赖很弱。

图69. 在不同宽度 $L_{y}$ 的 $L_{x} \times L_{y}$ 海森堡梯子中，最大x距离处（如文中所述对y方向间距进行了某种平均）的自旋关联函数随长度 $L_{x}$ 的变化。虚线形式为 $C ( x ) \propto 1 / x$

图69展示了宽度 $L_{\mathrm{y}} = 1 - 6$ 的梯子中长距离关联函数随 $x = L_{x} / 2$ 的变化。偶数 $L_{y}$ 和奇数 $L_{y}$ 之间的定性差异很明显，且对于奇数 $L_{x} > 1$ 的情况，当 $L_{x}$ 很大时，其行为与 $L_{x} = 1$ 系统非常相似，与形式 $C ( x ) \propto 1 / x$ 的偏差可以归因于对数修正（此处不作更详细分析）。对于偶数 $L_{x}$，当x很大时，关联函数呈指数衰减，但对于 $L_{\mathrm{y}} = 6$，已经可以看出一种不同的短距离行为开始出现。对于非常大的 $L_{y}$ 和x达到 $\sim L_{\mathrm{y}}$ 时，对于偶数 $L_{y}$ 和奇数 $L_{y}$，关联函数都应呈现二维形式，即 $C ( x )$ 应趋近于二维平方交错磁化强度的值 $m_{s}^{2} \approx 0.095$（图5中的数据有说明）。对于 $x \gg L_{y}$，关联函数应发生交叉，要么变为指数衰减（对于偶数 $L_{y})$，要么变为带有对数修正的 $1 / x$ 形式（对于奇数 $L_{y})$ 。在图69中可以看出，在 $x = 8$ 处，对于 $L_{\mathrm{y}} = 4 - 6$，关联函数几乎相同，但对于更大的x，关联函数仍随 $L_{y}$ 增加而增加。要清晰地观察到二维行为以及随后的交叉（进入偶数或奇数腿梯子的渐近行为），需要大得多的 $L_{y}$。

磁化率。 磁化率对温度的依赖性提供了一种方便的方法来提取偶数 $L_{y}$ 梯子中的自旋能隙 $\Delta$，同时研究预期渐近形式 (284) 对于奇数 $L_{y}$ 的适用性也很有趣。图70显示了 $L_{\mathrm{v}} = 1 - 6$ 的SSE结果，每个温度下的 $L_{x}$ 都足够大，以精确收敛到热力学极限。对于最低温度，使用了 $L_{x} = 1024$。对于固定的 $T_{\ast}$，在 $T / J > 0.2$ 时可以观察到随 $L_{y}$ 变化的收敛性，其中 $L_{\mathrm{v}} = 4 , 5$ 和 6 的结果几乎相同。收敛曲线对应于二维极限（在大的 $L \times L$ 晶格上的结果见图73）。偶数宽度和奇数宽度梯子之间的定性差异在低温下很明显，偶数 $L_{x}$ 的 $\chi$ 在某个依赖于 $L_{y}$ 的温度以下呈指数衰减至零。对于奇数 $L_{y}$，可以注意到 $T \to 0$ 的值随 $L_{y}$ 增加而减小。可以论证，对于小的 $L_{y}$，每个梯级的趋零温度磁化率大致应与 $L_{y}$ 无关[216]；因此，每个自旋的磁化率应大致按 $1 / L_{y}$ 标度。此处我们不对渐近形式 (284) 进行任何拟合，因为它仅对 $L_{y} \ge 3$ 的极低温度有效（并且除非使用极低的温度，否则可能需要高阶对数修正[216]）。注意，$\chi ( 0 ) \sim 1 / L_{y}$ 的形式与二维行为 $\chi_{2 D} ( 0 ) > 0$ 并不矛盾，因为梯子行为（对于偶数 $L_{y}$ 和奇数 $L_{y}$ 均成立）仅适用于低于某个 $T^{*} ( L_{y} )$ 的温度，该温度在 $L_{\mathrm{y}} \to \infty$ 时趋近于零。

图 70. 不同宽度 $L_{y}$ 的海森堡梯子磁化率的温度依赖性。长度 $L_{x}$ 对于每个温度都足够大，以精确表示 $L_{x} \to \infty$ 的极限。右图以更详细的尺度显示了 $L_{\mathrm{y}} = 2 - 6$ 的低温行为。由于使用了改进的磁化率估计量 (283)，统计误差不可分辨。

现在让我们分析偶数宽度梯子的自旋能隙 $\Delta$ 对 $L_{y}$ 的依赖性。由于基态是单重态（非磁性，且 $\langle \dot{M}_{z}^{2} \rangle = \bar{0}$），而第一激发态是三重态（磁性，且 $\langle M_{z}^{2} \rangle = 2 / 3$），因此人们大致预期低温磁化率呈指数行为，$\chi \sim \exp ( - \Delta / T )$。然而，由于在 ∆ 之上存在一整个连续的磁性态（对于 $L_{x} \to \infty$），这种行为会被修正。通过微扰分析弱耦合极限 [218]，得到了 2-腿梯子均匀磁化率的低温形式，

对于任何 $J_{y} / J_{x}$ 和偶数 $L_{y} ( \mathrm{for} L_{x} \infty )$，可以预期这种形式在 $T \to 0$ 时渐近成立。通过考虑 (285) 中 $\sqrt{T} \chi$ 的对数，可以从数值 $\chi ( T )$ 数据中提取出能隙，得到

因此，对 $- T \ln ( \sqrt{T} \chi )$ 与 T 拟合所得直线的截距即为能隙 ∆。图71展示了以此方式分析的数据。在所有情况下，低温下均呈现线性行为，直线拟合给出了图中列出的能隙。场论处理表明，自旋S梯子系统的能隙随 $L_{y}$ 的减小方式为 $\Delta \sim \mathrm{exp} ( - \pi S L_{y} )$ [220]，这与此处得到的结果（以及早期精度较低的计算[216]）吻合不佳。需要针对更大的 $L_{y}$ 进行计算（这完全可以实现），以更详细地研究此问题。

图71 根据方程(286)分析的偶数宽度梯子的磁化率。在低温线性区域的外推，以穿过拟合所用数据点的直线表示，给出了图中所示的自旋能隙。统计精度大致如所示的数字位数。

对奇数腿梯子进行比目前更大 $L_{y}$ 和更低温度的研究也很有意义，以提取(284)中参数 c 和 $T_{0}$ 的 $L_{y}$ 依赖性。另一种方法是将梯子的低能性质映射到具有长程相互作用[221]的单条链上，这种相互作用源于孤立的梯级一旦耦合，其 $S = 1 / 2$ 自由度并非完全局域在单个梯级上。局域化长度随 $L_{y}$ 增大，最终当 $L_{\mathrm{y}} \to \infty$ 时，应该会得到具有足够长程相互作用以产生奈耳序（如第4.3.3节所述）的有效模型。然而，这种映射方法尚未在大 $L_{y}$ 的梯子上得到检验，因此还应在更低温度下和针对更大 $L_{y}$ 研究完整的梯子系统。

5.3.3 二维中的长程序#

我们在第2.1节和第4.4.2节中已经讨论了二维海森堡模型反铁磁有序（奈耳）基态的性质。对子晶格磁化强度平方的最佳当前可用QMC结果[50]的外推，如图5所示，对于无限系统得到 $m_{s} = 0.30743 ( 1 )$，其中(1)表示前一位数字的统计误差。此结果与线性自旋波结果 $m_{s} = 0.3034$ 仅偏离约1%。高阶自旋波计算[40, 41, 42]给出 $m_{s} \approx 0.3070$ 。如第2.1节所述，这种与实际值的良好吻合可归因于量子涨落并不那么强，仅将子晶格磁化强度从经典值降低了约40%。其他重要的基态量，例如自旋刚度[222]、横向磁化率和自旋波速度，与自旋波理论同样吻合良好（但此处我们不进行具体数值的详细比较）。

在讨论一些通用的数据拟合问题之后，我们将在此处通过对 $T \to 0$ 收敛的SSE结果（通过对所有感兴趣的量使用如图65所示的检验来确认其正确收敛）进行有限尺寸外推，来计算能量和自旋刚度。除了通常使用的 $L \times L$ 晶格，我们还将考虑 $L_{x} = 2 L_{y}$ 的矩形 $L_{x} \times L_{y}$ 晶格。这可以对这种外推进行一致性检验。我们还讨论了有限温度磁化率，并使用足够大的 $L \times L$ 晶格来完全消除所考虑温度范围内的有限尺寸效应，将其外推到零温。最后，我们讨论了关联长度在 $T \to 0$ 时的发散，这是奈耳有序基态的另一个重要体现[5]。

之前已有大量论文展示了二维海森堡模型基态参数的外推结果，例如参考文献[223, 224, 205, 225, 226]。对于大多数实际应用而言，已有的精度已经足够。然而，由于该模型是量子磁性中最重要的原型系统之一，继续建立更精确的基准计算（例如，可用于测试其他方法[94]）仍然是有益的。下文展示的部分结果，特别是基态能量，代表了迄今为止最精确的计算。除了各种物理量的外推 $N = \infty , T = 0$ 值之外，它们的有限尺寸和温度修正也值得关注，因为基于场论方法[227]已经对这些修正进行了非常详细的预测，而数值检验有助于确定低能理论的有效范围。这些修正通常统计噪声更大，需要非常长的模拟才能精确确定。这里的讨论仅旨在展示其可能性。

数据拟合问题。在将有限尺寸结果外推到热力学极限时，了解基于解析计算所预期的尺寸修正形式是有帮助的。远离临界点时，维度为 $d > 1$ 的系统的尺寸修正通常采用逆系统长度 $1 / L$ 的多项式形式，但有时主导修正是 $\propto 1 / L^{a}$，其中 $a$ 是大于 1 的整数。如第 2.4.1 节所述，二维海森堡模型子晶格磁化强度的主导尺寸修正可以从奈尔态的线性自旋波理论中得到，并且是 $\propto 1 / L$。这个结果也是序参量为矢量这一事实的更一般性推论，如第 3.3 节所述。每个格点的基态能量 $E / N$ 的主导修正为 $\propto 1 / L^{3}$，由于这个幂次较高，即使基于较小的晶格[205]也相对容易获得高精度的能量。自旋刚度的主导修正为 $\propto 1 / L$ [228]。即使不知道主导幂次，只要数值精度足够高，通常也可以基于数据凭经验找到它。将数据拟合到多项式时，可能会发现线性项（以及可能的高阶项）的系数非常小，这很可能意味着它实际上应该严格为零。最好使用允许精确指定应包含 $1 / L$ 的哪些幂次的拟合程序。

对于一组给定的数据（这里是对一组晶格尺寸计算的一个物理量，但我们也将使用相同技术来拟合作为 $T$ 的幂次的温度相关量），应包含足够高的 $1 / L$ 幂次，以使拟合在统计上可靠，即对于大数据集，每个自由度的 $\chi^{2}$（$\chi^{2} / \mathrm{do} \dot{\mathrm{f}}$）应接近 1（如数据分析标准教材中所讨论的）。拟合参数的统计误差可以根据数据的误差棒计算。最安全的方法（不依赖于误差很小的任何假设，尽管通常误差应该很小）是进行大量拟合，在数据中添加高斯噪声（噪声的标准差等于相应的误差棒），然后计算所得参数值分布的标准差。注意，当增加拟合参数数量（这里是 $1 / L$ 的幂次）时，最佳拟合参数的误差棒通常会增大（且 $\chi / \mathrm{dof}$ 也可能变差）。因此，不应包含比将 $\chi^{2} / \mathrm{dof}$ 降低至接近 1 所需的更多的高阶项。然而，为了防止高阶修正可能残留的影响，可能需要在严格达到可接受 $\chi^{2} / \mathrm{dof}$ 所需参数数量的基础上，再增加一个参数。

图72. 纵横比$R = 1$（其中$N = L^{2}$）和2（其中$N = 2 L^{2}$）的二维海森堡晶格基态能量（左）和自旋刚度（右）的尺寸依赖性。对于$R = 2$，有限晶格上的自旋刚度在$x$和$y$方向上是不同的。图中显示了这两个值以及它们的平均值。曲线是文中讨论的多项式拟合。

基态能量与自旋刚度。 每格点能量$E / N$和自旋刚度$\rho_{s}$分别使用估计量(250)和(280)进行评估，其中能量已针对SSE模拟中每个键算符减去的常数$1 / 4$进行了调整。如5.2.4节所述，外推到热力学极限的自旋刚度值应乘以因子$3 / 2$，以说明SSE模拟中自旋旋转对称性未被破坏这一事实。

比较不同形状晶格的外推结果是检查数据拟合中包含的有限尺寸修正之外是否存在可探测的有限尺寸校正影响的好方法 [82, 94]。图72展示了$N = L \times L$（纵横比$R = 1$）和$N = 2 L \times L$（$R = 2$）晶格的结果，以$1 / \sqrt{N} \mathrm{~ ( = ~} 1 / L$ for $R = 1$ and $\propto 1 / L$ for $R = 2 )$为自变量作图。这些结果是使用高达$32 \times L$的逆温度$\beta$计算得到的，针对的是最大的系统（对于$R = 1$，$L$高达40；对于$R = 2$，$L$高达32）。最大的系统使用了数百CPU小时。误差棒太小，在图中不可见。以$32 \times 32$系统为例，基于$5 \times 10^{6}$次蒙特卡洛扫描得到的结果为$E / N = - 0.6695115 ( 8 )$和$\rho_{s} = 0.12606 ( 5 )$。其他系统的误差棒也具有相似的数量级。

对于$R = 1$，使用$L \geq 6$的系统和一个5阶多项式（不含线性项和二次项，如前所述，这些项预计不存在）给出$E / N = - 0.6694421 ( 5 )$，而对$L_{y} \ge 6$的$R = 2$数据进行的同阶拟合得到$E / N = - 0.6694422 ( 6 )$。这些结果在统计上完全一致，因此可以将它们的统计加权平均值$E / N = - 0.6694421 ( 4 )$作为基态能量的最终估计值。应注意，数值结果略微依赖于拟合中包含的晶格以及所使用多项式的阶数。一旦拟合在统计上可靠，这些波动应与统计误差一致。

现在讨论刚度，首先注意到，对于纵横比$R \neq 1$的有限晶格，存在两个刚度常数，分别对应在$x$和$y$方向施加的相位扭转。这两个值以及它们的算术平均值和$R = 1$的值均显示在图72的右面板中。可以看出，对于外推，$R = 1$的刚度和$R = 2$刚度的平均值比$R = 2$单独的$x$和$y$值表现更好，因此后者未包含任何拟合。对$R = 1$和$R = 2$（使用$L \geq 6$和$L_{x} \ge 8$数据）进行二次拟合分别给出$\rho_{s} = 0.12065 ( 4 )$和$\rho_{s} ( \infty ) = 0.12070 ( 6 )$。取这些统计一致值的平均值，并包含旋转因子$3 / 2$，得到无限系统刚度的最终估计值$\rho_{s} = 0.18100 ( 5 )$。

$T > 0$时的性质。由于默明-瓦格纳定理排除了二维海森堡模型在任何$T > 0$下出现磁有序的可能性，因此为与有序基态的行为一致，关联长度必然在$T \to 0$时发散。这一行为与经典二维海森堡模型相似，后者的关联长度呈指数发散。事实上，量子系统的长程关联和涨落可以通过重正化耦合[5]形式化地映射到经典系统，这种映射成立的低温区域被称为重正化经典区域（renormalized classical regime）。该区域的许多预言源自场论处理方法[5, 227, 84]。例如，关联长度和磁化率等物理量对温度的依赖关系已推导出极其详尽的结果。这些被视为在$T \to 0$极限下精确的预言形式，仅依赖于基态参数（如自旋刚度$\rho_{s}$和自旋波速度$c$）。

首先考察均匀磁化率，它在$T \to 0$时（对无限系统）应有如下形式[227]：

这里唯一依赖于模型的参数是基态自旋刚度，我们前面已经确定了它。为了将数值数据与这一结果进行比较，我们首先必须确保能够达到热力学极限。对于固定温度，我们可以研究系统尺寸对行为的影响，这里我们使用具有 $N = L^{2}$ 的方形晶格。图 73 的左图显示了几个晶格尺寸的结果。对于任何有限系统，由于单重态基态和谱中的有限尺寸能隙，磁化率在 $T \to 0$ 时消失。正如我们在第 4.4 节所讨论的，最低的激发是量子转子态，其能隙按 $1 / N$ 标度。对于有限系统，$T / J < 1 / L$ 的磁化率由这些激发主导。一旦 $T > 1 / L$，自旋波也可以被激发，因此可以推测，对于收敛到热力学极限，至少需要 $L > 1 / T$ 的系统尺寸。这里我们不详细研究收敛性，仅基于图 73 中的数据得出结论，$L = 256$ 的结果在 $T / J = 0.03$ 以下已安全收敛，可用于检验公式 (287)。

图 73. 使用多个 $L \times L$ 晶格计算得到的二维海森堡模型的均匀磁化率。右图以更详细的尺度显示了 $L = 128$ 和 256 的低 T 结果（两个数据集在统计误差内重合，表明已收敛至热力学极限）。右图中的实线是对 $L = 256$ 的数据（仅使用 $T / J \le 0.18$ 的数据点）按照形式 $\chi = a + b T + c T^{2}$ 的拟合，参数在文中讨论。误差棒太小而无法分辨，但对于低温下的 $L = 256$ 结果，典型值约为 $\approx 2 \times 10^{- 5}$。

解析预测的检验可以通过不同的方式进行。我们可以在这里使用从上述 $T = 0$ 数据的有限尺寸外推中提取的 $\rho_{s}$ 值，并检查 SSE 关于 $\chi ( T )$ 的结果与公式 (287) 的一致性，仅将 $\chi ( 0 )$ 作为拟合参数进行调整以获得最佳吻合。另一种方法是同时调整 $\rho_{s}$，其最佳拟合值可以与基态外推的结果进行比较。这里我们将采用一种不同的方式，该方式同时检验了 (287) 中的 $T$ 和 $T^{2}$ 修正项。通过将 SSE 数据拟合为形式 $\chi = a_{0} + a_{1} T + a_{2} T^{2}$，我们可以利用 (287) 中系数的形式，从线性项和二次项系数 $a_{1}$ 和 $a_{2}$ 中提取出相应的自旋刚度估计值 $\rho_{s} ( 1 )$ 和 $\rho_{s} ( 2 )$。如果拟合中使用的温度足够低，这些估计值应该彼此一致，并且与基态计算的结果一致。在较高温度下，包含 T 的更高次幂的修正项将变得重要，并应导致对二次形式的拟合不佳以及不同刚度估计值之间的不一致。

如果仅包含 $T / J \le 0.18$ 的数据，二次拟合效果良好，并能得到一致的指数值。该拟合如图 73 右图所示，从参数中提取的刚度常数分别为：由斜率得到的 $\rho_{s} ( 1 ) = 0.180 ( 3 )$ 和由二次修正项得到的 $\rho_{s} ( 2 ) = 0.179 ( 6 )$。这些值与基态结果 $\rho_{s} = 0.18100 ( 5 )$ 吻合良好（但噪声稍大）。固定形式 (287) 并提取 $\chi ( 0 )$ 和 $\rho_{s}$ 的最佳拟合值可以得到噪声更小的估计 [225]，但由于我们仍然在处理 $T = 0$ 量的修正项，要达到与我们在 $T = 0$ 时外推 $\rho_{s}$ 的有限尺寸值相当的精度并不容易。

图 74. 使用自旋结构因子定义（70）计算得到的二维海森堡模型在不同晶格尺寸下温度依赖的关联长度（绘制为温度倒数的函数），除以来自公式（72）的因子 $\sqrt{15 / 16}$。虚线是期望的渐近 $( L = \infty , T \to 0)$ 形式（289），参数为先前确定的 $\rho_{s} = 0.1810$ 和 $c = 1.660$。

根据上述拟合提取的 $T \to 0$ 磁化率为 $a_{0} = \chi ( 0 ) = 0.04378 ( 3 )$。正如我们在第4.4.2节中讨论的，在 $T = 0$ 对称性破缺态下的横向磁化率 $\chi_{\perp}$ 为该数值乘以 $3 / 2$，得到 $\chi_{\perp} = 0.06567 ( 4 )$。结合自旋刚度（spin stiffness），该数值使我们能够获得另一个基态参数：自旋波速度。我们在第3.5节中将自旋刚度常数作为一种弹性模量进行了讨论。这种将量子磁体长波特性视为弹性介质的方法也被称为流体动力学描述（hydrodynamic description）[229]。其中，横向磁化率类似于质量密度，我们可以根据下式将自旋波速度与自旋刚度和磁化率联系起来：

使用上述提取的 $\rho_{s}$ 和 $\chi_{\perp}$ 值可得速度 $c = 1.6602 ( 6 )$。回顾线性自旋波理论（第2.1.1节）给出的 $c = {\sqrt{2}}$。包含 $1 / S$ 修正的自旋波理论很好地捕捉了这两个数值之间的差异，并且 $c$ 的“重整化因子”也可以通过围绕伊辛模型（Ising model）的级数展开非常精确地计算出来[230]。

无限系统关联长度的温度依赖性也是已知的，包括对主导指数发散[5, 227]的修正：

我们已经确定了这里涉及的两个参数，即基态自旋刚度（spin stiffness）和速度，并将使用上述引用的数值，基于关联长度的“二阶矩”结构因子定义(70)，直接将解析预测与SSE结果进行对比。

由于关联长度是一个迅速发散的物理量，需要在非常大的晶格上计算才能收敛到低温下的热力学极限。晶格长度$L$应比$\xi$大数倍才能实现完全收敛。图74展示了$L = 256$以内若干晶格尺寸的结果。尽管过去已有对更大系统的研究[225, 226]，但这里$T / J \ge 0.25$时已很好收敛的结果表明，形式(289)非常准确地描述了数值数据。文献[225, 226]中讨论了小的偏差，这些偏差可理解为形式(289)的高阶修正所致。

5.3.4 二聚化系统中的量子相变#

我们现在研究图4(b)所示的二维柱状二聚化海森堡模型中的量子相变。基态的相变随耦合比$g = J_{2} / J_{1} > 1$发生，由量子涨落驱动，这对应于二聚体上单重态密度的增加，并在某个点导致$g \approx 1$时存在的奈耳序丧失。我们已经在图5中考察了该体系子晶格磁化强度的一些SSE结果，其行为表明在耦合比$g \approx 1.9$处发生了奈耳有序基态与非磁性基态之间的相变。在本节中，我们将应用临界点处的有限尺寸标度方法，我们在第3.3节中讨论过经典系统中的这一方法，并在第3.6节中推广到了量子系统。这里，我们首先分析在足够低的温度下计算的若干物理量的临界行为，以获取基态临界行为，然后讨论$T = 0$量子临界点（quantum critical point）在非零温度下的后果。在此之前，我们先花几段文字讨论一下非均匀耦合模型下SSE算法的实现。

SSE方法。SSE算法在海森堡模型中的一个吸引人之处在于，耦合强度仅在对角更新中进入接受概率(265)和(266)，其中耦合常数$J$被吸收进逆温度$\beta = J / T$中。对于非恒定的耦合，我们只需为每个键$b$定义$\beta_{b} = J_{b} / T$。在伪代码27的实现中，需要使用某种方案来识别插入尝试中生成的键$b$所对应的耦合。在具有大量不同耦合（甚至可以是随机的，每个键都不同）的系统中，耦合（以$\beta_{b}$的形式）应作为列表存储。对于仅有两种不同耦合的系统（例如这里考虑的二聚化晶格），最好对键进行排序，以便快速确定耦合，例如，$b \leq N / 2$和$b > N / 2$分别对应$J_{2}$和$J_{1}$。注意，对于算符移除，与代码27中均匀耦合的情况不同，我们现在也需要提取当前键$b = S ( \boldsymbol{p} ) / 2$，以确定概率。

$T \to 0$ 的有限尺寸标度。在分析基态时，我们可以像上一节标准均匀2D海森堡模型那样，在足够低的温度下模拟系统，以使得所有感兴趣的物理量都收敛到$T = 0$的结果。这种方法在多种二聚化系统中均有讨论，参见例如文献[85, 88]（以及许多更早的工作）。另一种方法是在逆温度$\beta = L^{z}$下研究系统，其中$z$是动力学临界指数（我们在3.6节讨论过）[231]。这种方法的动机如下，通过推广有限尺寸标度假设（64）：在量子系统中，标度函数$f ( \xi / L )$应被一个含有两个参数的函数$f ( \xi / L , \xi_{\tau} / L_{\tau} )$所取代，其中虚时维度上的关联长度与空间关联长度$\xi$的关系为$\xi_{\tau} \sim \xi^{z}$（这定义了动力学指数），而系统在虚时方向上的长度为$L_{\tau} = c / T \sim \beta$（其中$c$是速度）。如果我们选择$\beta \propto L^{z}$，那么标度函数可以写为$f [ \xi / L , ( \xi / L )^{z} ]$，这只是一个单一参数$\bar{\xi} / L$的函数。因此，有限尺寸标度方法可以完全像3.3.2节中讨论的经典系统那样使用。如果我们对每个L取$\beta \to \infty$极限（实际计算中取足够大的有限$\beta$以收敛到该极限），情况也是如此，因为此时$\xi_{\tau} / L_{\tau} \to 0$，标度函数中又只剩一个参数$\xi / L$。

图75. 二聚化海森堡模型的Binder累积量（左图）和自旋刚度（x方向）乘以系统长度（右图）。不同L曲线交叉点趋于耦合比$g$的临界值。误差棒远小于符号大小。

已有充分证据表明，在二聚化海森堡模型中$z = 1$，因此我们将使用$\beta = L$的系统。这使得我们能够研究比$\beta \to \infty$极限情况更大的系统，尽管在事先并不明确哪种方法最终更好，因为主导有限尺寸标度行为的修正项可能不同。这里我们使用的L最大到$L = 128$。我们还将通过研究与z相关的物理量，明确检验$\beta = L$的系统是否表现出与$z = 1$一致的行为。

我们首先通过考察在$g_{c}$处应呈现尺寸无关性的物理量来定位临界耦合。图75显示了Binder累积量和自旋刚度对g的依赖关系，其中自旋刚度乘以L以补偿预期的量子临界标度形式$\rho_{s} \sim 1 / L$，该形式来源于将经典形式（99）中的$d$替换为$d + z = 3$。

Binder累积量的定义遵循（77），其中分量数$n = 3$。但需要注意的是，（77）的定义使用了（75）中的完整标量积$m^{2} = {\dot{\mathbf{m}}} \cdot{\mathbf{m}}$，而这里使用SSE方法我们只计算了z分量期望值$\langle m_{z}^{2} \rangle$和$\langle m_{z}^{4} \rangle$（非对角分量更难以评估[232]）。通过对经典3D单位矢量的z分量cos(Θ)进行角度积分，可以方便地找到联系这些量的几何因子，得到$\langle \overline{{m}}^{2} \rangle = 3 \langle m_{z}^{2} \rangle$和$\langle m^{4} \rangle = 5 \langle m_{z}^{2} \rangle$。对于定位临界点，这些因子不起作用，我们也可以使用定义为$R_{2 z} = \langle m_{z}^{4} \rangle / \langle m_{z}^{2} \rangle^{2}$的简易Binder比值。

图76. 从Binder累积量、自旋刚度常数和关联长度的(L, 2L)交叉点提取的二聚化海森堡模型的尺寸依赖临界耦合。曲线显示了形式为$g_{c} ( L ) = g_{c} ( \infty ) + a / \bar{L^{\omega}}$的拟合结果。

由于二聚化晶格不具有 $90^{\circ}$ 旋转对称性，x 和 y 方向的刚度常数不同。虽然数值上确实有很大差异，但它们在临界点附近的标度行为却非常相似，不过 [x 方向刚度大约大 2 倍——二聚体沿 x 方向取向，如图 4(b) 所示]。图 75(b) 中仅显示了 x 方向的刚度。

在图 75 中，对于 $U_{2}$ 和 $\rho_{s} L$ 确实可以看到曲线交叉，并且在对于小尺寸 $L$ 系统（例如尺寸为 L 和 2L 的系统），交叉点经过一些显著的漂移后，似乎都在两种情况下收敛到大致相同的值。注意，$U_{2}$ 和 $\rho_{s}$ 的交叉点从相反的方向趋近于 $g_{c}$，这对于界定临界值范围很有用 [85, 88]。可以通过将适当阶数的多项式拟合到某些数据点，并重复该过程多次并添加高斯噪声来计算误差条，从而在数值上定位交叉点。图 76 显示了针对 Binder 比值、x 和 y 方向刚度常数以及 x 和 y 方向关联长度 [使用定义 (70) 计算] 进行此类过程的结果。图中还显示了形如 $g_{c} ( L ) = g_{c} ( \infty ) + a / L^{\omega}$ 的拟合结果。该形式很好地描述了所有 $L \ge 10$ 的数据（图中显示的尺寸）。所有外推得到的 $g_{c}$ 值都落在 [1.9094, 1.9096] 范围内，如果排除 $L = 10$ 的数据，该范围进一步缩小。所有量的指数都在 $2 \sim 2.5$ 范围内（$U_{2}$ 的指数最大）。将所有这些外推中得到的五个值视为独立的统计数据，得到 $g_{c} = 1.90948 ( 4 )$ 作为临界点的最终估计值。这与近期使用 $L$ 最大为 64 的晶格上相同量的 $T \to 0$ 数据得到的估计值 $g_{c} = 1.9096 ( 2 )$ 吻合良好（但误差棒更小）。交叉点偏移指数也吻合得很好。

存在一个极限位置 $g_{c}$ 使得 $L \rho_{s}$ 曲线交叉，这并不能证明 $z = 1$。还应该检查 $L \rho_{s}$ 的值是否行为良好，即，$( g , \rho_{s} L )$ 平面中是否真正形成了一个交叉点。类似于图 76 中针对尺寸为 L 和 2L 晶格在交叉点处的 $\rho_{s} L$ 所做的图证实了这一点。在更大的耦合窗口内检查 Binder 累积量时，无法找到任何负值的迹象，而负值将是一级相变（如第 3.4 节所述）中相共存的标志。因此，到目前为止的标度行为支持一个 $z = 1$ 的连续量子临界点。

现在我们可以继续进行数据坍缩拟合，以寻找关联长度指数 $\eta$，正如我们在第3.3节中对经典系统讨论的那样。这项工作已经在二聚化海森堡模型中完成，并且包含了标度修正 [85, 88]，结果与预期的三维经典海森堡普适类非常吻合。不过，我们这里仅讨论临界关联函数 (53) 中出现的指数 $\eta$，其中我们应再次将 $d$ 替换为 $d + z = 3$。交错结构因子 $S ( \pi , \pi )$ 是关联函数 (269) 的空间积分（对二维晶格求和），而用于交错磁化率 $\chi ( \pi , \pi )$ 的久保积分 (271) 则对应一个三维空间-时间积分（二维晶格求和加上虚时积分）。使用上述临界形式的关联函数，并以系统尺寸 $L$ 作为截断，执行这些积分给出 $S ( \pi , \pi ) \sim L^{d - z - \eta}$ 和 $\chi ( \pi , \pi ) \sim L^{d - \eta}$。我们将检验临界系统的这种标度行为，并使用上面找到的 $g_{c}$ 值来提取指数 $\eta$。与其在这一点上进行新的模拟（这一点仅已知在统计误差范围内），不如在现有数据集内执行多项式插值。然后，我们还可以轻松检查在偏离 $g_{c}$ 最佳估计值（正负一个误差线）的点上的行为，以检验拟合得到的 $\eta$ 对临界点位置的敏感性。分析临界行为的最简单方法是对 ln(S) 和 $\ln ( \chi )$ 相对于 $\ln ( L )$ 拟合一条直线（如同该模型 $S ( \pi , \pi )$ 的研究[88]中所做的那样）。总会有一些对标度行为的修正，如果数据误差较小，直线拟合只能用于较大的晶格。对于本文所使用的数据集，只有当使用系统尺寸 $L \ge 48$ 时，才能获得统计上可接受的线性拟合。因此，我们还将包含次主导修正，并假定以下形式

其中，预期次主导指数远小于主导指数（甚至可能为负值）。原则上可以对两个量进行联合拟合，固定两者（与常数不同）的相同值，但这里我们将独立拟合数据集。

为了更清楚地看到次主导修正的作用，我们绘制了图 77 中的 $S(\pi,\pi)$ 和 $\chi(\pi,\pi)$，并除去了主导的 $L$ 和 $L^2$ 因子。渐近行为 $L \to \infty$ 在两种情况下均为 $\sim L^{-\eta}$，其中 $\eta$ 预计很小。当前经典 3D 海森堡普适类的最佳估计是 $\eta = 0.0375(5)$ [113]。对 (290) 式的拟合给出：$S(\pi,\pi)$ 的 $\eta = 0.029(2)$，$\chi(\pi,\pi)$ 的 $\eta = 0.020(3)$。原则上，(290) 的形式当然只在 $g_c$ 严格成立，但在实际中，对于有限范围的系统尺寸，它在真实临界点附近某窗口内能很好地拟合数据。如图 77 所示，在 $g_c$ 加减一个误差棒处的两个量的插值值与中点处的值有明显偏差，但三种情况的拟合在统计上均可接受。上述统计误差主要来源于临界耦合的不确定性。(290) 中的次主导指数为：$S(\pi,\pi)$ 的 $\omega = -0.2(2)$，$\chi(\pi,\pi)$ 的 $\omega = 0.6(2)$。

图 77. 在估计的临界耦合比 $g_{c} = 1.90948(4)$ 附近交错结构因子（左）和磁化率（右）的有限尺寸标度。已除去了对应 $\tan z = 1$ 和 $\eta = 0$ 的 $L$ 的幂次。剩余渐近尺寸依赖应由 $\eta$ 的实际值决定。对 (290) 形式的拟合给出：由 $S(\pi,\pi)$ 得 $\eta = 0.029(3)$，由 $\chi(\pi,\pi)$ 得 $\eta = 0.020(4)$。拟合基于 $L \geq 8$ 的数据，但所得曲线也通过了 $L = 4$ 和 6 的点。虚线显示了在 $g = 1.90948$ 且不考虑标度修正时的行为（即保留拟合中获得的包含修正的其他参数值）。

图 77 还展示了忽略次主导修正后（保持其他参数为包含修正时的拟合值）的拟合函数。显然，修正相当显著，完全导致了 $L \approx 10$ 处两条曲线的极大值。当然，在没有次主导修正的情况下，对较大晶格可以得到更好的拟合。如上所述，鉴于所用数据的误差棒相当小，只有 $L \ge 48$ 的系统可以纳入这样的拟合。即便如此，忽略的修正仍会有一些影响。如果不包含修正，$\eta$ 的值确实会略低一些。包含修正后，可以纳入所有数据 $(L \geq 4)$，但为保险起见，此处引用的拟合及图中展示的结果仅包含 $L \geq 8$ 的数据。这些拟合仍然非常接近地通过了 $L = 4$ 和 6 的数据点，这进一步证实了所用函数形式的质量。

这里获得的数值以及文献[88]中的结果，均略低于（相差数个误差棒）上文引用的最佳经典三维海森堡值[113]。这些微小差异很可能是由进一步的标度修正导致的，但仍有必要对多种二聚化海森堡模型进行更高精度的量子蒙特卡洛计算（使用更大的晶格、更密集的晶格尺寸网格，并减小计算量的误差棒），以确凿确立与经典临界指数的吻合。鉴于 $\mathrm{图 4 (c)}$ 中所示的交错二聚体模型（以及某些其他二聚化模式）似乎显示出与预期的三维海森堡临界指数存在微小但统计上显著的偏差[89, 233, 234]，这一点尤为有趣且重要。这些系统可能并不属于最初提出的不同普适类[89]，因为偏差也可能源于异常大的标度修正[233, 234]。无论偏差的根本原因如何，其起源仍不清楚，有待进一步检验（相关研究正在进行中[234]）。

图78. 一般二聚化二维海森堡模型耦合-温度平面中的交叉行为（由粗曲线指示）。量子临界点（圆圈）控制 $T \geq 0$ 量子临界（QC）区域内的行为，直至晶格效应变得重要的温度（即 $\bar{\xi} \sim 1$ 时）。对于远离临界点的耦合，存在一个低温交叉，进入重整化经典（RC）或量子无序行为。对于二维系统，仅在 $T = 0$ 时存在奈耳序。分隔不同区域的曲线厚度反映了交叉并非尖锐，而是在显著的温度窗口内发生。

$T > 0$ 时的量子临界标度。量子临界性最显著和最重要的方面之一是，$T = 0$ 临界点的性质（普适类）也能强烈影响 $T > 0$ 的系统，通常可达相当高的温度[84]，并且即使系统不严格具有对应 $T = 0$ 临界点的耦合（但接近该点，我们将在下文更精确地说明这一含义）也是如此。因此，临界涨落和标度行为可能比经典系统[142]更普遍地成为量子系统的特征，在经典系统中，临界行为通常只能在非常接近临界点处观察到。

在量子系统中，当温度是主导能量尺度时，可以预期出现临界行为。在二聚化二维海森堡模型中，相变奈耳侧和非磁侧的能量尺度分别是自旋刚度 $\rho_{s}$ 和单重态-三重态激发能隙 ∆。因此可以预期，对于 $g < g_c$ 和 $g > g_c$，量子临界性的表现在 $T \gg \rho_{s}$ 和 $T \gg \Delta$ 时出现。在连续场论中，量子临界行为没有上限温度，而对于晶格哈密顿量，当温度超过某一阈值，自旋因热涨落而变得基本独立时（连续场论无法描述这种行为，因为其中不存在独立的离散自旋），必然会出现交叉，进入另一种行为。与高温极限的交叉可以定义，例如，依据当 $T \gg J$ 时获得的居里（独立自旋）磁化率 $\chi 1/(4T)$。正如我们将在下面看到的，量子临界磁化率是完全不同的。

综合以上讨论，对于二维二聚化海森堡模型，可以预期图78所示的一种通用的"交叉图"（cross-over diagram）[5, 84]。在高温交叉线以下的温度范围内，存在三个不同的"区域"（regimes），其中体系的行为由问题中对应的最大能量尺度 $\rho_{s} , \Delta$ 或 $T$ 所支配。在重整化经典区域（RC regime）中，自旋刚度（spin stiffness）占主导，关联长度按照(289)式随 $T \to 0$ 发散，正如我们在第5.3.3节中针对标准均匀海森堡模型明确看到的那样。在量子无序区域（QD regime）中，当 $T \to 0$ 时，关联长度趋于一个有限值，在 $T \approx \Delta$ 以下 $\xi$ 大致保持常数。在量子临界区域（QC regime）中，当 $T \to 0$ 且位于 $g_{c}$ 时，关联长度应按以下关系发散：

这种行为可以通过将 $d$ 维量子系统的路径积分映射到经典 $( d + z )$ 维系统来理解。在此映射中，量子系统的耦合参数 $g$ 对应 $( d + z )$ 维系统的温度 $T$。我们可以把量子系统的温度标度问题看作是虚时维度上的有限尺寸标度问题，其中有效系统的厚度是有限的： $L_{\tau} \propto 1 / T$ （即，我们保持 $T$ 不变，只改变有效系统的厚度）。根据动力学指数的定义，在 $g_{c}$ 处我们有

如果时间维度的厚度 $L_{\tau} = \infty$ ，那么在该点 $\xi$ 和 $\xi_{\tau}$ 都会发散，但在有限 $L_{\tau}$ （但空间大小 L 为无穷大）的情况下，我们应在标度公式中用 $L_{\tau} \propto 1 / T$ 替换 $\xi_{\tau}$ （这完全类似于经典有限尺寸标度中当 $L < \xi$ 时的做法），由此(292)式即给出(291)式。

基于这些思想，通过重整化群方法和 $1/N$ 展开（其中 $N$ 是 SU(N) 对称自旋的分量数，对于物理上的 $S = 1 / 2$ 自旋， $N = 2$）在连续场论 [ (2+1) 维非线性 σ 模型] 中进行了定量且极其详细的研究 [84]，已推导出许多物理量的量子临界标度形式。要点在于，即使当 $g \neq g_{c}$ 时，也能在图77所示的 QC 区域内观察到标度行为（并带有一定的修正，修正大小取决于所考察的量）。偏离 $g_{c}$ 后，体系会交叉进入 RC 或 QD 行为。这里我们只考察两个重要的物理量；除了关联长度，我们还分析均匀磁化率，这或许是在偏离 $g_{c}$ 后 QC 行为最稳健的量。预测的形式是，在 $g_{c}$ 点处磁化率随 $T$ 线性变化，偏离 $g_{c}$ 后则有一常数偏移：

其中 $c$ 是自旋波速度。常数 $a$ 并非精确已知，但基于 $1/N$ 展开中主导项计算出的值被认为接近实际值 [84]。线性行为应仅在 $T \to 0$ 时严格适用，因为存在 $T$ 的更高次幂（带有未知前因子），而这些项并未包含在(293)中。

我们首先研究在 $g_{c}$ 处（使用 $g_{c} = 1.9095$，在如上确定的临界点误差范围内）磁化率的有限尺寸收敛性，以确保我们能够可靠地获得反映低温下热力学极限的结果。图79左图显示了多个2的幂次系统尺寸（直至

图 79. 在临界点 $(g = 1.9095)$ 处，使用不同晶格尺寸计算得到的二聚化系统磁化率除以温度。左图中在小系统上观察到的峰是一个有限尺寸特征，该峰随 $1/L$ 向 $T = 0$ 移动。右图中的实线是对 $L = 512$ 数据 $( T / J_{1} \leq 0.20 )$ 的立方多项式拟合，该拟合给出式(293)中的前因子 $a / c^{2} = \bar{0.0916} ( 1 )$。作为比较，右图中也显示了 $g = 1.9090$ 时尺寸收敛的结果（最低的曲线，标记为 +）。在此温度范围内，在该耦合处也观察到了量子临界行为，该耦合点略处于相变的非磁性侧（因此最终当 $T \to 0$ 时 $\chi / T \to 0$）。

$L = 256$）的 $\chi / T$ 的温度依赖性（画图观察 $\chi / T$ 比 $\chi$ 本身更容易）。正如预期，结果在高温下收敛迅速，有限尺寸效应大约在 $T \sim 1 / L$ 时开始出现（由于动力学指数 $z = 1$，这在意料之中）。图79的右图以更细致的尺度展示了低温结果，其中也包括了 $L = 512$ 的结果。基于对不同尺寸结果的比较，可以得出结论：使用 $L \leq 512$ 的晶格，可以研究 $T / J_{1} = 0.03$ 以下的热力学极限（甚至可能更低一些）。

图79的右图还展示了低温下对 $\chi / T$ 的立方多项式拟合。对渐近线性形式(293)的修正相当显著，导致 $\chi / T$ 从 $T / J_{1} \approx 0.3$ 附近的最小值增加到最终的 $T \to 0$ 值，增幅约为 $8\%$ [这可以看作(293)中有效的温度依赖前因子 $a / c^{2}$]。原则上，外推得到的 $T = 0$ 值 $a / c^{2} \approx 0.0916$ 可用于提取自旋波速度，但由于前因子 $a$ 并非精确已知 [84]，这个估计可能不太精确。此处获得的结果的另一个用途是，将其与通过其他方式获得的 $c$ 的精确估计结合使用，这将有助于检验前因子的近似计算结果 [84]。然而，这超出了本文讨论的范围。

图80(a)展示了在$g = g_{c}$以及临界点两侧一定距离的两个值处$\chi$的温度依赖性。在$T \approx J_{1}$处有一个宽峰，对应于向最终高温居里行为的交叉。在峰值以下，这些近临界系统均表现出近似T线性的磁化率，与公式(293)形式一致，其中对于$g > g_{c}$有$b < 0$，对于$g < g_{c}$有$b > 0$。在更低的温度下，会出现向RC或QD行为的交叉，但由于温度依然过高，在此处并不明显。在QD区域，对于低于能隙∆的温度，$\chi \to 0$呈指数形式。RC形式也与QC形式一样是线性的，但斜率在交叉区域内发生变化。这种交叉实际上可以在标准的二维海森堡模型中看到，对应$g = 1$，它可能看起来离临界点太远了。尽管如此，在图73的左图中，仍然可以在$T / J \in \left[ 0.3 \sim 0.5 \right]$范围内看到$\chi$的近似线性行为，然后才在更低的温度下进入RC行为。该窄窗口内的$\chi$斜率与由已知自旋波速度以及(79)式中近似计算的前因子a得到的斜率非常吻合，这表明该行为实际上是远离量子临界点的量子临界行为的一种体现[84, 225]（尽管这种高度吻合可能在某种程度上是偶然的，因为正如我们上面所总结的那样，在刚好$g_{c}$处，存在对纯线性形式的重要修正，其温度远低于$g = 1$时的QC窗口）。

图80. 临界$( g = 1.9095 )$及近临界系统中磁化率（左）和关联长度（右）的温度依赖性。实线显示了$\chi$的渐近量子临界$T$线性行为（来自图79中拟合的$T = 0$值）以及$\chi_{x}$的$1/T$线性行为[遵循带有小负常数修正的形式(292)]。

图80(b)展示了在图80(a)相同近临界耦合下x方向的关联长度（在$g_{c}$处，该值比y方向关联长度大约30%）。虽然对于渐近形式(292)（其中$z = 1$）来说，行为非常线性，并带有一个小的常数修正，但略微偏离临界点的系统的结果在$T / J_{1} \approx 0.3$以下显著偏离线性。由于对于$g \le g_{c}$，关联长度是一个发散量，因此其对临界形式的修正远大于不发散的均匀磁化率。

5.3.5 J-Q 模型中的Néel-VBS相变#

通过使用第2.4.3节中介绍的J-Q模型，可以研究量子相变，在这类相变中，不仅反铁磁长程序消失，而且在价键固体形成时，非磁性态中还会有不同的对称性被破缺。表面上看，价键固体可能与二聚化海森堡模型的非磁性态颇为相似，因为在这两种情况下，系统都会展现出由强关联和弱关联最近邻对构成的图案。然而，这些态在根本上是不同的，因为在“人为”二聚化的海森堡模型中，哈密顿量本身破缺了平移对称性，并且强关联键也对应于强耦合的自旋对。相比之下，在J-Q模型中，哈密顿量遵守晶格的所有对称性，而平移对称性是在价键固体态中被自发破缺的。因此，这种价键固体的量子涨落，尤其是在相变点附近，是非常不同的（并且有趣得多）。

在我们将要考虑的四方晶格上，以及在我们选择的相互作用类型下，VBS可以形成四种等价的图案；因此，破缺的对称性是$Z_{4}$。我们将首先讨论具有四自旋耦合的标准J-Q模型[17]（此后术语J-Q模型将特指这一特定情况）。在这种情况下，Néel-VBS相变似乎是连续的，并展现出所提出的[16, 235]（且仍有争议的[236]）退禁闭量子临界点的许多特征。此情况下的VBS态很可能是柱状的，但该态的确切性质被VBS态中涌现的U(1)对称性所掩盖（在该模型中，VBS态总是近临界的，具有VBS序参量的大幅涨落）。这是所谓的退禁闭量子临界点的显著特征之一。尽管J-Q模型的行为在许多方面与该理论一致，但Néel-VBS相变也有一些方面仍未得到解释——尽管这未必与理论相冲突，因为详细的分析计算非常困难[235, 237]，并且基于描述该相变提出的场论，SU(2)自旋系统中退禁闭量子临界点的许多性质根本无法定量知晓。我们将在这里考察量子临界标度行为的一个强修正（可能为对数形式）的例子，退禁闭量子临界点理论并未预言这种修正。

我们还将讨论另一种多自旋耦合的一些结果，这种耦合会导致交错的VBS图案[如图4(c)中的二聚体图案]。这种情况下的相变显然是一级的[109]。我们将讨论两个模型中相变性质截然不同的原因。

J-Q模型的SSE实现。我们首先讨论使用SSE方法研究J-Q模型的一些实现问题。该模型最初是使用一种在价键基（我们将在第6节简要讨论）中构建的投影QMC技术进行研究的。然而，为其实现SSE方法也很容易[238]，最初的实现使用了“有向环”更新（一种广义的环算法，我们也将于第6节总结）。也可以设计一种J-Q模型的SSE算法，它与各向同性海森堡模型使用算符环的标准算法几乎相同。为了说明这一点，我们可以将具有四自旋相互作用的J-Q哈密顿量，即式(23)，用海森堡模型SSE算法中使用的对角和非对角键算符(255)和(256)表示：

与海森堡模型类似，[b] 表示一条含两个相互作用自旋的键 $[ i ( b ) , j ( b ) ]$ ，而 [bc] 表示两条平行键 $[ i ( b ) , j ( b ) ; k ( c ) , l ( c ) ]$ ，对应于 (23) 中的单重态投影算符乘积 $S_{i j} S_{k l}$ 。我们将这些键按图8所示排列，(294)中的求和遍历晶格上的所有平移，并同时包含水平和垂直键方向。然而，该方案独立于单重态投影算符的排列方式，且推广到 Q 项中任意数量的键也是平凡的。唯一的约束是我们必须避免符号问题，只要 $J > 0$ 且 $Q > 0$（考虑到(294)中的负号，这在能量上有利于 J 和 Q 项所包含键上的单重态），就能满足这一条件。关于不存在符号问题的讨论基于[238]中的子晶格旋转，也可使用第5.2节中用于二分海森堡模型的简单算符计数论证来证明。

图81. J-Q 模型中顶点的一些例子。仅展示了与算符对应的部分允许的顶点腿自旋态（空心和实心圆圈）。空心和实心条分别表示对角和非对角键算符。J 顶点与海森堡模型中的相同。式(294)中的 Q 项是两个键算符的乘积，展开后包含对角和非对角因子的所有组合。允许的环路仅通过 Q 顶点情况下的其中一个算符因子，如图所示，环路段在所有允许的腿对处均有展示。翻转环路可导致允许的算符和自旋态的任意组合。

我们现在有带四条腿的 J 顶点以及带八条腿的 Q 顶点，如图81所示。Q 顶点可视为两个 J 顶点连接而成，包含了来自(294)中 Q 项的四个算符所产生的对角和非对角部分的所有可能组合。显然，我们可以沿用处理海森堡模型的方法，通过结合对角更新和环路更新来更新算符字符串和存储的状态。这里的关键仍然在于所有 J 和 Q 顶点的矩阵元相同（值分别为 $J / 2$ 和 $Q / 4$），这意味着不改变顶点类型（J 或 Q）的环路更新总是可以接受的。对于 Q 顶点，满足此约束的环路在同一算符因子处进入和退出，如图81中的环路段所示。

在对角更新中，插入和移除 J 和 Q 对角算符。插入对角算符的最简单方法是从(294)中所有可能的单键 [b] 和双键 [bc] 实例中完全随机选择。水平和垂直键 [b] 各有 N 个，水平和垂直键对 [bc] 也各有 N 个，总共可供选择的有 4N 种情况。当前状态中，被选算符作用的键或键对上的自旋必须反平行，若满足此条件，接受概率是对式(265)的简单修改，用 $\beta_{J} = J / ( 2 T )$ 或 $\beta_{Q} = Q / ( 4 T )$ 替换 $\beta$，具体取决于插入算符的类型。键数 $N_{b}$ 被替换为键和键对的总数，即 4N。相同的修改也适用于移除概率 (265)。如果比值 $Q / J$ 远偏离于1（在 VBS 态和相变点处正是如此），最好在生成键或键对时就已考虑这一比值（但即使是朴素的随机算符生成方式实际上也运行得很好）。注意，在我们感兴趣的参数范围内，Q 远大于 J，因此最好以 Q 为单位定义温度，即在程序中设置 $Q = 1$。

图82. J-Q模型交错磁化强度的Binder累积量。左图展示了$J / Q = 1$以下的结果，其中$U_{2}$趋近于1，表明存在奈尔态。对于较小的$J / Q$，累积量随系统尺寸增大而消失，右图更详细地展示了这一点，证明当$J / Q < 0.045$（低于交叉点，这些交叉点在临界点处汇聚）时存在非磁态。对于所有$J / Q$值，累积量保持正值，表明不存在相共存，即基于这些结果，该相变是连续的。

环路更新与我们为海森堡模型开发的更新方法相同，只是简单扩展了顶点链表（在海森堡模型的图57中有所说明），现在每个顶点有八个元素，而海森堡模型情况下只有四个。尽管J-顶点只有四条腿，最实用的做法是为所有顶点分配八个存储元素，对于J-顶点仅使用前四个（并将未使用的元素填充为负数，这样在环路更新中永远不会访问到它们）。

J-Q模型中的连续量子相变。对J-Q模型的初步研究[17]证实了对于小$J / Q$值存在VBS基态，且存在一个符合$z = 1$量子临界点的Néel–VBS相变。这些结果很快在$T > 0$量子临界区[238, 239]的SSE研究中得到确认。后来的一项世界线QMC研究声称存在一级相变的证据[240]，其他体系的研究[236]也支持了这种观点。然而，近期在非常大晶格（最大$L = 256$）上的计算并未发现任何一级行为的迹象[110]。这违反了“朗道”规则，根据该规则，破坏不相关对称性的两个态之间的直接相变通常应该是不连续的，但这与退禁闭量子临界点理论[16]相符。然而，在某些量中发现了显著且意外的标度修正（这在[240]中可能被误认为是一级相变的迹象）。这里我们首先讨论支持单一连续Néel–VBS相变的一些证据。

Binder累积量是分析连续相变和一级相变的一个有用量。正如我们在第3.4节中，通过图20的例子所讨论的，一级相变处的相共存会导致累积量在相变点出现负的发散峰。图82展示了使用SSE方法，在$\bar{\beta = {Q} / {T}} = L$（基于先前表明$z = 1$的工作）条件下计算的J-Q模型的反铁磁Binder累积量$U_{2}$的结果。对于$J / Q > 0.05$，$U_{2}$随着L增大而增大并趋向于1，表明系统具有长程奈尔序。对于较小的耦合比，累积量减小到0，没有出现负峰迹象。交叉点行为良好，表明临界点位于$( J / Q )_{c} \approx 0.045$。

图83. J-Q模型自旋关联长度（左图）和二聚体关联长度（右图）的有限尺寸标度。注意，$\xi_{s}$和$\xi_{d}$的交叉点随着晶格尺寸增大而向不同方向移动，相互靠近。

接下来，我们检查关联长度，该长度使用结构因子定义（70）对自旋-自旋和二聚体-二聚体关联进行计算。自旋（奈尔）序参量的波矢为$( \pi , \pi )$，而柱状VBS对应于$( \pi , 0 )$或$( 0 , \pi )$，后者被用作定义（70）中的参考点$q = 0$。为了计算二聚体关联，我们使用定义（25），但仅包含键算符（26）中的对角关联，即

尽管基于旋转对称定义和$z$分量定义的关联函数的有限尺寸标度性质可能不同[238]，但从这些函数提取的关联长度在临界点处仍应都表现出有限尺寸标度$\sim L$（直至标度修正项）。

我们可以使用标度关联长度$\xi_{s} / L$（自旋）和$\xi_{d} / L$（二聚体）来检验两个序参量是否在同一点消失，或者是否存在可能两个不同的相变。图83显示了原始的SSE数据，这些数据同样是在$\beta = L$条件下生成的。不同系统尺寸的曲线相互交叉，尽管即使对于最大的晶格，也存在一些水平方向和垂直方向的漂移。虽然小晶格下$\xi_{s} / L$和$\xi_{d} / L$的$J / Q$交叉值相距甚远，但随着系统尺寸增大，它们彼此靠近，其方式与单一临界点一致。

图84展示了从L和2L数据中提取的若干不同量（上述讨论的量以及自旋刚度）的交叉点。将这张图与二聚化海森堡模型对应的图76进行比较是很有趣的。显然，J-Q模型的临界点收敛速度要慢得多。虽然所有量在一定系统尺寸范围内都几乎与$1 / L$呈线性关系，但对于最大系统，行为开始发生变化，这暗示最终将以$L$的某种更高次幂收敛到共同的临界点$( J / Q )_{c} \approx 0.045$。回想一下，对于二聚化海森堡模型，行为呈$\sim 1 / L^{\omega}$，其中$\omega$在$2 \sim 2.5$范围内，并且这种行为在相当小的晶格$( L \approx 8 )$上就已适用。由于收敛缓慢，很难基于这些单独的交叉点估计值来精确估计J-Q模型的临界点。然而，图84中有一个显著特征：$\xi_{s} / L$和$\xi_{d} / L$的交叉点形成的曲线看起来非常对称，以相同速率但从相反两侧趋近于表观渐近值。因此，这两个估计值的平均值几乎完全不依赖于系统尺寸（如图中所示）。基于这些平均值，人们可以获得比最初预期好得多的临界点估计。基于四个最大L点（它们在统计误差内完全一致）的结果是$( Q / J )_{c} = 0.04498 ( 3 )$。

图84. 基于L和2L系统的曲线交叉，使用多个不同量——Binder累积量$U_{2}$、尺寸标度自旋刚度$\rho_{s} L$以及自旋和二聚体关联长度$\xi_{s} / L$和$\xi_{d} / L$——得到的尺寸依赖临界耦合估计值。$\xi_{s} / L$和$\xi_{d} / L$估计值的平均值（以实心方块显示）几乎不表现出尺寸依赖性。这些点在四个最大L处的平均值给出了临界点估计值$( J / Q )_{c} = 0 . {\bar{0}} 4498 ( 3 )$，如水平线所示。

交点缓慢收敛表明底层物理量中存在大的标度修正。值得注意的是，在图84分析的各量中，Binder累积量的交点具有最弱的尺寸依赖性，这点在对比图82（累积量）和图83（关联长度）的原始数据时也很明显。由于Binder累积量涉及两个相似量的比值，这些量中的异常标度修正可能会部分相互抵消。两个关联长度交点中的有限尺寸修正几乎完全相互抵消的原因完全不清楚。这可能表明近临界的VBS态和Néel态之间存在某种隐藏的对称性——或许是退禁闭量子临界性固有的某种对偶性[这适用于理论的U(1)版本[16]，但在相关的SU(2)变体中并不明确，而后者恰是J-Q模型相关的版本]。注意，不仅是交点位置具有对称性，事实上图83中的关联长度曲线本身看起来也相当对称。将其中一条曲线关于临界点反射，$J / Q - ( \bar{J} / Q )_{c} \to - [ J / \bar{Q} - ( J / \bar{Q} )_{c} ]$，最大系统的自旋和二聚体关联长度几乎完全重合，即，这些关联长度的标度函数本身看起来也是相同的。

图85. J-Q模型的均匀磁化率$\chi / T$，在临界值$( J / Q )_{c} =$ $0.04498 ( 3 )$附近三个耦合比下，分别以线性（左）和对数（右）温度尺度作图。实线是在$J / Q = 0.045$（该值在临界值的误差范围内）处对形式$\chi / T = a + b \overline{{\ln}} ( Q / T ) + c T^{2}$的拟合。SSE计算使用了足够大的晶格（最低温度下晶格大小达$L = 512$），以获得代表热力学极限的结果。

$T > 0$时的量子临界磁化率标度。在常规量子临界点处，人们会预期均匀磁化率具有渐近线性温度依赖性，正如我们在图79的二聚化海森堡模型中所见，那里$\chi / T$比值中也可见到更高次T幂形式的修正。在J-Q模型中，对预期量子临界形式的修正要强得多。如图85所示，似乎存在一个乘法对数修正，即渐近$T \to 0$临界行为可能具有形式$\chi \sim T \ln ( 1 / T )$。结果分别在线性和对数温度尺度上作图，以清晰传达其行为如何不同于图79中二聚化模型的行为。J-Q模型中的行为似乎与常规的高次幂律修正不相容，而对数项却很好地描述了该行为。

根据当前的退禁闭量子临界性理论，并未提出对数修正。所有在非紧致 $\mathrm{C P^{1}}$ 场论框架内的解析计算，都是通过将该系统推广至 SU(N) 并对其大 N 展开来进行的（该场论被推广至 $\mathrm{C} \dot{\mathsf{P}}^{N - 1}$，其中 $N = 2$ 是自旋分量的物理数目）[16, 16, 237]。对于 $N = 2$ 的情况，可能会出现不同的特征。另一方面，对该场论格点版本的蒙特卡洛计算[241, 242, 236]并未指出任何反常的标度修正（尽管有的计算只发现了一级相变[236]）。由于这些计算彼此不一致，$N = 2$ 理论的情况仍然悬而未决。在相关的费米子规范场论中，磁化率及其他量的对数修正确实会出现[243]。如我们在第4.3节和第5.3.1节所讨论的，在场论语言中，对数修正可以作为边际算符的结果出现，但其他解释也是可能的。原则上，这些修正可能并非真正的对数形式——它们也可能源于近似边际（但最终不相关）的算符。相应的、具有非常小指数的幂律修正，在可达到的系统尺寸下可能与对数形式无法区分。虽然这种情形对某些量是可能的，但对于均匀磁化率来说，这将非常不寻常，因为通常二维反铁磁体中所有对量子临界行为的修正都是 $T$ 的整数幂次[84]。

图86. 展示了价键对的局域涨落如何导致柱状态（a）逐渐变为方块态（d）。在（b）和（c）中，少量键对被翻转。当方块（键叠加）的排列如（d）所示时[四种等价的平移图案之一]，所有键对的最大数量的局域水平和垂直涨落可以发生。这些方块上的垂直键对的比例 $p \in [ 0 , 1 ]$ 对应于一个VBS角 $\phi \in [ 0 , \pi ]$，该角度使用 x 和 y 二聚体序参量 (296) 定义。在这些图中，角度 $\phi \in [ \pi , 2 \pi ]$ 对应于将垂直键向上移动一步。

VBS中的涌现 U(1) 对称性。退禁闭量子临界性理论的一个显著特征是，除了正常的关联长度 $\xi$（在奈尔态中可以取为二聚体关联长度 $\xi_{d}$，在VBS态中可以取为自旋关联长度 $\xi_{s}$）之外，在VBS相中还存在另一个长度尺度 $\Lambda$。这是激发态 $S = 1 / 2$ 自旋子（它们是VBS中的缺陷，一种涡旋[245]）被禁闭为 $S = 1$ “三子”时的距离尺度（即 $\Lambda$ 是自旋子束缚态的大小）。相同的长度尺度也应与一种涌现的 U(1) 对称性相关联，这种对称性联系着一类连接两种不同VBS有序图案（柱状态和方块态）的特殊量子涨落。这些“角度”VBS涨落以及两类VBS态如图86所示。角度涨落在小于 $\Lambda$ 的长度尺度上发生（即在小于 $\Lambda$ 的长度尺度上，系统处于两种不同VBS态的叠加态）。激发态的自旋子可以彼此自由移动到距离 $\Lambda$ 之外，但在更大的距离上，它们由于VBS长程序所介导的有效势而变得禁闭（在这方面，可以注意到与粒子物理学中某些规范理论中禁闭/退禁闭的有趣联系[18]）。当接近临界点时，$\Lambda$ 应作为关联长度的幂次发散，即 $\Lambda \sim \xi^{a}$，其中指数 $a > 1$。

VBS序参量中近似 U(1) 对称性的存在，可以通过检查对应于具有 x 和 y 方向键的柱状序的序参量 $D_{x}$ 和 $D_{y}$ 的分布 $P ( D_{x} , D_{y} )$ 来检验，定义如下：

对于给定的蒙特卡洛构型，可以在固定虚时（SSE中的传播态）下，使用第5.2.5节讨论的改进估计量（improved estimator）来计算关联函数，从而评估 $D_{x}$ 和 $D_{y}$ 。该改进估计量是自旋旋转不变的，通过考虑环构型与价键基（valence bond basis）[244, 31, 50] 之间的等价性（我们将在第6节简要讨论），可以最好地理解这一点。事实上，只要温度足够低以至于模拟仅采样单重态扇区（singlet sector），等时关联函数的改进环估计量就完全等价于算符在特定价键态中的期望值。

图87. 对于 $L = 128$ 系统，在 $J / Q = 0$ (左) 和 0.045 (右) 下计算得到的 VBS 序参量分布 $P ( D_{x} , D_{y} )$ 。直方图显示在区域 $| D_{x} | , | D_{x} | \leq D_{\mathrm{cut}}$ 内，其中 $D_{\mathrm{cut}}$ 针对两个耦合比选择以覆盖权重显著的区域。更亮的特征对应更高的概率密度。这里的绝对标度并不重要，重要的是分布的圆形形状。由于模拟难以使大系统的 VBS 角度发生涨落，这些直方图中的权重并非完全均匀分布[既不服从完美的 $Z_{4}$ 对称性，也不服从 U(1) 对称性]。

虽然 $P ( D_{x} , D_{y} )$ 的计算可以采用 SSE 方法，但下面将要讨论的结果实际上是通过价键基（valence bond basis）中的基态方法生成的（在其最新表述 [50] 中，该方法实际上与 SSE 方法非常相似，如第6节所述）。这里并不严格需要旋转不变的估计量 [109]，也可以仅使用 (296) 式中的自旋 z 分量，正如我们在关联函数方程 (295) 中所做的那样。然而，在价键基中，自然要计算旋转不变的量。请注意，$P ( D_{x} , D_{y} )$ 不是一个算符期望值（即，它不是正常的物理可观测量），而仅仅是在特定基下进行的模拟中单个测量的分布。然而，其对称性性质反映了系统 VBS 序参量的实际对称性。

我们已经看到了 $P ( D_{x} , D_{y} )$ 的一个经典类比：阻挫伊辛模型中条带序参量（stripe order parameter）的分布 $P ( M_{x} , M_{y} )$ ，其结果如图19所示。在这种情况下，序参量在转变温度及其以下时的四重对称性很明显，而到非常高的温度（图中未显示）时，由于系统不同部分短程序的独立高斯涨落，分布变成具有完全旋转对称性[即 U(1) 对称性]的中心峰。类似地，在 VBS 中，可以预期 $P ( D_{x} , D_{y} )$ 在奈尔态（仅存在短程 VBS 关联）中具有 U(1) 对称性，并在 VBS 态内部变为四重 $( Z_{4} )$ 对称性。然而，从图87（显示了比文献 [17, 240, 108] 中更大系统的近期结果）可以看出，实际情况并非如此。在 VBS 相内部，$J / Q = 0$ 时，分布呈现环状，中心没有权重，也没有四重对称性的迹象。这表明 VBS 序参量的幅度 $D$（$D^{2} = D_{x}^{2} + D_{y}^{2}$）已经形成，$\langle D^{2} \rangle > 0$，但其角度对称性尚未建立；也就是说，尽管 VBS 序参量的幅度涨落相对较小，但单个分量 $D_{x}$ 和 $D_{y}$ [它们定义了 $( D_{x} , D_{y} )$ 平面中 VBS 的角度 $\phi$] 的涨落却非常大。在靠近临界点处，如图中的 $J / Q = 0.045$，分布如预期那样在中心达到峰值，但同样不存在 $Z_{4}$ 对称性。

在连续场论中，负责形成 VBS 的算符是危险无关的 [16, 245, 246]，这意味着尽管它是一个具有 $Z_{4}$ 对称性的微扰，并导致 $Z_{4}$ 对称性破缺的 VBS 态，但这个微扰并不影响临界点。退禁闭量子临界点理论具有固有的 U(1) 对称性，并且在 VBS 算符存在的情况下，这种对称性保持不变。J-Q 模型没有显式的 U(1) 对称性——这里固有的反而是正方晶格离散的 $Z_{4}$ 旋转对称性，并且人们通常期望在该晶格上形成的 VBS 会表现出 $Z_{4}$ 对称性，正如第 3.4 节讨论的伊辛模型的条纹序一样。相反，我们看到了一个 U(1) 对称的序参量分布 $P ( D_{x} , D_{y} )$，根据退禁闭量子临界点理论，如果系统小于退禁闭长度尺度 $\Lambda$，就应该会看到这种现象（因此，根据图 87 的结果，即使在 $J / Q = 0$ 时，$\Lambda$ 也必须很大，$\Lambda > 128$）。如果增加系统尺寸，我们最终应该会在 $L > \Lambda$ 时观察到四重对称分布（已在一些相关系统中看到 [108, 247, 248]）。柱状态对应于在 (0, D) 和 ( D, 0) 处的峰值。

如上所述并在图 86 中说明，涌现的 U(1) 对称性是柱状和 plaquette VBS 模式之间局部涨落的结果 [245]，其方式使得所有与幅度 D（弱涨落）兼容的 $D_{x}$ 和 $D_{y}$ 可能值都被同等（或几乎同等——完美的 U(1) 对称性仅适用于临界点）采样。这与基础量子力学中的一个类比相似：一个存在于略微不对称双势阱中的粒子。如果势阱之间的势垒非常高，那么粒子就局域在较深的势阱中；但如果势垒不是很高（相对于动能而言），那么粒子就会在两个势阱之间涨落，并且波函数在势垒内部也有权重。在 VBS 中，势阱对应于柱状态和 plaquette 态，其角度分别为 $n \pi / 2$（柱状）和 $n \pi / 2 + \pi / 4$（plaquette），其中 $n = 0 , 1 , 2 , 3$，它们之间的涨落对应于在这些点之间的连续角度上出现非零概率。在退禁闭理论中，柱状和 plaquette VBS 态几乎是简并的，并在临界点处变为精确简并。当接近临界点时，它们之间的势垒也会减小。在双势阱类比中，势阱变得更加简并，并且它们之间的势垒减小。那么，在有限系统中（直到 $L \sim \Lambda$），柱状态和 plaquette 态之间将存在显著的涨落，这意味着 $P ( D_{x} , D_{y} )$ 在所有角度上都有权重，并且当 $L > \Lambda$ 时，对称性会从 U(1) 交叉到 $Z_{4}$。这同样适用于在无限大（或非常大）系统中长度为 L 的子系上计算的序参量（这已在具有 q 重序参量并表现出 U(1) 到 $Z_{q}$ 对称性交叉的经典模型中研究过 [249, 250]）。

由于在迄今为止研究过的最大晶格中，J-Q 模型尚未清晰地展现出 $Z_{4}$ 对称性，因此关于 VBS 本质的问题仍然悬而未决。然而，在一个 Q 项由三个平行单重态投影算符乘积组成的 J-Q 模型中，VBS 序更加稳健，并且可以清晰地看到对应于柱状态的 $Z_{4}$ 对称分布 [108]。对标准 J-Q 模型的 SU(N) 推广（$N = 3$ 和 4）进行的模拟也显示出了柱状 VBS 基态 [108]。最有可能的情况是，标准的 SU(2) J-Q 模型也具有柱状态，并且在近期（正在进行的）针对 $L = 64$ 的长程模拟中确实有一些微弱的迹象表明这一点。

对于系统尺寸 $L \gg \Lambda$，VBS分布应具有四重对称性，通过研究从U(1)对称性到 $Z_{4}$ 对称性的转变过程，原则上可以确定将退禁闭长度标度Λ与关联长度 $\xi$ 联系起来的指数 a。经典交叉现象的实例（这也启发了理论[16]）在文献[249, 251, 250, 252]中得到了研究。一种针对对 $Z_{4}$ 对称性敏感的量的有限尺寸标度方法可用于提取指数 a [252]。在标准J-Q模型中，尚无法进行此类分析，因为在当前可访问的晶格中，向 $Z_{4}$ 对称序参量的交叉现象不够清晰。然而，该分析已在具有柱状六自旋相互作用的J-Q模型中完成，结果表明 $a \approx 1.2$。

其他临界指数。除了动力学指数 $z = 1$ 之外，其他临界指数已在几项QMC研究[17, 238, 239, 108, 110]中进行了分析。最值得注意的结果是 $\eta_{s}$（临界自旋-自旋关联函数的指数）异常大，$\eta_{s} \approx 0.3$，这与退禁闭量子临界点理论定性一致。基于 $\mathrm{CP}^{\dot{N} - 1}$ 场论 $1 / N$ 展开中的 $N = \infty$ 值 $\eta_{s} = 1$（即不同于传统的平均场值 $\eta_{s} = 0$），论证了 $\mathbf{A} \ \mathbf{\tilde{\Sigma}} \mathrm{large}^{3} \ \eta_{s}$ [16]。$1 / N$ 修正难以计算，但总体上支持异常大的指数[237]（尽管无法可靠地将结果扩展到 $N = 2$）。在SU(N)系统[108]的QMC计算中得到的二聚体-二聚体指数 $\eta_{d}$ 的行为，似乎遵循了理论[246, 253]中预测的随N的标度行为（其中 $\eta_{d}$ 随N发散）。

鉴于在图85中分析的磁化率以及自旋刚度[110]中存在显著的标度修正（可能为对数形式），一个重要的问题是计算出的指数是否也受到影响。到目前为止，这一点并不明显，因为使用合理的次主导指数值的标准数据塌缩方法效果很好[108]。目前基于 $T > 0$ [238] 和 $T = 0$ [108] 的QMC计算对 $\eta_{s}$ 的最佳估计值彼此吻合。另一方面，此处及文献[110]中讨论的结果将临界点置于比以前略高的 $J / Q$ 处，并且包含对数修正后分析得到的关联长度指数[110]也显著小于其他估计值[238, 108]。因此，最好使用几个不同的量和更大的晶格重复对所有指数的研究（这是一项正在进行的工作）。

从对杂质响应的研究中也发现了J-Q模型中存在对数修正的证据。通常，在二聚化模型的传统量子临界点处，向系统中引入一个空位并研究由此产生的磁化不平衡的空间分布，可以观察到该分布随系统尺寸的标度行为[254, 255]。在J-Q模型中，这种标度分析再次需要引入一个乘性的对数修正[255]。

对其他方法在J-Q模型中的应用。J-Q模型令人感兴趣且令人困惑的一个方面是，似乎很难用分析性多体技术正确捕捉其性质，即使是那些能够相当准确地定位二维阻挫$\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$海森堡模型中Néel-VBS相变的方法也是如此。从柱状二聚体态出发的均场处理得到的临界$J / Q$值与QMC结果相去甚远[256]。这种方法可以通过更好地考虑plaquette上的局部涨落来改进，这改善了临界点的值，但似乎导致了一级相变[256]。团簇均场计算随团簇尺寸的收敛性非常差[257]。这些结果表明存在异常强的非局域量子涨落[256]，这很难用从小团簇出发的局域方法或围绕固定二聚体图案的常规涨落来捕捉。造成这些难以捕捉VBS态的原因可能是涌现的U(1)对称性，这使得同时获得正确的长程行为（很可能是柱状序）以及更短长度尺度上柱状序和plaquette序之间的强涨落变得困难。

图88. 在此处研究的$\mathbf{J} {\cdot} \mathbf{Q}_{3}$模型中，三个单重态投影算符按交错图案排列。包含了投影算符乘积的所有不同取向（如图所示）和平移。

交错J-Q模型中的一级相变。测试涌现U(1)对称性与连续Néel-VBS相变之间联系的一种方法是构建一个模型，其中负责旋转VBS角的局域涨落被抑制。实现这一点的直觉来自Rokhsar-Kivelson (RK) 量子二聚体模型[258, 259]，该模型可以被视为一个由短价键主导的极端非磁性系统的有效模型（其中内部单重态结构也被忽略——键构型被视为正交态）。正方晶格上的RK哈密顿量可以写为$H_{\mathrm{RK}} = \nu V - k K$，其中$V$是对角（势能）算符，计算可翻转plaquette【平行键对，与图86(b,c)中的完全一样】的数量，而$K$是一个非对角（动能）项，用于翻转这样一个键对。该模型在$k = \nu$处有一个临界点，它将$\nu < k$时的plaquette VBS态【类似于图86(d)中的态】和$\nu > k$时的交错VBS态【键图案完全如图4(c)所示】分开。当$\nu \to k^{-}$时，plaquette态被量子涨落连续破坏，而交错态（有四个对称性相关的等价态）没有涨落，因为它没有可翻转的plaquette。因此，在$\nu \to k^{+}$时的相变是一级的。

RK模型的这个简单图像表明，实际自旋模型中的交错VBS也应具有强烈抑制的局域涨落，因此不应与涌现U(1)对称性相关联。由于局域涨落（以及由此的大尺度涨落）被抑制，它和奈尔态之间的相变应该是一级的。然而，这个图像并不完整，因为考虑到一个合理的自旋哈密顿量将远非二聚体模型，并且构成基态的价键构型也应包含一部分更长的键（存在这些键时，始终可能存在一些局域涨落），在交错VBS态中显然必须存在一些涨落。因此，值得在基态可以从奈尔反铁磁体调到交错VBS的模型的模拟中检验这一图像。通过在J-Q框架内以交错方式排列Q项的单重态投影算符（而不是在$2 \times 2$ plaquette内，或像[108]那样以三列或更多投影算符排列），可以轻松构造此类模型。然而，结果表明，由两个交错投影算符组成的Q项不足以破坏奈尔态。不过，使用三个按图88排列的投影算符（我们将其耦合称为$Q_{3}$），可以稳定一个交错VBS。我们在此简要总结近期对交错 $\mathbf{J} {\cdot} \mathbf{Q}_{3}$ 模型[109]进行SSE研究时发现的Neel-VBS一级相变证据。

图89. 交错 J $\mathrm{\partial}_{-} \mathrm{Q}_{3}$ 模型[109]中一级相变的证据。子晶格磁化强度的Binder累积量 (a) 呈现出一个随系统尺寸增大的负峰。(b) 观察到滞后效应：在耦合比缓慢增加和减少的模拟中，大系统会表现出不同的性质，能量在某个点上交叉，该点可视为依赖于尺寸的转变点。对不同系统尺寸的这些点进行外推，得到 $( Q_{3} / J )_{c} = 1.1933 ( 1 )$，在左图中以垂直虚线标出。所有结果均在逆温度 $Q_{3} / T = L$ 下通过SSE模拟获得。

一些结果如图89所示。奈尔序参量的Binder累积量呈现出一个随系统尺寸增大的负峰。正如第3.4节所讨论的（并且结果可与图20中经典阻挫伊辛模型的结果进行比较），这是一级相变的有力指示。对于 $L > 12$，数据变得非常嘈杂，因为模拟在两个共存态之间转换的效率低下。对于 $L > 24$，系统可能会完全陷入亚稳态。在VBS一侧接近相变处，若模拟在平衡过程中恰好进入类奈尔构型，它将在亚稳态奈尔态中停留很长时间（反之在相变另一侧亦然）。可以利用这种俘获现象以其他方式研究相变。通过在 $Q_{3} / J$ 紧密间隔的值上进行一系列模拟（从VBS或奈尔相内部深处开始，并且每次新模拟都从上次模拟的最后构型开始），可以跟踪基态进入相变“错误”一侧的相应亚稳态。从两次此类模拟（$Q_{3} / J$ 增加或减少）中计算出的能量如图89所示。两个分支相互交叉，这清楚地显示了相变的一级性质。能级交叉很可能是两个具有相同量子数的态之间避免的交叉，但这里的能级间距应该非常小。该间距应与共存点处两个态之间的隧穿时间成反比（但请注意，模拟时间与系统的实时动力学之间的关系尚不清楚）。

两个能量分支的交叉点可以看作临界耦合的一个有限尺寸定义。它随系统尺寸的变化大致为 $1 / L^{3}$，这与动态指数 $z = 1$ 的一级相变一致 [260]。外推得到的相变点为 $Q_{3} / J = 1.1933 (1)$。在相变点，VBS序参量约为其最大值（在完美的交错价键态中）的75%，这个值相当大，因此该相变被归类为强一级相变。