正方光学超晶格中超冷玻色子系统的量子相变

杨硕,姜颖

(上海大学 理学院,上海 200444)

0 引言

人工光晶格超冷原子系统,由于其纯净的环境和精确的可控性[1],近年来成为研究量子多体物理的一个重要的平台[2-4],并且得到了广泛的关注和长足的发展。对于该系统中一个重要的课题,也即量子相变问题[5-7],通常有几种不同的研究方法。其中最早也是最广泛被采用的方法是平均场方法[7-8],也即解耦合近似;另一种是强耦合展开法[9]。在与相应的蒙特卡洛数值模拟结果[10]比较后不难发现,平均场计算结果偏低而强耦合展开法计算结果偏高。与以上方法不同,近年来基于场论和瑞利-薛定谔微扰论而发展起来的一种新的解析方法,即有效势方法[11],不仅可以用来精确计算简单正方光晶格中超冷玻色系统的超流态(SF)到莫特绝缘态(MI)量子相变相图[11],同时也适用于具有更为复杂的晶格结构的系统。利用此方法,我们精确计算了三角光晶格、六角光晶格及笼目晶格(Kagomé)上的量子相变相图[12],所得解析结果与相应蒙特卡洛数值结果相比误差不超过10%。除了高精度外,该有效势方法原则上是可以用来进行任意高阶微扰近似的解析计算。基于此理论和微扰论的加藤(Kato)表示[13]而提出的过程链(process chain)技术[14-15]使得我们可以借助于计算机来计算任意阶微扰近似下的相图,并通过外推法得到与蒙特卡洛数值结果几乎重合的极为精确的量子相变相图。

随着人工光晶格冷原子实验手段的不断进步,越来越多的研究不再局限于简单的正方光晶格系统,而将精力投放到了具有更加复杂的晶格几何结构[16-18]或更加复杂的相互作用[19-20]或是多组分[21]的光晶格冷原子系统。这其中一个重要且有趣的问题就是人工光学超晶格超冷玻色系统[22-28]量子相变问题。人工光学超晶格可以使用多波长激光束实现[29],对于正方光学超晶格玻色超冷原子系统,它可由以下这个推广了的Bose-Hubbard哈密顿量来刻画[30-31],

(1)

Fig.1 Illustration of square optical superlattice.The potential barrier between sublattice A and B is Δμ图1 正方光学超晶格示意图,两个子晶格A和B的势阱差为Δμ

有意思的是上述有效势方法不仅适用于均匀系统,并可进行推广而适用于多子晶格结构系统的解析研究[30,32-34]。利用推广的有效势理论,借助于过程链技术,我们已给出了Δμ=0.5 U情形下正方光学超晶格玻色冷原子系统的高精度的量子相变相图[32]。然而,对于任意势阱深度差的情形下的全域相图,还需要进一步研究。之前,研究人员基于平均场理论已详细探讨了光学超晶格玻色冷原子系统的量子相变[24],并给出了平均场下的全域相图。在本文中,借助于推广的有效势理论和过程链技术,我们将细致计算高阶近似下光学超晶格超冷玻色系统的全域量子相变相图,并进一步通过外推法给出无穷阶的结果。我们的计算结果可以为后续的蒙特卡洛数值模拟和相关实验提供参考。

1 非微扰哈密顿量基态分布

研究表明,对于Bose-Hubbard模型,在Bogoliubov近似下,即将相互作用项当作微扰进行计算,并不能得到Bose-Hubbard系统中的量子相变[35]。为了研究系统的量子相变,我们需把跃迁项当作微扰项来进行处理。为了后续的讨论,首先对(1)式中哈密顿量的非微扰部分,即

(2)

所对应的基态特性进行分析。式中U和Δμ的竞争将导致非微扰基态在μ,Δμ不同区域呈现不同的量子相。我们用A、B两个子晶格上的填充数来表示相应的本征态|nA,nB〉,不难看出,相应的本征能量为

(3)

其中

(4)

本征态|nA,nB〉能成为基态的条件是

(5)

于是,我们就有

(6)

这就给出了如图2所示的非微扰基态相图。当nA=nB时该区域即为莫特绝缘态,而当nA≠nB时该区域即为密度波态(DW)。

Fig.2 Phase diagram of spinless bosons in a square superlattice described by the Hamiltonian (2)实线表示相界,(nA,nB)表示子晶格A、B中玻色子的占据数图2 正方光学超晶格中玻色子体系非微扰哈密顿量基态分布图

2 有效势方法及相界方程

当加入跃迁项后,由于跃迁项和相互作用项不对易,这将导致量子涨落,并随着跃迁振幅的增强最终导致由局域态(MI/DW)到超流态的量子相变的发生。与之前的工作类似,为了探究光学超晶格中玻色系统的量子相变并且解析地确定相应的相界边界,我们将在光学超晶格Bose-Hubbard哈密顿量中加入外源项J=(JA,JB)T[30]以破坏其U(1)对称性,即

(7)

并将跃迁项和外源项同时视作微扰[11]。很显然,系统基态的能量必可写成跃迁振幅t和外源强度J(J†)的泰勒级数展开。由于系统的基态的局域性,系统自由能表达式中的J,J†必须成对出现。基于以上的分析,我们就可以写出自由能的展开表达式为

F(J†,J,t)=Ns[F0(t)+J†C2(t)J+J†J†C4(t)JJ+…],

(8)

其中

(9)

且

(10)

Ns表示总格点数。定义超流序参量为[30,36-37]

(11)

不难看出,序参量Ψ和外源J之间满足以下关系式

(12)

其中i∈A,B。据此对式(8)中的自由能做Legendre变换

Γ(Ψ,Ψ†,t)=F/Ns-Ψ†J-ΨJ†+…,

(13)

就得到了系统的有效势

Γ(Ψ,Ψ†,t)=F0(t)-Ψ†A2(t)Ψ+… .

(14)

我们看到,这个有效势的形式正是朗道φ4-理论的形式。

从以上Legendre变换式,我们看到

(15)

这说明没有外源的实际体系中的超流序参量对应于上面有效势的极小值点。当有效势极小值点出现在超流序参量Ψ=0处,则体系处于局域态(MI/DW),而当有效势极小值点出现在非零超流序参量处时,体系则处于超流态。根据朗道的二级相变理论,其间的相变点就发生在DetA2(t)=0处。由上面的计算可知A2和C2互为逆矩阵,也即体系量子相变相界由1/Det|C2(t)|=0来确定。我们知道,

(16)

因此,寻求相变点就是在寻求以上幂级数的收敛半径。于是我们就得到了相界方程的n-阶微扰计算结果

(17)

3 过程链算法及相图计算

(18)

而

(19)

将体系的能量修正改写成

(20)

对于αl的限制条件为

(21)

式中

(22)

由于本征态间的正交性,所以

(23)

从上面可以知道{αl}中至少有两个值为0,考虑到求迹中算符的轮转对易性,把其中的一个α为0的Sα移到式(20)的最左边,此时当αn+1≠0时,可以证明能量的修正为0,所以我们总是可以得到一个两端α均为0的加藤表示形式

(24)

将S0的表达式代入,有

(25)

(26)

观察上式可以发现每一种{αl}的组合都代表着一个从基态|m〉出发,经过若干不同的中间态然后回到基态的过程,为了表示方便,我们用{αl}的序列来代替上式求和中的项

(α1α2…αn-1)=〈m|VSα1VSα2…Sαn-1V|m〉,

(27)

这些项称为加藤项,不同的数字序列的排列对应着不同的加藤项[14]。

式(27)中每一个αl都可以取0～(n-1)中的任意值,遍历所有可能的αl的取值我们就能得到所有的加藤表示。考虑到式(23)的限制,其实可以在遍历的时候增加下面的条件,去掉对能量修正无贡献的部分:①若数字序列的两端一端为零一端不为零,则根据式(23),这一项的能量修正为零,可以舍弃;②若数字序列的两端均为零,这一项的能量可以直接保留;③若数字序列的两端均不为零,根据式(23),利用轮转对易性对数字序列进行整理,若在整理过程中出现一端为零的情况,则把含零的项按(22)式进行展开。这样我们就最终得到了首尾两端均为零的数字序列,分别用左括号和右括号代替首尾两端的零得到了所有能量修正不为零的加藤项。当然,为了简化计算,我们需要将相同能量修正的加藤项进行合并,最后保留不同能量修正的加藤项以及它的权重。

如何找出能量修正相同的加藤项并对他们进行合并?首先找到加藤项中所有的0,将0变为“-)(”,这样每一个括号内的部分就变成了所谓的基本矩阵元[14],将这些基本矩阵元按包含元素数目从小到大进行排序,然后比较相等元素数目的基本矩阵元的每一个元素的数值大小,同样按从小到大进行排序,这样我们就在不改变其对应的能量修正的前提下对加藤项的形式做了转换,从而有助于我们能寻找相同能量修正的加藤项。以加藤项(02031021)和(20021031)为例来展示这种变换:

(0,2,0,3,1,0,2,1)=-()(2)(3,1)(2,1)=-()(2)(2,1)(3,1)=(02021031),

(28)

(2,0,0,2,1,0,3,1)=-(2)()(2,1)(3,1)=-()(2)(2,1)(3,1)=(02021031).

(29)

保留完每个加藤项的最简形式,接下来对相同形式的加藤项进行合并计算权重。所有的计算过程均可以数值化进行计算,而且由于我们并没有对体系的性质做过多的限制,所以最简的加藤项一旦生成,可以适用所有的微扰计算。我们以四阶微扰为例再来完整地演示整个加藤项的生成及化简过程,见表1。

加藤项其实就是从基态出发经过一系列中间态然后又回到基态的过程。具体地来说,对于我们的问题,加藤项包含在任一格点产生一个粒子、一个格点上消灭一个粒子以及连接这两格点间的最近邻跃迁。这些步骤的排列并不唯一，但要保证最后一步结束后系统恢复到原初非微扰态。每一个算符排列就代表一个微扰过程,称为过程链(process chain)[14],而加藤项就是这些相应的过程链的和。

过程链算法的基本思路是:将相关过程数字化以方便我们利用计算机进行编程,然后考虑系统的对称性对跃迁图进行简化,保留剩余的跃迁图及其相应的权重。对于一个确定的跃迁图,需要考虑其包含的所有算符的排列,不同算符的排列对应着不同的加藤项,将算符的排列与加藤项一一对应然后对所有的排列求和就得到了该跃迁图相应的能量修正。

表1 四阶加藤项的生成和简化

我们可以如下进行过程的数字化:用0,1,2,3四个数字来标记右,上,左,下四个跃迁过程。任意选择跃迁的起始点表示外源的流入,因为系统最终要恢复到初始状态,所以跃迁箭头的终点即为外源的流出点,这样我们只需考虑中间的跃迁过程并将每一步跃迁过程用相应的数字标出,即完成了对一个过程的数字化描述。(注意,这里的数字序列与之前所讲的加藤项的数字表示毫不相干,请注意不要造成混淆。)例如,如图3所示的过程就可由(03120)来表示。观察该数字序列不难发现,为了生成所有的数字序列,我们可以按照数值从小到大,一位位地向上递增的方式得到所有可能的跃迁图,然而不同的数字序列可能会表示相同的跃迁图(这里是指不考虑算符的排列,单单从形式来看)。对于图3,我们既可以用(03120)来表示,也可以用(02031)来表示,而我们最终只需要其中的一个,不然将会造成重复计算,所以当如上得到所有跃迁图的数字序列之后,还需对这些表示做筛选处理,保留每一个图对应的其最小表示。

Fig.3 Different number order represent the same hopping graph图3 不同数字序列(03120或02031)代表相同的跃迁图

何谓最小表示[14]?仍以图3为例,可以看到从初始点出发,只有一种选择向右(标记为0),然后就要面临两种选择:向上(标记为3)或者向左(标记为2),人为地,我们选取数值小的路径,向左移动回到初始点(标记为2),然后向右移动(标记为0),之后向上(标记为3),最后向下(标记为1)。这样就得到了这个图所谓的最小表示(02031),而(03120)这一项将会被扔掉。

然而,即使去掉了所有重复的跃迁图,留下的最小表示的跃迁图的数目依然十分庞大,计算的效率依然很低,需要我们进一步从这些跃迁图中找出对能量修正贡献相同的部分进行合并,保留其中的一项并统计相同图的数目作为该图的权重予以保留。首先我们考虑跃迁图的旋转对称性,从任一跃迁图出发,利用点群的D4对称性对其进行操作,将操作之后的结果做最小表示后与剩下的跃迁图进行比对,若是发现相同的图,则在剩余的结果中去掉这一图,原图的单图重复次数增加一,在完成所有的变换后,根据群论中的拉格朗日定理,得到我们权重的计算公式:权重=变换次数/单图重复次数,可以证明这一权重必为整数。

进一步的研究中我们发现,还存在着这样的一些跃迁图:它们的数字序列并不相同,旋转对称也不重复,但是却对能量修正有着同样的贡献,如图4所示。

Fig.4 Topological invariance of hopping graph括号表示格点坐标,数字代表拓扑序列,两个跃迁的拓扑序列均为12345图4 跃迁图的拓扑不变性

将这一情况定义为跃迁图的拓扑不变性,我们可以根据这一特性对跃迁图的数目做进一步的化简。需要说明的是,这一部分所有的简化操作都是为了提高我们的计算效率,并不要求最为精简的算法,也不会对我们的能量修正的值造成任何影响。考虑到拓扑操作的特殊性,我们对跃迁过程中经过的格点按从小到大的数字进行标记,重复的格点不赋予新的标记值。我们可以看到这两个对能量修正具有相同贡献的跃迁图拥有完全相同的拓扑序列。

通过以上操作,我们就能够挑出所有不同的跃迁图像并确定出它们的权重值,接下来只需将跃迁图与加藤项相对应便可以得到数值解。需要说明的是,每一个跃迁图中都包含着若干的算符,而这些算符的作用顺序可以是任意的,这就要求我们在具体的计算中要对所有的算符进行全排列,针对每一种全排列的结果,找出与之相匹配的加藤项,然后将这些结果全部相加才是一副跃迁图完整的能量修正结果。

Fig.5 Full-size 2nd-6th order phase transition diagram of an ultracold Bose system in a square superlattice(绿色、黄色、蓝色分别表示2、4、6阶的计算结果),相界之下区域表示不同的MI态/DW态,之上为超流态图5 正方光学超晶格中超冷玻色系统的全尺寸2～6阶相变图像

Fig.6 (a)Logarithm of versus the order n in different cases;(b)Extrapolation scheme for determining the phase boundary:tc are plotted versus 1/n,and extrapolated linearly to 1/n=0图关于阶数n的线性关系曲线;(b)相界tc关于1/n的线性拟合结果,将曲线外推到1/n=0处得到无穷阶时的相界值

Fig.7 Full-size infinite order phase transition diagram of an ultracold Bose system in a square superlattice相界之下区域表示不同的MI态/DW态,之上为超流态图7 外推到无穷阶微扰极限下正方光学超晶格中超冷玻色系统的全域相图

4 结论

文章通过分析正方光学超晶格中超冷玻色子体系的非微扰哈密顿量基态分布图像,很直观地得到了不同参数区间下玻色子体系所处的相,该体系在μ和Δμ的不同区域分别处于莫特绝缘态或密度波态。在此基础上,利用推广的有效势方法,得到了体系从莫特绝缘态/密度波态到超流态的相界曲线方程,并利用瑞利-薛定谔微扰论的加藤表示,借助于过程链的计算方法,得到了体系超出平均场的高阶微扰下的全域相图。基于所得到的高阶微扰结果,利用线性拟合外推技术,我们进一步得到了无穷阶微扰下的系统全域相图。虽然目前还没有相关的蒙特卡洛数值结果可以对比,但从过往对其他系统研究的结果以及上述的整个计算过程看,我们有理由相信,我们的结果有相当的精确度，可为后续的对该系统的蒙特卡洛数值研究及相关实验提供参考。