Dynamical Polarizability of Graphene with Spatial Dispersion
DDynamical polarizability of graphene with spatial dispersion
Tao Zhu, ∗ Mauro Antezza,
2, 3 and Jian-Sheng Wang Department of Physics, National University of Singapore, Singapore 117551, Republic of Singapore Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Universit´e de Montpellier, F-34095 Montpellier, France Institut Universitaire de France, 1 rue Descartes, F-75231 Paris, France
We present a complete theory of electric polarizability of graphene using different theoreticalapproaches. From Kubo’s linear response formalism, we give a general expression of frequencyand wave-vector dependent polarizability within the random phase approximation. Four theoreticalapproaches have been applied to the single-layer graphene and their differences are given by theband-overlap of wavefunctions. By comparing with the ab initio calculation, we discuss the validityof methods used in literature. Our results show that the tight-binding method is as good as thetime-demanding ab initio approach in calculating the electric polarizability of graphene. Moreover,due to the special Dirac-cone band structure of graphene, the Dirac model reproduces results of thetight-binding method for energy smaller than 3 eV. For doped graphene, the intra-band transitionsdominate at low energies and can be described by the Lindhard formula for two-dimensional electrongases. At zero temperature and long-wavelength limit, all theoretical methods reduce to a long-waveanalytical formula and the intra-band contributions agree to the Drude polarizability of graphene.Effects of electrical doping and temperature are also discussed. This work may provide a solidreference for researches and applications of the screening effect of graphene.
I. INTRODUCTION
As the first isolated single-layer material, graphene hasattracted intense interest for its unique electronic struc-tures [1–3]. Due to the unusual gapless linear band struc-ture with a vanishing density of state at the Fermi en-ergy, screening properties in graphene differ significantlyfrom conventional two-dimensional (2D) materials andbecome one of the most important quantities for prospec-tive applications in many fundamental physics [4]. Thedynamical electric screening can be charactered by a fre-quency and wave-vector dependent polarizability Π( q , ω )which describes the response of induced charge density toa time-dependent effective potential. The quantities suchas electrical conductivity, dielectric function, electricalsusceptibility, reflection coefficient, have simple relationswith the electric polarizability that usually involve thepure Coulomb interaction v . For example, the electricaldielectric function has the form ε = 1 − v Π with its zeroscorrespond to the poles of the response function that de-scribe the plasmon modes. In the momentum space, the2D electrical conductivity has the relation σ = iω Π /q that can be obtained from the linear response electro-dynamics. These quantities have been proved to play apivotal role in developing a variety range of physical ap-plications especially for the near-field energy transport[5–12] and the Casimir effects [13–18].In general, the polarizability of graphene is nonlo-cal that depends on both the frequency and the wave-vector. Other factors like temperature and electricaldoping (chemical potential) can also play an importantrole. Many theoretical approaches have been developedto study the screening effect of graphene, for both po- ∗ [email protected] larizability [19, 20] and conductivity [21–26]. In par-ticular, Wunsch et al [19] gives an analytical formulaof dynamical polarizability with the Matsubara tech-nique at zero temperature, which agrees with the work ofHwang and Das Sarma [20]. In the long-wavelength limit q →
0, Falkovsky and Varlamov [23, 24] obtained sim-ple analytical expressions of conductivity of graphene viaGreen’s functions approach. Moreover, Klimchitskaya etal [25, 26] performed the quantum electrodynamics in(2+1)-dimensional space-time to get the components ofpolarization tensor and conductivity of graphene. How-ever, most of the theoretical treatments are based onthe Dirac model that the energy dispersion is assumedto be linear around the K points of the Brillouin zone.This assumption is only applicable to the nearest π -bandof graphene with energy less than a few electron volts.Moreover, with finite electrical doping, the Drude modelhas been widely used to calculate the intra-band contri-butions to the heat transfer and the Casimir force [5, 13].A detailed discussion on the validity of these theoreticaltreatments is still lacking.In this article, we present a comprehensive derivation,comparison, and discussion of different methods to cal-culate the frequency- and wave-vector-dependent electricpolarizability of graphene. We study four theoretical ap-proaches, say, density functional theory (DFT) based abinitio method, nearest-neighbor tight-binding method,the Dirac model, as well as the Lindhard formula fortwo-dimensional electron gas (2DEG) [27, 28]. Fromthe linear response Kubo formula [29], we get a generalexpression of the electrical polarizability and then em-ploy it to graphene with different approximations. Theonly difference between different methods is given bya so-called prefactor that describes the band-overlap ofwavefunctions. Our results show that the tight-bindingmethod can produce comparable results to the ab initio calculation. When the energy is smaller than 3 eV, the a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Dirac model reproduces the results of the tight-bindingmethod. With finite electrical doping, intra-band tran-sitions become dominant at low energy and can be welldescribed by the Lindhard formula.
II. POLARIZABILITY
To get explicit response function from the microscopicelectronic structure, we consider a system under smalltime-dependent external perturbations that the Hamil-tonian has the form H = H + H ( t ), where H is theunperturbed term and H represents time-dependent in-teractions between the charge density ρ ( r ) and an exter-nal electric potential φ ( r , t ). Assuming the perturbationis switched on adiabatically, we can write the interactionterm as H ( t ) = (cid:90) d r ˆ ρ ( r ) φ ( r , t ) e ηt , (1)where ˆ ρ is the charge density operator and η is an in-finitesimal positive damping factor which ensures thatthe perturbation vanishes at t → −∞ . In the interac-tion representation, the charge density operator ˆ ρ ( r , t ) = e iH t/ (cid:126) ˆ ρ ( r ) e − iH t/ (cid:126) and the expectation value of the in-duced charge density δρ is a linear response of the exter-nal perturbation: (cid:104) δ ˆ ρ ( r , t ) (cid:105) = (cid:90) d r (cid:48) χ ( r , r (cid:48) , t − t (cid:48) ) φ ( r (cid:48) , t ) , (2)where χ is the electric susceptibility that is given by theKubo formula [29] χ ( r , r (cid:48) , t − t (cid:48) ) = − i (cid:126) Θ( t − t (cid:48) ) (cid:104) [ˆ ρ ( r , t ) , ˆ ρ ( r (cid:48) , t (cid:48) )] (cid:105) . (3) where Θ is the Heaviside step function that equals 1 for t > t (cid:48) and 0 otherwise.Here we are more interested in the electric polarizabil-ity Π which represents the linear response of the inducedcharge density to the total (external + induced) poten-tial of the system. The relation between Π and χ is givenby the Dyson equation χ = Π + Π vχ where v is the bareCoulomb interactions of 2D system [30, 31]. For non-interacting electron systems, we can consider only thebare bubble diagram of the Dyson equation so that theelectric polarizability in the random phase approximation(RPA) [32] writesΠ( r , r (cid:48) , ω ) = 2 e (cid:88) i,j ( f j − f i ) ψ j ( r ) ψ ∗ j ( r (cid:48) ) ψ ∗ i ( r ) ψ i ( r (cid:48) ) (cid:15) j − (cid:15) i − (cid:126) ω − iη , (4)where we transformed Eq. (3) into frequency domain andcalculate the expectation value of density operator com-mutator with independent particle wave-functions ψ andcorresponding energies (cid:15) of unperturbed Hamiltonian H .The index i, j denote different states, the factor 2 ac-counts for the spin degeneracy, and e is the electroncharge. The function f = (1 + e β ( (cid:15) − µ ) ) − is the Fermidistribution function with β = 1 /k B T and µ the chemicalpotential.For systems with a periodic crystal, it is more conve-nient to calculate the polarizability in momentum space.According to the Bloch theorem, the sum over states inEq. (4) can be replaced by a sum over states and k-points( i → n k , j → n (cid:48) k (cid:48) ) where k lies within the first Brillouinzone. Due to the lattice translation symmetry, we haveΠ( r + R , r (cid:48) + R (cid:48) ) = Π( r , r (cid:48) ), where R is any real spacelattice vector. It can be shown from Fourier transform ofEq. (4) that Π( q , q (cid:48) , ω ) is only nonzero if q and q (cid:48) differby a reciprocal lattice vector G and the quasi-momentumconservation holds, i.e. k (cid:48) = k + q modulo G . Writ-ing wave-function ψ n, k ( r ) = u n, k ( r ) e i k · r with the Blochfunction u n, k ( r ), we haveΠ G , G (cid:48) ( q , ω ) = 2 e N Ω (cid:88) n,n (cid:48) , k (cid:104) u n k | e − i G · r | u n (cid:48) , k + q (cid:105)(cid:104) u n (cid:48) , k + q | e i G (cid:48) · r (cid:48) | u n k (cid:105) f n (cid:48) , k + q − f n, k (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη , (5)where N is the number of k-points in the first Brillouinzone and r is the electron position operator. For 2Dmaterials, Ω is the area of the primitive cell and q is thein-plane Bloch wave-vector. III. MODELS FOR GRAPHENE
The general formula of polarizability Π in Eq. (4) andEq. (5) are summing over transitions from different statesand the key quantity is the unperturbed stationary eigen-functions and corresponding eigenvalues. In this sec- tion, we apply this general formula to graphene with spe-cific band structures obtained from different theoreticalframeworks and approximations. A. ab initio We can calculate the electronic band structure ofgraphene from the ab initio methods that solve the many-body Schr¨odinger equation in the self-consistent-field ap-proach. In principle, many first principles methods, suchas the GW method [33] and the time-dependent densityfunctional theory [31], can be used to obtain the excita-tion properties. In this work, we calculate the polariz-ability of graphene using DFT in the Kohn-Sham scheme[34] that the effective Hamiltonian writes H = − (cid:126) m ∇ + ( V ext + V H + V xc ) , (6)where V ext is the external potential, V H the Hartree po-tential, and V xc the exchange-correlation potential. The Kohn-Sham equation can be solved self-consistently toget the ground state wavefunctions which have the formof Slater determinant [35] of single-particle orbitals. Inthe RPA, the independent particle polarizability yieldsthe Adler-Wiser formula [36, 37]Π G , G (cid:48) ( q , ω ) = 2 e N Ω (cid:88) n,n (cid:48) , k (cid:104) φ n k | e − i ( q + G ) · r | φ n (cid:48) k + q (cid:105)(cid:104) φ n (cid:48) k + q | e i ( q + G (cid:48) ) · r (cid:48) | φ n k (cid:105) f n (cid:48) k + q − f n k (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη , (7) FIG. 1. Honeycomb lattice of graphene and its Brillouin zone.(a) Lattice structure of graphene, a = 1 . (cid:6) A the distancebetween two nearest neighbor carbon atoms; a and a arelattice vectors of the unit cell which consists of two differentsites A and B; δ i , i = 1 , , b and b are reciprocalvectors. The Dirac cones locate at the K and K (cid:48) points. where φ n k and (cid:15) n k are Kohn-Sham wavefunctions andeigenvalues, respectively. This formula reduces to Eq. (5)if writing the Kohn-Sham wavefunctions in the form ofBloch functions, i.e. φ n, k ( r ) = u n, k ( r ) e i k · r . Becausethe ground state electronic structure corresponds to zerotemperature, the Fermi function becomes a step functionthat f equals 1 for occupied states and 0 for unoccupiedstates.It is worth noting that Eq. (7) is not limited to theKohn-Sham-RPA scheme and further improvements arepossible. For example, one can substitute the Kohn-Sham eigenvalues (cid:15) n k by the quasiparticle particle en-ergies E n, k obtained from GW method (GW-RPA) asthe former has intrinsic problems of underestimating theenergy states. Moreover, the electron-hole interactionscan be taken into account by solving the Bethe-Salpeterequation [38] on top of the GW-RPA calculations [31]. B. Tight-binding
Figure 1 shows the hexagonal lattice structure ofgraphene and its Brillouin zone. The unit cell of graphene consists of two different sites A and B . Considering onlynearest-neighbor hopping, the tight-binding Hamiltonianfor graphene has the form [4] H = − γ (cid:88) (cid:104) i,j (cid:105) ( a † i b j + b † j a i ) , (8)where a i ( a † i ) annihilates (creates) an electron on site r i on sublattice A (similarly for b j , b † j on sublattice B). γ ≈ . (cid:104) i, j (cid:105) meansthe nearest neighbors are summed once for each pair.Due to the honeycomb lattice of carbon atoms, thematrix representation of Hamiltonian of graphene is di-agonal in the momentum space H = (cid:20) − γf ( k ) − γf ∗ ( k ) 0 (cid:21) , (9)where f ( k ) = (cid:80) i e i k · δ i with δ = a ( − , T , δ = a (1 / , √ / T , and δ = a (1 / , −√ / T as depictedin Fig. 1. The eigenvalues of the tight-binding Hamilto-nian can be written as (cid:15) n, k = nγ | f ( k ) | where n = ± − e − iϕ ( k ) = f ( k ) / | f ( k ) | , the cell normalized Bloch eigenstates ofgraphene become S n, k = 1 √ (cid:20) ne iϕ ( k ) (cid:21) . Because the charge distribution is homogeneous for thetwo sub-lattice of graphene, one can neglect the recipro-cal lattice vector G (i.e., the local field effect [39, 40]) andfinally write the tight-binding polarizability of grapheneas [7]Π( q , ω ) = 2 e N Ω (cid:88) n,n (cid:48) = ± , k S † n, k S n (cid:48) , k + q S † n (cid:48) , k + q S n, k f n (cid:48) , k + q − f n, k (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη . (10) C. Dirac model
We can numerically calculate the tight-binding polar-izability of Eq. (10) with proper sampling of k and q inthe first Brillouin zone. However, due to the special elec-tronic band structure, the energy dispersion of grapheneis linear around the vicinity of K and K (cid:48) points of theBrillouin zone. When q is small, we can expand the f ( k ) around the Dirac point and the matrix element of Eq. (9) becomes f ( k ) = − a ( k x − ik y ) / H = (cid:126) v F k · σ where σ i are the Pauli matrices. TheFermi velocity v F = 3 a γ/ (cid:126) has the value ∼ m/sthat does not depend on the energy or momentum. Theeigenvalues of the Dirac Hamiltonian can be expressed as (cid:15) n k = n (cid:126) v F | k | . Substituting the eigenvectors in Eq. (10)into Eq. (10) and taking into account the valley degen-eracy, we can further write the Dirac polarizability ofgraphene asΠ( q , ω ) = 4 e (cid:88) n,n (cid:48) = ± (cid:90) d k (2 π )
12 (1 + nn (cid:48) cos θ ) f n (cid:48) , k + q − f n, k (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη , (11)where we replaced the sum of k points in the first Bril-louin zone by an integral of k in the momentum spacebecause the Dirac Brillouin zone 2 π/a → ∞ . The an-gle θ = ϕ ( k + q ) − ϕ ( k ) and the factor 4 accounts forboth spin and valley degeneracies. In the Dirac model,the phase factor has the form e iϕ ( k ) = ( k x + ik y ) / | k | and thus θ is just the angle between k and k + q in themomentum space.It is worthwhile to emphasize that Eq. (11) is a com-plete expression that depends on wave-vector, frequency,chemical potential, as well as temperature. This formulacan be evaluated numerically by integrating the Brillouinzone with writing cos θ = k · ( k + q ) | k || k + q | (12)for all non-vanishing wave-vector of k and k + q .Besides the general expression of Dirac polarizabilityof Eq. (11), many analytical formulas based on the Diracmodel have been derived under certain conditions andapproximations. Firstly, at zero temperature, the Fermifunction becomes a simple step function that the conduc- tion band is empty ( f +1 = 0) and the valence band is fullyoccupied ( f − = 1). For intrinsic undoped graphene, be-cause of the gapless Dirac band structure, only inter-bandtransitions are possible. However, real graphene samplesare usually doped that both intra- and inter-band tran-sitions can occur. The doping effect is characterized bya nonzero chemical potential µ in the Fermi function.At zero temperature, Wunsch et al [19] reported a rel-atively simple analytical formula of the polarizability ofgraphene with arbitrary wave-vectors. Similar work hasbeen done by Hwang and Das Sarma [20] and their for-mula produce exactly the same numerical results as thatof Wunsch’s. Secondly, in the optical region, the mag-nitude of the wave-vector q → et al [25, 26] report an analytical formula of conductivityof graphene which agrees to previous reports.For the nonlocal ( q →
0) polarizability of grapheneat zero temperature ( T →
0) with ω > v F q fixed, theconsensus was reached and all reported analytical formu-las reduce to what we call the long-wave polarizability ofgraphene [19, 20, 23–26]Π ( T,q ) → ( ω ) = q e π ( (cid:126) ω + iη ) (cid:20) µ (cid:126) ω + iη + 12 ln | µ − (cid:126) ω − iη || µ + (cid:126) ω + iη | − i π (cid:126) ω − µ ) (cid:21) . (13)The first term in the square brackets of Eq. (13) givesintra-band transitions that proportional to the chemicalpotential. For graphene, the chemical potential can be written as µ = (cid:126) v F k F with the Fermi momentum k F = (cid:112) π | n | and | n | the 2D carrier density [20]. Given therelation σ = iω Π /q for 2D systems, the first term ofEq. (13) yields σ ( ω ) = iD/π ( ω + i Γ) which is exactly theDrude conductivity of graphene [41]. The scattering rateΓ = η/ (cid:126) and the Drude weight D = ( v F e / (cid:126) ) (cid:112) π | n | .The second and third terms represent inter-bandtranstions which become negligible when (cid:126) ω << µ .Nevertheless, for (cid:126) ω > µ , the inter-band transitions con-tribute and the corresponding conductivity of graphenehas the form of σ = e / (cid:126) which is the so-called univer-sal conductivity of graphene [25, 42]. D. Lindhard
We have shown that at finite electrical doping, thetransport property of graphene is governed by the intra-band transitions when (cid:126) ω << µ . With large carrierconcentrations, one can infer that the intra-band contri-butions to the polarizability of graphene correspond tothat of 2DEG. In this case, the wavefunction has the formof plane-waves so that the cell normalized Bloch func-tion u n, k become a constant and the overlapping term ofwavefunctions can be omitted. Then we have the Lind-hard polarizability of 2DEG as [27]Π( q , ω ) = 4 e (cid:88) n,n (cid:48) (cid:90) d k (2 π ) f n (cid:48) , k + q − f n, k (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη . (14)With the Dirac dispersion (cid:15) n k = (cid:126) v F | k | , one can nu-merically solve Eq. (14) to get the intra-band polarizabil-ity of graphene. Interestingly, Eq. (14) corresponding tothe Dirac polarizability of Eq. (11) with θ = 0 fixed thatimplies k + q has the same direction as k for intra-bandtransitions.For small wave-vector q , we can expand the numera-tor of Eq. (14) as q ∇ k f k and the difference of energiesin the denominator as (cid:126) qv k with v k = v F k / | k | is thecarrier velocity. Then we get the analytical formula ofpolarizability of graphene given by Svintsov and Ryzhii[43]Π( q , ω ) = 2 e π (cid:126) v F (cid:18) ln(1 + e βµ ) β (cid:19) (cid:18) s √ s − − (cid:19) , (15)where β = 1 /k B T and s = ( (cid:126) ω + iη ) / (cid:126) v F q .In Eq. (15), when T →
0, the first bracket term reducesto µ and when q →
0, the second bracket term reducesto 1 / s . In this case, we may writeΠ ( T,q ) → ( ω ) = µq e π ( (cid:126) ω + iη ) , (16)which is exactly the Drude term of polarizability asshown in Eq. (13). E. Summary
We have derived expressions of the polarizability ofgraphene from different methods. In the linear response scheme, the general formula of polarizability is given byEq. (5) that involves the eigenvalues and eigenfunctionsof unperturbed Hamiltonian. Due to the special elec-tronic properties of graphene, different theoretical frame-works and approximations can be applied at certain con-ditions. Up to now, we have shown four theoretical ap-proaches, i.e., ab initio , tight-binding, Dirac, and theLindhard formula for intra-band transitions. For all ofthese methods, we may write a unified formula of thepolarizability of graphene as:Π( q , ω ) = 2 e (cid:88) n,n (cid:48) (cid:90) d k (2 π ) F n,n (cid:48) k , q ( f n (cid:48) , k + q − f n, k ) (cid:15) n (cid:48) k + q − (cid:15) n k − (cid:126) ω − iη , (17)where the prefactor F n,n (cid:48) k , q denotes the band-overlap ofthe wavefunctions in Eq. (5). Neglecting the local fieldeffect, we can summarize the prefactor of each methodas:General form : (cid:104) u n k | u n (cid:48) , k + q (cid:105)(cid:104) u n (cid:48) , k + q | u n k (cid:105) Ab initio : (cid:104) φ n k | e − i q · r | φ n (cid:48) , k + q (cid:105)(cid:104) φ n (cid:48) , k + q | e i q · r (cid:48) | φ n k (cid:105) Tight-binding : S † n, k S n (cid:48) , k + q S † n (cid:48) , k + q S n, k Dirac : g v nn (cid:48) cos θ )Lindhard : g v (18)where g v = 2 is the valley degeneracy for graphene. IV. COMPARISON AND DISCUSSION
In this section, we compare the polarizability ofgraphene calculated from different models as discussedabove. We treat the parameter-free ab initio results asthe reference to discuss the validity of each method withdifferent wave-vectors and frequencies. Dependences onchemical potential and temperature are also discussed.It has been reported that the temperature dependence ofelectron transport in graphene is remarkably small in thetemperature range 0.03-300 K [44]. So we fix T = 300 Kfor all calculations except for the discussion on tempera-ture dependence. The damping factor η is set to 0 .
05 eV.We numerically solve Eq. (10), Eq. (11), and Eq. (14) tocalculate the tight-binding, Dirac, and Lindhard polar-izability of graphene, respectively. The nearest neighborhopping parameter γ is set to 2 . v F = 9 . × m/s. Other analyticalexpressions from the literature are also calculated andcompared.For the ab initio calculation, we start from the groundstate calculations by using DFT as implemented inQUANTUM ESPRESSO [45, 46]. Then the ab initio polarizability of graphene is calculated on top of theground state band structure by using the package Berk-erelyGW [47, 48]. We adapt the norm-conserving pseu-dopotential generated by the Martins-Troullier approach FIG. 2. Comparison of real (Π ) and imaginary (Π ) polariz-ability of graphene in different wave-vectors. Results of Wun-sch and Falkovsky are calculated from the analytical formulagiven in Ref. [19] and Ref. [24], respectively. [49] with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional in the generalized gradient approx-imation (GGA) [50]. A plane-wave basis set with 60Ry energy cut-off is used to expand the Kohn-Shamwave functions. We sample the first Brillouin zone bya 300 × × a = a = √ a = 2 . (cid:6) A. Toavoid interactions from the neighboring lattice in the z direction, a large lattice constant of a = 18 (cid:6) A is set tothe z direction of the unit cell. A. Wave-vector
We first consider the spatial dispersion of polarizabilityof graphene. In Fig. 2, we show the calculated polariz-ability of graphene obtained from different methods withdifferent wave-vectors. The frequency is fixed at 0 . µ = 0 . q < m − , the devi- ation of results from different methods gradually vanish.This agree to the fact that graphene is a Dirac materialand all formulas obtained from the Dirac model reduceto the long-wave formula when v F q << ω . Moreover, be-cause the frequency (cid:126) ω is not significantly greater than2 µ , the intra-band transitions dominate at small wave-vectors so that the Drude polarizability agree well withothers. At large wave-vectors, models with only intra-band transitions, i.e., Drude, Svintsov, and Lindhard, arefailed to describe the spatial locality of transport proper-ties of graphene. In particular, the Drude polarizabilitygrows as a positive parabola which presents an unrealscreening effect. However, when q < × m − , theDrude polarizability agree well with that of the Lindhardmodel. This suggests that one can safely use the Drudemodel to calculate the intra-band contributions in near-field heat transfer when wave-vector is less than 300 timesof ω/c where c is the speed of light. Both Lindhard andSvintsov polarizability saturate when q > × m − that implies the intra-band transitions is negligible fortransition with very large wave-vectors.Results calculated from Falkovsky’s formula and thelong-wave formula show a very similar character and theystart to deviate from other models when q > × m − .This is because both of them are based on the long-waveapproximation that is not supposed to be valid whenwave-vectors are large. On the other hand, the tight-binding results agree well with the ab initio results exceptfor a global shift. This shift originates from the deviationof tight-binding band structure from the ab initio calcu-lations at large wave-vectors. Moreover, the Dirac modelproduces the same result as the tight-binding model when q < × m − . This implies that the linear disper-sion of the Dirac-cone band structure extends up to 3 eV.Besides, polarizability calculated from the Wunsch’s for-mula coincides exactly with those of the tight-binding(and Dirac) model when q < × m − .Secondly, unlike the real part of polarizability definesthe screening of an electromagnetic wave in a medium,the imaginary part is responsible for the absorption of ra-diation and is proportional to the damping of oscillations.In Fig. 2(b), similar to the real part, models based on thelong-wavelength approximation (Drude, Long-wave, andFalkovsky) are only valid when q < m − and theyare monotonically decreasing with the increase of wave-vectors. The tight-binding, Dirac, and Wunsch modelgive almost identical values at all q values and they arevery close to ab initio results. Strong damping peaks areshown at q ∼ × m − , which is exactly the resonanceof (cid:126) v F q to the frequency (cid:126) ω = 0 . q > m − that indicates a weak absorptionat large wave-vectors. FIG. 3. Comparison of real (Π ) and imaginary (Π ) polariz-ability of graphene in different frequencies with q = 9 . × m − . Results of Wunsch and Falkovsky are calculated fromthe analytical formula given in Ref. [19] and Ref. [24], re-spectively. B. Frequency
Next, we discuss the dynamical screening propertiesof graphene. We first consider the case when the long-wavelength approximation is valid. It is shown in Fig. 2that when q < m − , the differences between eachmethod is not significant at small chemical potential andlow temperature. In Fig. 3 we show the polarizabilitiesof graphene calculated from different methods with wave-vector fixed at q = 9 . × m − and chemical potential µ = 0 . ab initio calculations, which further supports its validity at lowfrequencies. All Dirac-type models (Dirac, Long-wave,Falkovsky, and Wunsch) produce almost identical resultsand consistent with that of the tight-binding model.In Fig. 4, we compare polarizabilities of graphene ob-tained from different methods with wave-vector fixed at q = 2 . × m − and chemical potential µ = 0 . ∼ .
75 eV which resonant to the valueof (cid:126) v F q . When frequencies are smaller than 3 eV, the re- FIG. 4. Comparison of real (Π ) and imaginary (Π ) polariz-ability of graphene in different frequencies with q = 2 . × m − . Results of Wunsch is calculated from the analytical for-mula given in Ref. [19]. sults of tight-binding, Dirac, as well as Wunsch’s modelare consistent with each other. However, when frequen-cies exceed 3 eV, two prominent absorption peaks areshown in the imaginary part of the ab initio polariz-ability at 4 . . π band in the Dirac model is not appropri-ate at high energies. However, except for a global shift,the tight-binding model successfully reproduces the maincharacter of ab initio polarizability at high frequencies,agree to previously discussed. C. Chemical Potential
Besides the wave-vector and frequency, effects fromdoping (chemical potential) and temperature can alsoplay an important role in screening effects of graphene.Figure 5 shows the calculated polarizability of graphenewith different chemical potentials. The wave-vector isfixed at a small value of 9 . × m − and the fre-quency is fixed at 0 . abinitio method. In order to tune the chemical potential in ab initio method, one need to introduce extra charges inthe ground state calculations and it is unstable with largedoping. So, in Fig. 5, we show the chemical dependence ofgraphene polarizability obtained from other models. Aswe can see, both real and imaginary parts of polarizabil-ity have a linear relation to the chemical potential when µ >> (cid:126) ω . Results from the Dirac model coincide exactlywith that of the tight-binding model and they are asymp-totically approaching results from the Lindhard formula. FIG. 5. Comparison of real (Π ) and imaginary (Π ) polariz-ability of graphene in different chemical potentials. This implies that the intra-band transitions dominate thescreening effect at large chemical potentials. This agreeswith the fact that rich electrical doping makes the elec-tronic structure of graphene close to the 2DEG.On the other hand, when µ < (cid:126) ω , the inter-band tran-sitions contribute and the value of real part polarizabilitygradually decreases and saturates to a negative value atzero chemical potentials. In this case, the intra-bandtransitions vanish and the screening effect is from theinter-band transitions so that the system behaves as adielectric material.Similar behaviors are shown for results from modelswith the long-wave approximation. As shown, the Drudepolarizability is also linearly proportional to the chemicalpotential but with a different slope from the Lindhard re-sults. Because only intra-band transitions are considered,both the Drude and Lindhard model give a vanishingvalue at zero chemical potential. Again, results from thelong-wave formula asymptotically converge to the Drudepolarizability when chemical potential is large. D. Temperature
At last, we discuss the effect of temperature. The ab initio method calculates the ground state electronicstructure at zero temperature. As reported [44], thetemperature dependence is not significant for graphenepolarizability when
T <
300 K. In Fig. 6, we show thecalculated polarizability of graphene with temperatureranges from 200 K to 5000 K. The chemical potentialwas set to 0 . q = 9 . × m − and frequency (cid:126) ω = 0 . FIG. 6. Comparison of real (Π ) and imaginary (Π ) polariz-ability of graphene in different temperatures. polarizability has a smaller slope than that of Dirac andtight-binding results which implies that temperature de-pendence of intra-band transitions is not as sensitive asthat of the inter-band transitions. The imaginary partpolarizability, however, is nonmonotonic for results fromthe Dirac and tight-binding model. It increasing withtemperature when T <
T <
V. CONCLUSION
In this work, we comprehensively discussed the dynam-ical polarizability of graphene with four theoretical ap-proaches. We derive the general expression of polariz-ability from the linear response Kubo formula. Differentmethods are distinguished by a prefactor that representsthe band-overlap of wavefunctions of unperturbed Hamil-tonian. We discussed the validity of each method by com-paring their results with different wave-vector, frequency,chemical potential, and temperature. At finite doping,the Lindhard formula for 2DEG describes the intra-bandtransitions of graphene and it reduces to the Drude for-mula when both wave-vector and temperature go to zero.For models with inter-band transitions, the tight-bindingmethod produces similar results as the ab initio calcula-tion except for a global shift that originates from thedeviation of their band structure at high energies. More-over, the Dirac model is equally good as the tight-bindingmethod when energy is smaller than 3 eV. At zero tem-perature, all theoretical models reduce to a simple an-alytical long-wave formula in the long-wavelength limitand its intra-band term corresponds to the Drude con-ductivity. The intra-band transitions becomes dominantat large electrical doping and all theoretical model giveasymptotic results with large chemical potential. Thetemperature dependence of the real part polarizability ofgraphene is approximately linear while a non-monotonictemperature dependence is shown for the imaginary part.Our results may provide a solid reference for applicationsof the response and screening properties of graphene andthe methods can be employed to other Dirac materials with appropriate theoretical treatments.
VI. ACKNOWLEDGMENTS
We wold like to thank Jiebin Peng, Zhibin Gao, andZu-Quan Zhang for insightful discussions. This work issupported by a FRC grant R-144-000-427-114 and anMOE tier 2 grant R-144-000-411-112. [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Gregorieva, and A. A.Firsov, Science , 666 (2004).[2] A. K. Geim, and K. S. Novoselov, Nat. Mater. , 183-191(2007).[3] A. K. Geim, Science , 1530 (2009).[4] A.H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[5] A. I. Volokitin and B. N. J. Persson, Rev. Mod. Phys. , 1291 (2007).[6] J. -S. Wang and J. Peng, EPL , 155437(2017).[8] Z.-Q. Zhang, J. T. L¨u, and J.-S. Wang, Phys. Rev. B ,195450 (2018).[9] R. Yu, A. Manjavacas, and F. J. Garc´ıa de Abajo, Nat.Commun. , 2 (2017).[10] T. Zhu, Z. -Q. Zhang, Z. Gao, and J. -S. Wang, Phys.Rev. Appl. , 245437 (2017).[12] M. -J. He, H. Qi, Y. -T. Ren, Y. -J. Zhao, and M. An-tezza, Appl. Phys. Lett. , 263101 (2019).[13] G. L. Klimchitskaya, U. Mohideen, V. M. Mostepanenko,Rev. Mod. Phys. , 1827 (2009).[14] D. Drosdoff and Lilia M. Woods, Phys. Rev. B , 155459(2010).[15] V. Yannopapas and N. V. Vitanov, Phys. Rev. Lett. ,120401 (2009).[16] W. -K. Tse and A. H. MacDonald, Phys. Rev. Lett. ,236806 (2012).[17] C. Abbas, B. Guizal, and M. Antezza, Phys. Rev. Lett. , 126101 (2017).[18] B. Guizal and M. Antezza, Phys. Rev. B , 115427(2016).[19] B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J.Phys. , 318 (2006).[20] E. H. Hwang and S. Das Sarma, Phys. Rev. B , 205418(2007).[21] V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, Phys.Rev. Lett. , 256802 (2006).[22] V. P. Gusynin, S. G. Sharapov and J. P. Carbotte, Int.J. Mod. Phys. B , 4611 (2007).[23] L. A. Falkovsky and A. A. Varlamov, Eur. Phys. J. B , 195405 (2016). [26] G. L. Klimchitskaya, V. M. Mostepanenko, and V. M.Petrov, Phys. Rev. B , 235432 (2017).[27] G. D. Mahan, Many-Particle Physics , (Springer, NewYork, 2000).[28] J. Lindhard, K. Dan. Vidensk. Selsk. Mat. Fys. Medd. , 570 (1957),R. Kubo, M. Yokota, and S. Nakajima, J. Phys. Soc. Jpn. , 1203 (1957).[30] F. J. Dyson, Phys. Rev. , 1736 (1949).[31] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601-659 (2002).[32] H. Ehrenreich, M. H. Cohen, Phys. Rev. , 786 (1959).[33] F. Aryasetiawan, and O. Gunnarsson, Rep. Prog. Phys.
237 (1998).[34] W. Kohn, and L. Sham. Phys. Rev. , A1133 (1965).[35] J. C. Slater, Phys. Rev. , 413 (1962).[37] N. Wiser, Phys. Rev. , 62 (1963).[38] E. E. Salpeter and H. A. Bethe, Phys. Rev. , 1232(1951).[39] T. Zhu, P. E. Trevisanutto, T. C. Asmara, L. Xu, Y. P.Feng, and A. Rusydi, Phys. Rev. B , 235115 (2018).[40] S. G. Louie, J. R. Chelikowsky, and M. L. Cohen, Phys.Rev. Lett. , 155 (1975).[41] J. Horng, C. -F. Chen, B. Geng, C. Girit, Y. Zhang, Z.Hao, H. A. Bechtel, M. Martin, A. Zettl, M. F. Crommie,Y. Ron Shen, and F. Wang, Phys. Rev. B , 165113(2011).[42] M. Lewkowicz and B. Rosenstein, Phys. Rev. Lett. ,106802 (2009)[43] D. Svintsov and V. Ryzhii, Phys. Rev. Lett. , 219401,(2019).[44] Y.-W. Tan, Y. Zhang, H. L. Stormer, and P. Kim, Eur.Phys. J. Special Topics , 15–18 (2007)[45] P. Giannozzi et al. , J.Phys.:Condens.Matter , 395502(2009).[46] P. Giannozzi et al. , J.Phys.:Condens.Matter , 465901(2017).[47] M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390(1986).[48] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M.L. Cohen, and S. G. Louie, Comput. Phys. Commun. , 1269 (2012)[49] N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991).[50] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett., , 3865 (1996). [51] H. J. Monkhorst, and J. D. Pack, Phys. Rev. B , 5188(1976).[52] E. H. Hwang and S. Das Sarma, Phys. Rev. B79