- Email: [email protected]
Second-order envelope equation of graphene electrons
Solid State Communications 195 (2014) 95–101
Contents lists available at ScienceDirect
Solid State Communications journal homepage: www.elsevier.com/locate/ssc
Second-order envelope equation of graphene electrons Ji Luo Department of Physics and Institute for Functional Nanomaterials, University of Puerto Rico at Mayagüez, Mayagüez, PR 00681, USA
art ic l e i nf o
a b s t r a c t
Article history: Received 4 February 2014 Received in revised form 3 July 2014 Accepted 7 July 2014 by H. Akai Available online 12 July 2014
A treatment of graphene's electronic states based on the tight-binding method is presented. Like Dirac equation, this treatment uses envelope functions to eliminate crystal potential. Besides, a density-functionaltheory Kohn–Sham (KS) orbital of an isolated carbon atom is employed. By locally expanding envelope functions into second-order polynomials and by involving up to third-nearest atoms in calculating orbital integrals, the second-order envelope equation is obtained. This equation does not contain any experimental data except graphene's crystal structure, and its coefficients are determined through several kinds of integrals of the carbon KS orbital. As an improvement, it leads to more accurate energy dispersion than Dirac equation including the triangular warping effect and asymmetry for electrons and holes, and gives the Fermi velocity which is in good agreement with the experimental value. & 2014 Elsevier Ltd. All rights reserved.
Keywords: A. Graphene C. Electronic structure
1. Introduction One characteristic feature of graphene is that its electronic states are described by two envelope functions which vary gently at the unit-cell scale [1–3]. The envelope functions satisfy a Diraclike equation, which not only constitutes the basis of graphene's electronic properties, but also relates graphene physics to quantum electrodynamics [1,2]. Dirac equation was first obtained - -
according to a tight-binding (TB) based k U p method, where wave functions were supposed to have the specific form of Bloch functions of a Fermi wave vector multiplied by the factor -
-
-
- -
-
eið k k F Þ U r , with k , k F , and r respectively denoting a wave vector, a Fermi wave vector, and a position vector [4]. Envelope functions are simple because rapid oscillations within a unit cell have been separated out, and Dirac equation leads to a linear energy dispersion in the form of the Dirac cone [1,2,5]. Originally, TB calculation needs graphene's crystal potential for the hopping integrals of carbon 2pz orbitals, and if the eigenenergy of an isolated carbon atom is introduced, the atomic potential is also needed [6,7]. Since the crystal potential is difficult to obtain, it resorts to experimental or first-principles parameters to bypass the integrals [7,8]. The envelope-function method is rather different. Dirac equation contains only one parameter, the Fermi velocity vF . In definition vF is given by the momentum integral of carbon 2pz orbitals, and no potentials are involved [4,9]. This should have become another advantage of the envelopefunction method because empirical data could be avoided.
E-mail address: [email protected] http://dx.doi.org/10.1016/j.ssc.2014.07.007 0038-1098/& 2014 Elsevier Ltd. All rights reserved.
However, in practice vF is still obtained through experiments or first-principles calculations [7,8,10,11]. Simple as they are, Dirac equation and its results are valid only for low excited states whose wave vectors are near a Fermi one. For most excited states, even the TB method involving only nearest hopping integrals leads to an energy dispersion that deviates from the Dirac cone, and involving more integrals gives the dispersion that well accords with that obtained from first-principles calculations [7,12]. In addition to the energy dispersion, standard TB method also gives the electronic wave functions. However, such wave functions are expressed in terms of Bloch functions of a -
-
general k instead of k F . As a result, gentle envelope functions cannot be separated out to eliminate rapid oscillations within a unit cell. It is thus desirable that one can retain advantages of the envelope-function method and at the same time, improve its accuracy. In this work, we present a more general derivation of the equation satisfied by the envelope functions. Instead of presuming their simple form, we expand the envelope functions into series up to second-order terms for each unit cell. The resultant equation contains up to second-order derivatives of the envelope functions, and its coefficients are determined by several kinds of integrals of up to third-nearest carbon 2pz orbitals. These integrals depend only on the orbitals themselves and their distance, and no potential functions are involved. The carbon 2pz orbital function is obtained numerically as a Kohn–Sham (KS) orbital from densityfunctional theory (DFT) calculation [13,14]. As a result, the orbital integrals are calculated out and the second-order envelope equation is obtained without employing any experimental data except graphene's crystal structure. This equation leads to more accurate
96
J. Luo / Solid State Communications 195 (2014) 95–101
energy dispersion than Dirac equation. The deviation of the energy dispersion from the Dirac cone is now presented, such as the triangular warping effect and the asymmetry of the energy for electrons and holes. The Fermi velocity vF is also accurately obtained from the carbon KS orbital, and integrals of at least third-nearest orbitals are found necessary for this. In Section 2, we present the derivation of the second-order envelope equation. In Section 3, we present the numerical results of the equation. In Section 4, we discuss the results. Conclusions are presented in Section 5, and more details are presented in Appendix.
-
by x , y , and z . Two kinds of nonequivalent atoms are denoted by A and B. Two kinds of Dirac points K and K 0 are given by Fermi pffiffiffi pffiffiffi- wave vectors k F ¼ ð2π=3 3a0 Þð 3 x þ τ y Þ, with τ ¼ þ 1; 1 and a0 ¼ 0:142 nm being the carbon–carbon bond length. For an atom -
-
-
M at the position R M ¼ xM x þ yM y , we define -
-
-
¼ r RM; -
-
ϕM ð r Þ ¼ ϕð r M Þ;
-
ð1Þ
-
-
where r ¼ x x þ y y þ z z is the position vector and ϕ the 2pz orbital function of a carbon atom located at the origin. Originally graphene's electronic states are determined by the single-electron Hamiltonian: ℏ2 2 ∇ þ Vð r Þ; H^ ¼ 2m
ð2Þ
where m is the mass of an electron and V the periodical crystal potential including the average effect of electron–electron interaction. According to TB approximation, wave functions are expressed in terms of Bloch functions of the two-dimensional -
-
-
(2D) wave vector k ¼ kx x þky y , that is, - -
-
-
-
-
-
-
-
-
ψð r ; k Þ ¼ C 1 ð k Þ∑ ei k U R A ϕA ð r Þ þ C 2 ð k Þ∑ ei k U R B ϕB ð r Þ
ð3Þ
B
A
with C1, C2 being coefficients. Suppose -
-
-
ψ A ¼ eiπ=12 ∑ ei k F U R A ϕA ð r Þ;
-
-
-
ψ B ¼ e iπ=12 ∑ ei k F U R B ϕB ð r Þ
ð4Þ
B
A
-
are Bloch functions of a Fermi wave vector k F , where factors e 7 iπ =12 are introduced to simplify the results (see Appendix). Wave function (3) can always be expressed as - -
-
- -
- -
-
-
ψð r ; k Þ ¼ ψ 1 ð r ; k Þψ A ð r Þ þψ 2 ð r ; k Þψ B ð r Þ;
ψ1 ¼
ð5Þ
-
-
C 1 ð k Þ∑A ei k U R A ϕA ð r Þ -
ψ Að r Þ
;
ψ2 ¼
-
-
-
B atoms. Since ϕ is a localized function, in calculating the sums we include at most the third-nearest atoms of M and within this -
scope, we expand ψ j , j ¼ 1; 2 into power series at R M to secondorder terms, that is, -
-
-
-
-
-
-
-
-
: r M r M; - -
-
-
-
-
-
∇ψ j ð r ; k Þ ¼ ∇ψ j ð R M ; k Þ þ∇∇ψ j ð R M ; k Þ U r M ; - -
-
-
∇∇ψ j ð r ; k Þ ¼ ∇∇ψ j ð R M ; k Þ;
ð8Þ
- rMrM
where ∇∇ψ j and are dyads. As a result, a series of orbital integrals between two atoms M and N have to be calculated, which is R R R 3331 ϕM ϕN d r ; 1 ϕM ∇ϕN d r ; 1 r M ϕM ϕN d r ; ð9Þ R - R 331 r M ϕM ∇ϕN d r ; 1 r M r M ϕM ϕN d r : Of these integrals the first is a scalar one, the second and the third are vector ones, and the fourth and the fifth are dyad ones. For two nth-nearest atoms M and N with n ¼ 1; 2; 3, we define t MN as the unit vector along MN and n MN ¼ z t MN . For M ¼ N or n ¼ 0, we define t MN ¼ x and n MN ¼ y . According to the symmetry of ϕ (see Section 3), nonzero components of integrals (9) to be calculated are R 3sn ¼ 1 ϕM ϕN d r ; Z Z 33ð∇ϕN U t MN ÞϕM d r ; qn ¼ ð r M U t MN ÞϕM ϕN d r pn ¼ 1 1 Z 3un;1 ¼ ð r M U t MN Þð∇ϕN U t MN ÞϕM d r ; 1 Z 3un;2 ¼ ð r M U n MN Þð∇ϕN U n MN ÞϕM d r Z1 3un;3 ¼ ð r M U z Þð∇ϕN U z ÞϕM d r ; 1 Z 3vn;1 ¼ ð r M U t MN Þ2 ϕM ϕN d r ; 1 Z 3vn;2 ¼ ð r M U n MN Þ2 ϕM ϕN d r ; 1 Z 3vn;3 ¼ ð r M U z Þ2 ϕM ϕN d r : ð10Þ 1
C 2 ð k Þ∑B ei k U R B ϕB ð r Þ -
ψ Bð r Þ
-
-
:
ð6Þ
ℏ2 ℏ2 ðψ ∇2 ψ 1 þ ψ B ∇2 ψ 2 Þ ð∇ψ A U ∇ψ 1 þ ∇ψ B U ∇ψ 2 Þ ¼ εðψ A ψ 1 þ ψ B ψ 2 Þ: 2m A m
ð7Þ -
The crystal potential V disappears because eigen-energy of k F is set to be zero and corresponding Bloch functions approximately ^ ¼ 0 and H^ ψ ¼ 0. satisfy Hψ A
-
As a result, the triple integration of Eq. (7) multiplied by -
-
These two functions are envelope ones as long as they vary gently at the unit-cell scale. By substituting Eq. (5) into eigenequation H^ ψ ¼ εψ with ε being the eigen-energy, one obtains
-
eiπ=12 e i k F U R B ϕB ð r Þ, and conduct triple integration. According to Eq. (4), one has to calculate the following sums of integrals: R R 3∑N ei k F U ð R N R M Þ 1 ϕM ϕN ∇2 ψ j d r , ∑N ei k F U ð R N R M Þ 1 ϕM ∇ϕN U ∇ψ j d3 r , R 3and ∑N ei k F U ð R N R M Þ 1 ϕM ϕN ψ j d r , where N runs through all A or
as long as one defines -
-
ψ j ð r ; k Þ ¼ ψ j ð R M ; k Þ þ ∇ψ j ð R M ; k Þ U r M þ 12 ∇∇ψ j ð R M ; k Þ
The graphene is taken as xy-plane with x-axis parallel to one set of carbon–carbon bonds. Unit vectors of the axes are denoted
rM
-
M ¼ A or M ¼ B, multiply Eq. (7) by e iπ=12 e i k F U R A ϕA ð r Þ or
- -
2. Derivation
- -
To obtain equations of ψ 1 and ψ 2 , we choose a specific atom
B
e
-
iπ=12 i k F U R A
e
-
-
-
-
ϕA ð r Þ or eiπ=12 e i k F U R B ϕB ð r Þ leads to two equations -
-
of ψ j , j ¼ 1; 2 at R A and R B respectively [Eqs. (A9) and (A10) in Appendix]. Finally, for two atoms A and B within a unit cell, we -
-
-
-
-
-
express ψ j ð R M ; k Þ, ∇ψ j ð R M ; k Þ, ∇∇ψ j ð R M ; k Þ, j ¼ 1; 2, M ¼ A; B in -
terms of their values at the middle point R C by a series expansion similar to Eq. (8), that is, -
-
-
R C ¼ 12 ð R A þ R B Þ; -
-
-
r MC
-
-
¼ RM RC;
-
-
-
-
ψ j ð R M ; k Þ ¼ ψ j ð R C ; k Þ þ ∇ψ j ð R C ; k Þ U r MC -
-
-
-
þ 12 ∇∇ψ j ð R C ; k Þ : r MC r MC ;
J. Luo / Solid State Communications 195 (2014) 95–101 -
-
-
-
-
-
-
3. Numerical results
∇ψ j ð R M ; k Þ ¼ ∇ψ j ð R C ; k Þ þ ∇∇ψ j ð R C ; k Þ U r MC ; -
-
-
-
∇∇ψ j ð R M ; k Þ ¼ ∇∇ψ j ð R C ; k Þ:
ð11Þ
One finally obtains
εΨ þðεC x Dx Þ
∂Ψ ∂Ψ ∂2 Ψ þ ðεC y Dy Þ þ ðεC xx Dxx Þ 2 ∂x ∂y ∂x
þ ðεC yy Dyy Þ
∂2 Ψ ∂2 Ψ ∂2 Ψ ¼ 0; þ ðεC zz Dzz Þ 2 þðεC xy Dxy Þ ∂x ∂y ∂y2 ∂z
ð12Þ -
where Ψ ¼ ðψ 1 ψ 2 ÞT and all functions take their values at R C . Coefficients in eq. (12) are 2 2 matrixes τiq C y ¼ 32s σy;
a0 C x ¼ 3iq 2s σ x 2 σ z ;
C xx ¼ C yy ¼
C xy ¼
C zz ¼
v3 I2 ; 2s
Dzz ¼
q3 ¼ 1:798 10 11 m;
u0;1 ¼ 4:974 10 1 ;
u0;3 ¼ 4:924 10 1 ;
u1;1 ¼ 3:509 10
2
;
u1;2 ¼ 2:470 10 1 ;
u2;1 ¼ 1:400 10
1
;
u2;2 ¼ 8:794 10 2 ;
Dxy ¼
3τ ℏ ða0 p þ 2u2 Þσ x ; 4ms 2
ð13Þ
21
u3;2 ¼ 6:408 10 2 ; m2 ;
v0;3 ¼ 1:143 10 20 m2 ;
v1;1 ¼ 5:270 10 21 m2 ;
v1;2 ¼ 2:409 10 21 m2 ;
v0;1 ¼ 3:811 10
ℏ2 3ℏ2 ðs þ u1 ÞI 2 u2 σ y ; 2ms 4ms ℏ ðs þ 2u3 ÞI 2 ; 2ms
p3 ¼ 1:381 109 m 1 ;
m;
u3;1 ¼ 1:288 10 1 ;
3τiℏ2 p Dy ¼ σy 2ms
2
11
u2;3 ¼ 8:757 10 2 ;
ℏ2 3ℏ2 ðs þ u1 ÞI 2 ða0 p u2 Þσ y ; Dxx ¼ 2ms 4ms Dyy ¼
s2 ¼ 1:762 10 1 ; q1 ¼ 3:469 10
3τ ða0 q þ v2 Þσ x ; 4s
3iℏ2 p σx; Dx ¼ 2ms
We calculate the carbon 2pz orbital by using DFT program DMol3 employing the general gradient approximation [15,16]. The function ϕ is calculated out as the highest occupied molecular orbital of an isolated carbon atom. Numerical values of the function at lattice points with intervals 0:01 nm in a volume with size 2 2 2 nm3 are obtained. Contours of this orbital in the xz plane are presented in Fig. 1. The calculated ϕ is an even function of x or y, and an odd one of z. Because of this symmetry, some components of integrals (9) are zero. Numerical integration based on this orbital leads to the following results [in SI and see Eq. (10)]:
p1 ¼ 4:048 109 m 1 ;
a20 v1 3a0 q 3v2 þ I2 þ σy; 4s 8 4s 8s
v1 3v2 I2 þ σy; 4s 8s
97
v2;1 ¼ 4:330 10
21
v2;3 ¼ 3:259 10
21
v3;1 ¼ 3:805 10
21
2
m ;
v2;2 ¼ 1:086 10 21 m2 ;
2
m ; m2 ; v3;2 ¼ 8:401 10 22 m2 :
The corresponding coefficients in Eq. (12) are [in SI and see Eq. (13)] C x ¼ 5:319 10 11 iσ x 7:100 10 11 σ z
ðmÞ;
where I 2 is the 2 2 unit matrix, σ x , σ y , σ z are Pauli matrixes, and s, p, q, u1 , u2 , u3 , v1 , v2 , v3 are given by
C y ¼ 5:319 10
s ¼ 1 3s2 ;
q ¼ q1 q3 ;
C xx ¼ 2:056 10 21 I 2 8:598 10 22 σ y
ðm2 Þ;
u2 ¼ u1;1 u1;2 þ u3;1 u3;2 ;
C yy ¼ 4:576 10 21 I 2 þ 4:636 10 21 σ y
ðm2 Þ;
p ¼ p1 p3 ;
u1 ¼ 2u0;1 3u2;1 3u2;2 ; u3 ¼ u0;3 3u2;3 ;
v2 ¼ v1;1 v1;2 þ v3;1 v3;2 ;
v3 ¼ v0;3 3v2;3 :
11
C zz ¼ 1:758 10
v1 ¼ 2v0;1 3v2;1 3v2;2 ; ð14Þ
Eq. (12) is a set of equations of ψ 1 , ψ 2 at a set of discrete positions
-
R C . It can be regarded as the second-order differential equation satisfied by ψ 1 and ψ 2 , as long as they vary gently at the unit-cell scale. More details are presented in Appendix.
ð15Þ
τiσ y ðmÞ;
21
2
ðm Þ;
I2
C xy ¼ 1:305 10 20 τσ x Dx ¼ 1:036 10 28 iσ x Dy ¼ 1:036 10
28
ðm2 Þ; ðJmÞ;
τiσ y ðJmÞ;
Dxx ¼ 8:805 10 39 I 2 þ 1:872 10 39 σ y
ðJm2 Þ;
Dyy ¼ 8:805 10 39 I 2 9:229 10 39 σ y
ðJm2 Þ;
Dzz ¼ 1:546 10
40
2
I2
ðJm Þ;
Dxy ¼ 2:582 10 38 τσ x
ðJm2 Þ:
ð16Þ
Eq. (12) has solutions of the form ! c1 Ψ¼ eiðkx x þ ky y þ kz zÞ c2
ð17Þ
with kx , ky , kz being the wave-vector components. The eigenenergy ε is determined by the second-order secular equation
αε2 þ βε þ γ ¼ 0;
ð18Þ
where α, β, and γ are numerically given by (in SI)
α ¼ 1:000 þ6:323 10 21 ðk2x þ k2y Þ 3:516 10 21 k2z 2 3 þ 1:480 10 30 τkx ky 4:932 10 31 τky þ 3:487 10 42 kx 5:522 10 43 ky þ3:091 10 42 kz 4
4
4
1:435 10 40 kx ky 7:228 10 42 kx kz 1:609 10 41 ky kz 2 2
Fig. 1. Carbon 2pz orbital calculated from DFT program DMol3. The slice is in xz plane and has the size 0:6 0:6 nm2 , with the carbon nucleus located at the center. Values of contours are in atomic unit ð1 a:u: ¼ 2:598 1015 m 3=2 Þ.
2 2
β ¼ 6:586 10 39 ðk2x þ k2y Þ 3:093 10 40 k2z þ 5:828 10 48 τkx ky 1:943 10 48 τky 2
3
2 2
98
J. Luo / Solid State Communications 195 (2014) 95–101
Fig. 2. Energy wave-vector curves in sections perpendicular to kx ky plane, illustrating the deviation of the energy dispersion from the Dirac cone. The sections have the Dirac point K at the center and have different angles (marked in the figure) with respect to kx axis.
þ3:298 10 59 kx 4:995 10 60 ky þ 5:437 10 61 kz
4:768 10 76 kx ky 2:723 10 78 kx kz
5:237 10 58 kx ky 3:159 10 59 kx kz
2:723 10 78 ky kz
4
4
2 2
3:237 10 59 ky kz
2 2
4
2 2
2 2
2 2
ðJ2 Þ
ð19Þ
For kz ¼ 0, the energy dispersion is presented in Figs. 2 and 3.
ðJÞ
γ ¼ 1:074 10 56 ðk2x þ k2y Þ
4. Discussion
þ5:738 10 66 τkx ky 1:913 10 66 τky 2
3
þ7:401 10 77 kx 7:661 10 78 ky þ 2:391 10 80 kz 4
2 2
4
4
Second-order envelope Eq. (12) leads to more complicated energy dispersion and envelope functions than Dirac equation.
J. Luo / Solid State Communications 195 (2014) 95–101
99
Fig. 3. Wave-vector contours corresponding to different energy values. Dirac point K or K 0 is at the center of the contours. Triangular warping effect is illustrated. Contours for positive and negative energies demonstrate the asymmetry of the energy dispersion.
The energy dispersion now depends on wave-vector component kz as well, and the envelope functions on coordinate z as well. Components kx and ky are usually determined by the periodical boundary-conditions, yet a new quantum condition is needed to determine kz [17]. This may relate to some fundamental problems of quantum mechanics. The wave vector in Eq. (3) is a 2D one because graphene has 2D translational symmetry. The envelope functions in Dirac equation, however, are compulsively regarded as independent of z [4]. According to Eq. (6), in general envelope functions are functions of x, y, and z. For simplicity we consider states with kz ¼ 0, which can be compared with the results of Dirac equation. For kz ¼ 0 the energy dispersion is a curved surface in ε; kx ; ky space. Section curves in different vertical planes are presented in Fig. 2, demonstrating the deviation of the energy dispersion from the Dirac cone. A comparison with results in Ref. 12 indicates that the ε k curves accord with the first-principles calculations near the Dirac point. Coefficients α, β, γ in Eq. (18) are fourth-order polynomials of kx , ky , and kz . Graphene's reciprocal lattice has the constant b ¼ 4π =3a0 ¼ 2:95 1010 m 1 . We use polar coordinates kx ¼ k cos θ, ky ¼ k sin θ. According to Eq. (19), one can assume
α ¼ 1 unless k is comparable with b. If k o 0:01 b (k 108 m 1 or less), both the third-order and fourth-order terms in β and γ can be
neglected, and the energy dispersion is independent of θ. One has 3ℏ2 p ℏ2 u1 9pq 2 kþ 1þ þ 2 k : ε¼ 7 ð20Þ 2ms 2m s 2s The TB method involving nearest and second-nearest hopping R R 30 ^ ^ 0 31 ϕA HϕB d r , γ 0 ¼ 1 ϕA HϕA d r gives [6]
integrals γ 0 ¼
ε ¼ 7 32 γ 0 a0 k 94 γ 00 a20 k2 :
ð21Þ
We note that wave vector k in this work differs from that in Ref. [6] by the factor 2π . By comparing Eqs. (20) and (21), this work gives the nearest and second-nearest hopping integrals ℏ p γ 0 ¼ msa ¼ 3:036 eV; 0 2
γ 00 ¼
2ℏ2 u1 9pq 1 þ þ 2 ¼ 0:4530 eV: 2 s 2s 9ma0
ð22Þ
This provides a new way to determine TB parameters [2,8,11]. If k o 0:1 b (k 109 m 1 or less), the fourth-order terms in β and γ can be neglected and one has α ¼ 1:000, 3 2 2 pv2 9τ ℏ2 k þ 4ms β ¼ ℏm 1 þ us1 þ 9pq k sin 3θ; 2 qu2 þ 2 2s2
γ¼
9ℏ4 p2 2 9τℏ4 pu2 3 k þ k sin 3θ: 4m2 s2 4m2 s2
ð23Þ
100
J. Luo / Solid State Communications 195 (2014) 95–101
The factor sin 3θ in the k terms thus leads to the triangular warping of the energy dispersion. The expansion of the TB energy dispersion involving only nearest hopping integrals leads to the 2 factor sin 3θ in the k term [2]. This work indicates that more 3 accurately, this factor is in k term, which more accords with firstprinciples results [12]. Fig. 3 illustrates the triangular warping effect and the asymmetry of the positive and negative eigenenergies. 3
- -
It is interesting to compare this work with the usual k U p method [4]. If, instead of Eq. (8), one assumes - -
-
-
ψ j ð r ; k Þ ¼ ψ j ð R M ; k Þ;
- -
-
-
∇ψ j ð r ; k Þ ¼ ∇ψ j ð R M ; k Þ;
- -
∇∇ψ j ð r ; k Þ ¼ 0;
ð24Þ
then coefficients in Eq. (13) are zero except Dx and Dy. Eq. (12) reduces to Dirac equation ^ 3ℏp - σ U p Ψ ¼ εΨ ; 2ms
ð25Þ
^ where σ ¼ σ x x þ τσ y y and p ¼ ð iℏÞð x ∂=∂x þ y ∂=∂yÞ. Both Eqs. (20) and (25) lead to Fermi velocity
vF ¼
3ℏp ¼ 9:826 105 ms 1 ; 2ms
ð26Þ
the TB method to avoid empirical parameters but still retain its low computational cost.
5. Conclusions More accurate energy dispersion of graphene electrons can be obtained from a second-order envelope-function equation as an improvement of Dirac equation. DFT KS orbitals of isolated carbon atoms can play an important role in TB method by acting as Wannier functions, which endows KS orbitals with new significance. For graphene the improved TB method can totally avoid experimental data except the crystal structure. Further improvement is expected so that it can be applied to other systems.
Appendix: More details The symmetry of the orbital ϕ simplifies the calculation of integrals (9). For an atom M one has R R - 2 3R 2 331 ϕM d r ¼ 1; 1 ϕM ∇ϕM d r ¼ 0; 1 r M ϕM d r ¼ 0; Z ---3ϕM r M ∇ϕM d r ¼ u0;1 x x þu0;1 y y þ u0;3 z z ; 1 Z -- --3ϕ2M r M r M d r ¼ v0;1 x x þv0;1 y y þv0;3 z z ; ðA1Þ 1
which is in good agreement with the experimental value [10,11]. In Dirac equation, vF is defined according to the nearest momentum integral, that is, vF ¼ 3ℏp1 =2m. We note that p ¼ p1 p3 contains the third-nearest integral p3 , and s ¼ 1 3s2 the second-nearest one s2 . According to Eq. (15), p3 is comparable with p1 and 3s2 with s0 ¼ 1. If p3 and s2 are neglected, one has 3ℏp1 =2m ¼ 7:026 105 ms 1 . This explains the difference between the experimental value of vF and its definition according to the nearest momentum integral [4,6]. Even to determine vF from Dirac equation, at least the third-nearest integrals are necessary. More exact velocity of carriers can be calculated from -
the energy dispersion by v ¼ ð1=ℏÞ∇- ε, and the velocity difference of k
electrons and holes can be obtained from the asymmetry of the - -
energy dispersion. More accurate results in this work than the k U p method are achieved by expanding ψ 1 and ψ 2 to the second-order terms in Eq. (8) and by involving up to third-nearest orbital integrals. All numerical results in this work are based on the integrals of carbon 2pz orbitals, and no experimental data are involved except graphene's crystal structure. This is realized for two reasons. Firstly, the introduction of the envelope functions eliminates unknown crystal potential in Eq. (7). This is achieved due to graphene's unique crystal structure. Because of its π =3 rotational symmetry and the same carbon 2pz orbital for both A and B atoms, ^ by adjusting the zero point of the energy one has Hψ A ¼ 0 and ^ Hψ B ¼ 0 approximately and simultaneously. Secondly, the carbon 2pz orbital function is obtained from DFT calculation. There has been controversy about the significance of DFT KS orbitals. This work demonstrates that for carbon atoms in graphene, KS orbitals can be regarded as real single-electron wave functions in atomic systems and be used in TB wave functions. In the case of graphene electrons they constitute a good approximation of Wannier functions [18]. In this work, atomic orbitals are no longer symbols in TB method only for theoretical derivation. Instead they are used to numerically calculate the eigen-energy and envelope functions. Further improvement of this work may be the treatment of the crystal potential V. If, like the wave function, V can also be constructed in terms of atomic data calculated from DFT, this treatment may hopefully be applied to other systems such as semiconductors with direct band gaps [19–22]. Since atomic data can be calculated out as standard ones, this may ultimately enable
and for two nth-nearest atoms M and N with n ¼ 1; 2; 3, one has R 31 ϕM ϕN d r ¼ sn ; Z Z 33ϕM ∇ϕN d r ¼ pn t MN ; r M ϕM ϕN d r ¼ qn t MN ; 1 1 Z -3ϕM r M ∇ϕN d r ¼ un;1 t MN t MN þ un;2 n MN n MN þ un;3 z z ; 1 Z -- 3r M r M ϕM ϕN d r ¼ vn;1 t MN t MN þ vn;2 n MN n MN þ vn;3 z z : ðA2Þ 1
According to graphene's crystal structure, one can calculate out the following sums. For an atom A and its three nearest atoms B, one has - R 3∑ ei k F U ð R B R A Þ 1 ϕA ϕB d r ¼ 0; B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ∇ϕB d r ¼ eiπ=6 p1 ð x iτ y Þ; 2 1 B Z 3i 3∑ ei k F U ð R B R A Þ r A ϕA ϕB d r ¼ eiπ=6 q1 ð x iτ y Þ; 2 1 B Z 3∑ ei k F U ð R B R A Þ ϕA r A ∇ϕB d r
1
B
- 3i ¼ eiπ=6 ðu1;1 u1;2 Þð x þ iτ y Þð x þ iτ y Þ; 4 Z -
-
-
∑ ei k F U ð R B R A Þ B
¼
- -
1
3-
ϕA ϕB r A r A d r
- 3i iπ=6 e ðv1;1 v1;2 Þð x þ iτ y Þð x þiτ y Þ: 4
For an atom B and its three nearest atoms A, one has - R 3∑ ei k F U ð R A R B Þ 1 ϕB ϕA d r ¼ 0; A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ∇ϕA d r ¼ e iπ=6 p1 ð x þ iτ y Þ; 2 1 A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ϕA r B d r ¼ e iπ=6 q1 ð x þ iτ y Þ; 2 1 A Z 3i k F U ð R A R BÞ ϕB r B ∇ϕA d r ∑e
A
1
- 3i ¼ e iπ=6 ðu1;1 u1;2 Þð x iτ y Þð x iτ y Þ; 4
ðA3Þ
J. Luo / Solid State Communications 195 (2014) 95–101 -
-
-
∑ ei k F U ð R A R B Þ
Z
- -
1
A
-
3-
For an atom A and its six second-nearest atoms A , one has -
-
∑e A
0
-
-
-
∑ei k F U ð R A0 R A Þ A0
-
-
-
-
-
-
Z
0
i k F U ð R A0 R A Þ
∑e A
1 ϕA ϕA
0
0
1
ϕA ∇ϕA0 d r ¼ 0;
Eqs. (A3)–(A8) indicate that Pauli matrixes enter the envelope
-
3-
--
--
-
-
B
-
-
-
∑ei k F U ð R B0 R B Þ B0
-
-
-
∑ei k F U ð R B0 R B Þ B0
-
-
-
∑ei k F U ð R B0 R B Þ B0
R
1 ϕB ϕB
0
ðA5Þ
1
1
Z
1
at
3-
-
3-
1 --
2
--
--
--
¼ 32 ½ðv2;1 þ v2;2 Þð x x þ y y Þ þ 2v2;3 z z :
ðA6Þ
-
B
- -
1
3-
ϕA ϕB r A r A d r
- 3i ¼ eiπ=6 ðv3;1 v3;2 Þð x þ iτ y Þð x þ iτ y Þ: 4
ðA7Þ
For an atom B and its three third-nearest atoms A, one has - R 3∑ ei k F U ð R A R B Þ 1 ϕB ϕA d r ¼ 0; A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ∇ϕA d r ¼ e iπ=6 p3 ð x þ iτ y Þ; 2ℏ 1 A Z 3i 3∑ ei k F U ð R A R B Þ ϕA ϕB r B d r ¼ e iπ=6 q3 ð x þ iτ y Þ; 2 1 A Z 3∑ ei k F U ð R A R B Þ ϕB r B ∇ϕA d r
A
-
RB
according
to
1
- 3i ¼ e iπ=6 ðu3;1 u3;2 Þð x iτ y Þð x iτ y Þ; 4
2
2
2
- 3iℏ2 p 3iℏ2 u2 ∇ψ 1 U ð x þτi y Þ þ ∇∇ψ 1 : ð x τi y Þð x iτ y Þ 2m 4m v1 v3 ∂2 ψ ¼ sεψ 2 þ ε∇22 ψ 2 þ ε 22 4 2 ∂z - 3iq 3iv2 ε∇ψ 1 U ð x þ iτ y Þ ε∇∇ψ 1 : ð x iτ y Þð x iτ y Þ; þ 2 8
ðA10Þ
References
1
-
at
where ψ 1 , ψ 2 and their derivatives take values at R B .
- 3i ¼ eiπ=6 ðu3;1 u3;2 Þð x þ iτ y Þð x þ iτ y Þ; 4 Z -
series
-
- R 3∑ ei k F U ð R B R A Þ 1 ϕA ϕB d r ¼ 0; B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ∇ϕB d r ¼ eiπ=6 p3 ð x iτ y Þ; 2ℏ 1 B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ϕB r A d r ¼ eiπ=6 q3 ð x iτ y Þ; 2 1 B Z 3∑ ei k F U ð R B R A Þ ϕA r A ∇ϕB d r
∑ ei k F U ð R B R A Þ
power
For an atom A and its three third-nearest atoms B, one has
B
into
u1 2 ∇2 ψ 2 ℏ mu3 ∂∂zψ22 ℏ2ms∇2 ψ 2 ℏ2m
¼ 32 ½ðu2;1 þ u2;2 Þð x x þ y y Þ þ 2u2;3 z z ; Z - 3∑ei k F U ð R B0 R B Þ ϕB ϕB0 r B r B d r B0
-
ϕB ð r Þ and conducts the triple integration. The expan-
sion of ψ 1 , ψ 2 Eq. (8) leads to
ϕB r B ∇ϕB0 d r --
-
e
e
3-
--
R A . For a chosen atom B, one multiplies Eq. (7) by
iπ=12 i k F U R B
ϕB ϕB0 r B d r ¼ 0; -
-
-
ϕB ∇ϕB0 d r ¼ 0;
Z
ðA9Þ
where ∇22 ¼ ∂2 =∂x2 þ ∂2 =∂y2 , ψ 1 , ψ 2 and their derivatives take values
3-
d r ¼ 3s2 ;
Z
2
- 3iℏ2 p 3iℏ2 u2 ∇ψ 2 U ð x τi y Þ ∇∇ψ 2 : ð x þ τi y Þð x þ τi y Þ 2m 4m v1 v3 ∂2 ψ ¼ sεψ 1 þ ε∇22 ψ 1 þ ε 21 4 2 ∂z - 3iq 3iv2 ε∇ψ 2 U ð x τi y Þ þ ε∇∇ψ 2 : ð x þτi y Þð x þ τi y Þ; þ 2 8
For an atom B and its six second-nearest atoms B , one has -
2
0
0
2
2
¼ 32 ½ðv2;1 þ v2;2 Þð x x þ y y Þ þ 2v2;3 z z :
∑ei k F U ð R B0 R B Þ
-
u1 2 ∇2 ψ 1 ℏ mu3 ∂∂zψ21 ℏ2ms∇2 ψ 1 ℏ2m
--
--
1
--
-
-
¼ 32 ½ðu2;1 þ u2;2 Þð x x þ y y Þ þ 2u2;3 z z ; Z - 3∑ei k F U ð R A0 R A Þ ϕA ϕA0 r A r A d r A0
-
expansion of ψ 1 , ψ 2 into power series at R A according to Eq. (8) leads to
ϕA r A ∇ϕA0 d r --
ðA8Þ
e iπ=12 e i k F U R A ϕA ð r Þ and conducts the triple integration. The
3ϕA ϕA0 r A d r ¼ 0; -
1
1
3-
-
d r ¼ 3s2 ; 3-
1
- -
ϕA ϕB r B r B d r
equation due to the vectors x 7 iτ y , and their factors e 7 iπ =6 must be eliminated by e 7 iπ =12 in Eq. (4). For a chosen atom A, one multiplies Eq. (7) by
3-
Z Z
∑ei k F U ð R A0 R A Þ A
R
Z
- 3i ¼ e iπ=6 ðv3;1 v3;2 Þð x iτ y Þð x iτ y Þ: 4
ðA4Þ 0
-
-
A
- 3i ¼ e iπ=6 ðv1;1 v1;2 Þð x iτ y Þð x iτ y Þ: 4
i k F U ð R A0 R A Þ
-
∑ ei k F U ð R A R B Þ
ϕB ϕA r B r B d r
101
[1] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [2] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109. [3] L. Brey, H.A. Fertig, Phys. Rev. B 73 (2006) 235411. [4] D.P. DiVincenzo, E.J. Mele, Phys. Rev. B 29 (1984) 1685. [5] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197. [6] P.R. Wallace, Phys. Rev. 71 (1947) 622. [7] R. Saito, G. Dresselhaus, M.S. Dresselhaus, Physical Properties of Carbon Nanotubes, Imperial, London, 1998. [8] R. Saito, M. Fujita, G. Dresselhaus, M.S. Dresselhaus, Appl. Phys. Lett. 60 (1992) 2204. [9] J.W. McClure, Phys. Rev. 104 (1956) 666. [10] M.S. Dresselhaus, G. Dresselhaus, Adv. Phys. 51 (2002) 1. [11] R.S. Deacon, K.C. Chuang, R.J. Nicholas, K.S. Novoselov, A.K. Geim, Phys. Rev. B 76 (2007) 081406(R). [12] S. Reich, J. Maultzsch, C. Thomsen, P. Ordejon, Phys. Rev. B 66 (2002) 035412. [13] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [14] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [15] B. Delley, J. Chem. Phys. 92 (1990) 508; B. Delley, J. Chem. Phys. 113 (2000) 7756. [16] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [17] C. Kittel, Introduction to Solid State Physics, Eighth edition, John Wiley & Sons Inc., New York, 2005. [18] G. Wannier, Phys. Rev. 52 (1937) 191. [19] J.P. Perdew, M. Levy, Phys. Rev. Lett. 51 (1983) 1884. [20] L.J. Sham, M. Schluter, Phys. Rev. Lett. 51 (1983) 1888. [21] A.J. Cohen, P. Mori-Sanchez, W. Yang, Phys. Rev. B 77 (2008) 115123. [22] M.K.Y. Chan, G. Ceder, Phys. Rev. Lett. 105 (2010) 196403.
Contents lists available at ScienceDirect
Solid State Communications journal homepage: www.elsevier.com/locate/ssc
Second-order envelope equation of graphene electrons Ji Luo Department of Physics and Institute for Functional Nanomaterials, University of Puerto Rico at Mayagüez, Mayagüez, PR 00681, USA
art ic l e i nf o
a b s t r a c t
Article history: Received 4 February 2014 Received in revised form 3 July 2014 Accepted 7 July 2014 by H. Akai Available online 12 July 2014
A treatment of graphene's electronic states based on the tight-binding method is presented. Like Dirac equation, this treatment uses envelope functions to eliminate crystal potential. Besides, a density-functionaltheory Kohn–Sham (KS) orbital of an isolated carbon atom is employed. By locally expanding envelope functions into second-order polynomials and by involving up to third-nearest atoms in calculating orbital integrals, the second-order envelope equation is obtained. This equation does not contain any experimental data except graphene's crystal structure, and its coefficients are determined through several kinds of integrals of the carbon KS orbital. As an improvement, it leads to more accurate energy dispersion than Dirac equation including the triangular warping effect and asymmetry for electrons and holes, and gives the Fermi velocity which is in good agreement with the experimental value. & 2014 Elsevier Ltd. All rights reserved.
Keywords: A. Graphene C. Electronic structure
1. Introduction One characteristic feature of graphene is that its electronic states are described by two envelope functions which vary gently at the unit-cell scale [1–3]. The envelope functions satisfy a Diraclike equation, which not only constitutes the basis of graphene's electronic properties, but also relates graphene physics to quantum electrodynamics [1,2]. Dirac equation was first obtained - -
according to a tight-binding (TB) based k U p method, where wave functions were supposed to have the specific form of Bloch functions of a Fermi wave vector multiplied by the factor -
-
-
- -
-
eið k k F Þ U r , with k , k F , and r respectively denoting a wave vector, a Fermi wave vector, and a position vector [4]. Envelope functions are simple because rapid oscillations within a unit cell have been separated out, and Dirac equation leads to a linear energy dispersion in the form of the Dirac cone [1,2,5]. Originally, TB calculation needs graphene's crystal potential for the hopping integrals of carbon 2pz orbitals, and if the eigenenergy of an isolated carbon atom is introduced, the atomic potential is also needed [6,7]. Since the crystal potential is difficult to obtain, it resorts to experimental or first-principles parameters to bypass the integrals [7,8]. The envelope-function method is rather different. Dirac equation contains only one parameter, the Fermi velocity vF . In definition vF is given by the momentum integral of carbon 2pz orbitals, and no potentials are involved [4,9]. This should have become another advantage of the envelopefunction method because empirical data could be avoided.
E-mail address: [email protected] http://dx.doi.org/10.1016/j.ssc.2014.07.007 0038-1098/& 2014 Elsevier Ltd. All rights reserved.
However, in practice vF is still obtained through experiments or first-principles calculations [7,8,10,11]. Simple as they are, Dirac equation and its results are valid only for low excited states whose wave vectors are near a Fermi one. For most excited states, even the TB method involving only nearest hopping integrals leads to an energy dispersion that deviates from the Dirac cone, and involving more integrals gives the dispersion that well accords with that obtained from first-principles calculations [7,12]. In addition to the energy dispersion, standard TB method also gives the electronic wave functions. However, such wave functions are expressed in terms of Bloch functions of a -
-
general k instead of k F . As a result, gentle envelope functions cannot be separated out to eliminate rapid oscillations within a unit cell. It is thus desirable that one can retain advantages of the envelope-function method and at the same time, improve its accuracy. In this work, we present a more general derivation of the equation satisfied by the envelope functions. Instead of presuming their simple form, we expand the envelope functions into series up to second-order terms for each unit cell. The resultant equation contains up to second-order derivatives of the envelope functions, and its coefficients are determined by several kinds of integrals of up to third-nearest carbon 2pz orbitals. These integrals depend only on the orbitals themselves and their distance, and no potential functions are involved. The carbon 2pz orbital function is obtained numerically as a Kohn–Sham (KS) orbital from densityfunctional theory (DFT) calculation [13,14]. As a result, the orbital integrals are calculated out and the second-order envelope equation is obtained without employing any experimental data except graphene's crystal structure. This equation leads to more accurate
96
J. Luo / Solid State Communications 195 (2014) 95–101
energy dispersion than Dirac equation. The deviation of the energy dispersion from the Dirac cone is now presented, such as the triangular warping effect and the asymmetry of the energy for electrons and holes. The Fermi velocity vF is also accurately obtained from the carbon KS orbital, and integrals of at least third-nearest orbitals are found necessary for this. In Section 2, we present the derivation of the second-order envelope equation. In Section 3, we present the numerical results of the equation. In Section 4, we discuss the results. Conclusions are presented in Section 5, and more details are presented in Appendix.
-
by x , y , and z . Two kinds of nonequivalent atoms are denoted by A and B. Two kinds of Dirac points K and K 0 are given by Fermi pffiffiffi pffiffiffi- wave vectors k F ¼ ð2π=3 3a0 Þð 3 x þ τ y Þ, with τ ¼ þ 1; 1 and a0 ¼ 0:142 nm being the carbon–carbon bond length. For an atom -
-
-
M at the position R M ¼ xM x þ yM y , we define -
-
-
¼ r RM; -
-
ϕM ð r Þ ¼ ϕð r M Þ;
-
ð1Þ
-
-
where r ¼ x x þ y y þ z z is the position vector and ϕ the 2pz orbital function of a carbon atom located at the origin. Originally graphene's electronic states are determined by the single-electron Hamiltonian: ℏ2 2 ∇ þ Vð r Þ; H^ ¼ 2m
ð2Þ
where m is the mass of an electron and V the periodical crystal potential including the average effect of electron–electron interaction. According to TB approximation, wave functions are expressed in terms of Bloch functions of the two-dimensional -
-
-
(2D) wave vector k ¼ kx x þky y , that is, - -
-
-
-
-
-
-
-
-
ψð r ; k Þ ¼ C 1 ð k Þ∑ ei k U R A ϕA ð r Þ þ C 2 ð k Þ∑ ei k U R B ϕB ð r Þ
ð3Þ
B
A
with C1, C2 being coefficients. Suppose -
-
-
ψ A ¼ eiπ=12 ∑ ei k F U R A ϕA ð r Þ;
-
-
-
ψ B ¼ e iπ=12 ∑ ei k F U R B ϕB ð r Þ
ð4Þ
B
A
-
are Bloch functions of a Fermi wave vector k F , where factors e 7 iπ =12 are introduced to simplify the results (see Appendix). Wave function (3) can always be expressed as - -
-
- -
- -
-
-
ψð r ; k Þ ¼ ψ 1 ð r ; k Þψ A ð r Þ þψ 2 ð r ; k Þψ B ð r Þ;
ψ1 ¼
ð5Þ
-
-
C 1 ð k Þ∑A ei k U R A ϕA ð r Þ -
ψ Að r Þ
;
ψ2 ¼
-
-
-
B atoms. Since ϕ is a localized function, in calculating the sums we include at most the third-nearest atoms of M and within this -
scope, we expand ψ j , j ¼ 1; 2 into power series at R M to secondorder terms, that is, -
-
-
-
-
-
-
-
-
: r M r M; - -
-
-
-
-
-
∇ψ j ð r ; k Þ ¼ ∇ψ j ð R M ; k Þ þ∇∇ψ j ð R M ; k Þ U r M ; - -
-
-
∇∇ψ j ð r ; k Þ ¼ ∇∇ψ j ð R M ; k Þ;
ð8Þ
- rMrM
where ∇∇ψ j and are dyads. As a result, a series of orbital integrals between two atoms M and N have to be calculated, which is R R R 3331 ϕM ϕN d r ; 1 ϕM ∇ϕN d r ; 1 r M ϕM ϕN d r ; ð9Þ R - R 331 r M ϕM ∇ϕN d r ; 1 r M r M ϕM ϕN d r : Of these integrals the first is a scalar one, the second and the third are vector ones, and the fourth and the fifth are dyad ones. For two nth-nearest atoms M and N with n ¼ 1; 2; 3, we define t MN as the unit vector along MN and n MN ¼ z t MN . For M ¼ N or n ¼ 0, we define t MN ¼ x and n MN ¼ y . According to the symmetry of ϕ (see Section 3), nonzero components of integrals (9) to be calculated are R 3sn ¼ 1 ϕM ϕN d r ; Z Z 33ð∇ϕN U t MN ÞϕM d r ; qn ¼ ð r M U t MN ÞϕM ϕN d r pn ¼ 1 1 Z 3un;1 ¼ ð r M U t MN Þð∇ϕN U t MN ÞϕM d r ; 1 Z 3un;2 ¼ ð r M U n MN Þð∇ϕN U n MN ÞϕM d r Z1 3un;3 ¼ ð r M U z Þð∇ϕN U z ÞϕM d r ; 1 Z 3vn;1 ¼ ð r M U t MN Þ2 ϕM ϕN d r ; 1 Z 3vn;2 ¼ ð r M U n MN Þ2 ϕM ϕN d r ; 1 Z 3vn;3 ¼ ð r M U z Þ2 ϕM ϕN d r : ð10Þ 1
C 2 ð k Þ∑B ei k U R B ϕB ð r Þ -
ψ Bð r Þ
-
-
:
ð6Þ
ℏ2 ℏ2 ðψ ∇2 ψ 1 þ ψ B ∇2 ψ 2 Þ ð∇ψ A U ∇ψ 1 þ ∇ψ B U ∇ψ 2 Þ ¼ εðψ A ψ 1 þ ψ B ψ 2 Þ: 2m A m
ð7Þ -
The crystal potential V disappears because eigen-energy of k F is set to be zero and corresponding Bloch functions approximately ^ ¼ 0 and H^ ψ ¼ 0. satisfy Hψ A
-
As a result, the triple integration of Eq. (7) multiplied by -
-
These two functions are envelope ones as long as they vary gently at the unit-cell scale. By substituting Eq. (5) into eigenequation H^ ψ ¼ εψ with ε being the eigen-energy, one obtains
-
eiπ=12 e i k F U R B ϕB ð r Þ, and conduct triple integration. According to Eq. (4), one has to calculate the following sums of integrals: R R 3∑N ei k F U ð R N R M Þ 1 ϕM ϕN ∇2 ψ j d r , ∑N ei k F U ð R N R M Þ 1 ϕM ∇ϕN U ∇ψ j d3 r , R 3and ∑N ei k F U ð R N R M Þ 1 ϕM ϕN ψ j d r , where N runs through all A or
as long as one defines -
-
ψ j ð r ; k Þ ¼ ψ j ð R M ; k Þ þ ∇ψ j ð R M ; k Þ U r M þ 12 ∇∇ψ j ð R M ; k Þ
The graphene is taken as xy-plane with x-axis parallel to one set of carbon–carbon bonds. Unit vectors of the axes are denoted
rM
-
M ¼ A or M ¼ B, multiply Eq. (7) by e iπ=12 e i k F U R A ϕA ð r Þ or
- -
2. Derivation
- -
To obtain equations of ψ 1 and ψ 2 , we choose a specific atom
B
e
-
iπ=12 i k F U R A
e
-
-
-
-
ϕA ð r Þ or eiπ=12 e i k F U R B ϕB ð r Þ leads to two equations -
-
of ψ j , j ¼ 1; 2 at R A and R B respectively [Eqs. (A9) and (A10) in Appendix]. Finally, for two atoms A and B within a unit cell, we -
-
-
-
-
-
express ψ j ð R M ; k Þ, ∇ψ j ð R M ; k Þ, ∇∇ψ j ð R M ; k Þ, j ¼ 1; 2, M ¼ A; B in -
terms of their values at the middle point R C by a series expansion similar to Eq. (8), that is, -
-
-
R C ¼ 12 ð R A þ R B Þ; -
-
-
r MC
-
-
¼ RM RC;
-
-
-
-
ψ j ð R M ; k Þ ¼ ψ j ð R C ; k Þ þ ∇ψ j ð R C ; k Þ U r MC -
-
-
-
þ 12 ∇∇ψ j ð R C ; k Þ : r MC r MC ;
J. Luo / Solid State Communications 195 (2014) 95–101 -
-
-
-
-
-
-
3. Numerical results
∇ψ j ð R M ; k Þ ¼ ∇ψ j ð R C ; k Þ þ ∇∇ψ j ð R C ; k Þ U r MC ; -
-
-
-
∇∇ψ j ð R M ; k Þ ¼ ∇∇ψ j ð R C ; k Þ:
ð11Þ
One finally obtains
εΨ þðεC x Dx Þ
∂Ψ ∂Ψ ∂2 Ψ þ ðεC y Dy Þ þ ðεC xx Dxx Þ 2 ∂x ∂y ∂x
þ ðεC yy Dyy Þ
∂2 Ψ ∂2 Ψ ∂2 Ψ ¼ 0; þ ðεC zz Dzz Þ 2 þðεC xy Dxy Þ ∂x ∂y ∂y2 ∂z
ð12Þ -
where Ψ ¼ ðψ 1 ψ 2 ÞT and all functions take their values at R C . Coefficients in eq. (12) are 2 2 matrixes τiq C y ¼ 32s σy;
a0 C x ¼ 3iq 2s σ x 2 σ z ;
C xx ¼ C yy ¼
C xy ¼
C zz ¼
v3 I2 ; 2s
Dzz ¼
q3 ¼ 1:798 10 11 m;
u0;1 ¼ 4:974 10 1 ;
u0;3 ¼ 4:924 10 1 ;
u1;1 ¼ 3:509 10
2
;
u1;2 ¼ 2:470 10 1 ;
u2;1 ¼ 1:400 10
1
;
u2;2 ¼ 8:794 10 2 ;
Dxy ¼
3τ ℏ ða0 p þ 2u2 Þσ x ; 4ms 2
ð13Þ
21
u3;2 ¼ 6:408 10 2 ; m2 ;
v0;3 ¼ 1:143 10 20 m2 ;
v1;1 ¼ 5:270 10 21 m2 ;
v1;2 ¼ 2:409 10 21 m2 ;
v0;1 ¼ 3:811 10
ℏ2 3ℏ2 ðs þ u1 ÞI 2 u2 σ y ; 2ms 4ms ℏ ðs þ 2u3 ÞI 2 ; 2ms
p3 ¼ 1:381 109 m 1 ;
m;
u3;1 ¼ 1:288 10 1 ;
3τiℏ2 p Dy ¼ σy 2ms
2
11
u2;3 ¼ 8:757 10 2 ;
ℏ2 3ℏ2 ðs þ u1 ÞI 2 ða0 p u2 Þσ y ; Dxx ¼ 2ms 4ms Dyy ¼
s2 ¼ 1:762 10 1 ; q1 ¼ 3:469 10
3τ ða0 q þ v2 Þσ x ; 4s
3iℏ2 p σx; Dx ¼ 2ms
We calculate the carbon 2pz orbital by using DFT program DMol3 employing the general gradient approximation [15,16]. The function ϕ is calculated out as the highest occupied molecular orbital of an isolated carbon atom. Numerical values of the function at lattice points with intervals 0:01 nm in a volume with size 2 2 2 nm3 are obtained. Contours of this orbital in the xz plane are presented in Fig. 1. The calculated ϕ is an even function of x or y, and an odd one of z. Because of this symmetry, some components of integrals (9) are zero. Numerical integration based on this orbital leads to the following results [in SI and see Eq. (10)]:
p1 ¼ 4:048 109 m 1 ;
a20 v1 3a0 q 3v2 þ I2 þ σy; 4s 8 4s 8s
v1 3v2 I2 þ σy; 4s 8s
97
v2;1 ¼ 4:330 10
21
v2;3 ¼ 3:259 10
21
v3;1 ¼ 3:805 10
21
2
m ;
v2;2 ¼ 1:086 10 21 m2 ;
2
m ; m2 ; v3;2 ¼ 8:401 10 22 m2 :
The corresponding coefficients in Eq. (12) are [in SI and see Eq. (13)] C x ¼ 5:319 10 11 iσ x 7:100 10 11 σ z
ðmÞ;
where I 2 is the 2 2 unit matrix, σ x , σ y , σ z are Pauli matrixes, and s, p, q, u1 , u2 , u3 , v1 , v2 , v3 are given by
C y ¼ 5:319 10
s ¼ 1 3s2 ;
q ¼ q1 q3 ;
C xx ¼ 2:056 10 21 I 2 8:598 10 22 σ y
ðm2 Þ;
u2 ¼ u1;1 u1;2 þ u3;1 u3;2 ;
C yy ¼ 4:576 10 21 I 2 þ 4:636 10 21 σ y
ðm2 Þ;
p ¼ p1 p3 ;
u1 ¼ 2u0;1 3u2;1 3u2;2 ; u3 ¼ u0;3 3u2;3 ;
v2 ¼ v1;1 v1;2 þ v3;1 v3;2 ;
v3 ¼ v0;3 3v2;3 :
11
C zz ¼ 1:758 10
v1 ¼ 2v0;1 3v2;1 3v2;2 ; ð14Þ
Eq. (12) is a set of equations of ψ 1 , ψ 2 at a set of discrete positions
-
R C . It can be regarded as the second-order differential equation satisfied by ψ 1 and ψ 2 , as long as they vary gently at the unit-cell scale. More details are presented in Appendix.
ð15Þ
τiσ y ðmÞ;
21
2
ðm Þ;
I2
C xy ¼ 1:305 10 20 τσ x Dx ¼ 1:036 10 28 iσ x Dy ¼ 1:036 10
28
ðm2 Þ; ðJmÞ;
τiσ y ðJmÞ;
Dxx ¼ 8:805 10 39 I 2 þ 1:872 10 39 σ y
ðJm2 Þ;
Dyy ¼ 8:805 10 39 I 2 9:229 10 39 σ y
ðJm2 Þ;
Dzz ¼ 1:546 10
40
2
I2
ðJm Þ;
Dxy ¼ 2:582 10 38 τσ x
ðJm2 Þ:
ð16Þ
Eq. (12) has solutions of the form ! c1 Ψ¼ eiðkx x þ ky y þ kz zÞ c2
ð17Þ
with kx , ky , kz being the wave-vector components. The eigenenergy ε is determined by the second-order secular equation
αε2 þ βε þ γ ¼ 0;
ð18Þ
where α, β, and γ are numerically given by (in SI)
α ¼ 1:000 þ6:323 10 21 ðk2x þ k2y Þ 3:516 10 21 k2z 2 3 þ 1:480 10 30 τkx ky 4:932 10 31 τky þ 3:487 10 42 kx 5:522 10 43 ky þ3:091 10 42 kz 4
4
4
1:435 10 40 kx ky 7:228 10 42 kx kz 1:609 10 41 ky kz 2 2
Fig. 1. Carbon 2pz orbital calculated from DFT program DMol3. The slice is in xz plane and has the size 0:6 0:6 nm2 , with the carbon nucleus located at the center. Values of contours are in atomic unit ð1 a:u: ¼ 2:598 1015 m 3=2 Þ.
2 2
β ¼ 6:586 10 39 ðk2x þ k2y Þ 3:093 10 40 k2z þ 5:828 10 48 τkx ky 1:943 10 48 τky 2
3
2 2
98
J. Luo / Solid State Communications 195 (2014) 95–101
Fig. 2. Energy wave-vector curves in sections perpendicular to kx ky plane, illustrating the deviation of the energy dispersion from the Dirac cone. The sections have the Dirac point K at the center and have different angles (marked in the figure) with respect to kx axis.
þ3:298 10 59 kx 4:995 10 60 ky þ 5:437 10 61 kz
4:768 10 76 kx ky 2:723 10 78 kx kz
5:237 10 58 kx ky 3:159 10 59 kx kz
2:723 10 78 ky kz
4
4
2 2
3:237 10 59 ky kz
2 2
4
2 2
2 2
2 2
ðJ2 Þ
ð19Þ
For kz ¼ 0, the energy dispersion is presented in Figs. 2 and 3.
ðJÞ
γ ¼ 1:074 10 56 ðk2x þ k2y Þ
4. Discussion
þ5:738 10 66 τkx ky 1:913 10 66 τky 2
3
þ7:401 10 77 kx 7:661 10 78 ky þ 2:391 10 80 kz 4
2 2
4
4
Second-order envelope Eq. (12) leads to more complicated energy dispersion and envelope functions than Dirac equation.
J. Luo / Solid State Communications 195 (2014) 95–101
99
Fig. 3. Wave-vector contours corresponding to different energy values. Dirac point K or K 0 is at the center of the contours. Triangular warping effect is illustrated. Contours for positive and negative energies demonstrate the asymmetry of the energy dispersion.
The energy dispersion now depends on wave-vector component kz as well, and the envelope functions on coordinate z as well. Components kx and ky are usually determined by the periodical boundary-conditions, yet a new quantum condition is needed to determine kz [17]. This may relate to some fundamental problems of quantum mechanics. The wave vector in Eq. (3) is a 2D one because graphene has 2D translational symmetry. The envelope functions in Dirac equation, however, are compulsively regarded as independent of z [4]. According to Eq. (6), in general envelope functions are functions of x, y, and z. For simplicity we consider states with kz ¼ 0, which can be compared with the results of Dirac equation. For kz ¼ 0 the energy dispersion is a curved surface in ε; kx ; ky space. Section curves in different vertical planes are presented in Fig. 2, demonstrating the deviation of the energy dispersion from the Dirac cone. A comparison with results in Ref. 12 indicates that the ε k curves accord with the first-principles calculations near the Dirac point. Coefficients α, β, γ in Eq. (18) are fourth-order polynomials of kx , ky , and kz . Graphene's reciprocal lattice has the constant b ¼ 4π =3a0 ¼ 2:95 1010 m 1 . We use polar coordinates kx ¼ k cos θ, ky ¼ k sin θ. According to Eq. (19), one can assume
α ¼ 1 unless k is comparable with b. If k o 0:01 b (k 108 m 1 or less), both the third-order and fourth-order terms in β and γ can be
neglected, and the energy dispersion is independent of θ. One has 3ℏ2 p ℏ2 u1 9pq 2 kþ 1þ þ 2 k : ε¼ 7 ð20Þ 2ms 2m s 2s The TB method involving nearest and second-nearest hopping R R 30 ^ ^ 0 31 ϕA HϕB d r , γ 0 ¼ 1 ϕA HϕA d r gives [6]
integrals γ 0 ¼
ε ¼ 7 32 γ 0 a0 k 94 γ 00 a20 k2 :
ð21Þ
We note that wave vector k in this work differs from that in Ref. [6] by the factor 2π . By comparing Eqs. (20) and (21), this work gives the nearest and second-nearest hopping integrals ℏ p γ 0 ¼ msa ¼ 3:036 eV; 0 2
γ 00 ¼
2ℏ2 u1 9pq 1 þ þ 2 ¼ 0:4530 eV: 2 s 2s 9ma0
ð22Þ
This provides a new way to determine TB parameters [2,8,11]. If k o 0:1 b (k 109 m 1 or less), the fourth-order terms in β and γ can be neglected and one has α ¼ 1:000, 3 2 2 pv2 9τ ℏ2 k þ 4ms β ¼ ℏm 1 þ us1 þ 9pq k sin 3θ; 2 qu2 þ 2 2s2
γ¼
9ℏ4 p2 2 9τℏ4 pu2 3 k þ k sin 3θ: 4m2 s2 4m2 s2
ð23Þ
100
J. Luo / Solid State Communications 195 (2014) 95–101
The factor sin 3θ in the k terms thus leads to the triangular warping of the energy dispersion. The expansion of the TB energy dispersion involving only nearest hopping integrals leads to the 2 factor sin 3θ in the k term [2]. This work indicates that more 3 accurately, this factor is in k term, which more accords with firstprinciples results [12]. Fig. 3 illustrates the triangular warping effect and the asymmetry of the positive and negative eigenenergies. 3
- -
It is interesting to compare this work with the usual k U p method [4]. If, instead of Eq. (8), one assumes - -
-
-
ψ j ð r ; k Þ ¼ ψ j ð R M ; k Þ;
- -
-
-
∇ψ j ð r ; k Þ ¼ ∇ψ j ð R M ; k Þ;
- -
∇∇ψ j ð r ; k Þ ¼ 0;
ð24Þ
then coefficients in Eq. (13) are zero except Dx and Dy. Eq. (12) reduces to Dirac equation ^ 3ℏp - σ U p Ψ ¼ εΨ ; 2ms
ð25Þ
^ where σ ¼ σ x x þ τσ y y and p ¼ ð iℏÞð x ∂=∂x þ y ∂=∂yÞ. Both Eqs. (20) and (25) lead to Fermi velocity
vF ¼
3ℏp ¼ 9:826 105 ms 1 ; 2ms
ð26Þ
the TB method to avoid empirical parameters but still retain its low computational cost.
5. Conclusions More accurate energy dispersion of graphene electrons can be obtained from a second-order envelope-function equation as an improvement of Dirac equation. DFT KS orbitals of isolated carbon atoms can play an important role in TB method by acting as Wannier functions, which endows KS orbitals with new significance. For graphene the improved TB method can totally avoid experimental data except the crystal structure. Further improvement is expected so that it can be applied to other systems.
Appendix: More details The symmetry of the orbital ϕ simplifies the calculation of integrals (9). For an atom M one has R R - 2 3R 2 331 ϕM d r ¼ 1; 1 ϕM ∇ϕM d r ¼ 0; 1 r M ϕM d r ¼ 0; Z ---3ϕM r M ∇ϕM d r ¼ u0;1 x x þu0;1 y y þ u0;3 z z ; 1 Z -- --3ϕ2M r M r M d r ¼ v0;1 x x þv0;1 y y þv0;3 z z ; ðA1Þ 1
which is in good agreement with the experimental value [10,11]. In Dirac equation, vF is defined according to the nearest momentum integral, that is, vF ¼ 3ℏp1 =2m. We note that p ¼ p1 p3 contains the third-nearest integral p3 , and s ¼ 1 3s2 the second-nearest one s2 . According to Eq. (15), p3 is comparable with p1 and 3s2 with s0 ¼ 1. If p3 and s2 are neglected, one has 3ℏp1 =2m ¼ 7:026 105 ms 1 . This explains the difference between the experimental value of vF and its definition according to the nearest momentum integral [4,6]. Even to determine vF from Dirac equation, at least the third-nearest integrals are necessary. More exact velocity of carriers can be calculated from -
the energy dispersion by v ¼ ð1=ℏÞ∇- ε, and the velocity difference of k
electrons and holes can be obtained from the asymmetry of the - -
energy dispersion. More accurate results in this work than the k U p method are achieved by expanding ψ 1 and ψ 2 to the second-order terms in Eq. (8) and by involving up to third-nearest orbital integrals. All numerical results in this work are based on the integrals of carbon 2pz orbitals, and no experimental data are involved except graphene's crystal structure. This is realized for two reasons. Firstly, the introduction of the envelope functions eliminates unknown crystal potential in Eq. (7). This is achieved due to graphene's unique crystal structure. Because of its π =3 rotational symmetry and the same carbon 2pz orbital for both A and B atoms, ^ by adjusting the zero point of the energy one has Hψ A ¼ 0 and ^ Hψ B ¼ 0 approximately and simultaneously. Secondly, the carbon 2pz orbital function is obtained from DFT calculation. There has been controversy about the significance of DFT KS orbitals. This work demonstrates that for carbon atoms in graphene, KS orbitals can be regarded as real single-electron wave functions in atomic systems and be used in TB wave functions. In the case of graphene electrons they constitute a good approximation of Wannier functions [18]. In this work, atomic orbitals are no longer symbols in TB method only for theoretical derivation. Instead they are used to numerically calculate the eigen-energy and envelope functions. Further improvement of this work may be the treatment of the crystal potential V. If, like the wave function, V can also be constructed in terms of atomic data calculated from DFT, this treatment may hopefully be applied to other systems such as semiconductors with direct band gaps [19–22]. Since atomic data can be calculated out as standard ones, this may ultimately enable
and for two nth-nearest atoms M and N with n ¼ 1; 2; 3, one has R 31 ϕM ϕN d r ¼ sn ; Z Z 33ϕM ∇ϕN d r ¼ pn t MN ; r M ϕM ϕN d r ¼ qn t MN ; 1 1 Z -3ϕM r M ∇ϕN d r ¼ un;1 t MN t MN þ un;2 n MN n MN þ un;3 z z ; 1 Z -- 3r M r M ϕM ϕN d r ¼ vn;1 t MN t MN þ vn;2 n MN n MN þ vn;3 z z : ðA2Þ 1
According to graphene's crystal structure, one can calculate out the following sums. For an atom A and its three nearest atoms B, one has - R 3∑ ei k F U ð R B R A Þ 1 ϕA ϕB d r ¼ 0; B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ∇ϕB d r ¼ eiπ=6 p1 ð x iτ y Þ; 2 1 B Z 3i 3∑ ei k F U ð R B R A Þ r A ϕA ϕB d r ¼ eiπ=6 q1 ð x iτ y Þ; 2 1 B Z 3∑ ei k F U ð R B R A Þ ϕA r A ∇ϕB d r
1
B
- 3i ¼ eiπ=6 ðu1;1 u1;2 Þð x þ iτ y Þð x þ iτ y Þ; 4 Z -
-
-
∑ ei k F U ð R B R A Þ B
¼
- -
1
3-
ϕA ϕB r A r A d r
- 3i iπ=6 e ðv1;1 v1;2 Þð x þ iτ y Þð x þiτ y Þ: 4
For an atom B and its three nearest atoms A, one has - R 3∑ ei k F U ð R A R B Þ 1 ϕB ϕA d r ¼ 0; A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ∇ϕA d r ¼ e iπ=6 p1 ð x þ iτ y Þ; 2 1 A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ϕA r B d r ¼ e iπ=6 q1 ð x þ iτ y Þ; 2 1 A Z 3i k F U ð R A R BÞ ϕB r B ∇ϕA d r ∑e
A
1
- 3i ¼ e iπ=6 ðu1;1 u1;2 Þð x iτ y Þð x iτ y Þ; 4
ðA3Þ
J. Luo / Solid State Communications 195 (2014) 95–101 -
-
-
∑ ei k F U ð R A R B Þ
Z
- -
1
A
-
3-
For an atom A and its six second-nearest atoms A , one has -
-
∑e A
0
-
-
-
∑ei k F U ð R A0 R A Þ A0
-
-
-
-
-
-
Z
0
i k F U ð R A0 R A Þ
∑e A
1 ϕA ϕA
0
0
1
ϕA ∇ϕA0 d r ¼ 0;
Eqs. (A3)–(A8) indicate that Pauli matrixes enter the envelope
-
3-
--
--
-
-
B
-
-
-
∑ei k F U ð R B0 R B Þ B0
-
-
-
∑ei k F U ð R B0 R B Þ B0
-
-
-
∑ei k F U ð R B0 R B Þ B0
R
1 ϕB ϕB
0
ðA5Þ
1
1
Z
1
at
3-
-
3-
1 --
2
--
--
--
¼ 32 ½ðv2;1 þ v2;2 Þð x x þ y y Þ þ 2v2;3 z z :
ðA6Þ
-
B
- -
1
3-
ϕA ϕB r A r A d r
- 3i ¼ eiπ=6 ðv3;1 v3;2 Þð x þ iτ y Þð x þ iτ y Þ: 4
ðA7Þ
For an atom B and its three third-nearest atoms A, one has - R 3∑ ei k F U ð R A R B Þ 1 ϕB ϕA d r ¼ 0; A Z 3i 3∑ ei k F U ð R A R B Þ ϕB ∇ϕA d r ¼ e iπ=6 p3 ð x þ iτ y Þ; 2ℏ 1 A Z 3i 3∑ ei k F U ð R A R B Þ ϕA ϕB r B d r ¼ e iπ=6 q3 ð x þ iτ y Þ; 2 1 A Z 3∑ ei k F U ð R A R B Þ ϕB r B ∇ϕA d r
A
-
RB
according
to
1
- 3i ¼ e iπ=6 ðu3;1 u3;2 Þð x iτ y Þð x iτ y Þ; 4
2
2
2
- 3iℏ2 p 3iℏ2 u2 ∇ψ 1 U ð x þτi y Þ þ ∇∇ψ 1 : ð x τi y Þð x iτ y Þ 2m 4m v1 v3 ∂2 ψ ¼ sεψ 2 þ ε∇22 ψ 2 þ ε 22 4 2 ∂z - 3iq 3iv2 ε∇ψ 1 U ð x þ iτ y Þ ε∇∇ψ 1 : ð x iτ y Þð x iτ y Þ; þ 2 8
ðA10Þ
References
1
-
at
where ψ 1 , ψ 2 and their derivatives take values at R B .
- 3i ¼ eiπ=6 ðu3;1 u3;2 Þð x þ iτ y Þð x þ iτ y Þ; 4 Z -
series
-
- R 3∑ ei k F U ð R B R A Þ 1 ϕA ϕB d r ¼ 0; B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ∇ϕB d r ¼ eiπ=6 p3 ð x iτ y Þ; 2ℏ 1 B Z 3i 3∑ ei k F U ð R B R A Þ ϕA ϕB r A d r ¼ eiπ=6 q3 ð x iτ y Þ; 2 1 B Z 3∑ ei k F U ð R B R A Þ ϕA r A ∇ϕB d r
∑ ei k F U ð R B R A Þ
power
For an atom A and its three third-nearest atoms B, one has
B
into
u1 2 ∇2 ψ 2 ℏ mu3 ∂∂zψ22 ℏ2ms∇2 ψ 2 ℏ2m
¼ 32 ½ðu2;1 þ u2;2 Þð x x þ y y Þ þ 2u2;3 z z ; Z - 3∑ei k F U ð R B0 R B Þ ϕB ϕB0 r B r B d r B0
-
ϕB ð r Þ and conducts the triple integration. The expan-
sion of ψ 1 , ψ 2 Eq. (8) leads to
ϕB r B ∇ϕB0 d r --
-
e
e
3-
--
R A . For a chosen atom B, one multiplies Eq. (7) by
iπ=12 i k F U R B
ϕB ϕB0 r B d r ¼ 0; -
-
-
ϕB ∇ϕB0 d r ¼ 0;
Z
ðA9Þ
where ∇22 ¼ ∂2 =∂x2 þ ∂2 =∂y2 , ψ 1 , ψ 2 and their derivatives take values
3-
d r ¼ 3s2 ;
Z
2
- 3iℏ2 p 3iℏ2 u2 ∇ψ 2 U ð x τi y Þ ∇∇ψ 2 : ð x þ τi y Þð x þ τi y Þ 2m 4m v1 v3 ∂2 ψ ¼ sεψ 1 þ ε∇22 ψ 1 þ ε 21 4 2 ∂z - 3iq 3iv2 ε∇ψ 2 U ð x τi y Þ þ ε∇∇ψ 2 : ð x þτi y Þð x þ τi y Þ; þ 2 8
For an atom B and its six second-nearest atoms B , one has -
2
0
0
2
2
¼ 32 ½ðv2;1 þ v2;2 Þð x x þ y y Þ þ 2v2;3 z z :
∑ei k F U ð R B0 R B Þ
-
u1 2 ∇2 ψ 1 ℏ mu3 ∂∂zψ21 ℏ2ms∇2 ψ 1 ℏ2m
--
--
1
--
-
-
¼ 32 ½ðu2;1 þ u2;2 Þð x x þ y y Þ þ 2u2;3 z z ; Z - 3∑ei k F U ð R A0 R A Þ ϕA ϕA0 r A r A d r A0
-
expansion of ψ 1 , ψ 2 into power series at R A according to Eq. (8) leads to
ϕA r A ∇ϕA0 d r --
ðA8Þ
e iπ=12 e i k F U R A ϕA ð r Þ and conducts the triple integration. The
3ϕA ϕA0 r A d r ¼ 0; -
1
1
3-
-
d r ¼ 3s2 ; 3-
1
- -
ϕA ϕB r B r B d r
equation due to the vectors x 7 iτ y , and their factors e 7 iπ =6 must be eliminated by e 7 iπ =12 in Eq. (4). For a chosen atom A, one multiplies Eq. (7) by
3-
Z Z
∑ei k F U ð R A0 R A Þ A
R
Z
- 3i ¼ e iπ=6 ðv3;1 v3;2 Þð x iτ y Þð x iτ y Þ: 4
ðA4Þ 0
-
-
A
- 3i ¼ e iπ=6 ðv1;1 v1;2 Þð x iτ y Þð x iτ y Þ: 4
i k F U ð R A0 R A Þ
-
∑ ei k F U ð R A R B Þ
ϕB ϕA r B r B d r
101
[1] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [2] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109. [3] L. Brey, H.A. Fertig, Phys. Rev. B 73 (2006) 235411. [4] D.P. DiVincenzo, E.J. Mele, Phys. Rev. B 29 (1984) 1685. [5] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197. [6] P.R. Wallace, Phys. Rev. 71 (1947) 622. [7] R. Saito, G. Dresselhaus, M.S. Dresselhaus, Physical Properties of Carbon Nanotubes, Imperial, London, 1998. [8] R. Saito, M. Fujita, G. Dresselhaus, M.S. Dresselhaus, Appl. Phys. Lett. 60 (1992) 2204. [9] J.W. McClure, Phys. Rev. 104 (1956) 666. [10] M.S. Dresselhaus, G. Dresselhaus, Adv. Phys. 51 (2002) 1. [11] R.S. Deacon, K.C. Chuang, R.J. Nicholas, K.S. Novoselov, A.K. Geim, Phys. Rev. B 76 (2007) 081406(R). [12] S. Reich, J. Maultzsch, C. Thomsen, P. Ordejon, Phys. Rev. B 66 (2002) 035412. [13] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [14] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [15] B. Delley, J. Chem. Phys. 92 (1990) 508; B. Delley, J. Chem. Phys. 113 (2000) 7756. [16] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [17] C. Kittel, Introduction to Solid State Physics, Eighth edition, John Wiley & Sons Inc., New York, 2005. [18] G. Wannier, Phys. Rev. 52 (1937) 191. [19] J.P. Perdew, M. Levy, Phys. Rev. Lett. 51 (1983) 1884. [20] L.J. Sham, M. Schluter, Phys. Rev. Lett. 51 (1983) 1888. [21] A.J. Cohen, P. Mori-Sanchez, W. Yang, Phys. Rev. B 77 (2008) 115123. [22] M.K.Y. Chan, G. Ceder, Phys. Rev. Lett. 105 (2010) 196403.