Dynamical polarization of graphene at finite doping

The polarization of graphene is calculated exactly within the random phase approximation for arbitrary frequency, wave vector, and doping. At finite doping, the static susceptibility saturates to a constant value for low momenta. At $q=2 k_{F}$ it has a discontinuity only in the second derivative. In the presence of a charged impurity this results in Friedel oscillations which decay with the same power law as the Thomas Fermi contribution, the latter being always dominant. The spin density oscillations in the presence of a magnetic impurity are also calculated. The dynamical polarization for low $q$ and arbitrary $\omega $ is employed to calculate the dispersion relation and the decay rate of plasmons and acoustic phonons as a function of doping. The low screening of graphene, combined with the absence of a gap, leads to a significant stiffening of the longitudinal acoustic lattice vibrations.


I. INTRODUCTION
Recent progress in the isolation of single graphene layers has permitted the realization of transport and Raman experiments 1 which have stimulated an intense theoretical research on the properties of a monoatomic graphene sheet. Most of the unusual electronic properties can be understood in terms of a simple tight-binding approach for the π-electrons of carbon, which yields a gapless linear band structure with a vanishing density of states at zero doping. 2,3 The electronic band structure of an undoped single graphene sheet allows for a description of the electronic properties in terms of an effective field theory which is equivalent to Quantum Electrodynamics (QED) in 2+1 dimensions. Within such a framework important physical quantities such as the electron self-energy or the charge and spin susceptibilities can be calculated exploiting the special symmetries of the model. 4 Finite doping away from half-filling qualitatively changes the above description, since it breaks electron-hole symmetry and also the pseudo Lorentz invariance needed for the equivalence with QED in two space dimensions. One is thus forced to retreat to conventional condensed-matter techniques such as e.g. Matsubara Green's functions. This was recently done for the static susceptibility of graphene at finite doping 5 . Similar calculations for bulk graphite had been performed by Shung et al., 6 who investigated the dielectric function and the plasmon behavior. Studies of static screening in graphene based on the Thomas-Fermi approximation had also been considered. 7,8 The relation between the polarization and the transport properties of graphene has been recently discussed in Refs. 5,9,10. Thermo-plasma polaritons in graphene are discussed in citeVafek.
In this paper, we calculate the dynamical polarization within the random phase approximation (RPA) for arbitrary wavevector, frequency, and doping. We discuss the method of calculation in the next section. We present in Section III results for static screening, where we analyze the Friedel oscillations induced by a charged or magnetic impurity, comparing our results with those of the two-dimensional electron gas (2DEG). Section IV discusses the plasmon dispersion relation and lifetime. In Section V we analyze the screening of the longitudinal acoustic modes by the conduction electrons. Finally, Section VI presents the main conclusions of our work. Some of the more technical aspects of the formalism are explained in the Appendix.

II. RPA CALCULATION
Within the effective mass approximation, and focusing on one of the two unequal Kpoints, the Hamiltonian of an hexagonal graphene sheet is given in the Bloch spinor repre- Here φ k = k x + ik y and ψ k = (a k , b k ) T , where a k and b k are the destruction operators of the Bloch states of the two triangular sublattices. We have also introduced the Fermi wavevector k F which is related to the chemical potential µ via k F = µ/hv F . The Fermi velocity v F = 3at/2h is determined by the carbon-carbon distance a = 1.42Å and the nearest neighbor hopping energy t = 2.7 eV resulting in v F ≃ 9 × 10 5 m/s. 11 We note that the effective Hamiltonian given above is valid only for wave vectors k ≪ Λ, where Λ ≃ 8.25 eV is a high-energy cutoff stemming from the discreteness of the lattice. 11 The quantity of interest for many physical properties is the dynamical polarization, since it determines e.g. the effective electron-electron interaction, the Friedel oscillations and the plasmon and phonon spectra. In terms of the bosonic Matsubara frequencies ω n = 2πn/β, it is defined as 12 where A denotes the area. The average is taken over the canonical ensemble and the density operator is given by the sum of the density operators of the two sub-lattices ρ = ρ a + ρ b .
This amounts to working in the long-wavelength limit. To first order in the electron-electron interaction, we obtain with E ± (k) = ±hv F k − µ the eigenenergies, n F (E) = (e βE + 1) −1 the Fermi function, and g S = g V = 2 the spin and valley degeneracy. A characteristic difference between the polarization of graphene and that of a 2DEG is the appearance of the prefactors f ss ′ (k, q) coming from the band-overlap of the wave functions 5,6 f ss ′ (k, q) = 1 2 1 + ss ′ k + q cos ϕ |k + q| , where ϕ denotes the angle between k and q.
At zero temperature, the Fermi functions yield simple step functions. We define the following retarded function by replacing iω n → ω + iδ where g ≡ g S g V . The +(−) sign corresponds to intra(inter)-band transitions and D is a general upper limit.
For µ = 0, the retarded polarization thus reads For µ > 0, i.e., for nonzero electron doping, the retarded polarization has an additional term In the Appendix we give additional details of the calculation of Eq. (5). Here we summarize these expressions in terms of two complex functions, F (q, ω) and G(x), defined as From now on it is always assumed that ω > 0, noting that the polarization for ω < 0 is obtained via P (1) (q, −ω) = P (1) (q, ω) * . Equations (6) and (7) are then rewritten in the following compact form with P (1) and Equations ( Two limits of the polarization are of particular importance: (i) The long wavelength limit q → 0 with ω > v F q fixed, which is relevant for optical spectroscopy and for the plasma dispersion. (ii) The static case ω = 0 with q arbitrary, which is relevant for the screening of charged or magnetic impurities.
For the first case we obtain And for the second case we recover previous results 5 where G < (x) ≡ −iG(x) is the real function (for |x| < 1) given in Eq. (A2). Note that for 2µ ≤hv F q ≤ Λ the absolute value of the polarization is linear in q, and acquires rapidly the behavior of P (1) 0 . By contrast, the polarization of the ordinary 2DEG decreases with increasing wave vector for q > 2k F . Furthermore, and also in contrast to the 2DEG, where the first derivative of the static polarization is discontinuous at q = 2k F , doped graphene has a continuous first derivative athv F q = 2µ and a discontinuous second derivative. Table I summarizes the ω dependence of ImP (1) for ω → 0 for a 2DEG and for doped graphene. While the dependence is the same for q = 2k F , a subtle difference appears at q = 2k F . This difference is intimately related to the nonanalytic behavior (as a function of q) of the static polarization at q = 2k F . It leads to e.g. a different power law decay of the electron screening, as will be shown in subsection III. Table I also indicates that, depending on the order in which the limits q, µ → 0 are taken, the low-frequency behavior may be that of an insulator, a metal, or a hybrid between the two.
The polarization is a continuous function of q and ω except for the square-root divergence which F (q, ω) shows at ω = v F q. Figure 1 shows real and imaginary parts of the polarization given by Eq. (9). We note that ImP (1) (q, ω) = 0 forhv F q <hω < 2µ −hv F q or 0 <hω < hv F q − 2µ, and negative otherwise, while ReP (1) (q, ω) < 0 for ω < v F q.
The divergence at ω = v F q vanishes in the self-consistent RPA result of the polarization given by 13 Here v q = e 2 /2κ 0 q denotes the in-plane Coulomb potential in vacuum. We note that so that the self-consistent polarization has a real and finite value at   Figure 2 shows real and imaginary parts of the self-consistent polarization. While the singularity at ω = v F q is absent, a new singularity appears in ReP RPA (q, ω) at ω ∝ √ q, which reflects the existence of plasmons, as will be discussed in subsection IV.

III. STATIC SCREENING
An external charge density n ext (r) = Zeδ(r) is screened by free electrons due to the Coulomb interaction. This results in the induced charge density δn(r) Here ǫ(q, 0) ≡ lim ω→0 ǫ(q, ω). Within the RPA approximation 6 The effective dielectric constant ǫ 0 includes high energy screening processes. We take ǫ 0 ≃ 2.4. 5 There are two contributions to the induced charge density. A non-oscillating part comes from the long-wavelength behavior of the polarization. This contribution is obtained within the Thomas-Fermi (TF) approximation and is given by where α ≡ e 2 /4πκ 0h v F ≃ 2.5. We note here that in a 2DEG the TF contribution also decays as r −3 . 15 The second contribution to the long distance behavior is oscillatory and comes from the non-analyticity of the polarization athv F q = 2µ. 14,15 However, and quite importantly, in graphene the non-analyticity results from a discontinuity occurring only in second derivative, the first derivative being continuous. This leads to an oscillatory decay which contrasts with the behavior of a 2DEG, where Friedel oscillations scale like δn(r) ∝ cos(2k F r) r −2 . This difference has been previously noted in Refs. 8,16, where however TF screening was not considered. We emphasize here that, because the TF contribution is of the same order of magnitude as the oscillatory part, it is essential to consider both of them. Furthermore we note that, while the TF contribution to screening given in Eq. (17) is independent of the dielectric constant ǫ 0 , the amplitude of the oscillatory part given in  remarkable general property can be clearly appreciated in Fig. 3 for the particular case of ǫ 0 = 2.4: The induced density δn(r) does not change sign, but oscillates around a finite offset.
The polarization also determines the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction energy between two magnetic impurities as well as the induced spin density due to a magnetic impurity, both quantities being proportional to the Fourier transform of P (1) (q, 0). 17 For a magnetic impurity at r the induced spin density δm(r) is Note that the magnetic response is proportional to the bare polarization, since the spin-spin interaction is mediated via the short ranged exchange interaction (in contrast to the long ranged Coulomb interaction in the case of charge response). That is why the TF contribution is missing and the induced spin density oscillates around zero. Figure 4 shows the numerically calculated P (1) (r) in units of k 3 F /hv F and illustrates this behavior. Specifically in the long wavelength limit we obtain which is clearly seen in the inset of graphene we recover the monotonous r −3 decay obtained in Ref. 18.
Finally we note, that in contrast to the behavior of charge screening, where the induced charge scales like 1/µ at large distances, the envelop function of the spin density is independent of the chemical potential at large distances.

IV. PLASMONS
The plasmon dispersion is determined by solving for ǫ(q, ω p − iγ) = 0, where γ is the decay rate of the plasmons. For weak damping, the plasmon dispersion ω p (q) and the decay rate γ are determined by 13 Solutions to the first equation only exist for ReP (1) > 0, which is the case only for finite doping and ω > v F q. Furthermore, a stable solution requires ImP (1) = 0, as is the case for region 1B of Fig. 6. Using the low-q expansion of the polarization given in Eq. (12), and neglecting the logarithmic correction, we obtain with α defined below Eq. (17). The √ q behavior of the plasmon dispersion also appears in the 2DEG.
Here P el , ǫ el are the polarization and the dynamical dielectric function of the electrons as given in Eqs. (9) and (16), while P ion , ǫ ion are the corresponding quantities for the ions. In the calculation of P ion (q, ω) we assume that the ions have a quadratic energy dispersion where M denotes the ion mass. The ionic charge density is two positive charges per unit cell, so that the Fermi wave vector of the ions is k ′ F ≃ (8π/A c ) 1/2 , where A c = 3 √ 3a 2 /2 denotes the area of the hexagonal unit cell in real space. 11 This value of k ′ F is exact for field effect doping and approximate for chemical doping.
In the calculation of P ion we assume that, for all relevant frequencies, we can take ω ≫ hk ′ F q/M. In the case of acoustic phonons this assumption is equivalent to v s ≫hk ′ where v s denotes the sound velocity. This relation is fulfilled for all meaningful dopings as shown below. In this regime the ion polarization is real and given by where we defined E 0 =h 2 /MA c ≃ 7 × 10 −5 eV. We note that E 0 is of the order of the ion confinement energy.
We may estimate the dispersion and the decay rate of the acoustical phonons by inserting Eq. (23) into Eq. (21). For finite doping µ > 0, the acoustic phonons at long wavelengths lie in the region 1A of Fig. 6 In this regime the phonon dispersion is easily obtained: To be consistent with the precondition ω < v F q the expression in the square root has to be smaller than one. This sets a lower limit to the values of the chemical potential for which Eq. (25) is valid. The sound velocity v s and the decay rate γ may be derived from Eqs. (21) and (25) in the limit of low q. We obtain with ξ = 4πE 0 /gµ. Some additional remarks on the validity of Eq. (25) go in place. We have already said that, for q → 0, it only applies provided ξ < 1. We also wish to note that the acoustical phonons are only well-defined if their frequency is much larger than their decay rate, i.e. if ω ph /γ ≫ 1. Since ω ph /γ = 1 for ξ = 4/5, we conclude that the notion of acoustic phonons is justified for ξ ≪ 1, which corresponds to µ ≫ 4πE 0 /g ≃ 2.2×10 −4 eV. A detailed discussion of acoustical phonons for ξ > 1 is left for future work.
We finish this section by estimating the sound velocity of a typical graphene sample.
Assuming a concentration of "conduction-band" electrons of n el = 10 10 − 10 12 cm −2 we get µ ≃ 10 −2 − 10 −1 eV. We note that this corresponds to ξ ≃ 2 − 20 × 10 −3 , so that On the other hand, in a semiconductor with a gap larger than the typical acoustic phonon frequencies, the electrons follow adiabatically the ions. In that case, the lattice vibrations can be described as oscillations of neutral particles.

VI. CONCLUSIONS
In this article we have derived a compact and closed expression for the dynamical polarization of graphene within the RPA approximation. The obtained result is valid for arbitrary wave vector, frequency, and doping. As particular cases, we have derived the longwavelength limit q → 0 and the static limit ω → 0. We have employed the RPA polarization The dynamical polarization has been used to calculate the dispersion relation and the decay rate of plasmons and acoustic phonons. Like in the 2DEG case, the plasmon frequency shows a √ q-behavior in the long wavelength regime. We have determined the region in the (q, ω) plane where the plasmon is stable, as well as the decay rate in the regime where it is not.
The dispersion of acoustical phonons has been shown to be strongly dependent on the chemical potential. In particular, we have found that the sound velocity approaches the Fermi velocity at low doping. However, the same limit shows an increase in the decay rate of acoustic phonons due to electron-hole pair excitation.
Although we have focused on applications of the RPA calculation to the case of doped graphene, some aspects of the low-frequency, long-wavelength dynamics of pure graphene appear to be intriguing and worth studying further.
Note added. When this work was about to be submitted, we became aware of a related paper. 21 Overlapping results in the two papers are in agreement.

Imaginary part
The imaginary part of the functions χ ± D (q, ω) defined in Eq. (5) has the following form dk α=± α I αβ (k, q, ω) , The ϕ-integration yields which is always real. The final k-integration can now simply be performed. We obtain for In order to present the result for µ > 0, we introduce the real functions f (q, ω), For the additional term at finite doping given by Eq. (7), we obtain in the language of Fig. 6 Im∆P

Real part
The Kramers-Kronig relation valid for the retarded function ReP  This integral is calculated directly such that the Kramers-Kronig relation is not needed, here. The ϕ-integration yields where J α (k, q, ω) is given by We thus get in the language of Fig. 6 Re∆P (1) (q, ω) = − gµ 2π