Tight Binding Approximation（TBA）

緊束縛模型計算石墨烯能帶

Lattice vectors：

a_{1}=\frac{3a}{2}i+\frac{\sqrt{3}a}{2}j ,\quad a_{2}=\frac{3a}{2}i-\frac{\sqrt{3}a}{2}j

where

a

is the lattice constant。 The corresponding reciprocal vectors：

b_{1}=\frac{2\pi}{3a}\hat{k}_{x}+\frac{2\sqrt{3}\pi}{3a}\hat{k}_{y},\quad b_{2}=\frac{2\pi}{3a}\hat{k}_{x}-\frac{2\sqrt{3}\pi}{3a}\hat{k}_{y}

The lattice has two non-equivalent sites （marked by A and B） with the nearest hopping vectors：

\rho_{1}=\frac{a}{2}i+\frac{\sqrt{3}a}{2}j,\quad \rho_{2}=\frac{a}{2}i-\frac{\sqrt{3}a}{2}j,\quad \rho_{3}=-ai

The Bloch states are the normalized orthogonal basis and can be written in the form of Fourier transformation of Wannier basis ：

\psi_{nk}(r)=\frac{1}{\sqrt{N}}\sum_{l}e^{ik\cdot l}a_{n}(r,l)\quad

Wannier function

a_{n}(r,l)

is a function of

(r-l)

。 Here we only consider single band model thus the band index n will be omitted。 Considering two inequivalent A，B sites， we can write the Bloch states as two sets：

\varphi_{1}=\frac{1}{\sqrt{N}}\sum_{l}e^{ik\cdot R^{A}_{l}}\phi(r-R^{A}_{l}),\quad \varphi_{2}=\frac{1}{\sqrt{N}}\sum_{l

where N is the number of unit cells，

R^{A/B}_{j}

is the position vector of

A/B

atom in

Jth

unit cell。

\varphi_{1}

and

\varphi_{2}

are the Bloch states of A and B sites。 Because of the similar localization of wavefunctions we usually use the actual atom‘s orbital wavefunctions

\phi

to substitute the Wannier function

a_{n}(r-l)

。 Usually， this is a good approximation。 Except for the

sp_{2}

hybrid orbitals of

\sigma

bonds， every carbon atom has an extra

p_{z}

orbital for the

\pi

bond。 So here we use

2p_{z}

orbital wavefunctions of the carbon atom。 The wavefunction of the whole system can be written in the linear combination of atomic orbitals （LCAO）：

\psi(r)=c_{1}\varphi_{1}+c_{2}\varphi_{2}=\frac{1}{\sqrt{N}}\sum_{l,l

To solve the Schordinger equation：

H|\psi\rangle=E|\psi\rangle

we use the

|\varphi_{1}\rangle

and

|\varphi_{2}\rangle

to multiply the equation from the left， we have：

\\ \langle\varphi_{2}|H|\psi\rangle=c_{1}\langle\varphi_{2}|H|\varphi_{1}\rangle+c_{2}\langle\varphi_{2}|H|\varphi_{2}\rangle=E[c_{1}\langle\varphi_{2}|\varphi_{1}\rangle+\langle\varphi_{2}|\varphi_{2}\rangle]

\\ \langle\varphi_{1}|H|\psi\rangle=c_{1}\langle\varphi_{1}|H|\varphi_{1}\rangle+c_{2}\langle\varphi_{1}|H|\varphi_{2}\rangle=E[c_{1}\langle\varphi_{1}|\varphi_{1}\rangle+\langle\varphi_{1}|\varphi_{2}\rangle]

Assume that the overlap of different sites’ orbitals is zero：

\langle\varphi_{i}|\varphi_{j}\rangle=\delta_{ij}

。 So we have：

\\ c_{1}H_{11}+c_{2}H_{12}=c_{1}E, \quad c_{1}H_{21}+c_{2}H_{22}=c_{2}E

where

H_{ij}=<\varphi_{i}|H|\varphi_{j}>

and apparently

H_{12}=H_{21}^{*}

。 In order to get a nontrivial solution for the coefficients， the

H

and

E

must satisfy the secular equation：

\\ \begin{vmatrix} H_{11}-E \quad H_{12}\\ H_{21} \quad H_{22}-E \end{vmatrix} =0

If we look at the determinant we will find that it‘s equal to

det(H-E\cdot I)

which means we are solving the eigenvalues of

H

where

\\ H=\begin{bmatrix} H_{11} \quad H_{12}\\ H_{21} \quad H_{22} \end{bmatrix}

The

H

is written in two orthogonal Bloch basis ：

|\varphi_{1}\rangle

and

|\varphi_{2}\rangle

。 Therefore we have：

\\ E=\frac{1}{2} [H_{11}+H_{22}\pm \sqrt{[(H_{11}-H_{22})^{2}+4|H_{12}|^{2}]} ]

Next we compute the matrix elements：

H_{11}=\langle\varphi_{1}|H|\varphi_{1}\rangle=\frac{1}{N}\sum_{ll

since

\phi

is orthogonal and is the local orbital energy eigen-wavefunction。 Similarly：

H_{22}=\langle\varphi_{2}|H|\varphi_{2}\rangle=\frac{1}{N}\sum_{ll

H_{12}=\langle\varphi_{1}|H|\varphi_{2}\rangle=\frac{1}{N}\sum_{l,l

Where

J=\langle\phi(r-R_{l

with

H_{a}=-\frac{\hbar^{2}}{2m}\triangledown^{2}+v_{a}

。 Namely we separate the Hamiltonian into each individual single atomic Hamiltonian part

H_a

and the rest remaining potential part。

Considering that ：

\\ \langle\phi(r-R_{l

\\ \quad J=<\phi(r-R_{l

We call

J

as the overlap integral。 Assume that

J

is the same for all nearest interactions。 And if we only consider nearest interaction（ only when

|R^A_{l}-R^B_{l

with

\rho_i

being nearest hopping vector， can we get nonzero

J

）：

H_{12}=\frac{1}{N}\sum_{l,l

So we get：

E=\varepsilon_{p_{z}}\pm J\sqrt{3+2\cos(\sqrt{3}k_{y} a)+4\cos(\frac{3}{2}k_{x}a)\cos(\frac{\sqrt{3}}{2}k_{y}a)}

Set constant term

\varepsilon_{p_{z}} =0

， then we have：

E(k)=\pm J\sqrt{3+2\cos(\sqrt{3}k_{y} a)+4\cos(\frac{3}{2}k_{x}a)\cos(\frac{\sqrt{3}}{2}k_{y}a)} \\

Using mathematica to plot this band structure ：

緊束縛模型計算石墨烯能帶

J=2， kx，ky=［-4*（3）^（-1/2）pi/9a，+4*（3）^（-1/2）pi/9a］

Second Quantization of TBA

Now we repeat this process under the second quantization theory。

A state in Hilbert space can be written in Wannier basis or Bloch basis：

\Psi(r)=\sum_{l}C_{l}a(r-l),\quad \Psi(r)=\sum_{k}C_{k}\psi_{k}(r)

In the summation we omit the band index

n

and spin index

\sigma

。

H=\int \Psi^{\dagger}(r)h(r)\Psi(r)dr,\quad h(r)=-\frac{\hbar^{2}}{2m}\triangledown^{2}+V(r)

Again， if we only consider nearest hopping and use atom’s actual orbital to substitute the Wannier function then we have ：

H=\sum_{l,l

where

J=\langle\phi(r-R_{l}^{B})|V-v_{a}|\phi(r-R_{l}^{A})\rangle

。 Usually we first write

\hat{H}

in the form of

C^{+}_{l}

and

C_{l}

because it‘s easy to understand with clear physical meanings if we take

C^{+}_{l}

and

C_{l}

as production and annihilation operator。 Because of the Fourier transformation relation between

a(r-l)

and

\psi_{k}(r)

， we get：

C_{l}=\frac{1}{\sqrt{N}}\sum_{k}e^{ik\cdot l}C_{k}

We notice that

\sum_{l}C^{\dagger}_{l}C_{l}=\frac{1}{N}\sum_{l}\sum_{k,k

Therefore we have the following two equations：

\sum_{l}C^{\dagger}_{l}C_{l}=\sum_{k}C^{\dagger}_{k}C_{k}

\sum_{l}C^{\dagger}_{l}C_{l+\rho}=\sum_{k}e^{ik\rho}C^{\dagger}_{k}C_{k}

Therefore we can diagonalize the

\hat{H}

in the

C_{k}

basis：

\hat{H}=\sum_{k}(\varepsilon_{0}-J\sum_{\rho}e^{ik\rho})C^{\dagger}_{k}C_{k},\quad \hat{H}=E(k)\hat{n}_{k}

with the particle-number operator

\hat{n}_{k}=\sum_{k}C^{\dagger}_{k}C_{k}

and

E(k)=\varepsilon_{0}-J\sum_{\rho}e^{ik\rho}

。

As for the graphene， we have A&B nonequivalent sites thus we write

\hat{H}

in

a^{\dagger}_{l}

and

b^{\dagger}_{l}

：

\hat{H}=\varepsilon_{A}\sum_{l}a^{\dagger}_{l}a_{l}+\varepsilon_{B}\sum_{l}b^{\dagger}_{l}b_{l}-J\sum_{l,\rho_{i}}(a^{\dagger}_{l}b_{l+\rho_{i}}+b^{\dagger}_{l+\rho_{i}}a_{l})

Here we assume on-site energy

\varepsilon_{A}=\varepsilon_{B}=\epsilon_{pz}

。 Then we do Fourier transformation。

\hat{H}=\varepsilon_{A}\sum_{k}a^{\dagger}_{k}a_{k}+\varepsilon_{B}\sum_{k}b^{\dagger}_{k}b_{k}-J\sum_{k,\rho_{i}}(a^{\dagger}_{k}b_{k}e^{ik\cdot \rho_{i}}+h.c.)\\ =N\varepsilon_{pz}-J\sum_{k,\rho_{i}}(a^{\dagger}_{k}b_{k}e^{ik\cdot \rho_{i}}+h.c.)

when we drop the constant

N\varepsilon_{pz}

term， we can write the

\hat{H}

in the matrix form of Nambu basis which preserves the particle hole symmetry：

\hat{H}=\sum_{k}(a_{k}^\dagger,b_{k}^\dagger) \begin{pmatrix} 0 \quad H_{12}(k) \\H_{21}(k) \quad 0 \end{pmatrix}\begin{pmatrix}a_k\\b_k\end{pmatrix}

where

H_{12}(k)=-\sum_{\rho_{i}}Je^{ik\cdot\rho_{i}},H_{21}(k)=H_{12}(k)^{*}

，

The eigenvalue of matrix is obtained by

{\rm det}(\begin{bmatrix} 0 \quad H_{12}(k)\\ H_{21}(k) \quad 0 \end{bmatrix}-E(k)I)=0

Again if we only consider nearest neighbor hopping we will get：

\\ E(k)=\pm J\sqrt{3+2\cos(\sqrt{3}k_{y} a)+4\cos(\frac{3}{2}k_{x}a)\cos(\frac{\sqrt{3}}{2}k_{y}a)} \\

Dirac Cone and Massless Fermion

According to the rotation symmetry of graphene， we define new x，y axis to let the reciprocal vector in Brillouin zone become like this：

緊束縛模型計算石墨烯能帶

The contour plot of the energy dispersion for the lower band as a function of kx and ky。 The red dash lines mark the first BZ， which is a hexagon。

The the energy dispersion relation

E(k)

is written as ：

E(k)=\pm J\sqrt{3+2\cos(\sqrt{3}k_{x} a)+4\cos(\frac{3}{2}k_{y}a)\cos(\frac{\sqrt{3}}{2}k_{x}a)}

There are two kinds of K points in Brillouin zone ：

K=(\frac{4\sqrt{3}\pi}{9a},0),\quad K

The equivalent K points are connected by a reciprocal vector。 And in the K（k’） points， the

E(K)=0

：

E(K)=\pm J\sqrt{3+2cos(\frac{4\sqrt{3}\pi}{9a}\times\sqrt{3}a)+4cos(\frac{4\sqrt{3}\pi}{9a}\times\frac{\sqrt{3}a}{2})cos(\frac{3}{2}\times0)}

\quad \quad=\pm J \sqrt{3+2\times (-\frac{1}{2})+4\times(-\frac{1}{2})}=0

Near these K（K‘） points， the energy bands are linear functions of lattice momentum k。

Consider a small shift

q

around K points：

k_{new}=K+q

， then we have the expansion ：

E(K+q)=\pm J \sqrt{3+2cos(\sqrt{3}a(q_{x}+\frac{4\sqrt{3}\pi}{9a}))+4cos(\frac{\sqrt{3}}{2}a(q_{x}+\frac{4\sqrt{3}\pi}{9a}))cos(\frac{3}{2}q_{y}a)}

\quad \quad \approx \pm J \sqrt{2+3qa_{x}+\frac{3}{2}a^{2}q_{x}^{2}+(-2+\frac{3}{4}q_{x}^{2}a^{2}-3q_{x}a)\times(1-\frac{9}{8}q_{y}^{2}a^{2})}

\quad \quad=\pm \frac{3}{2}a\sqrt{q_{x}^{2}+q_{y}^{2}}+O(q^{2})=\pm \frac{3}{2}aq+O(q^{2})= \pm v_{F}\cdot q +O(q^{2})

where

v_{F}=\frac{3}{2}aq

is the Fermi Velocity。 If we omit the higher oder term， then we get

E(K+q)=\pm v_{F}\cdot q

緊束縛模型計算石墨烯能帶

Dirac Cone

As long as the time reversal symmetry and inversion symmetry is preserved， the Dirac cones always exist and are robust agains perturbation。

The

\hat{H}

can be written as ：

\hat{H}=v_{F}\sum_{k}(a_{k}^\dagger,b_{k}^\dagger)\begin{pmatrix} 0 \quad q_{x}-iq_{y} \\ q_{x}+iq_{y} \quad 0 \end{pmatrix}\begin{pmatrix}a_k\\b_k\end{pmatrix}

The kernel of the Hamiltonian is the matrix

H(K+q)=v_{F}\begin{bmatrix} 0 \quad q_{x}-iq_{y} \\ q_{x}+iq_{y} \quad 0 \end{bmatrix}

Using Pauli matrix

\sigma_{x}=\begin{bmatrix} 0\quad 1\\ 1\quad 0 \end{bmatrix}

\sigma_{y}=\begin{bmatrix} 0\quad -i\\ i \quad\quad 0 \end{bmatrix}

\sigma_{z}=\begin{bmatrix} 1\quad \quad 0\\ 0\quad -1 \end{bmatrix}

we get：

H(K+q)=v_{F}(q_{x}\sigma_{x}+q_{y}\sigma_{y})=v_{F}q\cdot \sigma

Notice that the electrons near Dirac Cone satisfy the massless Dirac equation：

i\hbar\frac{\partial}{\partial t}\psi=c\alpha\cdot (-i\hbar\triangledown+\beta mc^{2})\psi=c\alpha\cdot (-i\hbar\triangledown)\psi=H\psi

The

H=c\alpha\cdot (-i\hbar\triangledown)

and

H=v_{F}q\cdot\sigma

are of the same form。 The

\gamma_{i}\equiv-i\beta\alpha_{i}

and

\gamma_{i}

can be expressed using Pauli matrix。 For example， if we choose Pauli-Dirac basis， the 4x4 matrix

\gamma_{1}=\begin{bmatrix} 0 \quad-i\sigma_{1}\\ i\sigma_{i} \quad \quad0 \end{bmatrix}

。

For graphene， the lower band is filled and the upper band is empty， which is known as “half-filling”。 The “half” here means that the number of electrons

N_{e}

over the number of sites（ A and B）

N_{s}

is

\frac{N_{e}}{N_{s}}=\frac{1}{2}

To be clear， graphene is compound lattice： There are two sites in each unit cell。 So the number of electrons over the number of unit cells

\frac{N_{e}}{N_{u}}=1

。 However， one can tune the Fermi energy slightly away from the Dirac point by doping。

《固體理論》李正中（第二版）

《Advanced Quantum Mechanics》咯興林（第二版）

《Topological insulators Part III： tight-binding models》Sunkai

Phys620。nb（

2013）