关联量子模型的希尔伯特空间维度通常随模型所含粒子数呈指数增长, 如果不采取任何近似, 即使包含较少粒子数目的模型哈密顿矩阵也很容易超出计算机的内存空间而不能载入, 严格的数值求解更是无从谈起. 现有的关联量子模型哈密顿量近似求解手段包括变分方法^[1-6]、量子蒙特卡罗方法^[7-11]以及图展开方法^[12-15]等; 这些方法的有效性依赖于模型哈密顿量所包含的某些参数的取值范围: 在一定的参数区间范围内, 可以利用这些方法进行有效求解, 超出了这些范围, 采用这些方法求解所导致的误差往往无法控制而导致求解失败.

最近, Frame等^[16]提出了一种新方法求解关联量子模型哈密顿量: 特征矢量延拓. 其关键思路是: 尽管高维哈密顿矩阵的特征矢量在高维空间中, 当哈密顿矩阵随参数光滑变化时, 相应特征矢量的变化轨迹集中在1个低维的子空间中; 因此可以选取少量对应不同参数的特征矢量作为基矢来构建这个低维子空间, 将哈密顿量在此子空间中展开, 从而只需求解低维哈密顿矩阵. 他们将1个具有可调参数的$ 4\times 4\times 4 $晶格的Bose-Hubbard模型哈密顿量分别展开到1个由3个参数互相远离的特征矢量为基矢构成的三维子空间和1个由5个参数互相接近的特征矢量为基矢构成的五维子空间, 通过求解三维矩阵和五维矩阵均可获得高精度的近似解.

然而, 文献[16]中子空间基矢数量主要是以猜测的方式进行选取, 因此, 当模型粒子数目稍有增加, 比如从$ 4\times 4\times 4 $晶格增加到$ 6\times 6\times 6 $晶格, 获得高精度近似解需要的基矢是3个还是5个, 或者说轨迹子空间到底应该是三维还是五维, 亦或更高维, 文献中并没能给出明确答案. 确定或至少估计轨迹子空间维度大小是重要的: 一方面, 如果轨迹子空间维度大小随模型所含粒子数目的增加而迅速增加, 也就是采用特征矢量延拓方法所需的基矢数量迅速增加, 特征矢量延拓方法将很难奏效; 另一方面, 如果在特征矢量延拓计算中选用的基矢数目大于轨迹子空间的维度, 基矢有可能是线性相关的, 从而导致不合理的错误结果.

确定或估计轨迹子空间维度大小的困难在于: 尽管随参数光滑变化的高维哈密顿矩阵的特征矢量变化轨迹所在的子空间是低维的, 低维子空间的维度与多体模型所含粒子数目之间并不显含简单直接的函数关系. 本文系统性地研究了反铁磁Heisenberg链模型哈密顿矩阵的特征矢量随参数光滑变化的轨迹, 通过对特征矢量轨迹进行主成分分析画出其所在低维子空间, 指出采用特征矢量延拓求解反铁磁Heisenberg链模型哈密顿量所需基矢数目随其所包含自旋个数的增加而增加.

采用由$ L $个($ L $为偶数)自旋格点组成的, 具有周期边界条件的一维Heisenberg链模型(以下简称为模型), 如图1所示. 模型哈密顿量为

其中 ${\hat{S}}_{i}^{\alpha }$为格点的自旋$ 1/2 $算符, $ \alpha =x, y, z $; $ {J}_{i} $为相邻格点自旋之间的交换相互作用, 交替取值为$ J{\rm{c}}{\rm{o}}{\rm{s}} \theta , J{\rm{s}}{\rm{i}}{\rm{n}} \theta $, 当$ J=1 $时,

当参数θ∈[0°, 90°], ${J}_{i}\geqslant 0$, 模型为反铁磁Heisenberg链.

哈密顿量(1)的每个特征波函数可以对应于1个$ N $维矢量. 由于单个自旋的波函数可以由1个2维矢量表示, 两个维度可以分别对应于$ {S}^{z} $向上和向下的状态, 从而$ N\left(L\right)={2}^{L} $, 波函数矢量和哈密顿矩阵的维度随$ L $的增加而指数增加. 在本文中, 研究该模型的基态, 它对应的总自旋角动量$ z $分量$ {S}_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}}}^{z}=0 $, 因此自旋算符$ {S}^{z} $向上和向下的状态数都是$ L/2 $个, 所以波函数的维度$ N $减少为$N\left(L\right)= $$ {C}_{L}^{L/2}$. 例如 N(4) = 6, N(6) = 20, N(14) = 3432, N(16) = 12870, N依然随自旋个数$ L $的增加而迅速增加. 只有对自旋个数$ L $不太大的情形, 多体哈密顿量(1)的特征波函数和特征能量可通过精确对角化哈密顿矩阵的数值方法严格求解.

在本文中, 对$ L={\rm{4, 6}}, \cdots , 14 $的模型哈密顿量的基态矢量和能量进行了系统性研究. 对每一个$ L $, 将哈密顿矩阵的参数$ \theta $从0°到90°每隔1°连续变化, 通过对角化求解基态矢量和能量, 所获得的$ M $个基态矢量轨迹构成一个$ M\times N $的矩阵, 记为$ {\boldsymbol{X}}_{L} $, 其中$ M=91 $, $N={C}_{L}^{{L}/{2}}$.

我们的目标是通过对具有不同自旋个数$ L $的模型基态矢量轨迹矩阵$ {\boldsymbol{X}}_{L} $分别进行主成分分析^[17-20], 找出轨迹所在低维子空间的维度大小随自旋个数$ L $的变化关系, 从而用于指导利用特征矢量延拓法求解模型哈密顿矩阵的基态矢量和能量. 主成分分析是一种非监督学习算法, 被广泛应用于数据特征提取^[17-20]. 在轨迹矩阵$ {\boldsymbol{X}}_{{L}} $中, 每一行是1个基态波函数, 对应于1个$ N $维空间中的矢量; 根据特征矢量延拓法的思路^[16], $ {\boldsymbol{X}}_{L} $中的$ M $个基态矢量所构成的轨迹仅在部分基矢方向具有比较大的投影, 而且与其他基矢方向近似正交. 通过主成分分析找出这部分基矢, 并将原始高维基态矢量在其方向上投影, 从而可以在低维空间中直观而生动的展示轨迹图像.

由于矩阵$ {\boldsymbol{X}}_{L} $不是方阵, 首先对其进行变换$ {\boldsymbol{Y}}={\boldsymbol{X}}_{L}^{{\rm{T}}}\times {\boldsymbol{X}}_{L} $, 然后求解$ {\boldsymbol{Y}} $矩阵的特征值和特征矢量, 即

其中特征值$ {\lambda }_{n} $按从大到小的顺序排列${\lambda }_{1}\geqslant {\lambda }_{2}\geqslant $$ {\lambda }_{3}\geqslant \cdots \geqslant {\lambda }_{N}\geqslant 0$, 而对应的特征矢量${\boldsymbol{v}}_{1}, {\boldsymbol{v}}_{2}, {\boldsymbol{v}}_{3}, \cdots , $$ {\boldsymbol{v}}_{N}$是相互正交的, 称具有较大特征值$ {\lambda }_{n} $的特征矢量方向$ {\boldsymbol{v}}_{n} $为主成分方向, 原始高维波函数轨迹可以投影到相应主成分$ {\boldsymbol{v}}_{n} $方向进行展示, 即

其中$ {PC}_{n} $是波函数轨迹在第$ n\left(n={\rm{1, 2}}, 3, \cdots , N\right) $个主成分方向上的投影.

首先计算了$ L=4 $的模型基态矢量随$ \theta $变化的轨迹并将其投影到各主成分方向, 发现其仅在2个主成分方向$ {\boldsymbol{v}}_{1}, {\boldsymbol{v}}_{2} $的投影$ {PC}_{1}, {PC}_{2} $不为零, 其中

而在其他主成分方向的投影都小于$ {10}^{-8} $. 图2(a)用$ {\boldsymbol{v}}_{1}, {\boldsymbol{v}}_{2} $构成的2维空间中单位圆上的一段圆弧完整直观地展示了$ L=4 $的模型基态矢量随$ \theta $变化的轨迹, 说明尽管基态矢量是6维的, 其随$ \theta $变化的轨迹仅处于由前两个主成分方向所张成的2维空间中. 接着计算了$ L=6 $的模型基态矢量随$ \theta $变化的轨迹并将其投影到各主成分方向, 发现其不为0的主成分分量个数不再为2, 而是增加到了3. 图3(a)用这3个主成分方向构成的3维空间中单位球上的一段弧线完整直观地展示了$ L=6 $的模型基态矢量随$ \theta $变化的轨迹.

以上通过主成分分析方法获得了L = 4或6的模型基态矢量随$ \theta $变化的轨迹子空间的维度为2或3, 也就确定了采用特征矢量延拓方法求解具有其他角度$ \theta $的基态所需的基矢数目为2或3. 于是对$ L=4 $的模型哈密顿量分别在原6维空间($ N\left(4\right)= 6 $)中和2维轨迹子空间中求解基态能量, 记为$ {E}_{{\rm{4, 6}}} $和$ {E}_{{\rm{4, 2}}} $, 画在图2(b)当中. 毫不奇怪地, $ {E}_{{\rm{4, 2}}} $和$ {E}_{{\rm{4, 6}}} $几乎一样, 它们的最大绝对误差, ${\text{δ}} E= $$ \max\limits_{m=1, \cdots, M}\left|{E}_{{\rm{4, 2}}}\left({\theta }_{m}\right)-{E}_{{\rm{4, 6}}}\left({\theta }_{m}\right)\right| < {10}^{-15}$. 同样地, 对L = 6的模型哈密顿量分别在原20维空间 ($ N\left(6\right)= 20 $)中和3维轨迹子空间中求解基态能量, 并记为$ {E}_{{\rm{6, 20}}} $和$ {E}_{{\rm{6, 3}}} $, 如图3(b); $ {E}_{{\rm{6, 3}}} $和$ {E}_{{\rm{6, 20}}} $的最大绝对误差$\text{δ} E$也小于$ {10}^{-15} $. 然而, 确定采用特征矢量延拓方法求解L = 4或6的模型所需的基矢数目并不是非常有用, 因为利用计算机求解$ 6\times 6 $或$ 20\times 20 $的哈密顿矩阵并不非常困难. 我们的目标是能够确定或估计采用特征矢量延拓方法求解具有更大$ L $的模型所需的基矢数目. 随着模型中自旋数目$ L $的增加, 基态矢量的原空间维度$N\left(L\right)={C}_{L}^{ {L}/{2}}$迅速增长, 例如, $ N\left(16\right)= 12870 $, 利用计算机求解$ N\times N $的哈密顿矩阵变得困难, 因此需要找出轨迹子空间维度, 从而用于指导特征矢量延拓方法的求解计算.

尽管可以用主成分分析方法确定对L = 4或6的模型基态能量进行特征矢量延拓方法求解所需的基矢数目, 当模型变大以至于用计算机求解$ N\times N $的哈密顿矩阵变得非常困难甚至不可能时, 主成分分析方法将不再适合用于确定或估计采用特征矢量延拓方法求解所需的基矢数目; 因为首先需要进行$ N\times N $的哈密顿矩阵求解获得足够数量$ M $的基态矢量, 比如本文中, $ M= 91 $, 然后才能对这$ M $个矢量进行主成分分析. 为解决这个困难, 进一步研究了采用特征矢量延拓求解模型基态能量所需基矢数目与其所包含的自旋个数之间的关系.

对$ L=8 $的模型基态矢量随$ \theta $变化的轨迹进行主成分分析获得各主成分方向, 并分别在原70维空间($ N\left(8\right)= 70 $)中和由前$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $个主成分方向构成的子空间中求解基态能量, 记为$ {E}_{{\rm{8, 70}}} $和$ {E}_{8, {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}} $. 当$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}=6 $, $ {E}_{{\rm{8, 70}}} $和$ {E}_{8, {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}} $的最大绝对误差$\text{δ} E$小于$ {10}^{-7} $. 同样, 对$ L=10 $的模型, 当$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}=8 $, $ {E}_{{\rm{10, 252}}} $和$ {E}_{10, {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}} $的最大绝对误差$\text{δ} E$小于$ {10}^{-7} $; 对$ L=12 $的模型, 当$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}=8 $, $ {E}_{{\rm{12, 924}}} $和$ {E}_{12, {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}} $的最大绝对误差$\text{δ} E$小于$ {10}^{-7} $; 对$ L=14 $的模型, 当$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}=9 $, $ {E}_{{\rm{12, 924}}} $和$ {E}_{12, {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}} $的最大绝对误差$\text{δ} E$小于$ {10}^{-7} $. 明显地, 随着$ L $的增加, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $ 也在增加, 这个事实必须在采用特征矢量延拓进行计算时加以考虑. 将$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $作为$ L $的函数画在图4(a)中, 发现$ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $随$ L $的增加而迅速增长, 尽管其总是远远小于$ N $. $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $的快速增长对于特征矢量延拓方法的应用是非常不利的. 幸运的是, 当我们稍稍降低计算结果的精度要求, 比如从$\text{δ} E < {10}^{-7}$到$\text{δ} E < {10}^{-4}$, 尽管对$ L=4{\rm{和}}6, $ $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(L\right) $保持不变, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(8\right) $则从6降到了5, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(10\right) $从8降到了5, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(12\right) $从8降到了6, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(14\right) $从9降到了6. 从图4(a)可以看出, 如果采用$\text{δ} E < {10}^{-4}$, 当$L\leqslant 8$时, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}} $随$ L $的增加而快速增长, 而当$L\geqslant 8$时, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}{\rm{的}} $增长速度开始变得缓慢, 从而为特征矢量延拓方法的应用提供了可能性.

图4(a)清楚地展示了采用特征矢量延拓方法求解随$ \theta $变化的模型基态能量所需基矢数目随其自旋个数增加而增加的关系图像, 因此可用于指导求解具有更多自旋的模型基态. 通过延伸图4(a)中的关系曲线(图中的虚线), 可以预测, 如果采用$\text{δ} E < {10}^{-4}$, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(16\right)=6 $; 而如果采用$\text{δ} E < {10}^{-7}$, $ {N}_{{\rm{s}}{\rm{u}}{\rm{b}}}\left(16\right)=10 $. 于是分别计算了随$ \theta $变化的$ L=16 $的模型基态能量$ {E}_{{\rm{16, 6}}} $, $ {E}_{{\rm{16, 10}}} $, 如图4(b), 其中6维子空间的基矢由θ = 0°, 18°$,\cdots , $ 72°, 90°的基态矢量构成, 10维子空间的基矢由θ = 0°, 10°$,\cdots , $ 80°, 90°的基态矢量构成. 从图4(b)可以看出, $ {E}_{{\rm{16, 6}}} $和$ {E}_{{\rm{16, 10}}} $非常吻合, 计算可得$ {E}_{{\rm{16, 6}}} $和$ {E}_{{\rm{16, 10}}} $的最大绝对误差$\text{δ} E$小于$ {10}^{-4} $, 证实了图4(a)中关系曲线的预测能力.

本文通过采用主成分分析方法研究反铁磁Heisenberg链基态多体波函数矢量随参数变化的轨迹, 画出了采用特征矢量延拓方法求解其基态能量所需基矢数目随其所包含自旋个数增加而增加的函数曲线, 用于指导采用特征矢量延拓方法求解包含更多自旋的Heisenberg链模型哈密顿量. 将类似的方法应用于指导求解高维Heisenberg模型和其他关联量子模型以总结更一般的函数关系, 将在后续文章中报道.