解晓洁 孙俊松 秦吉红† 郭怀明

1)(北京科技大学物理系,北京 100083)

2)(北京科技大学,理论物理研究所,北京 100083)

3)(北京航空航天大学物理系,北京 100191)

1 引言

石墨烯具有较高的韧性和可弯曲性[1],作为一种独特的二维材料,其电子性能与机械形变密切相关.因此,可以采用对石墨烯施加应变的方法来研究它的电子性能[2].其中,应变工程在理论和实验领域都有广泛的研究.在实验中,研究人员通过在Pt(111)表面产生石墨烯纳米气泡,成功地测量到应变诱导下磁场强度超过300 T 的赝磁场[3].这一实验证明为探索超高赝磁场作用下二维材料的物理性质提供了新的基础.理论研究中,为了更好地研究石墨烯的特性,人们发展了几种应变方法: 弯曲石墨烯[4–6]、单轴应变[7,8]、三轴应变[9,10]和扭曲双层石墨烯[11].研究表明这些应变方法已成为调制石墨烯和其他二维电子材料相关量子特性的有力工具[12–14].

最近,应变诱导规范场被推广到非晶半导体[15]、蜂窝反铁磁体[16–18]和狄拉克半金属中[19].虽然蜂窝反铁磁体的朗道量子化是由应变诱导的非零赝磁场形成的,它的性质与石墨烯中的性质并不完全相同.具体来说,在蜂窝反铁磁体三轴应变诱导的非零赝磁场下,磁子赝朗道能级出现在能谱的上端并且是等间距的[18].通过类比石墨烯中的跳跃振幅并以与前者相同的方式修改交换耦合,将应变引入蜂窝反铁磁体中.应变诱导的规范场可能不像在石墨烯中那样作用于磁子.因此,用于设计石墨烯中均匀赝磁场的不同方法是否对蜂窝反铁磁体的自旋激励有相同的影响,是值得探讨的问题.

本文利用线性自旋波理论(LSWT)和量子蒙特卡罗模拟(QMC)研究了六角晶格反铁磁体在弯曲应变作用下的物理性质.首先利用线性自旋波理论分析了反铁磁体中赝朗道能级的形成,并观察到这些能级与石墨烯中所见能级的不同.然后研究了反铁磁体中磁序的演化.本文结构如下: 第2 节介绍研究的精确模型和采用的计算方法; 第3 节展示线性自旋波理论中的磁子赝朗道能级; 第4 节利用线性自旋波理论和量子蒙特卡罗模拟研究反铁磁序的演化; 第5 节给出了进一步的讨论和结论.

2 模型和方法

考虑具有自旋S的蜂窝量子反铁磁体,其哈密顿量可以表示为

其中〈i,j〉是最近邻配对;Jij是连接i,j位点的最近邻交换耦合,未加应变时格点间的交换作用一致Jij=J表现反铁磁性(J＞0);表示在位点i上自旋为1/2 的算符,位点i,j上的自旋算符满足对易关系(x,y,z循环).在没有应变的情况下,晶格格点i的位置矢量是ri,如图1,每个格点周围有3 个最近邻格点,格点i的3 个最近邻矢量δk(k=1,2,3)的值是:

图1 (a)弯曲应变下具有锯齿状边缘的内半径为R,外半径为R+W 的蜂窝扇形纳米带,纳米带在x 方向具有周期性,y 方向的宽度是Ly=6;(b)弯曲应变下宽度为Ly=200 的纳米带上边界处的晶格放大图.不同厚度的键表示不同的交换耦合值,数值标记在旁边Fig.1.(a) A honeycomb fan-shaped nanoribbon with jagged edges with an inner radius R and an outer radius R+W under bending strain.The nanoribbon is periodic in the x direction,and the width in the y direction is Ly=6.(b) Enlarged lattice at the upper boundary of a nanoribbon with a width of Ly=200 under bending strain.Keys of different thicknesses represent different exchange coupling values,and numerical values are marked next to them.

加入弹性应变后,晶格位置发生变化,从而影响最近邻交换耦合,使其偏离无应变下的值.形变后的晶格格点i位置矢量U(ri)=ri+u(ri)=((1/c+y)sin(cx),(1/c+y)cos(cx)−1/c),其中c=1/R表示应变强度,u(ri) 是平缓变化的位移函数,弯曲应变下u(x,y)=c(xy,−x2/2).新的最近邻矢量=U(ri+δk)−U(ri),连续极限下的键长变化量可写成=(δk·∂)U(ri);应变后的键长为.最近邻交换耦合以键长的指数函数形式变化为

其中β是格林艾森常数.这里考虑足够窄的纳米带上的较小弯曲,交换耦合可以近似表示为

当u(r) 是平缓的位移函数时,可以定义应变张量,应变张量产生赝矢势:

在弯曲应变下,εxx=cy,εyy=εxy=εyx=0;因此在x方向的赝矢势Ax=γcy/2.这个赝矢势产生一个均匀的赝磁场B=∇×A=.石墨烯中的狄拉克费米子在赝磁场下,形成相应的赝朗道能级.同样地,磁子也会受到赝磁场的影响,形成赝朗道能级.

考虑弯曲应变下的六角晶格纳米带,如图1(a)所示.扇形纳米带由A,B两套格子组成,每条zigzag 链上的A,B格子数量是N,格子总数是NA+NB=N2;纳米带宽度为W=,长度为.选取扇形纳米带底部中心为坐标原点,沿着y方向第j条zigzag 链上A格子格点的纵坐标为,随着y值的增加,交换耦合逐渐减小.对于有固定尺寸的弯曲蜂窝纳米带,能够研究的最大应变程度由相互作用的零值决定.上边界最高处的水平键预估会在最大应变强度下消失,因此有;纳米带宽度为W,得到最大应变强度.在最大应变强度下,赝朗道能级呈现最明显的特征.因此在后续的研究中,应变强度将取最大值.

下面的讨论中,在上述弯曲应变下使用线性自旋波理论和具有定向环更新的随机级数展开(SSE) 量子蒙特卡罗方法[20–24]研究等式(1)中的模型.SSE 方法将配分函数展开为幂级数,并且将求迹表示为对角矩阵元素的总和,定向循环更新使得模拟进行得很有效.纳米带沿x方向具有周期性,没有近似造成的系统误差,因此可以在一个离散构型空间中进行抽样来分析物理性质.在进行量子蒙特卡罗模拟时,设定温度的倒数为βt=200,该温度足够低可以获得基态性质.

3 线性自旋波理论

首先使用线性自旋波理论来研究应变模型(1)式的物理性质,通过Holstein-Primakoff (HP)变换[25]用玻色子算符代替自旋算符,在子晶格A,B上自旋算符写为

在线性自旋波近似下,自旋算符可以简化为

将简化后的自旋算符代入到哈密顿量中,保留双算符项,得到玻色子的紧束缚哈密顿量

在x方向做傅里叶变换

T矩阵有下列形式:

其中A+是NA×NA的矩阵;B−是NB×NB的矩阵.因为基矢包含玻色子算符,变换矩阵T保持玻色子的对易关系,T矩阵满足

因此,原基矢可以用对角化表象下的新基矢(αi,βi)表示

在线性自旋波理论下,可以计算得到局域磁化强度

晶格形变产生赝磁场并且诱导赝朗道能级,磁子能带中出现一系列平台.因此,在磁子态密度中出现一系列尖锐的峰是赝朗道能级形成的直接证据.利用线性自旋波理论获得的磁子能谱,可以计算磁子态密度

其中,δ函数近似为一个具有很小常数c0的窄高斯波包,在下面的计算中取展宽c0=3×10−3.在多体物理中,局域磁化率是一个与实验更相关的物理量,其定义为

在线性自旋波理论中,可以得到局域磁化率的表达式:

首先研究无应力的情况,图2(a)计算了磁子能谱,中间出现的独立分支与边界相关.我们发现对应的波函数主要分布在边界附近,说明了这些模式是边界态.磁子态密度如图2(c)所示,反铁磁体中的磁子激发呈现出线性色散的低能态.这种线性色散在系统的态密度中得到体现.在施加应变后,能级简并度降低,能谱变得更为扩展,如图2(b)所示.这种变化在动量kx=π 处最为明显,在没有应变的情况下,kx=π 处所有的能级解是简并的.高能范围的磁子能谱变得更平缓,同一能量上的态数目增多,态密度表现出明显的振荡行为.在高能附近,磁子态密度出现尖锐的尖峰,这是磁子能带中的一系列平台造成的.因此,态密度中尖锐峰的出现是磁子赝朗道能级出现的直接证据.

图2 锯齿状反铁磁蜂窝纳米带的磁子能谱(a)无应变情况;(b)弯曲应变情况;图(c)和图(d)是图(a)和图(b)两种情况对应的磁子态密度.在图(b)和图(d)中,应变强度为 c/cmax=1.图(d)中的红色虚线框内的曲线将放大显示在图3 中Fig.2.The magnon spectrum of the antiferromagnetic honeycomb nanoribbon: (a) Without strain; (b) under bending strain.Panels (c) and (d) are the magnon density of states corresponding to panels (a) and (b).In panels (b) and (d),the strain strength is c/cmax=1.The curve inside the red dashed box in panel (d) will be enlarged as shown in Fig.3.

图3(a)给出了图2(d)高能端(红色虚线框范围)的放大图,当施加弯曲应变时,磁子的赝朗道能级从能谱的顶端开始出现,这与石墨烯中赝朗道能级在低能端出现的情况有明显不同[26–28].然而,如图3(b)所示,通过拟合峰的位置,可以发现能量ωn与能级指数n的平方根成正比,这表明它们之间呈线性关系,这类似于石墨烯中狄拉克费米子的性质.磁子赝朗道能级的出现还可以在局域磁化率上反映出来.图3(c)和图3(d)给出了子晶格A和B上的局域磁化率,它们也形成了尖锐的峰,且峰的位置与态密度中峰的位置完全相同,进一步证实了弯曲应变下赝朗道能级的存在.

图3 (a)磁子能谱上端态密度;(b)赝朗道能级的能量ωn 作为能级指数n 的平方根的函数,实线反映了数据的线性拟合.图(c)和图(d)是子晶格A 和B 上的局域磁化率.在图(a),图(c)和图(d)中,应变强度采用c/cmax=1,线性尺寸Ly=200,Lx=20Fig.3.(a) The density of states at the upper end of the magnon spectrum; (b) the energy ωn of the pseudo-Landau level is a function of the square root of the energy level index n,and the solid line reflects the linear fit of the data.Panels (c) and (d) are the local magnetic susceptibility on the sublattices A and B.In panels (a),(c) and (d),the strain strength is c/cmax=1,the linear size is Ly=200,Lx=20.

4 反铁磁序的演化

海森伯自旋模型是一个典型的强关联模型,在二维空间以上并没有严格的解析解.在第3 节的线性自旋波理论中,将自旋算符表示成玻色子算符的(6)式和(7)式是精确的.但是得到的玻色子哈密顿量将包含四算符以上的相互作用项,解析上仍然无法处理,因此通常做线性自旋波近似(见(8)式和(9)式).该近似的物理意义是忽略自旋量子涨落,因此只能得到海森伯模型的定性物理性质.为了得到海森伯模型的精确解,需要采用完全考虑量子涨落的数值方法,比如本文采用的随机级数展开的量子蒙特卡罗数值方法(该算法详细细节见附录A).在下面的研究中温度倒数设置为βt=200,该温度足够低可以代表有限尺寸系统的基态性质.

进一步,研究了应变对反铁磁序的影响.在线性自旋波理论中,通过有限的磁化强度来确认反铁磁序的存在.由于弯曲的蜂窝纳米带在x方向上具有平移不变性,局域磁化只在扩展到整个带宽的单元格内变化(见图1).图4 给出了几种应变强度值下,单元格中局域磁化强度与格点位置之间的关系.在没有应变的情况下,边界附近的反铁磁序会受到扰动.当格点接近边界时,数值逐渐减小.然而,不管位置如何,局域磁化强度ms(i) 始终是有限的,且当远离边界时,大部分格点的数值趋于均匀.这意味着在存在开放边界的情况下,长程反铁磁序仍然能够保持.需要注意的是,边界最外层格点的局域磁化要比其附近的格点更大(详见图4(b)和图4(c)).这两个最外层格点代表了边界锯齿链的两个子晶格.尽管沿着一维链方向上是反铁磁的,但是这两个格点的磁矩不相等,这在zigzag 边界上形成了一个净铁磁矩.因此,在锯齿边界上形成了一个铁磁序,这在蜂窝纳米带上的哈伯德模型中也得到了揭示.

图4 (a)在应变强度c/cmax=0,0.5 和1 下,通过线性自旋波理论和量子蒙特卡罗模拟得到的局域磁化强度.图(b)和图(c)分别放大了图(a)中无应变和最大应变时靠近下边界的曲线.计算的格点范围延伸到400 (与此后的图表相同)Fig.4.(a) The local magnetization obtained by linear spin wave theory and quantum Monte Carlo simulations under c/cmax=0,0.5 and 1.Panels (b) and (c) magnify the curves near the lower boundary at no strain and maximum strain in panel (a),respectively.The calculated grid point range is extended to 400 (same as charts hereafter).

在施加应变后,局域磁化的数值随着应变强度的增加而逐渐减小.这是因为在反铁磁材料中,应变会影响交换耦合强度,而交换耦合决定了磁性的稳定性.在y方向上,交换耦合逐渐减小,因此离下边界较远的格点受到的影响更大.特别是在足够大的应变强度和靠近上边界时,局域磁化迅速减小,并在临界位置处变为负值,这意味着反铁磁序在此后消失.这是因为靠近上边界的位置,水平键的交换耦合非常小,系统可以被视为一系列孤立的一维垂直键.在此情况下,二维反铁磁序消失[29–32].相反地,靠近下边界的区域受到的影响较小.因为这里的交换耦合程度改变最小.因此,在该区域内,反铁磁序相对稳定.尽管线性自旋波理论忽略了自旋量子涨落,该量子涨落能量尺度正比于1/S,因此线性自旋波方法只有在自旋S取很大数值时才是精确的.本文研究的是自旋1/2 的海森伯模型,因此自旋波理论进行的近似较大.为了得到反铁磁序随应变演化的精确结果,需要使用大尺度的精确量子蒙特卡罗数值模拟[33,34]方法进行计算,该方法可以在计算机数值精度内准确地模拟海森伯模型基态的物理性质.

其中,通过对所有晶格格点j进行求和来计算局域磁化的数值;符号因子sgn(i,j)=1(−1)由i和j属于相同(相反)的子晶格来确定.图4 对比了相同应变强度下线性自旋波理论与量子蒙特卡罗模拟中的数值.可以看到,量子蒙特卡罗模拟得到的数值小于线性自旋波理论的预测值.此外,随着应变强度的增加,两种方法得到的数值差异也会增大.当没有应变时,量子蒙特卡罗模拟曲线缓慢上升,并在大部分区域趋于平缓.需要注意的是,在图4(a)中,量子蒙特卡罗模拟和线性自旋波理论的结果只有轻微不同,大部分数值之间的差异小于10%.这意味着HP 变换中使用的线性近似是相当准确的,这也得到了现有文献的证实.在取最大应变强度的0.5 倍的情况下,量子蒙特卡罗模拟和线性自旋波理论在上边界附近的结果明显不同.根据线性自旋波理论,局域磁化强度在靠近上边界时会迅速减小,然后趋于稳定; 而在更精确的量子蒙特卡罗模拟中,局域磁化在靠近上边界时出现了上翘的现象.因此,需要进一步讨论,以研究导致这种变化的原因.

如图5 所示,利用量子蒙特卡罗模拟分析了应变强度 0—1 区间的局域磁化强度,可以发现c/cmax为 0.4—0.7 时,上边界附近的曲线向上弯曲,局域磁化强度增大,当应变强度为最大应变强度的0.5 倍时,变化最为显著.

图5 在应变强度 c/cmax=0—1 下,量子蒙特卡罗模拟得到的局域磁化强度(a)0—400号格点上的数值;(b)250—400号格点上的数值Fig.5.The local magnetization obtained by quantum Monte Carlo simulations under c/cmax=0–1: (a) Values at grid points 0–400; (b) values at grid points 250–400.

接着研究了交换耦合随晶格位置的变化.J⊥表示竖直方向的交换耦合,J//表示沿zigzag 链方向的交换耦合;令J⊥=1,.图6表明,垂直交换耦合保持不变,水平交换耦合随着y值的增加逐渐减小,在靠近上边界时趋近于零.这些结果表明交换耦合几乎不影响上边界处局域磁化的上翘.

图6 在应变强度 c/cmax=0—1 下,作为zigzag 链指数函数的垂直和水平交换耦合Fig.6.Under c/cmax=0–1,the variation of horizontal bond and vertical bond with the lattice position.

最后计算了自旋关联强度,其定义为C(i,j)=.本研究中使用的模型由200 条锯齿状链组成,一个单元内有400 个格点.由于上翘只出现在系统上边界处,所以这里研究靠近上边界格点处的关联,如图7 所示,关注390—400 号格点,计算他们的水平关联平均值和垂直关联平均值.具体来说,水平关联平均值对应于zigzag 链上的第10 个格点的值.取390—-400 号中某一格点,计算该格点与380 到400 号格点的垂直关联的平均值,作为垂直关联平均值.结果表明,水平关联在上边界附近逐渐增加,而垂直关联基本不变.这些结果表明,局域磁化在上边界处的上翘在很大程度上受水平关联的影响.

图7 在c/cmax=0—1 范围内390—400 号格点上的(a)水平关联,(b) 垂直关联Fig.7.(a) Horizontal correlations,(b) vertical correlations at grid points 390–400 under c/cmax=0–1.

5 结论

本文利用线性自旋波理论和量子蒙特卡罗模拟研究了反铁磁蜂窝纳米带在弯曲应变下的赝朗道能级和反铁磁序的演化.在施加应变后,在能谱的上端形成了一系列平台.这些能级的间距与能级指数的平方根成正比.水平键的交换耦合强度在y方向上逐渐减小,在足够大的应变强度下,在纳米带上边界附近的水平方向的交换耦合减弱为零,系统解耦成孤立的垂直棒.两种计算方法都发现,局域磁化强度随y值的增加而逐渐减小,并且在靠近上边界处消失.这说明上边界附近的反铁磁序被弯曲应变破坏,虽然线性自旋波理论给出了反铁磁序演化的正确趋势,但它在相同应变强度下绘出的反铁磁区域更大,相比而言,量子蒙特卡罗结果更加准确.在过去的几年里,二维量子磁性材料领域取得了重大进展[36–39].应变工程作为一种控制和操纵磁性物态的重要方法,在二维磁性材料的研究中会有广泛的应用.本文研究的弯曲应变导致的新的物理现象一定会引起相关实验的兴趣.

目前人们已经发现一系列二维磁性材料,典型的材料包括三卤化铬 CrX3(X=Cl,Br,I)[40].在对这些二维磁性材料的研究中应变已经开始作为调控二维磁性材料磁性的一种有效方法.利用密度泛函理论研究二维单层材料CrI3的应变演化过程中,研究者发现应变可以诱导铁磁与反铁磁的磁态转变[41].与石墨烯类似,可以通过将单层二维材料放在衬底上,调控衬底来实现磁子赝朗道能级的应变.除了六角晶格反铁磁体,笼目晶格中非线性反铁磁中的应力效应也值得进一步研究,具体包括:非线性反铁磁序随应变演化,利用应变调控笼目晶格上的反铁磁序,以及磁子能谱中可能的赝朗道能级.最近的研究发现,具有层状笼目结构的Mn3Sn 材料中的非线性反铁磁可以由电流控制,这为研究笼目晶格的应力效应提供了额外的调控手段[42].

附录A

量子蒙特卡罗随机级数展开(SSE)方法就是通过对玻尔兹曼算子 e−βH做泰勒展开来构造适合蒙特卡罗采样的构型空间:

泰勒展开后来被更普遍地认为是广泛模型的精确量子蒙特卡罗 (QMC)算法的起点.指数算子的幂级数展开对于有限β和有限晶格收敛.选择一套基矢,配分函数可以写为

其中 {α}n对系统的所有基矢求和.以自旋为S=1/2 的海森伯模型为例,可以将哈密顿量写成对键算符求和的形式,键算符定义为

进一步可以将每个键对应的相互作用算符分成对角算符和非对角算符Ha,b:

其中a=1,2 分别代表对角和非对角算符,b代表键算符作用的键的位置.因此可以将原始的海森伯哈密顿量写成

将定义在键上的哈密顿量的形式代入到配分函数中,于是有

这个配分函数的形式就是SSE 算法的起点.其中Sn是哈密顿量中各个键上的对角算符和非对角算符乘积组成的算符序列,来源于Hn.当Sn这样一个算符序列作用在某一个基矢态后,会把这个基矢态变成其他的基矢态.根据定义的对角算符H1,b和非对角算符H2,b的形式可以知道,只有当这两种算符作用的键连接的两个自旋方向相反时才会对配分函数有贡献,因为Ha,b作用后非0 的只有4 种情况:

所以对于一个包含n个Ha,b的算符序列Sn,只有当这个算符序列作用在某个基矢态 |α〉后,能够回到原来的基矢的状态,这个算符序列作用在相应的基矢态时才会对配分函数有贡献,即

由于对配分函数有贡献的构型中每一个非单位算符会贡献一个1/2,因此每一个允许的构型对应的权重为

这样定义出每个构型对应的权重之后就可以按照经典蒙特卡罗的方法对符合条件的构型和算符序列按照它们对应的权重进行抽样.