MC-08 量子相变

本教程的目标是利用量子蒙特卡洛数据的有限尺寸标度，在一个二维自旋模型中定位并刻画量子相变。与经典相变不同，量子相变发生在 $T=0$ ，其驱动力是量子涨落而非热涨落。我们将用 Binder 累积量和自旋刚度来确定临界耦合，然后提取关联长度指数 $\nu$ 、反常量纲 $\eta$ 和动力学临界指数 $z$ ，从而确定该相变所属的普适类。事实证明该模型有两个量子临界点，我们对二者都将进行研究。

该模型由排布在正方晶格上的自旋 $\frac{1}{2}$ 海森堡梯子构成。每条梯子沿水平方向延伸：梯子内部的格点沿腿以耦合 $J_0$ 相连，沿横档以耦合 $J_1$ 相连。相邻梯子之间沿竖直方向以 $J_2$ 耦合：

o-J0-o-J0-o-J0-o   <- 腿
|         |
J1        J1        <- 横档（梯子内部）
|         |
o-J0-o-J0-o-J0-o   <- 腿
|         |
J2        J2        <- 梯子间耦合
|         |
o-J0-o-J0-o-J0-o   <- 腿
|         |
J1        J1        <- 横档（梯子内部）
|         |
o-J0-o-J0-o-J0-o   <- 腿

参数 L 给出每条腿上的格点数；W = L/2 给出行数，因此总共有 W/2 条梯子。在本教程中我们取 $J_0=J_1=1$ ，并改变梯子间耦合 $J_2$ （图示可参见 Wenzel and Janke, Phys. Rev. B 79, 014410 (2009) 的图 1）。尽管二维海森堡模型在有限温度下不存在相变（Mermin-Wagner 定理），但不同基态之间的相变可以发生在 $T=0$ 。

后台模拟

跨多个系统尺寸的有限尺寸标度运行，比下面针对单一系统的扫描耗时要长得多。现在就把它启动起来，让它在后台运行，同时你继续学习本教程的其余部分。

命令行

下载 parm8b 并运行：

parameter2xml parm8b
loop parm8b.in.xml &

Python

在继续之前，请在另一个终端中或以后台进程的方式运行 tutorial8b.py 的第一部分（设置部分与 pyalps.runApplication 调用）。

识别不同的相

我们先考察两个简单的极限：退耦合的梯子（ $J_2=0$ ）和各向同性的正方晶格（ $J_2=1$ ）。退耦合梯子的基态具有短程关联并表现出有限的自旋能隙：这是一个自旋液体相。相比之下，正方晶格表现出长程序，具有有限的交错磁化强度：这是一个反铁磁奈尔相。

探测这两个相的一个清晰办法是考察磁化率 $\chi$ 。在两种情形下都对一个 $8\times 8$ 的系统在一系列温度上进行模拟并比较结果。对于退耦合梯子，由于自旋能隙的存在，磁化率在低温下表现出激活行为；而在正方晶格上，当 $T\to 0$ 时磁化率趋于一个有限常数。请注意，在任何有限系统上，由于有限尺寸能隙的存在， $\chi$ 在足够低的温度下最终都会趋于零，但这不是我们这里关注的重点。你可以用参数文件 parm8a 在命令行上运行模拟：

parameter2xml parm8a
loop parm8a.in.xml

或者使用 Python 脚本 tutorial8a.py。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import pyalps.fit_wrapper as fw

parms = []
for j2 in [0.,1.]:
    for t in [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]:
        parms.append(
            { 
              'LATTICE'        : "coupled ladders", 
              'LATTICE_LIBRARY': 'lattices.xml',  # 自定义晶格定义的路径
              'MODEL_LIBRARY'  : 'models.xml',    # 自定义模型定义的路径
              'local_S'        : 0.5,
              'ALGORITHM'      : 'loop',
              'SEED'           : 0,
              'T'              : t,
              'J0'             : 1 ,
              'J1'             : 1,
              'J2'             : j2,
              'THERMALIZATION' : 5000,
              'SWEEPS'         : 50000, 
              'MODEL'          : "spin",
              'L'              : 8,
              'W'              : 4
            }
    )
    
input_file = pyalps.writeInputFiles('parm8a',parms)
pyalps.runApplication('loop',input_file)

对于 $J_2=0$ ，自旋能隙的数值可以由磁化率的有限温度行为估计得到（推导见 Phys. Rev. B 50, 13515 (1994)）：

\chi = \frac{A}{\sqrt{T}} \exp\!\left(-\frac{\Delta}{T}\right),

其中 $A$ 和自旋能隙 $\Delta$ 是拟合参数。对 $T\leq1$ 的数据进行拟合，提取自旋能隙的估计值。下面是在 Python 中执行该分析的一个示例：

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8a.task*.out.h5'),['Staggered Susceptibility','Susceptibility'])
susc1=pyalps.collectXY(data,x='T',y='Susceptibility', foreach=['J2'])

lines = []
gap_j0 = None
for data in susc1:
    pars = [fw.Parameter(1), fw.Parameter(1)]
    data.y= data.y[data.x < 1]
    data.x= data.x[data.x < 1]
    f = lambda self, x, pars: (pars[0]()/np.sqrt(x))*np.exp(-pars[1]()/x)
    fw.fit(None, f, pars, np.array([v.mean for v in data.y]), data.x)
    prefactor = pars[0].get()
    gap = pars[1].get()
    if data.props['J2'] == 0.0:
        gap_j0 = gap
    
    lines += plt.plot(data.x, f(None, data.x, pars))
    lines[-1].set_label('$J_2=%.4s$: $\\chi = \\frac{%.4s}{\\sqrt{T}}\\exp\\left(\\frac{-%.4s}{T}\\right)$' % (data.props['J2'], prefactor,gap))

拟合给出能隙 $\Delta$ 和前置因子 $A$ ；所得曲线如下所示。

plt.figure()
pyalps.plot.plot(susc1)
plt.xlabel(r'$T$')
plt.ylabel(r'$\chi$')
plt.title('spin gap $\\Delta \\approx$ %.4s' % gap_j0)
plt.legend()
plt.show()

对于 $J_2=0$ ，磁化率在低温处应急剧下降，拟合曲线应与数据吻合良好，给出自旋能隙 $\Delta \approx 0.5$ 。对于 $J_2=1$ ，当 $T \to 0$ 时磁化率应趋于一个有限常数，反映出正方晶格上不存在能隙。

定位相变

既然我们在 $J_2=0$ 和 $J_2=1$ 处识别出了两个不同的相，那么必然存在（至少）一个量子相变将它们分开。上面启动的后台运行会在耦合区间 $J_2 \in [0.2,0.4]$ 上，对系统尺寸 $L=8,10,12,16$ 、逆温度 $\beta=2L$ 进行扫描（完整参数集见 parm8b / tutorial8b.py）。它运行结束后，按下面的说明载入并分析结果。

逆温度的选择

这些模拟是在有限逆温度 $\beta = 2L$ 下进行的，而不是严格地在 $T=0$ 下进行。 $\beta \propto L$ 这一选择并非随意：在该量子相变处动力学临界指数为 $z=1$ ，这意味着时间和空间以相同的方式标度。取 $\beta = 2L$ 保证了路径积分的虚时长度随系统尺寸按相同比例增长，因而模拟采样到的是基态物理而非热激发。当你观察到 $\rho_s L$ 曲线交于同一点时，就隐含地验证了这一选择——只有当 $z=1$ 且 $\beta \propto L^z$ 时，这一交点才会干净利落地出现。

还要注意，圈算法的计算时间随 $\beta$ 线性增长：由于任何有限温度 QMC 算法的标度都不可能优于时空体积 $\beta L^d$ 的线性关系，因此即便在量子临界点处，这一标度也是最优的。

为确认结果不受有限温度影响，把 $\beta=2L$ 改成 $\beta=4L$ 并重复模拟——自旋刚度和 Binder 累积量应保持不变。也可以试试 $\beta=L/4$ ，观察当系统远离基态时结果如何变差。

交错磁化强度、Binder 累积量与自旋刚度

与经典蒙特卡洛教程中一样，我们用两个可观测量来确定相变位置。

第一个是交错磁化强度 $m_s$ 的 Binder 累积量， $m_s$ 是反铁磁相的序参量：

U_4 = \frac{\langle m_s^4\rangle}{\langle m_s^2\rangle^2}.

第二个是自旋刚度（Wenzel and Janke, Phys. Rev. B 79, 014410 (2009)），对该模型而言它比 Binder 累积量的交点具有更小的有限尺寸修正：

\rho_s = \frac{3}{4\beta} \langle w_x^2 + w_y^2\rangle,

其中 $w_x$ 、 $w_y$ 是世界线沿 $x$ 与 $y$ 方向的绕数。在量子临界点处 $\rho_s \propto L^{d-2-z}$ ，其中 $d$ 是空间维数， $z$ 是动力学临界指数。当 $z=1$ 时，组合量 $\rho_s L$ 在临界处是无量纲的，因此不同 $L$ 的曲线全都交于 $J_2^c$ 。两个可观测量都出现交点——而不是其中之一发散——表明这是一个连续相变而非一级相变。

用以下命令载入并绘制这些可观测量：

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm8b.task*.out.h5'),['Binder Ratio of Staggered Magnetization','Stiffness'])

binder=pyalps.collectXY(data,x='J2',y='Binder Ratio of Staggered Magnetization', foreach=['L'])
stiffness =pyalps.collectXY(data,x='J2',y='Stiffness', foreach=['L'])

for q in stiffness:
    q.y = q.y*q.props['L']

plt.figure()
pyalps.plot.plot(stiffness)
plt.xlabel(r'$J_2$')
plt.ylabel(r'Stiffness $\rho_s L$')
plt.title('coupled ladders')

plt.figure()
pyalps.plot.plot(binder)
plt.xlabel(r'$J_2$')
plt.ylabel(r'$g(m_s)$')
plt.title('coupled ladders')
plt.show()

不同 $L$ 的 $\rho_s L$ 曲线应从左侧的零值上升，并在某一点附近相交，给出 $J_2^c \approx 0.30$ – $0.31$ 的初步估计。 Binder 累积量曲线在同一区域相交，但有限尺寸修正更大，使得在这些系统尺寸下交点不够锐利。两个可观测量都相交（而不是其中之一发散）证实了该相变是连续的。

临界指数的估计

你已经得到了量子临界点 $J_2^c$ 的粗略估计。与经典情形一样，提取临界指数需要对 $J_2^c$ 作更精确的确定。

这可以通过在更精细的 $J_2$ 网格上运行更大的系统尺寸来实现，相应设置见 parm8d 和 tutorial8d.py。注意这些模拟非常耗费 CPU，因此留作练习。绘制不同系统尺寸下的 Binder 累积量 $U_4$ 和重标度后的刚度 $\rho_s L$ ；交点给出 $J_2^c$ 的更精确估计。为提取 $\nu$ ，考察这些量对 $J_2$ 的导数在 $J_2^c$ 处如何随系统尺寸标度。这些导数原则上可以在蒙特卡洛中直接测量，但对本教程而言，利用精细的 $J_2$ 网格作数值微分就已足够。

对两个量在不同系统尺寸下作数值微分，并把它们在 $J_2^c$ 处的值作为系统尺寸的函数画出。数据应按幂律标度：

\left.\frac{dU_4}{dJ_2}\right|_{J_2^c} \propto L\left.\frac{d\rho_s}{dJ_2}\right|_{J_2^c} \propto L^{1/\nu}.

用幂律拟合提取 $\nu$ 。两个量应给出彼此一致的估计，接近三维经典海森堡模型的值 $\nu \approx 0.71$ 。

练习： 与经典情形一样，你可以通过数据坍缩来直观判断估计值的质量。 $U_4$ 和 $\rho_s L$ 的标度形式与经典相变例子中 Binder 累积量的标度形式完全相同。

除 $z$ 和 $\nu$ 之外，最后一个独立指数 $\eta$ 可以由交错磁化率在临界点处的有限尺寸标度得到：

\chi_s(J_2^c) \sim L^{2-\eta}.

注意这里必须使用交错磁化率 $\chi_s$ ，即与序参量涨落相联系的磁化率，而不是均匀磁化率 $\chi$ 。把 $J_2^c$ 处的 $\chi_s$ 作为系统尺寸的函数画出，并从斜率提取 $\eta$ 。预期值为 $\eta \approx 0.034$ ，接近于零，这意味着交错磁化率在临界点处几乎按 $L^2$ 增长。

本教程中的两个量子相变都属于三维经典海森堡模型有限温度相变的普适类。

定位第二个量子临界点

当 $J_2$ 非常大时，梯子间耦合压过梯子内部的耦合 $J_0 = J_1 = 1$ ，系统重新组织成沿 $J_2$ 键的有效单重态二聚体。这同样是一个有能隙的自旋液体相，因此随着 $J_2$ 从 0 增大到很大的值，相图的结构为：自旋液体 → 奈尔反铁磁 → 自旋液体。因此必然存在第二个量子临界点 $J_2^{c_2}$ ，在该点奈尔序被破坏。

我们在参数区间 $J_2 \in [1.8, 2.1]$ 中重复有限尺寸标度分析，使用参数文件 parm8c 或脚本 tutorial8c.py，它们使用与 parm8b 相同的系统尺寸和 $\beta=2L$ ，但扫描更高的 $J_2$ 区间。

命令行

parameter2xml parm8c
loop parm8c.in.xml

Python

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np

parms = []
for l in [8,10,12,16]:
    for j2 in [1.8,1.85,1.9,1.95,2.,2.05,2.1]:
        parms.append(
            {
              'LATTICE'        : "coupled ladders",
              'LATTICE_LIBRARY': 'lattices.xml',
              'MODEL_LIBRARY'  : 'models.xml',
              'local_S'        : 0.5,
              'ALGORITHM'      : 'loop',
              'SEED'           : 0,
              'BETA'           : 2*l,
              'J0'             : 1,
              'J1'             : 1,
              'J2'             : j2,
              'THERMALIZATION' : 5000,
              'SWEEPS'         : 50000,
              'MODEL'          : "spin",
              'L'              : l,
              'W'              : l//2
            }
        )

input_file = pyalps.writeInputFiles('parm8c', parms)
pyalps.runApplication('loop', input_file)

用与之前相同的分析命令载入并绘制刚度 $\rho_s L$ 和 Binder 累积量，只需把文件模式中的 parm8b 换成 parm8c。

$\rho_s L$ 曲线应再次交于同一点，这次位于 $J_2^{c_2} \approx 1.91$ 附近。 Binder 累积量的交点也应出现在同一区域。第二个相变与第一个属于同一普适类——三维经典海森堡模型——因此同样的临界指数 $\nu$ 和 $\eta$ 适用。

问题

画出 $J_2=0$ 和 $J_2=1$ 时的磁化率。它们的低温行为有何不同？各自告诉了你关于基态的什么信息？
对 $J_2=0$ ，把 $\chi$ 拟合到 $A/\sqrt{T}\exp(-\Delta/T)$ 并提取自旋能隙 $\Delta$ 。你的估计与文献值（Phys. Rev. Lett. 73, 886 (1994); Phys. Rev. Lett. 77, 1865 (1996)）相比如何？
根据 $\rho_s L$ 曲线和 Binder 累积量曲线的交点，你对第一个量子临界点 $J_2^c$ 的估计是多少？
把 $\beta=2L$ 改成 $\beta=4L$ 并重复计算。自旋刚度和 Binder 累积量是否受到影响？再试试 $\beta=L/4$ ：结果是否有明显变化？这告诉了你温度在这些模拟中扮演什么角色？
由幂律标度 $dU_4/dJ_2|_{J_2^c} \propto L^{1/\nu}$ ，你得到的 $\nu$ 值是多少？
由交错磁化率的有限尺寸标度 $\chi_s(J_2^c) \sim L^{2-\eta}$ ，你得到的 $\eta$ 值是多少？
把你对 $\nu$ 和 $\eta$ 的估计与 Phys. Rev. B 65, 144520 (2002) 中报告的三维经典海森堡普适类的值作比较。
（进阶）用你估计的 $J_2^c$ 和 $\nu$ 对 $U_4$ 和 $\rho_s L$ 作数据坍缩。坍缩的质量对 $J_2^c$ 的精确取值有多敏感？
扫描 $J_2 \in [1.8, 2.1]$ 以定位第二个量子临界点。你的估计与 Wenzel and Janke, Phys. Rev. B 79, 014410 (2009) 的高精度结果相比如何？

MC-07 相变