第一性原理计算Co/h-BN界面上的微弱Dzyaloshinsky-M oriya相互作用∗

2018-06-19黄灿李小影朱岩潘燕飞樊济宇施大宁马春兰

物理学报 2018年11期

黄灿李小影朱岩† 潘燕飞‡ 樊济宇施大宁马春兰

1 引言

1962年,Skyrme[1]研究重子物理,命名了斯格明子(Skyrmion)这种受拓扑保护的粒子.斯格明子可以用拓扑数来表征而且不会被连续变化的外场所改变[1].1975年,Belavin和Polyakov[2]提出了一种自旋结构与斯格明子一致的二维海森伯自旋系统中的亚稳态,命名为磁性斯格明子.2006年,Bogdanov等[3]预言手征磁性导体MnSi中可能存在斯格明子晶体相.直至2009年,Pfl ciderer研究小组[4]在实验中观察到MnSi中含有磁性斯格明子.斯格明子是一种既具有拓扑性也具有手征性的纳米尺寸磁结构[3,5].

磁性斯格明子的产生机理是所谓的Dzyaloshinsky-Moriya相互作用(DMI)[6−10].1958年,Dzyaloshinsky[6]首先由晶体对称性提出DMI作用.1960年,Moriya[7]进一步地阐述了DMI,并且用来解释了铜氧超导体中的弱铁磁性.DM是这两位科学家名字的缩写,而DMI本质就是自旋轨道耦合和空间反演对称性破缺共同导致的磁性相互作用.根据DMI方向的不同,斯格明子可以分为布洛赫型(Bloch-type)斯格明子或奈尔型(Néeltype)斯格明子.由于DMI相对海森伯交换作用(Heisenberg interaction,HBI)较弱,只有在低温下才能有效果,而低温时较弱的热激发又会阻碍斯格明子的产生.

二维界面具有垂直方向的非中心对称,可能会产生稳定的奈尔型斯格明子[11],且二维材料与半导体技术发展能够更好地兼容.磁性和非磁性重金属之间的界面薄磁膜(例如Co/Pt,Co/Ir等)利用了界面来实现对称性破缺和重金属有效的自旋轨道耦合,引起了广泛的关注和研究.在技术上,已经利用自旋极化隧道电流手段在PdFe/Ir(111)薄膜上实现了单个斯格明子的写入和擦除[12].基于斯格明子的自旋电子器件,具有高稳定性、高存储密度、低电流驱动等特点.

越来越多的二维材料被用来研究产生斯格明子.其中,石墨烯(graphene)由于其独特的性质,包括具有单原子层厚度、无质量线性色散的电子结构以及长程自旋扩散[13−17],被广泛用来研究磁性相关的性质.石墨烯本身自旋轨道耦合比较弱,一般来说不会引起足够的DMI来影响手征磁性.但是在2017年,Yang等[18]通过自旋极化电子显微镜的实验研究和第一性原理计算表明,这种石墨烯诱导的DMI可以具有与重金属界面相同的量级:在石墨烯和Co(Co/Gra)的界面上每个原子的DMI可以达到1.14 meV.这项工作可以认为在Co/Gra界面上诱导的奈尔型斯格明子是稳定的[18],为斯格明子的研究开拓了新的领域.六方氮化硼(h-BN)薄膜,是一种类石墨烯蜂巢晶格结构薄膜,由交替的硼原子和氮原子组成.二维结构的h-BN薄膜,以其良好的电绝缘性(能带带隙近似6.0 eV)、低摩擦系数、高热稳定性和化学稳定性、可见光区内有良好的透过性和深紫外光区的光吸收性以及其与石墨烯极为接近的晶格结构(误差仅1.7%)吸引了广泛关注[19−22],在力学、电学、热学和光学等方面具有较为广泛的应用.

h-BN可以通过过渡金属掺杂成为良好的稀磁半导体,形成磁性半导体界面. 其中Co在h-BN上(Co/h-BN)也是一种研究得比较广泛的界面[23−25].本文通过第一性原理计算Co/h-BN界面DMI的大小.DMI的第一性原理计算刚刚开始,国内尚未有课题组开展类似工作.目前传播最广泛的第一性原理计算软件为VASP(Vienna ab-initio simulation package),国内外尚未见通过界面自旋螺旋各向异性的方法来计算磁性原子间各近邻DMI的相关报道.本课题组从VASP软件中挖掘出一套计算DMI的方法,本文采用此方法计算了Co/h-BN的各近邻间HBI和DMI参数,为以后各界面的DMI计算打下了坚实的基础.

2 计算方法与理论模型

采用基于密度泛函理论第一性原理赝势平面波方法[26]的VASP软件进行计算.VASP是利用赝势方法和平面波基组进行量子力学分子动力学从头计算的软件包.在VASP计算中,采用了VASP自带的Perdew-Burke-Ernzerhof势[27].VASP采用平面波基矢通过自洽迭代方式来求解Kohn-Sham方程,可以很容易计算出力与张量.计算体系中的原子再通过受力弛豫到基态上,力的收敛条件是10−1eV/nm.平面波截断能选取为500 eV,初始截断能设定为400 eV用来延展平面波函数,电子步自洽计算收敛条件为1.0×10−6eV.倒格空间的k点分布设置为25×25×1.采用平板周期模型(slab model)来模拟Co/h-BN薄膜:Co原子吸附在单层h-BN上.h-BN采用2×2超胞,层间距设置为1.5 nm.

过去绝大部分的VASP计算中,磁性原子间的磁交换能只考虑了最近邻原子间的作用J1,并没有得到其他更远近邻间HBI磁交换能Ji[28−31].2002年,Marsman和Hafner[32]首次通过广义布洛赫条件在VASP中计算了γ-Fe自旋螺旋(spin spiral)色散关系,得到了各个近邻的Ji,并用来解释γ-Fe磁性结构.本文尝试通过非线性计算[33],考虑自旋轨道耦合,利用广义布洛赫条件[34,35],计算自旋螺旋的色散关系,并通过色散关系的拟合得到各个近邻的HBI参数Ji以及DMI参数di.

具有DMI的海森伯模型系统的哈密顿量可以表示为

式中第一项为HBI,设置铁磁态为基态且能量为0,磁矩旋转导致的HBI能量变化是与铁磁态相比较的相对值;第二项为DMI,si(j)为归一化的磁矩.dij=dij(z×uij),其中uij为连接两个磁性原子i与j的基矢,z垂直于界面,方向从BN层指向Co层,dij为各近邻间DMI的大小.如果设置磁矩在垂直于二维界面的平面内变化,形成自旋螺旋,且dij的值不为零,那么DMI就打破了自旋螺旋的简并性,导致自旋螺旋色散关系的空间不对称性.本文通过VASP计算反方向的自旋螺旋色散关系E(q)与E(−q),其中q为自旋螺旋的方向矢量.则∆DMI(q)=E(q)−E(−q)为所计算体系2倍DMI能量随q的变化.一般而言,同一波矢q下,若DMI的绝对值越大,则dij的值越大,斯格明子越小,更便于实际应用.

3 结果与讨论

3.1 原子结构与电子结构性质

h-BN的晶格常数计算值为0.251 nm,与实验结果0.251 nm符合[36].该值接近局域密度近似(LDA)的计算结果0.250 nm[21],0.249 nm[19],而与广义梯度近似(GGA)计算结果一样[22].相应地,B—N键长为0.145 nm.Co原子在h-BN表面有4种吸附方式.如图1(a)所示,第一种是Co原子处在h-BN六角环结构的中间位置(hollow),简称h位;第二种是Co原子处在BN共价键的中间(bridge),简称b位;第三种是Co原子处在N原子的顶(top)位,简称tN;相应的第四种是tB位.从图2中可以发现,Co吸附位置从tN,tB,b到h,Co/h-BN体系的能量逐渐降低,即h位的结构最稳定,与之前的研究符合[21,22].本文选择h位结构进行计算.在2×2超胞中,吸附的Co原子磁矩朝上;在2×4结构中,吸附的两个Co原子一个磁矩安排朝上,另一个朝下,呈反铁磁结构,分别进行弛豫计算.铁磁与反铁磁结构的Co与h-BN层的垂直距离约分别为0.192 nm与0.191 nm,两者几乎一样.2×2超胞中,由于Co浓度高,Co原子间还保留有相互作用(在2×4结构中由于吸附了两个Co原子,所以浓度与在2×2超胞吸附一个Co原子一样),这导致垂直距离略大于4×4超胞中Co原子到h-BN层的垂直距离0.185 nm[22].在铁磁、反铁磁结构中,CoN与CoB键长分别为0.241,0.239 nm,也因此同分别样略长于4×4超胞中的0.237,0.235 nm.Co原子对BN键长影响不大,变化范围在0.001 nm之内.

图1 (a)Co/h-BN界面2×2超胞原子结构俯视图,蓝色大球、银色中球、绿色中球分别表示Co,N和B原子;α表示靠近Co的N或B原子,β表示远离Co的N或B原子;h,b和t分别表示Co吸附在h-BN的空位、桥位和顶位;Co原子周围的黄色曲面是自旋密度等高面;(b)结构图是将2×2原胞扩展成5×5倍后的俯视图,用来描述Co原子间各近邻关系,其中以中间的Co原子为中心原子,标记为0,其他Co原子与中心原子的近邻数也被标记在Co原子上;a 1和a 2为基矢,b1和b2为相应原胞基矢的倒格矢:Γ,M 和K 是第一个布里渊区的高对称k点;(c)平行于h-BN面并穿过Co原子中心的等自旋电荷密度图,(d)贯穿h-BN面中间的等自旋电荷密度图(c)与(d)等高线的值为ρn=0.001×n e/a30,其中n为图中等高线的标记值,a0为玻尔半径,e为电子电量;(e)反铁磁结构的等自旋电荷密度线图,图中标记的两个Co原子的磁矩排列相反,数字标注与(d)相同,粉色线是负等高线Fig.1.(a)Top view of Co/h-BN interface with 2×2 unit cell.Diff erent atoms is labeled:Big blue,middle silvery and Middle green ball is for Co,N and B atoms respectively.αandβpresents the sites of B or N atom close to or apart from Co atom respectively.The labeled h,b and t are for hollow,bridge,and top sites for Co atom adsorbed on h-BN sheet respectively.The yellow curved surface around Co atom is spin density;(b)extended 5×5 supercell of 2×2 unit cell is used to describe diff erent neighbors in Co/h-BN,as shown the yellow numbers on the atom.The labeled a 1 and a 2 is basis vector and b1 and b2 is reciprocal lattice vector;Γ,M and K are high symmetry k-points in fi rst Brillouin zone;(c)and(d)is the spin charge density contours plotted across Co atom and h-BN layer respectively,contour values are ρn=0.001×ne/a30,where n is the mark value of the contour,a0 is the Bohr radius,e is the charge of the electron;(e)is the spin charge density for antiferromagnetic structures,where magnetic moments are arranged anti-parallel to each other.Contour values are the same to(d),pink lines are for minus contour values.

图2 Co/h-BN中Co原子在不同吸附位置的能量Fig.2.Energy of Co/h-BN interface as a function of diff erent adsorption sites of Co.

图3 图1(a)中Co/h-BN各原子分态DOS图,黑色实线表示上自旋的DOS,红色实线表示下自旋DOS,能量为0处的紫色虚线代表费米面Fig.3.Partial DOS of the Co atom and selected N or B atom in Fig.1(a),zero energy indicates the position of the Fermi level shown as purple dashed line,positive and negative values of partial DOSare for the spin-up and spin-down channels,respectively.

Co原子的吸附给整个体系引入了磁矩,如图1(a)和图1(c)所示,整个体系的磁矩主要集中在Co原子上.在2×2超胞中,单Co原子吸附浓度达到了每原胞25%,Co原子间有较强的相互作用,整个体系的磁矩为2.77µB,稍小于4×4超胞中的3.00µB[22],其中µB为玻尔磁距.从图1(d)可以看到,主要是Co原子下方h-BN六角环内的N原子,即N_α受到了影响.N_α周围聚集着由Co诱导的自旋电荷密度.而其他原子旁边的自旋电荷密度非常小.自旋电荷密度在各原子中的分布可以从分态态密度图(density of states,DOS)中更清晰地看出来.如图3所示,Co原子的自旋分态DOS主要集中在−3.2—0 eV,而B与N的DOS主要分布在−3.5 eV以下,其中0指费米能级,这是由于两类元素的化学势不一样.但在Co原子DOS所在的能量范围内,B和N原子也还是被诱导出一定的DOS.Co上自旋通道与下自旋通道都有DOS分布,其中上自旋DOS被填满,下自旋DOS被半填满,费米能级附近也有DOS分布,整个体系显金属性.同样,从图3的绿色方框内可以看到,N_α(N_α是Co原子最近邻的N原子)的分态DOS受Co原子影响最大,而B_α(B_α是Co原子最近邻的B原子)在此范围内DOS分布很小.即Co原子磁矩间的相互作用主要是通过Co附近的N原子传递的.B_β(B_β是Co原子次近邻的B原子)的DOS要明显大于B_α,在费米能级上,N_α与B_β也具有一定的DOS.所以Co原子间的相互作用是主要通过Co-N_α-B_β网络来传递的.在反铁磁结构中,方向相反的磁矩也同样是通过相同网络来传递相互作用的,如图1(e)所示.N_β(N_β是Co原子次近邻的N原子)的磁矩则显示出一定的局域性.

3.2 自旋螺旋的计算与分析

图4 (a)离散点是Co/h-BN体系计算得到的自旋螺旋能量色散关系E(q);(b)E(q)与E(−q)之间的能量差∆DMI(q),插入部分是反方向自旋螺旋示意图,粉红色箭头表示磁矩方向,所选取的自旋螺旋中相邻磁矩间成90◦角Fig.4.(a)Scattered symbols are calculated energy dispersion E(q)of spin spirals in Co/Gra as a function of the spiral wave vector q;(b)∆DMI(q),gotten from the energy diff erences between E(q)and E(−q).In(a)and(b),lines are fi tted ones.Insert part is the spin spiral with opposite direction,where the pink arrows presents the direction of magnetic moments of Co atom and they perpendicular to each other between neighbors.

从3.1节可以得出Co原子磁矩间的作用网络,但各近邻间的HBI以及界面形成的DMI还需要通过自旋螺旋的方法计算得到.如图4(a)所示,本文通过VASP计算了自旋螺旋的能量与波矢q的色散关系E(q).自旋螺旋的波矢q选取了图1(b)所示的简约布里渊区的边界KΓ,ΓM,−K−Γ,−Γ−M.在HBI中,KΓ,ΓM与−K−Γ,−Γ−M是等价的.但DMI的作用使这种对称性发生了破缺,如图4(a)插入部分所示,当磁矩顺时针排布时,两个水平方向磁矩指向朝下即磁矩指向BN层的Co原子;当磁矩逆时针排布时,两个水平方向磁矩指向朝上即磁矩指向真空的Co原子.因此,E(q)与E(−q)的能量不完全一样.尽管在2×2超胞的Co/h-BN中,这种不同非常微小,但其能量差还是可以有规律地显示在图4(b)的离散点中.

E(q)与Co原子间的近邻数密切相关.如图1(b)所示,图中的中心原子Co设置为0,与中心原子最近邻的原子标为1,可以看到有6个最近邻原子.第二、第三近邻的原子分别标为2,3,也各有6个原子.然后将第四近邻的原子标为4,共有12个,这里并没有标全.即与一个Co原子作用的有6个最近邻,6个次近邻,6个第三近邻和12个第四近邻的Co原子.根据(1)式与h-BN六角格子的特殊结构,本文首先先推导出自旋螺旋的色散关系公式.

在三角结构中,选择一个磁矩作为零点(命名为S(0)),位置j处的磁矩位置用Rj=m a1+n a2表示,其中a1和a2为基矢.自旋螺旋的方向描述为q=q1b1+q2b2,其中b1和b2为倒格矢.在VASP计算中,所有的磁矩都设置在xz平面上,而S(0)沿着z轴方向.在广义布洛赫条件下,波矢为q的自旋波中,第j个相邻Co原子磁矩S(Rj)为

根据(1)式,只考虑到第四近邻,S(0)与第j近邻的磁矩S(Rj)之间的HBI之和分别如下:

总的HBI为

式中J1—J4为HBI参数,包含S(0)2.

DMI只考虑到第二近邻,S(0)和S(Rj)之间的DMI之和分别为:

总的DMI为

式中d1和d2为DMI参数,也包含S(0)2.

自旋螺旋E(q)的总能量包括HBI和DMI:

显然,EHBI(q)为q的偶函数,而EDMI(q)为q的奇函数.因此

图4(a)中,VASP计算的Co/h-BN自旋螺旋的色散关系E(q)与E(−q)应满足(4)式.实际上图4(a)通过最小二乘法拟合出的实线与计算值离散点符合得非常好.图4(b)所示的E(q)与E(−q)之差∆DMI(q)表示S(0)与其他磁矩之间的DMI之和的2倍,应满足(5)式.拟合出的实线与计算离散点趋势大致符合.

通过图4(a)可以看出,铁磁态Γ所对应的点能量最高.也就是说2×2超胞的Co/h-BN界面的磁有序态是非铁磁态.E(q)沿着ΓK,ΓM都快速下降,最稳定的自旋波态在K点,q=b1/3+b2/3.这种自旋波态实际上是三角晶格反铁磁[37−39],即相邻的3个Co原子的磁矩彼此成120◦角,是一种磁性几何阻挫系统.Mryasov等[34]是通过扩展超胞,排布几种磁矩构型,再计算得到各构型能量,得到体系的基态.这种方法比较消耗资源,并且只能模仿几种有限的磁矩构型,有时候并不能确定真正的基态.本文利用广义布洛赫条件,通过原胞就能得到体系自旋螺旋色散关系E(q),给出能量最低的自旋螺旋状态.在本文的体系中,最稳定的自旋波态在K点,但是M点的能量与K点相差不大,在2.8 meV内.即这种自旋波态不稳定,在稍高温度下,自旋螺旋形状就有可能发生改变,但铁磁态是最不可能出现的磁矩排布.

∆DMI(q)的趋势具有振荡特征,不是单调的正弦函数.由(3)式可知,叠加的两个正弦函数振幅应相可比拟.当然,∆DMI(q)本身就非常小,其最大值为0.1 meV,而K点E(q)的绝对值达到71.0 meV.所以在Co/h-BN结构里,DMI效应非常微弱,磁性原子间的相互作用主要是HBI.

通过最小二乘法拟合出的HBI系数J1—J4,DMI的系数d1和d2也显示在图5中.J1为−17.23 meV,其绝对值远大于J3的1.30 meV.而J2与J4均为0.02 meV,接近于0.起绝对主导作用的就是第一近邻的J1.J1为负值,这决定了Co原子间的耦合为反铁磁,这与图4(a)中E(q)在K,M点能量最低相呼应.在四方结构中,负的J1会导致G型反铁磁,即任意相邻磁矩方向排布都相反.但在六角结构中,由于最近邻原子数有6个,从对称性上来讲,不可能排布成任意两个磁矩都相反的磁有序态.所以只能排布成三角反铁磁或者M点型的反铁磁.在只有负的J1起主要作用的情况下,正如上段所叙述,三角反铁磁是最稳定的磁有序态,但是M点型的反铁磁能量也很接近三角反铁磁.

DMI的d1和d2的绝对值远小于0.1 meV,相比于J1,d1与d2更可以忽略不计.所以Co/h-BN是非DMI的界面.因此Co/Gra界面的DMI效应[18]并不是在所有二维六角环结构的界面上都会存在,在以后的研究中需要通过计算和实验一一甄别.当然,其接近于0的DMI,也有可以应用的地方.比如利用h-BN在空气中的稳定性,做一些强DMI体系的覆盖层,既可以保护原来的界面,又可以不对原来界面的DMI产生影响.

但是为什么Co/h-BN得到的DMI值与Co/Gra的DMI值相差很大呢?本文发现在2×2超胞中,Co到h-BN的距离为0.192 nm,远远大于Co/Gra中的0.156 nm.Yang等[30]指出最靠近Pt的Co原子层具有大的DMI值,次近邻的Co原子层DMI值就已经大幅下降.是否可以通过外界条件降低Co到h-BN的距离,从而调控DMI值呢?本文认为通过适当的措施(如加电场,尤其是从Co层指向BN层方向的电场),可以使得Co阳离子在电场的作用下靠近BN层,从而使得Co/h-BN表面产生DMI,这样就可以控制Co/h-BN界面在有DMI和无DMI之间来回切换,实现数据的存储.这将会是本课题组接下来进行的工作.

图5 拟合得到的HBI的J值和DMI的d值,其中J值考虑到第四近邻,d值考虑到第二近邻Fig.5.Fitted HBI J and DMI d,J is considered to the fourth neighbor and d is considered to the second neighbor.

4 结论

根据密度泛函理论的第一性原理方法计算了Co/h-BN的原子结构与电子结构,Co原子磁矩间的相互作用是主要通过Co-N_α-B_β网络来传递的.结合广义布洛赫条件,进一步计算了反方向的自旋螺旋能量与波矢的色散关系E(q)与E(−q).E(q)与E(−q)能量之差反映了Co/h-BN界面上下层之间空间对称性破缺引起的DMI的大小.HBI模型与DMI模型可以很好地拟合计算值.通过拟合得到了Co原子间各近邻的HBI参数J1—J4与DMI参数d1,d2. 在Co/h-BN中,J1为负值起完全主导作用,J3比J1小一个量级,其他参数接近于0.因此,Co/h-BN的基态是三角反铁磁,而基本上没有DMI.根据这种性质,h-BN可以作为其他DMI界面的覆盖层.Co/h-BN的DMI值与Co/Gra的DMI值相差很大的原因是,在2×2超胞中,Co到h-BN的距离为0.192 nm,远远大于Co/Gra中的0.156 nm.

[1]Skyrme T H R 1962 Nucl.Phys.31 556

[2]Belavin A A,Polyakov A M 1975 JETP Lett.22 245

[3]Rößler U K,Bogdanov A N,Pfl eiderer C 2006 Nature 442 797

[4]Mühlbauer S,Binz B,Jonietz F,Pfl eiderer C,Rosch A,Neubauer A,Georgii R,Böni P 2009 Science 323 915

[5]Yu X Z,Onose Y,Kanazawa N,Park J H,Han J H,Matsui Y,Nagaosa N,Tokura Y 2010 Nature 456 901

[6]Dzyaloshinsky I 1958 J.Phys.Chem.Solids 4 241

[7]Moriya T 1960 Phys.Rev.120 91

[8]Shu L,Chen Y G,Chen H 2002 Acta Phys.Sin.51 902(in Chinese)[殳蕾,陈宇光,陈鸿2002物理学报51 902]

[9]Cai Z,Lu W B,Liu Y J 2008 Acta Phys.Sin.57 7267(in Chinese)[蔡卓,陆文彬,刘拥军 2008物理学报 57 7267]

[10]Zhang Y L,Zhou B 2011 Acta Phys.Sin.60 120301(in Chinese)[张英丽,周斌 2011物理学报60 120301]

[11]Bode M,Heide M,von Bergmann K,Ferriani P,Heinze S,Bihlmayer G,Kubetzka A,Pietzsch O,Blügel S,Wiesendanger R 2007 Nature 447 190

[12]Romming N,Hanneken C,Menzel M,Bickel J,Wolter B,Bergmann K,Kubetzka A,Wiesendanger R 2013 Science 341 636

[13]Neto A H C,Guinea F,Peres N M R,Novoselov K S,Geim A K 2009 Rev.Mod.Phys.81 109

[14]Han W,Kawakami R K,Gmitra M,Fabian J 2014 Nat.Nanotech.9 794

[15]Roche S,Åkerman J,Beschoten B,Charlier J C,Chshiev M,Dash S P,Dlubak B,Fabian J,Fert A,Guimarães M,Guinea F,Grigorieva I,Schönenberger C,Seneor P,Stampfer C,Valenzuela O S,Waintal X,van Wees B 2015 2D Materials 2 030202

[16]Xu L,Dai Z H,Wang S,Liu B,Sun Y M,Wang W T 2014 Acta Phys.Sin.63 107102(in Chinese)[徐雷,戴振宏,王森,刘兵,孙玉明,王伟田2014物理学报63 107102]

[17]Xu L,Dai Z H,Sui P F,Wang W T,Sun Y M 2014 Acta Phys.Sin.63 186101(in Chinese)[徐雷,戴振宏,隋鹏飞,王伟田,孙玉明2014物理学报63 186101]

[18]Yang H X,Chen G,Cotta A A C,N’Diaye A T,Nikolaev S A,Soares E A,Macedo W A A,Schmid A K,Fert A,Chshiev M 2017 Ar Xiv 1704 09023

[19]PaciléD,Meyer J C,GiritÇ Ö,Zettl A 2008 Appl.Phys.Lett.92 133107

[20]Xu M S,Liang T,Shi M M,Chen H Z 2013 Chem.Rev.113 3766

[21]Guo S J,Dong S J 2011 Chem.Soc.Rev.40 2644

[22]Chen Q L,Dai Z H,Liu Z Q,An Y F,Liu Y L 2016 Acta Phys.Sin.65 136101(in Chinese)[陈庆玲,戴振宏,刘兆庆,安玉凤,刘悦林2016物理学报65 136101]

[23]Auwärter W,Muntwiler M,Greber T,Osterwalder J 2002 Surf.Sci.511 379

[24]Zhou Y G,Xiao-Dong J,Wang G,Xiao Y,Gao F,Zu X T 2010 Phys.Chem.Phys.12 7588

[25]Ma D W,Lu Z S,Ju W W,Tang Y N 2012 J.Phys.Condens.Matter 24 145501

[26]Kresse G,Joubert D 1999 Phys.Rev.B 59 1758

[27]Perdew J P,Burke K,Ernzerhof M 1996 Phys.Rev.Lett.77 3865