ED-03 Spectra

一维量子系统的能谱

在本教程中，我们将计算量子海森堡模型在多种一维格子上的能谱。主要工作将由 sparsediag 应用程序完成，它实现了 Lanczos 算法，一种迭代本征求解器。与前面几篇教程中的基态和能隙计算不同，这里我们在每个动量子空间中计算若干个最低本征值，从而构建出完整的、按动量分辨的激发谱 $E(k)$ 。仅仅通过观察低能谱随格子几何结构和耦合变化的形状，这是判断一个模型具有何种激发——单支色散磁振子、二自旋子连续谱、带能隙的三重子束缚态——的基本诊断手段。

下面的三个例子（链、梯子、孤立二聚体）都是同一个双腿梯子海森堡哈密顿量的特殊情形，

H = J_0 \sum_{\text{legs}} \mathbf{S}_i \cdot \mathbf{S}_j \; + \; J_1 \sum_{\text{rungs}} \mathbf{S}_i \cdot \mathbf{S}_j ,

其中 $J_0$ 是沿两条腿方向的耦合， $J_1$ 是横跨梯级的耦合；取 $J_0=J_1$ 得到各向同性梯子， $J_1=0$ 得到两条脱耦的链，而 $J_0=0$ 则得到 $L$ 个孤立的两格点二聚体。这一族模型及其极限之间的交叉行为，参见例如 E. Dagotto and T.M. Rice, Science 271, 618 (1996)。

参数

参数	含义	取值
`LATTICE`	内置几何结构	`chain lattice`（链）或 `ladder`（梯子/二聚体）
`MODEL`	量子自旋模型	`spin`
`local_S`	每个格点的自旋量子数	1/2
`J`	最近邻耦合（chain lattice）	1
`J0`	腿方向耦合（ladder lattice）	1（梯子）或 0（孤立二聚体）
`J1`	梯级耦合（ladder lattice）	1
`L`	线性尺寸（梯子的梯级数）	10–16（链），6–10（梯子/二聚体）
`CONSERVED_QUANTUMNUMBERS`、`Sz_total`	限制在 $S_z=0$ 子空间	`Sz`, 0

方法

我们使用 sparsediag，因为要构建 $E(k)$ 不仅需要基态，还需要每个动量子空间中若干个最低本征值；Lanczos 算法可以直接计算稀疏哈密顿量中任意指定数目的极端本征值，而无需对角化我们不关心的子空间。这里最大的子空间（链， $L=16$ ， $S_z=0$ ）维数为 $\binom{16}{8}=12870$ ，对于迭代式稀疏对角化来说完全可以胜任，而手工枚举则远远无法做到。

海森堡链

chain lattice（内置于 ALPS 格子库）是一个具有均匀耦合 $J$ 的、由 $L$ 个格点组成的周期性环：

    J     J     J           J
o-------o-------o--- ... ---o
0       1       2          L-1
|_________________________________|
                J   (bond L-1 -- 0, periodic)

从命令行准备并运行模拟

首先，我们来看一条具有海森堡耦合的 S=1/2 自旋链。参数文件 parm_chain 为具有 {L=10,…16} 个自旋的链，设置了 S_z=0 子空间的 ED 模拟。

LATTICE = "chain lattice", 
MODEL = "spin",
local_S = 0.5,
J = 1,
CONSERVED_QUANTUMNUMBERS = "Sz"
Sz_total = 0
{ L = 10; }
{ L = 12; }
{ L = 14; }
{ L = 16; }

使用下面这一系列命令即可运行对角化，然后用浏览器查看输出文件 parm_chain.out.xml。

parameter2xml parm_chain
sparsediag --write-xml parm_chain.in.xml

使用 Python 准备并运行模拟

要在 Python 中设置并运行模拟，我们使用脚本 chain.py。你可以在终端中用 python chain.py 来运行它。观察脚本的各个部分，我们可以看到在导入所需模块之后，输入文件是如何被准备为一个 Python 字典列表的。

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

parms=[]
for l in [10, 12, 14, 16]:
    parms.append({ 
        'LATTICE'                   : "chain lattice", 
        'MODEL'                     : "spin",
        'local_S'                   : 0.5,
        'J'                         : 1,
        'L'                         : l,
        'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
        'Sz_total'                  : 0
    })

接下来，输入参数被写入 XML 作业文件，然后运行 sparsediag 模拟。

input_file = pyalps.writeInputFiles('parm_chain',parms)
res = pyalps.runApplication('sparsediag',input_file)

为了绘制能谱，我们接着加载模拟产生的 HDF5 文件：

data = pyalps.loadSpectra(pyalps.getResultFiles(prefix='parm_chain'))

并将每个系统尺寸 L 的所有动量子空间的能量收集到一个 DataSet 中。为了得到一张美观的图，我们还从所有本征值中减去基态能量，并为每条能谱指定一个标签和线型。

spectra = {}
for sim in data:
    l = int(sim[0].props['L'])
    all_energies = []
    spectrum = pyalps.DataSet()
    for sec in sim:
        all_energies += list(sec.y)
        spectrum.x = np.concatenate((spectrum.x,np.array([sec.props['TOTAL_MOMENTUM'] for i in range(len(sec.y))])))
        spectrum.y = np.concatenate((spectrum.y,sec.y))
        spectrum.y -= np.min(all_energies)
    spectrum.props['line'] = 'scatter'
    spectrum.props['label'] = 'L='+str(l)
    spectra[l] = spectrum

现在可以将不同系统尺寸的能谱绘制到同一张图中：

plt.figure()
pyalps.plot.plot(spectra.values())
plt.legend()
plt.title('Antiferromagnetic Heisenberg chain (S=1/2)')
plt.ylabel('Energy')
plt.xlabel('Momentum')
plt.xlim(0,2*3.1416)
plt.ylim(0,2)
plt.show()

海森堡链绘制出的能谱如下图所示：

双腿海森堡梯子

只需对上面使用的输入参数做少量修改，我们就可以计算双腿海森堡自旋梯子的能谱。ladder 格子将 $2L$ 个格点排列成两行，每行 $L$ 个，通过腿方向的键 $J_0$ 和梯级键 $J_1$ 相连：

o---J0---o---J0---o---...---o     leg 0
|        |        |         |
J1       J1       J1        J1
|        |        |         |
o---J0---o---J0---o---...---o     leg 1

新的参数文本文件 parm_ladder 内容如下：

LATTICE = "ladder"
MODEL = "spin"
local_S = 0.5
J0 = 1
J1 = 1
CONSERVED_QUANTUMNUMBERS = "Sz"
Sz_total = 0
{ L = 6; }
{ L = 8; }
{ L = 10; }

我们只是将 “chain lattice” 替换为 “ladder”，并分别为腿方向和梯级方向定义了两个独立的耦合常数 J0、J1。除此之外，我们还缩小了线性系统尺寸 L，因为梯子共有 2L 个自旋。按照与之前完全相同的方式运行：

parameter2xml parm_ladder
sparsediag --write-xml parm_ladder.in.xml

Python 代码也需要做同样的修改，可以从这里下载：ladder.py

不同格子尺寸下海森堡梯子的能谱如下图所示：

孤立二聚体

如果我们将梯子腿方向的耦合设为 J0 = 0，就得到了 L 个孤立二聚体的能谱——每个梯级脱耦为一个独立的两格点单态/三重态问题，因此"能谱"在任意动量值下都退化为恰好两条能级（单态 $E=-3J_1/4$ ，以及平坦的三重态能带 $E=+J_1/4$ ）。这是通过参数文件 parm_dimers 实现的：

LATTICE = "ladder"
MODEL = "spin"
local_S = 0.5
J0 = 0
J1 = 1
CONSERVED_QUANTUMNUMBERS = "Sz"
Sz_total = 0
{ L = 6; }
{ L = 8; }
{ L = 10; }

以及 Python 脚本 dimers.py。运行方式相同：

parameter2xml parm_dimers
sparsediag --write-xml parm_dimers.in.xml

孤立二聚体的能谱如下图所示：

总结

打开或关闭单个耦合常数（ $J_0$ ）会使激发谱连续地从无色散的平带（孤立二聚体），经过具有束缚态/连续谱结构的两支能谱（梯子），演变为单链的无能隙二自旋子连续谱——这说明仅凭能谱的形状就可以揭示系统低能激发的本质。

问题

观察将不同系统尺寸的能谱组合在一起会如何形成漂亮的能带
在海森堡梯子的能谱中：辨认出连续谱和束缚态
链的能谱和梯子的能谱之间的主要区别是什么？
解释孤立二聚体的能谱
改变梯子中的耦合常数，观察能谱如何在之前讨论过的极限之间变化
附加问题：仔细观察不同系统尺寸下链的能谱：L/2 为偶数和为奇数的情形之间似乎存在差异。你能解释这一点吗？在 L 趋于无穷的热力学极限下会发生什么？

ED-02 Gaps ED-04 Criticality