本文主要引用自中科大郑奇靖博士的博文Plotly: Tight-binding Model for GrapheneCC BY 4.0)。同时此文可以作为教科书:Alexander Altland, Condensed Matter Field Theory(2nd edition)第二章第二节Applications of second quantization->Tight-binding systems的详细释读。本文不包含电子间相互作用部分。

紧束缚模型与Wannier函数介绍

讨论固体物理,非常重要的条件就是晶格的周期性。在讨论后续的模型之前,我们先已经知道,周期性的哈密顿量得到的波函数应当为布洛赫波(Bloch waves),也是近自由电子模型、紧束缚模型共同的出发点:

ψ𝐤n(𝐫)=𝚎𝚒𝐤𝐫u𝐤n(𝐫)\psi_{𝐤n}(𝐫)=𝚎^{𝚒𝐤\cdot 𝐫}u_{𝐤n}(𝐫)

其中𝐤𝐤在第一布里渊区中,nn是能带的标记,u𝐤n(𝐫)u_{𝐤n}(𝐫)拥有晶格周期性:

u𝐤n(𝐫+任意晶格矢量𝐑)=u𝐤n(𝐫).u_{𝐤n}(𝐫+\text{任意晶格矢量}𝐑)=u_{𝐤n}(𝐫).

可见布洛赫波弥散与整个晶体的空间,每个布洛赫波都有一个确定的参数𝐤𝐤。如果对于一个高迁移率(high mobility)导电电子的体系,我们常用的近似叫做近自由电子模型,直接认为导电电子所感受到的背景电势(晶体赝势)可以忽略。电子波函数就是平面波,u𝐤n=1u_{𝐤n}=1

但是对于晶格势能严重影响导电电子的体系,比如过度金属氧化物,我们用布洛赫波发展出了另一种从“局域电子”出发的模型:紧束缚模型(tight-binding systems)。

其实很容易理解:平面波弥散于整个空间,但是进行傅里叶变换,也就是不同的平面波积分,我们能够得到δ\delta函数(局域化);类似地,布洛赫波近似于平面波弥散在整个空间,如果我们对其进行傅里叶变换,应该得到一个虽然不是完美δ\delta函数,但也相当局域的函数。这个变换的结果就是Wannier瓦涅尔函数。可以容易想象,这种实质是傅里叶变换的纯数学操作自然不会丢失任何信息,也就是新的函数(Wannier函数)的正交性、完备性和变换前的布洛赫波一样。

数学的语言描述这种变换为:

(1)

ψ𝐑n>1N𝐤B.Z.𝚎𝚒𝐤𝐑ψ𝐤n>\left|\psi_{𝐑n}\right>\equiv\frac{1}{\sqrt N}\sum_𝐤^{\text{B.Z.}}𝚎^{-𝚒𝐤\cdot 𝐑}\left|\psi_{𝐤n}\right>

其中ψ𝐑n>\left|\psi_{𝐑n}\right>就是Wannier态,下标从布洛赫态的𝐤𝐤变成𝐑𝐑,因为是傅里叶变换嘛,很好理解。

反过来的变换为:

(2)

ψ𝐤n>1N𝐑𝚎𝚒𝐤𝐑ψ𝐑n>\left|\psi_{𝐤n}\right>\equiv\frac{1}{\sqrt N}\sum_𝐑𝚎^{𝚒𝐤\cdot 𝐑}\left|\psi_{𝐑n}\right>

需要注意,经过傅里叶变换,Wannier态在𝐤𝐤空间上是弥散的(因而在坐标空间局域),而且Wannier态不再是哈密顿量的本征态(布洛赫波是,相加后就不是了)。

还有一个可以想象的极限状态:如果我们的晶体中原子距离非常遥远以至于可以看作孤立原子,那么局域的Wannier函数应当会退化为原子中的电子波函数。


上面说了个概念。下面的内容我们要玩真的,来看Wannier函数下,我们如何计算体系的哈密顿量。数学的手段是二次量子化的操作,不太清楚的可以直接看Condensed Matter Field Theory(2nd edition)的第二章第一节。

晶体中单电子的Hamilton量为:

H^0=𝐩^22m+V(𝐫),\hat{H}_0=\frac{\hat{𝐩}^2}{2m}+V(𝐫),

在二次量子化的表象下,波函数成为了粒子数升降算符(aσ(𝐫)a^{\dagger}_\sigma(𝐫)),也要进入哈密顿量中:

(3)

H^0=d𝐫aσ(𝐫)[𝐩^22m+V(𝐫)]aσ(𝐫),\hat{H}_0=\int\mathrm{d}\,𝐫\,a^{\dagger}_\sigma(𝐫)\left[\frac{\hat{𝐩}^2}{2m}+V(𝐫)\right]a_\sigma(𝐫),

其中σ\sigma代表电子的自旋,此文用不到电子自旋,但是为了理论的完备,还是写出来。

在实空间中求哈密顿量本征态,得到的就是上面所示的布洛赫态ψ𝐤n>\left|\psi_{𝐤n}\right>,按照二次量子化,我们定义产生算符a𝐤σa^\dagger_{𝐤\sigma}(是aσ(𝐫)a^{\dagger}_\sigma(𝐫)的一个实例)使得

ψ𝐤n>=a𝐤σ0>,\left|\psi_{𝐤n}\right> =a^\dagger_{𝐤\sigma}\left|0\right>,

同时方便定义产生算符的共轭算符:湮灭算符a𝐤σa_{𝐤\sigma}。在这种二次量子化的表象下,Hamilton量被写为(所有单个电子求和)

(4)

H^0=𝐤ϵka𝐤σa𝐤σ.\hat{H}_0=\sum_{𝐤}\epsilon_ka^\dagger_{𝐤\sigma}a_{𝐤\sigma}.

这里已经用到了性质:布洛赫态是该哈密顿量本征态,因此有确定的能量ϵ𝐤\epsilon_𝐤

下面就是表象的变换,从布洛赫转换为Wannier,表象变换操作是大家在量子力学习题中应该见过多次的常用方法(运用完备性I=ψ𝐑n><ψ𝐑nI=\left|\psi_{𝐑n}\right>\left<\psi_{𝐑n}\right|),下面直接给出结果,不清楚的看Condensed Matter Field Theory书公式(2.23)前后:

(5)

a𝐤σ=1N𝐑𝚎𝚒𝐤𝐑iaiσ,aiσ=1N𝐤B.Z.𝚎𝚒𝐤𝐑ia𝐤σa^\dagger_{𝐤\sigma}=\frac{1}{\sqrt N}\sum_𝐑𝚎^{𝚒𝐤\cdot 𝐑_i}a^\dagger_{i\sigma}, \quad a^\dagger_{i\sigma}=\frac{1}{\sqrt N}\sum_𝐤^{\text{B.Z.}}𝚎^{-𝚒𝐤\cdot 𝐑_i}a^\dagger_{𝐤\sigma}

其中aiσa^\dagger_{i\sigma}就是我们需要的Wannier态的产生算符ψ𝐑in>=aiσ0>\left|\psi_{𝐑_in}\right>=a^\dagger_{i\sigma}\left|0\right>ii下标其实是简写了的𝐑i𝐑_i。将(5)代入(4)完成变换:

(6)

H^0=1Nii𝐤𝚎𝚒𝐤(𝐑i𝐑i)ϵkaiσaiσiiaiσtiiaiσ,\hat{H}_0=\frac{1}{N}\sum_{ii'}\sum_{𝐤}𝚎^{𝚒𝐤\cdot (𝐑_i-𝐑_{i'})}\epsilon_ka^\dagger_{i\sigma}a_{i'\sigma}\equiv\sum_{ii'}a^\dagger_{i\sigma}t_{ii'}a_{i'\sigma},

其中tii=N1ii𝐤𝚎𝚒𝐤(𝐑i𝐑i)ϵkt_{ii'}=N^{-1}\sum_{ii'}\sum_{𝐤}𝚎^{𝚒𝐤\cdot (𝐑_i-𝐑_{i'})}\epsilon_k。这就是目前不考虑其它相互作用,紧束缚模型得到的哈密顿量。

在我们使用紧束缚模型时,常常用到的近似是:tiit_{ii'}仅存在于最近邻原子𝐑i,𝐑i𝐑_i,𝐑_{i'}之间(因为Wannier函数局域,忽略更远近邻的影响)。

简单讨论一下紧束缚模型(6)公式的含义:

  1. 对于i=ii=i'也就对同一位置的电子,哈密顿量的对角位,tii=𝐤ϵkNϵt_{ii}=\frac{\sum_𝐤\epsilon_k}{N}\equiv \epsilon,理解上大致对应一个局域电子本身的能量;

  2. iii\neq i'最近邻电子之间,哈密顿量的非对角位,计算得到的tii=N1ii𝐤𝚎𝚒𝐤(𝐑i𝐑i)ϵkt_{ii'}=N^{-1}\sum_{ii'}\sum_{𝐤}𝚎^{𝚒𝐤\cdot (𝐑_i-𝐑_{i'})}\epsilon_k一般为一个负值(在紧束缚近似下,tii=t,t>0t_{ii'}=-t,\,t>0),理解为局域的电子向其最近邻位置有一定跳跃的可能,这种跳跃增加了局域电子的活动空间,从而降低了系统的能量。
    这里大家可能有点迷糊:tii=N1ii𝐤𝚎𝚒𝐤(𝐑i𝐑i)ϵkt_{ii'}=N^{-1}\sum_{ii'}\sum_{𝐤}𝚎^{𝚒𝐤\cdot (𝐑_i-𝐑_{i'})}\epsilon_k这个公式看起来非常的不直观,不知道它代表了什么样的物理意义,怎么就莫名奇妙约等于一个负的常数t-t了呢?如果我们换一种写法,不是用(4)的从布洛赫态展开哈密顿量,而是用(3)的更为一般的实坐标态下的哈密顿量写,就会更清楚。我们使用aσ(𝐫)aiσa^{\dagger}_\sigma(𝐫)\rightarrow a^\dagger_{i\sigma}的表象变换(推导过程与(5)原理一致):

    aσ(𝐫)=iψ𝐑i(𝐫)aiσa^{\dagger}_\sigma(𝐫)=\sum_i \psi_{𝐑_i}^*(𝐫)a^\dagger_{i\sigma}

    将(3)变为:

    H^0=iid𝐫ψ𝐑i(𝐫)aiσ[𝐩^22m+V(𝐫)]ψ𝐑i(𝐫)aiσ,\hat{H}_0=\sum_{ii'} \int \mathrm{d}\,𝐫\,\psi_{𝐑_i}^*(𝐫)a^\dagger_{i\sigma}\left[\frac{\hat{𝐩}^2}{2m}+V(𝐫)\right]\psi_{𝐑_{i'}}(𝐫)a_{i'\sigma},

    也就是:tii=<ψ𝐑iH^0ψ𝐑i>=tiit_{ii'}=\left<\psi_{𝐑_i}\right|\hat{H}_0\left|\psi_{𝐑_{i'}}\right>=t_{i'i}^*。这种量子力学中的矩阵元大家应该熟悉,体现的是最近邻两个Wannier函数的重合部分对能量的贡献,也应该理解为什么我们自然而然将它近似为常数。

紧束缚近似写为:

(7)

tii={ϵi=itii nearest neighbours0otherwiset_{ii'}=\left\{\begin{array}{ll}\epsilon & i=i'\\-t & i\,i' \text{ nearest neighbours}\\0 & \text{otherwise}\end{array}\right.

提醒一下,到这里我们讨论的全是哈密顿量H^0\hat{H}_0,本质是单电子哈密顿量,而没有考虑电子-电子的库伦相互作用。讨论相互作用的影响的紧束缚模型不在本文讨论范围内。

石墨烯结构

Honeycomb lattice of graphene where different colors are used to denote the two sublattices. The basis vectors of the unit cell are shown with black arrows.

This figure is generated by TikZ/LaTeX(From郑奇靖博文

石墨烯的蜂窝结构可谓非常经典了。只是需要说明下它的晶胞:原胞呈现菱形(灰色虚线),每个原胞内有两个碳原子A和B(红,蓝);每个原胞内有三条不同的最近邻连接线。𝐚1,𝐚2𝐚_1,𝐚_2为两个原胞矢量。

晶格矢量𝐑(x,y)=x𝐚1+y𝐚2𝐑(x,y)=x𝐚_1+y𝐚_2,其中x,yx,y为整数。

这种结构的成因:在石墨烯的平面内,sp2\text{sp}^2电子杂化,三条夹角120°的成键方向,产生每个碳原子连接三个碳原子的蜂窝结构。但是还剩下一个pz\text{p}_z轨道垂直于平面,各个原子的pz\text{p}_z轨道电子互相有一点的影响,产生π\pi键。而我们的紧束缚模型,正好可以很好的用于解释这个π\pi键电子的行为!

每个原胞中三条不同的最近邻线,以图中标注原子A起始为例,假设AB原子所在的晶胞位置为𝐑(0,0)𝐑(0,0)

  1. 向右连接本晶胞内𝐑(0,0)𝐑(0,0)的B,
  2. 向左上连接晶胞𝐑(0,1)𝐑(0,-1)中的B,
  3. 向左下连接晶胞𝐑(1,0)𝐑(-1,0)中的B。

这三个最近邻的作用应当是完全相同的,因为晶格的对称性。(旋转120°,240°完全一样)

其实你或者也可以从原子B起始,也写出3条最近邻线,这其实就是上面三条线晶格平移的结果:

  1. 向左连接本晶胞内𝐑(0,0)𝐑(0,0)的A,
  2. 向右下连接晶胞𝐑(0,1)𝐑(0,1)中的A,
  3. 向右上连接晶胞𝐑(1,0)𝐑(1,0)中的A。

从A开始或从B开始,二者等价。

紧束缚近似应用于石墨烯𝜋键

我们将(6)(7)应用于刚刚介绍的石墨烯,也就是计算石墨烯体系的能量。更为确切地说,是计算石墨烯(紧束缚近似下)的能带,即能量EE关于动量空间𝐤𝐤的函数。我们知道,其实Wannier态其实在动量空间是弥散的,所以最后的结果,我们不应该保留局域的Wannier态,而是仍然要傅里叶变换回布洛赫波的表象(拥有确定的𝐤𝐤可以作为自变量)。

注意一个小问题,石墨烯的原胞内有两个碳原子,需要区分讨论,因此(7)中本来只有一个的aiσa^\dagger_{i\sigma}这里也要分化出两个ai1σa^\dagger_{i1\sigma}ai2σa^\dagger_{i2\sigma},下标1,2分别代表同一个原胞ii中的原子A和B。两者对应的Wannier函数(各自的局域电子)自然也要求互相正交。(也就是说,对于每个晶格,并非只有一个Wannier态,而是2个原子分别有一个)Wannier函数变多了一倍,布洛赫波的个数自然也要扩大一倍:a𝐤1σa^\dagger_{𝐤1\sigma}a𝐤2σa^\dagger_{𝐤2\sigma}

由于这个变化,我们的哈密顿量的写法相对(6)要有小小的修改,看看ai1σa^\dagger_{i1\sigma}ai2σa^\dagger_{i2\sigma}怎么加入进去。整个原胞内的哈密顿量有2个原子+3条最近邻组合。首先看2个原子上的电子本身都拥有能量,对应着哈密顿量的对角项,也就是:

(8)

H^diag=iϵ(ai1σai1σ+ai2σai2σ)=ϵi𝐤𝐤1N𝚎𝚒𝐤𝐑ia𝐤1σ𝚎𝚒𝐤𝐑ia𝐤1σ+(12)=ϵ𝐤𝐤i1N𝚎𝚒(𝐤𝐤)𝐑ia𝐤1σa𝐤1σ+(12)=ϵ𝐤𝐤δ(𝐤𝐤)a𝐤1σa𝐤1σ+(12)=ϵ𝐤a𝐤1σa𝐤1σ+(12)\begin{aligned} \hat{H}_\text{diag}&=\sum_i\epsilon(a^\dagger_{i1\sigma}a_{i1\sigma}+a^\dagger_{i2\sigma}a_{i2\sigma})\\ &=\epsilon\sum_i\sum_{𝐤𝐤'}\frac{1}{N}𝚎^{-𝚒𝐤\cdot 𝐑_i}a^\dagger_{𝐤1\sigma}𝚎^{𝚒𝐤'\cdot 𝐑_i}a_{𝐤'1\sigma} + (1\rightarrow 2)\\ &=\epsilon\sum_{𝐤𝐤'}\sum_i\frac{1}{N}𝚎^{𝚒(𝐤'-𝐤)\cdot 𝐑_i}a^\dagger_{𝐤1\sigma}a_{𝐤'1\sigma} + (1\rightarrow 2)\\ &=\epsilon\sum_{𝐤𝐤'}\delta(𝐤'-𝐤)a^\dagger_{𝐤1\sigma}a_{𝐤'1\sigma} + (1\rightarrow 2)\\ &=\epsilon\sum_{𝐤}a^\dagger_{𝐤1\sigma}a_{𝐤1\sigma} + (1\rightarrow 2) \end{aligned}

其中第二行用到了(5)Wannier态和布洛赫态的转换;第三行用了i1N𝚎𝚒(𝐤𝐤)𝐑i=δ(𝐤𝐤)\sum_i\frac{1}{N}𝚎^{𝚒(𝐤'-𝐤)\cdot 𝐑_i}=\delta(𝐤'-𝐤)

剩下的还有3个最近邻作用,也就是最近邻hopping使得能量降低,对应的是哈密顿量中的非对角项。需要说明的是,虽说是3个最近邻作用,但是书写时我们经常把同一个作用的两写法A->B,与B->A都写出来,二者各乘0.5再加起来,这样写最为对称简洁。所以3个作用我们会写出对称的6项。

(9)

H^offdiag=it(ai1σai2σ+ai1σa(𝐑i+𝐑(0,1))2σ+ai1σa(𝐑i+𝐑(1,0))2σ)+it(ai2σai1σ+ai2σa(𝐑i+𝐑(0,1))1σ+ai2σa(𝐑i+𝐑(1,0))1σ)=ti𝐤𝐤1N[𝚎𝚒(𝐤𝐤)𝐑i+𝚎𝚒[𝐤(𝐑i𝐚2)𝐤𝐑i]+𝚎𝚒[𝐤(𝐑i𝐚1)𝐤𝐑i]]a𝐤1σa𝐤2σ+h.c.=t𝐤𝐤i1N𝚎𝚒(𝐤𝐤)𝐑i(1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚1)a𝐤1σa𝐤2σ+h.c.=t𝐤𝐤δ(𝐤𝐤)(1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚1)a𝐤1σa𝐤2σ+h.c.=t𝐤(1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚1)a𝐤1σa𝐤2σ+h.c.\begin{aligned} \hat{H}_\text{offdiag}=& \sum_i -t(a^\dagger_{i1\sigma}a_{i2\sigma} +a^\dagger_{i1\sigma}a_{(𝐑_i+𝐑(0,-1))2\sigma} +a^\dagger_{i1\sigma}a_{(𝐑_i+𝐑(-1,0))2\sigma} ) \\ &+\sum_i -t(a^\dagger_{i2\sigma}a_{i1\sigma} +a^\dagger_{i2\sigma}a_{(𝐑_i+𝐑(0,1))1\sigma} +a^\dagger_{i2\sigma}a_{(𝐑_i+𝐑(1,0))1\sigma} ) \\ =& -t\sum_i\sum_{𝐤𝐤'}\frac{1}{N}\left[𝚎^{𝚒(𝐤'-𝐤)\cdot 𝐑_i}+𝚎^{𝚒[𝐤'\cdot(𝐑_i-𝐚_2)-𝐤\cdot 𝐑_i]}+𝚎^{𝚒[𝐤'\cdot(𝐑_i-𝐚_1)-𝐤\cdot 𝐑_i]}\right]a^\dagger_{𝐤1\sigma}a_{𝐤'2\sigma} +\text{h.c.}\\ =& -t\sum_{𝐤𝐤'}\sum_i\frac{1}{N}𝚎^{𝚒(𝐤'-𝐤)\cdot 𝐑_i}\left(1+𝚎^{-𝚒𝐤'\cdot𝐚_2}+𝚎^{-𝚒𝐤'\cdot𝐚_1}\right)a^\dagger_{𝐤1\sigma}a_{𝐤'2\sigma} +\text{h.c.}\\ =& -t\sum_{𝐤𝐤'}\delta(𝐤'-𝐤)\left(1+𝚎^{-𝚒𝐤'\cdot𝐚_2}+𝚎^{-𝚒𝐤'\cdot𝐚_1}\right)a^\dagger_{𝐤1\sigma}a_{𝐤'2\sigma} +\text{h.c.}\\ =& -t\sum_{𝐤}\left(1+𝚎^{-𝚒𝐤\cdot𝐚_2}+𝚎^{-𝚒𝐤\cdot𝐚_1}\right)a^\dagger_{𝐤1\sigma}a_{𝐤2\sigma} +\text{h.c.} \end{aligned}

其中h.c.\text{h.c.}代表了Hermitian conjugate(厄米共轭)。

综合(8)(9),我们可以把总哈密顿量对于每一个𝐤𝐤拆开写成H^(𝐤)\hat{H}(𝐤),并写成矩阵形式:

(10)

H^(𝐤)=[a𝐤1σa𝐤2σ][ϵt(1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚1)t(1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚1)ϵ][a𝐤1σa𝐤2σ]\hat{H}(𝐤)=\begin{bmatrix}a^\dagger_{𝐤1\sigma} & a^\dagger_{𝐤2\sigma}\end{bmatrix} \begin{bmatrix}\epsilon & -t\left(1+𝚎^{-𝚒𝐤\cdot𝐚_2}+𝚎^{-𝚒𝐤\cdot𝐚_1}\right)\\ -t\left(1+𝚎^{𝚒𝐤\cdot𝐚_2}+𝚎^{𝚒𝐤\cdot𝐚_1}\right) & \epsilon \end{bmatrix} \begin{bmatrix}a_{𝐤1\sigma} \\ a_{𝐤2\sigma}\end{bmatrix}

对其进行对角化,而本征值就是能量。求本征值的操作过于基础,此处忽略。直接得到结果:

E(𝐤)=ϵ±t1+𝚎𝚒𝐤𝐚1+𝚎𝚒𝐤𝐚1+𝚎𝚒𝐤𝐚2+𝚎𝚒𝐤𝐚2+(𝚎𝚒𝐤𝐚1+𝚎𝚒𝐤𝐚2)(𝚎𝚒𝐤𝐚1+𝚎𝚒𝐤𝐚2)=ϵ±t3+2cos(𝐤𝐚1)+2cos(𝐤𝐚2)+2cos(𝐤(𝐚1𝐚2)).\begin{aligned} E(𝐤) &= \epsilon \pm t\sqrt{ 1 + 𝚎^{-𝚒𝐤\cdot𝐚_1} + 𝚎^{ 𝚒𝐤\cdot𝐚_1} + 𝚎^{-𝚒𝐤\cdot𝐚_2} + 𝚎^{ 𝚒𝐤\cdot𝐚_2} + (𝚎^{-𝚒𝐤\cdot𝐚_1} + 𝚎^{-𝚒𝐤\cdot𝐚_2}) (𝚎^{ 𝚒𝐤\cdot𝐚_1} + 𝚎^{ 𝚒𝐤\cdot𝐚_2}) } \\ &= \epsilon \pm t\sqrt{ 3 + 2\cos(𝐤\cdot𝐚_1) + 2\cos(𝐤\cdot𝐚_2) + 2\cos(𝐤\cdot(𝐚_1-𝐚_2)) }. \end{aligned}

再为了方便后面绘图,代入xy坐标下的原胞矢量:

𝐚1=a2(3,3)𝐚2=a2(3,3),\begin{aligned} 𝐚_1 &= {a\over2}(3, \sqrt{3}) \\ 𝐚_2 &= {a\over2}(3,-\sqrt{3}), \end{aligned}

其中aa就是晶格常数。我们就可以把𝐤𝐤也写成分量,则有:

(11)

E(kx,ky)=±t3+2cos(3kya)+4cos(32kxa)cos(32kya).E(k_x, k_y) = \pm t \sqrt{ 3 + 2 \cos(\sqrt{3}k_y a) + 4 \cos({3\over2}k_x a)\cos({\sqrt{3}\over2}k_y a) }.

我们忽略了能量基准ϵ\epsilon

下面便是使用Plotly对能带的绘图: