ar
X
iv
:1
00
5.
35
92
v1
[
co
nd
-m
at.
me
s-h
all
]
20
M
ay
20
10
Accurate evaluation of the Green’s function of disordered
graphenes
W. Zhu1,2, Q. W. Shi1,2†, X. R. Wang2,3∗, X. P. Wang1, J. L. Yang1, J. Chen4,5, J. G. Hou1
1Hefei National Laboratory for Physical Sciences at Microscale,
University of Science and Technology of China, Hefei 230026, China
2Department of Physics, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong
3School of Physics, Shangdong University, Jinan, P. R. China
4Electrical and Computer Engineering,
University of Alberta, Alberta, Canada T6G 2V4 and
5National Research Council/National Institute
for Nanotechnology, Alberta, Canada T6G 2M9∗
(Dated: May 21, 2010)
Abstract
An accurate simulation of Green’s function and self-energy function of non-interacting electrons
in disordered graphenes are performed. Fundamental physical quantities such as the elastic relax-
ation time τe, the phase velocity vp, and the group velocity vg are evaluated. New features around
the Dirac point are revealed, showing hints that multi-scattering induced hybridization of Bloch
states plays an important role in the vicinity of the Dirac point.
PACS numbers: 81.05.ue, 71.55.-i, 71.23.-k
∗Electronic address:phsqw@ustc.edu.cn;phxwan@ust.hk
1

Introduction. −Graphene, a single layer of graphite, has been intensively studied in recent
years because of many intriguing transport properties [1, 2]. Examples include minimum
conductivity and linear carrier density dependence of conductivity [3, 4]. Electrons in an
ideal graphene are governed by the relativistic massless Dirac equation and exhibit a linear
dispersion relation in the vicinity of the Dirac point and zero density of state at the Dirac
point [5, 6]. Among many relativistic effects, the Klein paradox [7] is arguably one of the
most important effects that differentiates Dirac electrons [8] from the Schrodinger electrons
in disordered systems.
A deep understanding of Dirac electrons in disordered graphenes requires extended knowl-
edge on the self-energy function in order to extract such fundamental physical quantities as
the phase velocity vp, the group velocity vg, and the elastic relaxation time τe. However, it is
known [9] that accurate and reliable calculations are quite difficult and nontrivial. For this
reason, various approximations have been employed in different theoretical studies. So far,
almost all calculations [10, 11] concerning electron properties in disordered graphenes were
performed without fully considering disorder effects. The wave-nature of Dirac electrons is
more pronounced near the Dirac point because of very large electron wavelength there. It is
known that interference due to multi-scattering leads to weak localization and the Ander-
son localization in conventional disordered electron systems[12]. Quantum interference also
plays important roles in coherent wave propagation through quasi-random [13] and random
media[14, 15]. As it will be shown below, interference-induced the hybridization of Bloch
states is essential in understanding the diffusion properties near the Dirac point in realistic
disordered graphenes, although previous calculations [10, 11] only captured the essential
physics far away from the Dirac point.
In this Letter, we present a systematic method to exactly calculate Green’s function and
the self-energy of large size disordered graphenes. We extract accurate τe, vp, and vg values
from the spectral function A(k, E) derived from the self-energy function. In comparison
with the results from the self-consistent Born approximation (SCBA), it is found that τe
is overestimated by the SCBA in the strong disorder regime. We show that both vp and
vg deviate significantly from their unrenormalized values and exhibit substantial energy
dependence, especially near the Dirac point. The effective group velocity is larger than
the effective phase velocity, but significantly lower than its unrenormalized value when the
mixing of different Bloch states is dominant. Moreover, we generalize the Einstein relation
2

to calculate the conductivity of the disordered graphene.
Model-and-method. − pi-electrons of undoped graphene can be modeled by a tight-binding
Hamiltonian on a honeycomb lattice of two sites per unit cell, H0 = t
∑
<ij>
|i >< j| + h.c.,
where t = −2.7 eV is the hopping energy. The corresponding eigenvalues and eigenstates
of H0 near the Dirac point are [5, 17], respectively, Ek,± = ±~v0Fk and |k± >= (|kA >
±eiφ(k)|kB >)/2, where v0F is the unrenormalized Fermi velocity, and ~ is the Plank constant.
A and B stand for A- and B-sublattices. φ(k) is the polar angle of the momentum k,
and |kA(B) >= 1√NA(B)
∑
rA(B)
eik·rA(B) |rA(B) >, where rA(B) is the position vector of the
A(B)-lattice and NA(B) is the total A(B)-lattice points. The plus (minus) sign denotes the
conduction (valence) band. The Green function of clean graphene is, in a diagonal basis,
G0(k, E) = 1E+i0+−~v0F k |k+ >< k+ |+
1
E+i0++~v0F k
|k− >< k−|. A weak point-like disorder is
introduced through V =
∑
i|i >< i| with nimp randomly distributed impurity sites where
the on-site energy i of each impurity can take −V0 or V0 (measured in the unit of t) with
equal probability. A dimensionless parameter α = nimpV
2
0 Ac
2pi(~v0F )2
can be used to characterize
disorder strength. Here, Ac is the area of the unit cell. Our calculation shows that the
physical quantities such as self-energy are only determined by the parameter α[18].
The ensemble-averaged Green function is defined as G(k±, E) =
< k± | 1E+iη−H0−V |k± >, where the bar means the ensemble average. It can be cal-
culated by using the well developed Lanczos recursive method[19, 20]. In order to obtain
an accurate ensemble-averaged Green’s function near the Dirac point with high energy
resolution[21], a large sample containing N = Lx × Ly ' 6.0 millions carbon atoms
(2400×2400) is used in our simulation, where Lx(Ly) is the number of atoms in x(y)
directions. The periodic boundary condition is used in our simulation in order to reduce
the finite size effect. Thus, the wave vectors are kx = nx4pi/3aLx and ky = ny4pi/
√
3aLy,
where nx(y) is an integer and a is the lattice constant.
Self-energy function. − Fig. 1 shows the calculated real (a) and imaginary (b) parts of the
self-energy function for nimp/N = 10% and V0 = 0.5 (open squares) and 2.0 (open circles),
respectively. The self-energy function is defined as usual, where Σ(k, E) = G−10 (k, E) −
G−1(k, E) [23]. In principle, the self-energy function depends on energy E and wave vector
k. However, our simulation finds that the self-energy function is not sensitive to wave vector
k, while the one-particle Green’s function depends on both E and k. This is not surprising
3

since the scatterer size is much smaller than the electron wavelength near the Dirac point
so that the inhomogeneous structure of disordered graphene can be described well by an
effective homogeneous medium. This is also why the self-energy function is assumed to be
k-independent in many perturbative calculations. Our finding validates this assumption[23].
Our calculation should be compared with widely used results from the self-consistent
Born approximation (SCBA) that predicts [17, 24]
Σ(E) = {
−2E/α− i2Γ0, |E|  Γ0
−2α(E + ipiα|E|)ln|EcE | − ipiα|E|, |E|  Γ0
where Γ0 = Ece−1/α (Ec is the cut-off energy) As shown in Fig. 1, SCBA results agree well
with our exact self-energy function for very weak disorder (V0=0.5 and nimp/N=10%)[25].
When the disorder strength increases several times (V0=2.0 or α∼0.07), the perturbative
results can not capture the main features, especially near the Dirac point. The discrepancy
is obvious for ImΣ(0) as shown in Fig. 1(b). ImΣ(E) at the Dirac point is Γ0 ∼10−5t
(Ec =
√
3t [6]) from the SCBA, while our exact value is ImΣ(0) ∼10−2t. Thus, the true
broadening of states at or near the Dirac point is much greater than what has been predicted
by SCBA. This huge discrepancy can be attributed to the mixture of Bloch states caused
by the impurities. Furthermore, the level-repulsion effect pushes all energy level toward the
Dirac point so that the density of states at the Dirac point increases more in the presence of
impurities. Thus, the impurities make the imaginary part of the self-energy function (directly
associated with the density of states) at the Dirac point bigger. When the wavelength
becomes short and the quantum interference as well as the Bloch state mixing are less
important, ImΣ(E) are determined by the disorder scattering. The difference between our
exact simulation and that of the SCBA is small as shown in Fig. 1.
Spectral function. − The single-particle spectral function relates to Green’s function
through A(k±, E) = −ImG(k±, E)/pi. Fig. 2(a) is A(k+, E) for ky = 0, V0 = 1,
nimp/N=10%, and various kx (curves from the left to the right in the figure) ranging from
0.0 (nx=0) to 0.098 (nx=56) (in unit of a−1). In the absence of disorder, the spectral func-
tion A0(k±, E) is a delta function, reflecting that the wave vector k is a good quantum
number and has all its spectral weight precisely at the energy E = Ek±. In the presence of
disorder, the translational symmetry is broken and the spectral function is broadened due
to the disorder scattering effect. The widths of the spectral function are given by ImΣ(E)
4

-0.2 0.0 0.2
-40
-20
0
20
40
-0.2 0.0 0.2
-100
-80
-60
-40
-20
0
V
0
=0.5
V
0
=2.0
SCBA
SCBA
R
e
(
E
)

(
1
0
-
3
t
)
E(t)
(a)
(b)
E(t)

I
m
(
E
)

(
1
0
-
3
t
)




FIG. 1: (Color online) (a) Real part of self-energy as a function of energy. (b)Imaginary part of
self-energy as a function of energy. The open squares and open circles represent our numerical
calculations for the disorder concentration nimp/N = 10% and V0 = 0.5 and V0 = 2.0, respectively.
Dashed lines represent the SCBA results for the same disorder. Energy is measured in the units
of t.
that measures the elastic relaxation lifetime τe, where τe = ~−2ImΣ(E) . Therefore, the elastic
scattering relaxation time is akin to ImΣ(E), and τe around the Dirac point is mainly deter-
mined by the Bloch state mixing and the level repulsion effect[26]. Far away from the Dirac
point, the lifetime is mainly attributed to the disorder scattering. As shown in the inset of
Fig. 2, lifetime become shorter as the wave vector increases. Physically, this is because the
density of states ρ increases linearly with energy and the disorder scattering effects become
larger, which is qualitatively consistent with the prediction of SCBA[16, 17, 24].
Effective band velocity. − Dirac electron propagation velocities, including the group
velocity vg and phase velocity vp, are also greatly modified by the disorder effects. These
quantities relate to the shape of the dispersion relation that is the roots of E − E0(k) −
ReΣ(E) = 0[27, 28]. One can also extract the dispersion relation from the peak of the
spectral function A(k, E) for a given k. Fig 3 (a) shows our exact dissipation curve Eeff(k)
that is linear for very weak disorder (α < 0.01). However, when the disorder strength
becomes significantly large (α ∼ 0.1), the Eeff (k) becomes concave near the Dirac point,
indicating the reduction in both the group and phase velocities.
Fig. 3(b) shows k−dependence of the group and phase velocities obtained from vg =
∂Eeff (k)/∂k and vp = Eeff/k. For the very weak disorder α = 0.01 as shown in the figure,
5

0.0 0.1 0.2
0
500
1000
1500
0.00 0.01 0.02 0.03 0.04
1E-3
0.01
0.1
1
10
A
(
k
,
E
)
E(t)
n
imp
=10%
V
0
=1.0


V
0
=0.5
V
0
=1.0
V
0
=1.5
V
0
=2.0
E(t)

e
(
p
s
)
FIG. 2: (Color online) Single-particle spectral function A(~k+, E) plotted as a function of energy E
at several k-points (from left to right): k = 0.000, 0.014, 0.028, 0.042, 0.056, 0.070, 0.084, 0.098
(or nx = 0 ∼ 56) along the kx direction. The model parameters are nimp/N = 10% and V0 =
1.0. Inset: The energy dependence of single-particle relaxation time τe for V0 = 0.5 ∼ 2.0 and
nimp/N = 10%.
vg and vp are not too much different from the unrenomalized velocity v0F . Their values
are reduced by 5% in comparison with v0F . When the disorder strength increase several
times α = 0.07, the renormalized vg is higher than vp and both vg and vp are reduced by a
large percentage near the Dirac point. This shows that disorder not only renormalizes the
Dirac electron velocities, but also changes the linear dispersion relation. The fact of large
reductions of velocities at the Dirac point indicates that Dirac electrons near the Dirac point
are more sensitive to the disorder.
One can use the renormalization factor Z defined as Z = (1− ∂ReΣ(E,k)∂E )−1 [29] to measure
the effect of disorder on electronic structure. Its value equals the ratio vg/v0F . As shown
in Fig. 3(b), Z is very close to 1.0 for very weak disorder. In this regime, the Bloch state
is still a good starting point for understanding disordered graphene, while the transport
properties are expected to be described by the quasi-classical Boltzmann theory. When the
disorder strength increases several times, calculations shows that Z is much smaller than 1.0,
especially around the Dirac point. This unusual feature directly reflects that a Bloch states
around the Dirac point couples strongly with other nearby Bloch states. Therefore, the Bloch
states would not be a good starting point to perform perturbative calculations. In fact, the
similar feature has been observed in semiconductor alloys and has been termed non-Bloch
nature of alloy states[30]. This is why the SCBA method cannot produce accurate enough
6

-0.10 -0.05 0.00 0.05 0.10
0
5
10
15
20
25
30
35
40
V
0
=0.5
V
0
=1.0
V
0
=2.0
(
4
e
2
/
h
)
E (t)


-0.2 0.0 0.2
0.00
0.03
0.06
D
O
S
E(t)

FIG. 4: (Color online) Conductivity as a function of charge density. The model parameters are
nimp/N = 10% and V0 =0.5, 1.0 and 2.0. Inset: Density of states for the same parameters.
calculations find that σxx(0) takes the values 8.2e2/h, 5.9e2/h and 2.7e2/h for α = 0.0046,
0.018 and 0.07, respectively, and shows non-universal behavior.
Summary. − In conclusion, we studied point-like disorder effects on the one-electron
properties of graphene. The exact ensemble-averaged Green’s function is obtained from a
large-scale real-space calculation. Through the analysis of self-energy and spectral functions,
we conclude that the single-particle lifetime reduction and the linear dispersion relation are
modified by the hybridization of the Bloch states. Furthermore, we studied the diffusion
transport properties by using our exact self-energy and the Einstein relation. Our approach
is very general and robust, thus is applicable to many other disordered systems.
Acknowledgment.− This work is partially supported by NNSF of China (Nos.
10974187,10874165, and 50721091), by NKBRP of China under Grant No. 2006CB922000,
and by KIP of the Chinese Academy of Sciences (No. KJCX2-YW-W22). JC is supported
by the NRC and NSERC of Canada (No. 245680). XRW acknowledges the support of Hong
Kong RGC grants (#604109, RPC07/08.SC03, and HKU10/CRF/08-HKUST17/CRF/08).
[1] K. S. Novoselov, et. al, Science 306, 666 (2004)
[2] Y. Zhang, et. al, Nature, 438, 201 (2005).
[3] A. K. Geim, et. al, Nat. Mater. 6, 183 (2007).
[4] S. Das Sarma, et. al, Solid State Commun. 143, 1(2007).
[5] A. H. Castro Neto, et. al, Rev. Mod. Phys. 81,109 (2009).
8

[6] Z.F. Wang, et. al, Phys. Rev. B, 74, 125417(2006).
[7] M. Katsnelson, K. Novoselov, and A. Geim, Nat. Phys. 2, 620 (2006).
[8] Y.Y. Zhang, et. al, Phys. Rev. Lett. 102, 106401 (2009); Phys. Rev. B 78, 155413 (2008).
[9] E. N. Economou, Green’s Functions in Quantum Physics, Springer (2006).
[10] E. H. Hwang and S. Das Sarma, Phys. Rev. B, 77, 195412 (2008).
[11] M. I. Katsnelson, et. al, Phil. Trans. R. Soc. A 366, 195(2008).
[12] P. A. Lee, et. al, Rev. Mod. Phys. 57, 287 (1985).
[13] X. R. Wang, Phys. Rev. B 51, 9310 (1995); 53, 12035 (1996).
[14] M. C. W. van Rossum, et. al, Rev. Mod. Phys. 71, 313 (1999).
[15] M. P. van Albada, et. al, Phys. Rev. Lett. 66, 3132 (1991).
[16] K. Normura and A. H. MacDonald, ibid. 96, 256602 (2006); ibid. 98, 076602 (2007).
[17] T. Ando, J. Phys. Soc. Jpn. 75, 074716 (2006); N. H. Shon and T. Ando, J. Phys. Soc. Jpn.
67, 2421 (1998).
[18] Other uncorrelated short range impurity potential (such as Anderson disorder) produces sim-
ilar results.
[19] J. K. Cullum and R. A. Willoughby, Lancozs Algorithms for large symmetric eigenvalue prob-
lem, Birkhauser(1985); C. Lanczos, J. Res. Natl. Bur. Stand. 45, 255 (1950); R. Haydock, V.
Heine, and M. J. Kelly, J. Phys. C 5, 2845 (1972); S. D. Wu, et al., Phys. Rev. B 77, 195411
(2008).
[20] E. Dagotto, Rev. Mod. Phys. 66, 763 (1994); L. C. Davis, Phys. Rev. B 28, 6961 (1983); B.
Bauml, et.al ibid. 58, 3663 (1998); S. M. Anaage, et.al, ibid. 34, 2336 (1998).
[21] A small artificial cut-off energy η is used to simulate the infinitesimal imaginary energy in
Green function. This artificial cut-off energy may be viewed as a measure of electron-phonon
or electron-electron interactions in reality. We chose η = 0.01meV in this paper, which is much
less than the energy resolution in recent ARPES experiments [22]. See also W. Zhu et.al, Phys.
Rev. Lett. 102, 056803 (2009).
[22] A. Damascelli, et. al, Rev. Mod. Phys. 75, 473(2003).
[23] A. Altland and B. Simons, Condensed Matter Field Theory, Cambridge (2006).
[24] P. M. Ostrosky, et. al, Phys. Rev. B 74, 235443(2006).
[25] T. Fukuzawa, et. al, J. Phys. Soc. Jpn. 78, 094714(2009).
[26] M. Calandra and F. Mauri, Phys. Rev. B 76, 205411 (2007).
9

[27] G. D. Mahan, Many-Particle Physics, third edition, Kluwer Academic, New York (2000).
[28] K. Tarafder, et. al, Phys. Rev. B 74, 144204 (2006).
[29] C. H. Park, et. al, Phys. Rev. Lett. 99, 086804 (2007).
[30] L.-W. Wang, et. al, ibid. 80, 4725(1997); Y. Zhang, et. al, ibid. 101, 036403 (2008).
[31] E.A. Henriksen, et. al, ibid. 104, 067404 (2010); E.A. Henriksen (private communication).
[32] T. Ando, et. al, J. Phys. Soc. Jpn. 71, 1318 (2002).
10