- Email: [email protected]
Impurity levels in the electronic spectra of graphene
Superlattices and Microstructures 53 (2013) 55–62
Contents lists available at SciVerse ScienceDirect
Superlattices and Microstructures journal homepage: www.elsevier.com/locate/superlattices
Impurity levels in the electronic spectra of graphene A. Feher a, S. Feodosyev b, I. Gospodarev b, O. Kotlyar b,⇑, K. Kravchenko b, E. Manzhelii b, E. Syrkin b Institute of Physics, Faculty of Science, P.J. S˘afárik University, Park Angelinum 9, SK-04154 Kosice, Slovakia B.Verkin Institute for Low Temperature Physics and Engineering of National Academy of Sciences of Ukraine, 47 Lenin Ave., 61103 Kharkov, Ukraine
a
b
a r t i c l e
i n f o
Article history: Received 8 June 2012 Received in revised form 10 September 2012 Accepted 14 September 2012 Available online 13 October 2012 Keywords: Graphene Electronic spectra Impurity Local discrete level
a b s t r a c t We calculate and analyze, within a model that allows to obtain an analytical approximation using the method of Jacobi matrices, the electronic spectra of graphene with impurity atoms, especially with boron substitutional impurity. We present analytical expressions for the conditions of the existence of electronic local discrete levels due to the presence of a substitutional impurity in graphene. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction It is well known that graphene monolayers cannot exist as planar objects in the free state, because in the flat 2D-crystals the mean-square displacements of the atoms in the direction normal to the plane of the layer diverge even at T ¼ 0 (see, e.g., [1]). So we can study and practically apply only such a graphene, which is deposited on a certain substrate providing the stability of the plane graphene monolayers (see, e.g., [2–4]). Only small graphene flakes can be detached from the substrate and these flakes immediately acquire a corrugated shape [5]. For theoretical and experimental studying of the electronic properties of graphene various substrates are used (see, e.g., [6–8]). The presence of the substrate greatly increases the possibility to introduce various defects into graphene. Impurity atoms embedded in graphene may cause the appearance of impurity states outside the band of quasi-continuous spectrum. At low impurity concentrations (when impurity is considered as an isolated defect) these states appear in the form of local discrete levels (LDLs). Although such
⇑ Corresponding author. Tel.: +380 686143327; fax: +380 573403370. E-mail address: [email protected] (O. Kotlyar). 0749-6036/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.spmi.2012.09.004
56
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
levels in various quasiparticle spectra have been known and studied over 60 years, an adequate description is still absent, even in the harmonic approximation for sufficiently realistic models of the crystal lattice (even with an isolated defect). The dependence of appearance conditions and characteristics of the LDL on the parameters of a perfect lattice and defect was identified only in the most general terms. However, LDL been essential part of local density of states (LDOSs) may be used as source of information about the defect structure and force interactions in real crystals. Modern experimental methods (e.g. scanning tunneling microscopy (STM) technique) give a possibility to obtain LDOS of different atoms in graphene and ultra-thin carbon nanofilms [8,9]. To extract such information it is useful to have analytical expressions that relate the main characteristics of LDL with parameters of the defect and the host lattice. In this paper we present the results of our calculations and analysis of the characteristics of the electronic LDL for the substitutional impurities in graphene, especially for boron substitutional impurity, using analytical approximation based on the Jacobi matrices method. 2. Electronic spectrum of graphene with defects Graphene is a two-dimensional semiconductor with zero band gap. The fact that the charge carriers in graphene are formally described by the Dirac equation and not by the Schrödinger equation, is due to the symmetry of the crystal lattice of graphene, which consists of two equivalent carbon sublattices. Electronic subbands formed by the symmetric and antisymmetric combinations of wave functions in the two sublattices intersect at the edge of the Brillouin zone, which leads to a cone-shaped energy spectrum near the K and K 0 points of the first Brillouin zone. The electrons obey the linear dispersion law (in ordinary metals and semiconductors the dispersion law is parabolic). The electronic spectrum of graphene can be described by a strong coupling approximation, and it is sufficient to consider the interaction between nearest neighbors only (see, for example, [3,4,10–12]). The corresponding Hamiltonian is
b ¼ H
X
ei jiihij
i
X J ij jiihjj;
ð1Þ
i;j
where i and j are the labels of the nodes of the two-dimensional lattice, ei is the energy of an electron at node i, and J ij is the so-called overlap integral. Curve 1 in Fig. 1 shows the density of electronic states of graphene. These calculations are made by using the method of Jacobi matrices [13,14]. In a perfect graphene the local Green’s function (LGF)
Fig. 1. Electronic density of states of perfect graphene (curve 1) and the local density of states for an isolated boron substitutional impurity (curve 2). e0 ¼ 3J is the half-width of the quasi-continuous spectrum.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
57
^ 1 jii; Gðe; iÞ ¼ hij eI^ H (I^ is the identity operator) is coincide with the total
GðeÞ ¼ lim
N 1X ^ 1 jii hij eI^ H
N!1 N
i¼1
because of the physical equivalence of the atoms of both sublattices. Peculiarity of the density of states at e ¼ eðKÞ (the value eðKÞ corresponds to the Fermi energy eF in graphene) determines the behavior of the real part of the Green’s function near eF . For a wide class of perturbations caused by defects we can find by using the Lifshitz equation [15], quasilocalized states in the interval ½eðMÞ; eðMÞ (where M is the point of the first Brillouin zone). This equation, which determines the energy of these states, can be written as (see, e.g., [1,14])
ReGðeÞ ¼ Sðe; Kik Þ;
ð2Þ
b (Kik are matrix elements of where the Sðe; Kik Þ function is determined by the perturbation operator K this operator on defined basis). The local spectral densities .ðe; iÞ p1 Imlimc#0 Gðe ic; iÞ of impurity atoms are calculated in [4]. For an isolated substitutional impurity, different from the host lattice atom, with the energy of the impurity node ~e and the overlap integral J i0 ¼ eJ ð1 þ gÞJ, the function Sðe; ~e; gÞ has the form
Sðe; ~e; gÞ ¼
ð1 þ gÞ2 : ~e þ egð2 þ gÞ
ð3Þ
For an nitrogen impurity (according to [4]) ~e eðKÞ 0:525J and g ¼ 0:5. As shown in [16], the Eq. (2) has a solution for both interval ½eðMÞ; eðKÞ and interval ½eðKÞ; eðMÞ. The local density of states of the nitrogen substitutional impurity calculated in [4] has quasi-local maxima in both intervals. For an boron substitutional impurity g ¼ 0:5 [4], quasi-localized states are absent in interval ½eðMÞ; eðMÞ [4,17]. Fig. 2 shows the graphical solution of the Lifshitz Eq. (2) for a given impurity atom. In this case also the Lifshitz equation has no solutions in interval ½eðMÞ; eðMÞ (corresponding to the dependence SðeÞ, curves 2 in Fig. 2). The local density of states of the boron impurity (curve 2 in Fig. 1) has two ðÞ poles outside the band of quasi-continuous spectrum ed , which are called local discrete levels and which are also solutions of Eq. (2). As is clearly seen in Fig. 1 the area under the curve 2 is smaller (by the sum of the residues at these poles) than the area under the curve 1. Local discrete levels can be an important source of information about defective structure and force interactions in real crystals. To extract this useful information we should have analytical expressions that relate the main characteristics of LDL (primarily their energy) to the parameters of the defect and the host lattice.
Fig. 2. Graphical solution of Eq. (2) for boron substitutional impurity in graphene. Curve 1 is the real part of the Green function of ideal graphene, curve 10 is the ‘‘approximations of two moments’’ of this function. Curve 2 represents function SðeÞ for boron impurity (~e ¼ 0:525J; g ¼ 0:5).
58
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
Such expressions were obtained in [17] for localized vibrations in the phonon spectrum of a threedimensional crystal. In [17] authors proposed the analytical approximation of the basic characteristics of local vibrations based on the rapid convergence of the real part of the Green’s function outside the band of quasi-continuous spectrum using the method of Jacobi matrices [13,14]. 3. Analytical approximation of the characteristics of local discrete levels Let us, briefly, to the extent necessary to understand the use of the classification of the eigenfunctions of Hamiltonian (1), we present the basics of the method of Jacobi matrices. This method allows, without finding the dispersion laws, to calculate directly the local partial Green’s functions of the system, corresponding to the perturbation of one or more atoms. This disturbance is described by a socalled generating vector ~ h0 2 H, where H is the space of electronic excitations of atoms. Its dimension is qN, where N is the number of atoms in the system, and q is the dimension of the displacement of a P single atom ðq ¼ 1; 2; 3Þ. If, using the generating vector ~ h0 ¼ p1 p jji (p is the number of excited j¼1
b n~ h0 g1 atoms) and the Hamiltonian (1) we construct the sequence f H n¼0 , then the linear envelope covering the vectors of this sequence forms, in the H space, a cyclic subspace invariant to the operator b This subspace contains, within itself, all the atomic displacements generated by the vector ~ H. h0 . The
corresponding partial Green’s function is determined as a matrix element 1 b n~ b h0 g1 e Ib H j0i. In the basis f H n¼0 which is obtained by orthonormalization of
G00 ðeÞ GðeÞ ¼ h0j
the sequence the operator (1) is represented in the form of a tridiagonal Jacobi matrix (or J-matrix). This matrix has a simple spectrum, what greatly simplifies finding the partial Green’s functions and spectral densities. As can be seen, this method does not use explicitly the translational symmetry of the crystal, what makes it extremely effective for treating systems in which such symmetry is broken. The method of Jacobi matrices is particularly effective for treating systems with a simply connected quasi-continuous band of spectrum D. In this case, with increasing rank of the J-matrix ðn ! 1Þ its diagonal elements an converge to a corresponding to the middle of the bandwidth D, and nondiagonal elements bn converge to b corresponding to one-quarter of the bandwidth D. For the local Green’s function corresponding to the excitations of one or more atoms, which are determined by the generating vector ~ h0 , we get following expression using the J-matrix method
Qn ðeÞ bn1 Qn1 ðeÞK1 ðeÞ ^ 1~ h0 Þ ¼ lim Gðe; ~ h0 Þ ¼ ð~ h0 ; ½eI^ H ; n!1 P n ðeÞ bn1 P n1 ðeÞK1 ðeÞ
ð4Þ
where polynomials P n ðeÞ are determined by the following recurrence relation
bn P nþ1 ðeÞ ¼ ðe an ÞP n ðeÞ bn P n1 ðeÞ:
ð5Þ
Polynomials Qn ðeÞ determined by the same recurrence relation. The initial conditions is 1 P 1 ðeÞ ¼ Q0 ðeÞ ¼ 0; P 0 ðeÞ ¼ 1; Q1 ðeÞ ¼ b0 . K1 ðeÞ is the function corresponds to the LGF operator, with all elements of its Jmatrix being equal to their limit values a and b
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 K1 ðeÞ ¼ 2b ½e a ZðeÞ ðe a þ 2bÞðe a 2bÞ;
ð6Þ
ZðeÞ 0ða 2b eÞ þ i Hðe a þ 2bÞHða þ 2b eÞ Hðe a 2bÞ:
ð7Þ
The method of Jacobi matrices can treat, as a regular singular perturbation, a much larger number of perturbations of the quasi-particle spectrum due to the presence of various crystal defects, than traditional methods [13,14]. In addition, perturbations do not change the bandwidth of the quasi-continuous spectrum, and consequently, the asymptotic values of the elements of the J-matrix can be regarded as an asymptotically degenerate regular perturbation [17]. This type of perturbations covers virtually all perturbations of the quasi-particle spectrum caused by local defects. The calculations of dynamics and thermodynamics characteristics of such systems are performed, using the method of J-matrices, with the same accuracy as for the initial ideal system.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
59
In practice, it is usually possible to calculate the Jacobi matrix of the Hamiltonian of a finite rank. The expression
Qn ðeÞ bn1 Qn1 ðeÞK1 ðeÞ Gðe; ~ h0 Þ ¼ P n ðeÞ bn1 P n1 ðeÞK1 ðeÞ
ð8Þ
is called analytical approximation of LGF. Quantity qðeÞ p1 ImG00 ðeÞ is called the spectral density generated by the initial displacement ~ h0 . All dependencies in Fig. 1 and curves 1 in Fig. 2 were calculated by the formula (8), using the Jacobi matrix of the Hamiltonian (1) with rank n ¼ 600. If we count the energy from the Fermi energy level, then all diagonal elements of Jacobi matrices are zero (8an ¼ a ¼ 0; b ¼ e0 =2). A good accuracy of the approximations shown in figures is confirmed also by the fact that they show the nonanalyticity effects corresponding to the densities of states of systems with the dimension higher than unity (so-called van Hove singularities). In the vicinity of these singularities the expression (8) slowly converges to the true values of the real and imaginary parts of LGF. Curve 10 in Fig. 2 show the real parts of the LGF calculated by formula (8) with n ¼ 1. As can be seen in the band of the quasi-continuous spectrum these relationships have very little in common with curve 1. On the curve 10 in the interval ½eðMÞ; eðMÞ are absent non-monotonous parts and the logarithmic singularities at the edges of the band of quasi-continuous spectrum, characteristic for the 2D systems. However, outside the band of quasi-continuous spectrum (also in the area of intersection of the real part of LGF with curve 2) curves 1 and 10 practically coincide. If in the Lifshitz equation (2) we put instead of LGF its approximation (8) for n ¼ 1, the obtained solutions give the energies of LDL with quite high accuracy. Moreover, these solutions can be easily found analytically. The LGF approximation by formula (8) for n ¼ 1 was named as the approximation of two moments [17]. Indeed, from the orthonormality of the polynomials [13,14], it follows that
Z
e0
P1 ðeÞ.ðeÞde ¼ 0 ) a0 ¼
e0
Z
e0
e0
P21 ðeÞ.ðeÞde ¼
Z
Z
e0
e.ðeÞde ¼ M1 ;
e0
e0
e2 a20 2
b0
e0
.ðeÞde ¼ 1 ) b0 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M 2 M21 :
Finding characteristics of LDL is more convenient without using Eq. (2), looking for them as the e e; h~0 Þ ¼ ðh~0 ; ½ebI H b K ^ 1 h~0 Þ. In the approximation of poles of the LGF perturbed Hamiltonian Gð two moments for the subspace generated by the excitation of an impurity atom, we get
e e; ~ Gð h0 Þ ¼
1
e a0 b20 K 1 ðeÞ
ð9Þ
:
For the case of an isolated substitutional impurity
a0 ¼ ~e;
b0 ¼
pffiffiffi 1þg 3ð1 þ gÞJ ¼ pffiffiffi e0 ; 3
ð10Þ
from where we get
e e; ~ h0 Þ ¼ Gð
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 cÞe a0 þ ZðeÞc je2 e20 j RðeÞ
;
ð11Þ
where RðeÞ ¼ ð1 2cÞe2 2ð1 cÞa0 e þ a20 þ c2 e20 and
c b20 =2b2 ¼ 2ð1 þ gÞ2 =3:
ð12Þ
Local discrete levels are poles (9), i.e. the roots of RðeÞ are ðÞ d
e
¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc 1Þa0 c a20 þ ð2c 1Þe20 2c 1
:
ð13Þ
60
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
ðÞ e eÞ are called the intensities of LDL and they determine the se¼eðÞ Gð Residues at these poles l0 ¼ re d relative LDL ‘‘amplitude’’ on the impurity atom. The condition that the intensity differs from zero defines the existence region of LDL. In this case
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ca0 þ ðc 1Þ a20 þ ð2c 1Þe20
a0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hð 1 þ cÞ; e0 2 2 ð2c 1Þ a0 þ ð2c 1Þe0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ca0 þ ðc 1Þ a20 þ ð2c 1Þe20 a0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Hðc 1 Þ: e0 2 2 ð2c 1Þ a0 þ ð2c 1Þe0
ðþÞ 0
l
¼
ðÞ 0
l
ð14Þ
In [17] it was shown that
b 1~ hn Þ ¼ P m ðeÞQn ðeÞ þ P m ðeÞP n ðeÞGðe; ~ Gmn ðe; ~ h0 Þ ¼ ðh~m ; ½ebI H h0 Þ:
ð15Þ
This implies that the damping of LDL, i.e. the decay of its intensity with increasing distance from the impurity atom (i.e. with the increase of n) follows the equation ðÞ 2 ðÞ 2 ðÞ e e lðÞ n ¼ rese¼eðÞ G nn ðeÞ ¼ P n ðed Þrese¼eðÞ GðeÞ ¼ P n ðed Þl0 : d
ð16Þ
d
Using the method of mathematical induction we can prove that
ðÞ d Þ
P n ðe
2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3n a20 þ ð2c 1Þe20 a0 pffiffiffiffiffiffi 5 : ¼ 2c4 e0 ð2c 1Þ
ð17Þ
and
l
ðÞ n>0
ðÞ n
¼ 2c½q
l
ðÞ 0 ;
qðÞ
2qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 32 a20 þ ð2c 1Þe20 a0 5 6 1: ¼4 e0 ð2c 1Þ
ðÞ
ð18Þ
Thus the intensities ln form with increasing n an infinitely decreasing geometric progression with ðþÞ ðÞ denominators qðÞ 6 1. Value qðÞ ¼ 1 if ed ¼ e0 or ed ¼ e0 . Summing these progressions, we see P1 ðÞ that l ¼ 1, that is the formation of each LDL is the formation of one quasiparticle outside n¼0 n the band of the quasi-continuous electron spectrum.
Fig. 3. Regions of the existence of discrete levels for substitutional impurities in graphene.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
61
Fig. 4. The dependencies on g of basic characteristics of LDL on the substitutional impurity in graphene. The value ~e ¼ 0:525J corresponds to boron [4].
4. Analysis of the characteristics of the local discrete level on the boron subsitutional impurity Formulas (10)–(14) give simple analytical expressions for the conditions for the existence of local discrete levels due to the presence of a substitutional impurity in graphene. Regions of the LDL exisðþÞ ðÞ tence are for ed (in this case ~e > e0 ð1 cÞ ¼ e0 ½1 2ð1 þ gÞ2 =3) and for ed (in this case ~e < e0 ðc 1Þ ¼ e0 ½2ð1 þ gÞ2 =3 1) shown in Fig. 3. ~ and g. The absence of LDL is posIt is seen that such levels exist in a very widep range ffiffiffiffiffiffiffiffi of variables peffiffiffiffiffiffiffiffi sible only in a narrow range of values ~e at g < 3=2 1. If g > 3=2 1 at least one local discrete level exists. Realistically, the lines delineating the area of existence of LDL in graphene in the plane ð~e; gÞ must pass through the origin of coordinates, since ReGðeÞ ! 1 for e ! e0 . However, since this divergence is logarithmic, for any appreciable splitting of LDL from the boundary of quasi-continuous spectrum there is a certain threshold. Curves 1 and 10 in Fig. 2 merge at j e j> e0 . Fig. 4 shows, for the value ~e 0:525J (for boron, from [4]), the dependences of energies, the intensities of LDL at the boron impurity and the damping parameters of the value g that characterizes the change in the overlap integral of the boron impurity atom. The solid lines show the characteristics of LDL, calculated using the approximation of the Green’s function (9), i.e. according to the analytical formulas (13), (14) and (18). The open circles correspond to the results of the numerical calculations of ðÞ ðÞ dependences ed ðgÞ and l0 ðgÞ, using the Green’s function in the form of (8), calculated by the Jacobi matrices of the n ¼ 600 rank. ðÞ It is seen that at the threshold values ed ¼ e0 of the LDL formation, intensities vanish, and the ðÞ parameters of damping are equal to the unit. Further increase of j ed j is accompanied by the increase ðÞ of l0 and by the strengthening of the containment level.
5. Conclusions Good agreement between the results of numerical calculation of the LDL characteristics using Jacobi matrices of high rank, and their analytical description by the Green’s function (9) gives the possibility of extraction of the defect parameters from the known characteristics of LDL. Experimental measurement (e.g. by scanning tunneling microscopy) of values should allow, using (13) and (14), the determination of the parameters and this might lead to a significant advance in creating nanomaterials with predetermined spectral characteristics. As follows from (14) and can be seen from Fig. 4, ðþÞ ðÞ with increasing the intensity of LDL l0 ! l0 ! 1=2. That is, the impurity levels cannot be completely localized on the impurities, but they also appear in the spectra of carbon atoms surrounding an impurity atom. This greatly increases the probability of experimental detection of such levels, even
62
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
at low concentrations of impurities. Experimental study of the local discrete levels in graphene with boron impurities using STM techniques at temperatures near 0.3 K is in progress now. Acknowledgments This work is the result of the project implementation: Research and Education at UPJŠ – Heading towards Excellent European Universities, ITMS project code: 26110230056, supported by the Operational Program Education funded by the European Social Fund (ESF). References [1] A.M. Kossevich, The Crystal Lattice (Phonons Solitons Dislocations), Wiley VCH Verlag Berlin GmBH, Berlin, 1999. [2] K.S. Novoselov, A.K. Gein, S.V. Morozov, D. Diang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197–201. [3] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109–155. [4] N.M.R. Peres, F.D. Klironomos, S-W. Tsai, J.R. Santos, J.M.B. Lopes dos Santos, A.H. Castro Neto, EPL 80 (2007) 67007–67017. [5] J.C. Meyer, A.K. Geim, M.I. Katsnelson, K.S. Novoselov, T.J. Booth, S. Roth, Nature 446 (2007) 60–63. [6] L. Zhao, K.T. Rim, H. Zhou, R. He, T.F. Heinz, A. Pinczuk, G.W. Flynn, A.N. Pasupathy, Solid State Comm. 151 (2011) 509–513. [7] E. Stolyarova, K.T. Rim, S. Ryu, J. Maultzsch, P. Kim, L.E. Brus, T.F. Heinz, M.S. Hybertsen, G.W. Flynn, PNAS 104 (2007) 9209– 9212. [8] J. Xue, J. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, B.J. LeRoy, Phys. Rev. Lett. 108 (2012) 016801–1016801-5. [9] L. Zhao, R. He, K.T. Rim, T. Schiros, K.S. Kim, H. Zhou, C. Gutiérrez, S.P. Chockalingam, C.J. Arguello, L. Pálova, D. Nordlund, M.S. Hybertsen, D.R. Reichman, T.F. Heinz, P. Kim, A. Pinczuk, G.W. Flynn, A.N. Pasupathy, Science 333 (2011) 999–1003. [10] Yu.V. Skrypnyk, V.M. Loktev, Phys. Rev. B 73 (2006) 241402-1–241402-4. [11] Yu.V. Skrypnyk, V.M. Loktev, Low Temp. Phys. 34 (2008) 818–825. [12] C. Bena, S.A. Kivelson, Phys. Rev. B. 72 (2005) 125432-1–125432-7. [13] V.I. Peresada, Condens. Matter Phys. 4 (1968) 172–210. [14] V.I. Peresada, V.N. Afanas’ev, V.S. Borovikov, Sov. J. Low Temp. Phys. 1 (1975) 227–232. [15] I.M. Lifshits, Dokl. Akad. Nauk SSSR (in Russian) 48 (1945) 83–86. [16] A. Feher, E. Syrkin, S. Feodosyev, I. Gospodarev, K. Kravchenko, in: Physics and Applications of Graphene: Theory (ed. S. Mikhailov), InTech, Rijeka, 2011, pp. 93-112. [17] O.V. Kotlyar, S.B. Feodosyev, Low Temp. Phys. 32 (2006) 256–269.
Contents lists available at SciVerse ScienceDirect
Superlattices and Microstructures journal homepage: www.elsevier.com/locate/superlattices
Impurity levels in the electronic spectra of graphene A. Feher a, S. Feodosyev b, I. Gospodarev b, O. Kotlyar b,⇑, K. Kravchenko b, E. Manzhelii b, E. Syrkin b Institute of Physics, Faculty of Science, P.J. S˘afárik University, Park Angelinum 9, SK-04154 Kosice, Slovakia B.Verkin Institute for Low Temperature Physics and Engineering of National Academy of Sciences of Ukraine, 47 Lenin Ave., 61103 Kharkov, Ukraine
a
b
a r t i c l e
i n f o
Article history: Received 8 June 2012 Received in revised form 10 September 2012 Accepted 14 September 2012 Available online 13 October 2012 Keywords: Graphene Electronic spectra Impurity Local discrete level
a b s t r a c t We calculate and analyze, within a model that allows to obtain an analytical approximation using the method of Jacobi matrices, the electronic spectra of graphene with impurity atoms, especially with boron substitutional impurity. We present analytical expressions for the conditions of the existence of electronic local discrete levels due to the presence of a substitutional impurity in graphene. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction It is well known that graphene monolayers cannot exist as planar objects in the free state, because in the flat 2D-crystals the mean-square displacements of the atoms in the direction normal to the plane of the layer diverge even at T ¼ 0 (see, e.g., [1]). So we can study and practically apply only such a graphene, which is deposited on a certain substrate providing the stability of the plane graphene monolayers (see, e.g., [2–4]). Only small graphene flakes can be detached from the substrate and these flakes immediately acquire a corrugated shape [5]. For theoretical and experimental studying of the electronic properties of graphene various substrates are used (see, e.g., [6–8]). The presence of the substrate greatly increases the possibility to introduce various defects into graphene. Impurity atoms embedded in graphene may cause the appearance of impurity states outside the band of quasi-continuous spectrum. At low impurity concentrations (when impurity is considered as an isolated defect) these states appear in the form of local discrete levels (LDLs). Although such
⇑ Corresponding author. Tel.: +380 686143327; fax: +380 573403370. E-mail address: [email protected] (O. Kotlyar). 0749-6036/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.spmi.2012.09.004
56
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
levels in various quasiparticle spectra have been known and studied over 60 years, an adequate description is still absent, even in the harmonic approximation for sufficiently realistic models of the crystal lattice (even with an isolated defect). The dependence of appearance conditions and characteristics of the LDL on the parameters of a perfect lattice and defect was identified only in the most general terms. However, LDL been essential part of local density of states (LDOSs) may be used as source of information about the defect structure and force interactions in real crystals. Modern experimental methods (e.g. scanning tunneling microscopy (STM) technique) give a possibility to obtain LDOS of different atoms in graphene and ultra-thin carbon nanofilms [8,9]. To extract such information it is useful to have analytical expressions that relate the main characteristics of LDL with parameters of the defect and the host lattice. In this paper we present the results of our calculations and analysis of the characteristics of the electronic LDL for the substitutional impurities in graphene, especially for boron substitutional impurity, using analytical approximation based on the Jacobi matrices method. 2. Electronic spectrum of graphene with defects Graphene is a two-dimensional semiconductor with zero band gap. The fact that the charge carriers in graphene are formally described by the Dirac equation and not by the Schrödinger equation, is due to the symmetry of the crystal lattice of graphene, which consists of two equivalent carbon sublattices. Electronic subbands formed by the symmetric and antisymmetric combinations of wave functions in the two sublattices intersect at the edge of the Brillouin zone, which leads to a cone-shaped energy spectrum near the K and K 0 points of the first Brillouin zone. The electrons obey the linear dispersion law (in ordinary metals and semiconductors the dispersion law is parabolic). The electronic spectrum of graphene can be described by a strong coupling approximation, and it is sufficient to consider the interaction between nearest neighbors only (see, for example, [3,4,10–12]). The corresponding Hamiltonian is
b ¼ H
X
ei jiihij
i
X J ij jiihjj;
ð1Þ
i;j
where i and j are the labels of the nodes of the two-dimensional lattice, ei is the energy of an electron at node i, and J ij is the so-called overlap integral. Curve 1 in Fig. 1 shows the density of electronic states of graphene. These calculations are made by using the method of Jacobi matrices [13,14]. In a perfect graphene the local Green’s function (LGF)
Fig. 1. Electronic density of states of perfect graphene (curve 1) and the local density of states for an isolated boron substitutional impurity (curve 2). e0 ¼ 3J is the half-width of the quasi-continuous spectrum.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
57
^ 1 jii; Gðe; iÞ ¼ hij eI^ H (I^ is the identity operator) is coincide with the total
GðeÞ ¼ lim
N 1X ^ 1 jii hij eI^ H
N!1 N
i¼1
because of the physical equivalence of the atoms of both sublattices. Peculiarity of the density of states at e ¼ eðKÞ (the value eðKÞ corresponds to the Fermi energy eF in graphene) determines the behavior of the real part of the Green’s function near eF . For a wide class of perturbations caused by defects we can find by using the Lifshitz equation [15], quasilocalized states in the interval ½eðMÞ; eðMÞ (where M is the point of the first Brillouin zone). This equation, which determines the energy of these states, can be written as (see, e.g., [1,14])
ReGðeÞ ¼ Sðe; Kik Þ;
ð2Þ
b (Kik are matrix elements of where the Sðe; Kik Þ function is determined by the perturbation operator K this operator on defined basis). The local spectral densities .ðe; iÞ p1 Imlimc#0 Gðe ic; iÞ of impurity atoms are calculated in [4]. For an isolated substitutional impurity, different from the host lattice atom, with the energy of the impurity node ~e and the overlap integral J i0 ¼ eJ ð1 þ gÞJ, the function Sðe; ~e; gÞ has the form
Sðe; ~e; gÞ ¼
ð1 þ gÞ2 : ~e þ egð2 þ gÞ
ð3Þ
For an nitrogen impurity (according to [4]) ~e eðKÞ 0:525J and g ¼ 0:5. As shown in [16], the Eq. (2) has a solution for both interval ½eðMÞ; eðKÞ and interval ½eðKÞ; eðMÞ. The local density of states of the nitrogen substitutional impurity calculated in [4] has quasi-local maxima in both intervals. For an boron substitutional impurity g ¼ 0:5 [4], quasi-localized states are absent in interval ½eðMÞ; eðMÞ [4,17]. Fig. 2 shows the graphical solution of the Lifshitz Eq. (2) for a given impurity atom. In this case also the Lifshitz equation has no solutions in interval ½eðMÞ; eðMÞ (corresponding to the dependence SðeÞ, curves 2 in Fig. 2). The local density of states of the boron impurity (curve 2 in Fig. 1) has two ðÞ poles outside the band of quasi-continuous spectrum ed , which are called local discrete levels and which are also solutions of Eq. (2). As is clearly seen in Fig. 1 the area under the curve 2 is smaller (by the sum of the residues at these poles) than the area under the curve 1. Local discrete levels can be an important source of information about defective structure and force interactions in real crystals. To extract this useful information we should have analytical expressions that relate the main characteristics of LDL (primarily their energy) to the parameters of the defect and the host lattice.
Fig. 2. Graphical solution of Eq. (2) for boron substitutional impurity in graphene. Curve 1 is the real part of the Green function of ideal graphene, curve 10 is the ‘‘approximations of two moments’’ of this function. Curve 2 represents function SðeÞ for boron impurity (~e ¼ 0:525J; g ¼ 0:5).
58
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
Such expressions were obtained in [17] for localized vibrations in the phonon spectrum of a threedimensional crystal. In [17] authors proposed the analytical approximation of the basic characteristics of local vibrations based on the rapid convergence of the real part of the Green’s function outside the band of quasi-continuous spectrum using the method of Jacobi matrices [13,14]. 3. Analytical approximation of the characteristics of local discrete levels Let us, briefly, to the extent necessary to understand the use of the classification of the eigenfunctions of Hamiltonian (1), we present the basics of the method of Jacobi matrices. This method allows, without finding the dispersion laws, to calculate directly the local partial Green’s functions of the system, corresponding to the perturbation of one or more atoms. This disturbance is described by a socalled generating vector ~ h0 2 H, where H is the space of electronic excitations of atoms. Its dimension is qN, where N is the number of atoms in the system, and q is the dimension of the displacement of a P single atom ðq ¼ 1; 2; 3Þ. If, using the generating vector ~ h0 ¼ p1 p jji (p is the number of excited j¼1
b n~ h0 g1 atoms) and the Hamiltonian (1) we construct the sequence f H n¼0 , then the linear envelope covering the vectors of this sequence forms, in the H space, a cyclic subspace invariant to the operator b This subspace contains, within itself, all the atomic displacements generated by the vector ~ H. h0 . The
corresponding partial Green’s function is determined as a matrix element 1 b n~ b h0 g1 e Ib H j0i. In the basis f H n¼0 which is obtained by orthonormalization of
G00 ðeÞ GðeÞ ¼ h0j
the sequence the operator (1) is represented in the form of a tridiagonal Jacobi matrix (or J-matrix). This matrix has a simple spectrum, what greatly simplifies finding the partial Green’s functions and spectral densities. As can be seen, this method does not use explicitly the translational symmetry of the crystal, what makes it extremely effective for treating systems in which such symmetry is broken. The method of Jacobi matrices is particularly effective for treating systems with a simply connected quasi-continuous band of spectrum D. In this case, with increasing rank of the J-matrix ðn ! 1Þ its diagonal elements an converge to a corresponding to the middle of the bandwidth D, and nondiagonal elements bn converge to b corresponding to one-quarter of the bandwidth D. For the local Green’s function corresponding to the excitations of one or more atoms, which are determined by the generating vector ~ h0 , we get following expression using the J-matrix method
Qn ðeÞ bn1 Qn1 ðeÞK1 ðeÞ ^ 1~ h0 Þ ¼ lim Gðe; ~ h0 Þ ¼ ð~ h0 ; ½eI^ H ; n!1 P n ðeÞ bn1 P n1 ðeÞK1 ðeÞ
ð4Þ
where polynomials P n ðeÞ are determined by the following recurrence relation
bn P nþ1 ðeÞ ¼ ðe an ÞP n ðeÞ bn P n1 ðeÞ:
ð5Þ
Polynomials Qn ðeÞ determined by the same recurrence relation. The initial conditions is 1 P 1 ðeÞ ¼ Q0 ðeÞ ¼ 0; P 0 ðeÞ ¼ 1; Q1 ðeÞ ¼ b0 . K1 ðeÞ is the function corresponds to the LGF operator, with all elements of its Jmatrix being equal to their limit values a and b
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 K1 ðeÞ ¼ 2b ½e a ZðeÞ ðe a þ 2bÞðe a 2bÞ;
ð6Þ
ZðeÞ 0ða 2b eÞ þ i Hðe a þ 2bÞHða þ 2b eÞ Hðe a 2bÞ:
ð7Þ
The method of Jacobi matrices can treat, as a regular singular perturbation, a much larger number of perturbations of the quasi-particle spectrum due to the presence of various crystal defects, than traditional methods [13,14]. In addition, perturbations do not change the bandwidth of the quasi-continuous spectrum, and consequently, the asymptotic values of the elements of the J-matrix can be regarded as an asymptotically degenerate regular perturbation [17]. This type of perturbations covers virtually all perturbations of the quasi-particle spectrum caused by local defects. The calculations of dynamics and thermodynamics characteristics of such systems are performed, using the method of J-matrices, with the same accuracy as for the initial ideal system.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
59
In practice, it is usually possible to calculate the Jacobi matrix of the Hamiltonian of a finite rank. The expression
Qn ðeÞ bn1 Qn1 ðeÞK1 ðeÞ Gðe; ~ h0 Þ ¼ P n ðeÞ bn1 P n1 ðeÞK1 ðeÞ
ð8Þ
is called analytical approximation of LGF. Quantity qðeÞ p1 ImG00 ðeÞ is called the spectral density generated by the initial displacement ~ h0 . All dependencies in Fig. 1 and curves 1 in Fig. 2 were calculated by the formula (8), using the Jacobi matrix of the Hamiltonian (1) with rank n ¼ 600. If we count the energy from the Fermi energy level, then all diagonal elements of Jacobi matrices are zero (8an ¼ a ¼ 0; b ¼ e0 =2). A good accuracy of the approximations shown in figures is confirmed also by the fact that they show the nonanalyticity effects corresponding to the densities of states of systems with the dimension higher than unity (so-called van Hove singularities). In the vicinity of these singularities the expression (8) slowly converges to the true values of the real and imaginary parts of LGF. Curve 10 in Fig. 2 show the real parts of the LGF calculated by formula (8) with n ¼ 1. As can be seen in the band of the quasi-continuous spectrum these relationships have very little in common with curve 1. On the curve 10 in the interval ½eðMÞ; eðMÞ are absent non-monotonous parts and the logarithmic singularities at the edges of the band of quasi-continuous spectrum, characteristic for the 2D systems. However, outside the band of quasi-continuous spectrum (also in the area of intersection of the real part of LGF with curve 2) curves 1 and 10 practically coincide. If in the Lifshitz equation (2) we put instead of LGF its approximation (8) for n ¼ 1, the obtained solutions give the energies of LDL with quite high accuracy. Moreover, these solutions can be easily found analytically. The LGF approximation by formula (8) for n ¼ 1 was named as the approximation of two moments [17]. Indeed, from the orthonormality of the polynomials [13,14], it follows that
Z
e0
P1 ðeÞ.ðeÞde ¼ 0 ) a0 ¼
e0
Z
e0
e0
P21 ðeÞ.ðeÞde ¼
Z
Z
e0
e.ðeÞde ¼ M1 ;
e0
e0
e2 a20 2
b0
e0
.ðeÞde ¼ 1 ) b0 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M 2 M21 :
Finding characteristics of LDL is more convenient without using Eq. (2), looking for them as the e e; h~0 Þ ¼ ðh~0 ; ½ebI H b K ^ 1 h~0 Þ. In the approximation of poles of the LGF perturbed Hamiltonian Gð two moments for the subspace generated by the excitation of an impurity atom, we get
e e; ~ Gð h0 Þ ¼
1
e a0 b20 K 1 ðeÞ
ð9Þ
:
For the case of an isolated substitutional impurity
a0 ¼ ~e;
b0 ¼
pffiffiffi 1þg 3ð1 þ gÞJ ¼ pffiffiffi e0 ; 3
ð10Þ
from where we get
e e; ~ h0 Þ ¼ Gð
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 cÞe a0 þ ZðeÞc je2 e20 j RðeÞ
;
ð11Þ
where RðeÞ ¼ ð1 2cÞe2 2ð1 cÞa0 e þ a20 þ c2 e20 and
c b20 =2b2 ¼ 2ð1 þ gÞ2 =3:
ð12Þ
Local discrete levels are poles (9), i.e. the roots of RðeÞ are ðÞ d
e
¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc 1Þa0 c a20 þ ð2c 1Þe20 2c 1
:
ð13Þ
60
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
ðÞ e eÞ are called the intensities of LDL and they determine the se¼eðÞ Gð Residues at these poles l0 ¼ re d relative LDL ‘‘amplitude’’ on the impurity atom. The condition that the intensity differs from zero defines the existence region of LDL. In this case
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ca0 þ ðc 1Þ a20 þ ð2c 1Þe20
a0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hð 1 þ cÞ; e0 2 2 ð2c 1Þ a0 þ ð2c 1Þe0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ca0 þ ðc 1Þ a20 þ ð2c 1Þe20 a0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Hðc 1 Þ: e0 2 2 ð2c 1Þ a0 þ ð2c 1Þe0
ðþÞ 0
l
¼
ðÞ 0
l
ð14Þ
In [17] it was shown that
b 1~ hn Þ ¼ P m ðeÞQn ðeÞ þ P m ðeÞP n ðeÞGðe; ~ Gmn ðe; ~ h0 Þ ¼ ðh~m ; ½ebI H h0 Þ:
ð15Þ
This implies that the damping of LDL, i.e. the decay of its intensity with increasing distance from the impurity atom (i.e. with the increase of n) follows the equation ðÞ 2 ðÞ 2 ðÞ e e lðÞ n ¼ rese¼eðÞ G nn ðeÞ ¼ P n ðed Þrese¼eðÞ GðeÞ ¼ P n ðed Þl0 : d
ð16Þ
d
Using the method of mathematical induction we can prove that
ðÞ d Þ
P n ðe
2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3n a20 þ ð2c 1Þe20 a0 pffiffiffiffiffiffi 5 : ¼ 2c4 e0 ð2c 1Þ
ð17Þ
and
l
ðÞ n>0
ðÞ n
¼ 2c½q
l
ðÞ 0 ;
qðÞ
2qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 32 a20 þ ð2c 1Þe20 a0 5 6 1: ¼4 e0 ð2c 1Þ
ðÞ
ð18Þ
Thus the intensities ln form with increasing n an infinitely decreasing geometric progression with ðþÞ ðÞ denominators qðÞ 6 1. Value qðÞ ¼ 1 if ed ¼ e0 or ed ¼ e0 . Summing these progressions, we see P1 ðÞ that l ¼ 1, that is the formation of each LDL is the formation of one quasiparticle outside n¼0 n the band of the quasi-continuous electron spectrum.
Fig. 3. Regions of the existence of discrete levels for substitutional impurities in graphene.
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
61
Fig. 4. The dependencies on g of basic characteristics of LDL on the substitutional impurity in graphene. The value ~e ¼ 0:525J corresponds to boron [4].
4. Analysis of the characteristics of the local discrete level on the boron subsitutional impurity Formulas (10)–(14) give simple analytical expressions for the conditions for the existence of local discrete levels due to the presence of a substitutional impurity in graphene. Regions of the LDL exisðþÞ ðÞ tence are for ed (in this case ~e > e0 ð1 cÞ ¼ e0 ½1 2ð1 þ gÞ2 =3) and for ed (in this case ~e < e0 ðc 1Þ ¼ e0 ½2ð1 þ gÞ2 =3 1) shown in Fig. 3. ~ and g. The absence of LDL is posIt is seen that such levels exist in a very widep range ffiffiffiffiffiffiffiffi of variables peffiffiffiffiffiffiffiffi sible only in a narrow range of values ~e at g < 3=2 1. If g > 3=2 1 at least one local discrete level exists. Realistically, the lines delineating the area of existence of LDL in graphene in the plane ð~e; gÞ must pass through the origin of coordinates, since ReGðeÞ ! 1 for e ! e0 . However, since this divergence is logarithmic, for any appreciable splitting of LDL from the boundary of quasi-continuous spectrum there is a certain threshold. Curves 1 and 10 in Fig. 2 merge at j e j> e0 . Fig. 4 shows, for the value ~e 0:525J (for boron, from [4]), the dependences of energies, the intensities of LDL at the boron impurity and the damping parameters of the value g that characterizes the change in the overlap integral of the boron impurity atom. The solid lines show the characteristics of LDL, calculated using the approximation of the Green’s function (9), i.e. according to the analytical formulas (13), (14) and (18). The open circles correspond to the results of the numerical calculations of ðÞ ðÞ dependences ed ðgÞ and l0 ðgÞ, using the Green’s function in the form of (8), calculated by the Jacobi matrices of the n ¼ 600 rank. ðÞ It is seen that at the threshold values ed ¼ e0 of the LDL formation, intensities vanish, and the ðÞ parameters of damping are equal to the unit. Further increase of j ed j is accompanied by the increase ðÞ of l0 and by the strengthening of the containment level.
5. Conclusions Good agreement between the results of numerical calculation of the LDL characteristics using Jacobi matrices of high rank, and their analytical description by the Green’s function (9) gives the possibility of extraction of the defect parameters from the known characteristics of LDL. Experimental measurement (e.g. by scanning tunneling microscopy) of values should allow, using (13) and (14), the determination of the parameters and this might lead to a significant advance in creating nanomaterials with predetermined spectral characteristics. As follows from (14) and can be seen from Fig. 4, ðþÞ ðÞ with increasing the intensity of LDL l0 ! l0 ! 1=2. That is, the impurity levels cannot be completely localized on the impurities, but they also appear in the spectra of carbon atoms surrounding an impurity atom. This greatly increases the probability of experimental detection of such levels, even
62
A. Feher et al. / Superlattices and Microstructures 53 (2013) 55–62
at low concentrations of impurities. Experimental study of the local discrete levels in graphene with boron impurities using STM techniques at temperatures near 0.3 K is in progress now. Acknowledgments This work is the result of the project implementation: Research and Education at UPJŠ – Heading towards Excellent European Universities, ITMS project code: 26110230056, supported by the Operational Program Education funded by the European Social Fund (ESF). References [1] A.M. Kossevich, The Crystal Lattice (Phonons Solitons Dislocations), Wiley VCH Verlag Berlin GmBH, Berlin, 1999. [2] K.S. Novoselov, A.K. Gein, S.V. Morozov, D. Diang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197–201. [3] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109–155. [4] N.M.R. Peres, F.D. Klironomos, S-W. Tsai, J.R. Santos, J.M.B. Lopes dos Santos, A.H. Castro Neto, EPL 80 (2007) 67007–67017. [5] J.C. Meyer, A.K. Geim, M.I. Katsnelson, K.S. Novoselov, T.J. Booth, S. Roth, Nature 446 (2007) 60–63. [6] L. Zhao, K.T. Rim, H. Zhou, R. He, T.F. Heinz, A. Pinczuk, G.W. Flynn, A.N. Pasupathy, Solid State Comm. 151 (2011) 509–513. [7] E. Stolyarova, K.T. Rim, S. Ryu, J. Maultzsch, P. Kim, L.E. Brus, T.F. Heinz, M.S. Hybertsen, G.W. Flynn, PNAS 104 (2007) 9209– 9212. [8] J. Xue, J. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, B.J. LeRoy, Phys. Rev. Lett. 108 (2012) 016801–1016801-5. [9] L. Zhao, R. He, K.T. Rim, T. Schiros, K.S. Kim, H. Zhou, C. Gutiérrez, S.P. Chockalingam, C.J. Arguello, L. Pálova, D. Nordlund, M.S. Hybertsen, D.R. Reichman, T.F. Heinz, P. Kim, A. Pinczuk, G.W. Flynn, A.N. Pasupathy, Science 333 (2011) 999–1003. [10] Yu.V. Skrypnyk, V.M. Loktev, Phys. Rev. B 73 (2006) 241402-1–241402-4. [11] Yu.V. Skrypnyk, V.M. Loktev, Low Temp. Phys. 34 (2008) 818–825. [12] C. Bena, S.A. Kivelson, Phys. Rev. B. 72 (2005) 125432-1–125432-7. [13] V.I. Peresada, Condens. Matter Phys. 4 (1968) 172–210. [14] V.I. Peresada, V.N. Afanas’ev, V.S. Borovikov, Sov. J. Low Temp. Phys. 1 (1975) 227–232. [15] I.M. Lifshits, Dokl. Akad. Nauk SSSR (in Russian) 48 (1945) 83–86. [16] A. Feher, E. Syrkin, S. Feodosyev, I. Gospodarev, K. Kravchenko, in: Physics and Applications of Graphene: Theory (ed. S. Mikhailov), InTech, Rijeka, 2011, pp. 93-112. [17] O.V. Kotlyar, S.B. Feodosyev, Low Temp. Phys. 32 (2006) 256–269.