利用能级交叉和二聚化算符表征海森堡反铁磁系统的量子相变

Level Crossing and Dimerization Operator for Characterization of Quantum Phase Transition in the Heisenberg Antiferromagnet

摘要

用逆迭代计算阻挫自旋\mathrm{1}/\mathrm{2}海森堡反铁磁J₁-J₂方格子的基态和低激发态及相应的本征能量。在总自旋量子数的子空间中利用精确对角化和正交化波算符方法，得到哈密顿量\mathit{H}、总自旋量子数S、总自旋z分量量子数S_{\mathrm{T}}^z的共同本征态。利用能级交叉和二聚化算符，分析了中间相量子相变点，讨论了自旋液态存在的区域。计算结果表明，按照无量纲参数g=J_{\mathrm{2}}/J_{\mathrm{1}}，无序相分布在g=\mathrm{0.41}~\mathrm{0.71}，其两侧g=\mathrm{0.00}~\mathrm{0.40}和g=\mathrm{0.72}~\mathrm{0.99}分别为\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}相和条纹相。低激发自旋单态S=\mathrm{0}和自旋三重态S=\mathrm{1}的能级交叉将中间相分隔成g=\mathrm{0.41}~\mathrm{0.49}和g=\mathrm{0.50}~\mathrm{0.71}两个区域。在g=\mathrm{0.50}~\mathrm{0.52}的区域显示有低激发态二度简并的自旋单态存在，通过二聚化算符分析其自旋液态特征。

Abstract

We present numerical inverse-iterative result for the ground state and the low-lying excitations and corresponding eigen-energies for the frustrated Spin-\mathrm{1}/\mathrm{2}J₁-J₂ Square Heisenberg Antiferromagnet. Using exact diagonalization and orthogonal wave-operator in the subspace of total spin, we obtain the common eigen-states among Hamiltonian \mathit{H}, total spin number S, and its z-component S_{\mathrm{T}}^z. Using level crossing and dimerization operator, we analyze the quantum phase transition point in the intermediate region, and discuss spin liquid regime as well. As revealed by the result, in the metric of dimensionless parameter g=J_{\mathrm{2}}/J_{\mathrm{1}}, the nonmagnetic phase spreads in g=\mathrm{0.41}~\mathrm{0.71}, of which on the left and right sides are Neel and stripe phases for g=\mathrm{0.00}~\mathrm{0.40} and g=\mathrm{0.72}~\mathrm{0.99}, respectively. The low-lying level crossing between the spin singlet S=\mathrm{0} and triplet S=\mathrm{1} separates the intermediate region as two regimes g=\mathrm{0.41}~\mathrm{0.49} and g=\mathrm{0.50}~\mathrm{0.71}. Doubly degenerated spin singlets present in the low-lying excitations for g=\mathrm{0.50}~\mathrm{0.52}. Further, we use dimerization operator to analyze the characteristics of spin liquid in this regime.

Graphical abstract

0　 引 言

学者们对二维方格子阻挫自旋1/2海森堡模型的研究兴趣源于高温铜基超导体的提出^［1］。阻挫的反铁磁耦合对可以形成自旋液态（spin liquid，SL）^［2~5］，比如共振价键态（resonating valence bonds，RVB）^{［4，6~8］}可以通过掺杂的方式转变为超导态。海森堡模型作为一个二维系统原型，无论是在理论研究还是数值计算领域都具有很重要的意义。特别是顺磁相中精细的量子相变，它的本质仍在广泛讨论中，引起了许多研究者的参与兴趣^{［7，9~17］}。

该模型中间相的形成有多种猜想，比如可能是由于自旋单态、共振价键体或是自旋液态引起的格点体系自发对称破缺导致的，进而推测中间相可能存在栅格态（columnar dimer）、块价键体、有能隙或者无能隙的自旋液态^{［4~6，17］}。该模型目前最新的研究表明，中间相可以分成两个结构不同部分，并且中间可能存在无能隙的自旋液态相将这两个部分分开，而不仅仅是一个相变点^{［4，6，17］}。

受限于尺度效应，目前中间相很多问题都处在争议和讨论之中。一些数值方法本身引入的一些先入为主的设定，也是进一步深入理解中间相本质的很大障碍。虽然目前没有一个天才的方法可以一次解决该模型所有方面的问题，但对于特定的方面，可以找到相适应的方法进行有效的研究。一个相当不错的方法就是，通过数值方法计算能级或能谱来确定中间相内精细相变的位置。由于能级交叉点对尺度的依赖是平滑的，可以较可靠地外推到热力学极限^{［10，12，14，15，17］}。

为了确定体系结构，聚合化算符^{［6~10，12，13］}也是一个实用且有效的工具。二聚化算符是最简单的结构化算符，涉及的格点数最少，而且它是一个标量，可以很方便地用来确定中间相的二聚化结构成分。

本文主要研究周期性边界条件下阻挫自旋1/2海森堡反铁磁J₁-J₂方格子的基态和低激发态。根据能级交叉和二聚化算符，确定中间相的相变点和能态结构，并寻找自旋液态。本文组织如下：第1部分报道能态计算、能级交叉，还有忠实度、磁结构因子、磁化强度和自旋关联等，给出相应的相变点位置；第2部分通过结构化算符表征自旋液态结构。

1　 基本能态结构

二维海森堡模型是一个典型的多体问题，其哈密顿量可以写为

（1）

式中所有记号全部按照约定习惯：{\mathit{S}}_{\mathit{i}}是格点i处自旋\mathrm{1}/\mathrm{2}算符；对和\ll i,k\gg的求和分别对最近邻（nearest-neighbor，NN）和次近邻（next-nearest-neighbor，NNN）进行，格点间每个键计算一次。为简便计，取J₁为单位\mathrm{1}，反铁磁模型下J₂>0为耦合参数，或者取无量纲数g=J₂/J₁作为耦合参数。

对于格点数N=\mathrm{16}的方格子体系，通过正交化波算符^［17］配合迭代法，得到体系的基态和第一激发态，发现在激发态出现能态简并。为了处理简并能态，进一步通过三分法和共轭梯度法进行逆迭代，在低激发态能级间作地毯式能量扫描，确定体系的能级结构。然后在本征能量下用逆迭代计算过完备的基矢组，利用正交化方法取总自旋量子数S子空间正交完备基，将总自旋z分量在S子空间中精确对角化，得到哈密顿量\mathit{H}、总自旋量子数S、总自旋z分量量子数S_T^z的共同本征态。

依据参数g划分能态，计算结果通过表\mathrm{1}给出。基态\mathrm{G}\mathrm{S}全部为单态；第一激发态\mathrm{L}\mathrm{E}\mathrm{V}\mathrm{1}在\mathrm{0.41}处分开成两个部分，\mathrm{0.00}~\mathrm{0.40}的区域属于自旋三重态，\mathrm{0.41}~\mathrm{0.99}的区域属于自旋单态；第二激发态\mathrm{L}\mathrm{E}\mathrm{V}\mathrm{2}将\mathrm{L}\mathrm{E}\mathrm{V}\mathrm{1}的两个区域进一步分成5个更精细的区域，\mathrm{0.00}~\mathrm{0.21}属于自旋五重态，\mathrm{0.22}~\mathrm{0.40}属于自旋单态，\mathrm{0.41}~\mathrm{0.49}属于自旋三重态，\mathrm{0.50}~\mathrm{0.71}属于自旋单态，\mathrm{0.72}~\mathrm{0.99}属于自旋三重态。特别地，根据{\mathbb{Z}}_{\mathrm{2}}几何自旋液态的拓扑二重简并特征^{［3，5，17］}，发现在\mathrm{L}\mathrm{E}\mathrm{V}\mathrm{2}中区域\mathrm{0.50}~\mathrm{0.52}存在两套简并的自旋单态相互正交，以它为基矢可以形成一个二维流形，流形上的每一个态都是自旋单态，它很可能就是将中间相分成栅格态和块价键体^［18］的{\mathit{\mathbb{Z}}}_{\mathrm{2}}自旋液态；区域\mathrm{0.72}~\mathrm{0.99}存在两套简并的自旋三重态，它则是由于条纹相在激发态C_{\mathrm{4}}自发对称破缺的消失造成的，条纹不再自发选择x或z方向中的一个，而是都有可能且可以任意叠加。

格点数N=\mathrm{16}体系的两个自旋单态和一个最低的自旋三重态相对于基态的能隙随着无量纲参数g的变化过程在图\mathrm{1}中进行了展示。自旋单态E_{\mathrm{1}}\left(S=\mathrm{0}\right)和三重态E_{\mathrm{1}}(S=\mathrm{1})在临界点g_{\mathrm{c}\mathrm{1}}=\mathrm{0.41}的位置发生交叉，如图\mathrm{1}(\mathrm{a})所示，而自旋单态E_{\mathrm{2}}(S=\mathrm{0})和三重态E_{\mathrm{1}}(S=\mathrm{1})在临界点g_{\mathrm{c}\mathrm{2}}=\mathrm{0.71}的位置发生交叉，如图\mathrm{1}(\mathrm{d})所示，由g_{\mathrm{c}\mathrm{1}}和g_{\mathrm{c}\mathrm{2}}之间的区域给出中间相的位置。在中间相区域，自旋单态E_{\mathrm{1}}\left(S=\mathrm{0}\right)相对基态的能隙随着g的增加衰减很快，如图\mathrm{1}(\mathrm{a})中对应的曲线E_{\mathrm{1}}\left(S=\mathrm{0}\right)所示。自旋单态E_{\mathrm{2}}(S=\mathrm{0})与三重态E_{\mathrm{1}}(S=\mathrm{1})在临界点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}发生交叉，如图\mathrm{1}(\mathrm{b})所示，并在g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.53}的位置发生交叠，图\mathrm{1}(\mathrm{c})所示，这两个临界点将中间相分成三个精细区域。在图\mathrm{1}(\mathrm{b})中两曲线交叉点的右侧，区域g=\mathrm{0.50}~\mathrm{0.52}有两组简并的自旋单态，可能是具有二度拓扑简并的自旋液态相，而在g=\mathrm{0.53}~\mathrm{0.71}的区域只有无简并的自旋单态。

在中间相区域，自旋液态相的左边和右边分别是栅格态和价键体。需要说明的是，低激发态可以反映基态的性质^{［14，17］}，并且基态相变可以通过给定总自旋量子数子空间内部或者不同总自旋量子数子空间正交化能级间的交叉来表达。根据这种交叉模式，在热力学极限下，分别可能出现三种情况，如图\mathrm{2}所示：(\mathrm{a}) 无能隙的自旋液态^［4］，上下能级有一个重叠的平台；(\mathrm{b}) 无能隙的一个临界点^［18］，上下能级的顶点重叠；(\mathrm{c}) 有能隙的自旋液态^［5］，上下能级的谷峰位置有一段重叠区域。

图\mathrm{1}中能级交叉是从能量角度给出的，还可直接从量子态角度考虑，利用基态及低激发态的忠实度^［19］

F\left(g,\delta \right)=\left|〈\Psi \left(g\right)|\Psi \left(g+\delta \right)〉\right|（2）

如图\mathrm{3}(\mathrm{b})所示，基态在{g_{c\mathrm{2}}}^{\text{'}}=\mathrm{0.63}的位置忠实度有一个向下的峰，此处对应一级相变，在热力学极限下，忠实度仅在一点发生跳变，其值为\mathrm{0}，而其余值为\mathrm{1}。如图\mathrm{3}(\mathrm{a})所示，第一激发态在{g_{c\mathrm{2}}}^{\text{'}}=\mathrm{0.63}的位置忠实度有一个向下的凹陷，激发态在一定程度上反映了基态一级相变的信息。图\mathrm{3}(\mathrm{a})中第一激发态的忠实度在g_{c\mathrm{1}}=\mathrm{0.41}的位置发生跳变，并断裂成明显的两段线条，因为发生了自旋单态E_{\mathrm{1}}\left(S=\mathrm{0}\right)和三重态E_{\mathrm{1}}\left(S=\mathrm{1}\right)的能级交叉，如图\mathrm{1}(\mathrm{a})所示，而在基态g_{c\mathrm{1}}=\mathrm{0.41}处的忠实度是连续的，如图\mathrm{3}(\mathrm{b})所示，所以此处对应二级相变。如图\mathrm{3}(\mathrm{c})所示，第二激发态在g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}和g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.53}，忠实度发生跳变，将中间相断开为三段线条，因为发生了自旋单态E_{\mathrm{2}}(S=\mathrm{0})与三重态E_{\mathrm{1}}(S=\mathrm{1})在临界点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}的交叉和g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.53}的交叠，分别如图1（b）和（c）所示，这对应于中间相内部的精细结构相变。

对于中间相与\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}相和条纹相之间的相变点，其他研究者用不同的处理方法进行了计算。簇耦合方法（coupled cluster method，CCM）给出g_{c\mathrm{1}}=\mathrm{0.454},g_{c\mathrm{2}}=\mathrm{0.588}^［16］，有限尺度精确对角化方法给出g_{c\mathrm{1}}=\mathrm{0.40},g_{c\mathrm{2}}=\mathrm{0.65}^［10］，Lanczos精确对角化方法给出g_{c\mathrm{1}}=\mathrm{0.35},g_{c\mathrm{2}}=\mathrm{0.66}^［14］，密度矩阵重整化群（density matrix renormalization group，DMRG）给出g_{c\mathrm{1}}=\mathrm{0.41},g_{c\mathrm{2}}=\mathrm{0.62}^［5］，蒙特卡洛（Monte Carlo， MC）给出g_{c\mathrm{1}}=\mathrm{0.4},g_{c\mathrm{2}}=\mathrm{0.6}^［4］。本文的结果g_{c\mathrm{1}}=\mathrm{0.41}与其他方法的结果相近，说明通过第一与第二激发态的能级交叉可以较好地反映二级相变点的位置，而g_{c\mathrm{2}}=\mathrm{0.71}与其他方法的结果出入比较大，说明当有限尺度下基态与第一激发态能级交叉没有表现出来时，用第二和第三激发态的能级交叉来指示一级相变点，存在较大的尺度依赖性。

对于中间相内的精细相变点，簇耦合方法给出g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49},g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.58}^［16］，蒙特卡洛方法只识别出一个相变点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.50}^［4］，密度矩阵重整化群方法识别出g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.46},g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.52}^［5］。本文的结果g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49},g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.53}所能确定自旋液态区域相比其他方法要小，而且在N=\mathrm{16}体系下表现为能隙很小的{\mathbb{Z}}_{\mathrm{2}}自旋液态。中间相内的相变是目前争议最大的问题，不仅对相变点的位置有所保留，甚至中间相有没有自旋液态，自旋液态有没有能隙都在激烈的讨论中。产生这些分歧的主要原因是计算体系的尺度不同，以及伴随处理这些体系近似方法的差异。本文计算的是有限尺度N=\mathrm{16}下体系表现出来的特征，并进一步计算了N=\mathrm{18}和\mathrm{20}尺度下能级交叉表现出来的趋势。

采用Marshall-Peierls态比重方法^{［11，14］}，计算了基态在三个样本态g\text{'}=\mathrm{0.00,0.50,0.99}上的投影大小\mathrm{S}\mathrm{W}，并在图\mathrm{4}中作出了三条对应的态比重\mathrm{S}\mathrm{W}随参数g变化的曲线。发现反铁磁Neel态g\text{'}=\mathrm{0.00}在的区域维持着很高的比重(>\mathrm{80}\mathrm{\%})，在g=\mathrm{0.50}~\mathrm{0.60}的区域Neel态的比重随着g的增大迅速衰减，当g到达\mathrm{0.71}的时候Neel态的比重已经很低。相反，条纹态g\text{'}=\mathrm{0.99}在g>\mathrm{0.71}的区域维持着高比重(>\mathrm{90}\mathrm{\%})，在g=\mathrm{0.50}~\mathrm{0.71}的区域条纹态的比重随着g的减小迅速衰减，当的时候条纹态的比重已经非常低。特别地，自旋液态g\text{'}=\mathrm{0.50}的比重在g=\mathrm{0.5}附近达到峰值，往左随着g的减小自旋液态的比重逐渐衰减到\mathrm{80}\mathrm{\%}以下，而往右随着g的增大自旋液态的比重迅速衰减到\mathrm{20}\mathrm{\%}以下。

在给定量子态下，对于总自旋量子数S，对应总自旋算符的平方可以表示为

〈{{\mathit{S}}_{\mathrm{T}}}^{\mathrm{2}}〉=〈{\left(\sum _i{\mathit{S}}_{\mathit{i}}\right)}^{\mathrm{2}}〉=S\left(S+\mathrm{1}\right)（3）

以及总自旋算符z方向分量为

〈S_{\mathrm{T}}^z〉=〈\sum _i{S}_i^z〉=\sum _i〈S_i^z〉（4）

格点自旋的有效磁矩〈S_i^z〉和总自旋z分量量子数S_{\mathrm{T}}^z可以相互验证，结果表明二者是一致的。自旋单态格点x,y和z方向有效磁矩均为零，随着S_{\mathrm{T}}^z绝对值的增长，z方向有效磁矩也随之增长。

两个自旋格点间的关联函数可以表示为C_{ij}=〈{\mathit{S}}_i·{\mathit{S}}_j〉，根据哈密顿量的表达式，能量可以通过自旋关联表达出来

（5）

因为该模型中能量的本质就是自旋在两种耦合模式下的关联相互作用，所以能级交叉也可以通过自旋关联的交叉来反映，如图\mathrm{5}所示。对于基态，自旋关联的值随距离平方的增加震荡衰减，如图\mathrm{5}(\mathrm{b})所示，进一步作出自旋关联的绝对值随着距离的变化曲线，如图\mathrm{5}(\mathrm{c})所示，基本上呈指数衰减。在参数g空间中，自旋关联的交叉给出临界点，如图\mathrm{5}(\mathrm{a})所示：最近邻r^{\mathrm{2}}=\mathrm{1}与次近邻r^{\mathrm{2}}=\mathrm{2}的自旋关联交叉于{g_{c\mathrm{1}}}^{\text{'}}=\mathrm{0.63}的位置，给出一阶相变点的位置，这与其他方法^{［4，5，10，14，16］}给出的一阶相变点比较吻合；阻挫相互作用下，在相变点{g_{c\mathrm{1}}}^{\text{'}}的左边\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}序占优，而在相变点{g_{c\mathrm{1}}}^{\text{'}}右边条纹序占优；该体系中最远关联r^{\mathrm{2}}=\mathrm{8}与次远关联r^{\mathrm{2}}=\mathrm{5}有两个交叉点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}和{g_{c\mathrm{1}}}^{\text{'}}，在临界点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}处两种关联均接近于零；在r^{\mathrm{2}}=\mathrm{8}与r^{\mathrm{2}}=\mathrm{5}谷峰交叠区域的左侧，二者基本上等大反号，而在右侧g>{g_{c\mathrm{1}}}^{\text{'}}，随着g的增大，次远关联r^{\mathrm{2}}=\mathrm{5}的数值很平滑地趋于零，最远关联r^{\mathrm{2}}=\mathrm{8}反而较快地增长到超过g=\mathrm{0.00}时的水平；次近邻r^{\mathrm{2}}=\mathrm{2}与次远r^{\mathrm{2}}=\mathrm{5}的自旋关联在临界点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.54}处发生交叉，在交叉点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}的左侧，随着g的增大，中程r^{\mathrm{2}}=\mathrm{4}与次近邻r^{\mathrm{2}}=\mathrm{2}的自旋关联一致地下降，而在交叉点的右侧中程r^{\mathrm{2}}=\mathrm{4}自旋关联不降反升，与最远关联r^{\mathrm{2}}=\mathrm{8}保持一致。

根据关联函数，引入波矢\vec{k}，通过傅里叶变换可以得到磁结构因子^{［5，10，14］}

M_N^{\mathrm{2}}\left(\vec{k}\right)=\frac{\mathrm{1}}{N}\sum _{i,j}{C}_{ij}{\mathrm{e}}^{i\vec{k}·\left({\vec{r}}_i-{\vec{r}}_j\right)}（6）

通过图\mathrm{6}可以发现，{g_{c\mathrm{2}}}^{\text{'}}=\mathrm{0.63}与g=\mathrm{0.99}态的自旋关联在傅里叶变换下得到的结构因子，差别可以明显地表现出来，而g_{c\mathrm{2}}=\mathrm{0.71}与g=\mathrm{0.99}态的结构因子，差别不明显。说明傅里叶变换会将态的局域关联细节平滑掉，进而更加展示出整体性特点。对比峰值位于{\vec{k}}_{\mathrm{1}}=(\pi ,\pi )的\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}态g=\mathrm{0.30}的结构因子，可以发现g=\mathrm{0.51}时闭合的峰环打开，形成明显的谷沟，这说明两个态属于不同的类型，后者的定域性比前者弱，流动性强于前者，某种程度上支持了自旋液态的猜想。g=\mathrm{0.57}处的结构因子在{\vec{k}}_{\mathrm{1}}附近形成了一个微凸的平台，它往{\vec{k}}_{\mathrm{2}}=(\pi ,\mathrm{0})方向下降平缓并逐渐稳定在一个较高的水平，相反地容易看出{\vec{k}}_{\mathrm{3}}=(\mathrm{0,0})处存在明显值为0的低谷，这表明往波矢{\vec{k}}_{\mathrm{3}}方向的衰减速度比较快。g=\mathrm{0.63}时，中间平台{\vec{k}}_{\mathrm{1}}的结构强度已经转移到了波矢{\vec{k}}_{\mathrm{2}}的位置，最终波矢{\vec{k}}_{\mathrm{1}}处形成一个具有四重对称的鞍点，这对应具有C_{\mathrm{4}}对称性的块价键体。在临界点g_{c\mathrm{2}}=\mathrm{0.71}处，波矢{\vec{k}}_{\mathrm{2}}的位置所获得的结构强度已经形成明显的峰值，结果导致一个C_{\mathrm{4}}对称的鞍点{\vec{k}}_{\mathrm{1}}逐渐转变成位置分布具有C_{\mathrm{4}}对称的4个鞍点{\vec{k}}_{C\mathrm{4,1}}=\left(\mathrm{\pi }/\mathrm{2},\mathrm{\pi }/\mathrm{2}\right),{\vec{k}}_{C\mathrm{4,2}}=(\mathrm{3}\mathrm{\pi }/\mathrm{2},\mathrm{\pi }/\mathrm{2}),{\vec{k}}_{C\mathrm{4,3}}=(\mathrm{\pi }/\mathrm{2,3}\mathrm{\pi }/\mathrm{2})，{\vec{k}}_{C\mathrm{4,4}}=(\mathrm{3}\mathrm{\pi }/\mathrm{2,3}\mathrm{\pi }/\mathrm{2})。随着g的继续增大，中间平台\mathrm{0.5}强度的等高线开始和四角\mathrm{0.5}强度的等高线交汇，鞍点{\vec{k}}_{C\mathrm{4},i}(i=\mathrm{1,2},\mathrm{3,4})转变为峰值间衔接的不解析点，即当条纹态g=\mathrm{0.99}选择了峰值位置波矢{\vec{k}}_{\mathrm{2}}=\left(\mathrm{\pi },\mathrm{0}\right)\mathrm{或}\left(\mathrm{0},\mathrm{\pi }\right)中的一个时，就很难转变成另外一个，这解释了条纹相的自发C_{\mathrm{4}}对称破缺。

\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}相和条纹相分别在{\vec{k}}_{\mathrm{1}}=(\mathrm{\pi },\mathrm{\pi })和{\vec{k}}_{\mathrm{2}}=(\mathrm{\pi },\mathrm{0})出现最大峰值。分别计算它们磁序结构的磁化强度m_{\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}}^{\mathrm{2}}和m_{\mathrm{S}\mathrm{t}\mathrm{r}}^{\mathrm{2}}，有

m_{\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}}^{\mathrm{2}}=\frac{\mathrm{1}}{N}{M}_N^{\mathrm{2}}\left(\mathrm{\pi },\mathrm{\pi }\right)（7）

m_{\mathrm{S}\mathrm{t}\mathrm{r}}^{\mathrm{2}}=\frac{\mathrm{1}}{N}{M}_N^{\mathrm{2}}\left(\mathrm{\pi },\mathrm{0}\right)（8）

图\mathrm{7}根据m_{\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}}^{\mathrm{2}}与m_{\mathrm{S}\mathrm{t}\mathrm{r}}^{\mathrm{2}}的交叉点给出一级相变临界点的位置{g_{c\mathrm{2}}}^{\text{'}}=\mathrm{0.63}，与图\mathrm{5}中自旋关联所给的临界点一致。\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}序强度从g_{c\mathrm{1}}=\mathrm{0.41}附近开始随着g的增大迅速减小，在g_{c\mathrm{2}}=\mathrm{0.71}之后逐渐趋于零，这说明\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}序经过阻挫不断增大的无序致有序过程后衰减得迅速且彻底。条纹序强度反而由于无序致有序过程从g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}附近开始随着g的增大迅速增长，在g_{c\mathrm{2}}之后增长开始变得平缓，之后稳定在较高水平。这说明随着g的增大，中间相的存在为激烈的无序致有序过程产生条纹序提供了一种过渡。

通过计算基态能量及其对参数g的一、二阶导数，如图\mathrm{8}(\mathrm{a})所示，基态能量在中间相达到极大值，能量的一阶导数在{g_{c\mathrm{2}}}^{\text{'}}=\mathrm{0.63}处有一个陡峭的下降，能量的二阶导数在{g_{c\mathrm{2}}}^{\text{'}}处有一个极小值，对应一阶相变点。在热力学极限下基态能量在{g_{c\mathrm{2}}}^{\text{'}}处不连续，所以{g_{c\mathrm{2}}}^{\text{'}}是一个在有限体系下没有显现出来的基态能级交叉点。图\mathrm{8}(\mathrm{b})给出N=\mathrm{16}体系的基本能级结构和能隙，可以看到各能级间的交叉以及能隙的零点位置。第一激发态相对于基态的能隙在g_{c\mathrm{2}}附近达到极小值接近于零，而且基态能量相对于参数g的二阶导数在{g_{c\mathrm{2}}}^{\text{'}}处给出一个极小峰，这两个点给出了有限尺度下一阶相变点可能存在的大致范围{g_{c\mathrm{2}}}^{\text{'}}~g_{c\mathrm{2}}。根据分级平均场^［20］和簇密度矩阵内嵌^［21］得到热力学极限下临界点{g_{c\mathrm{2}}}^{\text{'}}附近基态能量不解析，蒙特卡洛方法^［4］得到临界点{g_{c\mathrm{2}}}^{\text{'}}附近存在块价键态和条纹态的能级交叉，表明此处发生的是一阶相变。图\mathrm{9}通过计算N=\mathrm{18,20}这样更大体系的能级发现，当体系到达\mathrm{20}个格点的时候，就可以比较准确地看到一阶相变的能级交叉。

本文使用的计算工具是第一作者自主编写的\mathrm{C}程序：迭代算法的核心是抽象算符作用和矢量化，避开存储超大型矩阵（{\mathrm{2}}^{\mathrm{32}}个双精度矩阵元，对于N=\mathrm{16}体系），哈密顿量的作用过程通过独立的翻转逻辑函数来实现，翻转操作对象为数组化的态矢量；逆迭代算法的核心是共轭梯度法，对于海森堡模型中算符两两作用的二次收敛问题，其精度和计算效率十分可靠。基于共享的原则，作者将整个程序包放在\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h}\mathrm{G}\mathrm{a}\mathrm{t}\mathrm{e}科研意识形态网站（https：//www.researchgate.net）上开放阅读交流。

在计算精度方面，我们采用迭代前后变分能量E_k=〈{\psi }_k|\mathit{H}|{\psi }_k〉的相对误差\frac{E_{k+\mathrm{1}}-E_k}{E_k}进行精度控制，误差限设置为{\mathrm{10}}^{-\mathrm{8}}（一般精度模式）和{\mathrm{10}}^{-\mathrm{11}}（高精度模式）。程序中还做了随机误差验证，其判据为随机误差矢量\left|\mathrm{e}\mathrm{r}\mathrm{r}\right.〉=\frac{\mathit{H}}{E_k}\left|{\psi }_k\right.〉-\left|{\psi }_k\right.〉，结果表明其模||\mathrm{e}\mathrm{r}\mathrm{r}||的波动范围为{\mathrm{10}}^{-\mathrm{15}}~{\mathrm{10}}^{-\mathrm{13}}。为了防止哈密顿量\mathit{H}将\left|\mathrm{e}\mathrm{r}\mathrm{r}\right.〉翻转到一个能量贡献起主要作用的方向而产生较大的误差，程序还计算了交叠误差〈{\psi }_k|\mathit{H}|\mathrm{e}\mathrm{r}\mathrm{r}〉，通过评估发现其波动范围为{\mathrm{10}}^{-\mathrm{12}}~{\mathrm{10}}^{-\mathrm{13}}是合理的。

对于计算效率，一般精度模式下每个能态的计算耗时\mathrm{3}分钟，高精度模式下每个能态的计算耗时\mathrm{1}小时。一般精度模式主要用于三分法中对能级作地毯式扫描，精确计算能态则要用到高精度模式。

2　 中间相结构的表征

由于相互作用的传递，不仅不同距离的格点某种程度地关联起来，而且一些格点簇也会呈现出关联，进而产生不同的局域化结构。二聚化算符^［4，5］可以定义为一种沿着\zeta =x或z方向的键作用算符{D_i}^{\zeta }={\mathit{S}}_i\cdot {\mathit{S}}_{i+\zeta }，是一种最简单的结构化算符，并且是一个标量，使用十分方便。类似于用自旋算符描述自旋关联一样，同样可以利用二聚化算符描述格点簇的二聚化关联

C_D\left({D_i}^{\alpha },{D_j}^{\beta }\right)=〈{D_i}^{\alpha }{D_j}^{\beta }〉-〈{D_i}^{\alpha }〉〈{D_j}^{\beta }〉（9）

类比自旋格点的结构因子，还可以构造二聚化结构因子^［4，17］

M_D^{\alpha \beta }\left(\vec{k}\right)=\frac{\mathrm{1}}{N}\sum _{i,j}{\mathrm{e}}^{i\vec{k}·\left({\vec{r}}_i-{\vec{r}}_j\right)}{C}_D\left({D_i}^{\alpha },{D_j}^{\beta }\right)（10）

取x方向二聚化得到格点簇的聚合强度

m_{D,x}^{\mathrm{2}}=\frac{\mathrm{1}}{N}{M}_D^{xx}\left(\vec{k}\right)（11）

图\mathrm{10}给出4个典型参数g处的二聚化结构因子，其峰值出现在{\vec{k}}_{\mathrm{2}}=(\mathrm{\pi },\mathrm{0})或者{\vec{k}}_{\mathrm{3}}=(\mathrm{0,0})，其中自旋液态的峰出现在{\vec{k}}_{\mathrm{2}}和{\vec{k}}_{\mathrm{3}}两个位置，且在波矢{\vec{k}}_{\mathrm{2}}处的峰值更高，块价键态和条纹态的峰出现在{\vec{k}}_{\mathrm{3}}。

通过图\mathrm{11}分别给出波矢{\vec{k}}_{\mathrm{2}}和{\vec{k}}_{\mathrm{3}}的二聚化格点簇聚合强度，静态二聚化聚合强度{\vec{k}}_{\mathrm{3}}曲线从g=\mathrm{0.50}附近开始随着g的增大迅速增长，在g=\mathrm{0.65}附近达到极大值，随后下降但仍处在一个较高水平。这是由于在\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}序主导的态中，沿水平和垂直方向的自旋二聚化聚合强度较弱，而在条纹序主导的态中，二聚化可以沿水平或垂直方向固化聚合，所以二聚化聚合强度仍处于较高水平。静态二聚化聚合强度反映的是类条纹（stripe-like）的二聚化序，它的二聚化簇同向排列，有类似条纹的形态。动态二聚化聚合强度{\vec{k}}_{\mathrm{2}}曲线随着g的增加平滑增长，在g=\mathrm{0.50}附近达到极大值，随后有一个陡峭的下降并逐渐趋于零。这是由于引入波矢{\vec{k}}_{\mathrm{2}}=(\mathrm{\pi },\mathrm{0})的二聚化运动对\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}序和条纹序在特定方向的排列会造成破坏，而自旋液态中二聚化单态更加有利于波矢{\vec{k}}_{\mathrm{2}}的传播。这种动态二聚化序是类\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}（Neel-like）的二聚化序，它的二聚化簇间隔反向，有类似\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}自旋的形态。{\vec{k}}_{\mathrm{2}}和{\vec{k}}_{\mathrm{3}}曲线表明，在中间相g=\mathrm{0.50}附近的区域，动态和静态二聚化聚合强度都取得较大值，且强度值也很接近。在这一区域，间隔反向和间隔同向的二聚化簇都可以显著地共存叠加，说明二聚化簇的排列方向并不沿特定方向，可以自由地流动，其有效二聚化强度可以在\mathrm{0}和峰值之间变化，表现了自旋液态的特点。在g=\mathrm{0.50}的左侧，两类二聚化簇聚合强度随着g的减小一致地衰减，对照表\mathrm{1}的结果可以发现，基态倾向于由单态激发到三重态，二聚化簇的稳定性较差。而在g=\mathrm{0.50}的右侧，类条纹的二聚化序随着g的增加，很快地获得了支配性的优势，当类\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}的二聚化序在g=\mathrm{0.71}之后衰减到接近\mathrm{0}的时候，类条纹的二聚化序便凝聚为条纹排列的自旋。

3　 结 语

考虑到格点系的哈密顿量\mathit{H}(g)是无量纲参数g的函数，我们追踪\mathit{H}(g)的基态能量随着g的演化函数。对于有限体系，基态能量对参数g是光滑的、解析的。出现违反这一规则的情况，主要是由于g耦合了两个相互对易的部分。在给定g的邻域内，当近邻项H_{\mathrm{0}}和次近邻项H_{\mathrm{1}}变化不大时，能量曲线在该g值附近可以近似用一个切线E_H=A_{H_{\mathrm{0}}}+g{B}_{H_{\mathrm{1}}}来拟合，其中A_{H_{\mathrm{0}}}和B_{H_{\mathrm{1}}}分别是与H_{\mathrm{0}}和H_{\mathrm{1}}相关的常数参量。通过体系总自旋量子数可以看出，虽然本征能量随着g改变，但是本征函数却和g无关，这样便会出现能级交叉。一个激发态在这个临界点同时也是基态，这种交叉破坏了基态能量对g的解析性，形成奇异点。由于哈密顿量是对一个体系相互作用的一般描述，所以不仅仅是基态能量，低激发态能量对g的解析性也会由于能级交叉的出现而产生奇异点。

一个在有限体系下没能表现出来的能级交叉，随着格点体系的增大，能级交叉将会越来越显著，最后在无限大体系下形成一个尖锐的临界点。热力学极限下，我们可以找出任何基态能量不解析的奇异点。这些点就是体系的量子相变点，它们可以是有限体系下已经出现的能级交叉点，也可能是有限体系下没有表现出来的能级交叉点。

本文计算了自旋\mathrm{1}/\mathrm{2}海森堡反铁磁J₁-J₂方格子的基态和低激发态及相应的本征能量，并在总自旋量子数S子空间内精确对角化，进一步得到哈密顿量\mathit{H}、总自旋量子数S、总自旋z分量量子数S^z的共同本征态。根据能级交叉和二聚化算符，近似给出中间相内的相变点g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{1}}=\mathrm{0.49}和g_{\mathrm{c}\mathrm{i}\mathrm{n}\mathrm{2}}=\mathrm{0.53}。计算结果表明，在中间相g=\mathrm{0.41}~\mathrm{0.71}的左右两侧分别为\mathrm{N}\mathrm{e}\mathrm{e}\mathrm{l}相(g=\mathrm{0.00}~\mathrm{0.40})和条纹相(g=\mathrm{0.72}~\mathrm{0.99})；中间相被低激发的自旋单态S=\mathrm{0}和自旋三重态S=\mathrm{1}的能级交叉分隔成了g=\mathrm{0.41}~\mathrm{0.49}和g=\mathrm{0.50}~\mathrm{0.71}两个区域，并且在g=\mathrm{0.50}~\mathrm{0.52}的区域发现低激发态拓扑二重简并的{\mathbb{Z}}_{\mathrm{2}}几何自旋液态。

参考文献

[1] MANOUSAKIS E. The spin-1/2 Heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides ［J］. Reviews of Modern Physics， 1991， 63（1）： 1. DOI： 10.1103/RevModPhys.63.1 .

[2] KNOLLE J， MOESSNER R. A field guide to spin liquids［J］. Annual Review of Condensed Matter Physics， 2019， 10： 451-472. DOI： 10.1146/annurev-conmatphys-031218-013401 .

[3] HUANG C Y， LIN F L. Topological order and degenerate singular value spectrum in two-dimensional dimerized quantum Heisenberg model ［J］. Physical Review B， 2011， 84（12）： 125110. DOI： 10.1103/PhysRevB.84.125110 .

[4] MORITA S， KANEKO R， IMADA M. Quantum spin liquid in spin 1/2 J₁-J₂ Heisenberg model on square lattice： Many-variable variational monte carlo study combined with quantum-number projections ［J］. Journal of the Physical Society of Japan， 2015， 84（2）： 024720. DOI： 10.7566/JPSJ.84.024720 .

[5] JIANG H C， YAO H， BALENTS L. Spin liquid ground state of the spin-1/2 square J₁-J₂ Heisenberg model ［J］. Physical Review B， 2012， 86（2）： 024424. DOI： 10.1103/PhysRevB.86.024424 .

[6] CAPRIOTTI L， SORELLA S. Spontaneous plaquette dimerization in the J₁-J₂ Heisenberg model ［J］. Physical Review Letters， 2000， 84（14）： 3173. DOI： 10.1103/PhysRevLett.84.3173 .

[7] DAGOTTO E， MOREO A. Exact diagonalization study of the frustrated Heisenberg model： A new disordered phase ［J］. Physical Review B， 1989， 39（7）： 4744. DOI： 10.1103/PhysRevB.39.4744 .

[8] MAMBRINI M， LÄUCHLI A， Poilblanc D， et al. Plaquette valence-bond crystal in the frustrated Heisenberg quantum antiferromagnet on the square lattice ［J］. Physical Review B， 2006， 74（14）： 144422. DOI： 10.1103/PhysRevB.74.144422 .

[9] DAGOTTO E， MOREO A. Phase diagram of the frustrated spin-1/2 Heisenberg antiferromagnet in 2 dimensions ［J］. Physical Review Letters， 1989， 63（19）： 2148. DOI： 10.1103/PhysRevLett.63.2148 .

[10] SCHULZ H J， ZIMAN T A L. Finite-size scaling for the two-dimensional frustrated quantum Heisenberg antiferromagnet ［J］. Europhysics Letters， 1992， 18（4）： 355. DOI：10.1209/0295-5075/18/4/013/pdf .

[11] RETZLAFF K， RICHTER J， IVANOV N B. Analysis of the ground-state wave function of the J₁-J₂ quantum Heisenberg antiferromagnet ［J］. Zeitschrift für Physik B Condensed Matter， 1993， 93（1）： 21-31. DOI: 10.1007/bf01308803 .

[12] DARRADI R， DERZHKO O， ZINKE R， et al. Ground state phases of the spin-1/2 J₁-J₂ Heisenberg antiferromagnet on the square lattice： A high-order coupled cluster treatment ［J］. Physical Review B， 2008， 78（21）： 214415. DOI： 10.1103/PhysRevB.78.214415 .

[13] ISAEV L， ORTIZ G， DUKELSKY J. Hierarchical mean-field approach to the J₁-J₂ Heisenberg model on a square lattice ［J］. Physical Review B， 2009， 79（2）： 024409. DOI： 10.1103/PhysRevB.79.024409 .

[14] RICHTER J， SCHULENBURG J. The spin-1/2 J₁-J₂ Heisenberg antiferromagnet on the square lattice： Exact diagonalization for N= 40 spins ［J］. The European Physical Journal B， 2010， 73（1）： 117-124. DOI： 10.1140/epjb/e2009-00400-4 .

[15] RICHTER J， DARRADI R， SCHULENBURG J， et al. Frustrated spin-1/2 J₁-J₂ Heisenberg ferromagnet on the square lattice studied via exact diagonalization and coupled-cluster method ［J］. Physical Review B， 2010， 81（17）： 174429. DOI： 10.1103/PhysRevB.81.174429 .

[16] RICHTER J， ZINKE R， FARNELL D J. The spin-1/2 square-lattice J₁-J₂ model： The spin-gap issue ［J］. The European Physical Journal B， 2015， 88（1）： 1-6. DOI： 10.1140/epjb/e2014-50589-x .

[17] WANG L， SANDVIK A W. Critical level crossings and gapless spin liquid in the square-lattice spin-1/2 J₁-J₂ heisenberg antiferromagnet ［J］. Physical Review Letters， 2018， 121（10）： 107202. DOI： 10.1103/PhysRevLett.121.107202 .

[18] SENTHIL T， VISHWANATH A， BALENTS L， et al. Deconfined quantum critical points ［J］. Science， 2004， 303（5663）： 1490-1494. DOI： 10.1126/science.1091806 .

[19] YOU W L， LI Y W， GU S J. Fidelity， dynamic structure factor， and susceptibility in critical phenomena ［J］. Physical Review E， 2007， 76（2）： 022101. DOI： 10.1103/PhysRevE.76.022101 .

[20] QIN J， JIE Q， FAN Z. An efficient implementation of the hierarchical mean-field theory with large block ［J］. The European Physical Journal B， 2014， 87（12）： 1-8. DOI： 10.1140/epjb/e2014-50054-0 .

[21] FAN Z， JIE Q. Cluster density matrix embedding theory for quantum spin systems ［J］. Physical Review B， 2015， 91（19）： 195118. DOI： 10.1103/PhysRevB.91.195118 .