第4章 · 精确对角化方法#

通过精确对角化哈密顿量，可以完整了解一个量子自旋系统——有了本征态，任何静态或动态量都能计算。原则上，对于一个有限量子系统，通过构建哈密顿量矩阵并对其进行数值对角化，可以精确计算出所有本征态。然而在实践中，这种精确对角化研究仅限于相当小的晶格，即几十个自旋，因为基矢大小随自旋数量呈指数增长（对于 $S = 1/2$ 的情况为 $2^{N}$ 个态）。因此，在得出关于热力学极限的结论时必须极为谨慎，这甚至可能不可行（如果可用晶格太小而无法容纳无限大系统的物理特性）。尽管如此，从精确对角化研究中获得的见解仍然非常有用，它们本身有价值，同时也作为对其他计算方法的补充。小晶格的精确结果对于检验例如量子蒙特卡洛程序的正确性也是不可或缺的。此外，精确对角化方法为学习量子力学的许多重要方面提供了具体途径，特别是多体态的对称性质。

分块对角化。给定一个哈密顿量 H，精确对角化计算的第一步是选择一个基矢，在该基矢下表达哈密顿量及其他感兴趣的算符。一个自旋-1/2 系统的工作基矢通常由单自旋态 $\uparrow_i$ 和 $\downarrow_i$ 组成，其中 $i = 1, \ldots, N$（量子化方向通常取为 z 方向）。然而，由于希尔伯特空间随系统中自旋数量以 $2^{N}$ 增长，应尽可能利用对称性首先将哈密顿量简化为分块对角形式，如图 25 所示。在这种方案中，自旋态借助适用的对称操作进行组合和排序。这些块对应于具有与对称性相关的不同守恒量子数的态，例如由晶格平移对称性导致的晶体动量守恒，或者总自旋 z 分量的守恒（磁化强度）。这些块可以相互独立地对角化，从而大幅降低计算成本。除了减少计算量之外，直接获取量子数对于对激发进行分类也非常有用。

图 25. 分块对角化的示意图。在原始基矢下，哈密顿量没有明显结构（左）。通过构造由守恒量子数标记的态，矩阵分解成若干块（阴影方块之外的所有矩阵元为零），这些块可以相互独立地对角化（中）。应用另一个对称性（守恒律），这些块可以进一步分解成更小的块（右），由两个不同的量子数标记，等等。

有些对称性相对容易利用（例如磁化强度的守恒），而另一些则需要更多工作，并导致更复杂的基矢态（例如动量态）。一些原则上可用的对称性通常不被实施，因为计算中的实际复杂性可能超过其带来的好处。总自旋守恒就是一个例子。

在精确对角化中使用对称性可以通过群论语言来讨论 [143]。然而，这种形式化并非必需（而且常常令人困惑），这里采用更实用的方法，不涉及群论术语。在处理复杂晶格时，群论实际上非常有用，但其形式化的力量或许可以通过先对简单晶格采用较不形式化的方法，深入理解对称操作和分块对角化之后，才能更好地体会。这里我们考虑一维链和二维简单正方晶格。

本节概览#

我们的论述首先将以$S=1/2$海森堡链为具体例子，在4.1节展开。我们首先在4.1.1节讨论其对称性，并介绍一种使用整数的比特位来表示基态的计算机表示方法。在4.1.2节中，我们在伪代码中使用这种表示法来构建不考虑任何对称性的完整哈密顿量，然后利用磁化强度守恒进行分块对角化。4.1.3节讨论动量态，4.1.4节讨论宇称（反射对称性）的使用，4.1.5节讨论自旋反转对称性。4.1.6节通过一些结果展示了完全对角化以获得有限温度热力学性质的过程。为了计算更大尺寸链的少量本征态，将在4.2节中介绍Lanczos方法。4.3节将展示该技术在几个一维系统中的效用。4.3.1节讨论海森堡链的基态和低能激发。4.3.2节研究在存在阻挫性次近邻相互作用时发生的二聚化相变，4.3.3节则考虑一个具有长程相互作用的扩展模型。4.4节简要描述了二维系统的对角化。文中将讨论动量态和其他正方晶格的对称性，并讨论与该系统奈尔基态相关的计算结果。

4.1 海森堡链的对角化#

我们将研究$S = 1 / 2$反铁磁海森堡链，其哈密顿量为

其中，为了下面将明了的理由，这里将自旋标记为$i = 0 , \ldots , N - 1$比较实用。当我们考虑动量态时，将假设周期性边界条件$\mathbf{S}_{N} = \mathbf{S}_{0}$，但在此之前边界条件是任意的。

4.1.1 态与对称性的表示#

晶格变换。我们使用标准记号$| S_{0}^{z} , \ldots , S_{N - 1}^{z} \rangle$表示基态，其中从左到右的$S_{i}^{z}$始终对应于编号为$0 , 1 , 2 , \dots$的晶格位点上的自旋态，无论下标i的顺序如何。因此，如果我们写出一个态$| S_{1}^{z} , S_{0}^{z} \rangle$，除非$S_{0}^{z} = S_{1}^{z}$，否则它不同于$| S_{0}^{z} , S_{1}^{z} \rangle$。前者可以指通过置换算符P交换两个自旋后从后者得到的态，该算符可以在操作上定义为影响位点指标；$P | S_{0}^{z} , S_{1}^{z} \rangle = | \bar{S}_{1}^{z} , S_{0}^{z} \rangle$，例如，$P | \uparrow \downarrow \rangle = | \downarrow \uparrow \rangle$。将此推广到N自旋晶格变换、反射、平移和旋转（在二维和三维中），它们都是根据自旋指标的置换来定义的。例如，对于一个周期链，我们将平移算符定义为将自旋循环地向“右”移动一步；

这对应于在ket中的每个位置i处将自旋指标减一（模系统尺寸N）：$S_{i}^{z} \leftarrow S_{i - 1}^{z}$。当用向上和向下自旋（分别用$\uparrow$和$\downarrow$表示）写出特定态时，例如$\mid \uparrow \downarrow \uparrow \downarrow \cdots \rangle$，通常不需要指标（并且，为了记号的紧凑，箭头之间也不需要用逗号隔开）。

具有周期性边界条件的海森堡哈密顿量（106）在平移下是不变的，即，它与$T$对易：$[ H , T ] = 0$。因此，我们可以构造动量态$| \Psi ( k ) \rangle$，根据定义，它们是平移算符的本征态，

这里允许的动量为 $k = 2 n \pi / N$，其中 $n = 0 , \ldots , N - 1$，这是由 $T^{N} = 1$ 得出的。具有不同 $k$ 的态形成哈密顿量各自可单独对角化的区块。实际中如何构造动量态以及如何在计算机程序中使用它们，将在第4.1.3节讨论。

海森堡哈密顿量也与反射（宇称）算符对易，我们以推广已考虑过的二自旋置换的方式，根据自旋指标变换 $i N - 1 - i$ 来定义该算符：

对于 $P$ 的本征态，$T | \Psi ( p ) \rangle = p | \Psi ( p ) \rangle$，其中 $p = \pm 1$，因为 $P^{2} = 1$。我们将使用 $T$ 和 $P$ 进行分块对角化，但正如我们在第4.1.4节详细讨论的，它们不能总是同时使用，因为 $[ T , P ] = 0$ 仅在希尔伯特空间的一个子空间中成立。对于具有开边界的系统，$T$ 未被定义，但可以使用 $P$。

自旋量子数。由于哈密顿量是自旋旋转不变的，其本征态也必须是总自旋平方 $\mathbf{S}^{2}$ 的本征态，其中

对于一个本征态，我们有 $\mathbf{S}^{2} | \psi ( S ) \rangle = \mathbf{S} \cdot \mathbf{S} | \psi ( S ) \rangle = S ( S + 1 ) | \psi ( S ) \rangle$。这里存在潜在的混淆，因为相同的符号 $S$ 既用于表示单个自旋的自旋大小（即，$S_{i} = S$），也用于表示多体态的总自旋。上下文应始终能说明含义（而且我们在这里只考虑 $S_{i} = 1 / 2$）。

由于总自旋守恒，我们知道态可以形成由量子数 $( S , m_{z} )$ 标记的区块，其中 $m_{z} \in \{- S , - S + 1 , . . . , S \}$ 是沿量化轴方向的总磁化强度，

如果我们使用动量态，每个 $k$ 区块当然也会分裂成更小的区块 $( k , S , m_{z} )$，因为 $[ T , \mathbf{S} ] = 0$。使用 $m_{z}$ 进行分块对角化很容易，但实现总 $S$ 的守恒通常很繁琐（除非自旋数量很少），因此很少使用（尽管价键基中的 $S = 0$ 态有时被使用[144]）。我们将结合晶格对称性来使用 $m_{z}$ 守恒。注意，$m_{z}$ 守恒比总自旋 $S$ 守恒更普遍——即使我们在哈密顿量中给不同的项赋予不同的前因子，例如，给哈密顿量(106)中的伊辛（对角）项，$m_{z}$ 仍然守恒，尽管总自旋 $S$ 不守恒。因此，本节讨论的所有技术也可以直接应用于此类各向异性模型。

对于特殊情况（且最重要的情形）$m_{z} = 0$（对于偶数 $N$），我们可以利用自旋空间中所有可能旋转的一个离散子集进行分块对角化；即自旋反转对称性，即关于翻转所有自旋的不变性。这可以通过一个我们称为 $Z$ 的算符来形式化定义：

对于该算符，我们同样有$Z^{2} = 1$，其特征值为$z = \pm 1$。由于$Z$与P和$T$都对易，它可以与这些算符一起使用来进一步将H分块对角化，我们将在第4.1.5节中进行。

图26. 顶行显示一个$N = 8$的自旋态s及其作为整数s的前八位的表示。注意，我们将自旋$i = 0 , \ldots , N - 1$从左到右标记，而比特的常规标记是从右到左的二进制数。要翻转两个自旋i和j，可以使用一个掩码f（中间行）进行按位异或（XOR）运算。f的位i和j被设置为1，所有其他位为$0$（即$f = 2^{i + j}$）。底行显示了XOR运算的结果（Fortran 90中的$\mathbf{ieor} ( s , f )$）。

总自旋算符$\mathbf{S}^{2}$可以写成类似于所有自旋之间具有相等相互作用的海森堡模型的形式：

构造这个算符的矩阵形式（用于计算量子数 $S$）与构造哈密顿量矩阵非常相似。

自旋态的比特表示。$S = 1 / 2$ 模型很特殊，因为自旋可以直接用整数的比特值 0 和 1 在计算机中表示。我们将利用这一特性。比特通常从 0 开始编号，这也是我们如此给自旋编号的原因。为了引用整数 s 的比特 $i = 0 , \ldots , 31$（或长整数的 $i = 0 , \ldots , 63$），我们将使用记号 $s [ i ]$。因此，对于 N 个自旋，一个基矢态 $| S_{0}^{z} , \ldots , S_{N - 1}^{z} \rangle$ 在计算机中用整数 s 表示，其中 $s [ i ] = \stackrel{\smile} {S_{i}^{z}} + 1 / 2 \stackrel{\cdot \cdot} {\mathrm{for}} i = 0 , \dots , N - 1$，且 $s [ i ] = 0$ 对于 $i > N - 1$。

大多数计算机语言都有检查和操作比特的函数。海森堡哈密顿量 (106) 的非对角项会翻转两个自旋。在伪代码中，我们将使用一个比特函数 $\mathbf{fiip} ( s , i , j )$ 来实现此操作，该函数翻转表示状态的整数 s 的比特 i 和 j（即 0 和 1 互换）。其具体实现方式取决于所用的语言。一种可能性是使用带掩码的按位异或操作，如图 26 所示。在 Fortran 90 中，此操作可以实现为 ieor(s, 2^{i+j})。稍后，我们还需要实现状态的各种对称性变换的函数：平移、反射和自旋反演。

使用标准的 4 字节整数，用单个整数的比特表示法最多适用于 $N = 32$。使用长整数，可以将方案扩展到 $N = 64$。这远超出了实际中精确对角化技术能处理的最大尺寸，除非磁化扇区 $2 m_{z} = n_{\uparrow} - n_{\downarrow}$ 足够大，使得块尺寸 $N ! / ( n_{\uparrow} ! n_{\downarrow} ! )$ 变得可处理。然而，磁化 $m_{z} = 0$ 和其他低 $m_{z}$ 扇区通常是主要关注点。

在讨论构建基矢集和哈密顿量矩阵的算法时，我们将从最简单的方法开始，即不使用任何对称性，然后实现磁化守恒。在第 4.1.3-4.1.5 节中，我们将纳入更多对称性。哈密顿量矩阵的实际对角化以及利用其本征态计算物理可观测量将推迟到第 4.1.6 节讨论。

4.1.2 哈密顿量的计算机生成#

如果不利用任何守恒律，哈密顿量就是一个 $2^{N} \times 2^{N}$ 的矩阵。我们只需使用数字 $a = 0 , 1 , . . . , 2^{N} - 1$ 来标记基矢态。这些整数的比特值直接对应于自旋态。确定哈密顿量矩阵的对角贡献 $\langle a | S_{i}^{z} S_{1 + 1}^{z} | a \rangle = \pm 1 / 4$ 只需检查对应的比特对 $a [ i ] , a [ i + 1 ]$（相同或不同），而非对角算符 $( S_{i}^{+} S_{i + 1}^{-} + \bar{S_{i}^{-}} S_{i + 1}^{+} ) \bar{/ 2}$ 作用于一个满足 $a [ i ] \neq a [ i + 1 ]$ 的状态 $| a \rangle$ 上，会产生状态 $\left| b \right.$，其中两个比特被翻转。相应的矩阵元为 $\langle b | H | a \rangle = 1 / 2$。对于 $a [ i ] = a [ i + 1 ]$ 的状态，当然不存在非对角矩阵元。以下伪代码生成完整的哈密顿量矩阵 H：

do a = 0,2$\mathbb{N}$ 1
do i = 0,N 1
j = mod(i + 1, N) {2}
if (a[i] = a[ j]) then
H(a, a) = H(a, a) + 1/4
else
H(a, a) = H(a, a) - 1/4
b = flip(a, i, j); H(a, b) = 1/2
endif
enddo
enddo

这里 $j = \mathbf{mod} ( i {+} 1 , N )$ 是 $i$ 的“右”最近邻，模函数处理了周期性边界条件。在 Fortran 90 中，可以用布尔函数 btest $( a , i )$ 来检验 $a [ i ] = a [ j ]$，以此检查 $a$ 的第 $i$ 位。哈密顿量中的每个键算符对应于一个非对角矩阵元，而对角元则有 $N$ 个不同的贡献。其他感兴趣的算符对应的矩阵，当然可以通过类似的方式构造，即根据涉及的 $S_{i}^{z}$ 和 $S_{i}^{\pm}$ 算符的任意组合来检查和翻转比特位。

使用固定磁化板块。 进一步提升复杂度，我们接下来实现磁化守恒。我们希望构造作用在具有给定 $m_{z} = ( n_{\uparrow} - n_{\downarrow} ) / 2$ 的所有态上的板块哈密顿量。这样的态共有 $M = N ! / ( n_{\uparrow} ! n_{\downarrow} ! )$ 个，我们需要一个列表来记录它们。这个列表中的顺序将用作板块态的标签，即，我们将列表中的第 $a$ 个态记为 $| a \rangle$。然后我们还需要一个整数列表 $\left\{s_{a} \right\}$，其比特位代表第 $a$ 个态的自旋构型。之后，我们必须搜索列表 $s_{a}$ 以找到给定态整数 $s$ 的位置标签 $a$，因此将列表排序是实用的；$s_{a} < s_{a + 1}$。当明确提及自旋时，我们有时会用符号 $\vert s_{a} \rangle$ 代替 $| a \rangle$，即 $| s_{a} \rangle = | s_{a} [ 0 ] - 1 / 2 , . . . , s_{a} [ N \dot{-} 1 ] - 1 / 2 \rangle$。根据上下文可以明确内部标签是指 M 个态的列表中的位置，还是指包含自旋信息的整数。

为了构建态列表，我们遍历整数 $\mathbf{\bar{\Sigma}}_{s} = 0 , . . . , 2^{\bar{N}} - 1$，并检查置位比特数（态中自旋向上的数目 $n_{\uparrow}$）是否对应于目标扇区；$n_{\uparrow} = m_{z} + N / 2$。在初始化态计数器 $a = 0$ 之后，我们可以使用以下伪代码生成具有给定磁化的所有态：

现在我们有了一个大小为 M 的基，以整数 $s_{a} , a = 1 , \ldots , M$ 的形式存储。

为了构建哈密顿量，我们遍历标签 $a = 1 , . . . , M$，并像之前一样使用位运算作用于对应的态整数 $s_{a}$。当对 $ {\left| {s_{a}} \right.}$ 的非对角运算导致另一个态 $\left| {s^{*}} \right. = \left| {s_{b}} \right.$ 时，我们必须找到态列表中整数 $s^{*}$ 的位置 $b$。由于这个列表是有序的，我们可以通过二分搜索来完成，平均需要 $\sim \log_{2} ( M )$ 步。这种搜索通过一系列括号化进行，每一步我们都可以通过检查范围中点处的态来将 $b$ 的可能范围减半，初始括号为 $b_{\operatorname*{min}} = 1$ 和 $b_{\operatorname*{max}} = M$。以下子程序用于找到态整数 $s^{*}$ 的位置 $b$；

subroutine findstate $\cdot ( s^{*} , b )$ $b_{\operatorname*{min}} = 1 ; b_{\operatorname*{max}} = M$ {4} do $b = b_{\operatorname*{min}} + ( b_{\operatorname*{max}} - b_{\operatorname*{min}} ) / 2$ if $( s^{*} < s_{b} )$ then $b_{\operatorname*{max}} = b - 1$ elseif $( s^{*} > s_{b} )$ then $b_{\operatorname*{min}} = b + 1$ else exit endif enddo

这里整数除以 2 应视为标准方式，即对于奇数 $i$，$i / 2$ 等于 $( i - 1 ) / 2$。当基态 $s_{b}$ 等于目标态 $s^{*}$ 时，循环退出。也可以使用更高效的搜索过程，例如所谓的哈希表（hash-table）[145]。然而，二分法简单得多，并且在大多数情况下足够快。

使用子程序 $\{4 \}$，可以通过对代码 2 进行以下修改，构造出哈密顿量中源自对态 $| a \rangle$ 中自旋对 $( i , j )$ 进行操作的部分；

如果只是想对非常小的晶格快速得到结果，或者系统是非周期的（如开链或具有随机耦合的系统，此时动量不守恒），那么以这种形式构造哈密顿量并对其直接对角化（依次对所有所需的$m_{\bar{z}}$分块进行对角化处理）可能就足够了。但对于平移不变（周期性）系统的深入研究，值得实现额外的对称性来进一步将固定的$m_{z}$分块块对角化。

4.1.3 动量态#

接下来，我们构造由式(107)定义的平移算符$T$的本征态。平移N步后，自旋将回到其原始状态。因此，$T^{N} = 1$，这意味着本征值为$\mathrm{~ e}^{i k}$，其中一组N个不等价动量可以选取为，

这里晶格常数取为1。动量态可以使用一个参考态$| a \rangle$（z分量的一个基矢态）及其所有平移态来构造；

可以很容易地验证（通过对求和指标进行周期性边界条件允许的移位），对此态作用平移算符(107)得到$T | a ( k ) \rangle = {\mathrm{e}}^{i k} | a ( k ) \rangle$，这正是动量态的定义。

为了构造给定k（通常也给定$m^{z}$）的动量基，需要找到一组代表态，使其构成一组完备且可归一化的正交态。显然，要使两个态$| a ( k ) \rangle$和$| b ( k ) \rangle$正交，相应的代表态必须满足：对所有$r$，$T^{r} | a \rangle \neq | b \rangle$。因此，在平移态集合$| a ( r ) \rangle = T^{r} | a \rangle$（$r = 0 , \ldots , N - 1$）的所有态中，只能选一个作为代表态。由于标记采用位表示，实际中通常选择整数$a ( r )$最小的那个作为代表态（通过执行所有平移来确定）。

归一化与排除的代表态。 如果所有平移态${\cal T}^{r} | a \rangle$都是不同的，则式(115)中的归一化常数就是$N_{a} = N$。然而，有些参考态的周期性小于$N$，这会影响归一化。一个态的周期性定义为满足下式的最小整数$R_{a}$：

如果$R < N$，那么式(115)的求和中会出现同一态的多个副本，需要相应修改归一化常数。此时也可以将式(115)中的求和限制为$r = 0 , \ldots , T_{a} - 1$，但在对态进行形式化操作时，保留全部N项（无论$R_{a}$如何）更为实用。

动量基的一个重要方面是，代表态的周期性必须与动量相容，才能使式(115)成为有效的态。这种相容性与可归一化性有关。在式(115)的求和中，与代表态$| a \rangle$相关的相位因子的和为：

此时的归一化常数为：

因此，如果$F ( k , R_{a} ) = 0$，则无法使用参考态$| a \rangle$定义具有动量$k$的态。因此，对于给定的$| a \rangle$，允许的动量是那些$k R_{a}$为$2 \pi$倍数的动量，即：

对于允许的动量，式(115)中的归一化常数为：

本章新增术语：

当参考态不为其自身的任何非平凡平移时，有 $R_{a} = N$ 和 $N_{a} = N$。再次注意，给定的参考态只能出现在一个基矢态中。于是动量态显然是正交归一的：$\langle \bar{b} ( k^{\prime} ) | a ( k ) \rangle = \delta_{a b} \delta_{k k^{\prime}}$。还要再次注意，不与给定动量相容的参考态不会出现在该组块态中。构建动量基的一个重要方面是检查潜在代表态与该动量的相容性。

哈密顿矩阵。接下来，我们构建动量基中的哈密顿矩阵。这里考虑的周期性海森堡哈密顿量具有平移不变性，由 $N$ 个键算符 $\mathbf{S}_{i} \cdot \mathbf{S}_{i + 1}$ 组成。将所有对角项合并在一起，并分别考虑非对角项是比较方便的。为了简化后续的形式，我们相应地定义算符：

因此 $\begin{array} {r} {H = J {\sum_{j = 0}^{N}} {H_{j}}} \end{array}$。我们现在设 $J = 1$。我们需要求出当 $H$ 作用于动量态 (115) 时所得到的状态。由于 $[ H , T ] = 0$，我们可以写为

并且我们只需要将哈密顿算符 $H_{j}$ 作用于参考态。每次操作我们都会得到一个不同的态，或者在对角 ($j = 0$) 的情况下，得到相同的态。无论哪种情况，我们都可以写为 $H_{j} | a \rangle = h_{j} ( a ) | b_{j}^{\prime} \rangle$，其中 $h_{j} ( a)$ 是由 (121) 或 (122) 得出的矩阵元，并且为了简化符号，我们不在 $\left| b_{j}^{\prime} \right.$ 中包含任何表明它也依赖于 $| a \rangle$ 的显式指示符。$\left| b_{j}^{\prime} \right.$ 上的撇号用于表明这个新态不一定是用于定义基矢的参考态之一，因此，不应直接基于它来写动量态。假设 $\left| b_{j}^{\prime} \right.$ 与动量相容，则必然存在一个参考态 $\left| b_{j} \right.$，通过若干次平移与之相关：

利用这个关系，我们有

这里再次简化了符号，没有明确表明 $l_{j}$ 依赖于实际状态 $\left| b_{j} \right.$（因此也依赖于 $| a \rangle$）。我们现在可以将 (123) 写为

通过对求和中的指标进行移位，并注意到 $\left| b_{j} \right.$ 可能具有与 $| a \rangle$ 不同的归一化因子（周期性），我们得到

我们现在可以简单地提取出哈密顿算符 $H_{j}$ 的矩阵元：

严格来说，对于非对角键算符 $( j > 0 )$，这不是唯一的矩阵元，因为单独来看它们并不是平移不变的算符。因此，它们在不同动量态之间也有矩阵元。这里我们考虑对所有 $j$ 求和，求和后只有对 $k$ 对角的元素才保留下来。

哈密顿量中可能有若干项对同一个矩阵元 (128) 有贡献，因为有可能（事实上非常可能） $H_{j} | a \rangle \propto T^{- l_{i}} | b_{i} \rangle$ 且 $H_{j} | a \rangle \propto T^{- l_{j}} | b_{j} \rangle$，其中 $\left| b_{i} \right. = \left| \bar{b_{j}} \right.$ （但 $l_{i} \neq l_{j} $）。还要注意，对于作用于两个平行自旋的非对角算符，(125) 中的矩阵元 $h_{j} ( a )$ 为零。还应牢记，即便 (125) 中 $h_{j} ( a ) > 0$，但如果 $\left| b_{j} \right.$ 与动量不兼容，则 $| b_{j} ( k ) \rangle$ 也可能不存在。此时动量态的矩阵元 (128) 就不存在。考虑到这些在后续实现中需要处理的注意事项，(128) 式指定了哈密顿量的所有非零矩阵元。对于这里考虑的具体海森堡模型，我们可以将 $h_{j} ( a )$ 替换为对角和非对角项的实际值。

图27. 8自旋态$| s \rangle = | \uparrow \downarrow \downarrow \uparrow \uparrow \downarrow \downarrow \uparrow \rangle$（顶行）的比特表示，及其平移变换$T^{r} | _{s} \rangle$，其中$r = 1 , \ldots , R ,$，$R = 4$为周期。整数$s ( r )$对应比特的r次连续循环左移。代表为$s ( 1 ) = 51$，即该集合中的最小整数。

对角矩阵元得到

对于其他平移不变算符（这些算符应在动量基中使用）的矩阵元，如果算符传递的动量为0（即，一系列相同的局域算符$O_{i}$，$i = 0 , \ldots , N - 1$之和），则可以完全相同的方式得到。对于传递动量$q \neq 0$的算符$A_{q}$（例如$S_{q}^{z}$，即$S_{i}^{z}$的傅里叶变换），其过程的区别仅在于，当作用$A_{q} | a ( k ) \rangle$时，所得到的基态具有动量$k + q$（在这种情况下，必须存储所涉及的两个动量的基组）。

在程序中构建动量态基组。 我们现在再次使用标签a来指代基组中态的位置，并以整数位$s_{a} , a = 1 , \ldots , M$的形式存储相应的自旋。一个动量态由其代表态$ {\left| s_{a} \right.}$定义，根据(115)，通过平移算符T从该代表态生成完整的态。在这种情况下，我们无法先验地知道基组大小M，因为兼容代表态的数目取决于动量。如果所有态的所有平移都是唯一的，并且没有与动量不兼容的态，那么态的数量将等于$\big [ N ! / \big ( n_{\uparrow} ! n_{\downarrow} ! \big ) \big ] / N$。然而，由于许多态具有小于N的周期性，并且有些态与$k$不兼容（除非$k = 0$），因此这只是一个近似的基组大小。$k = 0$动量区块总是最大的，因为没有与该动量不兼容的态。对于大系统，不允许出现的态和周期$< N$的比例很小，因此上述对基组大小的近似是相当好的。

要在计算机程序中实现动量基组，我们首先需要决定如何选择代表态——原则上，态$| s \rangle$的代表态可以是与其通过平移相关联的态组$\{| s ( r ) \rangle = T^{r} | s \rangle , r = 0 , \ldots , N - 1\}$中的任何一个。图27以一个8自旋态为例，展示了其在位表示下的所有唯一平移。我们必须构建一个代表态列表，并且能够轻松识别与任意自旋态相对应的代表态。如前所述，在使用位表示时，自然地选择平移后的态$| s ( r ) \rangle$（其中整数$s ( r )$最小）作为代表态。

为了生成基组，我们遍历整数$s = 0 , \ldots , 2^{N} - 1$，并且如代码3所示，只处理那些对应选定磁化强度的整数。判断一个态是否是有效的新代表态所需的操作通过位运算来执行。当平移s时，如果某个$s ( r > 0 ) < s = s ( 0 )$，那么$| s \rangle$的代表态已经在基组中，不应再次使用[如图27所示，原始整数$s = s ( 0 ) = 153$，但平移$s ( 1 ) = 51 < s ( 0 )$]。如果所有$s ( r ) \geq s ( 0 )$，则s的代表态尚未出现在列表中。然而，根据(119)，由于其周期性可能与该动量不兼容，s可能仍然是不被允许的。

在下面的伪代码段中，动量仅由乘以(114)中$2 \pi / N$的整数$\in \{- N / 2 + 1 , \ldots , N / 2 \}$表示，我们在此称之为k。我们将定义一个子程序，用于检查态整数s是否为新的有效代表态，并为有效的代表态返回其周期性R。为了对态$| s \rangle$进行平移，我们使用其位的循环排列。我们假设这是通过一个函数cyclebits( i , n )完成的，该函数对整数i的前n位执行向左的循环排列（如图27所示）。在Fortran 90中，可以使用内建函数ishftc( s , d , n )精确完成此操作，该函数将s的前n位循环右移d步。根据(107)，T的定义对应于d = 1。我们的态检查子程序可以实现如下：

subroutine checkstate(s, R)
R = 1; t = s {6}
do i = 1, N
t = cyclebits(t, N)
if (t < s) then
return
elseif (t = s) then
if (mod(k, N/i) = 0) return
R = i; return
endif
enddo

这里整数R首先初始化为1，除非输入状态整数s是所有平移整数中最小的且周期性与动量相容（此时返回时R等于周期性），否则该值将被返回（并随后用作不允许代表的指示符）。现在我们可以轻松构建M个基态的列表；

构建哈密顿量时需要用到周期性列表$R_{a}$。

构建哈密顿量矩阵。 为生成哈密顿量，我们遍历代表态$\vert s_{a} \rangle , a = 1 , . . . , M$。对于每个代表态，检查所有最近邻对$( i , j )$对应的比特位。对角矩阵元的处理与代码5完全相同。对于非对角部分，需要进行少量修改。翻转两个自旋后，得到的态$\left| s^{\prime} \right.$通常不是代表态，需要使用式(124)找到它。矩阵元(130)还需要用于将$\left| s^{\prime} \right.$变换到代表态$\left| s_{b} \right.$的平移次数$l_{j}$。我们将这些任务实现为一个子程序，其内容与构建基态时使用的子程序checkstate（代码$\{6\}$）较为相似。给定状态整数s，我们以所有可能的方式平移其比特位，并存储到目前为止找到的最小整数r（潜在代表态）及对应的平移次数l；

subroutine representative $( s^{\prime} , r , l )$
r = s′; t = s′; l = 0
do i = 1,N 1
t = cyclebits(t, N)
if (t < r) then r = t; l = i endif
enddo

{8}

找到代表态r后，需要在列表$\left\{s_{a} \right\}$中定位其位置b。这与之前相同，通过子程序findstate（代码$\{4\}$）实现。但注意，由于动量施加的周期性约束，潜在代表态r可能实际不在列表中。因此，对findstate稍作修改，若列表中无元素r，则返回$b = - 1$，使用代码4中if...endif语句后的以下代码段：

对于每次自旋翻转，现在可以用以下代码替换代码5中的倒数第二条语句，将贡献添加到哈密顿量：

哈密顿量是复数的，除非$k = 0$或$N / 2$（实际动量为0或$\pi$）。

4.1.4 反射对称性与半动量态#

接下来，除了平移态${\cal T}^{r} | a \rangle$，我们还将考虑由式(109)定义的反射（宇称）算符P生成的态。算符T和$P$不对易，因此初看起来似乎无法构造同时为两者本征态的态。然而，如下文所示，在k=0的动量块中，这些算符是对易的。此外，我们将构造对任意k也是宇称本征态的半动量态。这些态（及哈密顿量）是实值的，与标准复数动量态形成对比。

具有宇称的态。 考虑动量态(115)的以下推广；

其中 $p = \pm 1$。显然，这是一个动量为 k 的态 [即它满足式 (108)]，但它是否也是宇称为 $p$ 的 P 的本征态？我们可以通过对 P 进行显式操作来检查，利用 $P^{2} = \overline{{1}}$、$p^{2} = 1$ 以及关系式 $P T = T^{- 1} P$：

这并不完全具有式 (131) 的形式，除非 $k = 0$ 或 $k = \pi$，因为此时 $\mathrm{e}^{i k r} = \mathrm{e}^{- i k r}$（即动量 k 与 k 等价，并且与该态相关联的方向性不存在）。因此，在这两种特殊情况下，宇称和平移不变性可以同时用于分块对角化，并且 $| a ( k , p ) \rangle$ 确实是一个具有宇称 p 的动量态（换言之，在动量 $k = 0$ 和 $k = \pi$ 的子空间中 $[ \check{T} , P ] = 0$）。

半动量态和宇称。 除了 $k = 0 , \pi$ 这些特殊情况外，宇称不能用于进一步对 H 的动量分块进行分块对角化。然而，我们可以通过混合具有相反的动量 k 的态，以一种不同的方式将宇称与动量结合使用。我们考虑这些态的和与差：

其中 $\sigma = \pm 1$，并且为了方便起见，我们引入了一个函数 $C_{k}^{\sigma} ( r )$：

这里我们将式 (133) 中的 k 称为半动量（semi-momentum）。注意它不是一个守恒的量子数，只是一个指示符，表明动量态是如何组合成半动量态的。严格来说，特殊值 $k = 0$ 和 $k = \pi$ 仍然是常规的（晶体）动量（并且对于这些 k 而言 $\sigma = - 1$ 态不存在），只有 $0 < k < \pi$ 范围内的 k 才应称为半动量。这里我们将所有 $0 \leq k \leq \pi$ 同等对待（需要注意的是，只使用了第一布里渊区的一半，因为另一半对应于相同的态）。对于特殊情况，归一化常数有所不同：

其中我们引入了因子

在实际应用中，$g_{k}$ 在这里并不重要，因为我们只考虑 H 的矩阵元，而它们是 k 对角的。因此，在矩阵元中出现的归一化常数比例中，g 因子会相互抵消。

检查重叠，具有相同 k 但不同 $\sigma$ 的态是正交的：

并且半动量基的其他正交归一化要求 $\langle a^{\tau} ( k^{\prime} ) | a^{\sigma} ( k ) \rangle = \delta_{\sigma \tau} \delta_{k k^{\prime}}$ 也可以很容易地验证。半动量态的优势在于，与复动量态相比，它们对所有 k 都是实值的。

在处理半动量态时，以下等式（标准的三角恒等式）很有用：

其中 $C^{\pm}$ 表示 $\sigma = \pm 1$。利用式 (139) 很容易看出，作用在半动量态上的哈密顿量会混合 $\sigma = \pm 1$ 的态。由 $H_{j} | a \rangle = h_{j} ( a ) T^{- l_{j}} | b_{j} \rangle$ 可得：

从中可以提取哈密顿矩阵元并写为

引入宇称。由于哈密顿量（141）在 方向 上并非对角，半动量块中的状态数是传统动量块中的两倍，因此这样使哈密顿量实化（通过使状态实化）似乎并不能获得比复动量态更多的优势。然而，我们也可以在半动量态中引入宇称，通过定义

这也是一个半动量态，因为反射分量 $P | a^{\sigma} ( k ) \rangle =$ $\sigma | P a^{\sigma} ( k ) \rangle$ 是通过使用反射代表态 $P | a \rangle$ 代替 $| a \rangle$ 得到的半动量态。同样容易验证，对于任何半动量，$( 142 )$ 是宇称算符的本征态；$P | a^{\sigma} ( k , p ) \rangle = p | a^{\sigma} ( k , p ) \rangle$，其中 $p = \pm 1$ 。这仅仅是因为在(142)中，算符 $( 1 + p P )$ 出现在平移算符 $T^{r}$ 之前，而在动量态(131)中，$( 1 + p P )$ 写在 $T^{r}$ 之后。因此，哈密顿量在 $p$ 上是对角的，每个 $( k , p )$ 块中的状态数大约是对应原始半动量 $k$ 块的一半。于是，我们回到了与传统动量态相似的块大小，但状态具有纯实数系数。

半动量态的正交性与归一化。考察状态(142)的正交性时，我们面临一个明显的问题：当存在 $m$ 使得 ${\cal T}^{m} P | a \rangle = \bar{| {a \rangle}}$ 时，状态 $| a^{+} ( k , p ) \rangle$ 和 $| a^{-} ( k , p ) \rangle$ 并不正交。此时我们可以利用(138)和(139)式，将保宇称半动量态(142)写成无宇称半动量态(133)的线性组合：

由此很清楚 $\langle a^{- \sigma} ( k , p ) | a^{\sigma} ( k , p ) \rangle$ 可能不为零。为了更仔细地考察这一点，我们首先确定(142)中的归一化常数 $N_{a}^{\sigma}$。注意我们在 $N_{a}^{\sigma}$ 上标了 $\sigma$ 以表示，保宇称半动量态的归一化常数与纯半动量态的归一化(135)不同，它可以依赖于 $\sigma$（除了对 $k , p$ 的隐式依赖外）。当对所有 $m$ 都有 $T^{m} P | a \rangle \neq$ $| a \rangle$ 时，$N_{a}^{\sigma}$ 的计算很简单。当存在 m 使得 ${\cal T}^{m} {\cal P} | a \rangle = | a \rangle$ 时，我们可以利用(143)以及纯半动量态 $| a^{\sigma} ( k ) \rangle$ 的正交归一性，得到

类似地，当存在 m 使得 ${\cal T}^{m} P | a \rangle =$ $| a \rangle$ 时，我们可以计算 $\sigma = \pm 1$ 状态之间的重叠，得到

这当然是在假设 $\sigma = \pm 1$ 两种状态都存在（即 $N_{a}^{\sigma} \neq 0$）的前提下得出的。因此，在 ${\cal T}^{m} P | a \rangle = | a \rangle$ 情况下，$\sigma = \pm 1$ 状态最多相差一个符号，我们只需在状态块中保留其中一个。为明确起见，若 $N_{a}^{+} \neq 0$，可选择 $| a^{+} ( k , p ) \rangle$，否则选择 $| a^{-} ( k , p ) \rangle$。在后续对半动量态进行一些形式化操作时，我们仍需注意非零重叠 (145)，这一点很快会看到。

半动量哈密顿量。 为了计算哈密顿量矩阵元，我们首先注意到，与普通动量态下的 (124) 式不同，当我们用哈密顿算符 $H_{j}$ 作用于 $| a \rangle$ 时，会得到一个态 $\left| b_{j}^{\prime} \right.$，该态通过若干平移和可能的反射与另一个代表态 $\left| b_{j} \right.$ 相关联：

因此，$H$ 作用于半动量态（或具有宇称的 $k = 0$ 动量态）的结果可以写成以下形式：

其中，代表态 $\left| b_{j} \right.$ 与 $H_{j} | a \rangle$ 相关：

利用关系式 (139)，我们可以将 (147) 写为：

如果对于所有 $m$ 都有 $T^{m} P | b_{j} \rangle \ne | b_{j} \rangle$，则 $\sigma = \pm 1$ 归一化常数之比为 1；否则可写为：

利用重叠 (145)，我们可以提取矩阵元。关于 $\sigma$ 为对角的部分是：

而非对角部分则是：

图28：一个态 $| s \rangle = | \uparrow \uparrow \downarrow \uparrow \uparrow \downarrow \downarrow \rangle$（左上）及其在比特表示下经平移算子T和反射算子P变换的所有状态。同时显示了对应的基10整数s及其变换结果——最小整数 $s = 27$ 对应代表态。对于这一特定态，有 $T^{5} P | s \rangle = | s \rangle$，即归一化常数(144)中 $m = 5$。

简并性。在常规动量基中，$\pm k$ 的态是简并的，但它们是正交态（除非 $k \neq 0 , \pi$）。在半动量基中，k态实际上是同一个态，最多相差一个符号（当 $\sigma = - 1$ 时），这就是半动量希尔伯特空间仅包含 $0 \leq k \leq \pi$ 的原因。简并性在这里转移到了宇称部分，使得 $p \pm 1$ 的态在 $k \neq 0 , \pi$ 时是简并的（但仍正交）。这可以通过在矩阵元(152)中将 $( p , \sigma , \tau )$ 替换为 $( -p , -\sigma , -\tau )$ 看出，该替换不改变 $\tau , \sigma = \pm 1$ 的矩阵元集合。因此，在计算中我们只需要考虑，例如，当 $k \neq 0 , \pi$ 时取 $p = + 1$（而当 $k = 0$ 时，我们只有 $\sigma = \tau = 1$，且 $p = \pm 1$ 部分不是简并的）。

构造半动量态基。推广我们在动量态中使用的约定，现在我们将一个态 $| s \rangle$ 的代表态选为 $| s ( r , q ) \rangle = T^{r} P^{q} | s \rangle$，使得对应的整数 $s ( r , q )$ 最小（其中 $r \in \{0 , . . . , N - 1 \} , q \in \{0 , 1 \} )$。现在，某个块中不能包含某个态，不仅可能因为它与动量不兼容，还可能因为反射量子数 $p$ 和非守恒态标签 $\sigma$ 施加的限制。因此，我们必须根据上文讨论的条件，确定应包含 $\sigma = \pm 1$ 半动量态(142)中的哪一个或全部。图28说明了一个态在比特表示下所有平移和反射的组合，该态与其自身反射的关系为 $| s \rangle = T^{m} P | s \rangle$，在这种情况下，基中应仅包含该态的 $\sigma = - 1$ 或 $\sigma = + 1$ 变体（并且该态与所有动量兼容，因为其周期 $R = N$）。

当同时需要 $\sigma = \pm 1$ 两种态时，我们将存储同一代表态的两个连续副本，这样我们可以继续使用列表 $s_{a}$ 中的位置作为块内每个态的标签。此外，我们还需要存储每个代表态的标签，以及将态反射回到自身的平移次数m，这在归一化和哈密顿量中都需要用到。

如前所述，我们将使用子程序checkstate来确定一个态是否是应添加到列表中的新代表态。除了新代表态的周期R之外，该子程序现在还在适用时给出反射平移数 $m$。我们可以使用 $m = - 1$ 作为标志，表示不存在使 $T^{m} P | s \rangle = | s \rangle$ 成立的m（否则 $m \geq 0$）。对反射态的操作是在执行完代码6中的指令后进行的，首先反射态整数s中的自旋比特并将其存储为t，其中 $t [ i ] = \bar{s} [ N - 1 - i ]$，$i = 0 , \ldots , N - 1$。这一反射由函数reflectbits(s, N)完成。然后，对得到的整数进行与之前相同的平移。因此，我们将代码6修改如下：

{subroutine checkstate(s, R, m) t = reflectbits(s, N); m = 1 do i = 0,R 1 if (t < s) then R = 1; return elseif (t = s) then m = i; return endif t = cyclebits(t, N) enddo}

{11}

这里...表示代码6中的所有操作。需要注意的是，上述代码中的循环从 $i = 0$ 开始，并且在比较 t 与 s 之后执行平移，因为一个状态可以是自身的反射（在这种情况下 $m = 0$）。此外，由于我们已经确定了状态的周期性R，在此循环中只需进行到 $i = R - 1$。函数 reflectbits 需手动实现；通常没有对应的内部函数。在Fortran 90中，可以使用函数 btest(i, b)、setbit(i, b) 和 clrbit(i, b) 来检查、设置和清除整数 i 的单个位 b。

通过上述检查的状态仍然不一定是 $\sigma = \pm 1$ 的有效代表元，因此必须进一步检查，以确定应将代表元的 $\sigma = \pm 1$ 副本中的哪一个（或两个）添加到基中。以下代码段根据上述讨论的标准，存储将包含在区块中的 $\sigma = \pm 1$ 态的状态信息：

call checkstate(s, R, m) do = 1 (若 k = 0 或 k = N/2，仅执行 = +1) if (m = 1) then if (1 + pcos(ikm2 /N) = 0) R = 1 if (σ = -1 且 1 - p cos(ikm2 /N) ≠ 0) R = -1 endif if R > 0 then a = a + 1; s_a = s; R_a = R; m_a = m endif enddo

{12}

回忆这里 m = 1 表示不存在 m 使得 $T^{m} P | s \rangle = | s \rangle$，因此两个 状态都应包含。另一方面，如果存在这样的 m，则只包含其中一个状态——如果其归一化常数 (144) 非零，我们选取 $\sigma = 1$。上述代码通过两个连续的 if 语句实现这一点（其中对无效代表元设置 R = 1）。我们必须存储 ，由于周期性 R 为正且非零，我们可以将其与 一起以 R 的形式存储在一个列表 $R_{a}$ 中。当然，也可以将 R、、m 更紧密地打包成一个整数。

构建半动量哈密顿量。对于一般的半动量，同一个代表元可能在列表 $s_{a}$ 中出现一次或两次，而对于特殊情况 $k = 0 , \pi$，每个代表元总只有一个 $( \sigma = 1 )$ 副本。我们这里构造可以处理两种情况的代码，但对于两个特殊动量可以编写更简单的代码。在通过循环状态索引 a 以及对应的代表元 $s_{a}$ 来构建哈密顿量时，我们首先检查列表中的前一项和后一项。如果前一项 $s_{a - 1} = s_{a}$，则跳过此代表元，因为它已经在之前的循环过程中处理过了。如果下一个代表元 $s_{a + 1}$ 等于 $s_{a}$，则将相同代表元的数量存储为 $n = 2$，否则设置 $n = 1$。哈密顿量构建的这一部分由以下代码完成：

do a = 1,M if (a > 1 且 $s_{a} = s_{a - 1} )$ then 跳过至下一个 a elseif $( a < M$ 且 $s_{a} = s_{a + 1} )$ then n = 2 else n = 1 endif enddo

{13}

这里的 skip 命令会跳过循环的当前迭代，直接进入下一次迭代，而...代表循环的主体部分，我们接下来讨论其主要特征。

首先考虑对角矩阵元 $H ( a , a )$，对于给定的 a，我们先对 $s_{a}$ 中位配置的所有对角贡献求和，并将结果存储为 $E_{z}$。考虑到代表元出现一次或两次的可能性 $( n = 1$ 或 2)，我们可以使用以下循环来相应地分配对角矩阵元：

接下来讨论非对角操作。遍历所有格点对 $(i, j)$，如果比特 $s_{a}[i] \neq s_{a}[j]$，则翻转自旋，得到如代码10所示的状态整数 $s^{\prime}$。再次需要执行所有对称操作以确定哈密顿量矩阵元（151）和（152）中出现的变换指标 $l$ 和 $q$，从而找到对应的代表态 $| r \rangle = T^{l} P^{q} | s^{\prime} \rangle$。实现此操作的子程序是代码8的简单修改：

subroutine representative $( s^{\prime} , r , l , q )$
t = reflectbits(s′, N); q = 0
do i = 1,N 1
t = cyclebits(t, N)
if (t < r) then r = t ; l = i; q = 1 endif
enddo

{15}

这里的 ... 代表代码8中头标之后的部分。调用此子程序后，我们可以在列表 $\{s_{a}\}$ 中搜索代表元 r。这与之前的方法相同，使用子程序 findstate $( r , b )$。

由于哈密顿量比标准动量基更复杂，我们还需要定义函数 helement $( a , b , l , q )$，它根据（151）或（152）返回矩阵元。该函数输入两个代表元（在矩阵中的位置）的标签 a 和 b，以及由 representative 返回的变换指标 $l, q$。矩阵元当然还依赖于动量 k 和宇称量子数 $p = \pm 1$，这里没有明确表示出来。根据（151）和（152）实现该函数是直接的，这里不列出任何代码。

与代码段13中的状态整数 $s_{a}$ 类似，这里也需要考虑代表元 $s_{b}$ 在代表元列表中可能出现一次或两次（下面用 $m = 1$ 或 2 表示）。如果存在两个副本，我们不知道子程序 findstate 返回的是第一个还是第二个位置。因此，再次需要检查与返回位置相邻的代表元列表位置。以下代码段完成所有这些操作，然后分配矩阵元：

$s^{\prime} = \mathbf{flip} ( s_{a} , i , j )$ call representative $( s^{\prime} , r , l , q )$ {16} call findstate(r, b) if $\left( b \geq 0 \right)$ then if $( b > 1$ and $s_{b} = s_{b - 1} )$ then $m = 2 ; b = b - 1$ elseif $( b < M$ and $s_{b} = s_{b + 1} )$ then m = 2 else m = 1 endif do j = b,b + m 1 do i = a,a + n 1 $H ( i , j ) = H ( i , j ) -$ + helement(i, j, l, q) enddo enddo endif

这段代码替代了纯动量基下简单得多的四行代码段 {10}。它的优点是矩阵 H 现在变成了实数值。

4.1.5 自旋翻转对称性#

自旋翻转对称性，其算符 Z 定义在（112）中，可用于减小磁化强度扇区 $m^{z} = 0$ 中的哈密顿量块大小。然而，由于磁化强度守恒，在 $| m_{z} | > 0$ 扇区中使用这种对称性不会带来任何好处（因为这样的基将由 $\pm | m_{z} |$ 状态的混合组成，对称和反对称组合的块大小将与原始块相同）。在不守恒 $m_{z}$ 的模型中，例如横场伊辛模型，自旋翻转可以对所有状态加以利用。我们用 z 表示 Z 的本征值。由于 $Z^{2} = 1$，因此再次有 $z = \pm 1$。

与上一节考虑的宇称算符不同，（112）中的自旋翻转算符 Z 与平移算符对易，$[ T , Z ] = 0$。因此，在所有动量扇区中，关联的量子数 z 与 k 一起守恒。在磁化强度扇区 $m_{z} = 0$ 中，我们可以通过使用以下形式的状态，将动量块始终分成两个更小的块：

其中归一化常数很容易求得为：

对 $| \alpha \rangle$ 进行哈密顿量操作，结果为

我们得到矩阵元

这对任意 k 均成立。和往常一样，对于 $j = 0$，上式简化为 $h_{0} ( a )$。

具有半动量态的自旋翻转对称性。我们也可以考虑结合自旋翻转对称性的半动量态；

这些态同时是 $Z$ 和 $P$ 的本征态。

在计算归一化和构造哈密顿量时，我们现在必须考虑五种不同类型的参考态，具体取决于哪些对称性操作组合能将参考态变换为自身；

这里，不等式应对所有 m 成立，等式则对某些 m, n 成立。计算遵循前一节的过程，我们仅列出结果。归一化常数为：

对于情形 $^{2 ) , 4 )}$ 和 5)，应由一对 $\sigma = \pm 1$ 的态中仅选取一个包含在基组中，因为这样一对内的两个态仅仅相差一个因子。为了明确，我们可以再次选取 $\sigma = 1$（如果这能使 $N_{a}^{\sigma} > 0$），否则选取 $\sigma = - 1$。现在，将哈密顿量算符作用于一个代表元后，会得到一个新的代表元 $\left| b_{j} \right.$，其形式为

H 的对角线矩阵元为

而在 $\sigma$ 上非对角的矩阵元为，

这些表达式看起来可能相当复杂，但需要注意的是，它们完全适用于一维系统，而不仅仅限于此处考虑的海森堡链。一旦它们在某个模型中得到实现和测试，代码可以很容易地复用于其他系统。

程序实现；组合所有对称性。在具有 $m_{z} = 0$ 的半动量基（semi-momentum basis）中纳入自旋反转对称性，其中态 $| a^{\sigma} ( k , p , z ) \rangle$ 由 (157) 式定义，我们需要增强代码段 11 中的子程序 checkstate。它不仅要返回对应于代表在反射 P 下的平移数 $T^{m_{p}} P | s \rangle = | s \rangle$，还要返回涉及 Z 和 $P Z$ 的对称关系中的平移数：$\dot{T^{m_{z}}} Z | s \dot{\rangle} = | s \rangle$ 和 $T^{m_{p z}} P Z | s \rangle = | s \rangle$。为了在位表示中实现自旋反转，我们定义一个函数 invertbits(s, N)，它将所有位 $s [ i ]$ 翻转：$1 - s [ i ] , i = 0 , N - 1$。这可以通过不使用逐位操作的简单方式实现，例如 $s = 2^{N - 1} - s$（当 N 小于所用整数位宽时；否则必须使用不同的公式）。代码 6 和 11 可以轻松修改，以同时检查反转态和反射-反转态，从而确定 $m_{z}$ 和 $m_{z p}$。

表 1. 磁化 $m_{z} = 0$ 时，具有不同宇称和自旋反转量子数 $( p , z )$ 的 $k = 0$ 态块的大小

如果修改后的 checkstate 返回一个潜在的代表 $R \neq - 1$，那么需要进一步检查该态是否满足作为 $\sigma = \pm 1$ 的一个或两个态代表的其他条件。现在，一个允许的代表根据其对称性质（如 (158) 式总结的那样）属于五类不同的之一。类别 $c \in \{1 , . . . , 5 \}$ 可以由 checkstate 提供的平移整数 $m_{p} , m_{z} , m_{z p}$ 确定。确定 c 之后，三个平移整数可以缩减为两个；即归一化常数 (159) 中的 m 和 n。类别 c 可以连同 m、n 打包成一个整数，例如以 $m + n ( N + 1 ) + c ( N + 1 )^{2}$ 的形式，并与周期性和索引（像之前一样打包为 R）一起存储。

为了构造哈密顿量，需要进一步修改子程序 representative（其主要部分已在代码段 8 和 15 中描述），使其除了返回 q 次反射和 l 次平移外，还返回一个索引 g，对应于将自旋翻转生成的态 s 变换为相应代表所需的 0 或 1 次自旋反转；即 $| r \rangle = T^{l_{j}} P^{q_{j}} Z^{g_{j}} | s \rangle$。然后哈密顿量的构造与之前完全相同，只需修改函数 helement，将 g 作为代码段 (5.15) 的一个参数。这是直截了当的，此处无需列出任何代码。

态块大小的示例。我们对于在一维系统中应用对称性的讨论到此为止。对于海森堡模型，总自旋也是守恒的，但在基中纳入该对称性要复杂得多。原则上，对于总自旋单重态，可以使用价键基（valence-bond basis）[144] 来实现，但实际上很少这样做，因为由此产生的哈密顿量非常稠密，不像此处实现的对称性所得到的稀疏矩阵（对于一个 $M \times M$ 矩阵，非零元素的数量与 $M \times N$ 成正比，且 $N \ll N$）。这种稀疏性在 4.2 节讨论的兰佐斯（Lanczos）计算中将非常有用。

表1显示了在利用所有对称性时，对于几种不同链长N，希尔伯特空间分块的大小。此处动量$k = 0$，并给出了所有宇称和自旋反转量子数组合$( p = \pm 1 , z = \pm 1 )$下的态数。最大的分块始终是完全对称的那一个。随着系统尺寸增大，分块大小的相对变化减小，对于较大的N，所有尺寸都能很好地近似为$N ! / [ 4 ( \frac{N} {2} ! )^{2} N ]$。对于$k = \pi$，分块的大小与表中$k = 0$的情况大致相同；对于其他动量，由于$\sigma = 1$和$\sigma = -1$的态都是允许的，分块的大小大约是两倍。

4.1.6 期望值与热力学#

现在我们已经准备好对角化哈密顿量并计算物理可观测量。这里我们首先考虑完全对角化，这意味着我们计算H的所有本征值和本征态。原则上，非奇异的$M \times M$矩阵H的本征值$\lambda_{n} , n = 1 , . . . , M$可以通过求解久期方程得到：

其中det[]表示行列式，I是单位矩阵。随后可以通过求解线性方程组得到本征向量$\nu_{n}$：

然而，由于对于大型矩阵，久期方程非常复杂，这种方法在实践中并不使用。大多数数值矩阵对角化方法基于迭代搜索一个酉矩阵U，使得：

其中E是对角本征值矩阵，$E_{n n} = E_{n}$。对于复厄米矩阵H，酉矩阵U的逆$U^{- \bigtriangledown}$是其复共轭矩阵的转置：$U^{- 1} = U^{* T}$。如果矩阵是实对称的，则有$\boldsymbol{U}^{- 1} = \bar{\boldsymbol{U}}^{T}$。

对角化矩阵U的列包含了H的本征向量。这可以通过将(165)式左边乘以U看出，得到$H U = U H$。由于E是对角矩阵，该等式右边$U$的第n列被第n个本征值（即矩阵元$E_{n n} = \lambda_{n} $）相乘。在左边乘以H时，U的第n列给出了$H U$的第n列$( H U )_{n}$，即$H U_{n} = E_{n n} U_{n}$，因此H的本征向量$\nu_{n}$被识别为$\nu_{n} = U_{n}$。

某个可观测量（算符）A在H的一个本征态中的期望值由对应矩阵变换到能量基后的对角元给出：

人们通常关心所考虑分块的基态（假设本征值从最低到最高排列），特别是所有分块中能量最低的那个实际基态。对于热平均，需要所有对称分块的所有本征值和/或本征态：

表2. 动量 $k = 0$ 的16位点链在宇称和自旋反转量子数 $( p = \pm 1 , z = \pm 1 )$ 块中的八个最低能量

这里 $\beta = T^{- 1}$ 是逆温度（单位采用 $k_{B} = 1$ ），指标j集体表示大小为 $M_{j}$ 的块的不同量子数 $m_{z} , k , p , z$。

矩阵对角化是一种标准的线性代数运算，许多软件库中都提供了成熟的子程序。因此，我们在此不讨论矩阵对角化过程的具体内部机制。通常，对角化程序会以升序返回一个包含特征值的向量，以及一个由相应特征向量组成的矩阵U。

对角线化一个 $M \times M$ 矩阵所需的运算次数通常与 $M^{3}$ 成比例，而存储所需的内存约为 $\sim M^{2}$ （即使对于稀疏矩阵，中间步骤通常也无法保持稀疏性）。这严重限制了实际中可以完全对角化的矩阵大小。目前，在单个工作站上，可以较为轻松地处理 $M \approx 10^{4}$ 的矩阵，而在超级计算机上可以处理略大几倍的矩阵。查看表1，可以清楚地看到，对于超过 $N \approx 20$ 的海森堡链，实际上无法进行完全对角化。对于旨在计算更大系统的基态以及可能的一些激发态的计算，可以采用其他方法，例如第4.2节讨论的Lanczos方法。

总自旋。了解能量本征态的总自旋S是很有用的。由于我们在对哈密顿量进行块对角化时没有利用S的守恒性，因此必须利用这些态来计算 $\mathbf{S}^{2}$ 。如前所述，该算符在形式上等价于一个具有长程相互作用的海森堡哈密顿量，其在(113)中显式写出。我们可以使用构建哈密顿量时开发的略有修改的程序来获得 $\mathbf{S}^{2}$ 的矩阵。然后我们用 $U$ 对其进行变换，之后我们就可以像(166)中那样，将 $S ( S + 1 )$ 表示为对角矩阵元。

让我们看一下在 $k = 0$ 动量扇区中，一个16位点链获得的一些结果。表2中列出了每个对称扇区 $\left( p , z \right)$ 的八个最低能量，以及计算得到的这些态的自旋。最低能量是一个 $( p = 1 , z = 1 )$ 的单重态。系统的第一激发态实际上并不在 $k = 0$ 扇区；它是一个 $k = \pi$ 的三重态，具有 $( p = - 1 , z = - 1 )$ ，能量为 $E = - 6.87210668$。

表中值得注意的一个有趣特征是，在每个对称扇区中，所有态的自旋要么是偶数，要么是奇数。这是由于自旋反转对称性。一个具有N个位点和总自旋S的 $m_{z} = 0$ 态可以写成由 $N / 2 - S$ 个单重态 $( | \uparrow_{i} \downarrow_{j} \rangle - | \downarrow_{i} \uparrow_{j} \rangle ) / \sqrt{2}$ 和S个三重态 $( | \uparrow_{i} \downarrow_{j} \rangle + | \downarrow_{i} \uparrow_{j} \rangle ) / \sqrt{2}$ 构成的线性组合[48]。

由于单重态在自旋反转下是奇的，而三重态是偶的，因此对于一个N位点系统，自旋S态的自旋反转量子数$z = +1$当且仅当$N/2$和S同为偶数或同为奇数，否则$z = -1$。因此，即使我们不计算S值，也能根据z值获得关于S的某些有限信息。低能态通常也具有低自旋，即$S = 0$或1，这些态可以仅通过z来区分。但需要注意的是，给定$(p, z)$扇区中的最低态并不一定具有最低的（奇或偶）自旋，如表2中最低的$(p = -1, z = 1)$能级所例证的，其自旋为$S = 2$，而非$S = 0$。然而，这并非一个低能态，因为在$k = 0$扇区中有11个能量更低的状态（并且在其他动量扇区中还有更多）。$k = 0$扇区中能量第二低的状态也具有$S = 2$，但在其他动量扇区中有四个$S = 1$的状态的能量低于它。

磁化率和比热。作为热力学的一个例子，让我们计算海森堡链的两个重要性质：比热C和磁化率$\chi^{z}$。它们可以根据以下公式进行评估：

两者都可以利用(167)方便地导出。在$\chi$的定义中，$h$是z方向的外磁场，它通过哈密顿量中增加的一项$-h m_z$耦合。这里我们只考虑零场情况，因此$\langle m_z \rangle = 0$。

C和$\chi$都是特殊的物理量，因为它们只依赖于能谱而不依赖于态。对于海森堡模型，由于总自旋S守恒，我们可以只考虑$m_z = 0$的能级来计算这些量。一个自旋为S的态具有$(2S+1)$重简并，磁化强度$m_z = -S, ..., S$。因此，对于C，我们只需对$m_z = 0$能级求和(168)，并用这个简并因子加权，即可得到对所有$m_z$的全平均。对于$\chi_z$，我们可以用同样的方法，同时利用在自旋S多重态中，(169)中的平均$\langle m_z^2 \rangle = \langle \mathbf{S}^2 \rangle / 3 = S(S+1)/3$这一事实。

情况稍微复杂之处在于，可能存在与哈密顿量明显对称性无关的“偶然”简并。如果这些简并态具有不同的S，那么对角化过程会给这些态带来一些随机的（取决于具体对角化方式）混合自旋。原则上，此时应在简并子空间中对角化$\mathbf{S}^2$以获得本征自旋值。然而在实践中，这种偶然简并非常罕见。对于小尺寸海森堡链，它们只出现在$k = \pi$, $p = 1$, $z = 1$的扇区。对于$N=10$，存在两个能量$E=0.5$的此类态；对于$N=12$和14，没有偶然的混合S简并；对于$N=16$，在能量2J处有四个混合的能级。由于简并数量如此之少，且能量如此之高，我们可以忽略它们。

图29说明了C和$\chi$的结果如何随系统尺寸增大而收敛到热力学极限。计算到$N=16$时，结果在$T/J \approx 0.25$以下已很好地收敛。高温下收敛更快是自然的，因为在$T \to \infty$极限下，自旋在统计上是退耦合的，系统表现为一组独立的自旋。因此，在$T \to \infty$时，每个自旋的磁化率具有单自旋居里形式：$\chi \sim 1/(4T)$。对于C，通过精确计算2-位点系统的能量（单个自旋能量恒定，因此$C=0$），求导，并仅保留$T^{-1}$展开的首项，可以得到一个解析的高温形式。这给出当$T \to \infty$时，$C \sim 3/(13T^2)$。高温极限也在图29中示出。对于C，该形式在$T/J \approx 1.5$以下都能很好地描述其行为，而对于$\chi$，只在更高温度下才吻合良好。

图29. 长度为$N=4, 8, 12$和16的海森堡链的比热和磁化率随温度的变化。细曲线表示主导的高温形式。

当$T \to 0$时，对于有限系统，C和$\chi$都指数趋近于零，这是由于基态与激发态之间存在有限能隙（因此当$T \to 0$时，只有基态对热平均有贡献），并且对于$\chi$，当N为偶数时，海森堡链的基态总是单重态（而当N为奇数时，基态是$S=1/2$的二重态，因此$\chi$按$T^{-1}$发散）。在无限系统中，$C \to 0$，但$\chi$趋近于一个非零的常数。

请注意，这里以及我们之后将遇到的许多情况中，在讨论 $T = 0$ 的性质时，考虑取极限 $T \to 0$ 和 $N \to \infty$ 的顺序非常重要。对于有限（偶数）$N$ 的链，其磁化率在 $T \to 0$ 时总是趋近于零。磁化率开始指数下降到零的温度取决于有限尺寸的自旋能隙，在海森堡模型的情况下，这个能隙 $\propto N^{- 1}$（我们将在第 4.3.1 节中进一步讨论并用数据说明）。因此，对于给定的 $N$，可以预期在直到某个 $T / J \propto N^{- 1}$ 的温度范围内，结果基本与热力学极限一致（在图29中也可以大致看到这一点），并且如果先取 $N \to \infty$ 的极限再取 $T \to 0$ 的极限，那么磁化率在 $T = 0$ 时保持非零。

在低温下，无法使用精确对角化结果来可靠地外推海森堡链在热力学极限下的性质。对于磁化率，早在几十年前就已经有人尝试使用特定的函数形式 [147] 对精确对角化结果 [146] 进行外推——其结果被称为 Bonner-Fisher 磁化率。此外，该外推还利用了已知的 $T \to 0 ( N = \infty )$ 时 $\chi$ 的精确值，该值来自基态的精确 Bethe 拟设解。Bonner-Fisher 曲线常被用于分析准一维反铁磁体（弱耦合链）的实验结果。这里最重要的特征是在 $T / J \approx 0.6$ 处最大值的位量和形状，如图29所示，这个值在较小的系统中已经很好地收敛了。现在已经有能力使用其他计算方法来研究非常长的海森堡模型的各种性质（我们将在第 5.3.1 节中讨论），因此可以更可靠地研究低温极限。此外，通过结合 Bethe 拟设 [7] 的场论研究，也获得了本质上精确的结果。这些计算显示了一个反常特征；趋近 $T \to 0$ 极限的过程是对数形式的，因此很难用有限 $N$ 的计算来重现。因此，Bonner-Fisher 曲线的低温形式是不正确的。对于实验来说，最重要的特征是那个宽峰。通常，在更低的温度下，超出海森堡链的效应（例如，链间耦合、自旋各向异性或无序）会出现。一些准一维材料的自旋磁化率实际上遵循海森堡链的 $\chi$ 一直降到 $\approx J / 50$ 的温度，在这个温度下对数行为非常显著 [8]。

4.2 Lanczos 方法#

正如我们在前一节中所看到的，对于自旋数超过 $\approx 20$ 的 $S = 1 / 2$ 系统，对哈密顿量（或各个区块）进行完全对角化变得耗时到不可行。对于更高的 $S$，情况当然更糟。如果我们只限于研究基态，并且可能还包括一些低能激发态，那么通过使用 Krylov 子空间技术，例如 Lanczos 方法，我们可以将可处理的系统尺寸扩大大约一倍。

4.2.1 Krylov 子空间#

Krylov 子空间是全希尔伯特空间的一个子空间，其构造方式使得感兴趣的哈密顿量 H 的低能态在其中可以被很好地近似。考虑一个任意态 Ψ，例如，其矢量元素 $\Psi ( i ) , i = 1 , \ldots , M$ 随机生成的态，所在空间是定义 H 的 M 维希尔伯特空间，并将该态按照哈密顿量的本征态 $\left| {\psi_{n}} \right.$ 展开（按本征值递增顺序排列，此处标记为 $n = 0 , 1 , \ldots , M - 1 )$ 。我们用哈密顿量的幂作用于这个选定的态；

如果幂次 Λ 很大，对应最大模本征值 $\lvert E_{\mathrm{max}} \rvert$ (即 $\operatorname*{max}=0$ 或 $\operatorname*{max}=M-1$) 的态将在求和式中占主导，前提是展开系数 $c_{\operatorname*{max}} \neq 0$。因此，对态多次作用哈密顿量会投影出本征值为 $E_{\mathrm{max}}$ 的本征矢量。如果我们希望确保得到基态 $\left| \psi_{0} \right.$ (当 $\Lambda \to \infty$ 时精确得到，有限 Λ 时近似得到)，可以用 $(H-c)^{\Lambda}$ 代替 $H^{\Lambda}$，其中 c 是一个足够大的正数，以确保 $\left| E_{0} - c \right| > \left| E_{M-1} - c \right|$。在下文中，我们将假设如果需要这样的常数，它已被吸收到 H 中。

虽然当 Λ 足够大时，式(170)保证能产生基态，但一种更有效的构造当 $\Lambda \to \infty$ 时趋近基态的状态的方法，不仅是考虑 $H^{\Lambda} | \Psi \rangle$，而且要考虑由态集合 $H^{m} | \Psi \rangle, m=0,...,\Lambda$ 张成的整个希尔伯特子空间。这些态可以通过对初始态 $\left| \Psi \right.$ 连续作用 H 来逐个构造。在这个子空间中，存在一个最佳线性组合的向量来近似基态（最小化能量期望值），找到它的方法是在生成的 $\Lambda+1$ 个向量子空间中对 H 进行对角化。除了在相对较小的 Λ（通常比整个希尔伯特空间的大小 M 小多个数量级）下投影出基态外，这种方法还可以精确地再现一些低能激发态。

通过对初始态多次作用 H 获得的希尔伯特子空间称为 Krylov 子空间。我们这里讨论的是量子力学中的哈密顿量，但 Krylov 子空间方法当然也更普遍地适用于本征值问题，并在科学和工程的许多领域中广泛使用。

4.2.2 Lanczos 基与哈密顿量#

在 Lanczos 方法[148]中，利用 Krylov 子空间态的线性组合构造了一组正交基，使得在此基中写出的哈密顿量是三对角的。在下面将要描述的标准方法中，首先构造一组正交但未归一化的基 $\{| f_{m} \rangle \} , m=0,...,\Lambda$。随后将这些态归一化，得到一组正交归一集 $\{| \phi_{m} \rangle \}$，其中哈密顿量呈现特别简单的三对角形式。我们还将讨论一种略有不同的方法，直接生成归一化基 $\{| \phi_{m} \rangle \}$。

Lanczos 基的生成。 基 $\left\{\left| f_{m} \right. \right\}$ 的构造从一个任意归一化态 $\left| f_{0} \right.$ 开始，该态的唯一要求是它与 H 的基态不正交（如果目标是寻找基态）。随机生成的 $\left| f_{0} \right.$ 应满足此条件，但我们也可以从一个已知与基态有较大重叠的向量开始。下一个态由下式给出：

其中常数 $a_{0}$ 应使得 $\left| f_{1} \right.$ 与 $\left| f_{0} \right.$ 正交。为此，我们来检查两个态之间的重叠：

这里我们引入了归一化常数和哈密顿量对角矩阵元的符号，这些符号也将用于后续的态：

为了使重叠(172)为零，我们必须选择 $a_{0}=H_{00}/N_{0}$。接下来的 $m=2$ 的态用前两个态表示为：

其中 $a_{1}$ 和 $b_{0}$ 可以选择使得 $\left| f_{2} \right.$ 与 $\left| f_{0} \right.$ 和 $\left| f_{1} \right.$ 都正交。利用从式 (171) 和 (174) 得到的 $H | f_{0} \rangle$ 和 $H | f_{1} \rangle$，重叠积分计算如下：

因此，适当的系数是 $a_{1} = H_{11} / N_{1}$ 和 $b_{0} = N_{1} / N_{0}$。对于所有后续迭代，式 (174) 推广为

并且使该态与 $\left| f_{m} \right.$ 和 $| f_{m - 1} \rangle$ 正交的系数是我们已找到的 $a_{0} , a_{1}$ 和 $b_{0}$ 表达式的推广：

这可以通过直接计算重叠积分来轻松验证。接下来需要证明，使用这些系数得到的态 $\left| f_{m + 1} \right.$ 也与所有 $k < m - 1$ 的先前态 $\left| f_{k} \right.$ 正交。在归纳证明中，我们可以利用所有先前生成的态彼此正交这一事实，得到

其中我们还使用了式 (176) 形式的 $H | f_{m - k} \rangle$。

这个迭代过程 (176) 一直持续到 $m = \Lambda$，其中 Λ 可以根据计算量的某个收敛准则自动、动态地确定，我们将在下文中讨论。首先我们需要哈密顿量矩阵元。

Lanczos基中的哈密顿量。 构造出一组 $\Lambda + 1$ 个 Lanczos向量 $| f_{m} \rangle , m = 0 , \ldots , \Lambda$ 后，我们可以构建该基中的哈密顿量。利用式 (176) 得到的 H 作用于一个基态的表达形式：

因此非零矩阵元为

归一化的基态为

根据式 (177)，$b_{m} = N_{m + 1} / N_{m}$。因此，非零矩阵元为

此矩阵是三对角矩阵，可以使用比通用对角化算法（许多线性代数子程序库中都有）更快的特殊方法进行对角化。然而，主要优势并不在于该矩阵更容易对角化——它通常规模不大，即便使用通用对角化程序，花费在该部分计算上的时间也可以忽略不计。其优势在于该矩阵的构建相对迅速，特别是当矩阵 $H$ 是稀疏矩阵时，只需考虑其非零元素（在作用于态时即时存储或生成）。有时，当将方法应用于非常大矩阵时，过程仅需存储三个大规模态矢量 $\left| f_{m} \right.$ 也是很有用的。在计算谱函数（动力学关联）时，三对角化也非常方便，例如参考文献 [43, 143] 中的讨论。

注意，当基态用 $m = 0 , \ldots , \Lambda$ 标记时，三对角矩阵的大小（实际基的大小）是 Λ，而不是 $\Lambda + 1$，因为在 (179) 式中当 $m + 1 = \Lambda$ 时，最后一个为用于 (183) 式生成的对角矩阵元是 $a_{\Lambda - 1}$。态 $\left| f_{\Lambda} \right.$ 不需要构建；只需要项 $H | f_{\Lambda - 1} \rangle$ 来计算最后一个系数 $a_{\Lambda - 1}$。

由于对角化兰佐斯哈密顿量矩阵的成本非常低，我们可以在每个新的基态生成后执行对角化，并跟踪本征值的演化。我们可以根据合适的收敛准则停止该过程，例如，期望的本征值在迭代之间变化小于某个容差。收敛所需的基大小当然取决于所研究的系统，但通常对于量子自旋系统，Λ 在几十到几百的范围内就足够了。

使用归一化向量的替代公式。兰佐斯基的生成通常以非归一化态 $\left| f_{m} \right.$ 来讨论，正如我们上面所做的那样。然而，该过程的直接计算机实现有时会导致数值问题，因为归一化常数 $N_{m}$ 可能变得非常大（如果 $H$ 的本征值很大）。那么直接使用归一化态 $\vert \phi_{m} \rangle$ 可能会更好（作为将 $H$ 乘以合适因子的替代方法）。实际上，在实践中这个公式甚至更简单，额外的计算成本仅仅是每一步的归一化（与其他成本相比，这通常很小）。

我们从一个归一化态 $\left| \phi_{0} \right.$ 开始，并按以下方式生成第二个态：

这里 $N_{1}$ 是一个归一化常数，通过直接计算构建的态 $H | \phi_{0} \rangle - a_{0} | \phi_{0} \rangle$ 与自身的标量积确定。与 $\left| \phi_{1} \right.$ 的正交性再次要求 $a_{0} = H_{00}$。对于接下来的每一个态，可以很容易地证明

与所有先前的态正交。这里 $a_{m}$ 和 $N_{m}$ 的定义与 (177) 式中的先前定义不同；

其中 $\left| \gamma_{m} \right.$ 是归一化之前的生成态，如上面 (184) 式所示。哈密顿量矩阵元为

在这种表述中，所有存储的数值都是良态的。

简并态。需要注意的是，Lanczos方法无法从多重态中产生多个成员；对于一组简并态，只能得到这些态的特定线性组合（该组合取决于初始态 $\left| f_{0} \right. )$ ）。要理解其原因，我们再次考察态 $H^{\Lambda} | \Psi \rangle$ 的展开式(170)，并假设其中存在两个简并态 $\left| {\psi_{i}} \right.$ 和 $| \psi_{j} \rangle , E_{i} = E_{j}$。在展开式中，我们可以将这些态与其他项分离：

对于任意Λ，展开式都包含态 $| \psi_{j} \rangle$ 和 $\left| {\psi_{j}} \right.$ 的相同线性组合。因此，在由态集合 $H^{m} | \Psi \rangle , m = 0 , \dots , \Lambda - 1$ 张成的子空间中，无法得到这两个简并态的不同线性组合。当然，这也能推广到包含超过两个态的简并多重态。

正交性的丧失与恢复。当基矢数量Λ很大时，Lanczos过程通常会遇到数值不稳定性。在构造基矢组的过程中累积的舍入误差，最终会引入态之间的非正交性。这种数值误差会升级，并导致在Λ的某些狭窄范围内，连续地突然出现多个相同的本征值（回顾：Lanczos方案永远不会产生简并态）。我们将在下文看到一个例子。

当目标仅是获得基态时，正交性的丧失以及相同态的多份拷贝出现通常不是问题。然而，它会使激发态的计算变得复杂。为了解决这个问题，可以在基矢构造中添加一个步骤，即显式地使每个新构造的Lanczos向量与所有先前的基矢向量正交化。在这种稳定化过程的最简单实现中，所有Lanczos向量都被存储起来（在主存或辅存中）。当处理非常大的基矢组时，这可能会出现问题，但该方案非常简单。对于归一化的态，根据(184)构造的每个 $| \phi_{m + 1} \rangle$ 需要与所有先前的态正交化，按如下步骤进行：

对 $i = 0 , \ldots , m$ 依次执行。这使得研究数量多得多的激发态成为可能（原则上仅受限于计算机内存）。

除了在同一次运行中收敛多个激发态，也可以逐个瞄准激发态，每次都从一个与先前所有态正交的向量开始。此时，也需要针对这些态进行再正交化。

本征态与期望值。对三对角Lanczos哈密顿量进行对角化，得到本征值 $E_{n}$ 和本征向量 $\nu_{n} , n = 0 , \ldots , \Lambda - 1$。我们希望这些本征向量用原始基 $\{| a \rangle \}$ 表示，因为在该基下，我们能够计算感兴趣的算符 O 的矩阵元 $\langle b | O | a \rangle$。首先，Lanczos基态为

我们将哈密顿量的期望本征向量记为

三对角矩阵的前几个本征向量 $\nu_{n}$ 在Lanczos基下能够准确地表示哈密顿量的本征态（对于足够大的Λ，基本上是精确的）；

因此，我们想要构造的波函数系数由下式给出：

如果在基矢构建过程中所有Lanczos向量均已存储，则此变换可直接进行。若计算机内存充裕，应存储所有Lanczos向量；但通常计算已逼近计算机容量极限，此时只能存储最少信息，代价是计算时间延长。若基矢构建时未存储完整的Lanczos向量集，则需要以与之前相同的迭代方式重新生成它们，唯一的区别是系数 $a_{m}$ 和 $N_{m}$（如果使用未归一化Lanczos向量 $\left| f_{m} \right.$ 的方案，则 $b_{m}$ 也可用）已经获取，无需重新计算。我们可以利用本征向量 $\nu_{m}$ 在运行中根据(192)式变换这些态，逐步构建哈密顿量的一个或多个本征态。

生成一个或多个本征态后（假设这些态已按(190)式中系数 $\psi_{m} ( a )$ 的形式存储），某个算符 $\cal O$ 的期望值可通过先作用于此态来获得，得到一个我们称之为 $\vert \psi_{n}^{\cal O} \rangle$ 的未归一化态，

然后我们计算该态与原态 $\left| {\psi_{n}} \right.$ 的内积；

这种略显繁琐的期望值写法对应一个典型的计算流程：先计算状态向量 $\psi_{n}^{O}$，然后计算其与 $\psi_{n}$ 的标量积。

值得注意的是，在Lanczos计算中，我们经常处理四种不同的基：首先，我们有原始的单自旋态“计算”基。从中我们生成包含对称性的基，例如动量态（通常就是上面记为 $| a \rangle$ 的态）。然后我们构造Lanczos向量，它们是这些态的特殊线性组合。最后，对角化三对角矩阵（它是低能区域的有效哈密顿量），我们得到所需的能量本征态。为了用这些态进行计算，我们实际上反向进行基变换。(193)和(194)中的矩阵元最终是在单自旋的计算基中执行的（即，使用例如用于构造动量态的代表态）。接下来，我们将在实际程序实现中讨论所有这些流程。

4.2.3 Lanczos方法的编程#

虽然Lanczos方法可以应用于任何对称（或厄米）矩阵，但对于自旋哈密顿量，还有一个额外的优势：哈密顿量是稀疏矩阵。尽管哈密顿量（单个对称块）的尺寸可能非常大，但非零矩阵元的数量要少得多。对于N个格点上具有短程相互作用的模型，一个包含M个态的哈密顿量块大约有NM量级的非零元，对于大的M而言，这远小于总元素数 $M^{2}$ 。由于构造Lanczos基最耗时的部分是哈密顿量矩阵对状态向量的重复作用，以根据式(176)或(184)生成下一个基态，稀疏性使得时间节省巨大。同样，在内存节省方面也有显著优势。因此，哈密顿量的非零元素应以紧凑形式存储，或根据需要即时生成（当使用许多对称性时，这通常比甚至从磁盘读取元素还要耗时）。我们在此实现Lanczos方法，将哈密顿量以紧凑形式存储在主存中。

我们将使用式(184)和(183)生成每一步都归一化的Lanczos基态 $\{| \phi_{m} \rangle \}$ 。在计算机程序中，这些态以其向量分量 $\phi_{m} ( a )$ 的形式存储，Lanczos态由式(189)用这些分量给出。这里 $a = 1 , \ldots , M$ 标记我们工作基的态，在最大程度使用对称性的情况下，这将是动量或半动量态，例如 $| a \rangle = | k , p , z , m_{z} = 0 , a \rangle$ ，但为简单起见，我们不写出所有量子数。

基态的性质仅在哈密顿量（或要测量的某个算符）作用于基态时才会起作用，为此我们需要的矩阵元正是已经在4.1.3、4.1.4和4.1.5节中针对各种对称性实现讨论过的那些。我们将编写假设波函数系数为实数的伪代码；复状态所需的更改不言自明。

Lanczos基构造与本征态。为了讨论Lanczos程序的一般结构，我们首先假设已经实现了一个子程序 hoperation $( \phi , \gamma )$ ，它将哈密顿量作用于一个态向量； $| \gamma \rangle = H | \phi \rangle$ 。稍后，我们将描述这个子程序的实现，它是Lanczos程序中唯一需要指定模型和所实现对称性的部分。我们还使用一个子程序 normalize $( \phi , n )$ ，它首先计算 $n = \langle \phi | \phi \rangle$ ，然后重新缩放 $\phi$ ，使得返回时 $\langle \phi | \phi \rangle = 1$ 。

为了开始基的生成，我们首先将初始态 $\left| \phi_{0} \right.$ 加载为一个随机生成的归一化态。我们将分量记为 $\phi_{0} ( i ) , i = 1 , \ldots , M$ 且 $\langle \phi_{0} \vert \phi_{0} \rangle =$ $\begin{array} {r} {\dot{\sum}_{i} \phi_{0} ( i )^{2} = 1} \end{array}$ 。下一个Lanczos态根据(183)生成；

如果不想存储所有Lanczos向量，可以在实现迭代基生成时只循环三个向量 $\phi_{0}, \phi_{1}, \phi_{2}$，它们包含每一步所需的信息，如(184)所示：

如有可能，应存储所有态，以避免在后续计算期望值时重新生成（这会使计算时间加倍）。如果希望进行额外的重正交化（re-orthogonalization），以确保数值截断误差不会最终破坏Lanczos基，就必须存储所有态。此时，需要根据方程(188)插入一个额外的循环，与所有先前生成的态进行正交化：

完成代码段18或19后，三对角矩阵（式(187)）的内容已经就绪，接下来可以将这些信息输入对角化程序，该程序将给出特征值 $E_{n}$ 和特征向量（元素为 $\nu_{n} ( i )$，$n = 0 , \ldots , \Lambda - 1 , i = 1 , \ldots , M$）。需要再次指出，根据我们一直使用的约定，Lanczos基的规模（即矩阵的规模）为 $\Lambda$；在18和19中生成的最后一个态（其中 $m + 1 = \Lambda$）实际上并不需要（只需要最后一次迭代中生成的系数 $a_{\Lambda - 1}$）。

与花在基构造上的时间相比，对角化的计算量非常小。因此，我们完全可以在每次生成新的Lanczos向量后立即执行对角化。这样，我们可以监控能量随基规模的变化。当满足某个收敛准则时，就可以停止计算。例如，我们可以要求基态能量（或感兴趣的最高激发态）在步骤 $m-1$ 和 $m$ 之间的变化小于某个很小的数值。Lanczos方法通常能够将能量收敛到计算机的数值精度（并且应始终使用双精度）。正如我们将在下文看到的，其他性质的收敛速度比能量慢。

期望值。首先，考虑计算算符期望值（除了能量之外）的情况，但程序中没有存储所有Lanczos基向量。那么，为了能够将基态从Lanczos基变换回原始基（以便在态上进行"测量"），我们必须从第一次运行中相同的初始态出发，重复Lanczos基的构造。如果我们在初始化时使用了随机 $\phi_{0}$（如代码17所示），并且没有保存该态，我们可以通过使用与第一次构造时相同的种子初始化随机数生成器来重新生成它。无论如何，有了 $\phi_{0}$ 之后，我们可以像代码段18那样继续，但有一个小的简化：我们现在不必再计算系数 $a_{m} , n_{m}$，因为它们已经存在了。我们可以使用 $n_{m}$ 进行归一化，而不是调用18中的子程序normalize。在基的重新构造过程中，我们使用三对角矩阵对角化得到的基态向量 $\nu_{0}$（或对于激发态i，使用 $\nu_{n}$）来变换这些态，从而在一个向量中逐步构建出基态：

do m = 2,Λ - 1

call hoperation(1, 2)

enddo

至此，向量中的元素包含了基态波函数在所选基矢（通过hoperation实现）下的系数。如果我们存储了所有Lanczos向量，则无需此步骤，可直接对代码19中生成的状态进行变换。

计算仅依赖于自旋z分量且在基矢所用所有对称操作下均不变的算符期望值尤为简单。我们只需用代表态（representative states）对应的波函数系数平方（即该态的概率）对计算量加权即可。例如，所有距离$r = 0 , \ldots , N / 2$下的自旋-自旋关联函数$\langle S_{i}^{z} S_{i + r}^{z} \rangle$可通过如下方式得到：

do a = 1,M
  do r = 0, N/2
    do i = 0, N-1
      j = mod(i+r, N)
      if (s_a[i] == s_a[j]) then
        C(r) = C(r) + psi**2(a)
      else
        C(r) = C(r) - psi**2(a)
      endif
    enddo
  enddo
enddo
C = C / (4*N)

这里$\psi^{2} ( a )$是态$| a \rangle$的概率，当计入对称性后，该态由数值$s_{a}$表示。与前文相同，我们可以通过位测试确定自旋的相对取向，此处是相距r个格点的自旋。需要注意的是，尽管最终的自旋关联$\langle S_{i}^{z} S_{j}^{z} \rangle$仅依赖于$r = | i - j |$，但我们仍需对所有i进行上述平均，因为仅使用代表态（未经对称算符作用时，它们本身并不满足任何晶格对称性）。原则上，我们还应当对代表态的反射和自旋反转进行平均，但在完成平移平均后，自旋关联函数对于其他两种对称性也是不变的（但并非所有算符都如此）。另外注意，由于海森堡哈密顿量具有自旋旋转不变性，此处计算的关联函数$\langle S_{i}^{z} S_{i + r}^{z} \rangle$等于$\langle \mathbf{S}_{i} \cdot \mathbf{S}_{i + r} \rangle / 3$。

计算无法简单关联到z方向对角元的非对角算符期望值时，需要根据式(193)和(193)对态进行显式的操作，然后计算标量积。这原则上很容易——基本上与前面章节中构造哈密顿量的步骤相同——但比对角算符更耗时。

哈密顿量的紧凑存储。现在我们讨论Lanczos程序中使用的子程序hoperation(phi, gamma)的内部机制。它实现了如下操作：

这里正是所有对称性细节体现的地方。与完全对角化中构造哈密顿量的主要区别在于，我们不想将 $H$ 存储为完整的 $M \times M$ 矩阵，因为我们考虑的是基态数目 $M$ 可能高达数百万的计算（例如，表1中 $N=32$ 的情况）。因此，我们必须设计一种便捷的方式，只存储其非零矩阵元，其数量约为 $(195)$ 式中的 $NM$。

我们为紧凑的哈密顿量存储定义以下数据结构：对于 $(195)$ 式中的每个 $a = 1 , \ldots , M$，我们存储非零矩阵元 $\langle b | H | a \rangle$ 的个数 $e_{a}$。我们将这些非零矩阵元的位置 $b$ 作为连续整数存储在一个向量中，其元素为 $B ( i )$。在程序中，对于 $a=1$，方便使用 $i = 1 , \ldots , e_{1}$ 的元素；接着对于 $a=2$，使用 $i = e_{1} + 1 , \ldots , e_{1} + e_{2}$ 的元素，以此类推。向量 $B$ 所需的大小（非零元素个数）一开始是未知的，并且根据所使用的编程语言，可能需要用估计的大小来分配。我们可以将相应的非零矩阵元存储在一个浮点值向量 $H ( i )$ 中。我们还可以利用哈密顿量是对称矩阵这一事实，只存储其一个“三角”部分。

有了上述记号，执行运算 $(195)$ 现在变得非常简单：

subroutine hoperation $( \phi , \gamma )$ = 0; i = 0 {22} do a = 1,M do j = 1,e_a i = i + 1 (B(i)) = (B(i)) + H(i) (a) (a) = (a) + H(i) (B(i)) enddo enddo

计数器 $i$ 跟踪数据结构 $H$ 和 $B$ 中元素的位置。对角矩阵元，即 $B ( i ) = a$，在此过程中会被重复计数，因为哈密顿量矩阵的上三角和下三角的贡献都被加上了。因此，假设它们存储的值已被除以2。

矩阵元 $H ( i )$ 及其位置 $B ( i )$ 的生成方式与我们构造完整哈密顿量时类似，例如，如同使用具有宇称的半动量基中代码片段 $\{13 \} {-} \{16 \}$ 那样。我们还讨论了涉及自旋反转对称性的简单扩展。现在我们讨论当我们只想将非零元素加载到紧凑存储中时所需的修改。

在我们此处考虑的半动量基的情况下，由于同一个代表元可以在态列表中一次或两次出现 $( \sigma = \pm 1 )$，并且我们希望同时处理它们以避免不必要的重复任务，因此流程再次有些复杂。因此，我们遍历基态 $a$，并像代码13中那样确定相同代表元的个数 $n$。我们不能直接将矩阵元放入存储向量 $H$ 中，因为代码22的方案要求每一列的所有矩阵元连续存储。由于在代码 $\{16 \}$ 中，哈密顿量中的各个键算符是逐一处理的，当 $n=2$ 时，如果我们在生成每个矩阵元时直接存储，那么这两列就会混淆。因此，我们为当前正在处理的一列或两列使用临时存储，然后将其内容复制到最终存储的适当位置。将第一个 $(a)$ 列和可能的第二个 $( a + 1 )$ 列分别称为 $c = 1$ 和 $c = 2$，我们暂时将数据存储为 $B_{c} ( k )$ 和 $H_{c} ( k ) , k = 1 , \ldots , n_{c} .$，其中 $n_{c}$ 是该列的矩阵元个数。对角元（如上所述除以2）首先被输入到这些临时存储列表中；

对于非对角矩阵元，我们采用代码段16中的做法。但需要注意的是，当作用于给定基态 $| a \rangle$ 时，多个键操作可能导向同一个态 $\left| b \right.$ ，然而在紧凑存储中，我们不想将同一矩阵元的各个贡献单独存储。通过将代码16中 i 和 j 的循环替换为以下扩展版本，所有问题都得到了处理：

do i = a,a + n do j = max(b,i),b + m c = i a + 1 E = helement(i, j, l, q) do k = 1, n_c if B_c(k) = j exit enddo if (k > n ) then n = k; B (k) = j endif if ( j = i) then H_c(k) = H_c(k) + E/2 else H_c(k) = H_c(k) + E endif enddo enddo

{24}

这里，第二行中的 max(i, b) 确保只构造矩阵三角形 $j \geq i$ 中的元素。最内层关于 k 的循环用于检查是否已有存储的贡献对应于矩阵元 $\langle j \vert H \vert i \rangle$ 。我们假设当关于 k 的循环完成时（即未在 $k = n_{c}$ 之前退出），k 取值为 $n_{c} + 1$（这在许多计算机语言中都是如此）。因此，如果列表中不存在之前的位置 j，则 $k = n_{c} + 1$，这意味着应将其添加到列表中（此时列表大小为 k）。我们还假设所有 $H_{c} ( k )$ 初始设为零，这样每个新贡献 E 都可以加到 $H_{c}$ 的相应位置上。注意，虽然我们在此执行的是非对角操作，但在所使用的基组中，此类操作也可能产生对角矩阵元，此时需要将 E 除以 2。

在生成了当前代表态的所有矩阵元（即完成所有最近邻自旋对的循环）之后，我们将临时存储向量的内容复制到永久全存储中：

这里$n_{H}$是一个计数器，用于记录目前为止添加的非零矩阵元素总数。之后，如代码13所示，关于态标签$a$的循环结束。这样就完成了

图30. 24位点海森堡链在对称性扇区$k = 0 , p = 1 , z = 1$中，四个最低能级的能量（左）和总自旋（右）随Lanczos基尺寸收敛的示例。自旋S是通过假设$\langle \mathbf{S}^{2} \rangle \overset{} {=} S ( S + 1 )$提取得到的。

哈密顿量的构建。

对于其他算符，存储矩阵元素并不划算，因为我们通常不会多次重复使用可观测量算符。然而，哈密顿量会被反复使用，多达几百次。因此，使用代码段21，而不是每次都从头执行与提取每个哈密顿量矩阵元素相关的所有操作，可以带来非常显著的加速。

通过不将矩阵元素存储为双精度浮点数，而是基于整数表将其映射到实际数值，可以进一步压缩哈密顿量。唯一值的数量通常非常少（几十或几百个值，因此甚至可以用“短”整数作为指向实际值的指针），这可以节省一些内存，但代价是构建哈密顿量的时间会更长。然而，位置$B ( i )$必须存储为四字节整数。

当需要处理最大可模拟的系统尺寸的极限时，可能不得不将哈密顿量存储在磁盘上（并在执行代码22时依次读取其各部分），或者在不存储的情况下即时生成。后者本质上相当于每次用哈密顿量作用时，都基于代码24来执行操作。

4.2.4 Lanczos计算的收敛性#

Lanczos方法在使用足够多的Lanczos向量Λ时是基本精确的，通常这个数字并不需要非常大，大约在几十到几百的量级。应通过针对几个不同的Λ进行计算来检查收敛性，直到所关心的能量和期望值不再发生变化为止。基态收敛最快，能量比期望值收敛更快。例如，图30显示了一条N = 24链的结果。图中展示了在基态对称性扇区内，四个最低能级的能量和总自旋随Λ的变化。自旋量子数是根据方程(193)和(194)，通过对状态作用总自旋平方算符S²计算得出的，并假设S² = S(S + 1)来提取S（该假设仅在状态已收敛为该算符的本征态时有效）。可以看出，能量是单调收敛的，但其他量不一定如此，本例中S=2的状态就清晰地表明了这一点。收敛的细节当然取决于构建Lanczos基的初始态（本例中是一个随机态）。此例中，当Λ = 60时，所显示的四个能级（包括能量E和自旋S）都已收敛到小数点后十位以上，而基态在Λ = 30时就已经达到了这一收敛水平。当系统尺寸增大时，收敛会变慢一些，但对于这个特定模型，在可实际研究的最大系统尺寸下，收敛几个能级并无困难。

图31. (a) 16位点海森堡链在量子数 $( k = 0 , p = 1 , z = 1 )$ 下，最低四个能量随Lanczos基组大小的变化。由于正交性丧失，同一态的多份副本依次出现。(b) 对基组进行再正交化后得到的同一系统的五个最低态。

可以通过从接近基态的初始态出发来加速Lanczos计算的收敛。这类初始态可以通过多种方式构建，例如基于某种近似解析方法。但如果没有收敛问题，这样做可能不值得增加额外工作量。然而，如果要对哈密顿量中的某个参数进行一系列计算，后续计算可以从前一参数值的基态出发，该基态很可能与下一个基态有显著重叠。但应注意，如果初始态是基态的良好近似，它与前几个激发态的重叠会非常小，因此在这种计算中只有基态可能快速收敛。如果还需要激发态，可以通过从先前计算中得到的本征态的线性组合出发来启动后续计算，从而规避这一问题。

正交性的丧失与恢复。Lanczos基矢理论上应彼此完全正交，但数值截断误差会逐渐累积，最终导致对某些Λ时正交性的严重丧失。这表现为人为简并，即激发态"跌落"到更低态上。图31(a)展示了这一现象，其中描绘了16位点链在基态对称性扇区中四个最低能量随Λ的变化。可以看到，较高能量依次坍缩到紧邻的较低能量上，在即将坍缩的迹象首次出现后，仅经过几次迭代能级就完全简并。此时基组中包含了相同态的多份拷贝。

正交性的丧失乍看之下可能是一个严重问题，但在实践中却未必如此。基态的计算不会因存在多个副本而受到太大影响，只要确保在计算期望值之前对其进行正确归一化，而第一激发态通常（但并非总是）在降低到较低态之前就已经良好收敛。

小体系的非正交性通常比大体系出现得更早。这与此事实有关：小体系的低能态较少，因此 Lanczos 基组中的最低态可能在该子空间内变得“过完备”（如果 Lanczos 基组比总态数更大，则显然确实是过完备的，计算将无法进行）。

可以通过像方程（188）和代码19中实现的那样，加入显式的再正交化步骤来轻松补充 Lanczos 方法。其缺点是需要存储所有 Lanczos 基向量，这对于大型系统可能不可行。如果态的数量很大，该过程也可能变得耗时。如图31(b)所示，再正交化能够使更多的本征态收敛（仅受内存和时间限制）。

4.3 一维态与量子相变#

在本节中，我们讨论几个相关的基于Lanczos方法的1D计算以及最多至$N = 32$个自旋的有限尺寸标度方法。首先，我们研究标准海森堡链的临界基态，并分析激发态的一些性质。然后，我们引入阻挫，研究$\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$链中的二聚化相变。最后，我们在阻挫之外加入长程相互作用，在这种情况下，连续的二聚化相变演变为奈尔态（如果相互作用足够长程，即使在1D系统中也可能存在）与二聚化VBS态之间的一级相变。

4.3.1 海森堡链的基态与激发#

在这一部分中，我们进一步考察海森堡链的基态以及激发态性质。

尽管海森堡链存在精确的贝特假说解[51]，但其波函数非常复杂，在许多情况下，必须利用有限尺寸系统的数值计算来从中提取物理性质的信息[53, 54]。然而，对于有限链和热力学极限，基态能量[52]和低能激发[149]都可以精确计算。其他一些量也可以从非常大的链中提取[150]。除了使用精确解外，这类系统（即具有更广泛相互作用并保持系统对称性的系统）的许多性质都是基于渐近精确的低能场论描述而得知的；即具有拓扑耦合（或中心荷）c = 1 的韦斯-朱米诺-维滕非线性模型[56]。等价地，由于自旋旋转不变性，海森堡链也是描述一维自旋、费米子和玻色子相互作用系统的卢廷格液体态的一个特例[151, 152]。基于精确对角化和其他无偏有限格点技术（如QMC和DMRG）的数值研究，在指导和验证这些理论方面发挥了重要作用（例如，参见文献[153]和[7]）。要提取普适连续场论所无法描述的非普适（短程、高能）性质（包括实验可直接测量的量，例如文献[154, 155]中所讨论的），也需要进行计算研究。海森堡链的数值结果对于在研究其他自旋模型时，能够对推断技术进行严格的基准测试也极其重要，而这些模型的分析计算结果则知之甚少。

图32. 在自旋反转量子数$z = 1$（给出自旋$S = 0$或$2$的低能态）和$z = -1$（最低态总是给出$S = 1$）的扇区中，一个$N = 32$链相对于基态的最低激发色散关系。数据点下方的实线是des Cloizeaux和Person[149]利用贝特假说导出的精确无限大尺寸色散关系。

基态与低能激发。 对于大小为$N = 4 n$的海森堡链，其基态具有动量$k = 0$、宇称$p = 1$和自旋反转数$z = 1$；它是一个完全对称的单重态。对于$N = 4 n + 2$，最低能态反而位于完全反对称扇区；$k = \pi , p = -1 , z = -1 , S = 0$。

让我们来看海森堡链最重要的激发。对于不等于基态动量的$k$，有限系统的最低能态是三重态；因此，对于大小为$N = 4 n$和$4 n + 2$的链，$z$分别等于$-1$和$+1$。图32显示了一个32位点链的最低$z = +1$和$z = -1$态激发能$\omega ( k ) = E s ( k ) - E_{0}$的动量依赖性。最低的$z = +1$态总是具有$S = 1$，而最低的$z = -1$态则具有$S = 0$或$S = 2$。三重态能量非常接近由贝特假说[149]导出的精确$N = \infty$三重态色散关系$\omega_{S} ( k ) = \pi | \sin ( k ) | / 2$，特别是在接近$k = 0$处。在接近$k = \pi$时，偏差较大。$S = 0, 2$的激发能量通常略高，但接近$k / \pi = 0.85$的$S = 2$态实际上略低于三重态（可能是一个小$N$异常），而在$k / \pi = 2 / N$处则高得多。

对于大的$N \to \infty$，可以预期单重态和三重态激发将退化为简并态，其色散关系为$\omega_{0, 1} ( k ) = \pi | \sin ( k ) | / 2$。这可以理解为源于称为自旋子（spinons）的自旋1/2孤子状自由度的激发对。这些自旋子相互作用非常弱，在无限长链中表现为独立粒子，因此形成四个简并能级（由此可构成$S = 0$和1的态）。图32中绘制的色散关系只是自旋子激发连续谱的最低边缘，该连续谱可以使用贝特假说[53]进行详细计算，并且已在准一维反铁磁体[156]中通过实验观测到。

图33. 海森堡链的单重态能隙和三重态能隙随系统尺寸倒数的变化。插图显示了能隙乘以N后的结果，表明主导的$\sim 1/N$标度（动力学指数z = 1）存在乘性对数修正。

对于$N = 4n$链，最低激发是动量$k = \pi$处的三重态，其$p = -1, z = -1$；而对于$4n+2$链，最低三重态出现在$k = 0, p = 1, z = 1$处（即相对于基态的动量差始终为$\pi$，而p和z等于基态对应值的负值）。图33展示了这些最低三重态和单重态激发能量（单重态和三重态有限尺寸能隙）随系统尺寸倒数变化的有限尺寸标度。两者均表现出按$1/N$标度趋于零，对应动力学指数$z = 1$。这种标度形式存在源于对数修正的弱修正，在将能隙乘以N后（图33插图）可以清晰看到。基于连续场论方法[57, 153]已对此能隙标度做出详细预测。能隙标度的前因子直接与自旋子激发速度相关，但由于对数修正的存在，难以从有限尺寸数据中可靠提取该速度。速度也可以从接近$k = 0$和$\pi$处色散关系的线性部分（即$\omega_{S}(k)/k$或$\omega_{S}(k)/(\pi - k)$）提取，其中$k = 0$处的行为更易分析，如图32清晰所示。

自旋关联函数。 自旋关联函数是任何量子自旋系统的关键特征之一。基态中的关联是主要关注对象。对于自旋各向同性和平移不变的一维链，我们可以将其写为

对于任何参考位置 i。正如我们在第 2 节中讨论的，默明-瓦格纳定理排除了海森堡链中的反铁磁长程序。相反，在基态中，对于大 r，$C ( \boldsymbol r )$ 以 $( - 1 )^{r} / r$ 形式衰减，并附有一个乘性对数修正（源于场论描述中的边缘算符[57, 58, 59]）。

图34. 海森堡链基态中自旋关联函数绝对值乘以距离 r 的结果。展示了不同系统尺寸下作为距离函数的结果，以及在 r = N/2 处随链长 N 变化的结果。

图 34 展示了系统尺寸 $N = 16, 24$ 和 32 下 $C ( \boldsymbol r )$ 作为 r 函数的基态结果。结果已乘以 r，因此如果渐近形式为 $\sim 1 / r$，则对于大 r 应观察到与 r 无关的行为。同时乘以 $( - 1 )^{r}$ 也很有用，以抵消交错相位的振荡。注意 $C ( r ) ( - 1 )^{r}$ 中仍存在奇偶振荡。这种振荡在一维系统的各种关联函数中相当常见。它们随着系统尺寸的增大而减弱，但在本文研究的较小系统中仍然很强。Luttinger 液体理论也预测了这种振荡的形式；它们以 $1 / r^{2}$ 衰减[152]。

在图 34 中，可以观察到 $C ( r )$ 相对于 $1 / r$ 形式的明显偏离，这主要源于周期性边界条件，它增强了接近 $r = N / 2$ 处的关联。开放边界会导致更严重的有限尺寸效应（此外，还无法利用平移对称性进行块对角化）。对数修正也应该起一定作用。随着链长的增加，固定 r 处的关联函数收敛于无限 N 极限形式。然而，如图所示，结果近似收敛的最大 r 仅占 N 的很小一部分。对于 $N = 32$，可以可靠地说 $C ( \boldsymbol r )$ 仅对 $r \approx 4$ 以内的 N 收敛，这不足以说明长程行为。图中还展示了一种替代方法：研究固定 $r / N$ 处作为 N 函数的关联函数，其中 $r / N = 1 / 2$ 是检验长程行为最自然的选择。尽管周期性边界此时显著增强了关联，但 $C ( r = N / 2 )$ 的函数形式正比于无限尺寸的 $C ( \boldsymbol r )$ [尽管 $C ( r = N / 2 )$ 相对于 N 的整体前因子显然不同于实际 $N \to \infty$ 收敛的 $C ( \boldsymbol{r} )$]。基于图 34 中的数据，主导的 $1 / r$ 形式看起来相当合理，而剩余的增强与预测的对数修正一致。在第 5.3.1 节中，我们将基于对更长链的量子蒙特卡洛计算，进一步研究长程关联。

4.3.2 阻挫驱动的量子相变#

$S = 1 / 2$ 海森堡链的一个非常有趣的方面是，当在哈密顿量中加入次近邻相互作用时，它表现出多种基态相变（量子相变）；

本章新增术语：

该模型被称为$J_{1} {-} J_{2}$链，或马朱姆达-高希模型（以其首次全面研究该系统的作者命名[60]）。在此我们将讨论最重要（也最被充分理解）的相变，该相变发生在$J_{1}$和$J_{2}$均为反铁磁（正值）时。为方便起见，我们定义比值$g = J_{2} / J_{1}$。

对于$g < g_{c} , g_{c} \approx 0.2411$ [157]，系统处于与上一节讨论的$J_{2} = 0$海森堡链相同的相。基态是临界的，反铁磁关联按$1 / r$衰减，有限尺寸能隙按$1 / N$标度（两者均存在对数修正）。只有前因子（例如激发速度）和对数修正的强度作为g的函数发生变化。

对于$g > g_{c}$，系统处于一种完全不同的状态，具有指数衰减的自旋关联和一个在热力学极限下保持有限的三重激发能隙。基态具有长程二聚体序（是第2.2节中讨论的二维VBS态的一个更简单的一维版本）。最近邻键强度（$J_{1}$能量贡献）$\left. B_{i} \right. = \left. \mathbf{S}_{i} \cdot \mathbf{S}_{i + 1} \right.$在对称性破缺的无穷大尺寸态中呈现交替（周期二）形式$B_{i} = B + \delta ( - 1 )^{i}$（存在两个简并态，其振荡彼此反相）。与磁序不同，这种序在一维零温下是允许的，因为破缺的对称性（晶格平移对称性）是离散的。调制$\delta$在$g = g_{c}$处变为非零，此后随$g$增大而增加。自旋关联最初是交错的（在倒易空间中动量$k = \pi$处有峰值），但在$g \approx 0.52$处变为$k = \pi / 2 [ 158 ]$，并且对于$g > 1$，峰值可能随$g$连续变化（螺旋态）[159]。

具有铁磁（负值）$J_{1}$的模型也很有趣。在具有不同周期性的态之间存在若干相变[160]。这里我们将只考虑反铁磁情形，并且位于$g < 1$的区域。我们将使用兰佐斯方法研究进入VBS态的相变，并讨论这种有序态的性质。

马朱姆达-高希点。在讨论兰佐斯结果之前，值得指出的是，VBS序的存在可以在特殊点$g = 1 / 2$（马朱姆达-高希点）处精确证明，在该点基态非常简单[60]。在具有偶数N的环上，基态是双重简并的单重态乘积态，如图35所示。通过将哈密顿量作用于这些态，可以相当容易地证明它们是本征态（而证明它们是最低态则更为复杂[161, 162, 163]）。在这些破缺平移对称性的简并态$\left| \Psi_{A} \right.$和$\left| {\Psi_{B}} \right.$中，可以形成对称和反对称态，它们分别具有动量$k = 0$和$k = \pi$：

量子自旋系统的计算研究

图35 马朱姆达-高希点上简并基态的示意图。粗弧线表示单重态；交替键上这些单重态的乘积是$g = 1 / 2$处的精确基态。

这些是动量守恒的Lanczos计算中得到的态。因此，在有限周期系统上计算时，序参量$\langle B_{i} \rangle = E_{1}$是无特征的。与精确有限格点计算中一贯的情况相同，必须通过计算关联函数来观察对称性破缺（除非我们添加某种微扰来打破两种可能有序模式之间的对称性，但这本身也有其复杂性，因为需要取微扰消失的极限）。离开$g = 1 / 2$点后的基态比图35中简单的单重态乘积更为复杂，单重态配对存在涨落，但最近邻键上单重态的交替高低密度仍然存在。有限环上平移不变的$k = 0$态大致对应于如式(198)中的对称和反对称组合。然而，只有在$g = 1 / 2$处，这两个态对于有限$N$才是精确简并的。在其他情况下，只有在$N \to \infty$的极限下，这些态才会变得简并。

二聚体序参量。我们已经在2.2节讨论过VBS基态，并在2.4.3节讨论过一个如何通过二聚体（或键）关联来表征此类序的例子。现在研究式(25)中引入的二聚体关联函数，对于一维系统，其显式形式为

如果存在VBS序（即使对称性没有被破坏），对于大的r，此函数应当在“弱”值和“强”值之间交替。利用旋转对称性可以简化用Lanczos方法计算该关联函数的过程，并计算 $\langle ( \mathbf{S}_{i} \cdot \mathbf{S}_{i + 1} ) S_{i + r}^{z} S_{i + 1 + r}^{z} \rangle$，它是 $D ( r )$ 的 $1 / 3$。原则上，我们可以通过只使用z分量来定义序参量，使其更加简单：$\langle S_{i}^{z} S_{i + 1}^{z} S_{i + r}^{z} S_{i + 1 + r}^{z} \rangle$，该量不仅仅是(199)式的常数倍，但对于二聚化态仍然是一个有效的序参量。这里我们考虑完整的 $D ( r )$。

图36(a)展示了三种不同耦合下32位点体系的结果。在 $g = 0$ 时，即标准的海森堡模型，关联随r迅速衰减，并在 $r = N / 2$ 时非常小。在转变点 $g_{c} \approx 0.2411$（我们将在下面更精确地确定）处，关联明显更强。应用于自旋链的Luttinger液体理论预言，对于 $g \le g_{c}$，二聚体关联会以 $1 / r$ 的形式衰减（当 $g < g_{c}$ 时有更强的对数修正）[152]。当耦合比更大时，图36(a)显示，在 $g = 0.4$ 处的关联清楚地表明了长程序的存在。为了确认长程有序的存在，必须进行有限尺寸标度分析。就像我们在前一节中对海森堡链的自旋关联函数所做的那样，通常最好的做法是观察最长距离 $r = N / 2$ 处的关联随体系尺寸的变化。对于二聚体关联，即使是在无序相中，“裸”关联函数也不为零，因此必须减去对应于这个本底的常数。一种可能性是减去 $\langle B_{i} \rangle^{2}$，这是 $D ( r )$ 对于大r应收敛到的值。我们在本文中选择的另一种方法是采用差值 $D ( N / 2 ) - D ( N / 2 - 1 )$。只有当 $D ( r )$ 振荡时，即存在长程二聚体序时，该量在 $N \to \infty$ 时才会非零。图36(b)展示了三种不同体系尺寸的结果随耦合比的变化。在此图中，序参量在 $0.4 < g < 0.6$ 范围内似乎已很好地收敛，证实了长程有序确实存在。

图36. (a) N = 32 链在三种不同耦合比 g 下的键关联函数。(b) 三种不同体系尺寸的键关联函数随 g 的变化。实线和虚线分别表示临界耦合 g_c 和精确可解点 g = 1/2

自旋和二聚体关联的有限尺寸外推。为了定量地提取无限大体系的二聚体序参量，在大多数情况下有必要进行外推，除非体系尺寸足够大以至于不再有任何显著的尺寸依赖性[基于图36(b)，在靠近 $g = 1 / 2$ 时就是这种情况]。图37展示了代表性的g值（包括VBS相内部和外部）下，二聚体和自旋关联随体系尺寸倒数变化的结果。

图37右图展示的二聚体关联，在 $g < g_{c}$ 时具有 $1 / r$ 形式，并存在对数修正，但恰在 $g_{c}$ 处这些修正很小。对于 $g > g_{c}$，即VBS相内部，关联外推至非零值，此时渐近的 $N \to \infty$ 收敛应为指数形式（可在 $g = 1 / 2$ 处明确证明），但在靠近相变点时，实际上无法达到足以观测到此行为的体系尺寸。相反，越接近 $g_{c}$，行为在 $1/N$ 上本质上是线性的，如图中 $g=0.4$ 所示。然而，在 $g=0.45$ 处，可以看到行为实际上是非单调的，其大N行为与 $g=1/2$ 时类似，表现为从下方趋近于无穷大N值。如图36右图所示，当 $g > 0.6$ 时，二聚体序参量随 $g$ 增大而减小。在这种情况下，需要更大的体系才能进行恰当的外推，因此必须采用其他方法。人们相信，对于所有 $g > g_{c}$，体系都保持二聚化状态[158, 159]。

当系统进入VBS相时，自旋关联应从 $1 / r$ 形式变为指数衰减。为探测这一变化，关联函数 $C ( r = N / 2 )$ 在左图中乘以了 $N / 2$。如之前所见，在 $g = 0$ 时，存在一个对数修正，增强了长程自旋关联。在 $g = g_{c}$ 处，这些修正小得多——已知在二聚化相变点，主导的对数修正会消失[164]——尽管行为仍然不是纯粹的 $1 / r$ 形式（这只有在 $N \to \infty$ 极限下才能预期）。对于VBS相内部的 $g$，$r C ( r )$ 如预期那样衰减至零。在 $g_{c}$ 附近，为了观测到纯粹的指数衰减，当然必须采用大体系尺寸——要大于自旋关联长度（该关联长度在 $g_{c}$ 处发散）。

图37. 周期性链中最长距离处的自旋关联（左图）和键关联（右图）随体系尺寸倒数的变化。五条曲线的耦合比自上而下（左图）和自下而上（右图）依次为 $g = 0 , g_{c} , 0.40 , 0.45 , 0.50$，即用黑色圆圈显示的三组数据对应VBS相。

对于图37中使用的g值，自旋关联是交错的，即在倒易空间中峰值位于 $k = \pi$。当 $g \approx 0.52$ 时，峰位迅速变为 $\pi / 2$ [158]；对于 $g > 1$，自旋结构很可能演变为具有连续变化螺距的螺旋态[159]（如同经典版本模型的情况）。然而，自旋关联始终是指数衰减的（不同于经典的长程螺旋态）。

确定二聚化相变点。已知二聚化相变类似于经典二维XY模型的Kosterlitz-Thouless相变。与常规相变不同，序参量在 $g_{c}$ 处不遵循幂律，而是在 $g_{c}$ 附近呈指数小。因此，基于Lanczos计算所能达到的小体系，不可能提取出 $g_{c}$ 附近无穷大体系的序参量——由于交叉行为，即使在VBS相内部相当深处，这也十分困难。这样一来，确定二聚化相变的位置似乎非常困难。然而，有一种非常巧妙的方法，可以基于激发态能量间接地提取临界耦合[157]。

正如我们在第4.3.1节所见，纯海森堡链（$ g = 0 $）的最低激发态是一个三重态。另一方面，我们也知道在VBS态中，基态应该是二重简并的，并且这两个基态都必须是单重态。在精确可解点 $g = 1 / 2$ 处，对于任意N，简并性都是精确的，但偏离这个特殊点时，这两个态仅在N趋于无穷大的极限下才变得简并。简并性的趋近方式随N呈指数变化（下文的数据将说明这一点），而海森堡链的三重态激发能隙随 $1/N$ 趋于闭合，最低单重态也以 $1/N$ 的方式趋近基态（如图33所示）。既然如此，对于固定的N，必然存在某个耦合比 $g_{\mathrm{cross}} ( N )$，在该处单重态和三重态激发发生交叉。这个点可以作为依赖于系统尺寸的相变点定义，并且当 $N \to \infty$ 时，它应趋近于实际临界耦合；$g_{\mathrm{cross}} ( N \to \infty ) = g_{c}$

图38. (a) 16个格点链的最低三个能量作为耦合比$g$的函数。$k = \pi$单重态与三重态能级的交叉点$g_{cross}$可用作依赖于系统尺寸的临界（二聚化）耦合$g_{c}$的定义。竖线表示精确可解的Majumdar-Ghosh点，在该点单重态是简并的。(b) 交叉点关于$1/N^{2}$的图。曲线是对$1/N$的多项式拟合（无$\propto 1/N$项，因为主导修正项为$\propto 1/N^{2}$）

图38(a)展示了16个格点链的三个相关能级，并标出了相应的量子数（这些量子数适用于所有系统尺寸 $N = 4 n$）。注意两个单重态如何随着 $g \to 1/2$ 而变得简并。正如预期，单重态和三重态激发能级发生交叉；对于该系统尺寸，交叉点位于 $g \approx 0.242$。如图38(b)所示，交叉点作为 $N$ 的函数收敛得非常快——已知主导修正项与 $1/N^{2}$ 成正比[164]。基于这些结果，可以非常精确地确定 $g_{c}$。对 $8 \leq N \leq 32$ 范围内的交叉点进行高阶多项式拟合，得到 $g_{c} = 0.2411674 ( 2 )$，其中 (2) 表示由于在拟合中包含不同的系统尺寸范围以及变化多项式阶数时，外推值的波动所导致的最后一位数字的不确定性。

VBS相中的有限尺寸能隙。如前所述，为了在热力学极限下容纳对称性破缺，在VBS相内部，随着系统尺寸增大，单重态-单重态有限尺寸能隙将以指数方式快速闭合（Majumdar-Ghosh点除外，因为在该点简并性对任何尺寸都是精确的）。图39左图展示了对几个代表性的$g$值，能隙的指数标度行为。在$g_{c}$处，行为变为 $\propto 1/N$，这在此线性-对数图中并不明显，但在右图的对数-对数图中则更为清晰。右图还展示了三重态能隙随$1/N$的变化。在此可明显看出，在VBS相内部，外推值非零。这个有限能隙本质上对应于将一个单重态键激发到一个三重态所需的能量。注意，只有在Majumdar-Ghosh点，基态才是精确可解的，而三重态能隙在该点也依赖于尺寸。在$g_{c}$处，能隙按 $1/N$ 标度——图中展示了在依赖于尺寸的交叉点$g_{\text{cross}}$处的行为，它也按 $1/N$ 标度。

图39. 不同耦合比下，单重态（左图）和三重态（右图）能隙随系统尺寸的变化。这里$g_{\mathrm{cross}}$是与尺寸相关的耦合值，在此耦合下最低的单重态和三重态激发是简并的（因此$\Delta_{s} = \Delta_{t} )$。在$g_{\mathrm{cross}}$处的能隙按$1 / N$衰减（在右图中更清晰可见）。左图中$\Delta_{s}$相对于$N$采用了对数标度，以显示在有序的VBS相中，随着$N \to \infty$而变得简并的两个单重态之间的能隙呈指数级迅速闭合。右侧，三重态能隙相对于$1 / N$绘制，以显示其在VBS相中收敛至非零值。

4.3.3 具有长程相互作用的链#

虽然具有有限程相互作用的海森堡链中不可能实现长程自旋序，但长程相互作用使得在$T = 0$时能够形成磁序。在一个具有距离依赖耦合$J_{r} \propto ( - 1 )^{r - 1} / r^{\alpha}$[72]的系统中，会发生奈尔态与准长程有序（QLRO）态（即前面几节讨论的自旋关联按$1 / r$衰减的海森堡临界态）之间的相变。这里的符号有利于反铁磁序，并且不存在阻挫。当指数$\alpha < \alpha_{c} \approx 2$时，基态具有真正的长程奈耳序，而当$\alpha > \alpha_{c}$时，系统处于QLRO相。长程相互作用参数$\alpha_{c}$的临界值取决于耦合的细节，例如当所有其他$J_{r}$固定时最近邻耦合$J_{1}$的值，并且相变处的临界指数是连续变化的（这与QLRO相内恒定的指数形成对比）。

另一个具有长程相互作用的系统的有趣例子是Haldane-Shastry链[165, 166]，它具有阻挫相互作用$J_{r} = 1 / r^{2}$。它的临界基态类似于标准海森堡链，但在场论语言中，导致主要对数修正的边缘算符消失了[167]。因此，该系统恰好处于一个二聚化相变点，类似于上面讨论的$\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$链的情况。可以注意到，在Haldane-Shastry模型中$J_{2} / J_{1} = 1 / 4$，这非常接近$\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$链的临界比率$g_{c} = 0.2411$。因此，Haldane-Shastry链中距离$r = 2$以外的相互作用只起次要作用（但重要的是，它们实际上使该系统是精确可解的[165]）。

图40. 具有阻挫参数g和长程相互作用指数$\alpha$的链的基态相图，取自文献[73]。虚线表示连续相变，而粗实线表示一级相变。对于$g < 0$的细实线对应于文献[72]中研究的无阻挫相互作用的系统。在${\boldsymbol{\alpha}}^{- 1} = 0$（即$\mathbf{J}_{1} \mathbf{-} \mathbf{J}_{2}$链）处，VBS态中主要的自旋关联从$q = \pi$变为$\pi / 2$（在$g \approx 0.52$处）[158]。对于$\alpha^{- 1} > 0$，发生这种变化的曲线连接至所有相边界汇合的多临界点。

在存在长程相互作用的情况下，也可以实现奈尔态与VBS之间的一维直接量子相变。这里我们讨论文献[73]中引入的系统，它根据哈密顿量将无阻挫的长程相互作用和短程阻挫（在间隔$r = 2$处）结合起来：

其中耦合的距离依赖性由下式给出：

在此，通过将所有耦合（除 $J_2$）用 $J_{\Sigma}$ 归一化后，所有非阻挫相互作用 $\left| J_{r} \right|$ 之和为1，因此极限 $\alpha \to \infty$ 恰好对应具有 $g = J_2 / J_1$ 的 $J_1$-$J_2$ 链。这在比较两个模型时将很有用，并保证对 $N \to \infty$ 每个自旋的能量有限，即使对于 $\alpha < 1$ 也是如此。也可以将求和范围从 $r = N / 2$ 扩展到 $N / 2 < r < N$。然而，这不应影响 $\alpha > 1$ 时的相界和临界指数。

我们将研究当长程相互作用指数 α 的倒数 $\alpha^{-1}$ 从 0 增加时，作为阻挫强度 g 的函数发生的二聚化相变的演化。如图 40 中的半定量相图所示（基于 Lanczos 结果构建，下文将讨论），这个连续相变持续到 $\alpha \approx 2$，而对于更小的 α，它演变为奈尔态与一个既具有 VBS 序又在波数 $q = \pi / 2$ 处具有临界（或可能长程）自旋关联的状态之间的一级相变。该状态在相图中记作 $\mathrm{VBS} {+} \mathrm{QLRO} ( \pi / 2 )$。

二聚化相变的演化。我们已讨论过如何基于 $J_1$-$J_2$ 链的激发态能级交叉来研究二聚化相变，如图 38(a) 所示。如果 α 足够大，相同原理也适用于存在长程相互作用的情况。这在图 41 左图所示的 $\alpha = 4$ 的 16 自旋链中得以体现。对于 $N \to \infty$，在 VBS 相中，最低的 $k = 0$ 和 $k = \pi$ 单重态应该变得简并（这样才能形成对称性破缺的二聚化态）。图中可以看到 $g > 1 / 2$ 的区域存在一个非常接近简并的区域。这个近似简并区域（难以精确划定）随着 N 增加而非常缓慢地向更小的 g 扩展。如图 42 所示，单重态-三重态交叉点随系统尺寸收敛 [尽管比图 38(b) 中 $J_1$-$J_2$ 链显示的 $1 / N^{2}$ 收敛更慢]，并且可以可靠地给出 $\alpha > 2$ 时的 $\mathrm{QLRO} ( \pi ) {\mathrm{- V B S}}$ 相界 $g_{c} ( \alpha )$。

图 41. 长程指数 α 取两个值时，16 自旋系统的低能级；$\alpha = 4$（左）和 $\alpha = 1$（右）。左图标注了各态的自旋 S 和动量 k（注意，显示了两个具有 $S = 0 , k = 0$ 的能级）。右图插图示出了两个 $k = 0$ 单重态之间回避交叉的更多细节。这两个单重态在靠近时存在一个非常小的能隙，但远小于本图可分辨的尺寸。

能级的另一个有趣特征是，当 α 降至 $\approx 2$ 以下时，基态能量随 g 变化的宽极大值变得越来越尖锐。如图 41 右图所示，在 $\alpha = 1$ 时，由于与第二个 $k = 0$ 单重态发生了回避能级交叉，这个极大值已发展成一个尖锐尖顶。真实的单重态-三重态交叉也移到了同一区域。两个具有相同量子数的态之间的回避能级交叉，导致对于 $N \to \infty$ 基态能量对 $g$ 的导数出现不连续，这是一级相变的标志。需要注意的是，并非只有最低的两个单重态表现出这种回避能级交叉。其他低能态也在相变附近区域聚集，形成一个复杂的能级交叉级联——当 $N \to \infty$ 时，所有这些都应收敛到一级相变处的同一个点。该相变处各相的性质将在下文讨论。首先，让我们更详细地研究二聚化相变如何从连续演变为一级相变。

图42展示了能级交叉点$g_{\mathrm{cross}}$和基态能量最大值位置$g_{\mathrm{peak}}$对系统尺寸的依赖关系。在标准 Majumdar-Ghosh 阻挫链中，交叉点的尺寸修正$\propto 1 / N^{2}$ [如图38(b)所示]，这在大的$\alpha$值下也能观察到。对于较小的$\alpha$值，修正似乎为$\propto 1 / N$，但只要相变保持连续，对于大的 N，似乎很可能交叉到 $1 / N^{2}$ 行为。峰值位置则向相反方向移动。对于某些情况，当 $N \to \infty$ 时，$g_{\mathrm{cross}}$ 和 $g_{\mathrm{peak}}$ 应当重合。结果表明，在该点，$g_{\mathrm{cross}}$ 和 $g_{\mathrm{peak}}$ 的主导修正项都是$1 \mathord{\left/ {\vphantom{1 \ N}} \right. \kern - delimiterspace} N$。图42中展示了$\alpha = 2$时的线性拟合，此时外推得到的两个值之间仍存在一个小的间隙（如果线性标度仅在$g_{\mathrm{cross}} = g_{\mathrm{peak}}$这一特殊参数值下渐近成立，那么外推可能不完全准确）。在$\alpha = 1.7$处，相变是一级的，它们应该重合（此时渐近尺寸修正应是指数形式的）。

图42. 对于不同的长程相互作用指数$\alpha$，单重态-三重态交叉点$g_{\mathrm{cross}}$和基态能量最大值位置$g_{\mathrm{peak}}$对逆链长的依赖关系。两条线显示了将$g = 2.0$的数据外推至$N = \infty$的结果。

一级 Néel-VBS 相变。为了确认在$\alpha > 1.8$（近似）处存在具有不连续能量导数的回避能级交叉，我们在图43中以半对数坐标绘制了基态能量在其最大值处的二阶导数。这些结果基于在接近峰值的一个非常密集的网格点上计算能量，并可以可靠地拟合出多项式。对于$\alpha = 1.5$，从该多项式提取的二阶导数随 N 呈指数增长，这表明对于无限长链，能量曲线的斜率确实发生了不连续的变化。相反，在$\alpha = 3$时，二阶导数对大 N 减小。对于$\alpha = 2$，收敛到一个有限值似乎也是合理的，而$\alpha = 1.7$和1.8则似乎接近分隔这两种不同行为的临界线（在此临界线上，发散的形式与幂律一致）。该分析表明，连续的二聚化相变在$( g_{m} \approx 0.41 , \alpha_{m} \approx 1.8 )$处平滑地转变为一级相变。单重态-三重态交叉点向基态能量最大值移动，并在多临界点$\left( g_{m} , \alpha_{m} \right)$与之重合，超过该点后，它发展成一级奇异性。

为了分析参与一级相变的状态，接下来我们研究标准的自旋关联函数$C ( \boldsymbol r )$和在(199)中定义的二聚体关联函数$D ( r )$。我们还将研究这些关联函数的傅里叶变换；即静态自旋 (s) 和二聚体 (d) 结构因子：

图43. 在基态能量取最大值$g_{\mathrm{peak}}$处，基态能量关于阻挫参数 g 的二阶导数对系统尺寸的依赖关系。

如果在某个波矢 $q=Q$ 处存在长程序，则相应的结构因子与 N（对于大N）成正比，因此序参量可以定义为 $S_{s,d}(Q)/N$。这里，对于奈尔态和二聚化VBS态，$Q=\pi$。图44展示了在 $\alpha=1.5$ 时，$S_d(\pi)/N$ 和 $S_s(\pi)/N$ 对阻挫参数g的依赖关系。可以看到，在 $g\approx 0.42$ 处出现不连续性，这与一级相变情景一致（而在 $g=1$ 处，此处未显示，不连续性要尖锐得多）。奈尔序参量和VBS序参量在其各自相中都良好收敛到非零值，而在另一相中则很小（在该相中，当 $N\to\infty$ 时，它们应趋于零，但由于系统尺寸较小，无法清晰地观察到这一点）。图中还显示了 $q=\pi/2$ 处的自旋结构因子。有趣的是，它在VBS相中也变得很大，而在奈尔相中，随着N增大，它应衰减到零（基于此数据这是合理的，但同样，对于这些小系统无法清晰地看到）。基于这些结果，可以得出结论：在一级相变处获得的VBS态在 $k=\pi/2$ 处共存有长程磁序（即周期为四的螺旋序）。然而，考虑到系统尺寸较小，这并不能保证。这种行为也可能反映了自旋关联的缓慢幂律衰减。

图45展示了 $\alpha=1$ 时，在一级相变两侧的两个g值处的实空间自旋和二聚体关联。在 $g1$，几乎看不到随r的衰减）。在此耦合强度下，$D(r)$ 中没有结构，即没有VBS序。对于 $g>g_c$，显然存在VBS序，$D(r)$ 显示出特征性的交错图案。同时还观察到了强的周期四自旋关联，但尚不清楚这些关联是长程的，还是随着r增加而缓慢衰减到零。这种行为与 $1/r$ 衰减一致，但需要更长的链才能可靠地提取其行为。这里需要注意的一个有趣点是，尽管没有减去常数，但二聚体关联围绕零振荡。这与图36中看到的 $\mathbf{J}_1\mathbf{-}\mathbf{J}_2$ 链中围绕正值振荡的行为形成对比。

图44. 具有长程相互作用指数 $\alpha=1.5$ 的系统中，基于尺寸归一化的二聚体结构因子 $S_d(\pi)$（左图）和自旋结构因子 $S_s(\pi)$（中图）及 $S_s(\pi/2)$（右图）定义的VBS和磁序参量随阻挫强度g的变化。随着系统尺寸增大，在一级相变处形成不连续性。

回顾式(199)中二聚体关联函数的定义，这种行为是周期四自旋结构的另一种表现，即在键算符 $\mathbf{S}_i+\mathbf{S}_{i+1}$ 的奇分离处，测量的是铁磁键与反铁磁键之间的关联。虽然该VBS中磁结构的性质尚不完全清楚，但关联函数似乎表明存在一种键序倾向，形式为 $\left| s t s t s t \cdots \right.$，其中s和t分别代表具有高单重态密度和三重态密度的最近邻键。在图40中，该态被标记为磁性准长程序态，即VBS–QLRO（$\pi/2$）共存态，但事实上它也可能是VBS与长程磁序共存的态。还应指出，在理论上，该相在小晶格上可能对边界条件非常敏感，并且不能排除实际磁结构可能是螺旋状的，间距连续变化（如 $g>1$ 的 $\mathbf{J}_1\mathbf{-}\mathbf{J}_2$ 链中所示[159]）。

$\mathrm{VBS} {+} \mathrm{QLRO} ( \pi / 2 )$ 态应具有无能隙的自旋激发，无论自旋关联如何，只要它们不是指数衰减的。最低的三重态位于 $k = \pi / 2$。然而，基于小系统的数据很难证明其无能隙性质，因为系统尺寸 $N = 4 n$ 的能隙（及其他量）的尺寸依赖性表现出较大的奇偶 n 振荡（以及其他不规则尺寸效应）。在传统的 VBS 相中（图 40 相图的右下部分），最低的三重态位于 $k = \pi$，即使自旋关联（在此相中呈指数衰减）的峰值位于 $k = \pi / 2$。因此，原则上，最低的 $k = \pi$ 和 $k = \pi / 2$ 三重态之间的能级交叉可用于提取 VBS 和 $\mathsf{VBS} \mathrm{+ Q L R O} ( \pi / 2 )$ 相之间的边界。然而，交叉点的尺寸依赖性并不平滑，因此无法非常可靠地外推。占据主导的 $k = \pi$ 和 $k = \pi / 2$ 自旋关联之间的边界也未被精确提取。这种自旋关联的变化可能与向周期为四的态 [158] 的转变有关，尽管没有迹象表明 VBS 序发生了变化。

QLRO–Néel 相变。让我们回到图 41，关注能级谱的另一个有趣特征：对于小的 $g$，当 $\alpha = 4$ 时，最低的单重态激发具有动量 $k = \pi$，但当 $\alpha = 1$ 时，其动量为 $k = 0$。这些能级在 $g < g_{c}$ 时作为函数发生的顺序交换与 Néel–QLRO( ) 相变有关。能级交叉可用于非常精确地提取此相边界，直至 $g \approx 0.25$（而对于更高的 $g$，$N \to \infty$ 的外推变得困难）。虽然不完全清楚为什么 $k = 0$ 和 $k = \pi$ 激发单重态在此相变处发生交叉，但可以再次指出，在 $\mathbf{Q} \mathbf{L} \mathbf{R} \mathbf{O} ( \pi )$ 相中，低能激发来自两个退禁闭的自旋子。最低的（对于 $N = 4 n$）是位于 $k = \pi$ 的三重态。最低的单重态也位于 $k = \pi$，由于弱自旋子-自旋子相互作用，它与三重态之间存在一个小的有限尺寸能隙（图 33 展示了这些能级如何随 $N \to \infty$ 而简并的示例）。在 Néel 相中，自旋子不再退禁闭，低能谱的结构随之改变。

图 45. 在 $\alpha = 1$ 的 32 自旋链中，两个阻挫参数值下的自旋（左）和二聚体（右）关联。在 $g = 0.25$ 和 0.45 时，系统分别处于 Néel 和 $\scriptstyle \mathrm{V B S - Q L R O} ( \pi / 2 )$ 相。此情况下的一级相变发生在 $g \approx 0.39$。

正如我们将在下面第 4.4.2 节中讨论的那样，有限晶格上的二维 Néel 态具有类似于量子转子态的低能激发，这些激发的 $S = 1 , 2 , 3 , \dots$，并随着系统尺寸的增加与基态简并。对于偶数 S 和奇数 S，它们的动量分别为 $k = ( 0 , 0 )$ 和 $k = ( \pi , \pi )$。人们可能在一维 Néel 态中也期望这种量子转子"塔"。显然，不属于这种转子塔的最低 $k = \pi$ 单重态，其能量被推高很多，而 $k = 0$ 单重态取代它成为最低的单重态激发。

将能级交叉方法与文献[72]中无阻挫模型（$J_{1} = 1$，$J_{r > 1} = \lambda ( - 1 )^{r - 1} / r^{\alpha}$）的QMC计算结果进行比较是很有用的。将该模型（$\lambda = 1$）重新参数化，转换为(201)中使用的约定，得到了图40中$g < 0$的曲线。对Néel序参量的QMC数据进行有限尺寸标度，对于N最大至4096的情况，给出了$\alpha_{c} = 2.225 \pm 0.025$[72]。对于$N \leq 32$的情况，外推$k = 0 , \pi$单重态交叉点（此处其尺寸修正$\propto 1 / N^{\beta}$，$\beta \approx 1.50$）给出了一个略高（且可能更可靠）的值：$\alpha_{c} = 2.262 \pm 0.001$。分析交叉点处的单重态和三重态能隙，假设$\Delta \sim N^{- z}$，得到了动力学临界指数$z = 0.764 \pm 0.005$，与文献[72]的结果吻合良好。对于QLRO( )–Néel边界上的阻挫模型(201)，$1/N$的能隙标度在$g \mathrm{up} \approx 0.25$时给出$z \approx 0.75$，而对于更大的$g$值，由于大的标度修正以及随着多临界点$( g_{m} \approx 0.41 , \alpha_{m} \approx 1.8 )$的接近（作为参考，当$g = 0$时$\alpha_{c} = 2.220 {\pm} 0.005$，当$g = 0.2$时$\alpha_{c} = 2.170 {\pm} 0.01$），对于越来越大的系统而言不存在能级交叉，因此无法通过这种方式可靠地提取$\alpha_{c}$和$z$。

在多临界点处，可以利用单重态-三重态交叉点处的能隙和基态能量最大值来提取$z$，得到$z \approx 0.8$。因此，动力学临界指数在整个QLRO( )–Néel边界上（包括所有相边界交汇的多临界点）可能是一个常数，$z \approx 3 / 4$。

技术说明。从技术角度来看，我们在此研究的模型对于一维系统而言颇具挑战性，原因在于其长程相互作用。长程相互作用与阻挫的结合使得应用其他计算技术变得困难。尽管高效的QMC技术可以应用于无阻挫的长程相互作用系统[168, 72]，但由于符号问题（在第5节中讨论），在存在阻挫的$J_{2}$项时，这就不再可行了。另一方面，DMRG方法[28, 29]可以处理阻挫，但不易处理长程相互作用。尝试将DMRG，或基于矩阵乘积态[21]的相关技术应用于该系统，尤其是为了进一步研究假想VBS-QLRO(/2)态中的自旋结构，仍然会是很有意义的。

上述Lanczos计算利用了第4.1节中讨论的所有对称性。哈密顿量以紧凑形式存储在磁盘上，对于(201)中的每个相互作用距离$r$都有一个单独的文件。根据代码22，在与哈密顿量的每次操作中，这些文件被逐个读取和使用。对于$N = 32$，所需的Lanczos迭代次数通常为$60 \sim 80$次（对于较小的系统则更少）。所有基态都保存在主存中，每次迭代后都进行与所有先前态的重新正交化（如代码19所示）。

4.4 二维系统#

现在，我们以简单周期方格上的海森堡模型这一重要情况为例，讨论二维系统的精确对角化。与一维情况类似，我们可以利用动量守恒在固定磁化强度的扇区内对哈密顿量进行分块对角化。为了充分利用晶格对称性来进一步分割某些块，情况要稍微复杂一些，因为与哈密顿量对易的对称操作数量更多。从积极的一面来看，实现所有晶格对称性的回报是，对于高对称动量（可以利用晶格反射和/或旋转），对于给定的$N$，二维的块比一维的更小。实际上，这可能帮助不大，因为可以被对角化的$L \times L$晶格的最大线性尺寸$L$仍然小得可怜$( L = 6 )$。精确对角化研究仍然极具价值。

在一维问题中，我们使用了反射（宇称）算符$P$，并基于代表态$| a \rangle$及其反射态$( 1 + p P ) | a \rangle$构造了动量本征态（后来还加入了自旋反转）。虽然$P$不与平移算符对易，但我们发现对于两个特殊动量$k = 0 , \pi$，$k$和$p$仍然是好量子数。对于一般的$k$，$p$不是有效的量子数，但借助混合了不同$k$态的半动量态，我们可以利用$P$来实现实值表示。在二维中，原则上也可以为任意动量构造实值哈密顿量，但这在实践中比一维复杂得多，最终也不值得付出这种努力（因为它不会增加可研究的系统尺寸）$^5$。因此，这里我们将使用标准的复数动量态。

图46 $4 \times 4$周期方晶格上的格点和键的标记。在使用位表示时，习惯上将格点标号从0开始，此处根据$i = x + y L$，其中$x , y \in \{0 , \ldots , L - 1 \}$且$L = 4$。键根据与其相连的“左-下”参考格点$i$来标记；水平键和垂直键分别为$b_{x} = i$和$b_{\mathrm{y}} = L^{2} + i$。

定义晶格。 与一维类似，我们使用态符号$| S_{0}^{z} , \ldots , S_{N - 1}^{z} \rangle$，该符号与自旋的位表示配合使用很方便。对于一维系统，自旋指标直接对应于链的几何结构，但对于二维晶格，我们必须为格点建立一种标记约定。我们将考虑具有$N = L_{x} \times L_{y}$个格点的矩形晶格，其坐标为$( x , y )$，其中$x = 0 , \ldots , L_{x} - 1$且$y = 0 , \ldots , L_{y} - 1$。一种自然的选择是将格点编号为$i = 0 , \ldots , N - 1$，使得$i = x_{i} + y_{i} L_{x}$，如图46中一个$4 \times 4$系统所示。在不显式依赖晶格的符号表示中，我们可以将海森堡哈密顿量写为

其中$[ i ( b ) , j ( b ) ]$是由键$b$连接的两个最近邻格点。具体说明晶格只需创建这个格点对列表。

4.4.1 二维中的动量态#

我们再次定义自旋指标的平移，现在在$x$和$y$两个方向上进行，如图47所示，相应的算符$T_{x}$和$T_{y}$定义为

图47 用于二维方晶格的对称操作。$T_{x}$和$T_{y}$分别将自旋沿$x$和$y$正方向平移一个晶格常数。晶格是周期的；空心圆代表对边的格点。$P_{x}$和$P_{y}$是关于$x$轴和$y$轴的反射，$P_{d}$和$P_{e}$是关于两条对角轴的反射。对于偶数$L$（反铁磁体通常应使用偶数），$P_{x}$和$P_{y}$的轴线穿过晶格格点之间。反射$P_{d}$和$P_{e}$是关于连接远角的线，这些线穿过对角线上的格点（因此在反射下，这些格点上的自旋保持不变）。

其中平移后的自旋指标为

其中 $[ \gamma_{i} - 1 ]_{L_{\gamma}}$ 表示 $\gamma_{i} - 1$ 模 $L_{\gamma}$，即 $[ - 1 ]_{L_{\gamma}} = L_{\gamma} - 1$。

利用这些平移，基于某个代表态 $| a \rangle$ 的动量态定义为

其中可能的动量为

归一化常数 $N_{a}$ 取决于代表态的平移性质，即在 $L_{x} \times L_{y}$ 个平移操作作用于代表态 $| a \rangle$ 后，所得到的不同状态数目 $D_{a}$（例如，对于自旋呈棋盘格图案的状态，$D_{a} = 2$）。若在将 $| a \rangle$ 映射到自身的那些平移上，(207) 式中相位因子之和 $F_{a}$ 为零，则该代表态与动量不相容。对于相容的状态，归一化常数 $N_{a} = D_{a} | F_{a} | ^{2}$。在程序中计算这一点的最简单方法是直接完成所有平移操作，并在过程中累加 $D_{a}$ 和 $F_{a}$，而不是像在一维系统中那样使用显式公式（当然也可以使用更暴力的方法）。

哈密顿量矩阵的构造方式与一维情形相同。我们再次将哈密顿量拆分为对角部分 $H_{0}$ 和非对角键算符 $H_{j}$，如 (121) 和 (122) 式所示，不过此时 $j = 1 , \dots N_{b}$ 。将这些算符作用于代表态 $| a \rangle$ 后，我们可能同样需要对得到的状态进行平移，以获得对应于 $H_{j} | a \rangle$ 的新代表态 $\left| b_{j} \right.$，即

将其作用于动量态，我们得到

并且取与 $\langle b_{j} ( k_{x} , k_{y} ) |$ 的重叠积分，得到矩阵元

归一化常数（或它们可能数值的某些整数映射）应与代表态一起存储。

融入其他晶格对称性。 我们在此考虑一个 $L \times L$ 的正方形系统。这样的晶格有八个独立的旋转和反射，可以取为，例如，四个 $90^{\circ}$ 旋转以及关于水平轴或对角轴的反射。可以很容易地验证，额外的反射可以写成旋转和第一个反射的组合。看出这一点的最简单方法是考虑一个 $2 \times 2$ 阵列的所有旋转和反射：

这里第一行包含了第一阵列的所有旋转，第二行中这些阵列通过交换两列进行了反射。任何其他反射只会产生已属于上述八个集合中的阵列。元素的其他排列不对应于旋转和反射的组合（因此不对应于我们通常感兴趣的哈密顿量的对称性）。实践中，我们可以使用任意方便的反射集，或旋转与反射的组合，以生成这八个唯一的变换。

利用平移以及一组其他对称算符（我们暂时将它们合并为一个算符 $Q$，对应量子数 $\{\boldsymbol{q} \}$），可以定义一个动量态：

如上所述，对称操作的数量多于唯一变换的数量。我们可以选择对于给定动量最方便的那些。这里我们将使用关于 x 轴和 y 轴的反射 $P_{x}$ 和 $P_{y}$，而不是正方晶格的旋转。这些对应任何矩形 $( L_{x} \times L_{y} )$ 晶格的哈密顿量的对称性（注意 $180^{\circ}$ 旋转等价于 $P_{x} P_{y}$）。对于一般动量，反射不与平移对易，因此无法使用 $P_{x}$ 和 $P_{y}$ 进行进一步分块，但存在高对称动量，使得 $P_{x}$ 或 $P_{y}$（或两者）可以使用。对于 $\dot{k}_{x} = \pm k_{y}$，我们将改用关于两条垂直对角轴的反射 $P_{d}$ 和 $P_{e}$。所有反射定义在图 47 中。对于最特殊的动量 $\mathbf{k} = ( 0 , 0 )$ 和 $( \pi , \pi )$，且 $L_{x} , L_{y}$ 任意，我们还可以在 $P_{x} , P_{y}$ 之外使用 $P_{d} \ ( \mathrm{or} P_{e} )$。

为了在给定动量 $( k_{x} , k_{y} )$ 下实现最大块对角化，可以将 (213) 中的算符 $Q$ 选择如下：

(214)

其中所有反射量子数取值为 $\pm 1$。注意情况 (g) 中 $p_{x} = p_{y}$ 的限制，即对于 $p_{x} \neq p_{y}$ 不存在这样的态。在这种情况下，我们当然也可以使用选项 (f)，但使用 (g) 可以进一步将块分成两个。

为了证明上述定义的对称算符 $Q$ 确实在不同动量扇区内具有相应的守恒量子数，以下晶格变换之间的关系很有用：

剩下的对子是互易算符；$[ T_{x} , T_{y} ] = 0$，$[ T_{x} , P_{y} ] = 0$，$[ T_{y} , P_{x} ] = 0$，$[ P_{x} , P_{y} ] = 0$，$[ \bar{P_{e}} , \mathbf{\hat{\Pi}} {}^{} ] = 0$ 。我们只检查最复杂的情况，即(214)中的情况(g)。出于与第4.1.4节讨论的一维系统相同的原因，在这些特殊动量下，置换算符与平移算符对易。因此，我们只需研究$Q$的性质。显然，该态是$P_{d}$（即最左边的反射算符）的本征态。为了检查它是否也是$P_{x}$和$P_{y}$的本征态，我们利用(215)中的一些关系来验证$P_{x} Q = p_{x} Q$和$P_{y} Q = p_{y} Q$。经过一点代数运算，我们得到

以及$P_{y}$的类似结果。因此，如果$p_{x} = p_{y}$，那么对于特殊动量$(0, 0)$和$(\pi, \pi)$，该态确实也是$P_{x}$、$P_{y}$和$P_{d}$的本征态。

如果我们想研究所有动量扇区，利用反射对称性并不会带来太多好处，因为在大多数块中，我们必须在(214)(a)中使用通用选项$Q = 1$。这些对称性在研究基态时很有用，基态通常位于最大可分块扇区(g)中。反射最有用的方面可能仍然是它们能提供对态对称性方面的洞察（例如，用于对激发进行分类）。除了晶格对称性，自旋反转对称性可以在$m_{z} = 0$扇区中使用，这与一维情况完全相同。

我们已经考虑了反射对称性的组合，但我们也可以使用旋转和反射的组合。这将给出不同的量子数，但块的大小相同。在群论中，存在具有标准化名称的不可约表示，对应于各种晶格对称性。对于新情况，我们无需自行推导适用的对称性，而是可以查阅不可约表示及其对应的特征标表，并利用它们来构建具有选定不可约表示量子数的态。这在研究更复杂的晶格时非常有用。文献[143]讨论了这种更形式化的方法，并提供了链和方晶格的例子。

这里我们使用的对称性仅适用于沿方晶格x轴和y轴边缘的$L_{x} \times L_{y}$晶格。当$L_{x} = L_{y} = L$时（如果我们关注趋近无限二维晶格，这是最有意义的情况），我们只能取$L = 4$和$L = 6$（因为奇数$L$会给反铁磁体带来不期望的阻挫效应）。为了增加可访问的晶格尺寸$N$，我们还可以考虑“倾斜”晶格，其边缘不平行于方晶格轴[143]。

实现。为了研究$a \left. 6 \times 6 \right.$晶格，我们必须使用长整数（8字节）进行比特表示。为了节省时间，更好的方式不是遍历所有$2^{N}$个可能的态整数并筛选出那些对应特定磁化强度的态（如代码7中所示），而是采用一种从一开始就只构建具有给定$m_{z}$（即给定$n_{\uparrow}$个自旋向上数）的态的方案。当然，这也可以对一维系统进行。这种更复杂的特定$m_{z}$态的生成从整数$a = 2^{n_{\uparrow} + 1} - 1$开始，其中前$n_{\uparrow}$个比特为1，其余为0，然后执行：

do
call checkstate(a, pass)
if (pass) then M = M + 1; store representative information endif
c = 0
do b = 1,N
if (a[b 1] = 1) then
if (a[b] = 1) then
c = c + 1
else
a[0], . . . , a[c 1] = 1; a[c], . . . , a[b 1] = 0; a[b] = 1; exit
endif
endif
enddo
if (b = N) exit
enddo

表3. 磁化强度 $m_{z} = 0$ 的 $L \times L$ 格子 $( L = 4 , 6 )$ 上 $k = ( 0 , 0 )$ 态区块的尺寸 $M ( L )$，考虑了不同反射和自旋反转量子数。对于 $p_{x} = - p_{y}$ 的情况，不存在量子数 $p_{d}$。对于 $p_{x} = - 1 , p_{y} = 1$，区块结构与 $p_{x} = 1 , p_{y} = - 1$ 相同。

这里的内循环简单地搜索最低的比特位置 $b - 1$，在该位置，a 的一个置位（1）比特可以向左移动一步（即移到当前比特为0的位置b）。找到这样的位置后，所有低于该位置的之前已置位的比特（其数量由计数器c跟踪）被移到最低的位置 $( 0 , \ldots , c - 1 )$。这种比特演化正好对应于基2里程表（odometer）的数字从右向左进位的方式。checkstate 的内容取决于使用了哪些对称性，但与我们在1D情况下讨论的实现非常相似，除了在 (205) 中定义的平移比简单的循环比特置换更复杂，需要手动实现。

在构建哈密顿量矩阵元时所需的归一化常数，应由 checkstate 提供（如果一个代表态通过了测试）。计算归一化的最简单方法（而不是像1D情况那样使用形式化表达式）仍然是，对代表态执行所有对称操作，并为每一个将代表态映射到自身的对称操作，累加 (213) 中因子（包括复数动量相位以及来自反射量子数的正负号）的总和 $F_{a}$。该数值与 $| a \rangle$ 的非等价变换数量 $D_{a}$ 一起，给出了归一化 $\dot{N_{a}} = D_{a} | F_{a} | ^{2}$（如果代表态与量子数不兼容，则同样有 $N_{a} = 0$）。我们不再讨论如何实现基生成和哈密顿量矩阵构造的进一步细节，因为这些任务是对第4.1.3和4.1.4节讨论的一维实现的直接推广。$m_{z} = 0$ 扇区中的自旋反转对称性可以按照第4.1.5节讨论的方式实现。

区块尺寸示例。表3列出了Lanczos计算中关注的方格子系统在其中一个特殊动量 $\mathbf{k} = ( 0 , 0 )$ 下的区块尺寸，在该动量下可以利用最大数量的对称性。海森堡模型的基态位于此区块中（在所有其他量子数都等于1的子区块中）。对于 $L = 6$，即使最小的区块也包含超过1500万个态，而最大的区块（对于一般k）则大约大8倍。这种尺寸的区块仍然可以在Lanczos计算中处理（较小的区块在标准工作站上相当容易处理）。

为了测试基的构建，检查子区块尺寸之和是有用的，该和应等于区块中的总态数。对于 $L = 4$，$\mathbf{k} = ( 0 , 0 )$ 有822个态，这等于表3中所有反射区块尺寸之和。

4.4.2 奈尔态及其量子转子激发#

现在我们用海森堡模型的二维Lanczos计算结果加以说明。尽管这种方法可处理的小晶格不足以对热力学极限进行定量精确的外推，但这些计算确实展示了具有奈耳序的系统的一些重要方面（超出了我们在第2.1节自旋波理论框架中讨论的内容）。我们引入有限系统低能态的量子转子映射，并基于此讨论磁化率。我们还计算了子晶格磁化强度。

量子转子态的双自旋模型。如果存在反铁磁序，则子晶格A上的自旋主要取向相同，而子晶格B上的自旋主要取向与A相反。然而，如果自旋数$N$有限且对称性未被破坏，则由子晶格A自旋定义的整体方向并非固定。这种情况可以通过考虑单个子晶格A和B的自旋和来刻画[169]：

作为两个固定长度的自旋$S_{A} = S_{B} = N / 2$（更精确地说，对于大$N$，我们应写为$S_{A} = S_{B} = m_{s} N$，但常数无关紧要），如图48所示。假设这两个大自旋以最简单的旋转不变方式反铁磁耦合，即通过有效海森堡相互作用：

其中$\mathbf{S} = \mathbf{S}_{A} + \mathbf{S}_{B}$是总自旋。这里$\mathbf{S}_{A}^{2}$和$\mathbf{S}_{B}^{2}$只是与$N^{2}$成正比的常数，在讨论激发能量时可以忽略。然而，这些常数意味着耦合常数$J_{A B}$应$\propto 1 / N$，以使总基态能量$\propto N$。(218)的基态具有总自旋$S = 0$，激发态具有$S = 1, 2, \ldots$，能量比基态高出$J_{A B} S ( S + 1 ) / 2$。

我们定义一个新的与$N$无关的耦合$J_{\mathrm{eff}} = N J_{A B} / 2$，并将能量写为

这些激发态被称为量子转子态塔（tower of quantum rotor states）。当$N \to \infty$时，$S \ll \sqrt{N}$的态简并，它们可以组合成具有固定奈尔矢量方向的基态（类似于量子力学中局域化一个粒子所需的无限多动量态），从而实现对称性破缺，而对称性破缺正是自旋波理论的出发点。在热力学极限下，有序矢量的方向是固定的（因为与其旋转相关的时间尺度发散[169]），因此量子转子态实际上不会被涉及。它们在标准的自旋波计算（在第2.1节中讨论）中从一开始就被忽略，因为通过构造，序被锁定在z方向。通过考虑处于外磁场中的系统，可以将其调节到具有总磁化$S^{z} = S$的基态[170, 171]，从而在自旋波理论中仍能获取转子能量。在有限团簇中，转子态具有重要意义。

图48. 用两个大自旋$\mathbf{S}_{A} , \mathbf{S}_{B} ,$（对应于两个子晶格上的自旋总和）对旋转不变的奈尔矢量$\mathbf{m}_{s}$进行的有效描述。这些自旋之间存在有效的反铁磁耦合，导致了一个单重态基态和一个总自旋为$S = 1 , 2 , \dots$的量子转子激发态"塔"，其能量比基态高出$\mathbf{\bar{\Delta}} \Delta_{S} \sim S ( \bar{S} + 1 \bar{]} / N$。

对于给定的系统，如果我们能将(219)中的有效耦合$J_{\mathrm{eff}}$与某个依赖于转子激发的物理量联系起来，就可以确定它。一个明显的选择是均匀磁化率$\chi = d \langle m_{z} \rangle / d h$。在双自旋模型中计算它，当$T \to 0$时给出$\chi = 3 / J_{\mathrm{eff}}$。对于维度$d \geq 2$的有限团簇上的真实海森堡模型，当$T \ll 1 / L$时，应该由量子转子态主导，因为最低自旋波能量标度为$\propto 1 / L$（而量子转子态标度为$1 / L^{d}$）。因此，我们可以将具有奈耳基态的海森堡模型的有效量子转子态塔写为

其中$\chi$应在(先令)$N \to \infty$然后$T \to 0$的极限下评估。注意，这里的$I = ( 3 / 2 ) N \chi$扮演了转动惯量的角色，使得(220)与量子力学中刚性转子的能谱形成类比。

关系式(220)也可用作数值计算海森堡模型磁化率的一种方法；通过提取作为S函数的最低能量（对于小$S$，应适用量子转子映射）。更精确地说，小S能量给出了$\chi$的估计值，作为$N \to \infty , S \to 0$极限下量$\chi ( S , N )$的值

这里$E_{S}$表示总自旋为S的最低能量。注意，我们必须减去基态能量$( S = 0 )$，因为在两自旋有效模型中，我们只计算了相对于基态能量的激发能$\Delta_{S}$（而双自旋模型并不能精确给出基态能量）。只有当$S \ll L$时，我们才能预期观察到与S无关的行为，因为更高的转子态会受到双自旋模型未考虑效应的影响。

图49. 在$4 \times 4$和$6 \times 6$格点上海森堡模型的Lanczos结果。(a) 每个总自旋扇区最低态相对于基态的能量，经$N S ( S + 1 )$归一化，因此完美的量子转子激发应产生与S和N无关的值。(b) 周期性格点上所有可能距离r的自旋关联与距离的关系。虚线标出了已知的$r \to \infty$值（来自图5的QMC结果）。

与第4.1.6节和第4.3.1节讨论的低能一维量子数类似，在二维正方晶格上（其中所有偶数L有$N = 4 n$），量子转子态分别对应奇数S和偶数S的动量$(\pi , \pi)$和(0,0)。对于偶数S，最低态在(214)中的反射量子数为$p_{x} = p_{y} = p_{d} = 1$，而对于奇数S，相应的量子数为$p_{x} = p_{y} = - 1 , p_{d} = 1$。对于给定的S，最低态可以在磁化扇区$m_{z} = S$中获得（由于$m_{z} \neq 0$，这里无法使用自旋反转对称性）。$(S=0)$基态处于全对称扇区；动量为(0, 0)且$p_{x} = p_{y} = p_{d} = z = 1$。

$L = 4$和$L = 6$的Lanczos结果如图49(a)所示。对于(220)显然存在较大的修正，因为$\chi^{- 1} ( S , N )$随S显著减小并随N增大。对于固定的N，$S = 1$和$S = N / 2$之间的差异约为10%。考虑到$S \ll L$的极限实际上无法基于小晶格进行研究，这种量级的偏差并不令人惊讶。转子态已通过QMC计算在更大的晶格上进行了研究[171, 33]。对于小S和大N，最精确的计算给出$\chi^{- 1} ( S , \bar{N} ) \to 22.8$ [33]。

为了理解与转子图像的偏差，可以利用轻微非刚性量子转子的类比，这似乎很自然，因为双自旋模型是用固定长度的自旋定义的，而在实际系统中，子晶格自旋显然会波动（其方式可能取决于总自旋S）。通过在相互作用(218)中添加高阶项$( \bar{\bf S}_{A} \cdot{\bf S}_{B} )^{2}$等，或许可以考虑这些效应。然而，这种扩展双自旋图像的细节尚未得到研究。

横向磁化率（Transverse susceptibility）。对于$T = 0$的无限系统，自旋旋转对称性被破坏。然后我们可以考虑磁化率相对于奈尔矢量的横向和纵向分量$\chi_{\perp}$和$\chi_{\parallel}$。在双自旋模型中，显然$\chi_{\parallel} = 0$，这在海森堡模型中同样成立。在对称性未被破坏的大系统中，各分量被旋转平均，给出$\chi = ( 2 / 3 ) \chi_{\perp}$。在实际具有奈尔序的准二维系统（例如高温铜氧化物）中，由于各向异性和/或三维耦合，系统可以在$T > 0$时具有有限的有序温度，因此横向磁化率在实践中是相关的。不过，与有序相关的能量尺度很小，并且由孤立二维平面（以海森堡反铁磁体为模型）得到的$\chi_{\perp} = ( 3 / 2 ) \chi$在$T \to 0$的值，为这个量在更复杂系统中的估计提供了良好的参考。

子晶格磁化强度。首次尝试计算二维海森堡模型子晶格磁化强度的工作由 Oitmaa 和 Betts 于 1978 年使用 Lanczos 计算完成 [172]。当时，$6 \times 6$ 格点仍无法触及，基于更小格点大小的外推得到了过大的值。图49展示了在$4 \times 4$和$6 \times 6$格点上，所有可能间隔$r$下的自旋关联函数。该函数随$r$缓慢衰减，并且结果至少看起来合理地在趋近于一个非零值$\langle m_{s}^{2} \rangle > 0$，当$r \to \infty$（且$N \to \infty$）时。QMC 计算的结果（已在图5中展示）给出了图49中用虚线标出的$r \to \infty$结果。事实上，二维海森堡模型中长程有序是否存在一直存在争议，直到 20 世纪 80 年代末第一批可靠的 QMC 模拟完成 [173]。

可以尝试将 Lanczos 结果外推到$N \to \infty$，例如用于阻挫系统，在这些系统中其他计算极具挑战性 [98, 99]。然而，普遍来说必须对这些结果持谨慎态度，因为它假定了在非常小的格点上表现出的行为已经反映了最终的渐近行为。基于不同形状的“Betts 簇” [174] 的拟合过程可能会产生误导，即使尺寸依赖性是平滑的。存在一些例子，由于系统中存在某个有限长度尺度，当$L$超过该尺度后才能观察到渐近行为，因此在较大距离处发生了交叉行为 [82]。