跳至内容
教程资源
Exact Diagonalization
ED-05 ED Phase Transition

ED-05 ED Phase Transition

具有次近邻相互作用的海森堡链的临界点

在本教程中，我们将延续 ED-04 教程最后一部分关于海森堡链的讨论。我们将在哈密顿量中加入一个次近邻耦合项，

H=J1∑⟨i,j⟩Si⋅Sj  +  J2∑⟨⟨i,j⟩⟩Si⋅Sj, H = J_1 \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j \; + \; J_2 \sum_{\langle\langle i,j \rangle\rangle} \mathbf{S}_i \cdot \mathbf{S}_j ,

其中第一个求和遍历同一条周期链上的最近邻键，第二个求和遍历次近邻键。这就是著名的 J1J_1–J2J_2（或 Majumdar-Ghosh）链。

在 J2=0J_2 = 0 的极限下，该模型退化为临界海森堡链，可由贝特假设（Bethe ansatz）求解。在 J2/J1=0.5J_2/J_1=0.5 时，该模型也是精确可解的：此时体系处于 Majumdar-Ghosh 点，基态恰好是最近邻单态的直积（C.K. Majumdar and D.K. Ghosh, J. Math. Phys. 10, 1388 (1969)），

∣Ψ⟩=(∣↑⟩1∣↓⟩2−∣↓⟩1∣↑⟩2)(∣↑⟩3∣↓⟩4−∣↓⟩3∣↑⟩4)(∣↑⟩5∣↓⟩6−∣↓⟩5∣↑⟩6) |\Psi\rangle = \left(|\uparrow\rangle_1 |\downarrow\rangle_2 - |\downarrow\rangle_1 |\uparrow\rangle_2\right) (|\uparrow\rangle_3 |\downarrow\rangle_4 - |\downarrow\rangle_3 |\uparrow\rangle_4) (|\uparrow\rangle_5 |\downarrow\rangle_6 - |\downarrow\rangle_5 |\uparrow\rangle_6)

这个二聚化的、二重简并的基态与 J2=0J_2=0 时无能隙、平移不变的基态在性质上截然不同，这自然表明在某个中间的 J2/J1∈(0,1/2)J_2/J_1 \in (0,1/2) 处存在一个相变——数值上已知该相变发生在 J2/J1≈0.2411J_2/J_1 \approx 0.2411 附近。

在本教程的第一部分，我们将通过观察能谱——特别是不同对称子空间中的能隙——随耦合调节的变化，来定位临界点的位置。在第二部分，我们将重新考察临界链的 CFT 内容。可以解析地证明，临界点处的模型由与海森堡链相同的 CFT 描述，但导致对数消失的有限尺寸修正的边际算符权重为零，因此可以更精确地确定标度维数。

参数

参数含义取值
LATTICE带有最近邻和次近邻键的周期链nnn chain lattice
MODEL量子自旋模型spin
local_S每个格点的自旋量子数1/2
J最近邻耦合 J1J_11
J1次近邻耦合 J2J_2（注意：ALPS 的 nnn chain lattice 将该键参数命名为 J1，与上面哈密顿量中使用的 J2J_2 不同）第一部分在 [0,0.5][0,0.5] 范围内扫描，第二部分固定为 0.25
CONSERVED_QUANTUMNUMBERS、Sz_total分别解析单态（Sz=0S_z=0）和三重态（Sz=1S_z=1）子空间Sz, 0 和 1
NUMBER_EIGENVALUES向 sparsediag 请求的本征态数2（第一部分）或 5（第二部分）
L链长6, 8（第一部分）或 10, 12（第二部分）

格子

nnn chain lattice 在周期链的基础上扩展了第二组键，将每个格点与其次近邻相连：

        J          J          J
   o---------o---------o---------o---(periodic)
   0         1         2         3
    \___________________________/
        J1 (=J2)   next-nearest-neighbour bonds

这个内置格子是 ALPS 格子库中记载的其余格子之一，是前面教程中使用的普通 chain lattice 的扩展。

方法

与 ED-02 中一样，我们用 sparsediag 分别对角化 Sz=0S_z=0 和 Sz=1S_z=1 子空间，以得到单态和三重态能隙；这里最大的子空间（L=12L=12，Sz=0S_z=0，维数 924）对 Lanczos 算法而言微不足道。我们以足够密的间隔扫描 J2J_2 以观察到能隙交叉，从而定位临界点的有限尺寸估计值，然后在接近该交叉点的 J2/J1=0.25J_2/J_1=0.25 处再次对角化，请求更多的本征态以解析 CFT 能级塔。

那么首先，让我们绘制基态和第一激发态的能量，以及单态（Sz=0S_z = 0）和三重态（Sz=1S_z=1）子空间中的能隙。

使用命令行

parm5a 针对 L=6L=6 和 L=8L=8，在 Sz=0S_z=0 和 Sz=1S_z=1 两个子空间中都对 J2/J1∈[0,0.5]J_2/J_1\in[0,0.5] 进行了扫描：

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
NUMBER_EIGENVALUES=2
Sz_total=0
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }
Sz_total=1
{ L=6; J1=0.0 }
{ L=6; J1=0.1 }
{ L=6; J1=0.2 }
{ L=6; J1=0.3 }
{ L=6; J1=0.4 }
{ L=6; J1=0.5 }
{ L=8; J1=0.0 }
{ L=8; J1=0.1 }
{ L=8; J1=0.2 }
{ L=8; J1=0.3 }
{ L=8; J1=0.4 }
{ L=8; J1=0.5 }
parameter2xml parm5a
sparsediag --write-xml parm5a.in.xml

Python

为了自动化扫描和数据分析，我们使用脚本 tutorial5a.py。我们从常规的导入开始：

import pyalps
import pyalps.plot
from pyalps.dict_intersect import dict_intersect
import numpy as np
import matplotlib.pyplot as plt
import copy
import math

同样，我们使用 SzS_z 量子数，但现在将在不同的子空间 (Sz=0,1)(S_z=0,1) 中运行模拟。我们对系统尺寸 L=6,8L=6,8 进行计算，因为我们要寻找的效应在非常小的系统尺寸下就已经出现。

prefix = 'alps-nnn-heisenberg'
parms = []
for L in [6,8]:
    for Szt in [0,1]:
        for J1 in np.linspace(0,0.5,6):
            parms.append({
                'LATTICE'              : "nnn chain lattice",
                'MODEL'                : "spin",
                'local_S'              : 0.5,
                'J'                    : 1,
                'NUMBER_EIGENVALUES'   : 2,
                'CONSERVED_QUANTUMNUMBERS' : 'Sz',
                'Sz_total'             : Szt,
                'J1'                   : J1,
                'L'                    : L
            })

input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
# res = pyalps.runApplication('sparsediag', input_file, MPI=4)
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

这次的数据分析比之前的稍微复杂一些。特别是，我们将大量依赖分层数据集这一特性。要理解其中的物理，实际上只需查看基态和第一激发态就足够了——所以如果你觉得能隙的计算过于令人困惑，不必过于担心。

第一步，我们将给定 J1、L、Sz_total 组合下的所有能量合并并排序。首先，我们按参数——J1、L、Sz_total——分组。对 grouped 的每次循环中的元素因此将包含一组不同动量对应的数据集列表。我们将把它们合并，然后使用 dict_intersect 函数来确定合并后数据集的属性；该函数只是接收一个字典列表，并返回它们共有的部分。我们使用 numpy 的 argsort 函数获得将对 y 排序的索引列表；这样也可以相应地对 x 排序，尽管这可能并不是必需的。

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L', 'Sz_total'])
nd = []
for group in grouped:
    ally = []
    allx = []
    for q in group:
        ally += list(q.y)
        allx += list(q.x)
    r = pyalps.DataSet()
    sel = np.argsort(ally)
    r.y = np.array(ally)[sel]
    r.x = np.array(allx)[sel]
    r.props = dict_intersect([q.props for q in group])
    nd.append( r )
data = nd

接下来，我们需要把出现在 Sz=1S_z=1 子空间中的态从 Sz=0S_z=0 子空间中去除。我们按 J1、L 对它们分组，使得每一组都包含两个不同 Sz_total 子空间的能谱。然后我们使用函数 subtract_spectrum，它会从作为第一个参数传入的数据集中，移除同时也出现在第二个参数中的元素。作为可选参数，该函数还接受一个最大相对差。

grouped = pyalps.groupSets(pyalps.flatten(data), ['J1', 'L'])
nd = []
for group in grouped:
    if group[0].props['Sz_total'] == 0:
        s0 = group[0]
        s1 = group[1]
    else:
        s0 = group[1]
        s1 = group[0]
    s0 = pyalps.subtract_spectrum(s0, s1)
    nd.append(s0)
    nd.append(s1)
data = nd

现在，我们建立一个新的数据集列表（sector_E），其中只包含基态（‘gs’）或第一激发态（‘fe’）的能量。我们会把这个信息存储在属性 which 中。这将使我们随后可以使用 collectXY 函数，为每个 L 绘制基态和第一激发态能量随耦合变化的图像。

sector_E = []
grouped = pyalps.groupSets(pyalps.flatten(data), ['Sz_total', 'J1', 'L'])
for group in grouped:
    allE = []
    for q in group:
        allE += list(q.y)
    allE = np.sort(allE)
    
    d = pyalps.DataSet()
    d.props = dict_intersect([q.props for q in group])
    d.x = np.array([0])
    d.y = np.array([allE[0]])
    d.props['which'] = 'gs'
    sector_E.append(d)
    
    d2 = copy.deepcopy(d)
    d2.y = np.array([allE[1]])
    d2.props['which'] = 'fe'
    sector_E.append(d2)

sector_energies = pyalps.collectXY(sector_E, 'J1', 'Energy', ['Sz_total', 'which', 'L'])
plt.figure()
pyalps.plot.plot(sector_energies)
plt.xlabel('$J_1/J$')
plt.ylabel('$E_0$')
plt.legend(prop={'size':8})

最后一步，我们计算单态和三重态能隙。它们分别定义为体系最低态的能量与 a) 单态（Sz=0S_z=0）子空间中第一激发态、b) 三重态（Sz=1S_z=1）子空间中最低态之间的能量差。

grouped = pyalps.groupSets( pyalps.groupSets(pyalps.flatten(data), ['J1', 'L']), ['Sz_total'])

gaps = []
for J1g in grouped:
    totalmin = 1000
    for q in pyalps.flatten(J1g):
        totalmin = min(totalmin, np.min(q.y))
    
    for Szg in J1g:
        allE = []
        for q in Szg:
            allE += list(q.y)
        allE = np.sort(allE)
        d = pyalps.DataSet()
        d.props = pyalps.dict_intersect([q.props for q in Szg])
        d.props['observable'] = 'gap'
        print totalmin,d.props['Sz_total']
        if d.props['Sz_total'] == 0:
            d.y = np.array([allE[1]-totalmin])
        else:
            d.y = np.array([allE[0]-totalmin])
        d.x = np.array([0])
        d.props['line'] = '.-'
        gaps.append(d)

gaps = pyalps.collectXY(gaps, 'J1', 'gap', ['Sz_total', 'L'])

plt.figure()
pyalps.plot.plot(gaps)
plt.xlabel('$J_1/J$')
plt.ylabel('$\Delta$')
plt.legend(prop={'size':8})

plt.show()

输出数据

独立地对各子空间进行对角化，得到以下作为 J2/J1J_2/J_1（如上所述，在 ALPS 参数中标记为 J1）函数的单态和三重态能隙：

LLJ2/J1J_2/J_1单态能隙 /J1/J_1三重态能隙 /J1/J_1
60.01.30280.6847
60.11.03000.6688
60.20.76210.6455
60.250.63020.6302
60.30.50000.6119
60.50.00000.5000
80.00.95150.5227
80.10.75510.5039
80.20.55830.4800
80.250.46010.4663
80.30.36240.4516
80.50.00000.4045

对两种系统尺寸，单态能隙（起初较大，因为 J2=0J_2=0 时体系无能隙，最低激发只有在 L→∞L\to\infty 时才趋于零）在 J2/J1=0.2J_2/J_1=0.2 到 0.30.3 之间降到三重态能隙以下，交叉点几乎恰好在 J2/J1≈0.25J_2/J_1\approx0.25——即使在这样小的尺寸下也已经接近公认的热力学极限值 J2/J1≈0.2411J_2/J_1\approx0.2411。这个能级交叉正是进入二聚化相这一相变的有限尺寸标志；还请注意，单态能隙恰好在 Majumdar-Ghosh 点 J2/J1=0.5J_2/J_1=0.5 处消失，这与该点二重简并的基态是一致的。

具有次近邻耦合的海森堡链：CFT 归属

通过调节到阻挫 J1-J2 链中的一个临界点，可以显著减小有限尺寸修正。尽管耦合方式不同，可以证明该模型具有相同的连续临界场论，因此我们可以在这一极限下提取标度维数。关于这一点的详细讨论，请参阅参考文献 I. Affleck, D. Gepner, H.J. Schulz and T. Ziman, J. Phys. A: Math. Gen. 22, 511 (1989)。

与你上面得到的能谱相比较：你会发现与预期标度维数的对应关系要清晰得多，并且随着系统尺寸增大，收敛速度也快得多。

使用命令行

parm5b 将 J2/J1=0.25J_2/J_1=0.25 固定在上面找到的交叉点附近，并针对两种系统尺寸，在 Sz=0S_z=0 子空间中请求 5 个本征态：

MODEL="spin"
LATTICE="nnn chain lattice"
CONSERVED_QUANTUMNUMBERS="Sz"
local_S=1/2
J=1
J1=0.25
NUMBER_EIGENVALUES=5
Sz_total=0
{ L=10 }
{ L=12 }
parameter2xml parm5b
sparsediag --write-xml parm5b.in.xml

Python

脚本 tutorial5b.py 中使用的新参数为：

parms_ = {
    'LATTICE'              : "nnn chain lattice",
    'MODEL'                : "spin",
    'local_S'              : 0.5,
    'J'                    : 1,
    'J1'                   : 0.25,
    'NUMBER_EIGENVALUES'   : 5,
    'CONSERVED_QUANTUMNUMBERS' : 'Sz',
    'Sz_total' : 0
}
prefix = 'nnn-heisenberg'
parms = []
for L in [10,12]:
    parms_.update({'L':L})
    parms.append(copy.deepcopy(parms_))

脚本的其余部分与 ED-04 中海森堡链的 CFT 分析完全一致：写入并运行输入文件，加载本征态测量结果，按 Sz=0S_z=0 子空间中最低两个态之间的能隙重新标度，并叠加预期的标度维数：

input_file = pyalps.writeInputFiles(prefix,parms)
res = pyalps.runApplication('sparsediag', input_file)
data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix=prefix))

E0 = {}
E1 = {}
for Lsets in data:
    L = pyalps.flatten(Lsets)[0].props['L']
    allE = []
    for q in pyalps.flatten(Lsets):
        allE += list(q.y)
    allE = np.sort(allE)
    E0[L] = allE[0]
    E1[L] = allE[1]

for q in pyalps.flatten(data):
    L = q.props['L']
    q.y = (q.y-E0[L])/(E1[L]-E0[L]) * (1./2.)
spectrum = pyalps.collectXY(data, 'TOTAL_MOMENTUM', 'Energy', foreach=['L'])

for SD in [0.5, 1, 1.5, 2]:
    d = pyalps.DataSet()
    d.x = np.array([0,4])
    d.y = SD+0*d.x
    spectrum += [d]

pyalps.plot.plot(spectrum)
plt.legend(prop={'size':8})
plt.xlabel("$k$")
plt.ylabel("$E_0$")
plt.xlim(-0.02, math.pi+0.02)
plt.show()

输出数据

以与上面相同的方式，对 J2/J1=0.25J_2/J_1=0.25 处独立对角化得到的 Sz=0S_z=0 子空间能谱进行重新标度，得到：

LL能级 1能级 2能级 3能级 4、5
1000.500（定义值）0.9471.401
1200.500（定义值）0.9601.432

将其与 ED-04 中普通海森堡链（J2=0J_2=0）在相同尺寸下的数值相比较：对于应收敛到 1 的能级，L=10L=10 时为 0.880，L=12L=12 时为 0.857。而在 J2/J1=0.25J_2/J_1=0.25 处，同一能级已经达到 0.947 和 0.960——既更接近目标值，又随 LL 增大朝正确的方向移动，这正是一旦调节掉了 ED-04 中导致对数修正的边际算符之后所预期的行为。

总结

比较阻挫 J1J_1–J2J_2 链的单态和三重态能隙，通过有限尺寸能级交叉，仅用 L=6,8L=6,8 两条链就能将二聚化相变定位到与公认值 J2/J1≈0.2411J_2/J_1\approx0.2411 相差仅百分之几的精度；在该点附近进行对角化还表明，ED-04 中辨认出的同一 c=1c=1 CFT 算符内容在这里出现得清晰得多，因为调节 J2J_2 消除了导致无阻挫链中缓慢对数有限尺寸修正的边际算符。

问题

  • 根据单态和三重态能隙表，通过在相邻数据点之间做线性插值，分别对 L=6L=6 和 L=8L=8 估计交叉点 J2/J1J_2/J_1。随着 LL 增大，这个估计值是否向 0.24110.2411 靠近？
  • 为什么恰好是单态能隙、而不是三重态能隙，在 Majumdar-Ghosh 点 J2/J1=0.5J_2/J_1=0.5 处消失？
  • 将 J2/J1=0.25J_2/J_1=0.25 处 L=10,12L=10,12 重新标度后的第三能级数值，与 ED-04 中 J2=0J_2=0 时的对应数值相比较。哪一个更接近 CFT 预言的 1？哪一个随 LL 收敛得更快？
ED-04 CriticalityED-06 Full Diagonalization