MC-06 量子王-朗道

本教程的目标是用量子王-朗道（QWL）算法计算量子自旋系统的热力学性质。 与在固定温度下模拟的路径积分量子蒙特卡洛方法不同，QWL 在配分函数展开阶数的空间中做随机行走，从而直接估计态密度。 这样，只需一次模拟就可以在整个温度范围内计算各种热力学可观测量——能量、比热、熵、磁化率——而无需重新运行。

量子海森堡自旋链的热力学

铁磁海森堡链

参数文件 parm6a 设置了 40 个格点的链上 $S=\frac{1}{2}$ 海森堡铁磁体的 QWL 模拟：

LATTICE="chain lattice"
MODEL="spin"
local_S = 1/2
L       = 40
CUTOFF  = 500
{J = -1}

CUTOFF 设定算法所抽样的配分函数展开阶数的上限；对于本教程所研究的温度范围内的 40 格点链，取 500 已经足够。

命令行

准备好输入文件并运行 qwl 程序之后：

parameter2xml parm6a
qwl parm6a.in.xml

用 qwl_evaluate 为所有热力学可观测量生成 XML 绘图文件：

qwl_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6a.task1.out.xml

这会为每个可观测量生成一个 XML 文件：

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_structure_factor.xml
parm6a.task1.plot.staggered_structure_factor.xml
parm6a.task1.plot.uniform_susceptibility.xml

用 plot2text 把数据提取为纯文本：

plot2text parm6a.task1.plot.energy.xml

工具 plot2xmgr（Grace）和 plot2gp（Gnuplot）会以相应的格式给出等价的输出。分析时推荐使用 Python，见下文。

Python

脚本 tutorial6a.py 设置并运行模拟，然后一次调用即可计算所有可观测量：

import pyalps
import matplotlib.pyplot as plt

data = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6a'),
                          DELTA_T=0.1, T_MIN=0.1, T_MAX=10.0)

for s in pyalps.flatten(data):
    plt.figure()
    plt.title("Ferromagnetic Heisenberg chain")
    pyalps.plot.plot(s)
    plt.show()

对于铁磁体（$J=-1$），你应当看到低温处有一个宽阔的比热峰（类似肖特基反常），以及在 $T\to 0$ 时发散的均匀磁化率，这与一维铁磁序在零温下形成的图像一致。

反铁磁海森堡链

要模拟反铁磁链，把 $J=1$ 而不是 $J=-1$。 参数见 parm6b，Python 脚本见 tutorial6b.py；其余一切与铁磁的情形完全相同。

对于反铁磁体，均匀磁化率在 $T\to 0$ 时保持有限（这是一维反铁磁体自旋液体基态的标志），而比热峰的位置和展宽方式则有所不同。

三维海森堡反铁磁体

模拟三维量子海森堡反铁磁体

参数文件 parm6c 设置了在含 $4^3=64$ 个格点的简单立方格子上 $S=\frac{1}{2}$ 海森堡反铁磁体的 QWL 模拟。 Python 脚本是 tutorial6c.py。 运行和计算的方式与上面的链完全相同。

交错结构因子 $S(\pi,\pi,\pi)$ 应当在 $T\approx 1$ 附近开始陡峭上升，标志着反铁磁关联的出现。 比热在相应位置出现峰值，而熵在同一温度区间内迅速下降。

用有限尺寸标度分析确定临界点

有限尺寸标度预言，临界点处的交错结构因子按 $S(L) \propto L^{2-\eta}$ 标度，其中 $\eta \approx 0.034$（三维经典海森堡普适类）。 把 $S(L)/L^{2-\eta}$ 对温度作图，不同 $L$ 的曲线应当在临界温度 $T_c$ 处相交。

参数文件 parm6d（或 tutorial6d.py）在两种系统尺寸（$L=4$ 和 $L=6$）下运行立方格子反铁磁体，并采用更大的 CUTOFF=1000 以保证低温下的精度。 运行结束后载入结果：

import pyalps
import matplotlib.pyplot as plt
import copy

results = pyalps.evaluateQWL(pyalps.getResultFiles(prefix='parm6d'),
                             DELTA_T=0.05, T_MIN=0.5, T_MAX=1.5)

提取每个系统尺寸的交错结构因子，按 $L^{-(2-\eta)}$ 重新标度，并用 $L$ 加以标注：

eta = 0.034
data = []
for s in pyalps.flatten(results):
    if s.props['ylabel'] == 'Staggered Structure Factor per Site':
        d = copy.deepcopy(s)
        l = s.props['L']
        d.props['label'] = 'L=' + str(l)
        d.y = d.y * pow(float(l), -(2 - eta))  # 按 L^{-(2-eta)} 重新标度
        data.append(d)

然后画出标度曲线：

plt.figure()
plt.title("Scaling plot for cubic lattice Heisenberg antiferromagnet")
pyalps.plot.plot(data)
plt.legend()
plt.xlabel('Temperature $T/J$')
plt.ylabel('$S(\pi,\pi,\pi)\, L^{-(2-\eta)}$')
plt.show()

不同 $L$ 的曲线应当在 $T_c \approx 0.946$ 附近相交。

思考题

• 比较铁磁（$J=-1$）与反铁磁（$J=1$）链的比热和均匀磁化率。差别在哪里最为显著？为什么两者在高温下几乎完全相同？
• 链在 $T=0$ 以及 $T\to\infty$ 时的熵各是多少？（如果 $L=40$ 的数据在低 $T$ 下噪声较大，可以改用 $L=8$ 重做。）这与热力学第三定律一致吗？
• 为什么均匀磁化率在铁磁体中于 $T\to 0$ 时发散，而在反铁磁体中却保持有限？
• 为什么三维反铁磁体的交错结构因子在 $T\approx 1$ 附近开始增大？还有哪些热力学特征伴随着这一变化？
• 重新标度后的结构因子曲线 $S(L)\,L^{-(2-\eta)}$ 是否交于同一个温度？你对 $T_c$ 的估计是多少？与文献值 $T_c = 0.946$ 相比如何？
• 怎样才能得到更精确的 $T_c$ 估计？需要在模拟中改变什么？
• 你认为立方格子上量子海森堡铁磁体的临界温度会与反铁磁体的相同吗？该如何检验？