海森堡模型的谱，到底有多靠谱
作者：孟子杨（中科院物理所副研究员）
靠谱，还是不靠谱，这是一个值得思考的问题。
这个问题不仅关乎哈姆雷特对于人生悖论的追问，其实在凝聚态物理学量子多体问题的动力学性质计算中，也是有一个让多少仁人志士扼腕叹息的，带有希腊悲剧意味的、形而上的命题。
所谓动力学性质，主要是指谱学行为，如关联电子系统中的准粒子 (quasiparticle) 能谱，如量子磁学系统中的自旋波磁振子 (magnon) 能谱。这类能量、动量依赖的谱函数，可以告诉人们量子多体系统的诸多本质信息，而且与现代凝聚态物理学的实验手段直接相关。比如角分辨光电子谱技术测量的就是固体材料的电子结构（即准粒子能谱），而中子散射技术测量的，就是量子磁学材料中的自旋波磁振子能谱。从准粒子谱中，可以看到费米面的形状，可以看到费米面上准粒子权重的分布，比如高温超导体中的费米弧 (Fermi arc)，还有拓扑绝缘体、半金属中的狄拉克锥 (Dirac cone)。从自旋波磁振子谱中可以看到系统的自旋排布，更可以看到自旋波磁振子的色散关系，从而确定自旋交换相互作用 J 的大小。再比如目前人们正在寻找的量子自旋液体材料中的连续谱和连续谱背后所蕴涵的分数化自旋子 (spinon) 元激发，这些谱学行为是寻找与确定拓扑序这种超越了朗道对称性破缺框架的物质形态的直接证据。
可见，谱是如此重要，谱学的结果一定要靠谱。但是，上面所说的，希腊悲剧意味的命题的意思就是，对于凝聚态量子多体问题的动力学性质的严格计算，从理论上讲相当困难。量子多体问题的动力学性质牵扯到多体系统的时间演化问题，对于具有指数多自由度的强关联系统，其静力学性质，如能量、序参量还有种种关联函数已经很难准确计算，而动力学性质的计算，即准确计算相互作用的指数多的自由度随着时间的演化，几近不可能完成的任务。这好比希腊悲剧中的英雄人物，任你有高贵的性格、一身的本事，命运总是无法抗拒，一如俄狄浦斯，一如安提戈涅。
在个别问题上，我们有严格可解的模型，但主要集中在一维量子系统。对于二维或者更高维的量子多体系统，解析的方法只能提供微扰论意义下的近似，比如量子磁学系统中的自旋波理论和它的高阶修正，真正严格的计算，还是需要发展数值计算方法。到目前为止，许多看似基本的问题，比如反铁磁海森堡模型的自旋激发谱，人们还没有得到全局性的认识。目前领域的前沿，就是逐步做到准确计算量子多体系统基本模型的动力学性质，然后和实验结果进行比对，验证计算结果是否靠谱。
中科院物理所理论室博士生秦彦齐，副研究员孟子杨，与北京计算科学研究中心邵慧博士，Stefano Chesi 副研究员，法国图卢兹大学 Sylvain Capponi 教授，美国波士顿大学 Anders Sandvik 教授，组成的研究团队，在计算量子多体问题动力学性质方面，进行了有益的尝试 [1]。他们将一套结合量子蒙特卡洛 (quantum Monte Carlo, QMC) 和随机解析延拓 (Stochastic analytic continuation, SAC) 的计算方法应用到二维方格子反铁磁海森堡模型的动力学行为研究之中，得到了海森堡模型 – 量子磁学系统的最基本模型 – 的完整能谱。如下文中图1、图2 所示，他们的结果既可以与最新的 Cu(DCOO)2.4D2O 材料中子散射的实验结果 [2] 进行对比，更揭示了海森堡模型的高能自旋能谱中蕴涵着分数化自旋子激发的痕迹，从而建立了海森堡模型与去禁闭量子临界行为的深层联系。这个发现说明，去禁闭量子临界行为和其伴随的分数化自旋子激发，其实在磁学系统的基本模型中也存在，它们是超越朗道对称性破缺和序参量的语言，是对于一大类量子相和量子相变问题更加完整的描述。
图 1. （上）二维方格子反铁磁海森堡模型与尼耳态。（中）Cu(DCOO)2.4D2O 材料的中子散射实验结果 [2]，这个系统很好地实现了二维反铁磁海森堡模型，其自旋激发谱与（下）中的量子蒙特卡洛+随机解析延拓的结果吻合，而且，在（中）和（下）里，系统的自旋激发谱在动量 (π,0) 附近都有连续谱出现。
这里，我们首先简单讨论一下量子蒙特卡洛+随机解析延拓 (QMC+SAC) 的这套方法的过程 [3-5]。
在量子蒙特卡洛计算中，要获得量子多体问题的实时关联函数并不容易，但虚时关联是可以严格计算的，比如对于物理可观测量 O，其虚时的关联函数为
其中
其中 S(ω) 即为谱函数，知道它便可以得到虚时关联函数 G(τ) 。但是这个问题的复杂性在于，我们掌握的信息，其实是通过量子蒙特卡洛计算所得的、带有统计误差的 G(τ) , 而为了得到可以和实验直接对比的能谱 S(ω) ，我们需要反解 Eq. (2)。这是数学物理中经常出现的、棘手的反问题 （这里是拉普拉斯变换的反过程）：知道 S(ω) 求 G(τ) 很容易，因为这是由因及果；而知道 G(τ) 求 S(ω) 很不容易，因为要由果及因，何况这个果，G(τ)，还是带着统计误差的。反解的存在性、唯一性和稳定性，都没有保证。这正是量子多体问题解析延拓困难的地方。为了解决这个问题，人们尝试了很多的办法，随机解析延拓是其中比较成功的一种。
随机解析延拓，顾名思义，就是把反解谱函数的过程转化成一个随机过程，进而运用蒙特卡洛的方法，在统计的意义上等到谱函数的分布和其最可几的期望值。
抛开技术上的细节，整个过程是这样的。我们假设谱函数在频率空间有如下形式，
其中参数 {ai , ωi} 刻画了谱函数的所有信息（Nω 多个 δ -函数的位置和高度），这样的假设涵盖了谱函数尽可能多的自由度。然后由 Eq. (2) 可以得到在这组参数下的虚时关联函数 G(τ) , 它与量子蒙特卡洛得到的准确的虚时关联之间的关系，可以用 一个 “goodness of fit” χ2 来描写，
其中 {Gi} 是量子蒙特卡洛得到的系综平均后的虚时关联函数值（即平均值），而
是量子蒙特卡洛数据之间的协方差矩阵，它包含了蒙卡数据之间的关联和统计误差等信息。N 是虚时点的个数，NB 是蒙特卡洛测量的 bin 数。
在如许构造之下，解析延拓的问题转换成了满足 Eq. (2) 和 Eq. (4) 的优化问题。如上所述，这是一个反问题，它不一定有简单的、确定性的最优解，对于这个问题的正确理解，是我们需要找到 { S(ω) } 的一个分布，这个分布的统计平均（系综平均）就扮演了那个最优解，而分布的误差则告诉我们可以在多大的置信区间中讨论这个最优解。行文至此，一切都变得明朗了，寻找一个正确的分布，按照其进行构型更新、抽样、系综平均，正是蒙特卡洛方法屡试不爽的法门。
所以，随机解析延拓的过程，本身就是一个蒙特卡洛的过程。如上所述，我们可以用 Eq. (4) 中的 χ2 来写下一个 Boltzmann 分布，
其中 P(S) 就是某个谱函数 S(ω) 构型的权重，权重有其对应的能量 χ2 和温度 θ , 这样解析延拓的问题就有了统计物理问题的味道，接下来的过程就和一般的、运用蒙特卡洛方法解决统计物理问题的过程没有差别了。具体而言，我们以 Eq. (3) 中的一套 {ai , ωi} 为一个构型，然后开始构型的更新，比如，
ωi → ωi + d
其中 d 是一个频率的随机改变量，然后按照 Eq. (5) 计算构型更新的 transition probability：如果构型的权重升高，更新一定接收；如果构型的权重降低，则以一定的概率（就是 exp(-χ2 / 2θ)）接收。如许的 Metropolis 抽样过程一直进行下去，我们就可以得到 S(ω) 的分布，其系综平均的结果，就是我们要寻找的谱函数。更多的细节，以及对于不同问题所采取的不同更新方法，读者请参见文献 [1,3-5]。
按照这样的方式，QMC+SAC, 一大类量子多体问题的动力学性质计算就可以进行下去了。步骤就是首先进行量子蒙特卡洛计算本身，得到高质量的系统的虚时关联函数，比如准粒子虚时关联函数（就是动力学格林函数）还有自旋关联函数，这一步的计算量比较大，往往需要在高性能计算平台上进行大规模并行计算。然后进行随机解析延拓，这一步的蒙特卡洛过程，相对简单，在本地就可以完成。这样得到的，就是靠谱的能谱。
在参考文献 [1] 中，我们研究了二维反铁磁海森堡模型的自旋能谱，主要结果总结在图1、图2之中。可以看到，系统的基态具有反铁磁长程序，自旋波磁振子 (magnon) 在动量点 (π, π) 没有能隙，这是所谓的 Goldstone mode，对应着海森堡模型自旋连续对称性的自发破缺，而且，磁振子在 (π, π) 的能谱，基本上就是一个 delta-函数，90% 以上的谱权重都在这个 delta-函数里。这个结果，与 Cu(DCOO)2.4D2O 材料的中子散射实验结果以及线性自旋波的解析计算，都是吻合的。
然而，更加有意思的是，实验和 QMC+SAC 的计算结果都发现，在动量点 (π,0) 附近，磁振子的位置，即第一激发态的能隙，明显小于线性自旋波的计算结果，如图2 （上）所示。而且，在自旋波磁子的 delta-函数之上，有相当宽的连续谱，这里和 (π, π) 点的情况相反，连续谱占据了近 60%的谱权重，而 delta-函数只有40%。在 (π,0) 的谱线峰之上，系统还具有足够展宽的谱权重。这些谱权重，是自旋波理论不能解释的。我们不禁要问，这 60% 的连续谱，到底是什么来历？
图 2. （上）随机解析延拓（QMC+SAC）计算出的磁振子色散关系（系综平均值）与Cu(DCOO)2.4D2O 材料中子散射实验结果的对比，吻合近乎完美。同时线性自旋波的计算结果 (Linear SWT)，对于反铁磁长程序的描述,在(π, π) 点附近没有问题，但是在 (π, 0) 点附近失效。这是因为自旋波的计算中，无法考虑分数化的自旋子激发。（下）随机解析延拓（QMC+SAC）计算出的能谱与Cu(DCOO)2.4D2O 材料中子散射实验结果的在(π,0) 和 (π/2,π/2) 点的细节对比。可以看到，在 (π, 0) 点，除了 delta-函数之外，高能连续谱占有相当大的谱权重 ；而在 (π/2, π/2) 点，主要的权重都集中在 delta-函数里。
为了更好的理解 (π,0) 附近的连续谱，我们仔细研究了系统在这里的动力学性质，并将海森堡模型和具有去禁闭量子临界点的 J-Q 的模型的动力学计算结果，做了认真地对比。关于去禁闭量子临界点与J-Q 模型，参见参考文献 [6, 7]，这里只介绍对比的结果。如图3（左）所示， 在J-Q 模型中，随着 Q/J 的增大，系统从反铁磁长程序的基态逐渐逼近去禁闭量子临界点，这里可以很清楚地看到，在 (π, 0) 点上的磁振子谱权重，随着 Q/J 的增大，急速地消失了，这表明即使在反铁磁长程序中，(π, 0) 点附近的 S=1磁振子，也在不断地分数化成 S=1/2 的自旋子。这个现象，之前没有被发现过，而我们 QMC+SAC 的计算刚好揭示了这个奇异的过程。这一发现帮助我们理解了海森堡模型在 (π, 0) 点附近超越自旋波理论的奇异行为：作为去禁闭量子临界点的先兆，整数自旋 (S=1) 的磁振子在此处分数化成了近乎去禁闭状态的半整数自旋 (S=1/2) 的自旋子。自旋子之间具有相互作用，自旋子和磁振子之间也有相互作用，这些相互作用改变了能谱在 (π, 0) 附近的形状，使其偏离线性自旋波，使其产生连续谱。从反铁磁长程序到去禁闭量子临界点这整个过程中，系统动力学行为的详细演化，我们正在做着进一步的研究。显而易见的是，分数化、去禁闭、超越朗道对称性破缺框架的凝聚态物理学相变理论新范式，其实就蕴藏在我们熟悉的基本模型之中，只是之前我们没有如 QMC+SAC 这样得力的工具，无法得到靠谱的能谱。
图 3. （左）QMC+SAC计算出的自旋波能谱在 (π, 0) 点和 (π/2, π/2) 点的行为。(a) 为磁振子谱峰的位置，(b) 为磁振子的谱权重。可见随着 Q/J 的增大（系统越来越靠近去禁闭量子临界点），(π, 0) 点的磁振子谱权重急速消失，系统的低能激发从整数自旋 (S=1) 的磁振子变成分数自旋 (S=1/2) 的自旋子。(π/2, π/2) 点则没有这样的变化。（右）QMC+SAC 计算出的J-Q 模型自旋能谱在布里渊区高对称线上的变化，可以看到随着 q=Q/(J+Q) 向着去禁闭量子临界点逼近 (qc~0.6)，(π, 0) 点连续谱越来越明显，显示出系统中磁激发的分数化行为、去禁闭行为越来越明显。
如开篇所述，量子多体问题的动力学信息的严格计算，是十分困难的问题，但是尴尬的是，动力学性质恰恰可以被现代凝聚态物理学实验手段直接观测。所以，理论和实验之间的差距变成了阻碍领域继续向前发展的，具有形而上意味的鸿沟。同时，凝聚态物理学中的新现象、新问题，如相互作用的拓扑物质形态、超越朗道对称性破缺框架的物质分类理论、去禁闭量子临界点、量子自旋液体中的涌现规范场、巡游电子量子相变和伴随的非费米液体行为，......，正在不断的出现，与之伴随的凝聚态物理学新范式正呼之欲出。对于这些现象的动力学性质的计算，即可以解释不断积累的实验现象，又可以推动理论的进一步发展，其重要性自不待言。再困难的哲学性问题，我们都要想办法克服。
这里讲述的 QMC + SAC 方法，能够在一些强关联模型中得到靠谱的能谱，应用前景十分广阔，比如在量子自旋冰模型中，人们讨论了很多年的规范场光子的激发谱，已经被 QMC+SAC 看到 [8]。希望，在不久的将来，靠谱，还是不靠谱，这个困扰着一代代诗人、哲人的问题，不再困扰我们。
[1] Nearly deconfined spinon excitations in the square-lattice spin-1/2 Heisenberg antiferromagnet,
Hui Shao, Yan Qi Qin, Sylvain Capponi, Stefano Chesi, Zi Yang Meng, Anders W. Sandvik,
Phys. Rev. X 7, 041072 (2017)
[2] Fractional excitations in the square-lattice quantum antiferromagnet,
B. Dalla, M. Mourigal, N. B. Christensen, G. J. Nilsen, P. Tregenna-Piggott, T. G. Perring, M. Enderle, D. F. McMorrow, D. A. Ivanov, H. M. Rønnow,
Nature Physics 11, 62–68 (2015)
[3] Stochastic method for analytic continuation of quantum Monte Carlo data,
Anders W. Sandvik,
Phys. Rev. B 57, 10287 (1998)
[4] Identifying the maximum entropy method as a special limit of stochastic analytic continuation
Kevin S. D. Beach,
arXiv:cond-mat/0403055
[5] Stochastic analytic continuation with various prominent spectral features,
Hui Shao and Anders W. Sandvik,
forthcoming
[6] Duality between the deconfined quantum-critical point and the bosonic topological transition,
Yan Qi Qin, Yuan-Yao He, Yi-Zhuang You, Zhong-Yi Lu, Arnab Sen, Anders W. Sandvik, Cenke Xu, Zi Yang Meng,
Phys. Rev. X 7, 031052 (2017)
[7] 西斯廷教堂中的对偶变换
http://mp.weixin.qq.com/s/XMZUlOI-CnXxkzZRyp48iQ
[8] Dynamics of topological excitations in a model quantum spin ice,
Chun-Jiong Huang, Youjin Deng, Yuan Wan, Zi Yang Meng,
arXiv:1707.00099
致谢
笔者感谢合作者，北京计算科学研究中心邵慧博士，对于这篇文字提出中肯的修改意见；感谢国家超级计算天津中心孟祥飞博士、赵洋工程师、菅晓东工程师等人对我们大规模蒙特卡洛计算所提供的资源、技术方面的有力支持。更想用这篇小文欢迎合作者，强关联电子系统数值计算领域的著名学者，Anders W. Sandvik 教授加入物理所，成为笔者的同事。可以想见，在量子多体问题大规模数值计算方面，更多靠谱的工作，将会从物理所不断涌现。
编辑：山寺小沙弥
近期热门文章Top10
↓ 点击标题即可查看 ↓
4. 换头手术真的成功了吗？
9. 这是一件占据你三分之一人生的大事，然而大多数人却对其背后的故事一无所知
10. 作死小能手再开脑洞：如果一个硬币大小的黑洞在你面前，会发生什么现象？