1.   引　言

过去十多年, 低维量子体系由于独特的结构和极强的可调性等特点, 成为凝聚态物理和材料科学研究的热点^[1,2]. 石墨烯是第一个由实验制备的稳定层状二维材料^[3]. 将两层石墨烯以不同角度堆叠可形成莫尔超晶格结构, 其性质对层数、层间的堆垛方式和衬底有极强依赖性^[4,5]. 两层石墨烯堆叠并旋转一定角度, 当转角减小到1.1°附近(此角度称为魔角), 体系的费米面附近出现平坦能带^[4]. 在平带中, 电子的运动速度趋于零, 具有极小的动能, 满足出现强关联的基本条件. 2018年, 实验研究人员在魔角石墨烯中观察到强电子-电子相互作用导致的莫特绝缘体和非常规超导态这两种关联量子物态^[6,7]. 随后, 通过门电压、转角、应变、垂直电场和外磁场等调控方式, 在石墨烯莫尔超晶格中实现了拓扑绝缘体^[8,9]、超导^[10]、铁电^[11]、铁磁^[12]以及量子反常霍尔 (quantum anomalous Hall, QAH) 效应^[13]等物理现象. 这些新奇现象迅速吸引了大批研究者的目光, 并发展出了一门新的科学分支, 即转角电子学 (twistronics)^[14].

最近实验和理论研究表明^[15-18], 当旋转双层石墨烯(twisted bilayer graphene, TBG)的转角小于1.2°时, 体系的晶格结构和电学性质发生明显的重构, 出现巨大的赝磁场^[19]. 随着转角逐渐减小, 晶格重构效应会进一步增强, 在极小角度的TBG中出现三角网格结构和一维拓扑态^[20]. 通常情况下, 将石墨烯放置在六方氮化硼(hexagonal boron nitride, hBN) 衬底上, 当二者晶向一致时, 会形成石墨烯/hBN莫尔超晶格. 研究发现hBN衬底会抑制魔角石墨烯中超导态的产生^[12,13,21]. 前期研究表明, 魔角石墨烯中超导态和绝缘态存在共生关系, 但在魔角石墨烯上叠加一层WSe₂时, 超导相和绝缘相由共生关系变成竞争关系, 且超导态能够在更小的转角下存在^[22]. 综上所述, 晶格弛豫和衬底会显著改变石墨烯莫尔超晶格的晶格结构和电子性质, 出现很多新的物理现象.

本文将围绕晶格弛豫和衬底对石墨烯莫尔超晶格性质调控的近期研究展开综述, 内容包括晶格弛豫对魔角TBG和魔角旋转三层石墨烯(twisted trilayer graphene, tTLG)的晶格结构和电学性质影响, 石墨烯和氮化硼莫尔超晶格中的能隙、超晶格狄拉克点的出现和反演对称性打破等现象.

2.   理论模型与方法

2.1   石墨烯莫尔超晶格的晶格结构

当单层石墨烯与另外一层石墨烯或其他层状二维材料以范德瓦耳斯力堆垛时, 在特定条件下会形成长周期的莫尔条纹(moiré pattern). 莫尔条纹的长度由两层材料的晶格失配和彼此之间形成的旋转角度 θ 决定. 具有严格实空间周期性的超晶格结构只存在于某些特定角度, 即所谓的“可公度”(commensurate)角度. 假设两层晶格结构的基矢分别为(${\boldsymbol{a}}_{\mathrm{t}1}$, ${\boldsymbol{a}}_{\mathrm{t}2}$)和(${\boldsymbol{a}}_{\mathrm{b}1}$, ${\boldsymbol{a}}_{\mathrm{b}2}$), 其对应倒格子基矢分别为(${\boldsymbol{k}}_{\mathrm{t}1}$, ${\boldsymbol{k}}_{\mathrm{t}2}$)和(${\boldsymbol{k}}_{\mathrm{b}1}$, ${\boldsymbol{k}}_{\mathrm{b}2}$). 这些基矢可以写成Park-Madden矩阵的形式^[23]:

$$\left(\genfrac{}{}{0pt}{}{{\boldsymbol{a}}_{\mathrm{t}1}}{{\boldsymbol{a}}_{\mathrm{t}2}}\right)=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\genfrac{}{}{0pt}{}{{\boldsymbol{a}}_{\mathrm{b}1}}{{\boldsymbol{a}}_{\mathrm{b}2}}\right) \text{, }$$

(1)

其中, 为了满足可公约条件, A, B, C, D均为有理数, 由莫尔条纹长度决定. 通常情况下, 样品的莫尔条纹可由扫描隧道显微镜(scanning tunneling microscope, STM)清晰地观测到. 为了精确求解(1)式中的矩阵元, 需要将超晶格结构的STM图进行(快速)傅里叶变换, 并定义莫尔超晶格的倒格子基矢为(${\boldsymbol{k}}_{\mathrm{m}1}$, ${\boldsymbol{k}}_{\mathrm{m}2}$), 则可公度性可用公度参数(i, j, k, l, m, n, q, r)描述为

$${\boldsymbol{k}}_{\mathrm{t}1}=i{\boldsymbol{k}}_{\mathrm{m}1}+k{\boldsymbol{k}}_{\mathrm{m}2} \text{, }$$

(2)

$${\boldsymbol{k}}_{\mathrm{t}2}=j{\boldsymbol{k}}_{\mathrm{m}1}+l{\boldsymbol{k}}_{\mathrm{m}2} \text{, }$$

(3)

$${\boldsymbol{k}}_{\mathrm{b}1}=m{\boldsymbol{k}}_{\mathrm{m}1}+q{\boldsymbol{k}}_{\mathrm{m}2} \text{, }$$

(4)

$${\boldsymbol{k}}_{\mathrm{b}2}=n{\boldsymbol{k}}_{\mathrm{m}1}+r{\boldsymbol{k}}_{\mathrm{m}2} .$$

(5)

由(1)式—(5)式可得:

$$\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)=\frac{1}{il-jk}\left(\begin{array}{cc}lm-jq& ln-jr\\ -km+iq& -kn+ir\end{array}\right) .$$

(6)

一旦莫尔条纹的周期(STM图中两个亮点间的距离, 如图2(a)所示)确定, 即可通过(2)式—(5)式求得一组公度参数, 且其解具有唯一性, 可进一步确定超晶格结构的应变$\varepsilon$和旋转角度$\theta$. 对于由两层六角蜂窝状晶格组成的超晶格结构, 假设两层的晶格矢量之比分别为${\boldsymbol{P}}_{1}=\dfrac{{\boldsymbol{a}}_{\mathrm{t}1}}{{\boldsymbol{a}}_{\mathrm{b}1}}$和${\boldsymbol{P}}_{2}=\dfrac{{\boldsymbol{a}}_{\mathrm{t}2}}{{\boldsymbol{a}}_{\mathrm{b}2}}$, 并且第二层相对于第一层的旋转角度分别为${\theta }_{1}=({\boldsymbol{a}}_{\mathrm{t}1}, {\boldsymbol{a}}_{\mathrm{b}1})$和${\theta }_{2}=({\boldsymbol{a}}_{\mathrm{t}2}, {\boldsymbol{a}}_{\mathrm{b}2})$. 则两层的晶格矢量可由推广的伍德记号(extended Wood’s notation)描述^[24]:

$$\left(\begin{array}{c}{\boldsymbol{a}}_{\mathrm{t}1}\\ {\boldsymbol{a}}_{\mathrm{t}2}\end{array}\right)=\left(\begin{array}{cc}{P}_{1}\left(\mathrm{cos}{\theta }_{1}+\dfrac{\mathrm{sin}{\theta }_{1}}{\sqrt{3}}\right)& \dfrac{2{P}_{1}}{\sqrt{3}}\mathrm{sin}{\theta }_{1}\\ -\dfrac{2{P}_{2}}{\sqrt{3}}\mathrm{sin}{\theta }_{2}& {P}_{2}\left(\mathrm{cos}{\theta }_{2}-\dfrac{\mathrm{sin}{\theta }_{2}}{\sqrt{3}}\right)\end{array}\right)\left(\begin{array}{c}{\boldsymbol{a}}_{\mathrm{b}1}\\ {\boldsymbol{a}}_{\mathrm{b}2}\end{array}\right) .$$

(7)

由(1)式和(7)式可求得参数${P}_{1}$, ${P}_{2}$, ${\theta }_{1}$和${\theta }_{2}$为

$${P}_{1}=\sqrt{{A}^{2}+{B}^{2}-AB} \text{, }$$

(8)

$${P}_{2}=\sqrt{{C}^{2}+{D}^{2}-CD} \text{, }$$

(9)

$${\theta }_{1}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{B\sqrt{3}}{2A-B}\right) \text{, }$$

(10)

$${\theta }_{2}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{C\sqrt{3}}{C-2D}\right) .$$

(11)

单轴异质应变(uniaxial heterostrain)和双轴异质应变(biaxial heterostrain)分别为

$${\varepsilon }_{\mathrm{b}\mathrm{i}}^{\mathrm{h}\mathrm{e}\mathrm{t}}=\sqrt{X-\sqrt{W}}-1 \text{, }$$

(12)

$${\varepsilon }_{\mathrm{u}\mathrm{n}\mathrm{i}}^{\mathrm{h}\mathrm{e}\mathrm{t}}=\sqrt{X+\sqrt{W}}-\sqrt{X-\sqrt{W}} \text{, }$$

(13)

其中$X=\dfrac{2\left({P}_{1}^{2}+{P}_{2}^{2}\right)+A}{3}$, $Y=\sqrt{\dfrac{4{P}_{1}^{2}{P}_{2}^{2}-{A}^{2}}{3}}$, $W={X}^{2}-{Y}^{2}$.

超晶格的原胞可能包含一个或多个莫尔晶胞. 若无特殊说明, 本文讨论的超晶胞原胞只包含一个莫尔晶胞. 以旋转角度为$\theta =6.01^{\circ}$的TBG为例, 图1为其结构示意图, 超晶胞的基矢${\boldsymbol{a}}_{\mathrm{m}1}$和${\boldsymbol{a}}_{\mathrm{m}2}$对应一组公度参数(i, j, k, l, m, n, q, r) = (11, 5, –5, 6, 11, 6, –6, 5), 推广的伍德记号为(${P}_{1}R{\theta }_{1}\times {P}_{2}R{\theta }_{2}$) = ($1 R6.01\times 1 R6.01$)(R是伍德记号标准形式中的R, 代表两层材料中对应基矢的旋转角度差, 即文中θ₁和θ₂), 一个原胞包含364个原子. 莫尔超晶格体系包含4种高对称结构, 分别是AA堆垛、AB堆垛、BA堆垛和畴壁(domain wall, DW)堆垛. 若TBG体系中不存在应变, 则莫尔条纹长度和转角满足简单关系式${{a}}_{\mathrm{m}1/\mathrm{m}2}=a/\left[2\mathrm{sin}(\theta /2)\right]$, 其中$a$是石墨烯的晶格常数.

2.2   低能有效模型

石墨烯莫尔超晶格的低能电子结构可由低能有效连续模型(low-energy effective continuum model)加以描述^[4,25,26]. 首先考虑多层转角石墨烯的层间跃迁项:

$$\hat{{\boldsymbol{U}}}=-\sum _{\alpha ,\beta }{t}_{\perp }^{\beta \alpha }({\boldsymbol{r}}_{\beta }-{\boldsymbol{r}}_{\alpha }){\hat{c}}_{\beta }^{†}{\hat{c}}_{\alpha } \text{, }$$

(14)

其中$\alpha =A, B$和$\beta =A', B\mathrm{'}$为上下层石墨烯的两套子晶格; ${\boldsymbol{r}}_{\alpha }$, ${\boldsymbol{r}}_{\beta }$分别是第一层$\alpha$轨道与第二层$\beta$轨道的位置; ${t}_{\perp }^{\beta \alpha }\left(\boldsymbol{r}\right)$是层间轨道的跃迁振幅. 在布洛赫表象下, 第一层$\alpha$子晶格上动量为$\boldsymbol{k}$的布洛赫态$|{\psi }_{\boldsymbol{k}}^{1, \alpha }\rangle$到第二层$\beta$子晶格上动量为$\boldsymbol{k}\boldsymbol{'}$的布洛赫态$|{\psi }_{\boldsymbol{k}\boldsymbol{'}}^{2, \beta }\rangle$的跃迁项可表示为^[26]

$$\begin{split}&{U}^{\beta \alpha }\left({\boldsymbol{k}}',\boldsymbol{k}\right)=\langle{\psi }_{{\boldsymbol{k}}'}^{2,\beta }|\hat{{\boldsymbol{U}}}|{\psi }_{\boldsymbol{k}}^{1,\alpha }\rangle = \\ \;&-\sum _{\boldsymbol{G},{\boldsymbol{G}}'}{\stackrel{~}{t}}_{\perp }^{\beta \alpha }\left(\boldsymbol{k}+\boldsymbol{G}\right){\mathrm{e}}^{-\mathrm{i}\left(\boldsymbol{G}\cdot{\boldsymbol{\tau }}_{\alpha }-{\boldsymbol{G}}'\cdot{\boldsymbol{\tau }}_{\beta }\right)}{\delta }_{\boldsymbol{k}+\boldsymbol{G},{\boldsymbol{k}}'+{\boldsymbol{G}}'}\text{, }\end{split}$$

(15)

其中, $\boldsymbol{G}={m}_{1}{\boldsymbol{k}}_{\mathrm{t}1}+{m}_{2}{\boldsymbol{k}}_{\mathrm{t}2}$, ${\boldsymbol{G}}'={m}_{1}'{\boldsymbol{k}}_{\mathrm{b}1}+{m}_{2}'{\boldsymbol{k}}_{\mathrm{b}2}$为上下两层的倒格矢; ${\boldsymbol{\tau }}_{\alpha }$, ${\boldsymbol{\tau }}_{\beta }$为子晶格原胞的位矢; ${\stackrel{~}{t}}_{\perp }^{\alpha , \beta }\left(\boldsymbol{q}\right)$为实空间跃迁振幅${t}_{\perp }^{\alpha , \beta }\left(\boldsymbol{r}\right)$的傅里叶变换,

$${\stackrel{~}{t}}_{\perp }^{\beta \alpha }\left(\boldsymbol{q}\right)=\frac{1}{S}\int \mathrm{d}\boldsymbol{r}{t}_{\perp }^{\beta \alpha }\left(\boldsymbol{r}\right){\mathrm{e}}^{-\mathrm{i}\boldsymbol{q}\cdot\boldsymbol{r}} .$$

(16)

(16)式表明, 当$\boldsymbol{q}$长度大于倒格子基矢长度$\left|{\boldsymbol{k}}_{\mathrm{t}1}\right|$时, 动量空间跃迁振幅${\stackrel{~}{t}}_{\perp }^{\beta \alpha }\left(\boldsymbol{q}\right)$会快速衰减, 故(15)式中有效的跃迁项只发生在第一层石墨烯中动量为$\boldsymbol{k}$的态与第二层动量为$\boldsymbol{k}$, $\boldsymbol{k}$+${\boldsymbol{k}}_{\mathrm{m}1}$和$\boldsymbol{k}$+${\boldsymbol{k}}_{\mathrm{m}1}$+${\boldsymbol{k}}_{\mathrm{m}2}$的态之间, 以及第二层石墨烯中动量为$\boldsymbol{k}$的态与第一层动量为$\boldsymbol{k}$, $\boldsymbol{k}{-\boldsymbol{k}}_{\mathrm{m}1}$和$\boldsymbol{k}{-\boldsymbol{k}}_{\mathrm{m}1}{-\boldsymbol{k}}_{\mathrm{m}2}$的态之间.

在长波近似下, 双层转角石墨烯的低能能谱可用六个动量值的哈密顿量来描述, 即第一层中的$\boldsymbol{k}$, $\boldsymbol{k}-\boldsymbol{k}_{\mathrm{m}1}$和$\boldsymbol{k}-\boldsymbol{k}_{\mathrm{m}1}-\boldsymbol{k}_{\mathrm{m}2}$和第二层的$\boldsymbol{k}, \boldsymbol{k}+{\boldsymbol{k}}_{\mathrm{m}1} \mathrm{和} \boldsymbol{k}+{\boldsymbol{k}}_{\mathrm{m}1}+{\boldsymbol{k}}_{\mathrm{m}2}$. 此时体系低能部分有效哈密顿量可约化为以($A, B, A\mathrm{'}, B\mathrm{'})$为基矢的4×4矩阵^[26,27]:

$${\boldsymbol{H}}^{\left(\xi \right)}=\left(\begin{array}{cc}{H}_{1}& {U}^{†}\\ U& {H}_{2}\end{array}\right) .$$

(17)

这里, ${H}_{l=\mathrm{1, 2}}$是l层内哈密顿量, 由外尔方程给出:

$${H}_{l}\left(\boldsymbol{k}\right)=-\hbar v \Big[R\Big(\pm \frac{\theta }{2}\Big)\left(\boldsymbol{k}-{\boldsymbol{K}}_{\xi }^{l}\right)\Big]\cdot(\xi {\sigma }_{x},{\sigma }_{y}) .$$

(18)

$U$是有效层间耦合, 由三项层间耦合组成:

$$U=\left(\begin{array}{cc}{U}^{{A}^\prime A}& {U}^{{A}^\prime B}\\ {U}^{{B}^\prime A}& {U}^{B^\prime B}\end{array}\right) .$$

(19)

根据(17)式中的哈密顿量, 可以得到TBG的低能能带色散关系.

在有效连续模型中, 晶格弛豫和应变对电子结构的影响可通过对(17)式中矩阵元的调整加以研究. 如文献[27]中讨论了石墨烯莫尔超晶格层间结构弛豫的影响, 由于石墨烯莫尔超晶格AA堆垛区域的层间结合能较AB和BA堆垛区域更高^[28-30], 莫尔超晶格自发弛豫后AA堆垛区域的层间距更大, 从而使(19)式中子晶格间的跃迁振幅不再相等来调整能带结构.

2.3   紧束缚模型与大尺度计算方法

紧束缚模型(tight-binding model)是凝聚态物理中广泛应用的理论模型, 对于周期性比较小的体系可以得到直观的电子结构, 同时也可以处理周期性很大或没有周期性的超大体系. 单电子近似下, 紧束缚模型的哈密顿量可以表示为

$$\hat{H}=-\sum \nolimits_{i,j}{t}_{ij}{{\hat{c}}_{i}}^{†}{\hat{c}}_{j}+\sum \nolimits_{i}{V}_{i}{{\hat{c}}_{i}}^{†}{\hat{c}}_{i} \text{, }$$

(20)

其中${t}_{ij}$表示i, j轨道间的跃迁振幅; ${V}_{i}$表示第i个轨道的占位能; ${{\hat{c}}_{i}}^{†}, {\hat{c}}_{i}$为第i个轨道上的产生、湮灭算符. 对于石墨烯体系, 为了研究结构弛豫对电子结构的影响, 可使用Slater-Koster近似计算轨道间的跃迁振幅^[31]:

$$\begin{split}&-t\left(\boldsymbol{R}\right)={V}_{\mathrm{p}\mathrm{p}\mathrm{\pi }}[1-{({{\boldsymbol{R}}}_{z}/\boldsymbol{R})}^{2}]+{V}_{\mathrm{p}\mathrm{p}\sigma }{({{\boldsymbol{R}}}_{z}/\boldsymbol{R})}^{2} \text{, }\\ &{V}_{\mathrm{p}\mathrm{p}\mathrm{\pi }} = {V}_{\mathrm{p}\mathrm{p}\mathrm{\pi }}^{0}{\mathrm{e}}^{-({R}-\tfrac{a}{\sqrt{3}})/{r}_{0}} , {V}_{\mathrm{p}\mathrm{p}\sigma } = {V}_{\mathrm{p}\mathrm{p}\sigma }^{0}{\mathrm{e}}^{-({R}-d)/{r}_{0}} , \end{split}$$

(21)

其中, $\boldsymbol{R}$为轨道间距离, ${{\boldsymbol{R}}}_{z}$为z方向的分量; ${r}_{0}$是跃迁衰减长度; ${V}_{\mathrm{p}\mathrm{p}\mathrm{\pi }}^{0}$是单层材料中最近邻轨道间跃迁能; ${V}_{\mathrm{p}\mathrm{p}\sigma }^{0}$是层间垂直方向上两个轨道的跃迁能^[32].

由于莫尔超晶格的晶胞包含原子数目过多, 使用对角化方法求其本征态在数值上会遇到较大的困难. 针对复杂的大尺度体系, 可采用基于含时薛定谔方程的紧束缚传播方法(tight-binding propagation method, TBPM)绕开对本征态的计算^[19,33-35]. TBPM是建立在对含时薛定谔方程的精确数值求解之上, 不对体系进行任何对角化处理, 可以对体系的各种物理性质进行快速而准确地计算, 且计算消耗和体系大小呈线性关系. 如运用Chebyshev多项式展开的方法, 可以从一个初始波函数直接准确地计算出其他任意时刻的波函数. 首先介绍在TBPM中对时间演化的处理方法. 根据含时薛定谔方程, 对于随机选取的初始态波函数$|{\psi }_{0}\rangle= \sum _{i}{c}_{i}|i\rangle$, 其时间演化态为

$$|{\psi }_{0}\left(t\right)\rangle={\mathrm{e}}^{-\mathrm{i}\hat{H}t/\hbar }|{\psi }_{0}\rangle .$$

(22)

在数值计算中, 使用有限大小的固定时间步长${\Delta }t$, 通过Chebyshev多项式对时间演化算符进行展开可得到^[33]

$$\begin{split}&|{\psi }_{0}\left(t+\Delta t\right)\rangle\\ =\;&\left[{\mathrm{J}}_{0}\left({\Delta }t\right){T}_{0} (\hat{H})+2\sum _{m=1}^{\infty }{\mathrm{J}}_{m}\left({\Delta }t\right){T}_{m}(\hat{H})\right]|{\psi }_{0}\left(t\right)\rangle \text{, }\end{split}$$

(23)

其中, ${\mathrm{J}}_{m}$为$m$阶第一类贝塞尔函数; ${T}_{m}$为$m$阶切比雪夫多项式, 可通过递推关系迭代求得, 即

$${T}_{m+1}(\hat{H})|\psi \rangle=-2\mathrm{i}\hat{H}{T}_{m} (\hat{H})|\psi \rangle + {T}_{m-1} (\hat{H})|\psi \rangle\text{, }$$

(24)

$${T}_{0} (\hat{H})|\psi \rangle=|\psi \rangle\text{, }~{T}_{1}(\hat{H})|\psi \rangle=-\mathrm{i}\hat{H}|\psi \rangle .$$

(25)

通过(23)式, 可求得一系列基于同一个初始态的时间演化波函数. 体系的电学、光学、输运和等离激元性质可以进一步通过波函数的含时演化得到. 下面以电子态密度(density of state, DOS)为例介绍如何通过时间演化方法计算体系物理性质:

$$D\left(E\right)=\sum _{n=1}\delta (E-{E}_{n}) .$$

(26)

利用性质:

$$\delta \left(E-{E}_{n}\right)=\frac{1}{2\mathrm{\pi }}\int _{-\infty }^{+\infty }{\mathrm{e}}^{\mathrm{i}(E-{E}_{n})t/\hbar }\mathrm{d}t \text{, }$$

(27)

通过对含时演化态与初始态的内积$\langle{\psi }_{0}|{\mathrm{e}}^{-\mathrm{i}Ht/\hbar }|{\psi }_{0}\rangle$进行傅里叶变换, 并对基于不同的随机初始波函数的结果进行统计平均^[36], 便可得到DOS:

$$D\left(E\right)=\underset{S\to \infty }{\mathrm{lim}}\frac{1}{S}\sum _{p=1}^{S}\frac{1}{2\mathrm{\pi }}\int _{-\infty }^{+\infty }{\mathrm{e}}^{\mathrm{i}Et/\hbar }\langle{\psi }_{0}|{\mathrm{e}}^{-\mathrm{i}\hat{H}t/\hbar }|{\psi }_{0}\rangle\mathrm{d}t .$$

(28)

类似地, 通过对不同时间的波函数本身进行傅里叶变换, 可以得到特定能量的准本征态(quasi-eigenstates):

$$|\psi \left(E\right)\rangle=\underset{S\to \infty }{\mathrm{lim}}\frac{1}{S}\sum _{p=1}^{S}\frac{1}{2\mathrm{\pi }}\int _{-\infty }^{+\infty }{\mathrm{e}}^{\mathrm{i}Et/\hbar }|{\psi }_{0}\left(t\right)\rangle\mathrm{d}t .$$

(29)

目前, TBPM方法已集成在开源软件TBPLaS (Tight-Binding Package for Large-scale Simulation)^[37]中. 这是一款通用的大型紧束缚模型计算软件, 不仅可以完成能带计算, 还可以应用到超大体系的各种电学、光学、输运和等离激元性质计算, 模拟体系的尺寸跨度达到近10个数量级, 可以模拟数十亿个原子(格点)的复杂量子体系, 非常适合莫尔超晶格体系的电学、光学、输运和等离激元性质的计算.

3.   石墨烯莫尔超晶格中的自发弛豫

旋转角度作为新的自由度, 可以有效调控旋转石墨烯体系的电学性质, 给石墨烯体系带来许多新奇的物理现象. 本节主要研究旋转石墨烯体系的弛豫现象以及弛豫对不同转角体系的晶格结构和电学性质的影响.

3.1   TBG莫尔超晶格

在TBG能带图中, K和K'点的中间位置有一个鞍点M, 对应DOS分布的极大值, 称为范霍夫奇点(van Hove singularity, VHS)^[5,38,39]. 实验和理论研究发现, TBG的VHS出现在低能部分, 且其能量位置随旋转角度减小逐渐接近费米能级^[5,38,39]. 例如, 当旋转角度减小到$1.16^{\circ}$时, VHS离费米能级的距离只有6 meV. 当旋转角度接近$1.05^{\circ}$, 体系在费米面附近出现平坦的能带, 对应的角度称为第一魔角^[4,39]. 魔角的定义主要有两种: 一是超结构倒空间中K和K'点的费米速度降为0; 二是体系在费米面附近的能带带宽最小, 并且最低能量电子态距离高能量电子态的带隙值最大^[40]. 理论研究发现, TBG有一系列魔转角, 且魔转角之间满足${\theta }_{n}^{*}={\theta }_{1}^{*}/n$, 其中${\theta }_{1}^{*}=1.05^{\circ}$是第一魔角^[4,39]. 由于理论模型和计算参数的差异, 第一魔转角的数值并不是定值, 而是出现在1.0°—1.5°之间. 2018年, 麻省理工学院Jarillo-Herrero课题组^[6,7]成功制备出魔角石墨烯, 并在体系中观察到强电子-电子相互作用导致的莫特绝缘体和非常规超导态这两种关联量子物态. 随后, 在这个强关联体系中又发现了铁电、铁磁、拓扑、QAH效应等一系列强关联现象^[8,11].

两层石墨烯以AA堆垛并旋转一定角度后, 体系中存在AA, AB, BA和DW四种高对称结构, 如图1所示. 其中AB和BA堆垛(又称Bernal堆垛)是能量最低的堆垛方式, DW堆垛能量次之, AA堆垛的结合能最高^[15]. 一方面, 不同堆垛的能量差异产生了面内应力, 使原子偏离平衡位置, 从而增加Bernal堆垛区域的面积, 减小层间相互作用能; 另一方面, 原子在面内移动产生的应力会阻止其进一步偏离原来位置. 基于以上原因, 莫尔超晶格中原子会偏离其原有位置, 发生结构形变, 因此, 晶格重构是层间能量达到最小值和位移产生应力的竞争结果^[15,41]. 理论研究发现, 当TBG的转角小于1.2°时, 体系发生明显的晶格弛豫^[15], 影响小角度TBG的晶格结构和电学性质^[17,18,42]. 对于魔角石墨烯体系, 晶格弛豫使平带的带宽由6 meV变为9 meV, 并且在Γ点打开带隙, 能隙值约为42 meV^[43]. Yoo等^[18]利用透射电子显微镜研究了小角度TBG的结构, 发现晶格弛豫使样品重构成一系列三角形的AB和BA堆垛区域, 同时, 在AB和BA交界处产生畴壁, 不同区域的畴壁形成一个三角网格. 畴壁上存在无能隙的一维拓扑电子态^[44], 可以增强畴壁区域的光电导^[45]. 结合STM探测手段、扫描隧道谱(scanning tunnel spectrum, STS)技术和理论分析, 目前已经确认TBG中晶格弛豫会引起自发的赝磁场^[19].

如图2(a)所示, TBG的STM形貌图展现出完美的六角对称性, 莫尔条纹长度通过测量亮点间距得到, 分别为$L1\approx L2\approx L3\approx 29.6$ nm. 随后求解(1)式—(6)式得到一组公度参数(I, j, k, l, m, n, q, r) = (139, 69, –69, 70, 139, 70, –70, 69), 并求得其旋转角度$\theta =0.48^{\circ}$. 在莫尔条纹区域不同的位置进行STS谱探测, 如图2(b)所示, AB堆垛区域的STS谱没有明显的特征峰, 呈现出AB堆垛双层石墨烯的性质. 另一方面, AA堆垛区域的STS谱展现出一系列等能量间距的特征峰, 与之前报道的结果截然不同^[38]. 早先的理论与实验研究表明, 在角度$1.4^{\circ} \leqslant \theta \leqslant 9.6^{\circ}$的TBG样品中, 其STS谱在费米能级附近出现两个VHS峰, 而不会出现一系列等能峰.

为了进一步探究STS中等能量特征峰背后的物理起源, Shi等^[19]对这一结构进行计算分析. STS谱反映样品表面的局域电子态密度(local density of states, LDOS). 图2(c)为采用紧束缚模型计算得到的LDOS. 对于非弛豫结构的TBG, AA堆垛区域的LDOS不存在等能量特征峰. 随后, 通过在体系中加入晶格弛豫修正, 新得到的理论结果与实验十分相符. 弛豫结构中原子在平面内和平面外的位移满足如下公式^[19]:

$$\begin{split}\Delta d\left(r\right)=\;&\Delta D\cdot \mathrm{sin}\left(\frac{\mathrm{\pi }r}{2{l}_{D}}\right)\left[1-\varTheta \left(r-{l}_{D}\right)\right]\\ &+\Delta D\cdot {\mathrm{e}}^{-{\left(r-{l}_{D}\right)}^{2}/{\sigma }_{D}}\varTheta\left(r-{l}_{D}\right) \text{, } \\ \Delta z\left(r\right)=\;&\pm \Delta Z\cdot {\mathrm{e}}^{-{r}^{2}/{\sigma }_{Z}}\left[1-\varTheta (r-{l}_{Z})\right] \text{, }\end{split}$$

(30)

其中, $r$是x-y平面内弛豫原子距离最近邻AA区域的长度; Θ是Heaviside阶跃函数; $\Delta D$和$\Delta Z$分别是平面内和平面外最大位移量, 其数值取决于旋转角度$\theta$^[17]; ${l}_{D}$, ${\sigma }_{D}$, ${l}_{Z}$和${\sigma }_{D}$是常数. 图2(g)显示的是原子沿路径M-N-P的面内和面外位移, 结果显示AA堆垛区域的原子形变最大, 与分子动力学弛豫结果一致^[17]. 在旋转石墨烯体系中, 晶格弛豫在体系哈密顿矩阵中引入两个势参数, 形变势能$V$和赝矢势$\boldsymbol{A}$, 其表达式为

$$\begin{split}&V=2{g}_{1}\frac{{(d}_{1}+{d}_{2}+{d}_{3})/3-d}{d} \text{, } \\ &{A}_{x}=c\frac{\beta {\gamma }_{0}}{d}\left({u}_{xx}-{u}_{yy}\right) \text{, } \\ &{A}_{y}=-c\frac{2\beta {\gamma }_{0}}{d}{u}_{xy} \text{, } \end{split}$$

(31)

其中, 石墨烯的形变势能${g}_{1}\approx 4$ eV; ${d}_{1}$, ${d}_{2}$和${d}_{3}$是弛豫结构中最近邻原子间距; 石墨烯的碳-碳键长$d=0.142$ nm; 常数因子$c\approx 1$, 取决于具体模型中的化学键长; Grüneisen参数$\beta \approx 2$; ${\gamma }_{0}$是平面内最近邻相互作用强度; ${u}_{ij}$是应变张量. 形变势能V以占位能的形式加入哈密顿矩阵, 赝矢势$\boldsymbol{A}$通过改变原子轨道间的跃迁能来修正体系哈密顿量. 对赝矢势求旋度可得到赝磁场$\boldsymbol{B}=\nabla \times \boldsymbol{A}$, 其绝对值近似为^[41]

$$\left|\boldsymbol{B}\right|\approx \frac{\hbar c}{\mathrm{e}}\left(\frac{2\beta }{3}\right)\frac{\sqrt{3{({d}_{3}-{d}_{2})}^{2}+{({d}_{2}+{d}_{3}-2{d}_{1})}^{2}}}{2{a}_{\mathrm{m}}{d}^{2}},$$

(32)

${a}_{\mathrm{m}}$是莫尔长度. 图2(d)和图2(e)是采用(31)式和(32)式计算的形变势能和赝磁场值. 对于转角$\theta =0.48^{\circ}$的TBG, 晶格弛豫在体系中产生不均匀的形变势能和赝磁场. 其中形变势能分布在AA堆垛区域, 最大值是0.2 eV. 类似地, 非均匀的赝磁场以AA堆垛中心为原点, 形成涡旋状磁场, 赝磁场最大值是6 T. AA堆垛区域的电子在赝磁场作用下出现了类似于真实磁场作用下的朗道量子化能级. 但是, 晶格弛豫产生的赝磁场没有打破时间反演对称性. 值得注意的是, 体系中的赝磁场是由于晶格自发弛豫产生的, 并且随着转角减小, 赝磁场值越大. 下面进一步讨论晶格弛豫对不同魔转角体系的原子结构和电子性质的影响.

魔角石墨烯的特别之处是在费米面附近出现平坦的能带^[4], 在平带中电子的运动速度趋于零, 具有极小的动能, 满足出现强关联的基本条件. 同时, 平带的电子态在实空间中高度局域在AA堆垛区域. 目前, 魔角双层石墨烯的实验和理论研究主要集中在第一魔角附近, 更高阶魔角的研究报道比较少. Huang等^[20]通过STM观察到旋转角度为0.25° (第四魔角附近)的双层石墨烯结构, 并利用STS探测到拓扑边界态. 在小角度TBG中, 晶格弛豫使样品重构成一系列AB和BA堆叠区域, 并在AB-BA堆叠的交界处产生畴壁. 在外加偏置电压下, 双层石墨烯可在狄拉克点打开能隙并获得非零的能谷陈数. 畴壁两侧AB和BA两种堆叠结构在某一特定能谷的能谷陈数符号相反, 导致畴壁上存在无能隙的一维拓扑电子态, 并形成一维导电通道^[20,46]. Xu等^[47]观察到电子在一维拓扑导电通道中产生的阿哈罗诺夫-玻姆效应(Aharonov-Bohm effect). 同时, 畴壁可以传输谷极化电流^[48]. 晶格重构会调节材料中的振动和电子结构, 导致电-声耦合行为发生变化^[42]. Lu等^[49]在第二魔角附近观察到一系列孤立的平带, 并通过强磁场下的输运性质得到进一步证实. 理论研究发现, 随着角度减小, 由于晶格重构, AA堆叠区域急剧减小, 从而导致高阶魔角体系中不存在孤立的平带^[50]. 结合分子动力学模拟和紧束缚模型, 在这里系统地阐述晶格弛豫对不同魔角石墨烯的结构、平带以及平带波函数在实空间分布的调控规律.

首先, 利用分子动力学模拟软件LAMMPS^[51]对不同魔角体系进行晶格弛豫, 然后采用(21)式将原子位移对层内和层间电子跃迁积分的修正加入紧束缚模型中. 通过调节(21)式的紧束缚参数, 将第一魔角对应的旋转角度值定为$\theta =1.05^{\circ}$. 在进行LAMMPS弛豫时, 层内相互作用和层间相互作用分别采用LCBOP^[52]和Kolmogorov-Crespi^[53]势参数. 若无特殊说明, 本文均采用相同的层内层间势参数对石墨烯莫尔超晶格进行晶格弛豫. 这里主要研究三类魔角体系: 第一魔角($\theta =1.05^{\circ}$)、第二魔角($\theta =0.53^{\circ}$) 和第三魔角($\theta =0.35^{\circ}$), 对应公度参数(I, j, k, l, m, n, q, r)分别是(63, 31, –31, 32, 63, 32, –32, 31), (125, 62, –62, 63, 125, 63, –63, 62)和(189, 94, –94, 95, 189, 95, –95, 94), 相 应晶胞所包含的原子数是11908, 46876和107164. 由于这些魔角体系晶胞包含的原子数众多, 无论是基于密度泛函理论的第一性原理计算还是基于对角化的紧束缚模型计算都难以直接进行. 因此, 采用TBPM方法来计算魔角石墨烯的LDOS和DOS在实空间的分布.

图3是不同魔角体系中原子在平面内和平面外的弛豫结果. 由于AA堆垛的结合能最高, 因此AA区域的原子发生明显偏移. 原子在平面内围绕着AA旋转, 形成“涡旋”结构, 并且上下层原子旋转方向相反. 同时, TBG在AA区域形成“突起”结构, 即原子朝面外移动. 因此, AA堆垛的层间距最大, AB堆垛的层间距最小. 这些“突起”结构即STM形貌图观察到的亮点. 随着角度减小, 例如在高阶魔角体系中, AA区域面积逐渐减小, AB和BA区域面积逐渐增大, 同时AB区域由圆形变为三角形, 并且在AB与BA区域间看到明显的畴壁. 特别是在第三魔角, 体系主要由一维拓扑畴壁连接不同的AA堆垛区域, 不同堆叠区域的畴壁形成一个三角网格, 三角网格内是交错排列的AB或BA堆垛结构. 从图3可以明显观察到, 由于晶格弛豫强度不同, 高阶魔角与第一魔角在结构上有显著的变化.

为进一步展示晶格弛豫对魔角体系电学性质的影响, 利用紧束缚模型计算不同魔角体系的DOS和LDOS. 如图4所示, 晶格弛豫对魔角体系的电子性质有显著影响. 和晶格重构强度的变化规律类似, 随着角度减小, 晶格弛豫对电子性质的影响越发明显. 首先, 我们研究晶格弛豫对DOS的影响. 对于第一魔角, 如图4(a)所示, 晶格弛豫导致费米面附近出现孤立的平带, 打开的带隙值约为42 meV. 对比非弛豫结构和弛豫结构在费米面附近的DOS分布, 发现晶格弛豫没有明显改变平带能量态密度. 理论上讲, 若第一魔角是$\theta =1.05^{\circ}$, 则第二魔角出现在$\theta =1.05^{\circ}/2=0.53^{\circ}$附近. 随后计算了转角$\theta =0.53^{\circ}$体系的DOS (图4(d)), 发现在非弛豫情况下, 费米面附近出现了两个特征峰. 体系弛豫后, 费米面附近只存在一个峰, 且峰的高度明显降低. 而第三魔角弛豫后的体系在费米面附近出现了一系列特征峰. 同时, 高阶魔角体系的能量态密度中没有出现能隙, 说明这些低能部分的平带都不是孤立的, 这一现象与第一魔角截然不同. 以上结果表明, 随着旋转角度减小, 体系电子性质受晶格弛豫影响会越发显著.

接下来, 考虑不同魔角体系中, 晶格弛豫对AA, AB和DW三种高对称区域LDOS的影响. 在魔角体系中, 不同区域电子密度分布有显著不同. 例如, 费米面附近的电子态主要局域在AA区域, 只有少部分分布在畴壁上. 对于高阶魔角的理想结构(非弛豫), AA区域态密度在高能部分出现一系列特征峰. 但是, 在弛豫结构中, 这些高能量的峰态几乎消失, 说明晶格弛豫减小了AA区域态密度对体系总态密度的贡献. 这些特征峰可以通过计算态密度在实空间的分布来进一步证实. 图4(b)和图4(c)是不同特征峰电子态在实空间的分布. 可以明显发现, 费米面附近平带的电子态主要局域在AA区域. 随着角度减小, AA区域局域电子态面积逐渐减小, 且晶格弛豫会进一步减小电子态局域面积. 对比非弛豫结构和弛豫结构的结果可看出, 晶格弛豫对电子态在实空间的分布有极大的影响. 以上计算结果说明晶格弛豫对高阶魔角的晶格构型和电学性质有更显著的影响, 这可能是实验中仍未发现高阶魔角超导特性的原因.

3.2   tTLG莫尔超晶格

实验和理论研究发现, 除TBG之外, ABC堆垛三层石墨烯^[54-57]、tTLG^[58-61]、双层-双层旋转石墨烯^[62,63]等结构同样能在石墨烯体系中引入平带, 并且对这些平带体系的研究中也发现了各种新奇的量子物态. 这一部分将重点介绍tTLG的平带以及这一平带体系中观察到的强关联量子物态. 相比于TBG, tTLG拥有更多可调节的自由度. 首先, tTLG有多种初始结构, 例如AAA堆垛、AAB堆垛、ABA堆垛和ABC堆垛等. 其次, tTLG的性质可由两个独立旋转角度${\theta }_{1}$和${\theta }_{2}$进行调控, 其中${\theta }_{1}$是第一层和第二层之间的旋转角度, ${\theta }_{2}$是第二层和第三层之间的角度. Wu等^[64]研究了不同初始结构的tTLG, 发现初始构型和旋转层的选取对其电学性质影响显著, 并在AÃA (初始结构是AAA堆垛, 旋转中间层)和ÃAB (初始结构是AAB堆垛, 旋转第一层) 两种构型中观察到了平坦能带. Carr等^[58]采用低能有效模型对tTLG-AÃA这种垂直镜对称性“三明治”结构进行研究, 在低能部分发现平带和狄拉克锥共存的现象. Park等^[65]发现tTLG-AÃA具备超导性能, 当中间层相对于其他两层以一定角度旋转时, 其超导性比双层结构更稳定. 同时Cao等^[66]研究了这种超导现象的自旋特性, 发现在高达10 T的磁场下体系仍表现出超导性, 是一种自旋三重态. Hao等^[67]以“魔角”依次旋转每层石墨烯来构建tTLG-AÃA结构, 观察到具有位移场可调谐的超导性, 其超导转变温度是2.1 K. 值得一提的是, 与传统的超导材料不同, tTLG中的超导性不是由常规的弱耦合Bardeen-Cooper-Schrieffer (BCS) 电子对机理引起的, 而是强烈依赖于电子的强耦合. Zhang等^[68]研究了tTLG的超莫尔晶格结构 (${\theta }_{1}={\theta }_{2}$), 发现体系在极低的电子密度(低于两个数量级)下仍表现出超导的传输特性. Turkei等^[69]使用STM直接表征魔角tTLG-AÃA的电子结构, 观察到体系具有基于莫尔晶格的强烈重构, 引入了一系列局部扭转角断层. 在非镜像对称堆叠结构中, 例如tTLG-ÃAB, 研究者观察到了非对称的关联绝缘态、铁磁态和类似超导的迹象^[8,70,71]. 随后, 研究tTLG-AÃA和tTLG-ÃAB两种构型中晶格弛豫对平带的影响.

如图5(a)所示, tTLG-AÃA构型的魔角$\theta = 1.35^{\circ}$(tTLG-AÃA-1.35体系), 是TBG的$\sqrt{2}$倍^[64]. 晶格弛豫后, 体系仍保持双重旋转对称性C_2z. 在AAA区域, 上层和下层的原子分别往面外移动. 从图5(d) DOS分布可以看出, 弛豫的体系在费米面附近打开能隙, 能隙值约为55.4 meV. 平带的电子态主要分布在AAA区域(图5(e)和图5(f)), 且中间层的DOS大于外层. tTLG-ÃAB可以看成由两部分组成, 下层是AB堆垛型双层石墨烯, 上层是相对于下层有一个转角的单层石墨烯. 因此, 第一层原子在AAB区域的z方向位移最大. 相比于tTLG-AÃA结构, 原子在z方向位移均有所减小. tTLG-ÃAB的魔角$\theta =1.05^{\circ}$ (tTLG-ÃAB-1.05体系) , 与TBG的魔角一致. 晶格弛豫后, 体系在费米面附近打开能隙, 约23.4 meV. 从图5(g)可以看出, 晶格弛豫使平带对应的DOS降低. 平带的电子态主要分布在AAB堆垛的上层和中间层. 以上结果表明, 在tTLG中, 晶格弛豫对由不同初始结构和不同旋转层构成体系的晶格结构和电子性质有不同的影响.

目前一些研究表明, 压强是调节旋转石墨烯性质的一种有效手段. Yankowitz等^[10]通过静水压力改变层间间距, 在旋转角度为$\theta =1.27^{\circ}$的TBG中诱导超导态, 观察到最大超导转变温度${T}_{\mathrm{c}}$约3 K. 这里, 通过具体计算阐述外部压强对tTLG物性的调控作用. 如图6(b)所示, tTLG-AÃA-1.89本身不存在平带, 且费米面附近电子态主要分布在中间层^[72]. 当对体系施加4 GPa的压强, 增强的层间相互作用使每一层费米面附近的电子态向AAA区域聚集, 从而诱导产生新的平带. 但是晶格弛豫后, 体系只打开6 meV的能隙, 远远小于tTLG-AÃA-1.35的能隙. 这个区别可由晶格弛豫后体系结构的不同重构来解释. 如图6(c)所示, 在4 GPa压强下, 晶格弛豫并未改变原子在z方向的分量, 仍类似于非弛豫结构, 使得AA堆垛与AB堆垛的层间距相同, AA区域的层间相互作用${\omega }_{0}$接近于AB的层间相互作用${\omega }_{1}$. 而对于tTLG-AÃA-1.35, ${\omega }_{1}/{\omega }_{0}$比值远大于1. 能隙的大小决定于${\omega }_{1}$与${\omega }_{0}$的相对大小, 具体来说, 正比于${\omega }_{1}/{\omega }_{0}$. 因此, 图6(d)中tTLG-AÃA-1.89加压后的能隙远小于tTLG-AÃA-1.35的能隙. 对于旋转角度为1.89°的TBG, 诱导产生平带所需的强压约为9 GPa, 约为tTLG的两倍^[73]. 上述结果进一步表明tTLG相对于TBG具有更加可调的晶体结构和电子性质.

4.   石墨烯-hBN异质结构中的自发弛豫

石墨烯上具有很多优异性能, 其应用的难点之一是高效地制备出高质量样品. 早期石墨烯的生长制备主要在SiO₂衬底表面, 但是由于其平整度不高, 制备的石墨烯样品会出现褶皱^[74]. 金属衬底上制备的石墨烯样品具有优异的力学性能, 但石墨烯低能部分的电子态会与衬底的金属态发生杂化, 从而影响样品的电子结构^[75]. 2010年Dean等^[76]在hBN衬底上制备了石墨烯样品, 显著提高了样品的力学和电学性能. hBN的晶体结构与石墨烯类似, 主要不同是hBN由氮原子和硼原子分别占据一个子晶格构成, 这种结构使其具有跟石墨烯完全不同的电子性质. hBN具有原子级平整的表面, 是一种绝缘体材料, 具有较大的能隙, 没有悬挂键, 光学声子的能量较高, 是目前制备高质量石墨烯样品常用的衬底.

随后, 研究者发现石墨烯和hBN会以范德瓦耳斯力堆垛形成异质结, 其电子性质与二者之间的旋转角度密切相关. 当旋转角度小于2°, 石墨烯表面会形成较大周期的莫尔条纹^[77-79]. 莫尔条纹可以视为施加在石墨烯上的周期势, 可对石墨烯的能带进行调控. 如在对石墨烯/hBN异质结输运性质的研究中, 研究者发现了其莫尔超晶格能级随着外界磁场的变化形成分形结构, 实现了Hofstadter提出的“霍夫斯塔特蝴蝶”(Hofstadter’s butterfly)能谱^[80-83]. 氮化硼衬底对单层石墨烯的性质调控主要有: 1) 在石墨烯狄拉克点打开一个能隙; 2) 在较高能量形成一系列超晶格狄拉克点. 超晶格狄拉克点的出现源于石墨烯上的莫尔超晶格势, 而石墨烯同hBN之间旋转角度越小, 超晶格狄拉克点同石墨烯狄拉克点的距离越近. 氮化硼/石墨烯/氮化硼形成的超莫尔三明治结构中同样存在超晶格狄拉克点^[84-87].

hBN的两个子晶格分别由两种类型的原子占据, 当石墨烯与hBN的晶向基本一致时, 理论上可以通过hBN的两个子晶格来调节石墨烯子晶格的占位能, 从而打破石墨烯平面内的反演对称性, 进而打开石墨烯的能隙. 同时氮化硼衬底会使石墨烯发生晶格重构, 可进一步提高能隙值^[88]. 目前报道的石墨烯/氮化硼莫尔超晶格(石墨烯与氮化硼晶向一致)的能隙, 不同的实验测量得到的数值并不一致. 例如, 输运测试中观察到的带隙值范围是10—30 meV^[82]; 磁光谱测得带隙值约为40 meV^[89]; 角分辨光电子能谱测得的带隙值是160 meV^[90]; 光电流谱测得的带隙值是14 meV^[91]; 而理论预测的能隙值范围是2—14 meV^[77,78,92]. 更有趣的是, 当魔角双层石墨烯与氮化硼对齐时, 体系中出现本征QAH效应^[13]. 这里可能是因为强相互作用将电子极化为具有单自旋/能谷分辨的孤立小能带, 每个小能带都有非零的陈数, 从而驱动了QAH效应的产生^[93,94]. 同时, 氮化硼衬底会影响魔角石墨烯中的超导态^[21,95]. 综上所述, 当石墨烯与氮化硼晶向一致时, 衬底对石墨烯体系的晶体结构和电子性质产生了显著的影响.

石墨烯与氮化硼的晶格结构存在1.7%的晶格失配, 可形成的最大莫尔条纹长度约为14 nm, 对应的晶胞包含12322个原子. 大多数理论研究都是从低能有效模型出发, 将衬底的效应以莫尔超晶格势的形式加入石墨烯的哈密顿量中, 忽略了晶格弛豫对体系的影响^[77,78], 或只考虑面外原子位移的变化^[88]. 低能有效模型可以有效解释部分实验现象, 但是模型过于简单, 不能考虑原子真实的形变. 接下来, 我们结合半经典分子动力学和紧束缚模型, 系统地阐述氮化硼衬底对单层石墨烯和魔角石墨烯晶格结构和电子性质的影响.

4.1   单层石墨烯-hBN异质结构

2015年, Slotman等^[41]利用分子动力学模拟软件LAMMPS^[51]对超晶格进行晶格弛豫, 石墨烯层内相互作用采用经验键序作用势^[96], 层间相互作用采用Kolmogorov-Crespi势. 对于石墨烯中碳原子与hBN中氮原子和硼原子间相互作用, 仍用Kolmogorov-Crespi势, 但C—B和C—N之间相互作用分别是C—C的0.6倍和2倍^[41]. 为了模拟块体衬底, hBN原子结构保持不变. 晶格弛豫会在体系中引入形变标势${V}_{D}$和矢势A, 这些效应都可以加入石墨烯紧束缚模型的哈密顿量中^[16]. 矢势A通过Slater-Koster近似改变原子间的跃迁振幅, 标势${V}_{D}$以占位能的形式加入, 其表达式是

$${V}_{D}\left({\boldsymbol{R}}_{\mathrm{i}}\right)={g}_{1}\frac{S\left({\boldsymbol{R}}_{\mathrm{i}}\right)-{S}_{0}}{{S}_{0}} \text{, }$$

(33)

其中, $S\left({\boldsymbol{R}}_{\mathrm{i}}\right)$是晶格弛豫后${\boldsymbol{R}}_{\mathrm{i}}$位置处原子的有效面积, ${S}_{0}=3\sqrt{3}{d}^{2}/4$是非弛豫结构中碳原子的有效面积, ${g}_{1}=4$ eV是石墨烯中屏蔽引入的形变势^[97].

图7(a)和图7(b)展示了转角为θ = 0°的石墨烯与氮化硼超晶格体系在结构弛豫后修正的紧束缚模型参数. 随后, 利用紧束缚传播方法计算体系的DOS和准本征态在实空间的分布情况. 在此部分计算中, 体系哈密顿矩阵(20)式只考虑碳原子的跃迁能和占位能. 如图7(c)和图7(d)所示, 对于转角θ = 0°的体系, 晶格弛豫后, 导带和价带都出现了超晶格导致的新的狄拉克锥, 这与实验结果相符合^{[80-82,98,99]}. 超晶格狄拉克点的位置取决于莫尔超晶格的倒格矢$\boldsymbol{G}$, 由${E}_{D}=\pm \hbar {v}_{\mathrm{F}}\left|\boldsymbol{G}\right|/2$给出, 其中${v}_{\mathrm{F}}$为费米速度^[100,101]. 因为倒格矢长度取决于θ, 超晶格狄拉克点的位置会随转角θ而变化(图7(d)). 结构弛豫效应会随θ的增加而减小, 这意味着对于较大θ, 超晶格狄拉克锥会消失.

图7(e)展示了不同能量准本征态在实空间的振幅分布. 由于hBN衬底打破了石墨烯子晶格的子晶格对称性, 我们将两套子晶格的振幅分布分别绘制. 如图7(e)所示, 当准本征态远离费米能量时, 出现了空间局域性. 对于能量接近费米能量的准本征态, 空间幅度的分布差异很小. 而对于能量更接近超晶格狄拉克点的准本征态, 空间分布表现出清晰的莫尔条纹. DOS中超晶格狄拉克点的出现和准本征态中的空间局域性特征表明单层石墨烯/hBN的电子结构受到了结构弛豫的强烈影响. 石墨烯在hBN衬底上的直流电导率的实验测量结果^[80-82]表明, 在空穴侧和电子侧的超晶格狄拉克点处, 电导率出现不对称的下降. 空穴侧电导率的降低更为显著, 其值甚至低于狄拉克点的最小电导率^[80].

4.2   TBG-hBN莫尔超晶格

为进一步研究氮化硼衬底对TBG晶体结构和电学性质的影响, 我们考虑TBG/hBN和hBN/TBG/hBN两大类体系, 并采用类似于石墨烯/hBN的势参数对体系进行晶格弛豫^[102]. 同时, 在紧束缚模型(21)式中加入氮化硼贡献的跃迁能和占位能, 其中硼原子、氮原子和碳原子的占位能分别是${V}_{\mathrm{B}}=3.34$ eV, ${V}_{\mathrm{N}}=-1.4$ eV和${V}_{\mathrm{C}}=0$, 跃迁能由Slater-Koster近似确定^[77].

对于TBG/hBN结构, 如图8所示, 当下层石墨烯与氮化硼衬底夹角小于2°时, 衬底会显著改变体系的晶体结构和电子性质^[102]. 随着转角${\theta }_{\mathrm{b}\mathrm{o}\mathrm{t}}$增大, AAA堆垛区域的涡旋状由三角形变为圆形. 晶格弛豫使下层石墨烯在ABC区域产生最大值约为60 meV的形变势能, 在AAA区域产生最大值为18 T的赝磁场. 而上层石墨烯的形变势能和赝磁场值都要小一个数量级. 因此, 在低能有效模型中通过在相邻石墨烯中引入有效周期势来考虑hBN衬底对TBG的影响是合理的^[103-105]. TBG/hBN哈密顿矩阵的占位能来自于两部分的贡献: 第一部分是氮原子和硼原子的占位能之差, 这种贡献打破石墨烯平面内的反演对称性, 是体系能隙的主要来源; 第二部分是由晶格弛豫引起的形变势能, 导致相邻子晶格的势能差, 会在体系局部产生能隙, 对体系能带有微调作用. 作用于K谷的赝磁场方向与K' 谷相反, 因此赝磁场没有打破体系的时间反演对称性.

如图8(d)—(f), 随着${\theta }_{\mathrm{b}\mathrm{o}\mathrm{t}}$的增大, hBN衬底对TBG能带影响减小, 在${\theta }_{\mathrm{b}\mathrm{o}\mathrm{t}}$接近3°时, 平带之间仍打开能隙. 当${\theta }_{\mathrm{b}\mathrm{o}\mathrm{t}} = 0.53^{\circ}$时, 能隙值约为25—30 meV. 对于无衬底的TBG, 层自由度、自旋和谷自由度没有相互耦合, 费米面附近的平带具有八重简并. 在TBG/hBN体系中, hBN衬底对下层石墨烯的影响更大, 打破TBG的层简并度. 同时, 上层K谷映射到下层K' 谷, 如图8(d)所示, 带边的平带出现分裂. 不同谷的能带可通过谷算子^[106]来分辨:

$${\hat{V}}_{z}=\frac{i}{3\sqrt{3}}\sum _{ \langle i,j\rangle ,s}{\eta }_{ij}{\sigma }_{z}^{ij}{C}_{i,s}^{†}{C}_{j,s} \text{, }$$

(34)

其中i, j表示对次近邻位点求和; ${\eta }_{ij}=\pm 1$分别表示顺时针和逆时针跃迁; ${\sigma }_{z}^{ij}$是定义在子晶格自由度的泡利矩阵. 谷算子的期望值范围从1到 –1, 分别对应K和K' 谷.

接着在TBG/hBN异质结表面再加另外一层氮化硼形成三明治结构, 如图9(a)所示, hBN/TBG/hBN体系中引入一个新的自由度${\theta }_{\mathrm{t}\mathrm{o}\mathrm{p}}$^[102]. 当${\theta }_{\mathrm{t}\mathrm{o}\mathrm{p}}={\theta }_{\mathrm{b}\mathrm{o}\mathrm{t}}$时, 上层和下层hBN晶格结构对齐, 体系中的平带打开约50 meV能隙, 对应赝磁场最大值约为36 T, 同时, 体系恢复八重简并. 由图9可以明显看出, 三明治结构的电学性质对hBN旋转角度高度敏感.

5.   总结与展望

目前, 石墨烯莫尔超晶格的实验和理论研究已取得一系列进展, 出现了很多新奇现象, 但还需要大量工作来揭示其背后的物理机制. 例如在实验和理论上, 需要更加深入探究魔角石墨烯体系中超导和关联绝缘态的机制. 高阶魔角体系中是否存在类似于第一魔角体系的物理现象? 旋转多层石墨烯体系中是否有更丰富、更新的新奇性质? 此外, 其他二维材料(如过渡金属硫族化合物、氮化硼和黑磷等^[2]) 所组成的莫尔超晶格引起了广泛关注. 如实验和理论已在过渡金属硫族化合物莫尔超晶格中发现了平坦能带^[107-109], 并观测到关联绝缘态^[110]、QAH效应^[111]、连续金属绝缘体转变^[112]、奇异金属态^[113]和铁电畴^[114]等. 相对于石墨烯构成的异质结, 其他二维材料的莫尔超晶格研究正处于起步阶段, 期待出现更多新奇现象. 与此同时, 晶格弛豫和衬底效应将会对其他二维材料的莫尔超晶格产生重要影响, 需要更多实验和理论研究去阐明相关的物理类机制.