One-parameter scaling at the Dirac point in graphene
Bardarson, J.H.; Tworzydlo, J.; Brouwer, P.W.; Beenakker, C.W.J.
Citation
Bardarson, J. H., Tworzydlo, J., Brouwer, P. W., & Beenakker, C. W. J. (2007). One-parameter
scaling at the Dirac point in graphene. Physical Review Letters, 99(10), 106801.
doi:10.1103/PhysRevLett.99.106801
Version: Not Applicable (or Unknown)
License: Leiden University Non-exclusive license
Downloaded from: https://hdl.handle.net/1887/71386
Note: To cite this publication please use the final published version (if applicable).
One-Parame ter Scaling at the Dirac Point in Graphene
J. H. Bardarson,
1
J. Tworzydło,
2
P. W. Brouwer,
3,4
and C. W. J. Beenakker
1
1
Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands
2
Institute of Theoretical Physics, Warsaw University, Hoz
˙
a 69, 00–681 Warsaw, Poland
3
Physics Department, Arnold Sommerfeld Center for Theoretical Physics, Ludwig-Maximilans-Universita
¨
t, 80333 Munich, Germany
4
Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853, USA
(Received 7 May 2007; published 5 September 2007)
We numerically calculate the conductivity of an undoped graphene sheet (size L) in the limit of a
vanishingly small lattice constant. We demonstrate one-parameter scaling for random impurity scattering
and determine the scaling function d ln=d lnL. Contrary to a recent prediction, the scaling flow
has no fixed point (>0) for conductivities up to and beyond the symplectic metal-insulator transition.
Instead, the data support an alternative scaling flow for which the conductivity at the Dirac point increases
logarithmically with sample size in the absence of intervalley scattering— without reaching a scale-
invariant limit.
DOI: 10.1103/PhysRevLett.99.106801 PACS numbers: 73.20.Fz, 73.20.Jc, 73.23.b, 73.63.Nm
Graphene provides a new regime for two-dimensional
quantum transport [1–3], governed by the absence of
backscattering of Dirac fermions [4]. A counterintuitive
consequence is that adding disorder to a sheet of undoped
graphene initially increases its conductivity [5,6].
Intervalley scattering at stronger disorder strengths enables
backscattering [7], eventually leading to localization and to
a vanishing conductivity in the thermodynamic limit [8,9].
Intervalley scattering becomes less and less important if
the disorder is more and more smooth on the scale of the
lattice constant a. The fundamental question of the new
quantum transport regime is how the conductivity scales
with increasing system size L if intervalley scattering is
suppressed.
In usual disordered electronic systems, the hypothesis of
one-parameter scaling plays a central role in our concep-
tual understanding of the metal-insulator transition
[10,11]. According to this hypothesis, the logarithmic de-
rivative d ln=d lnL is a function only of itself
[12]—irrespective of the sample size or degree of disorder.
A positive function means that the system scales towards
a metal with increasing system size, while a negative
function means that it scales towards an insulator. The
metal-insulator transition is at 0,
0
> 0. In a two-
dimensional system with symplectic symmetry, such as
graphene, one would expect a monotonically increasing
function with a metal-insulator transition at [13]
S
1:4 (see Fig. 1, green dotted curve).
Recent papers have argued that graphene might deviate
in an interesting way from this simple expectation. Nomura
and MacDonald [14] have emphasized that the very exis-
tence of a function in undoped graphene is not obvious,
in view of the diverging Fermi wave length at the Dirac
point. Assuming that one-parameter scaling does hold,
Ostrovsky, Gornyi, and Mirlin [15] have proposed the
scaling flow of Fig. 1 (black solid curve). Their function
implies that approaches a universal, scale-invariant value
in the large-L limit, being the hypothetical quantum
critical point of a certain field theory. This field theory
differs from the symplectic sigma model by a topological
term [15,16]. The quantum critical point could not be
derived from the weak-coupling theory of Ref. [15], but
its existence was rather concluded from the analogy to the
effect of a topological term in the field theory of the
quantum Hall effect [11,17]. The precise value of
is
therefore unknown, but it is well constrained [15]: from
below by the ballistic limit
0
1= [18–20] and from
above by the unstable fixed point
S
1:4.
In this Letter we present a numerical test first, of the
existence of one-parameter scaling, and second of the
FIG. 1 (color online). Two scenarios for the scaling of the
conductivity with sample size L at the Dirac point in the
absence of intervalley scattering. The black solid curve with two
fixed points is proposed in Ref. [15], the green dotted curve
without a fixed point is an alternative scaling supported by the
numerical data presented in this Letter. For comparison, we
include as a red dashed curve the scaling flow in the symplectic
symmetry class, which has a metal-insulator transition at
S
1:4 [13].
PRL 99, 106801 (2007)
PHYSICAL REVIEW LETTERS
week ending
7 SEPTEMBER 2007
0031-9007=07=99(10)=106801(4) 106801-1 © 2007 The American Physical Society
scaling prediction of Ref. [15] against an alternative scal-
ing flow, a positive without a fixed point (green dotted
curve in Fig. 1). For such a test it is crucial to avoid the
finite-a effects of intervalley scattering that might drive the
system to an insulator before it can reach the predicted
scale-invariant regime. We accomplish this by starting
from the Dirac equation, being the a ! 0 limit of the
tight-binding model on a honeycomb lattice. We have
developed an efficient transfer operator method to solve
this equation, which we describe before proceeding to the
results.
The single-valley Dirac Hamiltonian reads
H vp VxUx; y: (1)
The vector of Pauli matrices acts on the sublattice index
of the spinor , p i@@=@r is the momentum operator,
and v is the velocity of the massless excitations. The
disorder potential Ur varies randomly in the strip 0 <
x<L, 0 <y<W (with zero average, hUi0). This dis-
ordered strip is connected to highly doped ballistic leads,
according to the doping profile Vx0 for 0 <x<L,
Vx!1for x<0 and x>L. We set the Fermi energy
at zero (the Dirac point), so that the disordered strip is
undoped. The disorder strength is quantified by the corre-
lator
K
0
1
@v
2
Z
dr
0
hUrUr
0
i: (2)
Following Refs. [5,21], we work with a transfer operator
representation of the Dirac equation H 0 at zero en-
ergy. We discretize x at the N points x
1
;x
2
; ...;x
N
and
represent the impurity potential by Ur
P
n
U
n
yx
x
n
. Upon multiplication by i
x
the Dirac equation in the
interval 0 <x<Ltakes the form
@ v
@
@x
x
y
vp
y
z
i
x
X
n
U
n
yx x
n
x
y:
(3)
The transfer operator M, defined by
L
M
0
, is given
by the operator product
M P
L;x
N
K
N
P
x
N
;x
N1
K
2
P
x
2
;x
1
K
1
P
x
1
;0
; (4)
P
x;x
0
exp1=@x x
0
p
y
z
; (5)
K
n
expi=@vU
n
x
: (6)
The operator P gives the decay of evanescent waves
between two scattering events, described by the operators
K
n
. For later use we note the current conservation relation
M
1
x
M
y
x
: (7)
We assume periodic boundary conditions in the
y direction, so that we can represent the operators in the
basis
k
1
W
p
e
iq
k
y
ji;q
k
2k
W
;k 0; 1; 2... :
(8)
The spinors ji 2
1=2
1
1
, ji 2
1=2
1
1
are eigenvec-
tors of
x
. In this basis, p
y
kk
0
@q
k
kk
0
is a diagonal
operator, while U
n
kk
0
W
1
R
dyU
n
yexpiq
k
0
q
k
y is nondiagonal. We work with finite-dimensional
transfer matrices by truncating the transverse momenta
q
k
at jkjM.
The transmission and reflection matrices t, r are deter-
mined as in Ref. [19], by matching the amplitudes of
incoming, reflected, and transmitted modes in the heavily
doped graphene leads to states in the undoped strip at
x 0 and x L. This leads to the set of linear equations
X
k
kk
0
k
yr
kk
0
k
y
0
y; (9a)
X
k
t
kk
0
k
y
L
yM
0
y: (9b)
Using the current conservation relation (7) we can solve
Eq. (9) for the transmission matrix,
1 r
1 r
M
y
t
t
) t
1
hjM
y
ji: (10)
The transmission matrix determines the conductance G
4e
2
=hTrtt
y
, and hence the dimensionless conductivity
h=4e
2
L=WG. The average conductivity hi is
obtained by sampling some 10
2
–10
3
realizations of the
impurity potential.
Because the transfer matrix P has both exponentially
small and exponentially large eigenvalues, the matrix mul-
tiplication (4) is numerically unstable. As in Ref. [22], we
stabilize the product of transfer matrices by transforming it
into a composition of unitary scattering matrices, involving
only eigenvalues of unit absolute value.
We model the disorder potential Ur
P
N
n1
n
x
x
n
y y
n
by a collection of N isolated impurities dis-
tributed uniformly over a strip 0 <x<L, 0 <y<W. (An
alternative model of a continuous Gaussian random poten-
tial is discussed at the end of the paper.) The strengths
n
of
the scatterers are uniform in the interval [
0
,
0
]. The
number N sets the average separation d WL=N
1=2
of
the scatterers. The cutoff jkjM imposed on the trans-
verse momenta q
k
limits the spatial resolution
W=2M 1 of plane waves / e
iq
k
yq
k
x
at the Dirac point.
The resulting finite correlation lengths of the scattering
potential in the x and y directions scale with , but they are
not determined more precisely. The disorder strength (2)
evaluates to K
0
1
3
2
0
@vd
2
, independent of the corre-
lation lengths. We scale towards an infinite system by
increasing M / L at fixed disorder strength K
0
, scattering
range =d, and aspect ratio W=L.
This completes the description of our numerical method.
We now turn to the results. In Fig. 2 we first show the
PRL 99, 106801 (2007)
PHYSICAL REVIEW LETTERS
week ending
7 SEPTEMBER 2007
106801-2
dependence of the average conductivity on K
0
for a fixed
system size. As in the tight-binding model of Ref. [6],
disorder increases the conductivity above the ballistic
value. This impurity assisted tunneling [5] saturates in an
oscillatory fashion for K
0
1 (unitary limit [23,24]). In
the tight-binding model [6] the initial increase of was
followed by a rapid decay of the conductivity for K
0
* 1,
presumably due to Anderson localization. The present
model avoids localization by eliminating intervalley scat-
tering from the outset.
The system size dependence of the average conductivity
is shown in Fig. 3, for various combinations of disorder
strength and scattering range. We take W=L sufficiently
large that we have reached an aspect-ratio independent
scaling flow and L=d large enough that the momentum
cutoff M>25. The top panel shows the data sets as a
function of L=d. The increase of with L is approximately
logarithmic, hiconst 0:25 lnL, much slower than the
L
p
increase obtained in Ref. [5] in the absence of mode
mixing.
If one-parameter scaling holds, then it should be pos-
sible to rescale the length L
fK
0
; =dL such that the
data sets collapse onto a single smooth curve when plotted
as a function of L
=d. (The function f d=l
determines
the effective mean free path l
, so that L
=d L=l
.) The
bottom panel in Fig. 3 demonstrates that, indeed, this data
collapse occurs. The resulting function is plotted in the
inset. Starting from the ballistic limit [20]at
0
1=, the
function first rises until 0:6, and then decays to zero
without becoming negative. For >
S
1:4 the decay
/ 1= is as expected for a diffusive system in the sym-
plectic symmetry class. The positive function in the
interval (
0
,
S
) precludes the flow towards a scale-
invariant conductivity predicted in Ref. [15].
The model of isolated impurities considered so far is
used in much of the theoretical literature, whereas experi-
mentally a continuous random potential is more realistic
[14]. We have therefore also performed numerical simula-
tions for a random potential landscape with Gaussian
correlations [25],
hUrUr
0
i K
0
@v
2
2
2
e
jrr
0
j
2
=2
2
: (11)
The discrete points x
1
;x
2
; ...;x
N
in the operator product
(4) are taken equidistant with spacing x L=N, and
U
n
y
Z
x
n
x=2
x
n
x=2
dxUx; y: (12)
We take M, N, and W=L large enough that the resulting
conductivity no longer depends on these parameters. We
then scale towards larger system sizes by increasing L=
and W= at fixed K
0
. No saturation of with increasing K
0
is observed for the continuous random potential (as ex-
pected, since the unitary limit is specific for isolated scat-
0.3
0.6
0.9
1.2
1.5
1.8
1 10 100
〈σ〉 [4e
2
/h]
L/d
ξ/d
K
0
2.0 1.5 1.0 0.5 0.25
0.16
0.20
0.25
0.4
0.6
0.8
1.0
2.0
4.3
7.0
7.5
10
14
0.3
0.6
0.9
1.2
1.5
1.8
1 10 100 10
3
〈σ〉 [4e
2
/h]
L*/d
σ
β
0.3
0.0
1.50.5
FIG. 3 (color online). System size dependence of the average
conductivity, for W=L 4 (black and green or dark gray solid
symbols) and W=L 1:5 (all other symbols) and various com-
binations of K
0
and =d. The top panel shows the raw data. In
the bottom panel the data sets have been given a horizontal
offset, to demonstrate the existence of one-parameter scaling.
The inset shows the resulting function.
0.2
0.4
0.6
0.8
1
1.2
0.01 0.1 1 10 100 1000
〈σ〉 [4e
2
/h]
K
0
ξ/d
0.25
0.50
1.00
1.50
FIG. 2. Disorder strength dependence of the average conduc-
tivity for a fixed system size (W 4L 40d) and four values of
the scattering range.
PRL 99, 106801 (2007)
PHYSICAL REVIEW LETTERS
week ending
7 SEPTEMBER 2007
106801-3
terers [23,24]). Figure 4 shows the size dependence of the
conductivity—both the raw data as a function of L (inset)
as well as the rescaled data as a function of L
gK
0
L.
Single-parameter scaling applies for L * 5, where hi
const 0:32 lnL. The prefactor of the logarithm is about
25% larger than in the model of isolated impurities (Fig. 3),
which is within the numerical uncertainty.
In conclusion, we have demonstrated that the central
hypothesis of the scaling theory of quantum transport, the
existence of one-parameter scaling, holds in graphene. The
scaling flow which we find (green dotted curve in Fig. 1)is
qualitatively different both from what would be expected
for conventional electronic systems (red dashed curve) and
also from what has been predicted [24] for graphene (black
solid curve). Our scaling flow has no fixed point, meaning
that the conductivity of undoped graphene keeps increas-
ing with increasing disorder in the absence of intervalley
scattering. The fundamental question ‘‘what is the limiting
conductivity
1
of an infinitely large undoped carbon
monolayer’’ has therefore three different answers:
1
1= in the absence of any disorder [18,19],
1
1with
disorder that does not mix the valleys (this Letter), and
1
0 with intervalley scattering [8,9].
We thank C. Mudry and M. Titov for valuable
discussions. This research was supported by the Dutch
Science Foundation NWO/FOM, the European
Community’s Marie Curie Research Training Network
(Contract No. MRTN-CT-2003-504574, Fundamentals of
Nanoelectronics), and by the Packard Foundation.
Note added.—Since submission of this manuscript,
similar conclusions have been reported by Nomura,
Koshino, and Ryu [26].
[1] A. H. Castro Neto, F. Guinea, and N. M. R. Peres, Phys.
World 19, No. 11, 33 (2006).
[2] A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183
(2007).
[3] M. I. Katsnelson and K. S. Novoselov, Solid State
Commun. 143, 3 (2007).
[4] T. Ando, T. Nakanishi, and R. Saito, J. Phys. Soc. Jpn. 67,
2857 (1998).
[5] M. Titov, Europhys. Lett. 79, 17004 (2007).
[6] A. Rycerz, J. Tworzydło, and C. W. J. Beenakker,
Europhys. Lett. 79, 57003 (2007).
[7] E. McCann, K. Kechedzhi, V. I. Fal’ko, H. Suzuura, T.
Ando, and B. L. Altshuler, Phys. Rev. Lett. 97, 146805
(2006).
[8] I. L. Aleiner and K. B. Efetov, Phys. Rev. Lett. 97, 236801
(2006).
[9] A. Altland, Phys. Rev. Lett. 97, 236802 (2006).
[10] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57,
287 (1985).
[11] K. Efetov, Supersymmetry in Disorder and Chaos
(Cambridge University, Cambridge, 1997).
[12] We define the function in terms of the ensemble
averaged conductivity , measured in units of 4e
2
=h
(with the factor of 4 accounting for twofold spin and
valley degeneracies). This is the appropriate definition
for our system. For a more general definition of one-
parameter scaling, one needs to scale a distribution func-
tion of conductances [11].
[13] P. Markos and L. Schweitzer, J. Phys. A 39, 3221 (2006).
[14] K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 98,
076602 (2007).
[15] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Phys. Rev.
Lett. 98, 256801 (2007).
[16] S. Ryu, C. Mudry, H. Obuse, and A. Furusaki, arXiv:cond-
mat/0702529 [Phys. Rev. Lett. (to be published)].
[17] A. W. W. Ludwig, M. P. A. Fisher, R. Shankar, and G.
Grinstein, Phys. Rev. B 50, 7526 (1994).
[18] M. I. Katsnelson, Eur. Phys. J. B 51, 157 (2006).
[19] J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz, and
C. W. J. Beenakker, Phys. Rev. Lett. 96, 246802 (2006).
[20] We call
0
the ballistic limit because it is reached in the
absence of disorder, but we emphasize that it is a con-
ductivity—not a conductance. This is a unique property
(called ‘‘pseudodiffusive’’) of graphene at the Dirac point,
that its conductance scales / 1=L like in a diffusive
system even in the absence of disorder.
[21] V. V. Cheianov and V. I. Fal’ko, Phys. Rev. B 74,
041403(R) (2006).
[22] H. Tamura and T. Ando, Phys. Rev. B 44, 1792 (1991).
[23] D. V. Khveshchenko, Phys. Rev. Lett. 97, 036802 (2006).
[24] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Phys. Rev.
B 74, 235443 (2006).
[25] The Dirac equation with a delta-function correlated ran-
dom potential has a divergent scattering rate, see, e.g.,
A. A. Nersesyan, A. M. Tsvelik, and F. Wenger, Nucl.
Phys. B438, 561 (1995). Hence the need to regularize
the continuous potential model by means of a finite
correlation length .
[26] K. Nomura, M. Koshino, and S. Ryu, arXiv:0705.1607.
0
0.5
1
1.5
2
1 10 100 10
3
〈σ〉 [4e
2
/h]
L*/ξ
K
0
1.0
2.25
4.0
12.25
L/ξ
0
1
2
1 10 100
FIG. 4 (color online). System size dependence of the average
conductivity in the continuous potential model, for several
values of K
0
. The inset shows the raw data, while the data sets
in the main plot have a horizontal offset to demonstrate one-
parameter scaling when L * 5.
PRL 99, 106801 (2007)
PHYSICAL REVIEW LETTERS
week ending
7 SEPTEMBER 2007
106801-4
![FIG. 1 (color online). Two scenarios for the scaling of the conductivity with sample size L at the Dirac point in the absence of intervalley scattering. The black solid curve with two fixed points is proposed in Ref. [15], the green dotted curve without a fixed point is an alternative scaling supported by the numerical data presented in this Letter. For comparison, we include as a red dashed curve the scaling flow in the symplectic symmetry class, which has a metal-insulator transition at S 1:4 [13].](/figures/fig-1-color-online-two-scenarios-for-the-scaling-of-the-246cp5fm.webp)
