第3章 · 经典相变、蒙特卡洛模拟与有限尺寸标度#

量子相变的许多方面,以及如何基于数值有限格点数据对其进行分析,与经典相变非常相似。这里我们在经典相变的更简单背景下讨论这一通用形式体系,然后再转向量子系统的计算和数据分析。

在经典统计力学中,具有连续相变的系统原型是二维伊辛铁磁体,其精确解 [111, 112] 严格证明了这种相变的存在。该模型定义为

$$ E_{\sigma} = - J \sum_{\left. i j \right.} \sigma_{i} \sigma_{j} - h \sum_{i = 1}^{N} \sigma_{i} ,\tag{28} $$

这只是作为自旋 $\sigma_{i} = \pm 1$ 函数的(势)能。我们用 $\pmb{\sigma} = ( \sigma_{1} , \ldots , \sigma_{N} )$ 来整体表示整个自旋构型。若相互作用键 ij 限制在简单二维正方晶格上的最近邻,无穷大系统的精确临界温度为 $T_{c} / J = 2 / \ln ( 1 + \sqrt{2} ) \approx 2.269$ 。序参量为磁化强度 $\langle m \rangle = \langle \sigma_{i} \rangle$ 。它在 $T < T_{c}$ 且 $T \to T_{c}$ 时的渐近形式为 $m \propto | t | ^{\beta}$ ,其中 t 是约化温度 $t = ( T - T_{c} ) / T_{c}$ ,指数 $\beta = 1 / 8$ 。此类临界指数以及连续相变中标度行为的其他方面将在本节讨论。

伊辛模型的精确解非常特殊,通常人们以其他方式研究相变。临界指数在简单的平均场理论中已经出现,但其值通常不正确。然而平均场理论作为起点是必不可少的,这里我们针对伊辛模型进行概述。相变最重要的理论框架是重正化群,它解释了普适的(仅依赖于对称性和维度,而非相互作用细节)非平凡指数(即不同于一般平均场值的指数)是如何产生的(例如,参见 Cardy 的著作 [15])。为了以无偏方式计算临界指数和其他性质,通常使用蒙特卡洛模拟,其中自旋构型根据玻尔兹曼分布进行随机采样。$^3$ 本节概述蒙特卡洛模拟的基础知识(更多细节参见 Landau 和 Binder [115] 以及 Newman 和 Barkema [116] 的著作)。为了普遍地激发和说明有限格点计算,我们还将讨论如何从蒙特卡洛数据中提取热力学极限下的性质。

本节大纲。第3.1节讨论伊辛模型的标准平均场处理。第3.2节介绍蒙特卡洛方法,并用其说明对于大系统(形式上在无穷大系统大小的极限下),相变和对称性破缺如何实际发生。第3.3节回顾热力学极限下临界行为的关键方面(定义指数等),然后讨论有限尺寸标度假说,并通过分析二维伊辛模型的蒙特卡洛结果来演示其用法。第3.4节以一个阻挫二维伊辛模型为例,讨论一级相变以及用于探测不连续性的有限尺寸标度方法。第3.5节在XY(平面矢量)自旋模型的背景下引入自旋刚度(对应于固体的弹性模量)这一重要概念。通过三维和二维XY模型的蒙特卡洛结果说明自旋刚度的标度性质,其中后者不表现出进入有序态的正常相变,而是表现出进入低温临界态的不同的Kosterlitz-Thouless相变。第3.6节简要讨论如何将经典临界性概念推广到量子相变。

3.1 伊辛模型的平均场理论#


在平均场理论中,无限晶格子系统的环境被替换为一个外场,该外场代表子系统与环境之间的平均相互作用。子系统可以是单个自旋或由多个自旋组成的团簇。这里我们只考虑伊辛模型中最简单的单自旋计算情况,即单个自旋的无限环境被一个有效场取代。

为了方便,将哈密顿量(28)写成扩展形式,包含所有自旋之间(在任意晶格上)的一般相互作用$J_{i j}$,同时也包含一个外磁场:

$$ E_{\sigma} = - \frac{1} {2} \sum_{i = 1}^{N} \sum_{j = 1}^{N} J_{i j} \sigma_{i} \sigma_{j} - h \sum_{i = 1}^{N} \sigma_{i} .\tag{29} $$

注意因子$1 / 2$,这是因为求和中的每一对相互作用都被计算了两次(且$J_{i i} = 0$)。我们不对$J_{i j}$的符号和大小施加任何限制,但为简单起见,假设如果存在相变,有序态是铁磁性的。为了构造和论证单自旋近似,我们首先将(29)中涉及任意自旋$i$的所有相互作用归并在一起:

$$ E_{i} = - \sigma_{i} \left( \sum_{j} J_{i j} \sigma_{j} + h \right) .\tag{30} $$

注意这里求和前面没有因子$\frac{1} {\lambda}$,并且我们利用了$J_{i j} = J_{j i}$。

现在假设$h > 0$,因此$m = \langle \sigma_{i} \rangle > 0$。我们将在无场的情况下研究自发有序,最终令$h \to 0$。通过加减一个常数,可以将(30)中括号内的项写为

$$ \sum_{j} J_{i j} \sigma_{j} + h = m \sum_{j} J_{i j} + h + \sum_{j} J_{i j} ( \sigma_{j} - m ) .\tag{31} $$

在其最基本的表述中,平均场理论相当于忽略(31)中的第二个求和——即涨落项——之后我们就剩下一个易于求解的问题:单个自旋处于强度为$J_{s} m + h$的有效磁场中;

$$ E_{\sigma} = - ( J_{s} m + h ) \sigma ,\tag{32} $$

其中$J_{s}$是原始耦合的总和

$$ J_{s} = \sum_{j} J_{i j} ,\tag{33} $$

并且我们假设系统具有平移不变性,因此该总和与$i$无关。(32)中的磁化强度$m$在此阶段未知,将通过自洽条件$\langle \sigma \rangle = m$来确定。

当前表述下的平均场理论可以被证明是合理的,前提是所忽略的涨落平均而言远小于(31)中的其他项,即满足

$$ \delta_{m} = \frac{\left. \left| \sum_{j} J_{i j} ( \sigma_{j} - m ) \right| \right.} {J_{s} m + h} \ll 1 .\tag{34} $$

这里的 $\langle \rangle$ 表示自旋实际概率分布(玻尔兹曼分布)下的期望值,我们无法精确计算。原则上,涨落也可以在平均场理论框架内计算,作为一种内部自洽性检验。即使不做任何计算,我们也能粗略推断出 $\delta_{m}$ 较小的条件。显然,当存在显著有序时,即 $( 1 - m ) = ( 1 - \langle \sigma_{j} \rangle ) \ll 1$,$\delta_{m}$ 很小,因为此时大多数 $\sigma_{j} = 1 \approx m$。当 $h$ 较大时就是这种情况。对于 $h = 0$ 的情况,如果对称性自发破缺且 $m$ 接近 1(即存在相变时 $T \ll T_{c}$ 的情况),同样成立。此外,即使 $m$ 不是很大,只要对 $j$ 的求和涉及许多非零(且相对较大)的耦合常数 $J_{i j}$,$\delta_{m}$ 也可能很小——由于涨落的抵消,求和 $\textstyle \sum_{j} J_{i j} \sigma_{j}$ 通常(在大多数统计上重要的自旋构型中)会接近 $m \textstyle \sum_{j} J_{i j}$ 。在无限程均匀相互作用的极端情况下,对所有 $i , j$ 有 $J_{i j} = J / N$,且 $N \to \infty$(我们将 $J$ 视为有限常数,例如 $J = 1$,以得到有限能量),所有涨落完全抵消,$\delta_{m} = 0$(此时平均场理论成为精确理论)。这一结论对于无限维晶格上的短程相互作用同样严格成立。因此,一般而言,即使 $m$ 不大,对于高维系统和/或长程相互作用系统,可以预期涨落度量 $\delta_{m}$ 很小。这些正是预期平均场理论能够定量准确的适用条件。即使在定量不准确的情况下,平均场理论依然能在定性上提供有价值的见解。

现在让我们实际求解平均场问题 (32),即自洽条件 $\langle \sigma \rangle = m$ 下的单自旋问题 (32)。磁化强度为

$$ \langle \sigma \rangle = \frac{\sum_{\sigma} \sigma \mathrm{e}^{\sigma ( J_{s} m + h ) / T}} {\sum_{\sigma} \mathrm{e}^{\sigma ( J_{s} m + h ) / T}} = \mathrm{tanh} [ ( J_{s} m + h ) / T ] ,\tag{35} $$

因此自洽条件为

$$ m = \operatorname{tanh} [ ( J_{s} m + h ) / T ] .\tag{36} $$

该方程通常需要数值求解,可以通过对解进行连续区间定位来轻松完成。对于小的 $m$ 和 $h$,我们可以通过展开到 $m , h$ 的领头阶来进行解析处理。首先,当外场 $h = 0$ 时,$m = 0$ 是所有温度下的一个解。为了寻找其他可能的解,将 (36) 式展开到 $x = J_{s} m + h$ 的三阶项,tanh $( x ) = x - x^{3} / 3$,我们有

$$ m^{2} = 3 \frac{T^{2}} {J_{s}^{2}} \frac{J_{s} - T} {J_{s}} ,\tag{37} $$

图10. 伊辛模型的平均场解。(a) 温度高于、等于和低于 $T_{c}$ 时磁化强度与外场的关系。$h \to 0^{+}$ 和 $h \to 0^{-}$ 时的不连续性对应于 $h = 0$ 时的自发磁化强度。(b) 自发磁化强度的温度依赖关系。

从中我们可以识别出临界温度 $T_{c} = J_{s}$,低于该温度时磁化强度可以非零。$T \ll T_{c}$ 时的渐近行为为

$$ m = \left( 3 \frac{T_{c} - T} {T_{c}} \right)^{1 / 2} , \quad ( h = 0 , T < T_{c} ) .\tag{38} $$

我们还可以从 (36) 式得到 $T_{c}$ 处磁化强度对 $h \to 0$ 外场的依赖关系。在 $T = J_{s}$ 处保持三阶展开中 $h$ 的领头阶项,我们得到

$$ m = \left( 3 \frac{h} {J_{s}} \right)^{1 / 3} , \quad ( T = T_{c} ) .\tag{39} $$

同样有必要查看 m 的完整数值解。图 10(a) 展示了不同温度下 m 随外场变化的一些示例。对于 $T > T_{c}$,其行为在 $h = 0$ 处是解析的,而随着 $T T_{c}$ 会逐渐形成一个由 (39) 式描述的奇异性。在 $T_{c}$ 以下,行为是不连续的,这对应于 $m < 0$ 和 $m > 0$ 态之间关于 h 的一级相变。这种不连续性对应于零场下的自发磁化强度,其完整数值解绘制在图 10(b) 中,且当 $T T_{c}^{-}$ 时的渐近行为由 (38) 式给出。在这里,计算中 $h \to 0^{+}$ 或 $h \to 0^{-}$ 的作用是打破 m 解之间的简并性。$h = 0$ 时的简并也对应于恰好在一级相变处两个有序态的共存,我们将在第 3.4 节更一般地讨论这一点。

将 $T_{c} = J_{s}$ 的结果(其中 $J_{s}$ 由方程 (33) 给出)具体化,对于最近邻相互作用强度为 J 的二维伊辛模型,可得 $T_{c} = 4 J$。这远高于正确值 $T_{c} / J \approx 2.269$。可以预期 $T_{c}$ 会被高估,这是由于忽略了涨落(涨落自然会导致 $T_{c}$ 降低)。

(38) 和 (39) 式中的指数 $\beta = 1 / 2$ 和 $\delta = 1 / 3$,通常是在平均场理论中这些临界指数的值。其他指数也可以计算出来 [15]。尽管幂律行为确实是连续相变的普遍特征,但平均场值通常并不正确。如前所述,对于二维伊辛模型,根据 Onsager 的精确解,$\beta = 1 / 8$。在三维中,其在伊辛普适类中的值约为 $\beta \approx 0.33$,这是通过蒙特卡洛模拟和级数展开方法(以及精度稍低的场论方法和其他解析方法)确定的。平均场临界指数在四维及更高维度中是精确的(四维是伊辛模型的上临界维度——高于该维度,平均场临界指数变得精确)。我们将在后面更详细地讨论临界行为,届时将先发展蒙特卡洛模拟的技术细节以用于实际研究临界现象。

3.2 伊辛模型的蒙特卡洛模拟#

在蒙特卡洛模拟中,目标是生成一系列自旋构型 $\sigma ( 1 ) , \sigma ( 2 ) , . . . , \sigma ( K )$,它们代表来自玻尔兹曼分布的无统计偏差的样本,即,任意构型被选入样本的概率 $P ( \sigma )$ 应与温度 T 下的玻尔兹曼权重成正比:

$$ W_{\sigma} = \mathbf{e}^{- E_{\sigma} / T} ,\tag{40} $$

这里我们采用 $k_{B} = 1$ 的单位制。实际(适当归一化的)玻尔兹曼概率为 $P_{\sigma} = W_{\sigma} / Z$,其中 Z 是配分函数:

$$ Z = \sum_{\sigma} \mathrm{e}^{- E_{\sigma} / T} .\tag{41} $$

在蒙特卡洛模拟中,我们并不需要完整的配分函数,而只需要非归一化的权重 (40)。同样重要的是要注意,序列 $\sigma ( 1 ) , \ldots , \sigma ( K )$ 可以是相关的(我们将在下文中进一步讨论),但我们只需要概率 $P ( \sigma ) \propto W_{\sigma}$,目前无需担心诸如 $P ( \sigma_{1} , \sigma_{2} )$ 这样的联合概率。

Metropolis 算法。生成有效构型序列的最简单方法是使用 Metropolis 算法 [117]。这种模拟从任意自旋构型 $\sigma ( 1 )$(例如,随机生成的)开始。之后,根据一些简单的步骤,基于随机选择的自旋并以与期望分布相关的概率进行翻转,从先前构型 $\sigma ( k )$ 随机地获得每个后续构型 $\sigma ( k + 1 )$。我们只需要存储当前构型,并且从现在起省略“时间”索引 $k$。我们用 $\sigma^{- i}$ 表示翻转了第 $i$ 个自旋后得到的构型;$\pmb{\sigma}^{- i} = ( \sigma_{1} , \ldots , - \sigma_{i} , \ldots , \sigma_{N} )$。通常,我们将一次蒙特卡洛扫描定义为 $N$ 次这样的翻转尝试,因此在这个尺寸归一化的模拟时间单位内,平均会有 $\propto N$ 个自旋被翻转。一次蒙特卡洛扫描可以通过以下简单算法执行:$^4$

$$ \begin{array} {r l r} {{\mathbf{d o j = 1 , N}}} \\ & {} & {i = \mathbf{random} [ 1 , . . . , N ]} \\ & {} & {\mathrm{i f ~} ( \mathbf{random} [ 0 - 1 ] ) < W_{\sigma^{- i}} / W_{\sigma} ) \sigma_{i} = - \sigma_{i}} \\ & {} & {\mathrm{end} \mathbf{do}} \end{array}\tag{{1}} $$

这个 Metropolis 算法基于细致平衡原理。细致平衡原理是关于随机过程(某个任意构型空间中的马尔可夫链)的一个一般性定理,该过程应当能生成一个概率分布 $W$。粗略地说,这意味着采样构型的集合应随着构型数量的增加而趋近于分布 $W$,且与初始条件无关。用 $P ( A \to B )$ 表示当前构型为 $A$ 时“转移”到构型 $B$ 的转移概率,细致平衡原理指出,如果所有转移概率满足以下条件,就能生成期望的分布:

$$ {\frac{P ( A \to B )} {P ( B \to A )}} = {\frac{W ( B )} {W ( A )}} ,\tag{42} $$

其中 $A, B$ 是所有满足 $P ( A \to B ) > 0$ 的构型对。此外,采样应该是遍历的,即原则上,任何具有非零权重 $W ( C )$ 的构型 $C$ 都必须能够通过一系列从任意构型开始的移动,以非零概率到达。

在如代码 1 所实现的 Metropolis 算法中,转移概率由两个因子组成:

$$ P ( \sigma \to \sigma^{- i} ) = P_{\mathrm{select}} ( i ) P_{\mathrm{accept}} ( \sigma^{- i} ) ,\tag{43} $$

其中 $P_{\mathrm{select}} ( i )$ 是随机选择自旋 $i$ 的概率,在此处始终等于 $1 / N$。因此,(42) 中的完整转移概率可以用实际执行(接受)自旋翻转的概率 $P_{\mathrm{accept}}$ 代替。这个概率并不是唯一的。在 Metropolis 算法中,它被取为:

$$ P_{\mathrm{accept}} ( \sigma^{- i} ) = \operatorname*{min} \left[ \frac{W_{\sigma^{- i}}} {W_{\sigma}} , 1 \right] .\tag{44} $$

可以很容易验证这满足细致平衡条件(42)。概率在形式上必须 $\leq 1$,上述通过"取最小值"函数处理了这一点。在代码 1 中,将权重比与 [0, 1) 范围内的随机数进行比较,自动实现了相同效果。注意,玻尔兹曼因子的比值仅取决于与翻转候选自旋 $\sigma_{i}$ 相互作用的那些自旋(最简单情况下仅为最近邻自旋),因此可以快速计算。更物理地说,导致能量降低的自旋翻转总是被接受,而能量升高的构型仅以概率 $P_{\mathrm{accept}} ( \sigma^{- i} ) = \exp [ ( E_{\sigma} - \bar{E}_{\sigma^{- i}} ) / T ]$ 被接受。

图11. $L \times L$ 晶格($L=8$(上)和 $L=16$(下))在温度 $T/J = 2.2 \ : ( < T_{c} \approx 2.269 )$ 下,伊辛模型蒙特卡洛模拟产生的磁化强度时间序列。两种情况下初始构型均为完全极化 $m=1$,后续各点之间间隔 $N=L^{2}$ 次 Metropolis 自旋翻转尝试(构成一次蒙特卡洛扫描)。

Metropolis 算法在经过一段瞬态时间后,会得到正确的分布,该瞬态时间取决于初始构型、温度和系统尺寸。因此,在模拟中,应在"测量"物理可观测量之前丢弃一定数量的构型。经过这种平衡后,通常会在每次或每隔几次扫描后进行测量。

在铁磁伊辛模型的情况下,最重要的计算量是磁化强度。应使用所有自旋进行计算,以利用自平均性质来改善统计量。因此,我们通常使用

$$ m = \frac{1}{N} \sum_{i=1}^{N} \sigma_{i} .\tag{45} $$

这里我们只考虑外场 $h = 0$ 下的模拟。我们已经讨论过,哈密顿量的相关对称性(这里是离散的自旋反转对称性)在有限系统的模拟中不会被破坏。因此,对于伊辛模型,应该计算自旋反转不变的期望值,如 $\langle m^{2} \rangle$ 或 $\langle | m | \rangle$,以检测进入铁磁态的相变。

对称性破缺与有限系统。为了定性理解热力学极限下的对称性破缺在大晶格中如何实际体现,首先观察一下 m 的实际蒙特卡洛时间序列是有用的。图11展示了在低于 $T_{c}$ 的温度下,两个小系统的示例。可以看出,磁化强度在正负值之间波动,并且对于更大的系统,m 符号反转的典型时间更长。在更长的时间范围内绘制这个序列会使这一点更清晰,但值得注意的是,对于 $L = 16$ 的系统,系统接近 $m = 0$ 的时间明显比 $L = 8$ 的系统短得多,这当然与通常更长的反转时间直接相关。从这些结果可以清楚地看到,当 $T < T_{c}$ 时,m 值的分布在非零的正值和负值处出现峰值。而当 $T > T_{c}$ 时,分布则是中心位于 $m = 0$ 的单峰。这类分布的示例展示在图12中。用热力学的语言来说,这两种性质不同的分布可以理解为:在低温下,自由能 $F ( m ) = E ( m ) - T S ( m )$ 由内能 E 主导(在大 m 构型中 E 较低),而在高温下则受熵 S 主导(对于小 $| m |$ 而言 S 较大)。

图 12. 在临界温度 $T_{c}/J \approx 2.27$ 以下 $(T/J = 2.2)$ 和以上 $(T/J = 2.6)$, $L \times L$ 伊辛模型的磁化强度分布,其中 $L = 16$(左)和 64(右)。

低温下的双峰磁化强度分布表明,即使一个有限系统也可以被视为有序的,尽管在模拟中自旋反转对称性并未被破坏。似乎非常合理的是,典型的反转时间会随着 L 发散(并且通过多个 L 值的模拟数据和合适的反转定义不难验证这一点),那么对于大的 L,即使在非常长的模拟中,也不会发生反转。这个发散时间尺度的原因是,为了反转磁化强度,一系列局域自旋翻转必须使系统经过许多 $m \approx 0$ 的构型,而这些构型对于越来越大系统尺寸具有越来越高的能量(因此玻尔兹曼概率更低)。因此,对于大系统来说,当 $T < T_{c}$ 时,只能在具有固定 m 符号的构型子集中进行采样——随机过程实际上变得非遍历。对称性破缺和非遍历采样严格上仅在 $N = \infty$ 时体现,但实际上,在小于典型磁化反转时间的时间尺度上,对于大但有限系统也同样适用。这个时间尺度当然取决于自旋构型如何被热采样的细节——无论是在蒙特卡洛模拟中还是在实际磁体中——但任何局域采样方案下,它都随着 $N \to \infty$ 而发散。

通常,在蒙特卡洛模拟中,我们不会研究时间序列和物理量的完整分布(尽管有时这很有用)。通过计算 $\langle m^{2} \rangle$ 或 $\langle | m | \rangle$,我们不必担心反转的时间尺度。我们通常希望将有限 N 的结果外推到热力学极限,在热力学极限下 $| \langle m \rangle | = \langle | m \dot{|} \rangle = \langle m^{2} \rangle^{1 / 2}$(但请注意,对于有限 N,因为 m 分布中峰的有限宽度,$\langle m^{2} \rangle^{1 / 2} \neq \langle | m | \rangle$)。

自相关与统计误差。在计算期望值之前,我们必须讨论如何对蒙特卡洛数据进行统计分析。使用Metropolis算法生成的连续构型在统计上并不独立——只有间隔远大于自相关时间的若干次扫过的构型才是统计独立的。量$Q$的自相关函数定义为

$$ A_{Q} ( t ) = \frac{\langle Q ( i + t ) Q ( i ) \rangle - \langle Q \rangle^{2}} {\langle Q^{2} \rangle - \langle Q \rangle^{2}} ,\tag{46} $$

其中$t$和$i$表示模拟时间,通常以先前定义的蒙特卡洛扫过为单位,平均值是对参考时间$i$进行的。归一化使得$A_{\cal Q} ( 0 ) = 1$且$A_{\cal O} ( t \infty ) = 0$。渐近衰减是指数形式的,$A_{Q} ( t ) \sim{\bf e}^{- t / \tau_{Q}}$,可用于定义自相关时间$\tau_{Q}$。通常改用积分自相关时间,它也包含了通常占主导的非渐近行为的贡献:

$$ \tau_{Q}^{\mathrm{int}} = \frac{1} {2} + \sum_{t = 1}^{\infty} A_{Q} ( t ) .\tag{47} $$

这里我们不详细讨论自相关,只总结它们在决定计算量的统计精度("误差棒")中的基本作用。

伊辛模型中$m$的自相关时间大致对应于磁化强度反转的典型时间(如图11所示)。其他量,例如$m^{2}$和$| m |$,对$m$的符号不敏感,自相关时间更短(但通常仍远长于一次蒙特卡洛扫过)。长的自相关时间不会使计算的平均值产生偏差(即,即使$\tau_{Q} \gg 1$,在每次蒙特卡洛扫过后测量$Q$也不为错),只要总模拟时间远长于$\tau_{Q}$。然而,在计算统计误差时,自相关时间确实会(显式或隐式地)发挥作用。

为了计算统计误差,可以将一次模拟划分为$B$个组块(bin),每个组块包含$M$次蒙特卡洛扫过。对于某个量$Q$,先计算每个组块的平均值$\bar{Q}_{b} , b = 1 , \ldots , B$,然后根据下式计算最终平均值$\bar{Q}$和误差棒$\sigma_{Q}$(组块平均值平均值的一个标准差):

$$ \bar{Q} = \frac{1} {B} \sum_{b = 1}^{B} \bar{Q}_{b} , ~ \sigma_{Q}^{2} = \frac{1} {B ( B - 1 )} \sum_{b = 1}^{B} ( \bar{Q}_{b} - \bar{Q} )^{2} .\tag{48} $$

那么真实期望值$Q$的最终估计应表示为$\bar{Q} \pm \sigma_{Q}$。

对数据进行组块化的原因在于,根据中心极限定理,当$M$较大时,组块平均值的分布是高斯分布(这与单次测量值的分布不同,图12中磁化强度的分布清楚地体现了这一点),因此计算出的误差线具有明确且唯一的含义(例如,我们知道$\langle Q \rangle$落在$\bar{Q} ± \sigma_{Q}$一个误差线内的概率约为$2/3$)。这仅在组块长度$M$也远大于自相关时间时才成立,此时才能认为组块平均值在统计上是独立的。如果是这种情况,误差线应仅取决于总扫面次数;${\sigma_{Q}} \propto 1 / {\sqrt{M B}}$,其中比例因子为$\propto \sqrt{\tau_{Q}}$(当然也取决于单个测量值$Q$分布的具体形式)。没有必要显式计算$\tau_{Q}$。如果对组块是否足够长有任何疑问,可以通过使用相当多的组块(例如,在100-1000范围内)并将所有组块平均值保存到磁盘上来进行检查。之后可以在模拟结束后将数据重新合并成更长的组块,并检验$\sigma_{Q}$随组块长度的收敛性。将组块平均值保存到磁盘总是明智的做法,这不仅是为了分析误差线,也是为了便于稍后添加更多数据(以便在需要时改进结果)。

需要提到的一个重要点是,任何局部更新方案(例如Metropolis算法)的自相关时间在$T < T_c$ 或 $T > T_c$ 且$N \to \infty$时都会发散。我们之前在伊辛模型的磁化反转时间中看到了这样一个例子,但发散的自相关时间也会影响任何对$m$幅值涨落敏感的量(尤其是$m^{2}$和$| m |$),而实际上这正是大多数感兴趣量的情况。因此,在接近临界点的大系统中很难获得良好的结果。在许多情况下,包括伊辛模型,这个问题可以通过使用团簇算法[118, 119]来解决(或几乎解决),在该算法中,自旋团簇(其构造满足细致平衡条件)被集体翻转(而不是逐个翻转单个自旋)。虽然我们将在第5节中考虑一种用于QMC模拟的类似环-团簇算法,但我们在此不必讨论经典团簇算法;它们在文献中已有广泛描述(例如[115]和[116])。然而,应该指出,下面将要讨论的、针对自旋数高达$1024^{2}$的系统的伊辛结果,实际上是使用团簇算法获得的。如果在临界区域中,对如此大的系统使用Metropolis方案,将很难(在合理时间内)生成同等质量的数据。

3.3 有限尺寸标度与临界指数#

图13(a)显示了二维伊辛模型在几种晶格尺寸下,平方磁化强度$\langle m^{2} \rangle$的温度依赖性。随着$L$的增大,在已知的$T_{c}$附近形成了一个越来越尖锐的特征,并且看起来非常合理的结论是:对于$T > T_{c}$,$\langle m^{2} \rangle$在$N \to \infty$时为零;对于$T < T_{c}$,$\langle m^{2} \rangle$保持非零。在此图中,误差线不可见,因为它们远小于线宽。这些数据是在非常精细的温度网格上生成的,以产生看起来连续的曲线。

$\langle m^{2} \rangle$ 的尺寸依赖性如图13(b)所示,图中给出了三个不同温度下的数据。对于 $T < T_{c}$,它随着 $L$ 收敛到非零值,而对于 $T > T_{c}$,它则以 $L^{- 2}$ 的形式衰减到0。恰好在 $T_{c}$ 处,衰减遵循一个非平凡幂律:$\langle m^{2} \rangle \propto L^{- 2 \beta / \nu}$,其中 $\beta = 1 / 8$ 与 $T < T_{c}$ 时 $N = \infty$ 磁化强度的指数相同,而 $\nu$(关联长度指数,将在下文中讨论)对此模型等于1。这种行为是有限尺寸标度的一个例子。在临界点处的完美幂律行为仅在 $L \to \infty$ 时严格成立,而对于小系统,可能存在显著的标度修正(在此讨论的伊辛模型中,这些修正异常小)。通常我们不知道 $T_{c}$ 的精确值(这当然是进行模拟的原因之一),因此必须设计确定临界点的方法。有很多方法可以实现这一点,所有这些方法都以某种方式基于这样一个事实:在临界点处的大晶格上,序参量及相关量应表现为非平凡幂律(指数已知或未知)。

请注意,磁化强度的平方与自旋关联函数相关,

$$ C ( \mathbf{r}_{i j} ) = \langle \sigma_{i} \sigma_{j} \rangle ,\tag{49} $$

图13. (a) 多个 $L \times L$ 伊辛系统的磁化强度平方随温度的变化。(b) 在 $T_{c}$ 处、低于和高于 $T_{c}$ 时的有限尺寸依赖性。在 $T_{c}$ 处,行为是幂律 $m^{2} \sim L^{- 1 / 4}$,在此对数-对数坐标下对应斜率为 $- 1 / 4$ 的直线。对于 $T > T_{c}$,衰减到零的形式为 $L^{- 2}$,这也是 $T < T_{c}$ 时的收敛速度

其中$\mathbf{r}_{i j}$是自旋之间的距离。如果存在长程序,那么在无限晶格中当$r \to \infty$时,$C ( r ) \to \langle m^{2} \rangle$。对于有限系统,在$L \to \infty$的极限下,对于$C ( r_{\operatorname*{max}} )$同样成立,其中$r_{\mathrm{max}}$是周期性晶格上的最远距离,例如对于$L \times L$晶格,$r_{\mathrm{max}} = \sqrt{2} L$。如果不存在序,那么$C ( r ) \to 0$,在$T_{c}$处按照幂律衰减,在$T > T_{c}$时按照指数衰减。在$T_{c}$以下,无限系统中"连接"的关联函数,

$$ C^{*} ( r ) = C ( r ) - \langle m \rangle^{2} ,\tag{50} $$

当$r \to \infty$时指数衰减到零。$C ( \boldsymbol r )$和$C^{*} ( r )$的指数形式都由一个关联长度$\xi$表征,该关联长度在趋近$T_{c}$时从两侧发散。另外注意,$\langle m^{2} \rangle$可以精确地写成自旋关联的和;

$$ \langle m^{2} \rangle = \frac{1} {N} \sum_{\bf r} C ( {\bf r} ) .\tag{51} $$

$\langle m^{2} \rangle$ 的 $L^{-2}$ 收敛性(在 $T_{c}$ 上下均成立)与有限关联长度有简单的关系。

这里需要指出,$T_{c}$ 以下的自旋关联性质取决于我们所研究的序参量的对称性。对于矢量序参量,例如在海森堡模型(具有经典或量子自旋)中,对称性破缺系统的自旋分量平行(纵向)于序参量的关联呈指数衰减,与伊辛模型类似。然而,横向关联呈幂律衰减。这与序参量的连续对称性有关。在伊辛模型中,对称性是离散的——有序态破坏了自旋反转对称性,局域磁化涨落的自由能代价很大(与翻转畴的边界成正比)。这导致了有序相中关联函数(连接关联函数)呈指数衰减。另一方面,在海森堡模型中,存在无能隙的自旋波激发(这些激发是自旋在垂直于有序化方向上的激发,更一般地称为戈德斯通模)。由这些导致的涨落使得横向自旋关联呈幂律形式。因此横向关联长度形式上是无穷大的,并且在未明确破坏旋转对称性的有限格子计算中(因此计算出的 $\langle m_{s}^{2} \rangle$ 包含纵向和横向关联的贡献),$\langle m_{s}^{2} \rangle$ 的尺寸修正正比于 $1 / L$。我们在 2.4 节已经讨论过二维 $S = 1 / 2$ 海森堡反铁磁体基态中的这种行为,图 5(b) 非常清楚地显示了 $1 / L$ 修正。

3.3.1 标度性与临界指数#

为了更详细地讨论临界点附近的有限尺寸标度,我们首先需要回顾一些关于热力学极限下临界现象的基本知识。这里我们仅列出一些关键结果和定义;更多细节可参考任何关于临界现象的标准教材(例如,Cardy [15] 的著作)。

关联长度是相变和临界现象理论中最重要的概念之一。在无穷大系统中,当接近临界温度时,关联长度根据幂律发散;

$$ \xi \sim | t | ^{- \nu} .\tag{52} $$

从高于或低于 $T_{c}$ 靠近时,指数相同,但式 (52) 中的前因子在 $t \to 0^{+}$ 和 $t \to 0^{-}$ 时通常不同(然而,这些前因子之比——称为振幅比——通常是普适的)。在平均场理论中,$\nu = 1 / 2$。在 $T_{c}$ 处,尽管关联长度形式上是无穷大的,但系统尚未有序。相反,对于无穷大系统,恰在 $T_{c}$ 处的关联函数具有幂律形式,

$$ C ( {\bf r} ) \sim r^{- ( d - 2 + \eta )} ,\tag{53} $$

其中 $d$ 是系统的维度,$\eta$ 是另一个临界指数(也称为反常维度,因为它可以与临界点有序畴的分形维度相关联),其平均场值为 $\eta = 0$。长程序仅在略低于 $T_{c}$ 时出现,此时渐近长程关联趋于常数;$C ( r \to \infty ) = \langle m^{2} \rangle$。

我们已经讨论了对称性破缺态中当 $t \to 0^{-}$ 时序参量的行为;

$$ \langle{m} \rangle \sim | t | ^{\beta} .\tag{54} $$

我们还会对相应的磁化率感兴趣,其定义为

$$ \chi = \left. \frac{d \langle m \rangle} {d h} \right | _{h = 0} = \frac{N} {T} \big ( \langle m^{2} \rangle - \langle | m | \rangle^{2} \big ) ,\tag{55} $$

其中 $h$ 是与序参量耦合的场的强度,例如,在伊辛铁磁体的情况下,哈密顿量中有一项 $- h \textstyle \sum \sigma_{i}$。式 (55) 中的最后一个表达式明确显示了 $\chi$ 如何直接与序参量的涨落相关。磁化率在 $T_{c}$ 处发散;

$$ \chi \sim | t | ^{- \gamma} .\tag{56} $$

因此,当趋近临界点时,系统对耦合到序参量的场 $h$ 变得无限敏感,并且在 $T_{c}$ 处,线性响应形式 $\langle m \rangle = \chi h$ 不再成立。相反,在 $T_{c}$ 处,序参量依赖于(弱)场,其关系为 $h^{1 / \delta}$。根据方程 (39) 的结果,该指数的平均场值为 $\delta = 3$。

比热在 $T_{c}$ 处也是奇异的,

$$ C \sim | t | ^{- \alpha} ,\tag{57} $$

其中 $\alpha$ 可以为正或负。在平均场理论中,$\alpha = 0$。当 $\alpha < 0$ 时,不存在发散,只在 $T_{c}$ 处有一个尖点奇异性。在某些情况下,例如对于二维伊辛模型,$\alpha = 0$,但比热仍然存在弱对数发散。

我们以上遇到的临界指数 $\nu , \beta , \eta$ 等并非彼此完全独立。指数之间的关系由重整化群理论解释,在我们这里讨论的有序-无序相变类型中,该理论表明存在两个更基本的指数,可观测的物理指数可以用它们来表达 [15]。即使在这套终极相变理论出现之前(最完整的形式由 Widom 在 1960 年代中期给出),人们也已利用其他论证找到了指数关系,例如:

$$ \begin{array} {l} {\gamma = \nu ( 2 - \eta ) ,} \\ {\nu d = 2 - \alpha ,} \\ {\alpha + 2 \beta + \gamma = 2 .} \end{array}\tag{58} $$

这些关系对于检查指数数值计算的自洽性非常有用。涉及维度 $d$ 的指数关系称为超标度假定,其普适性不如其他关系。例如,它们在平均场理论中不成立,因此对于处于或高于上临界维度的任何系统也不成立。

3.3.2 有限尺寸标度假定#

有限尺寸标度理论 [120] 的基本假设是,当(无限系统的)关联长度 $\xi$ 与有限系统长度 $L$ 相当时,偏离无限大尺寸临界行为的情况就会发生。如果 $L \gg \xi$,则系统是有限的事实应该无关紧要,无限大尺寸的行为适用。另一方面,如果 $L \ll \xi$,那么最重要的长度尺度应该是 $L$,而非 $\xi$。为了理解这两个长度尺度是如何起作用的,用关联长度来表示感兴趣的物理量是有帮助的。考虑一个量 $Q$,它在 $T_{c}$(约化温度 $t 0$)处表现出幂律发散行为:

$$ Q \sim | t | ^{- \kappa} ,\tag{59} $$

例如磁化率(此时 $\kappa = \gamma$)。我们可以利用方程 (52) 将 $| t |$ 表示为关联长度的函数:

$$ | t | \sim \xi^{- 1 / \nu} ,\tag{60} $$

并由此将 $Q$ 写为:

$$ Q \sim \xi^{\kappa / \nu} .\tag{61} $$

这种形式应在 $\xi \ll L$ 时适用,但当 $\xi \sim L$ 时,有限晶格上的发散无法继续。因此,$Q$ 在有限晶格上能达到的最大值 $Q_{\mathrm{max}}$ 可以通过将方程 (61) 中的 $\xi$ 替换为 $L$ 得到,即:

$$ Q_{\operatorname*{max}} ( L ) \sim L^{\kappa / \nu} .\tag{62} $$

同样地,从方程 (60) 我们还可以推导出 $\xi$ 达到 $L$ 时约化温度的标度行为,这也就是量 $Q$ 达到最大值时的温度:

$$ \left. t_{\mathrm{max}} ( L ) \right. \sim L^{- 1 / \nu} .\tag{63} $$

需要注意的是,这只是一个比例关系,$L^{- 1 / \nu}$ 前的系数取决于所考虑的量。这种偏移也适用于非发散的量——任何在 $T T_{c}$ 处表现出奇异行为的特征都应该以 $L^{- 1 / \nu}$ 的速率发生偏移。如果某个量在 $T = T^{*} ( L )$ 处有一个明确的极大值或其他可区分的特征,那么这个温度就可以用作依赖于尺寸的临界温度(同样,这个温度并不是唯一的,而是取决于所考虑的量)。

有限尺寸标度律 (62) 和 (63) 源自一个更一般的有限尺寸标度假说 [120],该假说与无穷大系统的标度理论类似,最初是基于唯象考虑提出的,但后来利用重正化群理论推导出来。该假说认为,在热力学极限下,一个在 $T_{c}$ 处奇异可观测量,在接近 $T_{c}$ 时,其标度行为为 $L$ 的幂次乘以比值 $\xi / L$ 的一个非奇异函数。因此,任何奇异量(不一定发散)都应具有如下形式

$$ Q ( t , L ) = L^{\sigma} f ( \xi / L ) ,\tag{64} $$

利用 $\xi \sim | t | ^{- 1 / \nu}$,我们也可以将其写为

$$ Q ( t , L ) = L^{\sigma} g ( t L^{1 / \nu} ) .\tag{65} $$

这个标度律在临界点以上 $( t > 0 )$ 和以下 $( t < 0 )$ 都应成立。恰好在 $T_{c}$ 处,我们恢复尺寸标度 $Q ( 0 , L ) \sim L^{\sigma}$ 。为了与标准临界指数联系起来,我们可以利用这样一个事实:对于接近0的固定 $t$,随着系统增大,任何 $t \neq 0$ 的行为最终都必须由方程 (59) 给出;$\mathcal{Q} ( t , L \infty ) \sim | t | ^{- \kappa}$(对于奇异非发散量,$\kappa$ 为负值,例如,对于序参量,我们有 $\kappa = - \beta )$ 。为了得到这种形式,(65) 中的标度函数 $g ( x )$ 在 $x \to \infty$ 时必须渐进行为 $g ( x ) \sim x^{- \kappa}$ 。为了使 (65) 中的尺寸依赖性相互抵消,我们因此得出结论 $\sigma = \kappa / \nu$ ,即

$$ Q ( t , L ) = L^{\kappa / \nu} g ( t L^{1 / \nu} ) .\tag{66} $$

为了利用数值数据提取标度函数 $g ( x )$,我们可以定义

$$ y_{L} = Q ( t , L ) L^{- \kappa / \nu} , x_{L} = t L^{1 / \nu} ,\tag{67} $$

图14. 不同 $L \times L$ 晶格上伊辛模型磁化率 (55) 的蒙特卡洛结果。(a) 展示了温度依赖性,其中竖线表示 $T_{c}$ 。注意纵轴为对数坐标。(b) 中数据已使用伊辛指数的精确值 $\gamma = 7 / 4$ 和 $\nu = 1$ ,以及 $t = ( T - T_{c} ) / T_{c}$ 中 $T_{c}$ 的精确值进行了标度。

并针对不同系统尺寸绘制 $y_{L}$ 关于 $x_{L}$ 的图。如果标度假说成立,不同(大)系统尺寸的数据应落在同一条曲线上,该曲线即为标度函数(这称为曲线相互叠合);$g ( x ) = y_{L \to \infty} ( x )$ 。图14利用二维伊辛模型磁化率的蒙特卡洛数据说明了这一点。(a) 面板中峰的位置明显随着 $L$ 增大向已知的 $T_{c}$ 移动。根据上述程序对数据进行标度后,如 (b) 面板所示,在 $t = 0$ 附近曲线确实几乎相互叠合,但在远离临界点时,较小的系统会出现偏差。这些偏差源于对标度的修正,原则上可以用次主导指数来描述。

我们可以将标度形式 (66) 应用于关联长度本身,此时 $\kappa = \nu$ ,且 $L$ 标度与模型特定的指数无关。在普适类预先未知的情况下,这有助于通过曲线叠合 $\xi / L$ 数据来提取指数 $\nu$ ,而无需同时调整 (66) 中的另一个指数。

关联长度的实用定义。关联长度可以通过多种方式定义,不一定仅仅基于关联函数的渐近衰减(这通常难以可靠地提取)。一种常用且实用的关联长度定义基于关联函数的傅里叶变换,通常称为(静态)结构因子,

$$ S ( \mathbf{q} ) = \langle \sigma_{- \mathbf{q}} \sigma_{\mathbf{q}} \rangle = \sum_{\mathbf{r}} \mathrm{e}^{- i \mathbf{q} \cdot \mathbf{r}} C ( \mathbf{r} ) = \sum_{\mathbf{r}} \cos ( \mathbf{q} \cdot \mathbf{r} ) C ( \mathbf{r} ) ,\tag{68} $$

其中 $\sigma_{\mathbf{q}}$ 是单个自旋构型的傅里叶变换,

$$ {\sigma}_{\mathbf{q}} = \frac{1} {\sqrt{N}} \sum_{j} {\sigma}_{j} \mathrm{e}^{- i \mathbf{q} \cdot \mathbf{r}_{j}} .\tag{69} $$

我们用 $\mathbf{Q}$ 表示主导关联的波矢——对于铁磁体 $Q = 0$,对于二维反铁磁体 $\mathbf{Q} = ( \pi , \pi )$,等等。为简化符号,q 将用于表示偏离 Q 的波矢。那么 $q_{1} = 2 \pi / L$ 对应于最接近 $\mathbf{Q}$ 的波矢之一,例如 $\mathbf{Q} + ( 2 \pi / L ) \hat{\mathbf{x}}$ ,其中 $\hat{\mathbf{x}}$ 是倒易空间中 x 方向的单位矢量。

可以使用 $q = 0$ 和 $q_{1}$ 处的结构因子定义关联长度 $\xi_{a}$;

$$ \xi_{a} = \frac{1} {q_{1}} \sqrt{\frac{S ( 0 )} {S ( q_{1} )} - 1} .\tag{70} $$

可以证明,对于 d 维晶格,如果关联函数由奥恩斯坦-泽尼克形式(在平均场处理中得到 [15])给出,

$$ C_{O Z} ( \mathbf{r} ) \sim r^{- \frac{1} {2} ( d - 2 )} \mathrm{exp} ( - r / \xi ) ,\tag{71} $$

那么 $\xi_{a}$ 与该关联函数指数衰减中出现的原始关联长度 $\xi$ 的关系为

$$ \xi_{\mathrm{a}} = \xi \sqrt{\frac{( 1 + d ) ( 3 + d )} {8 d}} .\tag{72} $$

因此,对于$d = 1$和$3$,$\xi_{a} = \xi$,而二维情况是特殊的,有$\xi_{\mathrm{a}} = \xi (15 / 16)^{1 / 2}$(或者,可以说$d = 1, 3$是特殊情况,因为对于所有$d > 3$,该因子也都不等于1)。奥恩斯坦-泽尼克形式通常对无序相在$r \gg \xi$时成立[15]。短距离下与该形式的偏离意味着(72)并不严格成立,但无论短距离行为如何,当$\xi \infty$时该关系式严格成立。

在长程有序经典系统中,对于任意$L$,$\xi_{a}$通常在$T \to 0$时发散,因为基态中没有涨落(因此对于$q \neq 0$,结构因子为零)。为了从有序系统中移除关联函数非衰减部分的贡献,我们可以使用连通关联函数(50)。尽管对于有限系统,$\langle m \rangle^{2}$的定义并不唯一,但我们可以,例如,减去$C(\mathbf{r}_{\operatorname*{max}})$。另一种替代方法是,基于$q_{1}$和$q_{2} = 2 q_{1} = 4 \pi / L$处的结构因子,使用另一种关联长度定义

$ \xi_{\mathsf{b}} = \frac{1} {q_{1}} \sqrt{\frac{S ( q_{1} ) / S ( q_{2} ) - 1} {4 - S ( q_{1} ) / S ( q_{2} )}} 。\tag{73} $

可以证明,如果$C(\mathbf{r})$由奥恩斯坦-泽尼克关联函数给出,并且当$T < T_{c}$时还加上一个对应于长程有序的常数[因为这仅影响$S(0)]$,那么(72)对于该定义也成立。

图15展示了二维Ising模型的$\xi_{a} / L$和$\xi_{b} / L$的结果。注意,在$T \to 0$发散的$\xi_{a} / L$的情况下使用了对数坐标,而$\xi_{b} / L$是收敛的,因此在线性坐标上绘制。两个量在$T_{c}$处都表现出尺寸无关性(曲线相互交叉),但那里的值明显不同。这是因为关联函数的奥恩斯坦-泽尼克形式仅渐近地适用于$T > T_{c}$,并且没有理由认为$\xi_{a}$和$\xi_{b}$这两个定义在$T_{c}$处应该完全相同(尽管它们的值应该相关)。对于无序相中接近$T_{c}$的大系统,它们的值非常相似。交叉点(或$\xi_{b}$的峰值位置)可用于提取未知$T_{c}$的系统中的$T_{c}$。温度轴也可以像图14中那样进行缩放,以提取关联长度指数。

图15. 两个关联长度定义(式(70)和(73))的温度依赖性,对于$L \times L$ Ising模型通过尺寸L归一化。竖线表示$T_{c}$。

注意,对于所需的少量q点,可以使用蒙特卡洛模拟中生成的自旋构型的傅里叶变换(69)有效计算$S(\mathbf{q})$。由于结构因子是实值,我们有

$ S(\mathbf{q}) = \langle \mathrm{Re} \{\sigma_{\mathbf{q}} \}^{2} \rangle + \langle \mathrm{Im} \{\sigma_{\mathbf{q}} \}^{2} \rangle 。\tag{74} $

在模拟后计算完整的关联函数C(r)并对其进行傅里叶变换则要耗费更多时间。

Binder比和累积量。除了$\xi / L$之外,还有其他无量纲量在临界点处与尺寸无关,对于独立于临界指数值提取$T_{c}$很有用。最常用的可能是Binder比[121, 122]:

$ R_{2} = \frac{\langle m^{4} \rangle} {\langle m^{2} \rangle^{2}} 。\tag{75} $

在$T_{c}$处,幂律相互抵消,该比值与$L$无关(并且也是普适的),仅存在次主导的有限尺寸修正。通常,绘制不同系统尺寸的$R$对$T$的曲线会产生在$T_{c}$附近相交的曲线。通过定位成对系统尺寸(例如$L$和$2L$)的$R_{2}$交叉点,可以得到一个依赖于尺寸的临界点,该点通常比(63)中$L^{-1/\nu}$的偏移收敛得更快。可以将其理解为:在涉及两个系统尺寸的量中,主导修正相互抵消,从而只剩下一个根据更快的、高阶标度修正趋近于$T_{c}$的量。还可以基于$m$的其他幂次定义类似于(75)的比值,例如$R_{1} = \langle m^{2} \rangle / \langle | m | \rangle$。通过曲线交叉法定位$T_{c}$的方法同样可以应用于$\xi / L$。

Binder比还具有其他有趣的性质。在标量序参量的情况下(例如对于伊辛模型),Binder累积量定义为

$ U_{2} = \frac{3} {2} \left( 1 - \frac{1} {3} R_{2} \right) 。\tag{76} $

图16 $L \times L$伊辛模型的Binder累积量(76)。在(a)中可以看到,随着$L$增大,极限值$U_{2} \to 0$(对于$T > T_{c}$)和$U_{2} \to 1$(对于$T < T_{c}$)的趋近过程。在(b)中,靠近$T_{c}$(垂直线)的数据以更详细的尺度和更大的$L$绘制,以显示曲线的交叉。

在有序态中,当$N \to \infty$时,$U_{2} \to 1$,因为磁化强度分布$P ( m )$趋近于在$\pm \langle | m | \rangle$处的两个delta函数,因此$R_{2} \to 1$。相反,在无序相中,$m$的涨落围绕$m = 0$呈高斯分布(由中心极限定理得出,因为在大系统中,距离$\gg \xi$的区域内涨落是不相关的)。基于高斯积分,$R_{2} \to 3$且$U_{2} \to 0$。

将Binder累积量推广到n分量序参量(其中伊辛磁化强度对应$n = 1$,XY模型对应$n = 2$,等等)时,应记住(75)中的$m^{2} = \mathbf{m} \cdot \mathbf{m}$和$\overline{{m^{4}}}$对序参量的角度涨落不敏感。然而,对$| m |$的高斯分布在n维空间上进行积分以计算(75)中的平均值,会引入依赖于$n$的因子。因此,为了再现$U_{2}$的上述性质,必须将其定义一个一般序参量形式:

$ U_{2} = {\frac{n + 2} {2}} \left( 1 - {\frac{n} {n + 2}} R_{2} \right) 。\tag{77} $

图16展示了二维几种格点上$U_{2}$作为温度T函数的蒙特卡洛结果。随着系统尺寸的增大,它逐渐演化为一个在$T_{c}$处的阶梯函数,这一过程可以清晰看到。在这种情况下,所有曲线都在非常接近已知$T_{c}$的位置相交,反映出该模型中非常小的次主导修正。在其他系统中,交点往往会出现一些漂移,因此需要对交点进行仔细的外推(例如,基于尺寸为L和$2L$的数据集,或基于其他长宽比)。

有限尺寸标度在实践中的应用。我们简要讨论如何在实际中进行有限尺寸数据塌缩。涉及的参数数量(即$T_{c}$以及一个或两个临界指数,还有可能包括下面将讨论的次主导修正的指数)相当少,通常我们可以直接从原始数据观察和初步实验中获得这些参数的一些粗略知识,例如,通过定位非平凡的幂律行为(如图13所示),这通常就足以确定方程(62)和(66)中的指数比$\kappa / \nu$。对Binder累积量或$\xi / L$的分析可能已经能以足够的精度给出$T_{c}$,但检验其他拟合对其值的敏感性可能仍然有用。得益于现代计算机的强大能力,作为使用复杂的多维优化程序的替代方案,我们可以编写一个简单的暴力搜索计算机程序,在合适的有限网格上寻找最佳参数集。由一组参数产生的数据塌缩的质量,可以通过用一个单一的高阶多项式同时拟合所有可用的L的缩放数据点(方程(67)中定义的$\left( x_{L} , y_{L} \right)$)而得到的$\chi^{2}$值来量化。通常,对于不同的耦合和系统尺寸会有大量的数据点,并且由于缩放函数应该具有良好的行为,一个合理阶数(粗略指引为3到8阶)的多项式在数据可塌缩的窗口内应该能很好地工作。这个窗口的大小也必须调整,直到数据塌缩良好。例如,对于图14中的数据,x在$\left( -0.5, 0.5 \right)$范围内应该是合适的,尽管窗口也取决于分析中包含的系统尺寸以及误差棒(它决定了被忽略的次主导标度修正的敏感度)。

图17. (a) 二维伊辛模型磁化率的最佳拟合标度塌缩;与图14数据相同,但包含了更大的晶格,并调整了$t = (T - T_{c}) / T_{c}$中的$T_{c}$以及指数$\nu, \gamma$,以最小化相对于一个多项式形式(这里是四阶,如实线所示)的标度函数$g(x)$的$\chi^{2}$。该数据集的最优值为:$T_{c} / J = 2.26921 \pm 0.00002$,$\nu = 0.9985 {\scriptstyle \pm 0.0011}$,$\gamma = 1.750 {\scriptstyle \pm 0.002}$,这里误差棒(一个标准差)是通过多次向数据添加高斯噪声(幅度等于蒙特卡洛误差棒)重复拟合计算得出的。(b) 展示了减去拟合的标度函数$g(x)$后的数据,使得误差棒可见。该拟合在统计上是合理的,$\chi^{2} \approx 0.9$(每自由度)。

在这种拟合中,确定参数可靠误差棒并不总是容易的。除了纯粹的统计误差(可以通过例如多次向数据添加大小等于误差棒的高斯噪声并重复$\chi^{2}$最小化来确定)之外,还存在由于标度修正导致的系统误差,这部分难以估计。如果$\chi^{2}$在统计上合理(即每自由度接近1),通常可以假设被忽略的修正对参数的影响并未超出统计不确定性范围。

图17展示了一个数据塌缩的例子。利用 $L \in \{64 - 512 \}$ 在窗口 $| t | L^{1 / \nu} < 0.5$ (其中 $t$ 包含程序中调整的变量 $T_{c}$)内的数据,确定了二维伊辛模型的临界指数和 $T_{c}$。图注中列出的结果与已知的 $T_{c}$ 和二维伊辛指数吻合极好,且拟合在统计上也非常好。使用相同的拟合窗口,如果包含更小的格子,则无法获得统计上可接受的拟合,例如,包含 $L = 32$ 会得到 $\chi^{2} \approx 2$(每自由度),对于该拟合的自由度数而言略微偏大。在这种情况下,$T_{c} = 2.26924 \pm 0.00002$,比真实值高出约3个误差条,而指数则在两个误差条内正确。

对标度性的修正。作为检查,通常一个好主意是包含对主导标度形式的一些修正,有时为了获得良好拟合甚至必须这样做。最常用的方法是通过下式进行数据塌缩

$$ y_{L} = Q ( t , L ) L^{- \kappa / \nu} ( 1 + a L^{- \omega} )^{- 1} , ~ x_{L} = t L^{1 / \nu} ,\tag{78} $$

其中 $\omega$ 是次主导指数,a 是常数。这里 $x_{L}$ 标度(临界点的有限尺寸偏移)没有对标度性的修正,但原则上也可以包含这种修正 [85]。通常,当包含不使用修正时必须排除的晶格尺寸时,也能通过修正获得良好拟合。更大的系统尺寸范围可以部分补偿因包含更多参数而导致所有参数统计不确定性增加的事实。如果在包含和不包含修正的拟合中均获得一致的领先临界指数,那么就有理由确信结果是正确的。

3.4 一级相变#

我们在上一节讨论的标度性质适用于关联长度发散的连续相变。在一级(不连续)相变中,关联长度在相变点处保持有限,且序参量以及其他量呈现不连续跳跃。这些不连续性在无限系统极限下发展,通常遵循幂律,这些幂律也可以使用有限尺寸标度技术进行研究。与这些幂律相关的指数通常与系统维度有平凡的关系 [123, 124, 125]。例如,在一级相变中,比热随系统尺寸按 $L^{d}$ 发散,而在连续相变中,则按 $L^{\alpha / \nu}$ 发散(其中 $\alpha$ 通常非常小甚至为负)。临界点随系统尺寸的偏移按 $L^{- d}$ 标度,而连续临界点的偏移则为 $L^{- 1 / \nu}$。

尽管指数等于维度 $d$ 的有限尺寸标度原则上使得识别一级相变变得容易,但对弱一级相变的研究却很困难,因为其主导标度形式存在大的修正。此时,使用实际可得的系统尺寸,可能难以清楚地区分缓慢发展的不连续性与连续相变处的较弱奇异行为。

强一级相变的研究也很困难,原因则完全不同。蒙特卡罗模拟可能陷入亚稳态,此时计算出的量并不对应正确的热平均(另一方面,这与真实系统完全类似,对真实系统而言,亚稳性和滞后效应是一级相变的标志)。这甚至可能导致难以准确定位相变点。为了缓解这些问题,人们开发了多种多正则蒙特卡罗方法,这些方法在扩展系综中对构型进行采样,其中系统温度也在波动(变温或并行回火方法 [126, 127]),或者采用不同于玻尔兹曼概率的分布 [128, 129, 130](测量结果据此重新加权),这使得系统更容易探索构型空间。

3.4.1 相共存#

一级相变最重要的特征之一就是相共存。在相变点处(实际上在一个随系统尺寸增大而收缩至单一温度的窄窗口内)进行的蒙特卡罗模拟中,这一点表现为产生两类构型,分别对应于相变温度上下存在的两个不同相。前提是能够遍历地采样整个构型空间,如上述讨论,这在实践中并非总能实现。在小系统中,也会存在无法明确归因于其中一相的构型——它们对应于系统在两相之间转变所需的涨落(畴壁)。图18示意性地展示了一个伊辛序参量的情况。能够观察到三峰分布的转变窗口随系统尺寸增大而迅速变窄,并且这些峰演变为δ函数。当峰之间的权重变得非常小时,实际上可能无法再遍历地采样构型,因为系统会陷于仅对应其中一个峰的子空间内(这方式正与第3.2节讨论的对称性破缺完全类似)。

注意,图18中的序参量分布类型并不适用于图10中磁化曲线所示的伊辛模型场致一级相变。在该情况下,仅存在两个峰,当磁场调谐经过 $h = 0$ 时,权重在这两个峰之间转移。在 $h = 0$(且 $T < T_{c}$ )处的相共存对应于图12分布中的两个对称峰。通过考虑 $| m |$ 的一维分布,可以更清晰地看出这种情况与顺磁-铁磁相变(具有一般n分量序参量m)之间的相似性。当有序相和无序相共存时,该分布有两个峰,一个位于 $| m | = 0$,另一个位于 $| m | > 0$。再次注意,在连续相变中(图12),中央峰在 $T_{c}$ 以下连续分裂成两个峰,这与一级相变中在 $| m |$ 非零值处出现两个有序峰的情况形成对比。

对于一般的n维序参量,Binder 比率(式(77))仅取决于 $| m |$ 的分布。在一级相变处,它具有一个非常有趣的性质。可以利用例如理想化的共存分布 $P ( | m | ) = ( 1 - p ) G ( | m | ) + p \delta ( | m | - m^{*} )$ 来轻松检验,其中 $G$ 是归一化的“半高斯分布”(定义在 $| m | \geq 0)$ 上且 $0 < m^{*} \leq 1$ ,当有序相权重 $p$ 从0调谐到1时,并且如果高斯分布的宽度趋于零(对应于无穷大系统尺寸),$U_{2}$ 会表现出负发散。虽然完整的m分布当然包含更多信息,但定位负Binder累积量的窗口并检查其峰值是否发散,可能是分析一级相变的一种实用方法。

图18. 有限伊辛铁磁体在一级相变附近磁化强度分布的(示意性)演化。当温度调谐经过相变区域时,权重在中心峰(对应无序相)和两个非零磁化强度峰(两个铁磁态)之间快速转移。在无穷大系统尺寸极限下,这些峰变成δ函数。相变温度 $T_{c}$ 可以定义为这三个峰具有某些特定(但本质上任意)特征(如等权重或等峰高)的点。

3.4.2 阻挫伊辛模型(Frustrated Ising model)#

这里我们研究一个一级相变的具体例子,即哈密顿量如下的阻挫二维正方晶格伊辛模型

$$ E_{\sigma} = - J_{1} \sum_{\langle i j \rangle_{1}} \sigma_{i} \sigma_{j} + J_{2} \sum_{\langle i j \rangle_{2}} \sigma_{i} \sigma_{j} ,\tag{79} $$

其中两个耦合常数 $J_{1} , J_{2} > 0$(需注意系数前的不同符号)。第一项($\langle i j \rangle_{1}$ 表示最近邻自旋)是标准的伊辛铁磁体,而第二项($\langle i j \rangle_{2}$ 表示次近邻自旋,即 $2 \times 2$ 方格中的对角线,如图7所示)是反铁磁相互作用,并导致阻挫。对于耦合比 $g = J_{2} / J_{1} < 1 / 2$,系统的基态是铁磁态(完全极化)。在极限 $g \to \infty$下,系统简化为两个解耦的反铁磁体,具有四种条纹(或共线)基态,图7(c)展示其中一种。这些基态对所有 $g > 1 / 2$ 都成立。在 $g = 1 / 2$ 点处,基态是高度简并的。注意,由于 $E_{\sigma}$ 在棋盘子晶格之一的 $\sigma_{i} - \sigma_{i}$ 变换下具有不变性,$J_{1} < 0$ 的情形等价于 $J_{1} > 0$;因此在此情形下,铁磁相被反铁磁相取代,而条纹相不变。

阻挫伊辛模型已被研究很长时间,但其相图的某些特征仍存争议,或直到最近才得以解决[95, 96, 97]。当 $g < 1 / 2$ 时,相变属于标准伊辛类型,但由于较大的标度修正和较长的蒙特卡洛自相关时间,在接近 $g = 1 / 2$ 处难以研究。这里我们只考虑 $g > 1 / 2$,在耦合强度达到 $g^{*} \approx 0.8$ [97] 之前,高温顺磁相和低温条纹相之间存在一级相变,之后相变变为连续相变。我们分析使用标准单自旋Metropolis算法得到的结果(因为团簇蒙特卡洛方法在存在阻挫时无效)。温度将以 $J_{1}$ 为单位给出。

图19. 在耦合比 $g = 0.55$、尺寸 $L = 128$ 的二维阻挫伊辛模型的 $( m_{x} , m_{y} )$ 平面(全空间 $| m_{x} | \leq 1 , | m_{y} | \leq 1$)中的序参量分布。较亮特征对应更高的概率密度。各温度值(图中以 $J_{1}$ 为单位标注)处于一级相变区域,相共存(一个中心峰以及对应于x方向和y方向条纹的四个峰)在 $g = 0.772$ 时最为清晰

序参量分布。有序相的条纹可以沿x轴或y轴取向,对应的序参量为

$$ m_{x} = \frac{1} {N} \sum_{i = 1}^{N} \sigma_{i} ( - 1 )^{x_{i}} , m_{y} = \frac{1} {N} \sum_{i = 1}^{N} \sigma_{i} ( - 1 )^{y_{i}} ,\tag{80} $$

其中 $x_{i}$ 和 $y_{i}$ 是格点 i 的(整数)晶格坐标。我们首先考察序参量分布 $P ( m_{x} , m_{y} )$。如果在 $T_{c}$ 以下四个可能的有序态被等概率取样,则分布应呈现四重对称性,峰值位于 x 轴和 y 轴的负半轴和正半轴上。在顺磁相中应只有一个中心峰值。因此,一级相变处的共存应体现在该情形下存在五个峰值(而非图18中标量序参量的三峰分布)。图19展示了 $g = 0.55$ 时一个 $L = 128$ 的系统在转变窗口内三个温度下的结果。在所有情形下都能观察到非常明显的四重对称性。在 $T = 0.771$ 时,四个对称峰清晰地表明了有序相的存在,但存在强烈的涨落,这反映在分布向中心延伸的权重上。在 $T = 0.773$ 时,分布峰值位于中心,但权重也远延伸至有序区域。在这两个温度之间,即 $T = 0.772$ 时,可以同时辨别出四个有序峰和一个中心峰。所有这些图都显示出共存的特征;即使没有出现五个峰,有序构型和顺磁构型仍然具有显著的概率。远离狭窄的转变窗口,分布会迅速转变为具有明显有序或无序特征的分布。

Binder累积量。我们使用式(77)中的 $m^{2} = m_{x}^{2} + m_{y}^{2}$ 来分析 Binder 累积量。现在的问题是在该定义中应使用何种序参量维度 n。在高温下,大系统不同部分中 $m_{x}$ 和 $m_{y}$ 的涨落是不相关的。这意味着一个旋转对称(圆形)分布 $P ( m_{x} , m_{y} )$。引入这些与 n 相关的因子是为了使 $T \to \infty$ 时 $U_{2} \to 0$,因此在这种情况下应使用 $n = 2$。这在低温下也是正确的(以保证 $U_{2} \to 1$),因为无论序参量结构如何,在有序相中,当系统尺寸发散时 Binder 比 $R_{2} \to 1$。图20展示了两个耦合比率 $g = 0.55$ 和 $0.70$ 下作为温度函数的结果。对于 $g = 0.55$,在 $L \ge 8$ 时可以观察到负的 Binder 累积量,负峰变得非常狭窄,并且随着系统尺寸的增加明显发散。在 $g = 0.70$ 时,相变仍然是一级的,但间断性较弱,仅在尺寸约为 $L = 32$ 时才开始表现为共存和负的 Binder 累积量。在这两种情况下,蒙特卡洛模拟仍然是遍历的,但必须进行大量的更新步数(每种情形约 $10^{8}$ 步)才能获得平滑的曲线。

图20. 阻挫伊辛模型在耦合比率 $g = 0.55$(左)和 $g = 0.70$(右)下向条纹相转变时的 Binder 累积量。随着 L 增大而发展出的发散负峰是一级相变的明确信号。

间断性与有限尺寸标度。在那些不够强的一级相变中,如上文系统中的情况,精确提取物理可观测量的间断性大小可能并不容易。图21展示了一个例子,其中绘制了几个晶格尺寸下内能随温度的变化关系。虽然随着 L 的增大,间断性(潜热)明显发展,但很难精确确定最终跳跃会发生在哪两个能量值之间。这需要利用更大晶格进行仔细分析。

图21还展示了两个不同耦合比率下比热最大值随系统尺寸的变化。在 $g = 0.51$ 时,相变是相当强的一级相变,对于最大的晶格可以观测到符合预期的标度 $C_{\mathrm{max}} \sim L^{2}$。相比之下,在 $g = 0.55$ 时,直到 $L \approx 64$ 为止,行为似乎遵循着不同的幂律(指数接近 1.2)。然而,对于更大的晶格,结果开始显示出稍快的发散趋势,并且最终,对于非常大的晶格,可以预期在此情形下指数也会等于 $d = 2$。

3.5 自旋刚度与Kosterlitz-Thouless相变#

具有长程磁有序系统的一个重要方面是它表现出非零的自旋刚度。对于具有连续矢量自旋的系统(XY模型或海森堡模型),自旋刚度类似于固体的弹性模量。它也常被称为螺旋模量(helicity modulus)[131]。为了研究这一概念,我们这里考虑XY模型,其哈密顿量写为:

$$ H = - J \sum_{\langle i , j \rangle} \cos ( \Theta_{i} - \Theta_{j} ) ,\tag{81} $$

其中 $\Theta_{i} \in [ 0 , 2 \pi ]$ 是描述自旋 i 的角度。为简单起见,我们只假设最近邻相互作用,但推广到任意相互作用是直接的。在推导自旋刚度的表达式时,我们将考虑二维正方晶格(有时为了简化,使用一维链),但稍后也会研究三维立方系统。在二维 XY 系统中,自旋刚度是表征该系统所表现出的非常规(拓扑)Kosterlitz-Thouless相变的重要物理量,我们也将在此节简要讨论它。

3.5.1 自旋刚度的定义#

通俗地说,自旋刚度 $\rho_{s}$ 描述了有序自旋在响应于施加序参量方向调制的扰动时进行适应的趋势(与磁化率形成对比,后者衡量序参量在响应于施加在固定方向上的场时发生变化的趋势)。它在连续介质力学中类似于剪切模量(shear modulus),后者表征弹性物体形状变形的趋势(而压缩性对应于保持形状下的体积变化趋势)。自旋刚度的定义在 $T = 0$ 时最容易理解,我们首先考虑这种情况,然后再推广到 $T > 0$。

$T = 0$ 时的刚度。我们首先考虑一个在 x 方向具有开放边界、在 y 方向具有周期性边界的系统(稍后我们将推广到 x 方向也是周期性边界的情况)。图22 说明了当左右边界之间施加自旋角度 $\Phi$ 的整体变化时(例如,由于在边界列施加的强磁场),$T = 0$ 时的铁磁 XY 自旋如何进行调整以最小化能量。为了使能量最小化,每一列相对于下一列扭转一个角度 $\phi = \Phi / ( L_{x} - 1 )$,其中 $L_{x}$ 是系统在 x 方向上的长度(列数)。在 $T = 0$ 时,我们关心存在这种扭转时的能量,对于二维 XY 模型,它简单地由下式给出:

图22. 一个二维经典 XY 模型,在 x 方向具有相位扭转,由固定边界列的自旋处于相对角度 $\Phi$ 实现。为最小化能量,总扭转 $\Phi$ 被均匀分布,使得相邻列的自旋相对扭转角度 $\phi = \Phi / L_x$。

$$ E ( \phi ) = - J ( L_{x} - 1 ) L_{y} \cos ( \phi ) = E ( 0 ) + J ( L_{x} - 1 ) L_{y} [ 1 - \cos ( \phi ) ] .\tag{82} $$

对于小的 $\phi$,我们得到领先阶的 $E ( \phi ) - E ( 0 ) = ( J / 2 ) ( L_{x} - 1 ) L_{y} \phi^{2}$。受此结果启发,$T = 0$ 时的自旋刚度定义为:

$$ \rho_{s} = \frac{1} {N} \frac{\partial^{2} E ( \phi )} {\partial \phi^{2}} ,\tag{83} $$

对于大N(或者如果按相互作用x键的数量$L_{\mathrm{v}} (L_{x} - 1)$归一化,由于扭曲造成的能量代价分布在这些键上,则对任意N),我们有$\rho_{s} = J$。

通常考虑周期系统更为方便。为了推导这种情况下的自旋刚度表达式,首先在边界上施加一个相位扭曲。为简单起见,我们在一维链上计算,自旋角度为$\Theta_{x}, x = 0, \ldots, N - 1$,但该计算可以轻松推广到更高维度。每个键的相互作用能为$E_{x} = - J \cos ( \Theta_{x + 1} - \Theta_{x} )$,除了边界处,扭曲Φ对应于$E_{N} = - J \cos ( \Theta_{0} - \Theta_{N - 1} + \Phi )$。最小化能量的构型满足$\Theta_{x + 1} - \Theta_{x} = \delta$,其中$\delta$与x无关,从而得到总能量

$$ E ( \Phi ) = - ( N - 1 ) \cos ( \delta ) - \cos ( \Phi - [ N - 1 ] \delta ) ,\tag{84} $$

该式通过取$\delta = \phi \equiv \Phi / N$ 最小化,得到$E ( \phi ) = E ( 0 ) + J N [ 1 - \cos ( \phi ) ]$,根据(83)定义的自旋刚度再次为$\rho_{s} = J$。

考虑另一种在周期系统中扭曲自旋的方式也很有用,即在哈密顿量中引入一个扭曲场 $\Phi_{x} = x \phi$。对于一维系统,在此场存在的情况下,每个键的能量为$- J \cos ( \Theta_{x + 1} - \Theta_{x} + \Phi_{x + 1} - \Phi_{x} )$。为了正确处理边界,应将$\Phi_{x}$视为连续变量x的函数,它在$x = N$处从$N \phi$不连续地跳变到0。那么XY相互作用中出现的相位差应理解为

$$ \Phi_{x + 1} - \Phi_{x} = \int_{x}^{x + 1} \Phi_{x} d x = \int_{x}^{x + 1} \phi d x = \phi ,\tag{85} $$

该关系在边界处$( x + 1 = N )$也成立。除去因子,$\phi$类似于穿过环的通量。在此扭曲存在的情况下,当$\Theta_{x} = \Theta$与x无关时能量最小化,得到$E ( \phi ) = E ( 0 ) + J N [ 1 - \cos ( \phi ) ]$,与之前扭曲边界条件的情况相同。可以将扭曲场视为一种变换掉扭曲边界条件(旋转到存在边界扭曲的自旋参考系)的方法,这在实践中使许多计算更容易。从现在开始,我们将使用扭曲场。

$T > 0$ 时的刚度 现在我们考虑非零温度,此时自旋刚度定义为

$$ \rho_{s} = \frac{1} {N} \frac{\partial^{2} {\cal F} ( \phi )} {\partial \phi^{2}} ,\tag{86} $$

其中 $F ( \phi )$ 是在扭曲场(或等价地,扭曲边界条件)存在下的自由能,它与配分函数的关系为

$$ F ( \phi ) = - \frac{1} {\beta} \ln [ Z ( \phi ) ] .\tag{87} $$

当 $T \to 0$ 时,方程 (86) 显然退化为基态能量的定义 (83)。

在二维XY模型中施加一个仅依赖于x的扭曲场 $\Phi ( x , y ) = x \phi$,我们得到哈密顿量

$$ H ( \phi ) = - J \sum_{\langle i , j \rangle_{x}} \cos ( \phi + \Theta_{j} - \Theta_{i} ) - J \sum_{\langle i , j \rangle_{y}} \cos ( \Theta_{j} - \Theta_{i} ) ,\tag{88} $$

其中我们可以假设两个方向都采用周期边界条件(虽然原则上y边界可以是开放的)。我们可以通过使用标准三角恒等式简化对$\phi$的依赖关系,

$$ \cos ( \phi + \Theta_{j} - \Theta_{i} ) = \cos ( \Theta_{j} - \Theta_{i} ) \cos ( \phi ) - \sin ( \Theta_{j} - \Theta_{i} ) \sin ( \phi ) ,\tag{89} $$

我们将其展开到$\phi$的二阶项

在存在小$( \phi \to 0 )$扭曲的情况下的哈密顿量可以写为

$$ \begin{array} {r} {H ( \phi ) H ( 0 ) + \frac{1} {2} \phi^{2} H_{x} - \phi I_{x} ,} \end{array}\tag{91} $$

其中$H_{x}$是$\phi = 0$时哈密顿量的x键部分,而$I_{x}$是x晶格方向上的总自旋“流”:

当存在一个小扭曲时,配分函数可以写为

其中包含 $H_{x}$ 和 $I_{x}$ 的指数项被泰勒展开到所需阶数。我们现在可以将其写成对于 $\phi = 0$ 分布的期望值形式。到 $\phi$ 的二次项:

$$ Z ( \phi ) Z ( 0 ) [ 1 - {\textstyle{\frac{1} {2}}} \beta \phi^{2} \langle H_{x} \rangle + \beta \phi \langle I_{x} \rangle + {\textstyle{\frac{1} {2}}} \beta^{2} \phi^{2} \langle I_{x}^{2} \rangle ] .\tag{94} $$

由对称性 $\langle I_{x} \rangle = 0$,自由能 (87) 由下式给出:

$$ \begin{array} {r} {F ( \phi ) = F ( 0 ) + \frac{1} {2} \phi^{2} \big ( \langle H_{x} \rangle - \beta \langle I_{x}^{2} \rangle \big ) ,} \end{array}\tag{95} $$

由此我们可以推导出自旋刚度 (86) 的简单表达式:

$$ \rho_{s} = \frac{1} {N} ( \langle H_{x} \rangle - \beta \langle I_{x}^{2} \rangle ) ,\tag{96} $$

它可以通过蒙特卡洛模拟来评估。这个结果是针对 $x$ 晶格方向上的扭曲场。对于一个 $d$ 维各向同性系统,我们可以对所有等效方向进行平均,并将自旋刚度写为

$$ \rho_{s} = \frac{1} {N d} \left( \langle H \rangle - \beta \sum_{a = 1}^{d} \langle I_{a}^{2} \rangle \right) ,\tag{97} $$

其中指标 $a$ 对应于晶格方向 $a$ 上的电流。对于一个各向异性的 $d$ 维系统,所有晶格方向的刚度通常是不同的,即存在 $d$ 个不同的刚度常数,每个都由类似 (96) 的形式给出。

扭曲场定义的刚度可以很容易地推广到具有长程相互作用的系统。形式 (96) 和 (97) 仍然适用,其中电流 $I_{x}$ 正好如同哈密顿量中那样包含了所有相互作用的贡献。

与超流性和超导性的关系。XY 模型自旋刚度的一个非常有趣的方面是,它可以直接与超导体(或超流体,如 $\mathrm{^4He}$)的超流密度相关联 [131, 132]。像 XY 模型的磁化强度一样,超导体或超流体的序参量是 $U ( 1 )$ 对称的(对应于波函数的全局相位),而我们上面讨论的扭曲场直接类似于磁通量(在超导体的情况下)。非零刚度对应于超导体表现出的迈斯纳效应。旨在研究超薄膜超流体性质的 2D XY 模型的蒙特卡洛模拟在参考文献 [133, 134] 中讨论。

图23. (a) 不同尺寸 $N = L^{3}$ 格点上 3D XY 模型的自旋刚度随温度的变化。(b) 根据方程 (99) 的有限尺寸标度。不同 $L$ 的 $L \rho_{s}$ 曲线在临界温度处(对较大的 $L$ 渐近地)相互交叉。这些结果是使用蒙特卡洛簇算法 [119] 获得的,在所有情况下误差棒都小到不可见。

自旋刚度标度。自旋(或超流)刚度的临界标度性质最初是在超流体的背景下提出的 [131, 132]。在 $d > 2$ 的无限 $d$ 维系统中,已证明当 $T$ 从下方趋近 $T_{\mathrm{c}}$ 时,有

$$ \rho_{s} \sim ( T_{c} - T )^{( d - 2 ) \nu} ,\tag{98} $$

其中 $\nu$ 是标准关联长度指数。根据一般有限尺寸标度关系 (66),恰好处于 $T_{c}$ 时刚度的尺寸依赖性由下式给出:

$$ \rho_{s} \sim L^{2 - d} .\tag{99} $$

因此,自旋刚度与关联长度和 Binder 比一样,是确定临界点而不必调整任何未知(或未精确已知)指数的有用量。图 23 中的 3D XY 模型蒙特卡洛结果说明了这一点,其中在临界点处 $L \rho_{s}$ 与尺寸无关。

3.5.2 二维中的 Kosterlitz-Thouless 相变#

根据 Mermin-Wagner 定理 [38](如第 2 节所述),在 2D XY 模型中,在 $T > 0$ 时不可能转变到具有长程磁有序的相。值得注意的是,该系统展现出一种不同类型的相变,其中没有长程序形成,但自旋关联从指数衰减转变为幂律形式 [135, 136]。这种 Kosterlitz-Thouless(KT)相变本质上是拓扑的,是由于在温度 $T > T_{\mathrm{KT}}$ 时,自旋构型中非束缚涡旋(即拓扑缺陷)大量产生的结果。对于 $T \le T_{\mathrm{KT}}$,涡旋也存在(其密度在 $T \to 0$ 时趋于零),但它们都被束缚在涡旋-反涡旋对中,这些涡旋对没有净涡度(因此在自旋的粗粒化过程中消失)。自旋关联的幂律形式 $C ( r ) \sim r^{- \eta}$ 适用于所有 $0 < T \le T_{\mathrm{KT}}$。指数依赖于温度,在 $T_{\mathrm{KT}}$ 处值为 $\eta = 1 / 4$,而当 $T \to 0$ 时 $\eta \to 0$,因此真正的长程有序存在于 $T = 0$。

尽管对于 $0 < T \leq T_{\mathrm{KT}}$ 没有长程序,但自旋刚度在 KT 相中实际上非零——具有足够小指数的幂律关联足以支撑边界扭转的能量代价。刚度没有幂律形式的出现,而是有一个更显著的相变信号;在 $T_{\mathrm{KT}}$ 处从 $\rho_{s} = 0$ 到非零值的不连续跳跃。对应于 2D XY 模型的连续场论的重整化群计算给出了关于 KT 相变的非常详细的信息。最重要的结果之一是相变温度 $T_{\mathrm{KT}}$ 与恰好在此温度下的自旋刚度(即间断的大小)之间的严格关系(由 Nelson 和 Kosterlitz [139] 提出):

$$ \rho_{s} ( T_{\mathrm{KT}} ) = \frac{2 T_{\mathrm{KT}}} {\pi} .\tag{100} $$

对于有限晶格,$T_{\mathrm{KT}}$ 处的刚度以对数尺寸修正 [138] 趋近于无穷大尺寸的值:

$$ \rho_{s} ( T_{\mathrm{KT}} , L ) = \rho_{s} ( T_{\mathrm{KT}} , \infty ) \left( 1 + \frac{1} {2 \ln ( L ) + c} \right) ,\tag{101} $$

其中 $c$ 是一个系统依赖的参数。因此,可以说,在临界点处 $\rho_{s}$ 的一般有限尺寸标度律 (99) 成立,但当 $d = 2$ 时,对于主导阶的尺寸无关形式存在对数修正。

图 24(a) 显示了 2D XY 模型刚度随温度变化的蒙特卡洛结果。热力学极限中预期的跳跃作为尺寸的函数趋近得非常缓慢。这与以下事实有关:$T > T_{\mathrm{KT}}$ 的关联长度并非以幂律发散,而是按照指数形式:

$$ \xi \sim \mathrm{e}^{a / ( T - T_{\mathrm{KT}} )^{1 / 2}} ,\tag{102} $$

其中 $a$ 取决于晶格和相互作用的细节。因此,使用与我们在第 3.3.2 节标准临界点情况中相同的论点,$T_{\mathrm{KT}}$ 的有限尺寸偏移 [例如,使用 $\rho_{s}$ 下降最快的温度 $T^{*} ( L )$ 来定义] 由以下形式给出:

$$ T^{*} ( L ) - T_{\mathrm{KT}} \sim \frac{1} {\ln^{2} ( L )} ,\tag{103} $$

这比传统的幂律偏移 $T^{*} ( L ) - T_{c} \sim L^{1 / \nu}$ 慢得多。

使用标准有限尺寸标度假说 (64),并用 (101) 中的对数尺寸修正替换 $L^{\sigma}$(临界点处奇异量的尺寸依赖性),我们可以将 KT 相变处自旋刚度的尺寸和温度依赖性假设写为:

$$ \rho_{s} ( T , L ) = \left( 1 + \frac{1} {2 \ln ( L ) + c} \right) f ( e^{a / ( T - T_{\mathrm{KT}} )^{1 / 2}} / L ) ,\tag{104} $$

对尺度函数$f ( x )$的自变量取对数后,可将其用另一函数$g [ \ln ( x ) ]$表示为

$$ \rho_{s} ( T , L ) \left( 1 + \frac{1} {2 \ln ( L ) + c} \right)^{- 1} = g [ \ln ( L ) - a / ( T - T_{\mathrm{KT}} )^{1 / 2} ] .\tag{105} $$

图24. (a) 不同尺寸$L = 2^{n}$的格点上,二维XY模型自旋刚度温度依赖关系的蒙特卡洛结果(使用簇算法[119]获得)。随着$L \to \infty$,在$T_{\mathrm{KT}}$处出现一个不连续性,如已知转变温度处的竖线所示($T_{\mathrm{KT}} \approx$ 0.8933 [137])。根据Nelson-Kosterlitz关系(100),对于任何呈现KT相变的系统,恰好处于$T_{\mathrm{KT}}$时的自旋刚度必须落在所示的直线上:$\rho_{s} = 2 T / \pi$。(b) 根据组合$( T , L )$标度假说(105)的有限尺寸数据坍缩,使用已知的$T_{\mathrm{KT}} = 0.8933$,并选择两个参数$a = 1.5$和$c = 0.7$以获得良好的数据坍缩。竖线显示了根据Nelson-Kosterlitz关系预期的渐近$T T_{\mathrm{KT}} , L \to \infty$的值。

此处$T T_{\mathrm{KT}}$对应于自变量$z = \ln ( L ) - a / ( T - T_{\mathrm{KT}} )^{1 / 2} \to - \infty$。

在许多研究中,例如[133, 137],二维XY模型的KT相变温度以不同方式被提取出来。图24(b)展示了标度形式(105)的检验,使用了[137]中获得的$T_{\mathrm{KT}} = 0.8933$,并调整常数a和c,以使尺寸$L = 8 , 16 , \dots , 128$的系统数据(近似)最优地坍缩到一个共同的尺度函数上。该形式确实很好地描述了数据。这种图同时证实了KT相变的三个不同方面(这些方面在KT相变历史的不同阶段从理论上推导得出):指数发散的关联长度(102) [135, 136]、对数修正(101) [138]以及Nelson-Kosterlitz关系(100) [139]。

3.6 量子相变#

在本笔记的后续章节中,我们将研究几个量子相变的例子,这些相变发生在系统的基态中,并随某个模型参数变化[1]。我们之前讨论的经典连续相变和一级相变的标度性质仍然适用,但需要一些重要的扩展和修正。经典热相变的性质可以追溯到自由能中的奇异性。而在量子相变中,则是基态能量表现出奇异行为,这也在其他物理量中有所体现。这可以理解为源于量子系统的$T > 0$自由能,当$T \to 0$时,该自由能变为基态能量。这在经典情况下当然也是如此,但经典模型的基态通常不会随参数连续演化(正如我们在第3.4.2节阻挫伊辛模型的例子中看到的那样),因此表现为强一级相变。相比之下,正如第2.4节中已经讨论的,量子系统具有非平凡的基态,其零点涨落随参数变化而变化。由发散的量子涨落驱动的连续相变很常见。

正如我们在第5节讨论量子蒙特卡洛方法时将会看到的那样,使用路径积分可以将一个$d$维量子系统形式上映射为一个等价的$d+1$维经典统计力学问题(尽管有时其分布函数并非正定)。在新产生的“虚时”维度上,系统的大小为温度的倒数:$L_{\tau} = c / T$,其中c为速度。当$T > 0$时,这一维度是有限的,而空间维度可以是无限的。此时系统是$d$维的,具有一个“厚度”$L_{\tau}$。当$T \to 0$时,系统才严格具有$(d+1)$维性。在这种情况下,哈密顿量中的一个可调参数所扮演的角色与经典系统中的温度非常相似。有趣的是,改变这里的温度类似于在$L_{\tau}$上的有限尺寸标度[5],这可以用来推导出量子临界点附近的高温标度性质[84]。

在某些情况下,$d$维量子系统具有与$d+1$维同类型经典系统相同的低能性质。例如,第2.4节讨论的二聚化海森堡模型便是如此。这些模型的低能物理可以映射到一个3D经典海森堡模型,这通常是通过连续场论(如非线性σ模型[5, 84])来完成的(需要注意的是,二聚化仅仅是调节量子涨落强度的一种手段,在粗粒化的均匀有效模型中,它仅仅由一个场论的耦合常数,或是一个有效3D均匀经典模型中的温度来表示[140])。因此可以推断,由调节二聚化强度驱动的量子相变应属于3D经典海森堡模型温度驱动相变的普适类。在其他情况下,低能映射可能会给出一个与任何已知经典模型都不对应的有效系统,但仍然可以认为该系统对应于某个$(d+1)$维有效模型。

动力学临界指数。在许多情况下,例如上述二聚化海森堡模型,在映射到$d+1$维时产生的时间维度(在渐近意义下)与空间维度是等价的。那么,在量子相变处,此类系统的标度性质只需将临界关联函数(53)、超标度关系(如第二个方程(59))以及自旋刚度的标度形式(98)和(99)中的$d$替换为$d+1$即可得到。在另一些情况下,新维度中的关联可能与空间维度中的关联有本质区别。动力学指数$z$关联了空间和时间的关联所对应的幂律关系:例如,当某个参数$g$被调谐到其临界值时,空间关联长度$\xi$的发散形式为$\lvert g - g_{c} \rvert^{- 1 / \nu}$,则时间关联长度的发散形式为$| g - g_{c} | ^{z / \nu}$;在经典标度关系中,对于量子临界点,应将$d$替换为$d+z$。

动力学指数得名于它同时支配着波数为$q$的激发元的色散关系$\omega ( q ) \sim q^{z}$。这一点的一个重要方面是,通过设定$q \propto 1 / L$可得到有限尺寸激发能隙,其标度形式为$\Delta_{L} \sim 1 / L^{z}$。这一结果可直接用于数值计算,并且对依赖于激发谱的量(例如各种磁化率)的标度性质也有间接影响。提取动力学指数是量子相变计算研究中的一个重要方面。除此之外,连续和一级量子相变都可以使用与上述经典相变相同的有限尺寸标度方法进行分析。

需要指出的是,经典系统中也存在动力学指数,但这在平衡态统计力学中并不发挥作用(因为如果哈密顿量包含动能部分,与之相关的相空间概率分布会在期望值计算中相互抵消)。经典动力学指数取决于如何将动力学引入系统[141]。伊辛模型(Ising model)的Metropolis蒙特卡洛算法是经典动力学(称为Glauber动力学)的一个例子,临界点处自关联的长时间衰减由一个相关的动力学指数控制。然而,这种动力学并不决定平衡性质。在量子力学中,动力学是无法分离的,它是平衡性质的一个不可分割的组成部分,因为完整的哈密顿量始终参与其中,并同时包含了静态和动力学性质[142]。