自旋-玻色子模型的量子蒙特卡罗模拟：虫洞更新方法

使用虫洞更新的自旋-玻色子模型的量子蒙特卡洛模拟#

Manuel Weber

马克斯·普朗克复杂系统物理研究所，Nöthnitzer Str. 38, 01187 德累斯顿，德国

0（2021年8月30日收到；2022年3月31日修订；2022年4月5日接收；2022年4月15日发表）

我们提出了一种用于耦合到耗散玻色子浴的自旋系统的精确量子蒙特卡洛方法，该方法利用非局域虫洞更新来模拟由非对角自旋-玻色子耦合引起的延迟自旋翻转相互作用。该方法与随机级数展开密切相关，并将全局有向环更新的范围扩展到通过世界线构型的非局域移动。我们针对U(1)对称的双浴自旋-玻色子模型测试了该方法，其中自旋1/2粒子的非对角分量耦合到具有幂律谱的相同独立浴，并得到了临界相与局域相之间临界耦合的精确估计。我们的方法适用于任意空间维度中耦合到具有任意谱分布的玻色子模式的杂质系统和晶格模型。

I. 引言#

蒙特卡洛方法全局更新方案[1–4]的发展使得对出现在广泛强关联物理中的临界现象进行高精度数值研究成为可能。对于量子系统，虫洞算法和有向环算法[5–7]已成为分别模拟玻色子或自旋模型的先进技术。虽然这些方法随机采样配分函数的微扰展开，但虫洞/有向环更新沿着在扩展构型空间中构建的闭合环改变世界线构型。构建这些环的计算代价仅与系统大小和逆温度成线性关系，因此只要不存在符号问题，就能够在任何空间维度中进行高效模拟。因此，非常期望将这些算法扩展到无法以其他方式高效模拟的新模型类别。

计算多粒子理论中的主要挑战是模拟耦合到不遵守粒子数守恒的局域玻色子模式的量子系统。研究最多的例子包括具有局域耦合到声子的晶格模型。基于矩阵乘积态(MPS)的方法需要大键维来处理每个模式的无界玻色子希尔伯特空间，这使得它们比纯粹电子或自旋模型效率更低；然而，优化的基组可以带来显著改进[8]。当玻色子以第一量子化形式采样时，量子蒙特卡洛(QMC)模拟通常受到长自相关时间的困扰[9]。到目前为止，其二量子化形式中缺乏粒子数守恒阻碍了虫洞/有向环更新的简单高效表述（例外情况见文献[10,11]），使得低效的局域更新成为模拟玻色子的唯一可用选项[12,13]。如果哈密顿量在玻色子场中是二次型的，则可以避免直接玻色子采样的困难。此时，可以使用路径积分精确积分掉玻色子，从而在系统的自由度中获得延迟相互作用。最近将定向环算法推广到延迟相互作用[14]已被证明能够克服自相关问题，并已成功应用于一维[15,16]和二维[17]的声子耦合系统。

与玻色子模式的耦合在量子物理的许多领域都扮演着重要角色。例如，在固态系统中，玻色子激发以声子或磁振子等准粒子的形式出现；离子阱量子模拟器依赖于由离子振动模式中介的长程相互作用[18,19]；而在量子光学中，原子则与量子化的电磁场相互作用[20,21]。一个核心问题是当量子系统连接到具有无穷自由度的玻色子库时的耗散现象[22]。其最简单的实现之一是自旋-玻色子模型：一个二能级系统耦合到一个连续的玻色子模式。在该模型中，由外加磁场引起的库致局域化与退局域化之间的竞争会导致量子相变，而其临界性质直到发展出一种高效的团簇量子蒙特卡洛(QMC)方法，通过延迟相互作用模拟与库的耦合后才得以阐明[23]。一种变分矩阵乘积态(MPS)方法揭示了，将两个非对易自旋分量耦合到独立库可能会产生更复杂的相图[24]，其中包含一个由退相干机制中的受挫效应[27]引发的临界相[25,26]；然而，MPS方法向多库情形的推广受到局域玻色子希尔伯特空间大维度的限制。带有在位耗散的晶格模型已使用经典蒙特卡洛方法模拟长程相互作用[28,29]；完整的QMC模拟仅通过蠕虫算法应用于对角自旋-玻色子耦合[30]。到目前为止，对于更一般的自旋-玻色子相互作用，我们仍然缺乏能够同等高效地处理杂质模型和晶格模型的QMC更新方案。

量子杂质模型是动力学平均场理论的核心，其中到无穷大库的耦合模拟了与晶格上相邻粒子的相互作用效应[31]。对高效杂质求解器的需求推动了费米子系统连续时间QMC方法的发展，其中库的耦合以延迟相互作用的形式处理[32-34]。这些方法中最著名的一种基于杂化展开，并对费米子世界线构型进行采样[33]。该方法通过将费米子库传播子合并为一个行列式来避免负符号问题，该行列式仅允许局域蒙特卡洛更新，并最终导致反比于温度的三次方标度。在此，蠕虫采样仅用于获得格林函数的改进估计量[35]。该方法已被扩展以包含玻色子库[36]，用于求解玻色子杂质问题[37]、自旋-玻色子模型[38]或玻色-费米近藤模型[38-41]，但即使在没有纯自旋-玻色子模型中费米子行列式的情况下[38,41]，蒙特卡洛采样仍然局限于局域更新。例如，自旋-玻色子杂质模型的自洽求解对自旋玻璃[42-45]以及动力学平均场理论向自旋系统的扩展[46]都起着重要作用。玻色子模式也已被纳入相互作用展开的QMC方法中，该方法可应用于有限晶格上的费米子-玻色子模型[47-49]。

在本文中，我们提出了一种用于耗散自旋模型的连续时间量子蒙特卡洛方法。该方法结合了杂质求解器的优势以及为格点模型开发的全局虫洞/有向环更新。我们的方法基于具有全局有向环更新的随机级数展开[50]及其最近对延迟相互作用的推广[14]。与声子耦合系统不同，下面考虑的自旋-玻色子模型不守恒总$\hat{S}_{z}$量子数，这阻碍了虫洞/有向环更新的局域构造。为克服这一困难，我们引入了新的虫洞更新：虚时非局域相互作用中的两个算符可以作为虫洞/有向环的源和汇，从而允许通过世界线构型进行非局域移动。结果表明，构建有向环的规则与具有最近邻相互作用的常规自旋模型等价——只需将格点位指标替换为虚时变量——这使得可以轻松移植参考文献[7]中开发的方法。延迟相互作用的时间依赖性不影响环的构建，因为它在对角更新期间被采样[14]；对于任何谱分布的浴模式，浴传播子都可以高效采样。在下面的内容中，我们为杂质模型阐述了我们的量子蒙特卡洛方法，但虫洞更新可以平凡地扩展到有限格点上的自旋-玻色子模型，我们将在本文末尾讨论这一点。我们通过双浴自旋-玻色子模型展示了我们量子蒙特卡洛方法的效率。我们的方法达到了比先前量子蒙特卡洛方法[38,41]显著更低的温度，这使我们能够精确估计临界相和局域相之间的量子临界耦合，该估计与先前基于矩阵乘积态方法[24]的结果非常一致。未来，我们的量子蒙特卡洛方法将能够高效模拟各种与玻色子模式（如光子、声子、磁振子等）耦合的量子杂质和格点模型，这些模型出现于量子光学、离子阱模拟器或固态系统中。

本文组织如下。在第II节中，我们定义了自旋-玻色子模型；在第III节中，我们展示了如何从相互作用绘景推导出延迟相互作用；在第IV节中，我们介绍了我们的量子蒙特卡洛方法；在第V节中，我们展示了我们的结果；在第VI节中，我们讨论了方法的扩展；在第VII节中，我们进行了总结。

II. 自旋-玻色子模型#

我们考虑一个通用形式的自旋-玻色子模型

它可以分为自旋部分$\hat{H}_{\mathrm{s}}$、玻色子部分$\hat{H}_{\mathrm{b}}$和自旋-玻色子相互作用$\hat{H}_{\mathrm{sb}}$。目前，$\hat{H}_{\mathrm{s}}$描述一个局域自旋-$1 / 2$自由度$\hat{S}_{\ell} = \textstyle{\frac{1} {\gamma}} \hat{\sigma}_{\ell}$，其中$\hat{\sigma}_{\ell} , \ell \in \{x , y , z \}$是泡利矩阵（我们设$\bar{h} \overset{=} {=} 1$）。玻色子浴由谐振子的和给出，

其中$\hat{a}_{\mu}^{\dagger} ~ ( \hat{a}_{\mu} )$在频率为$\omega_{\mu}$的态$\mu$中产生（湮灭）一个玻色子。下面，上标$\mu$将指代一个连续的振子集合和浴模式的不同分量，但它也可以标记有限系统中的格点位置。自旋-玻色子相互作用

将浴耦合到一个自旋算符$\hat{\varrho}_{\mu} [ \hat{S}_{\alpha} ]$，该算符依赖于模型，并且也包含自旋-玻色子耦合常数。下面，我们定义两种符合通用形式类型的自旋-玻色子模型。

原始的泡利自旋-玻色子模型由一个两能级系统组成，该系统通过自旋的z分量耦合到玻色子浴，即

磁场 $h_{x}$ 在两个能级之间诱导了有限的隧穿，而通过具有模相关耦合常数 $\gamma_{q}$ 的连续浴的耦合会使自旋定域化。自旋-玻色子模型已被推广到与最多三个浴相互作用的情形，即

其中每个浴耦合到不同的自旋分量。我们设 $\gamma_{qx} = \gamma_{qy}$，并由于该模型与XXZ模型的相似性（这在后文引入延迟相互作用时会变得明显），将其称为XXZ自旋-玻色子模型。当 $\gamma_{qz} = 0$ 时，我们恢复为双浴自旋-玻色子模型；当 $\gamma_{qx} = \gamma_{qy} = \gamma_{qz}$ 时，则为玻色-近藤模型。在无磁场的情况下，前者具有U(1)对称性，而后者具有SU(2)对称性。关于U(1)对称情形的更多细节将在第V节讨论。关于杂质性质的综述可在文献[51]中找到。

量子光学中另一个经常研究的哈密顿量是Jaynes-Cummings模型[52]

其中量子化电磁场通过自旋翻转算符 $\hat{S}_{\pm} = \hat{S}_{x} \pm \mathrm{i} \hat{S}_{y}$ 与二能级系统相互作用。

我们考虑与连续浴模的耦合。杂质的物理性质完全由每个浴的谱函数决定，

我们假设 $J_{\ell} ( \omega )$ 具有指数为 $s$ 的幂律形式，

其中 $s = 1$ 对应欧姆浴，$0 < s < 1$ 对应亚欧姆区域。这里，$\alpha_{\ell}$ 是自旋-浴耦合常数，我们将在后文中使用，并选取频率截断 $\omega_{\mathrm{c}} = 1$ 作为能量单位。

III. 相互作用图像与延迟相互作用#

上一节引入的自旋-玻色子模型在玻色子算符中是二次的，因此可以精确积分掉浴。通常，我们使用路径积分形式推导系统自由度中的延迟相互作用。为避免自旋路径积分的复杂性，我们将展示延迟相互作用如何从相互作用图像中的微扰展开中产生。通过这种方式，我们将看到先前直接模拟玻色子的QMC方法所面临的困难是如何被克服的。我们的形式与用于费米子杂质模型的杂交展开[33,34]非常相似。

我们将整个系统的哈密顿量 $\hat{H} = \hat{H}_{0} + \hat{V}$ 分为无扰部分 $\hat{H}_{0}$ 和微扰部分 $\hat{V}$。配分函数的戴森展开为

其中所有算子 $\hat{V} ( \tau ) = e^{\tau \hat{H_{0}}} \hat{V} e^{- \tau \hat{H_{0}}}$ 按其虚时变量排序，即 $\beta \geqslant \tau_{1} \geqslant \tau_{2} \geqslant \dots \geqslant$ $\tau_{m} \geqslant 0 ( \beta = 1 / T$ 为逆温度；我们设 $k_{\mathrm{B}} = 1 )$ 为方便推导，我们引入时序算子 $\hat{\mathcal{T}}_{\tau}$，并将式 (9) 改写为

对于式 (1) 中的通用自旋-玻色子模型，我们确定

类似于费米子杂质模型中使用的杂化展开 [34]。在推导延迟自旋相互作用时，我们设 $\hat{H}_{\mathrm{s}} = 0$ 以简化符号表达。稍后，我们可以在不丧失一般性的情况下将 $\hat{H_{\mathrm{s}}}$ 包含在 $\hat{V}$ 中。为区分式 (11) 中的正则算子和伴随算子，我们引入了上标 c 及其相反的上标 ¯c。由此，相互作用展开变为

迹分解为自旋和玻色子各自的独立贡献。由于玻色子粒子数守恒，只有当 $m = 2 n$ 且产生与湮灭算子数目相等时，$\langle \bullet \rangle_{\mathrm{b}} = Z_{\mathrm{b}}^{- 1} \mathrm{Tr}_{\mathrm{b}} [ e^{- \beta \hat{H}_{\mathrm{b}}} \bullet ]$ 才给出非零贡献。玻色子迹可以通过 Wick 定理进一步简化。例如，考虑如下顺序排列算子的期望值：

这里，$S_{n}$ 是 n 阶对称群，$\pi \in S_{n}$ 是 n 个对象的排列。我们引入了自由玻色子传播子 $D ( \omega_{\mu} , \tau - \tau^{\prime} ) = \langle \hat{\mathcal{T}}_{\tau} \hat{a}_{\mu} ( \tau ) \hat{a}_{\mu}^{\dagger} ( \tau^{\prime} ) \rangle_{\mathrm{b}}$，其表达式为

并满足 $D ( \omega , \tau + \beta ) = D ( \omega , \tau )$。式(13)左边仅表示选择 $\{c_{1} , \ldots , c_{n} \}$ 的一种可能组合，使得产生算符和湮灭算符数量相等。总共有 $\binom{2 n} {n}$ 种组合。将式(13)代入式(12)，我们得到（定义 $\tau_{k}^{\prime} = \tau_{n + k} )$）

对所有排列的求和可以通过重新标记变量来求值，这在式(15)中引入了另一个因子 n!。最终，微扰展开可以改写为

对式(16)中得到的相互作用项进行重新指数化后，配分函数变为

因此，自旋-玻色耦合产生了一个延迟的自旋-自旋相互作用

该相互作用由自由玻色传播子(14)介导。注意，对于我们所选的 $\dot{H}_{0}$，自旋的时间演化是平凡的，因为 $\hat{H}_{0}$ 不包含任何自旋算符。然而，我们仍然需要 τ 标签来对微扰展开中的自旋进行时间排序，同时也用于玻色传播子的时间依赖关系。

最后，我们指定上一节引入的自旋-玻色模型的延迟自旋相互作用。对于 XXZ 自旋-玻色模型，我们得到

而 Jaynes-Cummings 哈密顿量则导致

这里我们已经对玻色子热库取了连续极限。我们明确看到，式(7)中定义的谱函数完全决定了热库对自旋子系统的影响。对于式(8)中的幂律谱，$D ( \omega , \tau )$的频率平均在虚时中诱导出一个长程相互作用，当$\omega_{\mathrm{c}} \tau \gg 1$时，该相互作用以$1 / \tau^{s \bar{+} 1}$形式衰减。由于XXZ自旋-玻色子模型的延迟相互作用在$\tau \tau^{\prime}$下是对称的，我们在式(19)中将$D ( \omega , \tau )$替换为对称化传播子$D_{+} ( \omega , \tau ) = [ D ( \omega , \tau ) + D ( \omega , \beta - \tau ) ] / 2$。

IV. 量子蒙特卡罗方法#

A. 蒙特卡罗方法的基本思想#

在描述我们的QMC方法的数学形式之前，我们首先想要介绍支撑该方法的主要思想。我们采用定向环算法[7]的框架，该算法最初是在SSE表示[50]中提出的，最近被推广到延迟相互作用[14]。在上一节中，我们已经证明，式(1)形式的一般自旋-玻色子模型可以映射到具有延迟相互作用的自旋问题。这样，我们就可以避免直接玻色子采样带来的困难，这种困难通常出现在粒子数不守恒的情况下。由于式(18)中得到的延迟自旋相互作用也不守恒总$\hat{S}_{z}$，我们将引入非局域虫洞更新，以利用定向环算法提供高效采样。

作为我们方法公式化的出发点，考虑式(17)中的配分函数，其中玻色子热库已被积分掉，从而得到延迟自旋相互作用。遵循文献[14]，我们在完整相互作用$\hat{\mathcal{H}}_{\mathrm{ret}}^{\mathrm{~ ~}}$中展开时间有序指数，即围绕平凡自旋部分$\hat{\mathcal{H}}_{0} = 0$进行相互作用展开，得到式(16)。这种在全哈密顿量中的展开是导致SSE表示的特征性特征，即使我们从相互作用表象出发也是如此。然后，自旋算符$\hat{S}_{\ell} ( \tau )$的时间演化变得平凡，其τ标签仅用于执行时间排序。对于等时（即时间无关）哈密顿量，所得展开中的虚时积分可以精确计算，从而精确映射到不具有显式时间依赖性的SSE表示[10,12]。在存在延迟相互作用的情况下（这是本文的重点），我们的蒙特卡罗方案必须显式地对虚时变量进行采样，以正确考虑玻色子传播子的时间依赖性。文献[14]表明，连续时间的精确采样仅在对角更新中引入（如下文所述），因此定向环算法的其余部分可以像在SSE表示中那样进行公式化。

我们首先讨论自旋算符的迹迹如何与原始SSE公式相关联。自旋算符的时序乘积对应于SSE表示中的算符序列。根据我们在第III节中的推导，时序不会导致任何负号是清晰的。自旋的迹 $\begin{array} {r} {\mathrm{Tr}_{\mathrm{s}} \bullet = \sum_{\alpha} \left. \alpha \right| \bullet \left| \alpha \right.} \end{array}$ 遍历初始态 $| \alpha \rangle$，该初始态通常选择在 $\hat{S}_{z}$ 基下。我们假设 $\hat{\mathcal{H}}_{\mathrm{ret}}$ 中的算符是非分支的，即 $\hat{\varrho}_{\mu_{p}}^{\bar{c}_{p}} ( \tau_{p} ) | \alpha_{p} \rangle \sim | \alpha_{p - 1} \rangle$ 唯一地将一个基态传播到另一个基态；此处我们按式(9)定义时序（传播指数递减），并设 $\left. \alpha_{0} \right. =$ $| \alpha_{2 n} \rangle = | \alpha \rangle$ 。初始态的虚时演化可以用世界线图像可视化（注意，世界线由传播的自旋态的时间依赖性定义，而非由算符定义）。由于传播态 $\textstyle | \alpha_{p} \rangle$ 仅在施加算符时发生变化，因此找到相互作用顶点的图形表示就足够了。$\hat{\mathcal{H}}_{\mathrm{ret}}$ 可能的顶点类型如图1所示。每个顶点由两个算符、玻色子传播子以及施加算符前后的自旋态组成，分别用两个黑色竖条、一条虚线和四条腿表示。对于我们的自旋-1/2模型，腿以实心（空心）圆表示，对应于 $\hat{S}_{z}$ 本征值为 $1 / 2 \left( - 1 / 2 \right)$ 的本征态。总的来说，对于式(19)和式(20)中的延迟相互作用，共有六种顶点类型：四个对角顶点（保持世界线构型不变）和两个非对角顶点（其中自旋态被翻转）。后者对应于 $\hat{S}_{-} ( \tau ) \hat{S}_{+} ( \tau^{\prime} )$ 和 $\hat{S}_{+} ( \tau ) \hat{S}_{-} ( \tau^{\prime} )$ 。包含两种顶点类型的典型世界线构型如图2(a)所示。关于相互作用顶点以及我们QMC方法构型空间的更正式定义，请参见第IV.B节。

图1：延迟相互作用顶点的图形表示。时间 τ 和 τ' 处的黑色竖条对应于一个顶点的两个算符，而实心（空心）圆表示施加算符前后、具有 $\hat{S}_{z}$ 本征值 $1 / 2 \ ( - 1 / 2 )$ 的自旋态。两个算符由一条虚线连接，代表自由玻色子传播子。对角顶点由顶点类型 $v = 1 \ – 4$ 标识，而非对角项 $\hat{S}_{-} ( \tau ) \hat{S}_{+} ( \tau^{\prime} )$ 和 $\hat{S}_{+} ( \tau ) \hat{S}_{-} ( \tau^{\prime} )$ 分别标记为 $v = 5$ 和 $v = 6$ 。

我们使用马尔可夫链蒙特卡罗采样来随机地评估对所有世界线构型的求和。为了构建一个遍历算法，我们需要两种类型的更新：对角更新和定向环更新。

图2：(a) 包含两个延迟顶点的世界线配置示意图。初始态|α>的虚时传播方向从左至右。由于仅当施加算符时传播态才会改变，因此只需标明每个顶点支路上的自旋配置即可。(a)与(b)之间对角顶点的添加或移除不会改变世界线配置。(c) 为初始化定向环更新，我们随机选取初始时刻并将环头传播至顶点入口支路。(d) 根据定向环方程确定的概率选择出口支路，此处环头利用自由玻色子传播子建立的非局域连接实现虫洞移动。(e) 重复此过程直至环头返回起始点。(f) 最终更新世界线配置，即翻转环连接的所有支路（本例中还包括初始态）。同时对角顶点可转化为非对角顶点，反之亦然。

对角更新的主要目的是改变平均展开阶数、用新顶点替换现有顶点，并——在原始SSE方法基础上——对虚时变量进行采样。由于对角算符不改变世界线配置，我们建议采用第IV.C节定义的梅特罗波利斯方案添加或移除对角顶点。图2(a)与2(b)分别展示了添加与移除此类顶点的示例。通过根据自由玻色子传播子对顶点算符间的时间差进行采样，可显著提升新顶点添加的接受率。每次对角更新块需要${\mathcal{O}} ( n )$次提议操作，以便用新顶点替换大量旧的顶点变量。

我们的蒙特卡洛方法的效率依赖于定向环更新[7]，这是一种全局更新方案，它沿着一个闭合环改变大段世界线的子集。该环的构造如下：首先，我们随机选取一个起始时刻，插入一对自旋翻转算符。其中一个算符被标识为环的头部，它可以作为世界线的不连续性在虚时中移动，而另一个算符则成为环的尾部，固定在起始时刻。由于自旋算符的时间演化是平凡的（$ \hat{H}_{0}$ 不含任何自旋项），头部将一直前进，直到它到达某个顶点中算符的腿部。在传统的定向环更新中，我们会根据局部细致平衡的定向环方程所确定的概率，选择一个局部顶点的出口腿。事实证明，这对于延迟相互作用的非局部顶点仍然成立。因此，在图2(c)中，通过左算符的左腿进入顶点的环头，可以选择反弹并反转传播方向、径直穿过左算符，或者从右算符退出。在图2(d)中，我们选择了最后一个选项，称之为虫洞更新，因为环头利用自由玻色子传播子建立的连接作为虫洞，在虚时中相距甚远的点之间瞬间隧穿——这一命名也类比于虫算法，其中环被称为虫[5]。一旦环从某个顶点的腿部退出，它将传播到下一个顶点的入口腿部，并重复相同的过程，直到环头返回到其尾部，如图2(e)所示。当环闭合时，我们翻转沿环的所有自旋构型，如图2(f)所示。因此，对角顶点可以转变为非对角顶点，反之亦然。由于定向环能够触及大量的顶点，它代表一种全局更新操作。

构建回路的计算开销与所触及顶点数量成线性关系。关于有向回路构造的更多细节见第IV.D节。

通过上述对角更新和有向环更新，我们的 QMC 算法采样了对角和非对角算符的完整图解展开，从而满足了遍历性。对角更新通过在任何虚时位置添加和移除对角顶点来采样展开阶数 $n$，而有向环更新则将对角算符转换为非对角算符，反之亦然。后一种转换能够达到所有具有非零权重的顶点类型，如其原始表述[7]所述。

总之，我们将用于杂质求解器以避免直接玻色子采样困难的延迟相互作用框架[34]与最初为格点模型开发的全局有向环更新[7]相结合。对于具有正定传播子的玻色子浴，虫洞更新的局部构造是可行的，而费米子浴则需要数值上昂贵的行列式计算以避免符号问题[33]。我们算法的表述与有向环算法的标准表述[7]密切相关。然而，虚时中的非局域性要求在实现上进行一些更改，我们将在附录中讨论。尽管我们是在有向环算法的框架中开发了虫洞更新，但其核心思想可以推广到基于虫/有向环更新的相关 QMC 方法。

B. 构型空间与顶点权重#

一个蒙特卡罗配置 $\mathcal{C} = \{n , \mathcal{C}_{n} , | \alpha \rangle \}$ 由展开阶数 $n$、有序顶点列表 $\mathcal{C}_{n} =$ $\{\nu_{1} , \ldots , \nu_{n} \}$ 以及本地 $\hat{S}_{z}$ 基下的初始态 $|\alpha\rangle$ 定义。每个顶点变量 $\nu = \{t_{\mathrm{int}} , a , \omega , \tau , \tau^{\prime} \}$ 包含一个频率 $\omega$ 和两个虚时值 $\tau$ 与 $\tau^{\prime}$。算符类型 $a$ 区分顶点的对角部分和非对角部分，而相互作用类型 $t_{\mathrm{int}}$ 标记不同的相互作用种类——后者仅在我们结合自旋-玻色子模型与额外相互作用时相关，如第 VI 节所述。我们自旋-玻色子模型的相互作用顶点由下式给出：

这里，我们对自由玻色子传播子和谱密度进行了重新缩放，

使得前者成为关于 $\tau$ 的概率分布，后者成为关于 $\omega$ 的概率分布。对于式 (8) 中的幂律谱，我们得到 $\mathcal{I} ( \omega ) = s \omega_{\mathrm{c}}^{- s} \omega^{s - 1}$。谱密度的归一化因子可以被吸收进算符部分 $\hat{h}_{a} ( \tau , \tau^{\prime} )$。

我们的蒙特卡罗方法对配分函数的微扰展开进行采样，$Z = Z_{\mathrm{b}} \bar{\sum_{\mathcal{C}} {W ( \mathcal{C} )}}$，如式 (16) 所示。一个配置的总权重由下式给出：

是单个顶点的权重。算符部分 $W [ \hat{h}_{a} ( \tau , \tau^{\prime} ) ]$ 的期望值完全由图 1 中六种顶点类型 $v$ 之一的权重 $W_{v}$ 决定。

对于 XXZ 自旋-玻色子模型，非对角算符部分由下式给出：

而对角部分变为：

这里，我们定义了耦合常数 $\lambda_{\ell} = 2 \alpha_{\ell} \omega_{\mathrm{c}} / s$，其中包含了谱函数的归一化。请注意，我们省略了在时刻 $\tau$ 和 $\tau^{\prime}$ 的单位算符，例如，$C$ 对应 $C \hat{\mathbb{1}} ( \tau ) \hat{\mathbb{1}} ( \tau^{\prime} )$。我们可以将磁场 $h_{z}$ 包含在延迟相互作用中，因为有 $d\omega$ $\begin{array} {r} {\mathcal{I} ( \omega ) \int d \tau^{\prime} P ( \omega , \tau -} \end{array}$ $\tau^{\prime} ) = 1$。我们在 $\hat{h}_{2} ( \tau , \tau^{\prime} )$ 中添加了一个常数 $C$ 以获得正权重。对于图1中定义的不同顶点类型 $v$，其权重由下式给出：

我们选择 $C \geqslant$ max $\begin{array} {r} {[ \frac{\lambda_{z}} {4} , \frac{| h_{z} |} {2} - \frac{\lambda_{z}} {4} ]} \end{array}$ 以确保权重为正。XXZ 自旋-玻色子模型的权重 $W_{v}$ 等价于具有最近邻自旋交换的铁磁 XXZ 模型的权重 [7]。

对于 Jaynes-Cummings 模型，我们有

(29)

顶点权重与式(27)中的类似，但 $\lambda_{z} = 0$ 且 $W_{5} = 0$。

C. 对角更新#

对角更新涉及使用 METROPOLIS-HASTINGS 算法添加或移除单个顶点 $\hat{h}_{2} ( \bar{\tau} , \tau^{\prime} )$。METROPOLIS 接受概率

由接受率决定

它依赖于式(23)中定义的蒙特卡洛权重 $W ( )$ 以及下面指定的提议概率 $T_{0} ( \mathcal{C} \to \mathcal{C}^{\prime} )$。

我们提出以概率密度添加一个新顶点：

其中只有第一时间变量随机选择，而频率和第二时间变量分别从 $\mathcal{T} ( \omega )$ 和 $P ( \omega , \tau - \tau^{\prime} )$ 中采样，如下文所述。注意，有 $n + 1$ 种可能性将新顶点插入到有序列表 $\mathcal{C}_{n + 1}$ 中。我们还包含了一个额外概率 $p_{t_{\mathrm{in}}}$，用于从一组相互作用项中选择相互作用类型 $t_{\mathrm{int}}$，详见第VI节；目前，$p_{t_{\mathrm{int}}} = 1$。对于随机选择一个对角顶点进行移除，我们得到：

其中 $n_{2}$ 是 $\mathcal{C}_{n}$ 中所有对角顶点的个数。利用式(23)中蒙特卡洛权重的比值，$W_{\alpha} ( \mathcal{C}_{n + 1} ) / W_{\alpha} ( \mathcal{C}_{n} ) = \mathcal{I} ( \omega ) P ( \omega , \tau - \tau^{\prime} ) W_{v}$ dω dτ $d \tau^{\prime} / ( n +$ 1)，我们得到 METROPOLIS 比率：

分别对应于添加和移除一个顶点。由于我们将 $\mathcal{T} ( \omega )$ 和 $P ( \omega , \tau - \tau^{\prime} )$ 包含在了建议概率中，因此它们从接受率中消失了。通过这种方式，我们确保了所有参数都具有较高的接受概率。

我们使用逆变换采样法分别从概率分布 $\mathcal{T} ( \omega )$ 和 $P ( \omega , \tau - \tau^{\prime} )$ 中抽取 $\omega$ 和 $\tau - \tau^{\prime}$。给定一个概率分布函数 $f ( x )$，我们计算其累积分布函数 $F ( x ) = \textstyle \int_{0}^{x} d x^{\prime} f ( x^{\prime} )$。如果我们选择一个均匀分布的随机数 $\xi \in [ 0 , 1 )$，那么 $F^{- 1} ( \xi )$ 将返回一个从 $f ( x )$ 中抽取的随机数。在我们的情况中，$\mathcal{T} ( \omega )$ 和 $P ( \omega , \tau - \tau^{\prime} )$ 的累积分布可以通过解析方法求逆。我们首先从式(8)的幂律谱中抽取 $\omega$，即

然后，我们使用选定的 $\omega$ 值从式(14)的玻色子传播子中抽取 $\tau - \tau^{\prime}$，即

为了从对称化传播子 $P_{+} ( \omega , \tau - \tau^{\prime} )$ 中抽取时间差，我们可以推导出一个类似的公式，或者以 $1 / 2$ 的概率将式(37)中的 $\tau - \tau^{\prime} \to \beta - ( \tau - \tau^{\prime} )$ 替换。注意，我们抽取 $\omega \in ( 0 , \omega_{\mathrm{c}} ]$ 和 $\tau - \tau^{\prime} \in [ 0 , \beta )$。

D. 有向环路和虫洞更新#

我们已经在第 IV.A 节中介绍了有向环路更新的概念及其向虫洞更新的扩展。接下来，我们将详细阐述其数学表述。

有向环路更新利用一个扩展的构型空间来连接两个需要大量局部改变才能关联的正则蒙特卡罗构型。对应的世界线构型仅在是否存在一条闭合环路（沿此环路翻转自旋）上有所区别。为了构造这种全局更新，我们需要满足细致平衡条件。文献[7]已证明，由于式(23)中的蒙特卡罗权重是可分解的，这种全局要求简化为每个顶点的一组局部条件，称为有向环路方程。文献[7]中的证明对顶点的结构没有任何假设，因此它也适用于我们的延迟相互作用。因为 $\mathcal{T} ( \omega )$ 和 $P ( \omega , \tau - \tau^{\prime} )$ 是每个顶点的全局前置因子，即它们不依赖于式(24)中的顶点类型，所以它们从有向环路方程中消去了[14]。因此，可以像原始方法[7]一样针对顶点权重 $W_{v}$ 来表述这些方程。在扩展构型空间中，每个顶点被分配了一个有向环路的入口和出口支路，相应的权重变为 $W_{v} ( l_{\mathrm{in}} , l_{\mathrm{out}} )$。有向环路方程为

图3：(a) 对于有向环路更新，自旋-玻色子模型的延迟顶点等价于晶格模型的最近邻顶点。(b) 用于求解有向环路方程的分配表与晶格模型的构建方式相同（详见正文说明）。这里，我们为每个顶点分配了一个有向环路段，并将有向环路更新期间会被翻转的顶点支路标为红色。

对应了细致平衡的局域版本以及概率守恒。这里，类型为 $\bar{v}$ 的顶点是通过对沿指定环路段的支路上的自旋进行翻转，从 $v$ 关联得到的。如果有向环路从支路 $l_{1}$ 进入一个顶点，则 $W_{v} ( l_{1} , l_{2} ) / W_{v}$ 给出了从支路 $l$ 离开的概率。

在解释如何求解有向环路方程之前，我们先重新审视图1中所示的顶点结构。对于有限晶格上的常规自旋模型（例如，海森堡或XXZ模型），每个顶点有四个支路。我们之前已经论证过，在求解有向环路方程时，顶点的非局域性并不起作用。我们发现，如果我们将 $\tau^{\prime}$ 处的子顶点置于 $\tau$ 处子顶点的右侧[如图3(a)所示]，那么图1中的顶点就对应到具有最近邻相互作用的常规自旋模型。因此，我们可以使用与常规自旋模型相同的技术来求解有向环路方程。此外，方程(27)中XXZ自旋-玻色子模型的顶点权重，等价于有限晶格上铁磁XXZ自旋模型的顶点权重。因此，对于这两种情况，有向环路方程具有相同的解。特别地，在SU(2)对称情况下，环路可以确定性地构造[6]。

为使内容完整，我们简要概述一下如何解析求解本文所考虑的自旋-玻色子模型的有向环路方程。在更一般的情况下，可以使用线性规划技术来求解有向环路方程[53]。考虑图3(b)中展示的赋值表。该表的构造方式是，每一行对应方程(39)，非对角元素通过方程(38)相关联。反弹概率在对角线上给出。对于图3(b)中选择的赋值，我们得到方程组

我们可以如下求解 $a , b ,$ 和 c：

我们的指导原则是最小化反弹移动权重 $b_{1} , b_{2} .$ 和 $b_{3}$ ，因为它们会降低算法的效率。如果我们代入XXZ自旋-玻色子模型的权重，我们得到一个显式解，

如果我们调整常数偏移量 C 以及反弹权重 $b_{1}$ 或 $b_{2}$ 中的一个，我们总能得到正定解。在没有外磁场 $h_{z}$ 的情况下，我们可以得到 $\lambda_{z} \leqslant \lambda_{xy}$ 的无反弹解。对于SU(2)对称情况，即 $\lambda_{z} = \lambda_{xy}$ ，我们可以选择 $C = \lambda_{z} / 4$ 来构造环路，使得入口支路精确决定出口支路，这与海森堡模型的情况相同[6]。必须对所有可能的赋值表求解有向环路方程。对于Jaynes-Cummings相互作用，可以推导出类似的解。

E. 观测量#

观测量的计算等同于标准的世界线量子蒙特卡罗方法。下面，我们仅简要概述最常见的观测量。

自旋系统的性质可以通过虚时关联函数获得

以及相应的磁化率

这里，$\Omega_{m} = 2 \pi m / \beta , m \in \mathbb{Z}$ 是玻色子松原频率（bosonic Matsubara frequencies）。我们还定义静态磁化率 $\chi_{\ell} =$ $\chi_{\ell} ( \mathrm{i} \Omega_{0} )$。这些可观测量在 z 方向上的分量不会改变世界线构型，可以直接从传播态 $\textstyle | \alpha_{p} \rangle$ 获得，而对角线外的分量则可以通过在定向环更新过程中追踪定向环的传播来获取。后者之所以可行，是因为环的头部和尾部被识别为自旋翻转算符——具体取决于模型，将它们识别为 $\hat{\sigma}_{x}$ 算符或 $\hat{S}_{+}$ 和 $\hat{S}_{-}$ 算符可能更方便。与 SSE 表示相比，我们可以直接访问每个顶点在虚时域的值，这简化了两种情况下的时移可观测量计算。关于对角和非对角测量的详细讨论，请参考文献[12,54]。

由于玻色子已被积分掉，因此无法直接从世界线构型中获取库（bath）的性质。然而，借助生成函数，可以通过高阶自旋-自旋关联函数推导出玻色子的可观测量，并最终从顶点的分布中恢复出来[55]。例如，$\langle n \rangle =$ $- \beta \langle \hat{H}_{\mathrm{s}} + \hat{H}_{\mathrm{sb}} \rangle$ 将平均展开阶数与自旋子系统能量关联起来。针对单个库频率的玻色子能量、玻色子传播子或比热的估计器已在文献[15,55]中推导出来。它们可以推广到连续库，但需要遵循文献[56,57]定义一个从离散频率到连续谱的额外映射，这给库的性质留下了一些不确定性。

V. 结果#

为了展示我们 QMC 方法的效率，我们考虑 $h_{\ell} = 0 , s = 0.8$ 的双库自旋-玻色子模型中临界相与局域相之间的量子相变，并将其作为自旋-玻色子耦合 α 的函数。对于这些参数，临界耦合已通过变分 MPS 方法确定，该方法基于具有对数离散化的玻色库（bosonic bath）的威尔逊链（Wilson chain）[24,58]。我们证明，我们的 QMC 方法能达到足够低的温度，通过自旋磁化率区分这两个相，并通过有限尺度分析精确估计临界耦合。

对于浴指数 $0 < s < 1$，微扰重整化群预言了一个中间耦合不动点，在该点杂质的部分屏蔽导致一个具有幂律自旋-自旋关联 $C_{xy} ( \tau ) \sim 1 / \tau^{1 - s}$ 的临界相 [25, 26, 59–61]。数值模拟揭示，临界相仅在 $s^{*} < s < 1 \left[ s^{*} \approx 0.76 \right]$ 和 $\alpha < \alpha_{\mathrm{c}}$ 时才稳定，超过该范围则形成满足 $\lim_{\tau \to \infty} C_{xy} ( \tau ) = m_{\mathrm{loc}}^{2} > 0$ 的局域矩 [24, 58]。关于相图、不动点结构以及双库自旋-玻色子模型临界性质的详细讨论见文献 [58]。本文发展的QMC方法使我们能够从有限温度趋近这两个相。图4(a)展示了序参量 $C_{xy} ( \beta / 2 )$ 在 $\alpha > \alpha_{\mathrm{c}}$ 时出现局域矩，而 $C_{xy} ( \beta / 2 )$ 在 $\alpha < \alpha_{\mathrm{c}}$ 时标度为零。临界相与局域相也可通过静态自旋磁化率 $\chi_{xy}$ 和 $\chi_{z}$ 的低温响应加以区分，如图4(b)所示。对于 $\alpha = 0$ 或 $T / \omega_{\mathrm{c}} \gg 1$，我们有 $\chi_{xy} = \chi_{z} = 1 / 4 T$。当温度足够低且处于局域相深处时，$\chi_{xy}$ 趋近居里定律 $\chi_{xy} = m_{\mathrm{loc}}^{2} / T$，其中局域矩 $m_{\mathrm{loc}}$ 有限。进入临界相后，$\chi_{xy}$ 明显偏离居里行为，并缓慢收敛至 $\chi_{xy} \sim T^{- s}$；我们未展示 $\alpha < 0.7$ 时的 $\chi_{xy}$，因为所有曲线均与 $\alpha = 0.7$ 时的 $\chi_{xy}$ 平行，且与其他曲线相交，降低了可区分性。另一方面，$\chi_{z}$ 也呈现幂律依赖，其指数随α增大而稳步减小，最终在局域相中变为零。尽管我们的QMC方法能达到 $T / \omega_{\mathrm{c}} = 10^{-7}$ 的低温，$C_{xy} ( \beta / 2 )$ 和 $\chi_{xy}$ 仍经历有限温度效应，导致缓慢收敛至预期的 $T \to 0$ 行为。

图4：(a) 不同温度下序参量 $C_{xy} ( \beta / 2 )$ 作为自旋-玻色子耦合α的函数。在局域相中，$\begin{array} {r} {\operatorname*{lim}_{\beta \to \infty} C_{xy} ( \beta / 2 ) = m_{\mathrm{loc}}^{2}} \end{array}$ 收敛到一个有限局域矩，而在临界相中 $m_{\mathrm{loc}} = 0$。(b) 不同α下自旋磁化率 $\chi_{xy}$ 和 $\chi_{z}$ 作为温度的函数。我们展示α ∈ $[ 0.7 , 1.5 ]$ 时的 $\chi_{xy}$ 和 α ∈ $[ 0.1 , 0.9 ]$ 时的 $\chi_{z}$，步长 $\Delta \alpha = 0.1$，如色条所示。黑线对应自由局域矩 $\chi = ( 4 T )^{-1}$，虚线表示临界相中的渐近行为 $\chi \sim T^{- s}$。这里 $s = 0.8$

我们的QMC方法足够强大，能够精确估计临界耦合。为此，图5(a)和5(b)展示了重标度磁化率 $T^{s} \chi_{xy}$ 和磁化率比值

两者分别对应。后者是关联比的推广形式，其倒数可关联于虚时方向关联长度的有限尺寸估计量[41]。当 $T \to 0$ 时，这两类可观测量在临界相内保持有限，而 $T^{s} \chi_{xy}$ 在局域相发散，$R_{xy}$ 则趋于零。我们可以通过对 $T^{s} \chi_{xy}$ 和 $R_{xy}$ 进行有限尺寸分析来估计临界耦合——此处逆温度扮演了系统尺寸的角色，但作用于虚时方向。对于每组温度数据 $\{T , T / 10 \}$，我们提取交叉点 $\alpha_{*} ( T )$ 并汇总于图5(c)。在足够低的温度下，预期交叉点遵循幂律行为 $\alpha_{*} ( T ) = \alpha_{\mathrm{c}} + A T^{e}$，由此可通过最小二乘法拟合估计 $\alpha_{\mathrm{c}}$。基于 $T^{s} \chi_{xy}$ 的交叉点，我们估计 $\alpha_{\mathrm{c}} = 0.76213 ( 6 )$；$R_{xy}$ 的 $\alpha_{*} ( T )$ 表现出非单调温度依赖，无法进行可靠拟合，但数据与 $T^{s} \chi_{xy}$ 的估计结果一致。从序参量比 $C_{xy} ( \beta / 2 ) / C_{xy} ( \beta / 4 )$ 也可获得一致交叉点，但较大的统计误差阻碍了对 $\alpha_{\mathrm{c}}$ 的精确估计。我们的临界耦合估计值与先前 MPS 结果 $\alpha_{\mathrm{c}}^{\mathrm{MPS}} = 0.765058 \ [ 58 ]$ 高度吻合；需注意后者由于 Wilson 链中玻色子浴的对数离散化而存在高达1%的系统误差[58]。考虑此 MPS 系统误差后，我们的 QMC 估计给出了更精确的临界耦合值。在后续研究[62]中，我们还确定了该量子相变处的关联长度指数，与 MPS 结果完全一致。

图5：有限温度分析：(a)重标度磁化率 $T^{s} \chi_{xy}$ 和(b)磁化率比值 $R_{xy}$ 随自旋-玻色子耦合α的变化。(c)图(a)和(b)中温度 $\{T , T / 10 \}$ 数据对交叉点 $\alpha_{*} ( T )$ 的温度依赖性。我们通过 $T^{s} \chi_{xy}$ 的幂律拟合 $\alpha_{*} ( T ) = \alpha_{\mathrm{c}} + A T^{e}$（红色虚线）估计临界耦合为 $\alpha_{\mathrm{c}} = 0.76213 ( 6 )$。该结果与文献[58]的变分 MPS 结果高度一致，后者具有高达1%的系统误差（阴影区域所示）。此处 $s = 0.8$

最后，我们再次强调，虫洞更新和玻色子传播子的高效采样使我们能够达到低至 $T / \omega_{\mathrm{c}} = 10^{-7}$ 的温度——这比先前仅能达到 $T / \omega_{\mathrm{c}} \approx 10^{-4}$ 的 QMC 研究低了数个数量级[38,41]。基于我们的算法发展，我们得以在图4中解析出临界相和局域相的特征低温行为，特别是 $\chi_{xy} \sim T^{- s}$ 与 $\chi_{xy} \sim T^{- 1}$ 之间的差异，并在图5中获得了临界耦合的精确估计。关于密切相关的 SU(2)对称自旋-玻色子模型的系统性研究将在文献[62]中呈现。定向环更新对于本文考虑的所有参数区间均高效运行，因为环的平均长度由 $\chi_{xy}$ 给出，因此随问题的特征温度尺度发散。对于最低温度和最强耦合，我们的展开阶数已达到约2000万个顶点。

VI. 其他相互作用类型的推广#

我们已为 XXZ 自旋-玻色子模型和 Jaynes-Cummings 模型引入了 QMC 方法。下面将讨论杂质模型如何包含其他相互作用类型，以及虫洞更新如何应用于晶格模型。

A. 杂质模型#

1. xy平面内的磁场#

对于之前讨论的杂质模型，在xy平面内施加磁场也可能是有意义的。我们可以将附加场分解如下：

由于式(19)和式(20)中的推迟自旋-玻色子相互作用总是包含一对 $\hat{S}_{+}$ 和 $\hat{S}_{-}$ 算符，因此微扰展开中的磁场项也必须成对出现，从而不会导致符号问题。

我们可以将xy平面内的磁场视为一种额外的相互作用类型。为了构建完整的更新过程，我们将对角项 $- C \int d \tau \hat{\mathbb{1}} ( \tau )$ 添加到相互作用顶点。方便地选择常数前因子等于 $h_{xy} / 2$，其中 $h_{xy} = \sqrt{h_{x}^{2} + h_{y}^{2}}$ 。这样，所有顶点权重都具有相同的绝对值，从而简化了有向环方程的解。如果有向环进入这些顶点之一，我们希望将单位算符变换为自旋翻转算符或反之，这会阻止环头继续在局部传播。然而，我们可以选择世界线构型中的一个随机磁场顶点（包括原始顶点），并从那里继续构建环。出射支路的选择必须与传播算符和世界线构型一致。当环头返回到其原始起始点时，环的构建就完成了。以这种方式包含磁场的优点是，我们可以从微扰展开中的非对角磁场顶点数量测量 $\langle h_{x} \hat{S}_{x} + h_{y} \hat{S}_{y} \rangle = \langle n_{x y , 1} \rangle / \beta$ 。

2. XYZ 自旋-玻色子模型#

我们在式(5)中对自旋-玻色子模型的讨论仅限于XXZ情况。在此需要提及的是，也可以模拟具有三个独立耦合 $\lambda_{\ell}$ 的完整问题。此时，式(25)中的非对角顶点必须替换为

新的项 $\hat{S}_{+} ( \tau ) \hat{S}_{+} ( \tau^{\prime} )$ 和 $\hat{S}_{-} ( \tau ) \hat{S}_{-} ( \tau^{\prime} )$ 必须被视为额外的顶点类型 v，即 $v = 7$ 和 $v = 8$。存在这些顶点类型时，有向环可以从顶点的所有四个支路退出，这使得有向环方程的解比以前稍微复杂一些。只要不在xy平面内施加附加场，式(47)中的第二项就不会导致负权重，因为在微扰展开中它必须总是成对出现。

3. 原始自旋-玻色子模型#

为完整起见，我们还想提及如何模拟式(4)中的原始自旋-玻色子模型。由于玻色子热库仅与自旋的z分量耦合，我们无法直接应用上述定义的虫洞更新。为了获得推迟的自旋翻转相互作用，我们可以在 $\hat{S}_{x}$ 基下构建QMC方法。这相当于使用Hadamard矩阵对自旋算符进行旋转（在不存在 $\hat{S}_{y}$ 的情况下交换 $\hat{S}_{x}$ 和 $\hat{S}_{z}$）。相互作用部分变为

并使得在 $C \geqslant | h_{x} | / 2$ 条件下蒙特卡罗权重为正。非对角项现在包含两种额外的顶点类型。因此，有向环可以退出所有四条顶点腿，并且有向环方程可以使用线性规划技术求解[53]。

对于自旋-玻色子模型，文献[23]中已经基于 $\hat{S}_{z}$ 算符之间的推迟相互作用提出了一种聚类算法。自旋-玻色子模型与横向磁场中的长程伊辛模型具有密切相似性，后者可以在SSE表示中进行高效模拟[63]。也可能将这种方法扩展到虚时中的长程相互作用。

B. 晶格模型#

我们为单个自旋耦合耗散玻色子浴发展出了QMC方法，但虫洞更新可以简单扩展到晶格模型。在最简单的情况下，每个晶格点耦合到一个独立的浴，因此相互作用顶点只需增加一个额外的格点指标i。我们的方法细节保持不变，只有在对角顶点提议添加时需要随机抽取一个 $i \in [ 1 , {\cal L} ]$，这会在式(32)中引入一个额外的 $1 / L$ 因子。为了耦合不同晶格点上的自旋，我们可以纳入任何可以在SSE表示中模拟且没有符号问题的相互作用顶点，例如最近邻之间的XXZ自旋交换耦合。额外的顶点 $\hat{H}$ 变换到相互作用表象为 $\begin{array} {r} {\hat{\mathcal{H}} = \int d \tau \hat{H} ( \tau )} \end{array}$。自旋算符获得虚拟时间变量，这些变量在对角更新期间被采样以确定其在世界线配置中的位置。对角更新必须使用第IV C节中的Metropolis方案来表述，而有向环更新的概率表保持不变。注意，不同相互作用类型的顶点结构不必相同，只要每种类型都可以单独制定环更新。最终，本文引入的虫洞更新使我们能够研究耗散对各种晶格模型的影响。虽然虫洞算法可以在对角算符耦合到玻色子浴时应用[30]，但我们使用有向环的公式允许模拟更一般的耦合，这些耦合也包括非对角项。海森堡链耦合欧姆浴的结果将在别处展示[64]。

我们针对在空间上局域但虚时上非局域的玻色子介导相互作用引入了虫洞更新。一般来说，只要蒙特卡罗权重为正，虫洞更新就可以应用于空间或时间上的长程相互作用。量子光学中会出现这样的例子，例如在Dicke模型[65]或离子阱模拟器[18]中，全局玻色子模式耦合到不同格点上的自旋——对于后者，需要根据具体情况评估符号问题的存在与否。然而，我们预期像Jaynes-Cummings-Hubbard模型这样的极化激元模型会在我们的公式中引入符号问题，因为积分掉玻色子跳跃项 $\hat{a}_{i}^{\dagger} \hat{a}_{j}$ 会导出一个非局域玻色子传播子，该传播子在 $i \neq j$[66] 时有负贡献；在这种情况下，直接模拟玻色子更好[67,68]。此外，我们可以设计在空间和时间上长程的相互作用，这些相互作用不一定与自旋-玻色子哈密顿量相关，但可以用虫洞更新进行模拟。

VII. 结论#

我们开发了一种精确的量子蒙特卡洛方法，该方法通过推迟自旋相互作用来模拟与耗散玻色子热库耦合的自旋系统。我们的方法利用了最初在 SSE 表示 [7] 中提出、近期被推广到推迟相互作用 [14] 的有向环更新。为了将推迟的自旋翻转相互作用纳入有向环更新的框架，我们引入了非局域的虫洞更新，该更新允许环在顶点的两个算符之间隧穿。有向环更新的表述仍和 SSE 表示 [7] 一样简单，因为推迟相互作用的时间依赖性已完全由对角更新 [14] 处理。我们的方法适用于任何谱分布的热库，在更新对角更新期间，可以高效地对该谱分布以及玻色子传播子的时间依赖性进行采样。本文提出的思想也适用于相关的量子蒙特卡洛方法，例如使用虫洞/有向环更新来采样世界线构型的虫算法 [5]。

我们针对双热库自旋-玻色子模型展示了我们量子蒙特卡洛方法的效率。与之前的量子蒙特卡洛方法 [38,41] 相比，我们的方法达到了显著更低的温度，这使我们得以识别临界相和局域相中自旋磁化率的特征性有限温度响应。为了估计这两个相之间的量子临界耦合强度，我们比较了不同可观测量在有限尺寸下的行为。我们估计的 $\alpha_{\mathrm{c}}$ 与之前变分矩阵乘积态研究的结果 [24,58] 吻合良好。

在未来的研究中，我们的量子蒙特卡洛方法将使我们能够确定与多个热库 [62] 耦合的自旋-玻色子模型的相图。悬而未决的问题包括这些模型的量子临界性质，以及它们与具有空间长程相互作用的自旋系统的关系。自旋-玻色子相互作用在更复杂的玻色-费米近藤模型中也扮演着重要角色，该模型是我们理解重费米子系统中出现的量子临界行为 [69] 的核心模型。对于本文考虑的局域自旋-玻色子耦合，虫洞更新可以与任何已在 SSE 表示中能够模拟的格点模型相结合，这将使我们能够研究耗散对有限系统 [64] 的影响。我们预期，我们的量子蒙特卡洛方法将能够模拟与更一般的玻色子模式耦合的格点模型，这些模式描述了量子光学或离子阱系统中的光与物质相互作用——自旋-玻色子模型可以在，例如，量子模拟器 [70] 中实现。总的来说，虫洞更新不仅适用于虚时方向的长程相互作用，也适用于空间方向的长程相互作用。因此，我们的算法发展开辟了一条在成熟的有向环算法框架内研究全新一类模型的途径。

致谢#

感谢 F. Assaad、M. Hohenadler、D. Luitz、F. Parisen Toldin 和 M. Vojta 进行了有益的讨论。

附录：实现#

有向环算法在 SSE 表示中的实现已有详细文献记载，关于必要数据结构的详细描述请参考

文献 [7,71]。下面，我们将快速概述在相互作用表示中由于虚时非局域性而需要进行的修改。

在 SSE 表示中，顶点被保存在一个固定长度的算符串中，在对角更新过程中会顺序遍历该算符串；在此过程中，单位算符与对角算符相互交换。在我们的公式中，顶点被保存在一个未排序的列表中，其中每个元素都包含所有顶点变量的结构。对于本文定义的自旋-玻色子顶点，每个顶点具有变量 $\nu = \{t_{\mathrm{int}} , v , \omega , \tau , \tau^{\prime} \}$；需要注意的是，保存顶点类型 v 比保存算符类型 a 更有用。在开始一组对角更新之前，我们遍历一次顶点列表，并建立一个新列表，其中包含在其上应用了非对角算符的格点 i 处的所有时间变量。对该列表排序后，我们可以从初始状态 |α 的知识访问世界线配置中每个位置的自旋配置。尽管排序算法对每个格点需要 (β ln β) 次运算，在数学上比有向环更新更昂贵，但对于我们模拟可达到的所有温度，环的构造仍然是最昂贵的任务。为了添加一个对角顶点，我们可以从排序列表中访问自旋配置，按照第 IV C 节所述绘制顶点变量，计算 METROPOLIS 接受概率，并最终将该顶点作为顶点列表的第 $n + 1$ 个元素添加。为了删除一个对角顶点，我们随机选取一个对角顶点，将其与顶点列表中保存的第 n 个元素交换，然后将计数顶点数的整数 n 减一。虽然对角更新会改变顶点数量，但我们可以使用固定长度的顶点列表，只考虑前 n 个元素。

[1] R. H. Swendsen and J.-S. Wang, Nonuniversal Critical Dynamics in Monte Carlo Simulations, Phys. Rev. Lett. 58, 86 (1987).

[2] U. Wolff, Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett. 62, 361 (1989).

[3] H. G. Evertz, G. Lana, and M. Marcu, Cluster Algorithm for Vertex Models, Phys. Rev. Lett. 70, 875 (1993).

[4] H. G. Evertz, The loop algorithm, Adv. Phys. 52, 1 (2003).

[5] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, Exact, complete, and universal continuous-time worldline Monte Carlo approach to the statistics of discrete quantum systems, J. Exp. Theor. Phys. 87, 310 (1998).

[6] A. W. Sandvik, Stochastic series expansion method with operator-loop update, Phys. Rev. B 59, R14157 (1999).

[7] O. F. Syljuåsen and A. W. Sandvik, Quantum Monte Carlo with directed loops, Phys. Rev. E 66, 046701 (2002).

[8] J. Stolpp, T. Köhler, S. R. Manmana, E. Jeckelmann, F. Heidrich-Meisner, and S. Paeckel, Comparative study of stateof-the-art matrix-product-state methods for lattice models with large local Hilbert spaces without U(1) symmetry, Comput. Phys. Commun. 269, 108106 (2021).

[9] M. Hohenadler and T. C. Lang, Autocorrelations in quantum Monte Carlo simulations of electron-phonon models, in Computational Many-Particle Physics, edited by H. Fehske, R.

定向环更新的实现与文献[7]非常相似。为了创建双向链接的顶点列表，我们首先再次遍历顶点列表，并为时间变量和子顶点标签分别创建独立的列表——后者是顶点编号和算子1或2标签的组合。我们为时间变量创建一个索引列表，用于对子顶点标签列表进行排序。由此，我们可以建立链接顶点列表。随后，定向环的构建遵循文献[7]中的描述。当我们在支腿 $l_{i}$ 处进入类型为 v 的顶点时，可以通过概率表确定出口支腿 $l_{e}$。即使对于虫洞移动，顶点变量和出口支腿的组合也能精确确定链接顶点列表中出口顶点的标签。在构建环的过程中，我们还会更新顶点列表中的顶点配置。请注意，正确测量非对角时间偏移关联函数要求我们为定向环随机选择一个起始时间。在完成一组定向环更新后，我们还会按照文献[7]的描述更新初始状态 $|\alpha\rangle$。

在蒙特卡罗预热过程中，对角更新中的提议次数和定向环更新中构建的环数量会进行调整——使得每个顶点平均至少被触及一次——之后保持不变。我们从随机选取的 $| \alpha \rangle$ 和空的顶点列表开始预热。为了加速预热过程，我们采用 beta 倍增方案。对于某些可观测量的测量，我们可能需要重新对部分顶点进行排序。

Schneider, 和 A. Weiße (Springer, Berlin, Heidelberg, 2008), 第 357–366 页。

[10] F. Assaad 和 H. Evertz, 自旋、声子和电子的世界线和行列式量子蒙特卡洛方法, 载于《计算多粒子物理》（Computational Many-Particle Physics）, H. Fehske, R. Schneider, 和 A. Weiße 编 (Springer, Berlin, Heidelberg, 2008), 第 277–356 页。

[11] H. Suwa, 《量子自旋-声子复杂系统的几何构造马尔可夫链蒙特卡洛研究》（Geometrically Constructed Markov Chain Monte Carlo Study of Quantum Spin-Phonon Complex Systems）(Springer, Japan, 2014)。

[12] A. W. Sandvik, R. R. P. Singh, 和 D. K. Campbell, 相互作用表象中的量子蒙特卡洛：应用于自旋-佩尔斯模型, Phys. Rev. B 56, 14510 (1997)。

[13] R. P. Hardikar 和 R. T. Clay, 半填充和四分之一填充下一维哈伯德-霍尔斯坦模型的相图, Phys. Rev. B 75, 245103 (2007)。

[14] M. Weber, F. F. Assaad, 和 M. Hohenadler, 延迟相互作用的定向环量子蒙特卡洛方法, Phys. Rev. Lett. 119, 097401 (2017)。

[15] M. Weber, F. F. Assaad, 和 M. Hohenadler, 佩尔斯链中的热和量子晶格涨落, Phys. Rev. B 98, 235117 (2018)。

[16] M. Weber, F. Parisen Toldin, 和 M. Hohenadler, 苏-施里弗-赫格模型中的竞争序和非传统临界性, Phys. Rev. Research 2, 023013 (2020)。

[17] M. Weber, 耦合到量子声子的蜂窝反铁磁体中的价键序, Phys. Rev. B 103, L041105 (2021)。

[18] D. Porras 和 J. I. Cirac, 用囚禁离子实现有效量子自旋系统, Phys. Rev. Lett. 92, 207901 (2004)。

[19] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, 和 N. Y. Yao, 用囚禁离子进行自旋系统的可编程量子模拟, Rev. Mod. Phys. 93, 025001 (2021)。

[20] M. O. Scully 和 M. S. Zubairy, 《量子光学》（Quantum Optics）(Cambridge University Press, Cambridge, 1997)。

[21] I. Bloch, J. Dalibard, 和 W. Zwerger, 超冷气体的多体物理, Rev. Mod. Phys. 80, 885 (2008)。

[22] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, 和 W. Zwerger, 耗散两能级系统的动力学, Rev. Mod. Phys. 59, 1 (1987)。

[23] A. Winter, H. Rieger, M. Vojta, 和 R. Bulla, 亚欧姆自旋-玻色子模型中的量子相变：使用连续虚时团簇算法的量子蒙特卡洛研究, Phys. Rev. Lett. 102, 030601 (2009)。

[24] C. Guo, A. Weichselbaum, J. von Delft, 和 M. Vojta, 单浴和双浴自旋-玻色子模型中的临界相和强耦合相, Phys. Rev. Lett. 108, 160401 (2012)。

[25] S. Sachdev, C. Buragohain, 和 M. Vojta, 近临界二维反铁磁体中的量子杂质, Science 286, 2479 (1999)。

[26] A. M. Sengupta, 涨落场中的自旋：Bose(+Fermi) Kondo模型, Phys. Rev. B 61, 4041 (2000).

[27] A. H. Castro Neto, E. Novais, L. Borda, G. Zaránd, 和 I. Affleck, 磁有序系统中的量子磁杂质, Phys. Rev. Lett. 91, 096401 (2003).

[28] E. Luijten 和 H. W. J. Blöte, 长程相互作用自旋模型的蒙特卡洛方法, Int. J. Mod. Phys. C 06, 359 (1995).

[29] P. Werner 和 M. Troyer, 耗散量子系统的团簇蒙特卡洛算法, Prog. Theor. Phys. Suppl. 160, 395 (2005).

[30] Z. Cai, U. Schollwöck, 和 L. Pollet, 识别相互作用自旋-玻色子模型中由浴诱导的玻色液体, Phys. Rev. Lett. 113, 260403 (2014).

[31] A. Georges, G. Kotliar, W. Krauth, 和 M. J. Rozenberg, 强关联费米子系统的动力学平均场理论与无穷维极限, Rev. Mod. Phys. 68, 13 (1996).

[32] A. N. Rubtsov, V. V. Savkin, 和 A. I. Lichtenstein, 费米子的连续时间量子蒙特卡洛方法, Phys. Rev. B 72, 035122 (2005).

[33] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, 和 A. J. Millis, 量子杂质模型的连续时间求解器, Phys. Rev. Lett. 97, 076405 (2006).

[34] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, 和 P. Werner, 量子杂质模型的连续时间蒙特卡洛方法, Rev. Mod. Phys. 83, 349 (2011).

[35] P. Gunacker, M. Wallerberger, E. Gull, A. Hausoel, G. Sangiovanni, 和 K. Held, 使用虫洞采样的连续时间量子蒙特卡洛, Phys. Rev. B 92, 155102 (2015).

[36] P. Werner 和 A. J. Millis, Holstein-Hubbard模型的高效动力学平均场模拟, Phys. Rev. Lett. 99, 146404 (2007).

[37] P. Anders, E. Gull, L. Pollet, M. Troyer, 和 P. Werner, 玻色子的动力学平均场理论, New J. Phys. 13, 075013 (2011).

[38] J. Otsuki, 连续时间量子蒙特卡洛中的自旋-玻色子耦合, Phys. Rev. B 87, 125102 (2013).

[39] J. H. Pixley, S. Kirchner, M. T. Glossop, 和 Q. Si, 赝能隙Bose-Fermi Kondo模型的连续时间蒙特卡洛研究, J. Phys.: Conf. Ser. 273, 012050 (2011).

[40] J. H. Pixley, S. Kirchner, K. Ingersent, 和 Q. Si, 赝能隙Bose-Fermi Anderson与Kondo模型中的量子临界性：费米子与玻色子诱导Kondo破坏的相互作用, Phys. Rev. B 88, 245111 (2013).

[41] A. Cai 和 Q. Si, 具有SU(2)对称性的Bose-Fermi Anderson模型：连续时间量子蒙特卡洛研究, Phys. Rev. B 100, 014439 (2019).

[42] A. J. Bray 和 M. A. Moore, 量子自旋玻璃的副本理论, J. Phys. C: Solid State Phys. 13, L655 (1980).

[43] S. Sachdev 和 J. Ye, 随机量子海森堡磁体中的无能隙自旋流体基态, Phys. Rev. Lett. 70, 3339 (1993).

[44] D. R. Grempel 和 M. J. Rozenberg, 量子随机海森堡顺磁体中的涨落, Phys. Rev. Lett. 80, 389 (1998).

[45] A. Georges, O. Parcollet, 和 S. Sachdev, 量子海森堡自旋玻璃的平均场理论, Phys. Rev. Lett. 85, 840 (2000).

[46] J. Otsuki 和 Y. Kuramoto, 量子自旋系统的动力学平均场理论：磁有序态解的检验, Phys. Rev. B 88, 024427 (2013).

[47] F. F. Assaad 和 T. C. Lang, 图解行列式量子蒙特卡洛方法：投影方案及其在Hubbard-Holstein模型中的应用, Phys. Rev. B 76, 035116 (2007).

[48] M. Hohenadler, F. F. Assaad, 和 H. Fehske, 一维半满带中电子-声子相互作用范围的影响, Phys. Rev. Lett. 109, 116407 (2012).

[49] M. Weber 和 M. Hohenadler, 二维Holstein-Hubbard模型：临界温度、Ising普适性与双极化子液体, Phys. Rev. B 98, 085405 (2018).

[50] A. W. Sandvik 和 J. Kurkijärvi, 自旋系统的量子蒙特卡洛模拟方法, Phys. Rev. B 43, 5950 (1991).

[51] M. Vojta, 杂质量子相变, Philos. Mag. 86, 1807 (2006).

[52] E. Jaynes 和 F. Cummings, 量子与半经典辐射理论的比较及其在束脉泽中的应用, Proc. IEEE 51, 89 (1963).

[53] F. Alet, S. Wessel, 和 M. Troyer, 量子蒙特卡洛模拟的广义定向环方法, Phys. Rev. E 71, 036706 (2005).

[54] A. Dorneich 和 M. Troyer, 使用随机级数展开获取大多粒子系统的动力学, Phys. Rev. E 64, 066701 (2001).

[55] M. Weber, F. F. Assaad, 和 M. Hohenadler, 费米子-玻色子晶格模型的连续时间量子蒙特卡洛方法：改进的玻色子估计量及应用于Holstein模型, Phys. Rev. B 94, 245138 (2016).

[56] R. Bulla, T. Pruschke, 和 A. C. Hewson, 赝能隙费米系统中的安德森杂质, J. Phys.: Condens. Matter 9, 10463 (1997).

[57] R. Bulla, N.-H. Tong, 和 M. Vojta, 玻色系统的数值重整化群及其在亚欧姆自旋-玻色子模型中的应用, Phys. Rev. Lett. 91, 170601 (2003).

[58] B. Bruognolo, A. Weichselbaum, C. Guo, J. von Delft, I. Schneider, 和 M. Vojta, 双浴自旋-玻色子模型：相图与临界性质, Phys. Rev. B 90, 245130 (2014).

[59] L. Zhu 和 Q. Si, Bose-Fermi Kondo模型中的临界局域磁矩涨落, Phys. Rev. B 66, 024426 (2002).

[60] G. Zaránd 和 E. Demler, Bose-Fermi Kondo模型中的量子相变, Phys. Rev. B 66, 024427 (2002).

[61] M. Vojta, C. Buragohain, 和 S. Sachdev, 二维反铁磁体和超导体中的量子杂质动力学, Phys. Rev. B 61, 15152 (2000).

[62] M. Weber 和 M. Vojta, SU(2)对称自旋-玻色子模型：量子临界性、不动点湮灭与对偶性, arXiv:2203.02518.

[63] A. W. Sandvik, 具有任意相互作用的量子Ising模型的随机级数展开方法, Phys. Rev. E 68, 056701 (2003).

[64] M. Weber, D. J. Luitz, 和 F. F. Assaad, 耗散诱导的序：耦合到欧姆浴的S = 1/2量子自旋链, arXiv:2112.02124.

[65] P. Kirton, M. M. Roses, J. Keeling, 和 E. G. Dalla Torre, Dicke模型导论：从平衡态到非平衡态，反之亦然, Adv. Quantum Technol. 2, 1800043 (2019).

[66] M. Weber, F. F. Assaad, 和 M. Hohenadler, 量子Su-Schrieffer-Heeger模型的激发谱与关联函数, Phys. Rev. B 91, 245147 (2015).

[67] P. Pippan, H. G. Evertz, 和 M. Hohenadler, 强关联晶格玻色子和极化激元的激发谱, Phys. Rev. A 80, 033612 (2009).

[68] M. Hohenadler, M. Aichhorn, S. Schmidt, 和 L. Pollet, Jaynes-Cummings-Hubbard模型的动力学临界指数, Phys. Rev. A 84, 041608(R) (2011).

[69] Q. Si, S. Rabello, K. Ingersent, 和 J. L. Smith, 强关联金属中的局域临界量子相变, Nature (London) 413, 804 (2001).

[70] D. Porras, F. Marquardt, J. von Delft, 和 J. I. Cirac, 束缚离子的介观自旋-玻色子模型, Phys. Rev. A 78, 010101(R) (2008).

[71] A. W. Sandvik, 量子自旋系统的计算研究, AIP Conf. Proc. 1297, 135 (2010).

阅读笔记#

一句话概括#

本文提出了一种结合延迟相互作用框架与定向环（directed loop）更新的连续时间量子蒙特卡洛方法，用于模拟自旋-玻色子模型。该方法的关键创新在于引入了虫洞更新（wormhole update）——利用自由玻色子传播子在虚时两个远距离顶点之间建立非局域连接，使环头可以瞬间隧穿，从而克服自旋翻转算符不守恒总$\hat{S}_z$导致的传统定向环无法局部构造的困难。相比仅能使用局域更新的先前QMC方法（最低只能达到$T/\omega_c \approx 10^{-4}$），该方法可达到$T/\omega_c = 10^{-7}$的极低温度，并实现了对双浴自旋-玻色子模型临界耦合$\alpha_c = 0.76213(6)$的高精度估计，与变分MPS结果$\alpha_c^\mathrm{MPS} = 0.765058$（系统误差$\sim 1\%$）高度吻合。

核心论证链#

1. 问题定位：自旋-玻色子模型中的$\hat{S}_\pm$算符不守恒总$\hat{S}_z$量子数，导致环头在世界线中前进时遇到延迟顶点中一对非对角算符$(\hat{S}_+,\hat{S}_-)$间隔任意远的虚时——中间无其他顶点可传播，传统定向环算法失效。现有QMC方法只能退回到低效的局域更新，温度下探能力受限（最低仅$T/\omega_c \approx 10^{-4}$）。

2. 消除玻色子自由度：将哈密顿量拆为$\hat{H} = \hat{H}_\mathrm{b} + \hat{H}_\mathrm{sb}$，在相互作用绘景中展开配分函数。利用玻色子部分在算符中二次型的特点，通过Wick定理精确积分掉玻色子浴，得到自旋自由度中由自由玻色子传播子$D(\omega,\tau-\tau')$介导的延迟相互作用$\hat{\mathcal{H}}_{\mathrm{ret}}$（式(18)）→ 避免了直接采样无界玻色子希尔伯特空间这一传统困难来源。

3. 顶点结构的等价映射：延迟自旋相互作用在微扰展开中产生由两个虚时点$(\tau,\tau')$和一个自由玻色子传播子构成的延迟顶点（图1），每个顶点有四个支路（两个来自$\tau$处的算符，两个来自$\tau'$处的算符）。如果重新标记$\tau'$处的子顶点在$\tau$处的子顶点右侧（图3(a)），则其拓扑结构与有限格点中最近邻自旋模型的顶点完全相同 → 顶点的非局域性不影响定向环方程的局部求解。

4. 权重等价性与环方程的直接套用：求出6种顶点类型的权重$W_v$（式(27)），发现XXZ自旋-玻色子模型的顶点权重恰好等于铁磁XXZ自旋模型的顶点权重 → 已有定向环方程的解析解（包括无反弹解的存在条件）可以直接套用，无需重新推导。

5. 虫洞更新的核心思想：环头从顶点一个算符的左支路进入后，不沿世界线逐点传播（此路径在延迟顶点中不存在），而是直接利用自由玻色子传播子建立的虚连线——该连接作为虫洞——瞬间隧穿到另一个算符的出口支路（图2(d)）。出口支路的选择概率由局部环方程决定 → 实现了虚时中任意远两点之间的全局连接，克服了总$\hat{S}_z$不守恒导致的大尺度构型阻隔。

6. 高效虚时和频率采样通过逆变换实现：对角更新中，使用逆变换采样从谱分布$\mathcal{I}(\omega)$和玻色子传播子$P(\omega,\tau-\tau')$中分别精确抽取$\omega$和$\tau-\tau'$（式(36)-(37)），然后将这些因子包含在提议概率中 → 它们在Metropolis接受率（式(34)）中完全抵消，保证了高接受概率。延迟相互作用的时间依赖性完全由对角更新处理，定向环的构建不受影响。

7. 相变判定与高精度临界耦合估计：在双浴自旋-玻色子模型($s=0.8, h_\ell=0$)中，计算重标度磁化率$T^s\chi_{xy}$和磁化率比值$R_{xy}$在不同温度下的交叉点$\alpha_*(T)$。对$\alpha_*(T)$进行有限温度标度拟合$\alpha_*(T) = \alpha_c + A T^e$ → 外推出$T\to 0$时的临界耦合$\alpha_c = 0.76213(6)$，统计误差显著低于MPS方法的系统误差($\sim 1\%$)。极低的$T/\omega_c = 10^{-7}$使$\chi_{xy} \sim T^{-s}$与$\chi_{xy} \sim T^{-1}$的低温渐近行为得以清晰区分。

实验参数详解#

批判性思考#

1. 顶点结构等价性的适用范围：本文将延迟顶点映射到最近邻顶点（图3(a)），该映射成立的前提是每个顶点仅涉及一个玻色子模式频率$\omega$的传播子。当推广到XYZ自旋-玻色子模型（式(47)）时，顶点类型增加到8种（新增$\hat{S}_+\hat{S}_+$和$\hat{S}_-\hat{S}_-$），环可从全部4条支路以不同权重退出 → 定向环方程失去解析解，必须用线性规划数值求解，计算代价显著上升。

2. 低温极限的实际瓶颈：尽管声称可达$T/\omega_c = 10^{-7}$，但文中披露最低温/最强耦合下展开阶数$n \sim 2\times 10^7$，且"环的构建仍然是最昂贵的任务"、"排序算法每格点需$O(\beta \ln \beta)$运算"。当方法推广到晶格模型（格点数$\mathcal{L}$）时，总顶点数$\propto \mathcal{L} \times n$，在$\mathcal{L} \sim 10^2$、$\beta \sim 10^7$下排序和环构建的复杂度是否仍可行未被论证。

3. 有限温度外推的不确定性：临界耦合估计依赖对交叉点$\alpha_*(T)$的幂律拟合$\alpha_*(T) = \alpha_c + A T^e$。$R_{xy}$交叉点呈现出无法用该公式拟合的非单调温度行为（图5(c)），表明外推形式可能不普适。另外拟合仅使用了$\{T, T/10\}$两组相邻温度的交叉数据点，不同温度区间的选取对外推结果的影响未被系统分析。

4. 符号问题的界限清晰但狭窄：方法要求所有顶点权重$W_v > 0$。这通过适当选择常数$C$（式(27)）和对角元权重正定的谱函数条件$J(\omega) > 0$实现。但对于Jaynes-Cummings-Hubbard模型中的玻色子跳跃项$\hat{a}_i^\dagger \hat{a}_j$，积分掉玻色子后非对角传播子在$i \neq j$时可能有负贡献 → 本文明确承认此局限（第VI B节最后一段），确认方法并非普适。

5. 玻色子可观测量恢复的额外误差：玻色子被积分掉后无法直接从构型获取库性质（如$\langle n \rangle$、比热等）。恢复这些信息需借助生成函数和从离散频率到连续谱的反问题映射（参考文献[56,57]），后者本质上是不适定的，会引入额外的系统误差。本文仅在理论层面讨论，未给出恢复后结果的精度验证。

局限性#

• 仅适用于正定玻色子传播子：方法依赖Wick定理对玻色子迹求值，要求自由玻色子传播子$D(\omega,\tau-\tau') > 0$。对于费米子库或包含玻色子跳跃$\hat{a}_i^\dagger \hat{a}_j$的系统（如Jaynes-Cummings-Hubbard模型），传播子可能出现负贡献导致符号问题（第VI B节末明确提及）。

• 仅验证了双浴模型单一参数点：数值结果仅针对双浴自旋-玻色子模型在$s=0.8$、$h_\ell=0$、无外场的条件下临界-局域相变。对于SU(2)对称的玻色-近藤模型、Jaynes-Cummings模型、有外磁场$h_x, h_y$的情形，论文虽给出了理论推广公式但缺乏QMC数值结果。

• 高效性依赖于无反弹参数区间：$h_z=0$且$\lambda_z \leq \lambda_{xy}$时可达无反弹解（环头确定性地选择出口），定向环更新效率最高。当参数偏离此区间（如$h_z \neq 0$或$\lambda_z > \lambda_{xy}$），必须引入非零反弹权重$b_1, b_2, b_3$（式(42)），环的平均长度因反弹而缩短，算法效率下降。

• 晶格模型扩展的计算复杂度未量化：向晶格模型的推广（第VI B节）使顶点数正比于$\mathcal{L} \times n$。定向环构建需为$\mathcal{L}$个格点分别建立双向链接列表并排序（$O(\mathcal{L} \beta \ln \beta)$），对于大晶格（$\mathcal{L} \sim 10^2–10^3$）结合极低温（$\beta \sim 10^7$），该开销是否能被接受未经评估。

• MPS外推误差对比不对称：论文指出MPS方法因对数离散化有1%系统误差，但QMC的外推$\alpha_*(T) = \alpha_c + A T^e$同样依赖拟合形式的选择和低温度点的精度。文中$R_{xy}$交叉点表现为非单调（图5(c)），无法进行外推拟合，暗示可能有一个未被捕获的误差来源。

关键公式速查#

• $$ \hat{H} = \hat{H}_{\mathrm{s}} + \hat{H}_{\mathrm{b}} + \hat{H}_{\mathrm{sb}},\quad \hat{H}_{\mathrm{b}} = \sum_{\mu} \omega_{\mu} \hat{a}_{\mu}^{\dagger} \hat{a}_{\mu}$$— 自旋-玻色子模型通用哈密顿量，式(1)-(2)
• $$ J_{\ell}(\omega) = 2\pi \alpha_{\ell} \omega_{\mathrm{c}}^{1-s} \omega^{s},\quad 0 < \omega < \omega_{\mathrm{c}} $$— 幂律谱函数（$s=1$欧姆，$0
• $$ D(\omega, \tau) = \frac{e^{-\omega \tau}}{1 - e^{-\beta \omega}},\quad 0 \leq \tau < \beta $$— 自由玻色子传播子，式(14)
• $$ \hat{\mathcal{H}}_{\mathrm{ret}} = - \iint_{0}^{\beta} d\tau d\tau' \sum_{\mu} \hat{\varrho}_{\mu}^{\dagger}(\tau) D(\omega_{\mu}, \tau-\tau') \hat{\varrho}_{\mu}(\tau') $$— 积分掉玻色子浴后的延迟自旋-自旋相互作用，式(18)
• $$ R(\mathcal{C} \to \mathcal{C}_{n+1}) = \frac{\beta W_v}{(n_2 + 1) p_{t_{\mathrm{int}}}} $$— 对角顶点添加的Metropolis接受率（$\mathcal{I}(\omega)$和$P(\omega,\tau-\tau')$已在提议概率中消除），式(34)
• $$ W_5 = W_6 = \lambda_{xy}/2,\quad W_1 = C + \lambda_z/4 - h_z/2,\quad W_2 = W_3 = C - \lambda_z/4,\quad W_4 = C + \lambda_z/4 + h_z/2 $$— XXZ自旋-玻色子模型六种顶点类型的权重，式(27)
• $$ \alpha_*(T) = \alpha_c + A T^{e} $$— 有限温度标度外推公式，用于从交叉点估计零温临界耦合，第V节
• $$ \alpha_c = 0.76213(6) $$— 双浴自旋-玻色子模型($s=0.8$)临界耦合的QMC估计值，第V节

术语对照#

延伸阅读#

• A. Winter et al., Phys. Rev. Lett. 102, 030601 (2009) — 最早使用聚类QMC方法研究亚欧姆自旋-玻色子模型的量子相变，与本文的延迟相互作用方法互补
• M. Weber, F. F. Assaad, M. Hohenadler, Phys. Rev. Lett. 119, 097401 (2017) — 延迟相互作用定向环QMC方法的最早提出（声子系统），本文的方法论基础
• C. Guo et al., Phys. Rev. Lett. 108, 160401 (2012) — 使用变分MPS方法研究双浴自旋-玻色子模型并发现临界相，本文的主要对比基准
• B. Bruognolo et al., Phys. Rev. B 90, 245130 (2014) — 对双浴自旋-玻色子模型相图和临界性质的系统性MPS研究（含Wilson链对数离散化的系统误差讨论）
• J. Otsuki, Phys. Rev. B 87, 125102 (2013) — 使用局域更新的连续时间QMC方法模拟自旋-玻色子耦合，本文在温度可达范围上的主要对比对象

延伸阅读#

• 利用虫洞更新进行自旋-玻色子模型的量子蒙特卡洛模拟 — Weber（2022）提出精确连续时间量子蒙特卡洛方法，通过非局域虫洞更新模拟耗散玻色子浴耦合的自旋系统。