好的，这是根据您的要求，对第7章内容进行完整翻译后的中文译文。

第7章 · 行列式方法#

本章介绍一种用于模拟晶格上相互作用电子的有限温度算法。由于该算法由 Blankenbecler、Scalapino 和 Sugar（1981；Scalapino and Sugar, 1981）提出，因此有时被称为 BSS 算法。该方法使用 Hubbard-Stratonovich 变换将相互作用电子问题转换为与一个虚时相关辅助场耦合的非相互作用问题。因此，它也被称为辅助场方法（auxiliary-field method）。我们在这里使用另一个名称——行列式方法（determinant method），因为到非相互作用电子问题的变换会产生作为统计权重的行列式，这很恰当。有限温度行列式算法是一种通用电子算法，能够计算各种局部可观测量和关联函数。关于零温行列式方法的讨论，请参考附录 I。

7.1 理论框架#

Feynman 和 Hibbs（1965）用构型空间中所有路径的积分来表述量子力学。在实时中，每条路径都会对积分贡献一个由该路径经典作用量决定的相位。两条路径可以发生相长或相消干涉。在经典极限下，只有稳相路径是重要的。实时量子动力学以其大量干涉路径为特征，对重要性采样构成了远超其能力的挑战。另一方面，统计力学涉及虚时路径积分。其积分贡献的幅值呈指数变化，但相位不变。因此，路径积分的贡献主要来自幅值较大的路径。量子蒙特卡罗方法的任务是识别这些重要路径并高效地对其进行采样。

在本章中，我们通过费曼路径积分来解决有限温度下的多电子问题经典化。由此产生的方法通常被称为行列式方法，因为路径的权重可以表示为行列式，这些看起来并不经典的权重，却非常强烈地暗示了费米子态的反称性。高效且稳定地采样这些权重需要特殊技术。我们从一个简要概述开始，以激发经典表示的一般形式和我们需采样的权重。

7.1.1 Hubbard-Stratonovich 变换#

如第5章所述，大多数量子蒙特卡洛算法通过将其矩阵元写为虚时路径积分来处理玻尔兹曼算符 $e^{-\beta H}$

其中 $\Delta \tau ~ = ~ \beta / N_{\tau}$ 。对于有限温度统计力学，配分函数

因此变为对所有在虚时上满足周期性的路径的积分。

在第5章中，我们还了解到，将哈密顿量分解为其各个项 $H = H_{1} + H_{2} + \cdot \cdot \cdot + H_{n}$ 并使用Trotter近似（5.12）

有助于评估哈密顿量指数算符的矩阵元。然后，我们选取一个方便的基矢，在这些基矢下指数算符要么是已知的，要么是容易近似的，从而表达出这些指数因子的矩阵元。我们现在增加一个额外的要求，即这种分解是将问题转化为只需要非相互作用问题矩阵元的第一步。为实现这一点，我们特别处理分解中描述电子之间直接相互作用的那些项。

作为一个示例，我们考虑一个杂质哈密顿量，其中存在一个同时占据杂质轨道的自旋向上和自旋向下电子之间的库仑相互作用。¹ 我们将其写为

¹ 我们将每个相互作用项写为 $U(n_{\uparrow} - 1/2)(n_{\downarrow} - 1/2)$ 而不是通常的 $U n_{\uparrow} n_{\downarrow}$ 。这是方便的选择，因为我们将在下面进行 Hubbard-Stratonovich 变换；$U(n_{\uparrow} - 1/2)(n_{\downarrow} - 1/2) = U n_{\uparrow} n_{\downarrow} - (U/2)(n_{\uparrow} + n_{\downarrow}) + U/4$ ，因此我们对单电子项进行了移位。由于 $H_{1}$ 包含了单电子（“hopping”）项，这种移位可以通过重新定义相互作用项来吸收到 $H_{1}$ 中。

其中 $n_{\sigma}$ 是杂质轨道上的电子数算符。现在，对相互作用项的关键操作是使用恒等式

其中参数 $\lambda$ 满足 $\cosh(\Delta \tau \lambda) = e^{\Delta \tau U/2}$ 。这个恒等式称为 （离散）Hubbard-Stratonovich 变换（Stratonovich，1958；Hubbard，1959）。它通过引入一个辅助伊辛场 $s$ 并对其求和，将库仑相互作用项 $U(n_{\uparrow} - 1/2)(n_{\downarrow} - 1/2)$ 替换为与自旋算符 $(n_{\uparrow} - n_{\downarrow})$ 耦合的项。 式（7.5）右侧的表达式与 $\Delta \tau$ 的三阶修正相同（练习7.1）。 额外的因子 $e^{- \Delta \tau U/4}$ 可以吸收到 $H_{1}$ 的指数化中。 同样，我们可以将 $U(n_{\downarrow} - 1/2)$ 和 $U(n_{\uparrow} - 1/2)$ 项吸收到 $H_{1}$ 中，使得我们可以写出

其中 $B(s)$ 是一个算符，其矩阵元由被非相互作用电子在有效势 $V(s)$ 中传播所决定。对于这种情况，$V(s) \propto s (n_{\uparrow} - n_{\downarrow})$ 。在实际应用中，我们需要处理的矩阵元是 $e^{- \Delta \tau H_{1}}$ 和 $B(s)$ 。 这些矩阵元可以是多体费米子虚时传播子的矩阵元，但虚时传播子是通过在有效场中单电子运动的积分来计算的。

对于具有多个相互作用的更普遍问题，我们可以重复这个过程。在Trotter分解（7.3）中，我们将每个相互作用项与一个辅助场（例如，自旋为 $\pm 1$ 的伊辛变量）配对，从而得到一个包含大量辅助变量的路径积分。在任意一个时刻 $\tau$，相互作用由一组辅助场变量 $\{s\}_{\tau}$ 来表征。 因此，总路径积分可以写为

其中 $B_{\ell}(\{s\}_{\ell})$ 是代表在时间片 $\ell$ 内、由该时间片上的辅助场 $\{s\}_{\ell}$ 所定义的有效势中传播的算符。这些传播子是自由电子的虚时传播子，但处于局域的有效势中。这些算符的迹可以通过Slater行列式的形式来计算，我们将在下面看到。在实践中，我们通常不计算求和（7.7）本身，而是通过对辅助场变量进行采样（通过蒙特卡洛方法）来计算诸如内部能量或关联函数等热力学量。该算法分为以下步骤：第一步是生成有效的单电子传播子；第二步是计算这些传播子的迹，这将得到配分函数或其他期望值。我们将在下面详细讨论这些步骤。

7.2 有限温度算法#

我们现在详细描述有限温度行列式方法。我们从设置哈密顿量开始。设 $H = H_{0} + H_{I}$ ，其中 $H_{0}$ 包含单电子项（如动能和电子-离子势），$H_{I}$ 包含电子-电子相互作用项。对于晶格模型，$H_{0}$ 通常描述电子在晶格格点之间的跳跃（hopping），$H_{I}$ 描述电子之间的排斥或吸引。我们的目标是计算热期望值

我们通过引入虚时路径积分来实现。将逆温度 $\beta$ 分为 $L$ 个时间片，每个时间片长度为 $\Delta \tau = \beta / L$ 。使用Trotter分解，我们有

现在，我们将Hubbard-Stratonovich变换应用于 $e^{- \Delta \tau H_{I}}$ 。 这会将每个相互作用项替换为辅助场的积分（或求和）。对于每个时间片 $\ell$ ，我们引入一组辅助场变量 $\{s\}_{\ell}$ ，使得

其中 $p(\{s\}_{\ell})$ 是某种权重因子（例如，对于离散Hubbard-Stratonovich变换是 $1/2$），而 $B(\{s\}_{\ell})$ 是一个描述在由 $\{s\}_{\ell}$ 决定的势中非相互作用电子传播的算符。因此，配分函数变为

通过对辅助场构型求和（或积分）并取其迹，我们得到了路径积分形式。关键在于，总的算符乘积 $e^{- \Delta \tau H_{0}} B(\{s\}_{\ell})$ 描述的是非相互作用单电子的传播。因此，它们在整个虚时周期 $0 \leq \tau \leq \beta$ 内的作用，可以用一个单一的有效单电子哈密顿量来表征，该哈密顿量是辅助场的函数。记这个乘积为

其中 $\{s\}$ 表示所有时间片上所有辅助场的完整集合。那么配分函数变为

其中 $p(\{s\}) = \prod_{\ell} p(\{s\}_{\ell})$ 。(注意，这里的 ${\mathrm{Tr}}$ 是多体算符的迹，但$U(\{s\})$是单电子算符的乘积。）

由于 $U(\{s\})$ 描述的是非相互作用多电子系统，其迹可以通过Slater行列式计算。具体来说，如果单电子哈密顿量作用在 $M$ 个单电子态上，那么描述 $N$ 个电子的多体系统的传播子 $U(\{s\})$ 的迹可以表示为

其中 $U(\{s\})$ 现在被视为一个 $M \times M$ 矩阵，代表单电子的虚时传播子。这个行列式出现在费米子系统的配分函数中，源于对费米子场实施反周期边界条件（antiperiodic boundary conditions）的路径积分。另一个常见形式是

其中 $B_{\ell}$ 是代表单个时间片内单电子传播子的 $M \times M$ 矩阵。完整的传播子 $U(\{s\})$ 由这些单时间片传播子的矩阵乘积给出。

因此，配分函数 $Z$ 可以写成：

观测量的热期望值类似地可以写成：

其中 $\langle A \rangle_{\{s\}}$ 是在给定辅助场构型 $\{s\}$ 下的自由电子的期望值。这些自由电子期望值也可以使用单电子传播子通过行列式和格林函数来计算。

在实际的蒙特卡罗模拟中，我们对辅助场构型 $\{s\}$ 进行重要性采样，其权重与 $p(\{s\}) \det[ I + U(\{s\}) ]$ 成正比。 这意味着我们需要评估对于每个试探构型更新的行列式。高效地更新行列式是算法成功的关键，通常通过使用基于快速更新的矩阵求逆引理来实现，例如Sherman-Morrison公式（或Woodbury公式）。我们将不再深入讨论这些技术细节，而是强调该算法的本质在于：

1. Trotter分解和Hubbard-Stratonovich变换，将相互作用问题映射到辅助场中。
2. 将费米子传播子的多体迹表示为单电子矩阵的行列式。
3. 使用蒙特卡罗方法对辅助场构型进行重要性采样。

该算法非常强大，因为它允许直接计算各种可观测量，包括单粒子格林函数、双粒子关联函数等。然而，由于行列式的正负号，它也会遇到费米子符号问题（fermion sign problem），尤其是在呈现受阻相互作用（frustrated interactions）的模型中。

7.3 格林函数的计算#

为了计算观测量的期望值，我们需要知道在给定辅助场构型 $\{s\}$ 下的格林函数。单粒子格林函数（也称为虚时格林函数）尤为有用，因为它可以让我们计算密度、动能以及各种关联函数。在这里，我们考虑等时格林函数：

其中 $c_{i}^{\dagger}$ ($c_{i}$) 是在格点 $i$ 上产生（湮灭）一个电子的算符，并且期望值是在给定辅助场构型下、温度为 $\beta$ 的统计系综中取的。对于一个非相互作用的系统，等时格林函数 $G$ 可以由传播子 $U(0, \tau)$ 和一个 $M \times M$ 的矩阵 $R(\tau)$ 给出，其中 $R(\tau)$ 将 $\tau=0$ 时的态传播到 $\tau$ 时的态。具体地说，对于一个具有 $N$ 个电子的系统（或更一般的巨正则系综），我们有

其中 $U(\beta, 0)$ 是从 $\tau=0$ 到 $\tau=\beta$ 的整个时间演进算符，对于这个构型它是已知的。在虚时 $0$ 和 $\beta$ 之间任何中间时间 $\tau$ 上的格林函数 $G(\tau)$ 也可以通过类似的关系得到。

在模拟中，我们当前有一个辅助场构型，并且需要评估格林函数。由于 $U(\beta, 0)$ 是许多小时间片的矩阵乘积，我们可以迭代地构建它。然而，每次蒙特卡洛更新都会改变该乘积中的一个因子 $B_{\ell}$ ，因此我们需要高效地更新格林函数。这通常通过以下方法完成：

• 使用 Sherman-Morrison 公式（或 Woodbury 公式）对矩阵进行秩更新。
• 或者，通过每次更新后对整个时间演进算符进行逐步乘积来重新计算格林函数。

第一种方法在大多数实现中更受欢迎，因为它高效且稳定。格林函数 $G$ 在计算各种观测量的期望值（如密度、动能、自旋-自旋关联函数等）时会被反复用到。

7.4 抽样和更新策略#

算法的核心部分是对辅助场构型进行蒙特卡罗抽样。通常，我们使用 Metropolis 类型的算法。在这个过程中，我们选择某个时间片 $\ell$ 上的某个辅助场变量 $s_{\ell}$ ，并尝试翻转它（对于离散伊辛场）或改变其值（对于连续场）。这个尝试的更新被接受或拒绝的判据基于配分函数的权重 $p(\{s\}) \det[ I + U(\{s\}) ]$ 的比值。

具体来说，如果 $W(\{s\}) = p(\{s\}) \det[ I + U(\{s\}) ]$ 是权重，并且我们提议从旧构型 $\{s\}$ 变为新构型 $\{s'\}$，则接受概率为

因为从 $\{s\}$ 到 $\{s'\}$ 的转变只影响少数的矩阵因子 $B_{\ell}$ ，所以权重比的数值计算，以及随后对格林函数的更新，都可以非常高效地完成。实际上，计算行列式的比值涉及到求解一个小型的矩阵问题，这可以通过每次被改变的辅助场变量的数量来定标。

这个算法的一个主要挑战是数值稳定性。随着虚时 $L$ 的增加（即 $\Delta \tau \to 0$），乘积 $U(\beta, 0)$ 中的矩阵可能会变得病态（ill-conditioned），因为指数上正的和负的特征值都会增长。这可能导致行列式和格林函数计算中的数值误差不断积累。为了保持数值稳定性，通常采用以下方法：

• 矩阵求逆的高斯消去法：对于小系统，直接对矩阵采取高斯消去法进行计算是稳定的。
• 周期性再正交化（periodic reorthogonalization）：每隔几个时间片，对时间演进算符 $U$ 的矩阵因子进行QR分解或类似操作，以保持其数值正交性。这是最常用的稳定技术之一。
• U 分解（U decomposition）：将传播子矩阵分解为上三角和下三角矩阵的乘积，以用于稳定计算。

通过结合这些稳定性技术，该算法可以模拟各种晶格模型，包括 Hubbard 模型、$t-J$ 模型和电子-声子相互作用模型。它已被广泛应用于研究强关联电子系统中的相图、超导性和量子临界行为。

7.5 测量可观测量#

一旦辅助场构型按照其权重分布被采样，我们就可以测量各种可观测量。在非相互作用参考系中，观测量期望值的计算通常归结为计算单粒子格林函数及其相关量。让我们列举一些重要的例子。

1. 内部能量 $\langle H \rangle$：能量可以用动能和势能的和来表示。对于动力项（hopping），能可以通过计算动能算符的期望值来获得，而这通常需要已知单粒子格林函数。势能项（如 Hubbard $U$）的期望值也可以通过格林函数计算。在给定辅助场构型下，内部能量可被表达为格林函数的函数。

2. 双粒子关联函数：例如密度-密度关联函数 $\langle n_{i}(\tau) n_{j}(0) \rangle$ 或自旋-自旋关联函数 $\langle S_{i}^{z}(\tau) S_{j}^{z}(0) \rangle$。这些更复杂的量可以通过虚时格林函数和Wick定理计算。由于系统在固定辅助场下是自由的，所以关于算符序列的任何期望值都可被分解为单粒子格林函数的乘积。

3. 静态结构因子：它是与密度或自旋相关的序参量，可用于研究不稳定性（如电荷密度波或自旋密度波）。通过傅里叶变换关联函数，我们可以在倒易空间中研究这些静态性质。

在测量过程中，我们通常收集辅助场构型对应的算符值，并对它们进行平均。由于采样是基于其统计权重进行的，算符样本的平均值就是热力学期望值，前提是没有符号问题。符号问题的出现是因为权重 $\det[ I + U ]$ 可能是负的（对于某些构型而言）。这会导致期望值的分母和分子都可能是两个相差很大的小量的商，造成统计误差巨大。当符号问题严重时，该算法的效率骤降。

7.6 局限性和符号问题#

虽然有限温度行列式方法功能强大，但它有重要的局限性。

• 符号问题：对于一般的相互作用模型，特别是那些在时间反演对称性下不保序的模型（如包含无自旋（spinless）费米子、或具有阻挫六角格点自旋的 Hubbard 模型），权重 $\det[ I + U ]$ 可能取负值。这在应用重要性抽样时带来了困难，因为负的权重无法作为概率密度。这个问题被称为费米子符号问题。当符号问题严重时，信噪比随系统尺寸和逆温度 $\beta$ 呈指数下降，使得对大系统或低温的模拟几乎不可能。

• 数值发散：如前所述，由于乘积 $U(\beta, 0)$ 可能导致非常大的数值动态范围，因此必须采用特殊技术（如周期性再正交化）来保证数值稳定性。对于非常低温或非常大的系统，计算资源的开销可能会非常高。

• 系统大小和温度约束：在计算机内存和时间的限制下，该算法通常用于研究中等系统尺寸（例如，$4 \times 4$ 到 $16 \times 16$ 格点）的模型，且逆温度 $\beta$ 上界也有限。这是由矩阵乘法的计算复杂度以及行列式更新的 $O(M^2)$ 或 $O(M^3)$ 开销所决定的。

尽管如此，该算法仍是研究强关联电子系统中最可靠的工具之一。它在 Hubbard 模型和类似系统中展示了反铁磁序、d-波超导等重要物理现象，并提供了在其他方法中无法直接访问的微观信息。

练习#

1. 练习 7.1：证明 Hubbard-Stratonovich 变换 (7.5) 在 $\Delta \tau U$ 的三阶以内成立。检查其 $O((\Delta \tau)^3)$ 误差的符号。

2. 练习 7.2：考虑具有 $M$ 个单电子态和 $N$ 个电子的非相互作用系统。证明巨配分函数可以写成 $Z = \det[ I + e^{-\beta (H_0 - \mu)} ]$，其中 $H_0$ 是以矩阵形式表示的单电子哈密顿量。将这个关系推广到具有辅助场的虚时路径积分形式。

3. 练习 7.3：推导一个简单单能级系统（$M=1$，$N$ 可以是 0,1,2）的配分函数的路径积分表示。显式计算相互作用部分 $U n_{\uparrow} n_{\downarrow}$ 的 Hubbard-Stratonovich 变换，并与直接对角化得到的精确结果进行比较。

其中 $H_{1}$ 是无相互作用问题的哈密顿量，即它是电子产生和湮灭算符的二次型。更深入地，这一表述还意味着 $\begin{array} {r} {H_{1} = \sum_{\sigma} H_{1}^{\sigma}} \end{array}$，且不同的 $H_{1}^{\sigma}$ 算符互相对易。通过 Trotter 近似 (5.12)，

其中 $n_{\sigma}$ 是自旋 $\sigma = \mathrm{\uparrow}$ 或 ↓ 的电子数算符。

将相互作用项与其余部分分离后，我们将该指数变换为辅助场集合上的积分。如果库仑相互作用 $U$ 是排斥性的，则将指数 (7.5) 的自变量写为

然后，我们进行以下 Hubbard-Stratonovich 变换 (Fetter and Walecka, 1971; Negele and Orland, 1988; Fulde, 1991; Enz, 1992)

从而将辅助场 $x$ 引入分析。Hubbard-Stratonovich 变换仅仅是平移高斯积分²的算符推广

变换的结果是辅助场变量上的高斯加权积分，其中被积函数是指数形式的算符，该算符与辅助场成正比，并且是电子产生和湮灭算符 $( n_{\sigma} ~ = ~ c_{\sigma}^{\dagger} c_{\sigma} )$ 的二次型。因为它是二次的，所以不显式表达电子间的相互作用。简而言之，该变换用耦合到杂质轨道净自旋 $n_{\mathord{\uparrow}} - n_{\mathord{\downarrow}}$ 的标量（玻色）场替代了杂质轨道中的排斥库仑相互作用。变换后，$U$ 为正值的效果是在该轨道上产生一个涨落的自旋矩。我们选择的变换形式确保辅助场始终为实数。若未如此选择，则会产生复值被积函数。

这个结果会有问题，因为我们希望将被积函数解释为概率密度。

如果 $U$ 是吸引性的，我们则写为

从而得到

在这里，辅助场耦合到一个电荷自由度（charge degree of freedom）。经过这样的变换后，$U$ 为负值的效果是在关联轨道中产生电荷涨落。

Hirsch（1983）提出了变换的另一种形式，它对于库仑相互作用表示为数算符乘积之和的格点模型尤其有用。Hirsch 变换引入了一个离散辅助场，在实际采样中它比连续场更高效。他的变换如下：

其中 $\cosh \alpha = \exp ( \Delta \tau \vert U \vert / 2 )$。同样，对于 $U > 0$，辅助场耦合到一个自旋自由度（spin degree of freedom）；而对于 $U < 0$，它则耦合到电荷。附录 G 中讨论了更多的 Hubbard-Stratonovich 变换。

无论是引入连续还是离散的辅助场，配分函数都变成了一个嵌套求和：

其中 $C = \{x_{1} , x_{2} , \ldots , x_{N_{\tau}} \}$ 表示虚时上所有可能的辅助场构型集合，并且

这里 $e^{- \beta H}$ 通过 Trotter 近似（7.3）离散化。该变换将每个时间片的相互作用分解为无相互作用项之和（每个自旋对应一项），从而产生一个有效的无相互作用哈密顿量。由于该哈密顿量的每个自旋部分都对易，其指数会分解为指数函数的乘积，每个自旋对应一组乘积。

求迹中的状态 $| \psi \rangle$ 通常是占据数基（occupation number basis）。在占据数基中，每个基态是自旋向上和自旋向下状态的直积：

其中 $\begin{array} {r} {n_{i \sigma} = 0 \mathrm{or} 1 , \sum_{i} n_{i \sigma} = N_{\sigma}} \end{array}$，且

为了符号紧凑，我们写作

来表示具有 $N_{\sigma}$ 个电子的基态，并使用

来表示具有 $\begin{array} {r} {N_{\mathrm{elec}} = \sum_{\sigma} N_{\sigma}} \end{array}$ 的直积态。我们看到，在这个基中，自旋状态也是可分解的。最终，

重要的是要注意，权重乘积中的自旋因子共享同一个辅助场构型。正是这种共享重建了相互作用。虽然这里讨论的杂质问题特别简单，但行列式方法将更一般、更复杂的哈密顿量表达为相同的形式。

通过将电子产生和湮灭算符中的二次项放入指数函数的参数中，Hubbard-Stratonovich 变换使我们能够在固定辅助场变量构型的情况下对电子自由度求迹。换句话说，它将给定问题转化为我们知道如何求解的问题。如果我们明确地对电子自由度求迹，我们会得到

以新权重（记为 $w_{C}^{\sigma}$，无参数）表示，这些权重可写成行列式形式。我们需要指定的正是这些权重。

尽管这是一个进展，但并非我们最终所需的结果。起始的配分函数 (7.2) 是在正则系综中，固定电子数 $N_{\mathrm{elec}}$ 的。

在有限温度下，使用巨正则系综来处理电子更为方便：

这里需要对所有可能的电子数求和。显然，我们可以对每个电子数重复刚才的分析，并为每个电子数找到一个行列式权重。我们还将看到，这些行列式的和是另一个行列式，即费米-狄拉克配分函数的行列式表示。

在下一小节中，我们将推导这个巨正则系综下的行列式权重。其核心是虚时单电子传播子。借助它，我们能够构造虚时多电子传播子，而利用这个传播子，我们又可以构建虚时单粒子格林函数。这一构造很重要：下一小节之后，我们将开始介绍行列式方法，其中采用辅助场构型的Metropolis采样。这种采样需要计算两个行列式的比值，以决定是否接受或拒绝辅助场的提议变化。为了高效地进行采样，我们需要使用特殊方法，而这些方法依赖于虚时单粒子格林函数。拥有这个格林函数还非常方便，因为它也有助于计算多体关联函数。行列式方法的一个优势在于，可以轻松计算多种多样的此类函数。

7.1.2 行列式权重#

从路径积分的角度来看，(7.1) 式中玻尔兹曼因子的矩阵元是由一系列矩阵乘法构成的，每个矩阵代表从一个时间步到下一个时间步的虚时演化。Trotter近似 (7.3) 和哈伯德-斯特拉托诺维奇变换 (7.6)、(7.8) 或 (7.9) 将这个传播的复杂性简化为自由电子态的传播。行列式方法数值框架的核心是那些描述多电子态从虚时 $\tau_{1}$ 到 $\tau_{2}$ 自由传播的矩阵，它们建立在一个矩阵B之上，而B的元素 $B_{i j} ( \tau_{2} , \tau_{1} )$ 则描述单个电子从 $\tau_{1}$ 时刻的状态 $j$ 传播到 $\tau_{2}$ 时刻的状态 $i$ 。对于给定的自旋 $\sigma$ 和辅助场构型，这个单粒子传播子定义为

其中 $\mathcal{T}$ 是编时算符，$H ( \tau )$ 是依赖于辅助场的无相互作用哈密顿量（$\textstyle ( \int d \tau H ( \tau )$ 是欧几里得作用量），|0>是真空态。这个传播子代表了电子产生和湮灭算符运动方程的解³

我们需要注意到 $c_{i} ( \tau )$ 与 $c_{i}^{\dagger} ( \tau )$ 并非厄米共轭关系，因为当前工作在虚时框架下。这意味着传播过程并非幺正的，因此最初正交归一的轨道将失去此性质。

描述一个 $N_{\sigma}$ 电子系统需要使用多电子传播子：

该表达式是一个矩阵的行列式，其矩阵元为单粒子传播子。我们在附录H中推导了这一恒等式，此处先专注于其物理含义。由于演化中的电子是独立的，我们预期多体传播子应为单粒子传播子的乘积。我们在时间 $\tau_{1}$ 将 $N_{\sigma}$ 个全同粒子注入轨道 $j_{1} , j_{2} , \ldots , j_{N_{\sigma}}$，但不知道在时间 $\tau_{2}$ 哪些特定粒子到达了轨道 $i_{1} , i_{2} , \ldots , i_{N_{\sigma}}$。该行列式对所有可能的到达方式进行求和，并配以适当的负号以反映全同电子交换为偶数次或奇数次。

现在回到巨正则系综中多电子配分函数的表达式。Hubbard-Stratonovich 变换（7.9）对所有辅助场的构型C进行求和，并结合 $\Delta \tau = \beta / N_{\tau}$ 下的 Trotter 近似（7.3），生成自旋分解式：

其中 $M _ { i } ^ { \sigma }$ 是定义每个自旋的辅助场依赖型非相互作用哈密顿量的矩阵（附录H）。对于内层求和，反复运用 (7.17) 可得：

在上述式中，$B_{i j}^{\sigma} = B_{i j}^{\sigma} ( \beta , 0 )$（附录 H）。通过直接展开容易验证，右边等于

因此配分函数 (7.18) 可简化为

其中

如所承诺，我们已证明权重为行列式。注意

是一串矩阵乘法（附录 H）。

7.1.3 单粒子格林函数#

单体等时格林函数 (7.22) 的矩阵元的标准多体物理定义为$^{4}$

其中期望值是对巨正则系综取的。为了将该定义与 (7.17) 联系起来，我们首先将作用量与算符 $c_{j \sigma}^{\dagger} c_{i \sigma}^{\vphantom \dagger}$ 在虚时某一瞬时的耦合

通过引入该相互作用，(7.15) 和 (7.20) 给出一个修正的权重，

其中 $Q$ 除 $Q_{j i} = 1$ 外均为零矩阵。由此表达式可得

我们现在利用恒等式

以及迹的循环性质 $\mathrm{Tr} A B C = \mathrm{Tr} B C A = \mathrm{Tr} C A B$ ，得到

因此，

这个格林函数在行列式方法中起着核心作用，因为它有两个不同的关键用途。第一个用途是辅助场采样。对于轨道 i 和虚时 τ，变量 $x_{i} ( \tau )$ 仅与 $n_{\mathord{\uparrow}} ( \tau )$ 和 $n_{\downarrow} ( \tau )$ 耦合，这些信息蕴含在 $G_{i i}^{\sigma} ( \tau , \tau ) = 1 - \langle n_{i \sigma} ( \tau ) \rangle$ 中。第二个用途是测量算符的期望值。模拟的非相互作用特性使得威克定理（Wick's theorem）（Negele and Orland, 1988）能够将这些期望值简化为格林函数乘积的和。

7.2 有限温度算法#

本章的方法最方便应用于诸如以下相互作用的电子模型：

其中库仑相互作用表示为数算符乘积的和（密度-密度相互作用）。⁵ 在 (7.25) 中，自旋依赖的对称矩阵元 $T_{i j}^{\sigma} ~ = ~ T_{j i}^{\sigma}$ 是跃迁积分，$c_{i \sigma}^{\dagger}$ 和 $c_{i \sigma}$ 是位点 i 上自旋为 $\sigma = \mathrm{\uparrow}$ 和 $\downarrow$ 的电子的产生和湮灭算符。库仑相互作用由与自旋无关的对称矩阵 $V_{i j} = V_{j i}$ 参数化，局域电荷密度为 $n_{i} = n_{i \uparrow} + n_{i \downarrow}$ 。这个哈密顿量非常通用，包括了诸如哈伯德模型和拓展哈伯德模型、单杂质和格点安德森模型，以及 Falicov-Kimball 模型等著名模型。

作为实例，我们采用排斥性哈伯德模型（repulsive Hubbard model），

其中 $U > 0$ 是局域库仑相互作用的强度，$\mu$ 是化学势，跃迁矩阵与自旋无关。虽然在哈伯德模型中，跃迁积分通常仅限于最近邻且值为 t，但并没有先验的必要性来假定这一点。在下文中，我们使用 N 表示格点位数，在占据数基矢下工作，并假设周期性边界条件。

7.2.1 矩阵表示#

使用 Trotter 近似 (7.3)，我们首先将单时间步的时间演化算符分解为两个因子

右边第一个因子中哈密顿量的两个部分在占据数表示中都是对角的，因此它们对易。因此，这个因子可以写为

对于右边的第一个乘积，我们现在对每个格点 i 应用离散 Hubbard-Stratonovich 变换 (7.9)。这个变换在辅助场上增加了一个求和，但将二次数算符放入了指数参数中，并允许我们将每个位点的指数分解为对易的自旋向上和自旋向下项。

对于给定的辅助场变量构型，我们有

其中

且

指数参数为 $\sigma \alpha x = \pm \alpha x$，取决于 $\sigma \ : = \mathrm{\uparrow}$ 或 $\downarrow$。参数 $\alpha$ 由双曲余弦 $\cosh ( \alpha ) = \exp ( {\textstyle{\frac{1} {2}}} \Delta \tau U )$ 确定。矩阵 $A^{\sigma} ~ ( 7.31 )$ 是对角矩阵，因为辅助场只耦合到局域电子自由度。在虚时演化中，$\exp ( \Delta \tau T )$ 的效果是扩散性的，$A^{\sigma} ( \tau )$ 像一个外散射势阱，而 $\exp ( \Delta \tau \mu )$ 是一个衰减因子。我们回顾一下，Trotter 近似 (7.27) 使得 $B^{\sigma}$ 成为一个非对称矩阵。

我们可以通过求 T 的特征值 $\lambda_{i}$ 和特征向量 $\nu_{i}$ 来计算 $\exp ( \Delta \tau T )$。这些特征向量构成将 T 对角化的相似变换矩阵 $V = ( \nu_{1} , \nu_{2} , \ldots , \nu_{N} )$ 的列。通过这一变换，

或者，我们可以通过棋盘分解 (checkerboard decomposition) (第 5.2.4 节) 来近似 $\exp ( \Delta \tau T )$，

其中 $T^{( i j )}$ 是稀疏矩阵，只有 $T_{i j}^{( i j )} = T_{j i}^{( i j )} = T_{i j}$ 非零，

因此

每个 $\exp ( \Delta \tau T^{( i j )} )$ 是一个稀疏矩阵，只有 $i i , i j , j i ,$ 和 $j j$ 元素与单位矩阵不同。这种构造跳跃矩阵近似指数的方法在 $T$ 本身稀疏时特别有用，例如当跳跃只存在于最近邻之间时。在这种情况下，稀疏性将矩阵-矩阵乘法中的 $N^{3}$ 次运算减少为 $N_{\mathrm{bonds}} N^{2}$，其中 $N_{\mathrm{bonds}}$ 是与格点 $i$ 相连的格点 $j$ 的数量。在大多数类哈伯德模型中，$N_{\mathrm{bonds}} \ll N$。然而，每步计算时间的减少必须与使用更小的 $\Delta \tau$ 以保持Trotter近似精度的需求相平衡，因此需要采取更多的虚时步骤。

7.2.2 Metropolis 采样#

确定了辅助场的统计权重（7.20）后，我们现在要对它们进行采样。这里，我们采用Metropolis算法⁶（第2.5.1节），在每个虚时，遍历格点，提议翻转伊辛型离散辅助场变量 $x$。我们根据比值 $\mathcal{R} = w_{C^{\prime}} / w_{C}$ 是否大于从均匀分布中选取的随机数 $\zeta$ 来接受或拒绝该提议。该提议仅通过（7.31）中定义的矩阵 $A$ 影响时间演化。对于在时间 $\tau$、格点 $i$ 处从 $x$ 到 $x^{\prime}$ 的翻转，

其中矩阵 $\Delta^{\sigma} ( i , \tau )$ 除 $\Delta_{i i}^{\sigma} ( i , \tau )$ 元素外全为零，该元素为

伴随而来的传播子的变化是

因此

这个表达式看起来可能有点难以计算，但利用两个标准的行列式关系，该表达式可以大大简化。一个关系是

另一个是

代数运算很直接。首先，我们将（7.33）写为

然后，我们利用关系（7.34）将该方程重写为

接下来，我们利用（7.35）写出

最后，利用

其中 $A \equiv B^{\sigma} ( \tau , 0 )$ 且 $B \equiv B^{\sigma} ( \beta , \tau )$，我们得到

其中

$G^{\sigma} ( \tau , \tau )$ 是等时单粒子格林函数（7.24）。由于 $\Delta^{\sigma}$ 的稀疏性，我们有

这是一个标量而非矩阵表达式，只要已知 $G_{i i}^{\sigma} ( \tau , \tau )$，其计算就极为简单。对于每次被接受的 $x_{i} ( \tau )$ 翻转，从（7.29）和（7.37）计算 $G^{\sigma}$ 需要相当大的计算量。然而，矩阵代数提供了一种技术，可以将计算量减少 N 倍。下面我们来推导这一技术。

在提议翻转 $x_{i} ( \tau )$ 的符号时，

其中我们用 $B^{\sigma}$ 表示 $B^{\sigma} ( \tau , 0 ) B^{\sigma} ( \beta , \tau )$ 的简写。电子算符的运动方程（7.16）给出

这是一对与 $B^{\sigma}$ 的传播子特性一致的关系。类似地，对于 $\tau^{\prime} > \tau$，我们通过以下方式将等时格林函数推进到更晚的时刻：

我们仍需展示，如果 Metropolis 提议被接受，如何高效地更新 $G^{\sigma} ( \tau , \tau )$。

解析地，⁷

其中 $\mathcal{R}^{\sigma} ~ ( \sigma = \uparrow , \downarrow )$ 是翻转前后的权重比。用紧凑的记法，

因此，在提议的翻转下，

用矩阵元表示，

（7.44）的意义变得更加清晰：如果我们知道 $G_{i j}^{\sigma}$ 并且接受了翻转，我们就能够以数量级为 $G^{\sigma}$ 元素数量（即 $N^{2}$）的算术运算次数，计算出辅助场新构型下的格林函数。这个数量比与矩阵求逆相关的 $\mathcal{O} ( N^{3} )$ 次运算少了 N 倍。这一节省相当显著，尤其是因为 $G^{\sigma}$ 的更新位于算法的内循环中。

7.2.3 算法#

现在我们讨论该算法及其相关的若干问题。整体蒙特卡洛算法结构与第 3.9 节讨论的相同。在初始化阶段，我们不仅将每个时间 τ 上每个格点 i 的辅助场变量 $x_{i}$ 设为 ±1 的值，还要构建所有时间步的矩阵 $B^{\sigma} ( \tau , \tau - \Delta \tau )$ 以及两个

格林函数矩阵在初始时刻的自旋分量。我们可以利用(7.37)构建这个格林函数，即先执行所有矩阵乘法，然后计算所需的矩阵逆运算。这些步骤看似非常直接，但如下文所述，它们隐藏着复杂性。

算法23描述了一次蒙特卡罗扫描。为了便于讨论扫描过程，我们为第 $i$ 个虚时步 $\tau_{i} = i \Delta \tau$ 引入更方便的记号 $G_{i} \equiv G ( \tau_{i} , \tau_{i} )$ 和 $B_{i}^{\sigma} \equiv$ $B^{\sigma} ( \tau_{i} , \tau_{i - 1} )$ 。利用这些记号：

算法23 有限温度行列式蒙特卡洛算法。 输入：$\overline{{N_{\tau} \mathrm{~} \times \mathrm{~} N}}$ 个辅助场变量 $x_{i , j}$ ，$\overline{{N \times N}}$ 矩阵 $B_{i}^{\sigma}$ ，以及 $N \times N$ 矩阵 $G_{1}^{\sigma}$ for i = 1 to $N_{\tau}$ do if mod(i, $n_{\mathrm{stab}} ) {=} 0$ then for $\sigma = \uparrow$ to ↓ do 用(7.46)计算 $G_{i}^{\sigma}$ ； end for end if for $j = 1$ to N do 用(7.43)计算 $\mathcal{R}^{\sigma} ( x_{i j} )$ ； 在 [0, 1] 内均匀随机抽取 ζ ； if $\mathcal{R}^{\sigma} > \zeta$ then $x_{i j} \gets - x_{i j}$ ； for $\sigma = \uparrow$ to ↓ do 用(7.45)更新 $G_{i}^{\sigma}$ ； end for end if end for if $i < N_{\tau}$ then for $\sigma = \uparrow$ to ↓ do $G_{i + 1}^{\sigma} \gets B_{i + 1}^{\sigma} G_{i}^{\sigma} [ B_{i + 1}^{\sigma} ]^{- 1}$ ； end for else for $\sigma = \uparrow$ to ↓ do $G_{1}^{\sigma} \gets B_{1}^{\sigma} G_{N_{\tau}}^{\sigma} [ B_{1}^{\sigma} ]^{- 1}$ end for end if end for return $G_{i}^{\sigma}$ 其中 $i = 1 , 2 , \dots , N_{\tau} .$

为了更新辅助场，我们使用 Metropolis 算法，并在每个虚时访问每个格点。如果提议的改变被接受，我们就使用 (7.44) 更新格林函数。访问完所有格点后，我们使用 (7.41) 将等时格林函数推进到下一个时刻，然后重新访问格点。在完成所有时刻的采样后，一次扫描（sweep）结束。我们注意到 (7.37) 意味着 $G^{\sigma} ( \beta , \beta ) ~ = ~ G^{\sigma} ( 0 , 0 )$，并且 (7.41) 将 $G^{\sigma} ( \beta , \beta )$ 滚动到 $G^{\sigma} ( \Delta \tau , \Delta \tau )$，这是必要的，因为虚时方向上具有周期性边界条件。当一个提议的改变被接受时，它会改变矩阵 $B_{i}^{\sigma}$ 的定义(7.30)中的矩阵 $A_{i}^{\sigma}$。在将 $G^{\sigma}$ 推进到下一个时刻时，我们使用具有新辅助场构型的 $B^{\sigma}$。

作为一种费米子算法，有限温度行列式方法通常存在符号问题（第5.4节和第11.1节）。当 $w_{C}^{\sigma}$ 为负而 $w_{C}^{- \sigma}$ 为正时，符号问题就会出现。从(7.43)中，我们看到对于给定的 x，当 $1 - G_{i i}^{\sigma}$ 为负时，$w_{C}^{\sigma}$ 为负。由于 $1 - G_{i i}^{\sigma} = n_{i \sigma}$，这种情况显然是不物理的。Gubernatis 和 Zhang (1994) 研究了这种情况的来源是否为计算格林函数时的数值精度不足。他们的结论是，虽然精度不足可能导致符号问题，但这并非问题的主要来源。他们指出，在某些情况下，Hubbard-Stratonovich 场的拓扑结构与费米子场的拓扑结构直接相关，并认为符号问题的根源是 Hubbard-Stratonovich 构型中产生的拓扑缺陷（零模）。我们还注意到，在某些情况下，例如低电子填充时，符号问题并不严重。

如果我们能够以无限数值精度执行所有算术运算，那么算法 23 效果很好。但使用有限精度算术时，对格林函数的重复更新和向前逐步推进会累积数值误差。为了克服这种精度损失，如算法 23 所述，我们定期使用 (7.46) 根据当前的辅助场值重新构建格林函数。一种策略是根据经验确定周期 $n_{\mathrm{stab}}$：幸运的是，当数值精度崩溃时，其发生往往相当突然。因此，我们可以尝试几个 $n_{\mathrm{stab}}$ 值，直到精度崩溃，然后选择一个低于该数值的值。

最严重的精度损失发生在对 (7.46) 中矩阵求逆的过程中。随着温度的降低，该矩阵会变得非常病态（ill-conditioned），使用标准矩阵求逆例程对其进行求逆得到的结果可能会完全丧失数值精度。在第7.4节中，我们将描述一些称为矩阵乘积稳定化方法（matrix product stabilization methods）的技术，这些技术能够实现这种求逆。简而言之，这些方法以分解形式累积 B-矩阵的乘积，从而在单位矩阵 I 与 $B^{\prime} {\mathrm{s}}$ 乘积相加时减少信息损失。这种分解形式易于求逆，但使用这些方法所需的编码量比算法 23 中所暗示的要更多。稳定化方法仅会增加少量的总计算成本。在算法中，当我们说“由 $(7.46)$ 得到 $G_{i}^{\sigma}$”时，我们假定这一步是使用矩阵乘积稳定化方法完成的。这些技术的描述虽然具有技术性，但却至关重要。由于其技术性质，我们将讨论推迟到本章后面（第7.4节）。使用这些方法对于精确估计可观测量（尤其是那些需要不同时格林函数的可观测量）也很重要。我们现在讨论可观测量估计的基础知识。

7.2.4 测量#

在无相互作用模型的解析计算中，Wick 定理（Fetter and Walecka, 1971; Negele and Orland, 1988; Fulde, 1991; Enz, 1992）将算符乘积的期望值转化为所有可能的产生湮灭算符成对缩并的乘积之和。这些缩并的期望值正是（无相互作用）格林函数的矩阵元。因此，Wick 定理将测量结果用模拟中的核心量——格林函数——表示出来。

作为在行列式方法中使用 Wick 定理的一个例子，我们考虑反铁磁结构因子：

其中 N 是格点数目，$n_{i \sigma}$ 是自旋为 σ 的电子在格点 i 上的数量，若格点 i 和 j 在同一个子晶格上，则 $( - 1 )^{i - j}$ 为 +1，否则为 -1。蒙特卡洛模拟通过对辅助场变量的 M 个独立样本进行平均，以估计 $S^{z z} ( \pi . \pi )$：

其中量子力学期望值对应于辅助场变量的特定构型 C。由于 Hubbard-Stratonovich 变换将问题在电子自旋上分解，不同自旋对之间的缩并，例如 $\langle c_{\sigma}^{\dagger} c_{- \sigma}^{\vphantom \dagger} \rangle_{C}$，为零。因此，

同自旋项会产生更多非零缩并：

我们注意到，$i = j$ 的情况简化为 $( 1 - G_{i i}^{\sigma} )$，这是因为 $\langle n_{i \sigma} n_{i \sigma} \rangle_{C} =$ $\langle n_{i \sigma} \rangle_{C}$。

在进行测量时，我们可以利用哈密顿量的对称性，一次性收集一个感兴趣的量的许多可能估计。这对于当前算法尤其实用，因为测量计算通常比一次蒙特卡洛扫描所需的计算机时间少得多。例如，要估计一个平移不变系统的动能，我们可以对所有格点位置 i、对称相关位置 $i + \delta$、电子自旋 σ 和虚时 τ 上的 $\langle c_{i , \sigma}^{\dag} ( \tau ) c_{i + \delta , \sigma} ( \tau ) \rangle_{C}$ 测量值进行平均。虽然单次测量不满足哈密顿量的对称性，但对多个构型平均后将满足这些对称性。对尽可能多的等价格点进行平均可以降低最终估计的方差。

模拟中的某些特征可能会人为地破坏哈密顿量的对称性。例如，在哈伯德模型中，交错磁化强度 (staggered magnetization)

其中 $\sigma^{\alpha}$ 是泡利自旋矩阵，这一项具有旋转不变性。由于它将辅助场与局域自旋的 z 分量 $n_{i \uparrow} - n_{i \downarrow}$ 耦合，Hubbard-Stratonovich 变换破缺了这一对称性，使得纵向估计量

比横向估计量（Hirsch, 1987）

的噪声大得多。通常更推荐使用后者。

7.3 Hirsch-Fye 算法#

在本节中，我们将从一般角度重新审视格林函数的更新方程。更一般的视角并不会为格点模型带来更高效的算法，但会引出一个有趣的更新变体，称为 Hirsch-Fye 算法。该算法对于杂质模型（impurity model）的模拟特别有用，它为前述非杂质问题中的更新提供了不同的视角。

通过 (7.45) 式，在给定的虚时 $\tau$ 上为一组 Hubbard-Stratonovich 场值 $x_1, x_2, \ldots, x_i, \ldots, x_N$ 更新格林函数 $G^{\sigma} ( x_{1} , x_{2} , \ldots , x_{i} , \ldots , x_{N} )$，会生成在同一虚时上仅改变一个场后的格林函数 $G^{\sigma} ( x_{1} , x_{2} , \ldots , x_{i}^{\prime} , \ldots , x_{N} )$。该过程精确地解决了在指定时刻改变单个位点上的场的问题。我们可以从无相互作用的格林函数开始，使用相同的更新方程，其中 $\Delta_{11}^{\sigma} ( i , \tau ) = \exp [ \sigma \alpha x_{1} ( \tau ) ] - 1$，生成仅在一个位点激活相互作用后的格林函数 $G^{\sigma} ( x_{1} )$；然后使用这个格林函数来生成在两个位点激活相互作用后的格林函数 $G^{\sigma} ( x_{1} , x_{2} )$，依此类推，直至所有位点上的相互作用都被激活。最终得到的格林函数将继承与 Hubbard-Stratonovich 场相关的虚时。更新所有位点上的 Hubbard-Stratonovich 场需要 $\mathcal{O} ( N^{3} )$ 次浮点运算，即对于 $N$ 个位点，每次更新需要 $\mathcal{O} ( N^{2} )$ 次操作。要对所有虚时切片完成此更新，则需要 $\mathcal{O} ( N_{\tau} N^{3} )$ 次操作。在第 7.2.2 节中，我们还讨论了将等时格林函数从一个虚时推进到另一个虚时，以及从等时格林函数生成不等时格林函数的程序，其复杂度也是 $\mathcal{O} ( N^{3} )$。

在一般情况下，对于离散 Hubbard-Stratonovich 场的配分函数的路径积分具有如下形式（Blankenbecler 等人，1981; Negele 和 Orland，1988）

其中 $X ( C ) = X^{\sigma} ( C ) X^{- \sigma} ( C )$ 是描述无相互作用电子动力学及其与辅助场耦合的矩阵表示。它的规模是 $N_{\tau} \times N$，但非常稀疏：

也就是说，它是一个 $N_{\tau} \mathrm{~} \times N_{\tau}$ 的分块矩阵，其中每个块都是一个 $N \times N$ 的矩阵。所有位点和所有时刻的格林函数 $\bar{G}^{\sigma}$ 正是 $X^{\sigma} \left( x \right)$ 的逆（Negele 和 Orland，1988）。虽然这些陈述并不直观，但通过一些简单的代数推导可以使其易于理解。例如，从 (7.47) 式出发，适当地选择分块行和分块列的各种置换，并利用以下矩阵分块求逆公式

其中

我们可以生成格林函数的全部等时（(7.24)式）和不等时（见下文(7.56)式）表达式。因此，利用 $X^{\sigma} \left( x \right)$ 的结构和稀疏性，可将（7.48）式求逆所需的 $( N_{\tau} N )^{3}$ 量级代价降低为之前推导算法中更高效率的 $N_{\tau} N^{3}$ 量级。

与(7.49)相关，我们还需注意一个有用的关系式

利用该关系式，容易证明 $X ( C )$ 的行列式与(7.19)、(7.20)和(7.21)式的等价性。

在 $N_{\tau} \times N$ 维的更大空间中，总格林函数的更新过程是怎样的？为推导此过程，我们首先从 $M^{\sigma}$ 矩阵的定义出发，并利用B矩阵的定义(7.30)式来定义

其中我们将化学势因子吸收到了跳跃矩阵 $T$ 中，并且由(7.31)式得到 $N \times N$ 对角矩阵

该定义将依赖于辅助场的项与不依赖辅助场的项分离开来。接下来，我们定义 $N_{\tau} \times N_{\tau}$ 分块对角矩阵

然后引入另一个定义：矩阵 $\hat{G}^{\sigma} \equiv \bar{G}^{\sigma} e^{V^{\sigma}}$ ，即

由最后一个表达式和算子恒等式

可以简明地证明：如果将 $\exp ( - V^{\sigma} )$ 更新为 $\exp ( - V^{\prime \sigma} )$ ，则更新后的格林函数 $\hat{G}^{\prime \sigma}$ 满足戴森方程（Dyson equation）

最后，将 $\hat{G}^{\sigma}$ 的定义代入该方程，得到感兴趣的格林函数 $\bar{G}^{\sigma}$ 的表达式

该方程推广了(7.39)式。新的戴森方程并非更新给定时刻的所有格点，而是更新所有时刻的所有格点，这是一个 $\mathcal{O} ( ( N_{\tau} N )^{3} )$ 量级的操作。

(7.51)式的通用性允许我们定义一个过程来更新给定格点在所有时刻的取值。这一能力对于杂质模型（impurity model）尤为有用，即一个格点在物理上与其他格点不同的问题；例如 (7.4)式。由(7.51)式，我们很容易确信：如果相互作用局域在某个特定格点（比如0）上，则通用更新方程可简化为仅针对该格点的格林函数在所有时刻的 $N_{\tau} \times N_{\tau}$ 矩阵方程：

按照与之前相似的代数运算，这个更新方程可化为与(7.45)相同的形式。这里的更新每次针对一个 $\tau$ 值进行。由于矩阵是 $N_{\tau} \times N_{\tau}$ 的，更新所有时刻所需的计算时间按 $N_{\tau}^{3}$ 标度增长。

这种杂质更新程序就是 Hirsch-Fye 算法。它最初是为杂质问题而发展的(Hirsch 和 Fye, 1986)，随后被推广用于晶格问题。该算法很容易推广到包含 $N_{c}$ 个杂质位点的团簇，相应算法的标度为 $( N_{c} N_{\tau} )^{3}$ 。当团簇达到一定大小时，“固定时刻的所有位点”流程变得更加高效。在连续时间杂质求解器（第8章）发展之前，Hirsch-Fye 算法被广泛用于解决单位点和团簇杂质问题。

关于 Hirsch-Fye 算法，有几点值得注意：尽管执行蒙特卡洛更新只需要杂质格林函数，但一旦有了这个格林函数，通用的更新方程告诉我们，如何利用杂质格林函数和没有杂质时的格林函数，来计算杂质与晶格中其他位点之间的空间矩阵元格林函数(Gubernatis et al., 1987)。这种能力，在空间维度上相当于找出所有位点在不同虚时时刻的格林函数。

在本节开始时，我们注意到，对于一个包含相互作用位点的晶格，如果已知无相互作用格林函数，就可以逐个激活相互作用。在杂质问题中，类似甚至更好的做法是可行的。为了构建无相互作用杂质格林函数，我们需要杂质位点上没有相互作用时的格林函数。一旦我们有了杂质格林函数，就可以激活 Hubbard-Stratonovich 场。之后，我们就不再需要其他的（空间非局域）格林函数。对于没有杂质时的格林函数，我们可以使用无限大晶格的格林函数（如果可用），或者为非常大的晶格计算它。简而言之，我们能够模拟杂质问题而无需担心有限尺寸效应。在下一章中，我们将讨论消除由虚时离散化引起的误差的方法。

7.4 矩阵乘积稳定化#

本节描述的矩阵方法是行列式方法的重要组成部分，但它们在本质上与多电子物理关系不大。它们主要用于计算等时格林函数。在这一计算中，B矩阵的乘积被累积并加到单位矩阵上，以便在对矩阵和求逆后仍保持足够的精度。⁸ 问题在于需要求逆的矩阵通常是病态的。

矩阵 M 的条件数 $\kappa$ 定义为 $\kappa = \| M \| \| M^{- 1} \|$ ，其中 $\lVert M \rVert$ 是某种矩阵范数(Golub 和 Loan, 1989; Meyer, 2000)。矩阵特征值绝对值的最大值就是其中一种范数。条件数衡量一个矩阵接近奇异的程度，即是否接近有一个或多个零特征值。如果 $\kappa = 10^{d}$ ，那么当矩阵被求逆时（例如，使用部分主元高斯消去法），$d$ 大致表示损失的十进制位数(Golub 和 Loan, 1989)。通过考虑 $U = 0$ 的情况并将hopping 限制为最近邻，我们可以感受到当前应用中的困难。对于正方晶格，最近邻hopping 矩阵 T 的特征值广为人知：$\lambda ( k_{x} , k_{y} ) = - 2 t [ \cos ( k_{x} ) + \cos ( k_{y} ) ]$ $( 0 \leq k_{x , y} \leq 2 \pi )$，因此 $- 4 t \leq \lambda ( k_{x} , k_{y} ) \leq 4 t$ 。当 $\mu = 0$ 时，

是一个条件数大约为$\exp (8 \beta t)$的矩阵。当温度约为$t / 5$时（这是电子带宽8t的一小部分），该条件数意味着损失17位十进制精度，这一数值恰好超出了双精度算术所能表示的范围。

在低温下对多电子系统进行热力学研究的基本困难在于，玻尔兹曼分布$e^{- \beta ( H_{1} - \mu N )}$赋予相对于化学势的最低能态以指数级大的权重。然而，泡利不相容原理阻止这些态被宏观占据，并迫使它们在被占据到某个“费米能”之前只能被单占据。对角化后的$B ( \beta , 0 )$的矩阵元恰好是1减去费米-狄拉克分布，

其中$\lambda$为每个特征值。正确地将1加到$\exp [ - \beta ( \lambda - \mu ) ]$上是确保电子具有正确费米-狄拉克分布的第一步。求$I + B ( \beta , 0 )$逆的过程，实际上就是将单位矩阵$I$正确地加到玻尔兹曼分布$\exp [ - \beta ( H - \mu N ) ]$上，以便矩阵和的求逆运算能够保持电子的费米-狄拉克特性。对于任意辅助场构型，单粒子传播子$B^{\sigma}$会衰减“高能”尺度，同时放大“低能”尺度。因此，提取化学势附近“中间能量”的信息是一项非常困难的数值任务，但又必须精确完成，因为正是这些能量态控制了有限温度的物理行为。如果没有将要描述的处理过程，模拟中费米子的特性就会丢失。

在求逆一个包含数值尺度差异极大的矩阵时，数值局限性完全不同于求逆单个数值。对于单个数值，双精度算术通常允许尺度在$10^{- 307}$到$10^{+ 309}$之间变化。针对求逆病态矩阵所设计的策略（Loh et al., 1989; Loh and Gubernatis, 1992; Loh et al., 2005）是对矩阵乘法进行分层，用$B$矩阵的一种表示将这些数值尺度放入对角矩阵中。我们可以精确地求逆这些矩阵，因为这只需对单个数值求逆。该过程还以避免混合不同尺度的方式来累积矩阵乘法。只有最后的矩阵乘法步骤会混合尺度并导致精度损失。剩余的精度足以保证测量误差至少与统计误差一样小。

所采用的方法是将矩阵分解为$UDV$的形式，其中对角矩阵$D$包含不同的尺度，而$U$和$V$的条件数相当良好且易于求逆。下面说明了因式分解形式的特征以及我们试图避免的问题：

这里我们使用字符 x 的尺寸来表示不同的数值量级。矩阵 U 和 V 具有单位量级。上述方程第一行中的三个矩阵表示我们在分解中想要实现的目标。第二行说明了乘法发生时的情形——最大的数值量级占据主导地位，而小量级的信息则丢失了。

另一方面，如果矩阵 M 的数值量级按列分层排列，我们可以稳定地对其进行分解，从而保留其小量级信息：

从左侧对 M 乘以 $U^{- 1}$ 仅会组合给定列中的元素，因此不会造成信息丢失。从右侧乘以 $V^{- 1}$ 会组合不同量级的列，但只要在较大量级的列被添加到较小小量级的列之前先将其进行缩小，就不会覆盖任何小量级信息。使用 $M V^{- 1} = ( V^{- T} M^{T} )^{T}$ 可以减少这种混合。

为了以稳定的方式计算 (7.29) 中多个矩阵的乘积，我们在整个计算过程中将多个量级分离开。作为说明，我们首先假设有 $B ( \tau , 0 ) = U D V$ 。如果我们希望得到 $\tau + \tau_{0}$ 处的传播子（其中 $\tau_{0}$ 是几个虚时步长的大小）⁹，我们可以使用标准矩阵乘法方法来扩展传播子。我们写出

从而得到 $B^{\sigma} ( \tau + \tau_{0} , 0 )$ 的 UDV 分解。我们所做的是将分层矩阵 $B^{\sigma} ( \tau + \tau_{0} ) U D$ 分解为 $U^{\prime} D^{\prime} V^{\prime}$ 。V 矩阵必须具有足够好的条件数，以便我们能够以稳定的方式将多个它们相乘。根据需要重复此过程，即可通过从左侧相乘来计算 B 矩阵 (7.21) 的长乘积链。

在矩阵乘积以分解形式组装完成后，我们通过以下步骤计算逆 $[ I + U D V ]^{- 1}$ ：

我们注意到，在形成和 $U^{- 1} V^{- 1} + D$ 之后，我们将其分解为 $U^{\prime} D^{\prime} V^{\prime}$ ，并分别对各个因子执行求逆运算。

上述方法的一种替代方案是执行两个部分乘积的分解

在第二个方程中，我们反转了 UDV 因子，因为我们打算通过从左侧相乘来构建部分乘积。现在我们有

示意性地，需要求逆的部分为

由于我们将对角矩阵保留在外侧，只有在最后一步，当我们将不同尺度的元素相加时，才会出现小尺度截断。我们通过先进行另一次 UDV 分解来对结果矩阵求逆。

如何创建 UDV 分解？有几种标准的数值矩阵方法可以使用，包括改进的 Gram-Schmidt、QR 和奇异值分解方法（Golub and Loan，1989；Meyer，2000）。$U D V^{T}$ 分解是奇异值分解方法的标准结果。标准的改进 Gram-Schmidt 和 QR 算法很容易重写以产生这一结果。改进的 Gram-Schmidt 方法虽然在数值稳定性上稍逊，但速度明显更快，通常也足够好。如何对各种矩阵求逆？对角矩阵的求逆非常简单。实际上，U 和 V 也是如此。改进的 Gram-Schmidt 方法产生一个列正交的 U 和一个单位上三角的 V。前者我们通过转置求逆，后者则通过对线性方程组 $V V^{- 1} = I$ 进行回代求解未知的 $V^{- 1}$ 来求逆，该逆矩阵也是一个三角矩阵。

Gram-Schmidt 正交化方法如下：给定一组向量 $a_{1} , a_{2} , \dots , a_{M}$（它们可能是一个矩阵的列），求出一组标准正交的向量 $\nu_{1} , \nu_{2} , \ldots , \nu_{M}$，它们张成相同的空间。这个思路非常简单：取 $a_{1}$，除以其长度 $\| a_{1} \| = \sqrt{a_{1}^{T} a_{1}}$ 得到 $\nu_{1}$。然后减去 $a_{2}$ 在 $\nu_{1}$ 方向上的分量。余量正交于 $\nu_{1}$。接着，将余量除以其长度得到 $\nu_{2}$。重复这个过程，减去 $a_{3}$ 在 $\nu_{1}$ 和 $\nu_{2}$ 方向上的分量。将余量除以其长度得到 $\nu_{3}$。继续这个过程，我们就得到了所需的标准正交集。不幸的是，如果任一 $a_{i}$ 接近线性相关，这个数值过程将不稳定。幸运的是，解决方法很简单：重新组织计算，使得一旦得到 $\nu_{i}$，就不再使用 $a_{i}$。算法如算法24所示。用一种更接近计算机语言的形式，算法25用伪代码描述了正交化和矩阵分解。

算法24 改进的 Gram-Schmidt 方法。
输入：M 个向量 a<sub>1</sub>, a<sub>2</sub>, . . . , a<sub>M</sub>。
for i = 1 to M do
$a_{i} a_{i} / a_{i} $；
for $j = i + 1$ to M do
$a_{j} a_{j} - ( a_{i}^{T} a_{j} ) a_{i}$；
end for
end for
return M 个标准正交向量 $a_{1} , a_{2} , \dots , a_{M} .$

算法25 改进的 Gram-Schmidt 方法：UDV 分解。 输入：N × N 数组 U。 for k = 1 to N do d(i) = 0； for i = 1 to N do $d ( k ) \gets d ( k ) + U ( i , k ) U ( i , k )$ end for $d ( k ) \sqrt{d ( k )}$； for i = 1 to N do $U ( i , k ) \gets U ( i , k ) / d ( k )$ end for V(k, k) = 1； for j = k + 1 to N do V(k, j) = 0； for i = 1 to N do $V ( k , j ) \gets V ( k , j ) + U ( i , k ) U ( i , j )$ end for for i = 1 to N do $U ( i , j ) \gets U ( i , j ) - V ( k , j ) U ( i , k )$ end for $V ( k , j ) V ( k , j ) / d ( k )$ end for end for return N × N 标准正交矩阵 U，N × N 对角矩阵 D 的 N 个对角元 d，以及 N × N 单位上三角矩阵 V。