引言

这是量子多体计算理论的学习笔记，主要参考资料是A. W. Sandvik的《Computational Studies of Quantum Spin System》¹，主要内容涉及量子蒙特卡洛方法理论基础和一些重要算法的介绍。

量子蒙特卡洛方法和随机级数展开方法

经典统计物理核心的任务就是确定系统的配分函数，量子力学则进一步丰富了配分函数的内涵，引入了密度算符和密度矩阵的概念，所以在量子系统下一般有配分函数的形式，即对系统虚时玻尔兹曼算符求迹：

$$\begin{align} \mathcal{Z}=\mathrm{Tr}\{e^{-\beta H}\}. \end{align}\tag{1}$$

于是，基于虚时演化的世界线方法（World Line Methods）在自旋系统和相关的格子系统中被发展出来，该方法便涉及对玻尔兹曼算符的Suzuki-Troutter分解。而早在上世纪六十年代初，Handscomb就在海森堡铁磁系统发展出了基于系统玻尔兹曼算符的级数展开方法，并准确得到了置换算符乘积的迹。因此，由于随机级数展开形式的引入，量子蒙特卡洛计算的幂级数方法便重获新生。这个推广不仅和计算机算法的特点契合，也显示了离散幂级数和虚时演化下的路径积分的密切关系。在一些更为特殊自旋相关的体系，例如海森堡模型，自旋部分可以用高效的循环(loop)集群更新，而这就是对经典的Swendsen-Wang集群算法在量子系统下的推广。关于loop概念的推广，已经在一些量子自旋和玻色系统开展了很深入的研究。
除了上述应用，在自旋阻挫系统量子蒙卡一般会遇到所谓符号问题，即非正定的路径积分或级数展开。然而这样的结果原则上仍然可以计算得出，所以也有一系列工作集中在高温下符号涨落的控制问题以及低温下由于这种符号混合造成的统计误差问题。
以上这些问题和模型已经显示了丰富的物理，并给量子蒙卡的前沿应用和研究提供了深刻的思路。

虚时演化的世界线方法和量子蒙特卡洛计算基础

引入虚时的概念是量子统计物理的重要一环，经典的量子力学演化算符是薛定谔时间演化算符的形式：

$$\begin{align} \mathcal{O}_S=e^{-\mathrm{i}\frac{Ht}{\hbar}}. \end{align}\tag{2}$$

联系到Eq.(1)，考虑$t=-\mathrm{i}\Delta_{\tau}$，令$\Delta_{\tau}=\beta/L$便能得到玻尔兹曼算符的形式。Eq.(1)求迹的过程相应地有

$$\begin{align} \mathcal{Z}=\mathrm{Tr}\{e^{-\beta H}\}=\mathrm{Tr}\left\{ \prod_{l=1}^{L}e^{-\Delta_{\tau}H}\right\}. \end{align}\tag{3}$$

求迹和可以写成矩阵对角元求和的形式，配分函数则有

$$\begin{align} \mathcal{Z}=\sum_{\alpha_0}\sum_{\alpha_1}\cdots\sum_{\alpha_{L-1}}\langle\alpha_0|e^{-\Delta_{\tau}H}|\alpha_{L-1}\rangle\cdots\langle\alpha_2|e^{-\Delta_{\tau}H}|\alpha_1\rangle\langle\alpha_1|e^{-\Delta_{\tau}H}|\alpha_0\rangle. \end{align}\tag{4}$$

因此这里可以看做是一个$d$维的量子系统映射到$d+1$维，增加的一个维度正是虚时$\tau$，而系统虚时演化的过程则体现在量子态的演化，即所谓“路径”或者路径积分下世界线上的一个位形：$\alpha_0\to\alpha_1\cdots\to\alpha_{L-1}\to\alpha_0$（构成位形空间）。考虑对虚时演化算符泰勒展开至一阶，Eq.(4)则有

$$\begin{align} \mathcal{Z}\approx\sum_{\{\alpha\}}\langle\alpha_0|1-\Delta_{\tau}H|\alpha_{L-1}\rangle\cdots\langle\alpha_2|1-\Delta_{\tau}H|\alpha_1\rangle\langle\alpha_1|1-\Delta_{\tau}H|\alpha_0\rangle. \end{align}\tag{5}$$

这里$\{\alpha\}$表示所有态的演化步骤，也就是所谓“位形(Configuration)”。$L=\beta/\Delta_\tau$项因子每一项系数都会产生$\Delta_{\tau}^2$的误差，因此对于固定的温度因子$\beta$配分函数的精度是$\Delta_\tau$，而每个展开因子的主导可忽略项仍然包含$H^2$，于是可以猜测离散化过程的误差，其标度是$N^2$。不过实际上准确的结果并不是简单地考虑$\langle H^2\rangle\propto N^2$，主要是因为真实的路径和Eq.(5)取线性近似的路径并不一样。所以为了取得不依赖$\beta$和$N$的误差项，在线性近似下我们至少需要$N\beta$阶的时间切片，即$\Delta_{\tau}\propto1/N$。

蒙特卡洛世界线更新

对于玻色系统，在晶格模型下仅考虑动能项，其Hamiltonian可以写作

$$\begin{align} H=K=-\sum_{\langle i,j\rangle}K_{ij}=-\sum_{\langle i,j\rangle}(a_j^+a_i+a_i^+a_j). \end{align}\tag{6}$$

其中$a_i^+,a_i$分别是格点为$i=1,\cdots,N$的生成湮灭算符，$\langle i,j\rangle$表示晶体最近邻格点。如果格点的在位排斥很强的话，我们可以将格点占据数设为$n_i=0,1$，即所谓“硬核玻色子”，然后仍然可以将该玻色系统赋予$S=1/2$的量子自旋，格点占据态和空态直接对应 $\uparrow$ 和 $\downarrow$。在不考虑相互作用的情况下，系统Hamiltonian等价于XY模型，$H_{ij}=-(S_i^+S_j^-+S_i^-S_j^+)$。考虑到Eq.(6)不包含化学势，利用正则系综，考虑系统具有固定数量玻色子$N_B$，则系统的磁性可以表示为$m=N_B/N-1/2$。

选用占据数表象，利用Eq.(5)，可以得到矩阵元$\langle\alpha_{l+1}|1+\Delta_\tau K_{ij}|\alpha_l\rangle$，这相当于把Hamiltonian表示成所有单独hopping项求和的形式。对于那些非零矩阵元，态$|\alpha\rangle$和$\alpha_{l+1}$满足两种可能的关系，要么$|\alpha_{l+1}\rangle=|\alpha_l\rangle$，要么$|\alpha_{l+1}\rangle=K_{ij}|\alpha_l\rangle$。对于后一种情况，一定代表在格点$i$会生成粒子，这是hopping动能$K_{ij}$从原来的空格点$j$作用的结果，反之亦然。

Fig.1 三玻色子硬核系统14格位20时间（虚时）切片的离散路径积分示意图。该一维时空玻色链满足时间和空间的周期性边界条件，图中值得注意的是，在线性切片算符$1-\Delta H$作用下，在时间切片上不能有超过一个粒子的hop，而在空间维度全同粒子由于周期性边界条件形成缠绕的图像，显示了粒子虚时演化的非局域特征。

因为这个条件对所有连续矩阵元都必须满足，即遍历$l=0,\cdots,L-1,0$，这就要求每个粒子都沿着一条世界线，在时间切片上构成一个循环(Fig.1)。而且虚时演化的方向与Eq.5求迹的过程一致，不过这些边界条件可以应用在玻色态，而不是严格意义上的世界线上。因为玻色子本身是全同粒子，它可以在虚时任意切片位置置换，并且仍然满足所要求的周期性条件。所以可以推知，对于一维硬核链，仅有的置换一定是系统中所有粒子的循环置换，因为粒子只考虑最近邻hopping的话，要满足hard-core的占据数条件。这就导致一个净的环形的粒子流(Fig.1右图)。同理，除此之外的其他不同置换可以出现在“软核(Soft-core)”系统。一维硬核的结果在二维plane或者三维bulk甚至在考虑超出最近邻hopping的一维链仍然只有世界线缠绕的图像会产生净的粒子流。而世界线缠绕系统的净次数又叫做缠绕数(Winding Number)，这个量是个拓扑量，它标志了系统的一种非局域性质。在更高维度，空间取向仍然存在“缠绕”，有对应的“缠绕数”。玻色系统缠绕数对应的是超流，具体而言是系统缠绕数的归一化方差。在自旋系统，缠绕数则和自旋刚度(spin stiffness)或者螺旋模量(helicity modulus)有关.

可以用一个适当的世界线位形权重因子来表示系统的配分函数，即

$$\begin{align} \mathcal{Z}=\sum_{\{\alpha\}}W(\{\alpha\}). \end{align}\tag{7}$$

如果仅仅考虑Eq.5的情况，玻色算符的矩阵元在硬核体系下总是1（占据数对应的简并度为1），所以这个路径权重就变成

$$\begin{align} W(\{\alpha\})=\Delta_\tau^{n_K}. \end{align}\tag{8}$$

这里$n_K$表示为动能跳跃的总次数（这就是动能事件，$K_{ij}$满足之前的条件反映在世界线上就是一次跳跃，同时产生一个$\Delta_\tau$，对应一个时间切片。），对应Fig.1世界线对角部分。这个结果可以简单推广到软核系统，对应的生成湮灭算符的项因子分别变成$\sqrt{n_i},\sqrt{n_i-1}$。在量子蒙卡模拟中，需要改变世界线以使体系的细致平衡满足权重因子$W(\{\alpha\})$的位形分布，这就涉及世界线的局域更新（见Fig.2）。

Fig.2 世界线的局域更新。这改变的是动能跳跃的次数，其中(a)只有两个时间切片，对角部分对应一次相互的跳跃，而(b)的相互跳跃被增加的一个时间切片分开。

这些局域移动如果只在固定缠绕数的区域遍历，并不会改变拓扑缠绕数，这是缠绕数的非局域定义。在更高维度，单独世界线的局域移动并不会引起需要包含在内的任意其他置换，涉及两条世界线的局域更新可以用于简化置换，但是并不能改变和周期性边界条件有关的缠绕数。

实验上我们更关注可观测量，所以考虑算符的热力学平均值

$$\begin{align} \langle\mathcal{O}\rangle=\frac{1}{\mathcal{Z}}\mathrm{Tr}\{\mathcal{O}e^{-\beta H} \}. \end{align}\tag{9}$$

把上式写成Eq.4时间切片的形式

$$\begin{align} \langle\mathcal{O}\rangle=\frac{1}{\mathcal{Z}}\sum_{\{\alpha\}}\langle\alpha_0|e^{-\Delta_\tau H}|\rangle\cdots\langle\alpha_2|e^{-\Delta_\tau H}|\alpha_1\rangle\langle\alpha_1|e^{-\Delta_\tau H}\mathcal{O}|\alpha_0\rangle. \end{align}\tag{10}$$

于是联系到之前蒙特卡洛抽样的权重因子

$$\begin{align} \langle\mathcal{O}\rangle=\frac{\sum_{\{\alpha\}}\mathcal{O}(\{\alpha\})W(\{\alpha\})}{\sum_{\{\alpha\}}W(\{\alpha\})}, \end{align}\tag{11}$$

把期待值写成可使用蒙特卡洛抽样的形式。不过上式成立的条件在于分子对应的路径和配分函数对应的路径是否相同。

对于占据数对角存在的量，例如密度关联$\langle n_in_j\rangle$，评估器退化成最平庸的形式，有$\mathcal{O}(\{\alpha\})=\mathcal{O}(\alpha_0)$，不管时间切片算符用何种近似，对角算符这种形式的评估总是有效的。又因为矩阵迹的循环特性，我们可以将算符插入到任意切片位置，从而得到所有时间切片的平均值，并利用

$$\begin{align} \mathcal{O}(\{\alpha\})=\frac{1}{L}\sum_{l=0}^{l-1}\mathcal{O}(\alpha_l) \end{align}\tag{12}$$

增强模拟过程的统计性。而且在不失重要统计信息的情况下，出于节约算力和时间只需要对Eq.12部分求和，例如对每$N$($L>N$)项时间切片求和。主要因为近邻的态其本征值相差不大。

非对角算符的期望值通常是非常复杂的，但是动能算符是个例外。除了容易估计以外，动能算符的期望值也是揭示路径积分物理图像非常有意义的量。对于在配分函数中时间切片算符的线性近似，我们可以利用Eq.10将$Ke^{-\Delta_\tau H}$近似表示。考虑一个位形$\{\alpha\}$使得一个特定的动能算符$K_{ij}$的评估有

$$\begin{align} K_{ij}(\{\alpha\})=\frac{\langle\alpha_1|K_{ij}|\alpha_0\rangle}{\langle\alpha_1|1-\Delta_\tau K|\alpha_0\rangle}. \end{align}\tag{13}$$

可见，仅当第一个时间切片存在格点$i$和$j$的世界线跳跃的时候上式分子才不为零，而此时$K_{ij}(\{\alpha\})=1/\Delta_\tau=L/\beta$。如果我们对Eq.10中所有的$K_{ij}$求平均，则有

$$\begin{align} \langle K_{ij}\rangle=\frac{\langle n_{ij}\rangle}{\beta}. \end{align}\tag{14}$$

所以总的动能的平均值可以用动能跳跃总数的平均值确定，即$\langle K\rangle=-\langle n_K\rangle/\beta$。当然利用经典热力学统计的公式$E=\frac{\partial \ln(\mathcal{Z})}{\partial \beta}$我们仍然可以获得动能的表达式。
Eq.14的意义在于它不仅可以用来计算动能，还携带了关于路径积分和世界线本质的最基本的信息。要问，在一般的世界线位形下可以期待存在多少这样的动能跳跃？在低温下$\beta\rightarrow\infty$，动能几乎不依赖于温度，而与晶格的大小成比例。因此根据Eq.14可以导出动能跳跃数目的预期大小和温标，即有

$$\begin{align} \langle n_K\rangle=-\langle K\rangle\beta\sim N\beta. \end{align}\tag{15}$$

这给出了世界线粗糙度的标度：对于较大的$L$，跳跃的时间间隔约等于$1/N$。这就是为什么我们期待在时间切片算符线性近似的离散路径积分中利用和$N\beta$相同量级的时间切片，因为这样可以避免随着$N$和$\beta$增长带来的误差。而较少的时间切片无法包含所需正确动能行为的事件数目。

Fig.3 连续时间世界线的位形。这里动能事件（格点跳跃）涉及在任意虚时区间$0\geq\tau<\beta$。

Fig.3显示了两种不同的世界线位形，左图满足世界线的时序周期性边界条件，所以贡献的世界线路径是量子配分函数的，而右图则生成湮灭事件（一对生成湮灭算符分别对应实心圆和空心圆。）说明它反映的是$a_i^+a_j$的期望值，这个是Eq.6的一部分，它的傅里叶变换正是系统的动量分布函数$\rho(k)=\langle a_k^+a_k\rangle$。因为这个原因，$a_i^+a_j$这样的非对角算符的期望值就不能利用Eq.11进行求解。如果我们考虑一个构成配分函数的世界线位形，必须满足$i,j$是最近邻表示，而$i,j$如果差了不止一个格位，构造的路径又不满足配分函数的要求，那么对所有$\mathcal{O}(\{\alpha\})$都有$W(\{\alpha\})=0$。

考虑具有势能项的玻色系统，即引入相互作用项，对应Hamiltonian $H=K+V$。类似的，可以将时间切片函数分离

$$\begin{align} e^{-\Delta_\tau H}=e^{-\Delta_\tau K}e^{-\Delta_\tau V}+O(\Delta_\tau^2). \end{align}\tag{16}$$

这里的误差小量是因为$K$和$V$不对易引起的（Glauber公式$e^{\lambda(\hat A+\hat B)}=e^{\lambda \hat A}e^{\lambda \hat B}e^{-\frac{1}{2}\lambda^2[\hat A\,\hat B]}$）。考虑$V$的本征表象，Eq.4时间切片算符的矩阵元则有

$$\begin{align} \left\langle\alpha_{l+1}\left|\mathrm{e}^{-\Delta_\tau H}\right| \alpha_l\right\rangle \approx \mathrm{e}^{-\Delta_{\mathrm{T}} V_l}\left\langle\alpha_{l+1}\left|\mathrm{e}^{-\Delta_\tau K}\right| \alpha_l\right\rangle. \end{align}\tag{17}$$

因此在相互作用下，Eq.8的权重因子变成

$$\begin{align} W(\{\alpha\})=\Delta_\tau^{n_K}\exp\left[-\Delta_\tau\sum_{l=0}^{L-1}V_l\right]. \end{align}\tag{18}$$

蒙特卡洛计算中并不是直接处理配分函数，仅需要位形权重的比值，用以计算世界线更新的接受概率。类似于Fig.2，在世界线上插入和移除动能事件完成一次最简单的序列更新，如果考虑虚时连续极限，这两个事件的概率分别是零和发散的结果。在这种背景下，基于“worm”和循环算法，处理连续极限问题的算法应用而生。下面基于局域更新的图像给出连续极限下蒙特卡洛权重问题的解决思路。

Fig.4 虚时连续极限下的事件更新。随机在世界线位形插入两个时间，用于接受概率的计算。

这里一开始可以设想一段平直连续的世界线线元，在施加$m$个时间切片，事件端点分别为$\tau_1,\tau_2$，即$\tau_2=\tau_1+m\Delta_\tau$。如果相邻格点之间没有世界线，引入两个时间随即改变世界线权重，沿用之前的纯动能Hamiltonian，则权重$W(\{\alpha\})=\Delta_\tau^2$。显然引入事件的插入方式有$m(m-1)/2$种，那么对于这个改变位形的子空间总的相对权重为$\Delta_\tau^2m(m-1)/2$，而原始位形对应的权重则为1。因此在连续极限下，原有位形引入这类事件之后权重因子就变成了$(\tau_2-\tau_1)^2/2$，这一定是个有限值。
除了插入事件仍然要考虑移除事件，Fig.3显示了二者的过程，(a)是插入事件，(b)是移除事件，$\tau_1,\tau_2$分别对应两个事件，当然如果$\tau_1=0,\tau_2=\beta$世界线上没有事件产生。产生的新事件分别命名为$\tau_a,\tau_b$，它们随机分布在$\tau_1$和$\tau_2$的之间的区域。可以针对任意分布在$[\tau_1\,\tau_2]$的已有事件（包括$\tau_{a’}\,\tau_{b’}$）集算总的权重，其相对权重一定为$(\tau_2-\tau_1)^2/2$，同样，如果区间内没有对应的事件权重则为1.
因此可以给出两个事件的接受概率

$$\begin{align} P_{\text {insert }}=\frac{\left(\tau_2-\tau_1\right)^2 / 2}{1+\left(\tau_2-\tau_1\right)^2 / 2}, \quad P_{\text {remove }}=\frac{1}{1+\left(\tau_2-\tau_1\right)^2 / 2}. \end{align}\tag{19}$$

Suzuki–Trotter分解

尽管晶格模型的量子蒙卡算法都是基于连续极限下的路径积分，历史上基于离散时间变量的Suzuki-Trotter近似却最先发展出来用于表示时间切片算符。Eq.16即是Suzuki-Trotter近似的一个例子，不过最一般的关系是

$$\begin{align} e^{\Delta(A+B)}=\left\{\begin{array}{l}e^{\Delta A}e^{\Delta B}+O(\Delta^2)\\ e^{\Delta B/2}e^{\Delta A}e^{\Delta B/2}+O(\Delta^3)\end{array}\right.. \end{align}\tag{20}$$

高阶近似通常过于复杂难以实际应用，不过原则上适当选择足够多的指数函数可以将任意阶的$\Delta$误差约去。如果更关心这些指数乘积的迹，比如在路径积分中，以上低阶近似满足

$$\begin{align} \operatorname{Tr}\left\{\prod_{l=1}^L \mathrm{e}^{\Delta B / 2} \mathrm{e}^{\Delta A} \mathrm{e}^{\Delta B / 2}\right\}=\operatorname{Tr}\left\{\prod_{l=1}^L \mathrm{e}^{\Delta A} \mathrm{e}^{\Delta B}\right\}. \end{align}\tag{21}$$

这里利用了求迹的循环性质。可见即便我们习惯上使用Eq.20第一个展开，对应的系数误差是$\Delta_\tau^2$量级，但实际的误差更小，即$\propto \Delta_\tau^3$。不过由于指数因子的数目是$L=\beta/\Delta_\tau$，所以在固定$N,\beta$的情况下，Eq.20这些指数乘积的总误差仍然是$\propto \Delta_\tau^2$。
另外，可以看出$L$足够大时，对于固定的时间切片误差并不会随着$N$的增大而增长，因此多数情况下可以保证$\Delta_\tau$并不依赖$\beta$和$N$（仍然要留心$\beta\to \infty$误差发散的问题）。

在路径积分的背景下，Suzuki-Trotter近似的实用性在于，即使不能直接估算配分函数的时间切片算符，如果能够找到Hamiltonian合适的分解，例如$H=H_A+H_B$,我们仍然可以利用该近似准确地估算$e^{-\Delta_\tau H_A}$和$e^{-\Delta_\tau H_B}$。一个典型的例子就是仅包含最近邻动能项和相互作用的一维晶格系统

$$\begin{align} H_A=\sum_{\text {odd } i} H_{i, i+1}, \quad H_B=\sum_{\text {even } i} H_{i, i+1}, \end{align}\tag{22}$$

通过这种奇偶分解，$H_A$和$H_B$中所有项都是相互对易的，即对于任意$i,m$有$[H_{i,i+1}\,H_{i+2m,i+2m+1}]=0$。于是在此基础上我们可以将二者分开处理

$$\begin{align} \mathrm{e}^{-\Delta_{\mathrm{r}} H_A}=\mathrm{e}^{-\Delta_{\mathrm{T}} H_{1,2}} \mathrm{e}^{-\Delta_\tau H_{3,4}} \ldots \mathrm{e}^{-\Delta_{\mathrm{T}} H_{N-1, N}} . \end{align}\tag{23}$$

由此类似Eq.4，我们可以将上式插入完备的投影算符，而其中的时间切片算符只涉及两个格点，矩阵元就可以有如下约化

$$\begin{align} \left\langle\beta\left|\mathrm{e}^{-\Delta_{\mathrm{r}} H_{i, j}}\right| \alpha\right\rangle=\left\langle n_{\beta i}, n_{\beta j}\left|\mathrm{e}^{-\Delta_{\mathrm{r}} H_{i j}}\right| n_{\alpha i}, n_{\alpha j}\right\rangle. \end{align}\tag{24}$$

最近邻格点相互作用的系统可以用$2\times2$的矩阵表示（$H_{ij}$也叫所谓的价键算符(The bond operator)，它保证了系统两格点粒子数的守恒），加上玻色态（自旋翻转，用$\alpha,\beta$表示）的两个自由度，所以上面这个矩阵元来自于$4\times4$的块矩阵，四个对角元在占据数表象下有$n_{\alpha i}=n_{\alpha j}=n_{\beta i}=n_{\beta j}$。如果相互作用不局限于最近邻格点，Hamiltonian的分解会更复杂。

直观地看，Eq.24的矩阵元用图的形式可视化对应的是Fig.5的图样。 棋盘连续两行构成一个时间切片，代表近似分解的两支；灰色格子代表非零的矩阵元，代表一个动能事件。和之前线性近似的模型不同，单独一个时间切片可以包含更多Hopping而不是只能容纳一个。

Fig.5 Suzuki-Trotter近似下一维世界线位形的棋盘分解图样。

Fig.5右边显示的六个模式对应下式

$$\begin{equation} \begin{aligned} &\left\langle\uparrow_i \uparrow_j\left|\mathrm{e}^{-\Delta_\tau H_{i j}}\right| \uparrow_i \uparrow_j\right\rangle=\left\langle\downarrow_i \downarrow_j\left|\mathrm{e}^{-\Delta_{\mathrm{r}} H_{i j}}\right| \downarrow_i \downarrow_j\right\rangle=+\mathrm{e}^{-\Delta_\tau / 4} \\ &\left\langle\uparrow_i \downarrow_j\left|\mathrm{e}^{-\Delta_\tau H_{i j}}\right| \uparrow_{i \downarrow j}\right\rangle=\left\langle\downarrow_i \uparrow_j\left|\mathrm{e}^{-\Delta_{\mathrm{r}} H_{i j}}\right| \downarrow_i \uparrow_j\right\rangle=+\mathrm{e}^{\Delta_\tau / 4} \cosh \left(\Delta_\tau / 2\right) \\ &\left\langle\downarrow_i \uparrow_j\left|\mathrm{e}^{-\Delta_\tau H_{i j} j}\right| \uparrow_i \downarrow_j\right\rangle=\left\langle\uparrow_i \downarrow_j\left|\mathrm{e}^{-\Delta_\tau H_{i j}}\right| \downarrow_i \uparrow_j\right\rangle=-\mathrm{e}^{\Delta_\tau / 4} \sinh \left(\Delta_\tau / 2\right). \end{aligned} \end{equation} \tag{25}$$

这里正是对海森堡反铁磁模型的描述，将之前讨论的玻色态$n_i$换成了自旋态$\uparrow,\downarrow$，得到海森堡交换关联项$H_{i,i+1}=\boldsymbol{S}_i\cdot\boldsymbol{S}_{i+1}$。值得注意的是，正如上式最后一项非对角矩阵元的符号，它是个负号，而在自旋阻挫系统由于同时具有正的权重和负的权重，就存在所谓“符号问题”。所以实际应用中世界线方法和类似的量子蒙卡方法主要用于研究双线性自旋系统和玻色模型。对于费米系统，世界线的置换仍然会引起符号问题，仅在一维系统且只考虑全局的循环置换才可能避免这个问题（通过选择对称或反对称周期边界条件）。另外，在更高维度，Suzuki-Trotter分解仍然可以适用局域更新的方式实现量子蒙卡的抽样。

而即便Suzuki-Trotter方法中时间切片的数目可以不必依赖系统的大小（因为不限制单位切片的事件数），量子蒙卡的计算量仍然与连续极限的公式一致，也就是用$N\beta$的占据数来标志（这个结果可以在Eq.15通过储存和操纵时间切片的事件获得）。出于这个原因，连续时间方法通常更为高效，而且不受离散化误差的影响。不过离散化方法在实际应用中仍然不可或缺。在一些考虑耗散的有效模型中，由于世界线存在含时的相互作用，很难取连续极限，离散化方法就不失为一种可行的方法。

级数展开表示

除了离散化时间切片方法或者连续极限的路径积分方法，仍然可以利用对玻尔兹曼算符泰勒展开构造可供蒙特卡洛抽样的位形空间。展开如下

$$\begin{align} \mathrm{e}^{-\beta H}=\sum_{n=0}^{\infty} \frac{(-\beta)^n}{n !} H^n. \end{align}\tag{26}$$

在特定态下，配分函数可以写作

$$\begin{align} \mathcal{Z}=\sum_{n=0}^{\infty} \frac{(-\beta)^n}{n !} \sum_{\{\alpha\}_n}\left\langle\alpha_0|H| \alpha_{n-1}\right\rangle \cdots\left\langle\alpha_2|H| \alpha_1\right\rangle\left\langle\alpha_1|H| \alpha_0\right\rangle. \end{align}\tag{27}$$

下标$\{\alpha\}_n$表明这是对$n$个态的求和，同Eq.5线性离散化的路径积分公式可比照来看。如果Hamiltonian采用Eq.6的玻色动能项，我们可以绘制出和Fig.1相似的世界线图样。二者主要的不同在于世界线切片的数目，即展开指数$n$是变动的，对于给定的$n$，每一个切片都存在粒子的hopping（仔细观察两个公式，线性离散化公式插入的完备基是针对时间切片算符的，所以沿着这个思路自然要考虑时间切片的连续极限。而级数展开公式插入的完备基是关于$n$的求和项，不需要人为考虑时间切片的大小，这是二者本质的区别。）。在这个表示下新的传播子维度和虚时演化下的不同，但是下标为$p=0,\cdots,n-1$的态和路径积分下的虚时密切相关。显然在只考虑动能项的情况下，Eq.27所有的矩阵元在任意位形下都等于-1，所以可以把展开式的负号抵消。所以蒙特卡洛的权重因子可以写作

$$\begin{align} W(\{\alpha\}_n)=\frac{\beta^n}{n!}. \end{align}\tag{28}$$

因为考虑路径积分量子化，玻色系统和非对角项没有阻挫的量子自旋系统下的权重因子都是正定的，因此可以利用正定条件实现一些技巧。首先导出系统内能

$$\begin{align} E=\frac{1}{\mathcal{Z}} \sum_{n=0}^{\infty} \frac{(-\beta)^n}{n !} \sum_{\{\alpha\}_{n+1}}\left\langle\alpha_0|H| \alpha_n\right\rangle \cdots\left\langle\alpha_2|H| \alpha_1\right\rangle\left\langle\alpha_1|H| \alpha_0\right\rangle. \end{align}\tag{29}$$

由于$H$额外增加了一个矩阵元，对给定$n$在Eq.27中实际上是对$n+1$个态的求和。上式所有位形都贡献于配分函数，但是权重因子则不同，所以为了匹配二者的位形可以把上式写成

$$\begin{align} E=-\frac{1}{\mathcal{Z}} \sum_{n=1}^{\infty} \frac{(-\beta)^n}{n !} \frac{n}{\beta} \sum_{\{\alpha\}_n}\left\langle\alpha_0|H| \alpha_n\right\rangle \cdots\left\langle\alpha_2|H| \alpha_1\right\rangle\left\langle\alpha_1|H| \alpha_0\right\rangle . \end{align}\tag{30}$$

如此扩展对$n$的求和不包含$n=0$（注意这里多了一个负号），使得Eq.27和Eq.30准确匹配。于是可以用$n/\beta$作为对系统能量的估算公式，在配分函数中依据对应权重的位形中抽样，系统能量有

$$\begin{align} E=- \frac{\langle n\rangle}{\beta}. \end{align}\tag{31}$$

当把$H$写成求和形式，即$H=-\sum_iH_i$，可以导出单项$H_i$的期望值

$$\begin{align} \langle H_i\rangle=\frac{\langle n_i\rangle}{\beta}. \end{align}\tag{32}$$

也就是出现在配分函数展开式中算符字符串出现的平均数（同样需要注意上下两式的符号。）。这个结果和Eq.14路径积分的表达式是等价的。所以这显示了两个公式中路径之间一一对应的关系，这种关系只在离散或连续路径积分并不涉及事件发生具体的时间点（切片时间点$\tau_l$），而只考虑动能事件在路径中发生顺序的时候有所体现。对于级数展开方法，主要不同在于系统的势能项和动能项处理的方式是一致的，这个在路径积分方法中则是利用Suzuki-Trotter近似分开处理。对于势能项所对应的“非事件”性质的影响，这本身并不改变世界线的图样，但是会影响修正路径权重的对角矩阵元的结果。而在路径积分方法中，则是以势能因子的形式在蒙特卡洛抽样中体现（Eq.11）。

从Eq.30中观察到$n=\infty$，但是在实际应用中仍然要考虑在$n_{\mathrm{max}}\propto N\beta$截断操作。这样做是非常合乎计算上的要求，满足物理上的需要。考虑对系统比热的评估，对Eq.31的内能求温度的导数

$$\begin{align} C=\langle n^2\rangle-\langle n\rangle^2-\langle n\rangle. \end{align}\tag{33}$$

当$T\to 0$，比热应该为零，即上式关于$n$分布的方差应该等于$\langle n\rangle$，所以很明显当该分布超过某个正比于$N\beta$的$n$次方时，它会以指数形式衰减。事实上，随机级数展开算法会基于该表示依据不同幂次$n$区域的相对权重自动对$n$的分布抽样。

通过对泰勒级数的截断，可以进一步地明晰路径积分方法和随机级数展开方法之间的关系。当截断发生在$n=L$，利用$L-n$个单位算符$I$在$n<L$的情况下通过增加所有幂次$H^n$，我们可以正式地构造一个固定大小的抽样空间。允许在$L$个算符乘积中插入所有可能的单位算符，我们能定义算符串$S=S_1,S_2,\cdots,S_L$，满足$S_i\in\{H,I\}$，并对所有序列求和。然后必须用等价项的数字$\begin{pmatrix}
L\\n
\end{pmatrix}$补偿权重，于是配分函数修正为

$$\begin{align} Z=\sum_S \frac{(-\beta)^n(L-n) !}{L !} \sum_{\{\alpha\}_L} \sum_{\left\{S_i\right\}}\left\langle\alpha_0\left|S_L\right| \alpha_{L-1}\right\rangle \cdots\left\langle\alpha_2\left|S_2\right| \alpha_1\right\rangle\left\langle\alpha_1\left|S_1\right| \alpha_0\right\rangle. \end{align}\tag{34}$$

这里$n$是指在算符串中元素$S_i=H$的数目。如果我们考虑之前的玻色系统，取连续极限$L\to\infty$，权重将约化为$(\beta/L)^n$，这和Eq.6在$L$个时间切片下的结果一样。在这个极限下，指标$p=0,\cdots,L-1$和虚时的关系有$\tau=p\beta/L$。而两者主要的区别在于全部权重$\beta^n(L-n)!/L!$的级数展开的结果对于$L\propto N\beta$是准确的，而在路径积分中为了避免离散化误差我们必须采用连续极限。对于有限的$L$,级数展开的指标$p$并不严格对应于虚时，但是表示虚时的一个分布。不过鉴于传播子指标$p$和虚时$\tau$之间的紧密关系，将$p$空间视作级数展开中时间维度也是合适的。
总的来看，随机级数方法以其易用性准确性相比于虚时极限的离散化表示更为有效，也仅仅在势能项比较大的系统中后者才体现出高效性（通常是因为级数展开算符$\langle n\rangle$远大于路径积分中动能事件数 \langle n_K\rangle$），所以对于量子自旋系统这种势能项和动能项有良好平衡的系统，基于随机级数的量子蒙卡方法有更多的优势。

海森堡模型的随机级数展开方法

位形空间

蒙特卡洛抽样

可观测量的评估

随机级数展开方法在1D和2D系统的应用

海森堡链

梯子系统

二维的长程序

二聚系统的量子相变

J-Q模型奈尔-价键固体相变

参考

1. Sandvik A W. Computational studies of quantum spin systems[C]//AIP Conference Proceedings. American Institute of Physics, 2010, 1297(1): 135-338.
2. N. Hatano and M. Suzuki, in Quantum Annealing and Other Optimization Methods, Edited by B. K. Chakrabarti and A. Das (Springer-Verlag, 2005)