{'fraction_non_alphanumeric': 0.03189344543582704, 'fraction_numerical': 0.011968085106382979, 'mean_word_length': 4.687362771407661, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Plasmon anomaly in the dynamical optical conductivity of graphene

We theoretically consider the effect of plasmon collective modes on the frequency-dependent conductivity of graphene in the presence of the random static potential of charged impurities. We develop an equation of motion approach suitable for the relativistic Dirac electrons in graphene that allows analytical high-frequency asymptotic solution in the presence of both disorder and interaction. We show that the presence of the acoustic plasmon pole (i.e. the plasmon frequency vanishing at long wavelengths as the square-root of wavevector) in the inverse dynamical dielectric function of graphene gives rise to a strong variation with frequency of the screening effect of the relativistic electron gas in graphene on the potential of charged impurities. The resulting frequency-dependent impurity scattering rate gives rise to a broad peak in the frequency-dependent graphene optical conductivity with the amplitude and the position of the peak being sensitive to the detailed characteristics of disorder and interaction in the system. This sample (i.e. disorder, elecron density and interaction strength) dependent redistribution of the spectral weight in the frequency-dependent graphene conductivity may have already been experimentally observed in optical measurements.


I. INTRODUCTION
Electronic transport and electrical conductivity of graphene is a subject 1 of great current interest in both fundamental physics and in applied physics and engineering. Recently, optical conductivity measurements have been used to identify the effects of electron-electron interactions of relativistic Dirac electrons 2,3 which have attracted a great deal of interest. In particular, the linear band dispersion as well as its gapless chiral 2D nature make graphene a unique laboratory system for studying electron-electron interaction effects and their interplay with disorder since the lack of Galilean invariance associated with the non-classical linear band dispersion leads to qualitative and quantitative phenomena not present in regular parabolic band metals and semiconductors. For example, graphene violates the wellknown Kohn's theorem, and both cyclotron resonance and optical conductivity are directly affected by electronelectron interaction effects in contrast to the corresponding parabolic band systems. In the context of engineering applications, strong broadband absorption and high tunability make monolayer graphene uniquely suited for a wide range of optoelectronic and plasmonic devices 4 . These prospects have motivated numerous measurements of the optical conductivity of monolayer graphene in the wide frequency range from the far infrared to the ultraviolet [5][6][7][8][9][10][11][12] . A detailed understanding of the optical properties of monolayer graphene is required for the progress in both of these directions. In this paper we consider only monolayer graphene within the effective linear gapless Dirac band model, and 'graphene' in this article refers exclusively to 'monolayer graphene'. Our goal here is to theoretically study the effects of electronic collective modes on the graphene dynamical conductivity treating disorder and interaction on the same footing.

A. Experimental motivation
The optical conductivity of undoped graphene is dominated by interband processes (arising from the particlehole excitations from the filled valence band to the empty conduction band) in the gapless Dirac spectrum. In the non-interacting case this gives rise to a universal frequency-independent value of the conductivity that can be expressed in terms of fundamental constants, a prediction 13,14 remarkably confirmed by the measurements 11,15 , establishing the basic gapless chiral massless Dirac band dispersion of graphene to be valid. At finite doping the interband processes are strongly suppressed by Pauli blocking in the partially filled conduction band (assuming positive Fermi level E F > 0, but obviously the same physics also applies for the hole-doped system) whereas the universal value of the conductivity is reached at higher frequencies ω 2E F where Pauli blocking effects become small, and the interband transitions again dominate the optical conductivity. The partially filled conduction band gives rise to the intraband response reflected in the Drude peak at ω ≈ 0 in the optical conductivity. In graphene electron-electron interactions are expected to give rise to relatively small yet observable corrections to the optical conductivity via the many-body renormalization of the electron spectrum and the excitonic effects near the absorption threshold ω = 2E F 1,3 . An important difference between graphene and ordinary parabolic band semiconductor systems in this context is that the linearly dispersing graphene electron and hole bands violate the usual Galilean invariance dominating the long wavelength optical response of the semiconductor systems where the electron-electron interaction, being a property of the relative coordinates of the individual electrons themselves, does not affect the long wavelength translationally invariant optical conductivity of the system since the latter is a specific property of the arXiv:1305.4940v1 [cond-mat.mes-hall] 21 May 2013 center of mass coordinates coupling to the external light in the long wavelength limit. The long wavelength optical properties in graphene, however, are explicitly affected by the interaction effects since graphene obeys Lorentz invariance which does not enable the usual separation of the system Hamiltonian into the center of mass and relative coordinates components.
In particular, the Drude weight in graphene has been predicted to be renormalized by the interaction effects 16 . Experimentally this effect is demonstrated by analyzing the parameters of the Drude peak. The latter is routinely phenomenologically fitted by a 2-parameter Lorentzian σ = iD/(π(ω + i/τ )) characterized by a width 1/τ due to the disorder scattering and the Drude weight D 2,3 . The latter has shown discrepancy with the free-electron value D = (ve 2 / ) √ πρ suggestive of electron-electron interactions playing a role.
Here v is the slope of the linear dispersion in graphene (i.e. the so called graphene velocity) and ρ is the electron density. On the other hand, the data reported in Ref. 17 suggest a non-Lorentzian shape of the Drude peak in graphene possibly reflecting a frequency-dependent relaxation rate 1/τ (ω) or the presence of other non-universal background at ωτ 1. More recently, additional measurements of the plasmon decay 18,19 , directly related to the optical conductivity, confirmed the discrepancies with the simple Lorentzian fits. The questions of the origin of the non-Lorentzian shape of the Drude peak and the actual value of the Drude weight extracted from the Drude response remain open.

B. Theoretical motivation
Several theoretical models of the optical response of graphene have been put forward that include the effects of disorder 13,20,21 , the band structure corrections due to the next-nearest neighbor hopping in the tight binding description 22 , excitonic effects near ω ≈ 2E F 23,24 , effects of phonons 25-30 relevant at high temperatures, plasmaron effects 31 near the Dirac point, the effect of strain 32 and the spectrum renormalization effects of electronelectron interactions 29,33 .
An effect omitted so far in the literature is the frequency dependence of the Lorentz invariant screening effect of the graphene relativistic electron gas. Formally, the dynamical screening of an electron gas is characterized by the inverse of the dynamical dielectric function, 1/ε(q, ω). The presence of a plasmon pole in the dynamical dielectric function, m[1/ε(q, ω)] ∼ δ(ω − ω pl (q)), results in the failure of the screening at the plasmon frequency ω ≈ ω pl (q) and therefore gives rise to an enhancement of the disorder scattering rate. As a result, the acoustic plasmon dispersion ω pl (q) (i.e. the plasma frequency going as the square-root of the wavevector at long wavelength) in the 2D electron gas (2DEG) gives rise to a strong frequency dependence of the scattering rate due to charged impurities [34][35][36][37] .
Despite the conceptual simplicity, a formal calculation of the frequency-dependent scattering rate in the electron gas in the presence of disorder and interaction is not straightforward. The standard method is to construct a diagrammatic perturbation theory for the current-current correlator relying on the "metallic" expansion parameter E F τ characterizing the disorder strength and then employing some sort of approximation scheme for the electron-electron interactions. Such a construction however requires the inclusion of an infinite series of diagrams. The summation of the infinite series results in an integral Bethe-Salpeter equation for the vertex function which is difficult and cumbersome to solve in the presence of both disorder and interaction even for the simplest approximation schemes.
In the case of the 2DEG with the parabolic dispersion the problem was successfully solved with a lot of work using the so-called 'memory function' approach 34 . This solution relied on a substantial simplification that is possible in the parabolic case due to the separation between the slow dynamics of the center of mass of the electron gas and the relative motion of electrons. As a result the dynamics of the 2DEG is described by an effective quantum Langevin equation for a macroscopic particle with a large mass (the mass of the electron gas) subject to a random force due to the disorder potential 38 . In this case the force-force correlator determines the response of the 2DEG to the external field. Electron-electron interactions enter the force-force correlator in a particularly simple way, only via the polarization operator, which greatly simplifies the calculation 34 .
In graphene the situation is more challenging. Due to the broken Galilean invariance the dynamics of the center of mass of the electron gas is coupled to the pseudospin degree of freedom 16,39 . This is reflected in the violation of Kohn's theorem in graphene. Therefore one has to take great care when applying the quantum Langevin equation approach to graphene. In particular, in the presence of electron-electron interactions the force-force correlator does not have the simple form of the parabolic case. Therefore the Langevin equation does not readily provide any advantage over the standard Bethe-Salpeter equation.
The situation is greatly simplified in the high frequency regime ωτ 1 where the summation of the infinite series of diagrams in powers of the disorder potential is not required. In this paper we adopt for the case of the Dirac electrons in graphene an equation of motion framework, developed in Refs. 40 and 41, which is suitable for the high frequency regime. This approach is somewhat similar to the derivation of the quantum kinetic equation 39,[42][43][44] , however, here we will employ an equal time correlator, the 'density matrix', as opposed to the quasiclassical Green function used in the former approach.

C. Summary of the calculation and results
In this paper we develop an equation of motion approach in order to calculate the high frequency asymptotic, 1/(ωτ ) 1, (accessible in current experiments 45 ) of the optical conductivity of graphene in the presence of electron-electron interactions and charged impurities. Interaction effects are included in the self-consistent mean field approximation which is equivalent to the random phase approximation (RPA) of the diagrammatic perturbation theory.
We find an enhancement of the disorder scattering rate at finite frequencies due to the presence of the graphene plasmon excitation in the spectrum. We make a quantitative prediction for the frequency-dependent correction to the optical conductivity of graphene due to the plasmon anomaly in the perturbative regime, 1/τ ω E F . The plasmon enhanced scattering gives rise to a broad peak in the frequency-dependent conductivity given by, The dimensionless coefficient c ≡ α/(k F d), depends on the Fermi wave vector k F and for typical values of parameters can be small c 1. Here d is the effective separation of the charged impurities from the graphene plane, and α ≡ e 2 /(κv ) is the dimensionless interaction strength (i.e. the so called graphene fine structure constant) and κ is the effective dielectric constant. The rapid increase of the scattering rate at ω cE F is limited by the flattening of the Coulomb potential of the impurities at the length scale d. Both the location ω * /E F = 3/2c and the magnitude eσ(ω * ) of the peak are sensitive to the characteristics of disorder and the Fermi level E F in the system. Therefore the effect is expected to be strongly sample dependent (i.e. impurity distribution, carrier density and interaction strength).
We use our perturbative high frequency results to make a qualitative prediction for the frequency-dependent optical conductivity at low frequencies ωτ 1. We find that in this regime the plasmon induced frequency-dependent scattering rate leads to a redistribution of the spectral weight in the optical conductivity, which may be substantial in the low density regime in strongly disordered samples. The redistribution of the spectral weight results in the non-Lorentzian shape of the Drude peak which has to be accounted for when extracting various system parameters from measurements of the Drude response.
The paper is organized as follows. In Sec. II we discuss the low energy Hamiltonian of monolayer graphene and the model of Coulomb disorder arising from charged defects. An equation of motion is derived in Sec. III and solved in Sec. IV and details given in Appendix B. The plasmon propagator is introduced in Sec. V. Asymptotic solutions of the equation of motion presented in Sec. VI are used to obtain quantitative predictions for the optical conductivity in Sec. VII. We conclude in Sec. VIII. Appendix A contains a description of the formal theory of the low-energy excitations in graphene that includes a band cut off and Appendix B gives some details of the derivation for the results in Secs. IV and VI.

II. MODEL SETUP
The unit cell of the hexagonal lattice of monolayer graphene contains two chemically equivalent carbon atoms (labeled A and B). The resulting degeneracy ensures the crossing of the conduction and valence bands at the corners of the hexagonal Brillouin zone (labeled K and K ). The low energy excitations (< 1eV ) are described by the two flavor Dirac Hamiltonian, with Σ x/y/z being the 4×4 matrices in the sublattice and valley space that form a Pauli matrix algebra Σ j , AKη and ψ AKη are the creation and annihilation operators for an electron on the sublattice A at the K point, characterized by momentum k and spin η. Below we assume that the Hamiltonian of graphene is diagonal in spin indexes (and a spin degeneracy of 2 as well as a valley degeneracy of 2 will be assumed throughout). Throughout the text we set ≡ 1 and restore it in the final answers presented in Sec. VII, unless explicitly stated. The full low energy Hamiltonian, includes the inter-electron Coulomb interaction, with V q ≡ 2πe 2 /(κq), characterized by a dimensionless ratio of the typical Coulomb energy to the typical kinetic energy α = e 2 κv (here ≡ 1) that is independent of the electron density. The full Hamiltonian also includes the potential of random charged impurities, that dominates the elastic disorder scattering in typical graphene devices fabricated on SiO 2 /Si substrate. Here q is a 2D momentum in the graphene plane, r i is the location of the ith impurity (assumed random) in the 2D plane parallel to the graphene layer, and d i is the distance between the graphene layer and the impurity in the direction perpendicular to the graphene plane. For simplicity, and with no loss of generality we assume that all impurities are located in a 2D plane parallel to the graphene layer located a distance d away from it. We will comment on the implications of the impurities being distributed in 3D later in the paper. We assume without loss of generality the distribution of impurities to be charge-neutral such that the average of the disorder potential over the random disorder configurations vanishes V (i) (r) = 0. We use the Born approximation for the disorder scattering rate such that the disorder scattering is fully characterized by the correlator of the disorder potential averaged over the random disorder configurations, where δ q,q stands for Kronecker symbol. Our model of uncorrelated random disorder can be straightforwardly generalized to correlated disorder scattering 46 if experimental information about disorder correlations is available. To keep the number of parameters a minimum, we assume the impurity-induced random disorder to be completely characterized by just two parameters, the impurity density ρ imp and their location with respect to the graphene layer d, which is the minimal possible model for Coulomb disorder.

III. EQUATION OF MOTION
We introduce an equal time correlator which will be called the 'density matrix' in this paper. The density matrix,ĝ k,k+q (t), is a 8 × 8 matrix defined in terms of the spinor creation/annihilation operatorsΨ † k in the basis of the Bloch wave functions labeled by the sublattice, valley and spin indices. The density matrix is diagonal in the spin indices which will give rise to a factor of η = 2 each time the trace is taken.
In the Heisenberg representation a time derivative of an operator −i∂ tĝk,k+q (t) is given by a commutator with the full Hamiltonian, Calculating the commutators of the density matrix with the Hamiltonian in Eq. (1) (see also Appendix A) we find, and, The interactions are included in the mean field approximation, where ρ(q) ≡ k Trĝ(k, k + q, t) is the Fourier transform of the spatial fluctuation of the electron density.
We include the external electric field in a gauge invariant formulation 47 which requires the monochromatic field E(t) = Ee iωt to enter as follows, where we also take the Fourier transform with respect to time. Here and throughout the text we assume that the frequency ω > 0 has an infinitesimal positive imaginary part which is set to zero at the end in the usual manner.

A. Diagonalized equation of motion
We introduce a projection operator, that projects on the subspace of positive and negative energy eigenstates of the Dirac Hamiltonian H 0 for s = 1 and s = −1, respectively. Here n k ≡ k/k. It can be verified directly that P 2 sn k = P sn k and P n k + P −n k = 1. Multiplying both sides of Eq. (10) by P sn k on the left and by P s n k+q on the right we diagonalize the free particle part Eq. (7) in the right hand side of the equation of motion Eq. (10) , where we used that P sn k vΣ · k = s k P sn k and k ≡ vk. We use an identity, which can be verified directly, in the left hand side of Eq. (10) multiplied by the projectors (as described above) to obtain, Psn k (ieE·∇ kĝk,k+q (ω))P s n k+q where all the terms containing perturbations are collected in the right hand side which will be more convenient for the construction of the perturbation theory.

IV. SOLUTION OF THE EQUATION OF MOTION
The solution of the matrix integro-differential equation (12) can be used to obtain the current response, where j = x, y. Here and throughout the text we take the trace over the sublattice, valley and spin indexes. In the linear response J i t = σ ij E j , i, j = x, y, and therefore it will be sufficient to keep the electric field to the first order in the solution of (12). The conductivity averaged over random disorder configurations is isotropic, nevertheless, it will be more convenient for us to discuss components of the current response and gives the results for the isotropic conductivity at the end.
Solving the equation (12) is still a daunting task as both the disorder potential and interactions have to be included self-consistently. In the present work we solve this equation only in the high frequency regime ωτ 1 where τ is the disorder scattering time. In this high frequency regime it is sufficient to keep the terms up to the 2nd order in the disorder potential. We treat electronelectron interactions in the self-consistent mean field approximation. This approach is equivalent to the standard random phase approximation (RPA) of the manybody perturnation theory which have been shown to provide an accurate description of the plasmon dispersion in graphene 48 . We neglect self-energy and vertex corrections due to Coulomb interactions in the plasmon dispersion because of the relatively weak electronelectron interaction strength in graphene. We emphasize that the long-wavelength square-root in wavevector plasmon dispersion in graphene is protected by current conservation, but the coefficient of this long-wavelength dispersion term is affected (weakly) by interaction effects in graphene (in contrast to parabolic band systems), because of the violation of Galilean invariance in graphene 49 .
We introduce a notation for the expansion of the density matrix, whereĝ (ij) corresponds to the ith order expansion in the electric field and the jth order in the disorder potential, and we drop the subscripts ofĝ k,k+q (ω) in Eq. (14). Note that due to the violation of the Kohn's theorem by the Dirac electrons, the external electric field couples to the homogeneous electron density,ĝ (10) k,k+q (ω) = 0, in contrast to the ordinary parabolic band spectrum 2D electron systems 41 .
In the following we find successive approximations of g k,k+q (ω) defined in (14) by recursively solving Eq. (12) in each order of the perturbation theory. In the zeroth order, free Dirac quasiparticles are described by the Fermi-Dirac distribution f s k = 1 1+e βs k , and the density matrix is proportional to the projector Eq. (11), where β = 1/(k B T ) is the inverse temperature and k B is the Boltzmann constant, and δ k,k+q stands for the Kronecker symbol.

A. First order terms
The linear response of disorder free graphene to an external electric field is described byĝ (10) k,k+q (ω) which satisfies,ĝ k,k+q (ω) = ss Psn k [Hee,ĝ k,k+q (ω)]P s n k+q We first solve Eq. (15) ignoring the Coulomb interaction.
In this case, in the first order in the electric field we obtain, g k,k+q (ω) = ieE 0 · where the sum is over s, s = ±1. From this we directly obtain that the electron density response to the homogeneous electric field, k,k+q (ω) = 0, vanishes in the linear order. Using Eqs. (16) and (13) we reproduce the universal value 13,14 of the conductivity of non-interacting disorder free graphene at the charge neutrality point (E F = 0), where we restore the Planck's constant. At finite doping Eq. (16) describes the optical absorption due to interband transitions with energies ω > 2E F as well as the disorderfree zero-frequency Drude peak 16 , where we restore Planck's constant.
In general, electron-electron interactions give rise to corrections to Eq. (16) reflected in the many-body renormalization of the Drude weight 16 and the renormalization of the universal conductivity σ E F =0 50 . Accounting for these renormalization effects would require including self-energies and vertex corrections in the expansion in the electron-electron interaction strength. We reiterate here that in the present work we neglect all such effects. We include the electron-electron interactions within the self-consistent mean field approximation, which in the first order in the electric field gives the density matrix correction that is proportional to the density response, which therefore vanishes. As a result, the interactions produce no change in Eq. (16) within the chosen approximation scheme in our model. The full inclusion of interaction-induced self-energy and vertex corrections in the graphene dynamical conductivity in the doped situation in the presence of impurity disorder is a formidable task which is well beyond the scope of our work where our interest lies in treating disorder and interaction on an equal footing, treating disorder perturbatively (i.e. highfrequency approximation) and interaction at a mean-field level (i.e. RPA which correctly incorporates the plasmon collective mode in the high-frequency conductivity).
The static density fluctuation induced by the disorder potentialĝ (01) k,k+q satisfies Eq. (12) with the first order terms in disorder included and the terms proportional to the electric field omitted in the right hand side, where s, s = ±1, and we used Eqs. (8) and (9). Taking the trace over the spinor indices and summing over k on both sides of Eq. (17) we solve it for the density ρ (01) (q, 0) ≡ k Trĝ k,k+q (0), Here we introduced the RPA dielectric function, 1/ε(q, ω) ≡ 1/(1 − V −q χ(q, ω)), with the free electron polarization operator defined as 51,52 , η = 2 is the spin degeneracy. Going back to Eq. (17) we obtain the first order correction to the density matrix due to the disorder potential screened by the interacting electron gas, g k,k+q (0) = where the sum is over s, s = ±1.

B. Second order terms
The contribution to the static density fluctuations in the second order in disorderĝ (02) k,k+q vanishes after the averaging over the disorder configurations. This is because the disorder averaged value ĝ (02) k,k+q = ξδ q,0 is homogeneous in space. The charge neutral distribution of disorder does not give rise to non-zero homogeneous corrections to the density, therefore ξ = 0. This conclusion can also be confirmed by showing that in the second order in the disorder potential Eq. (17) does not have non-trivial static homogeneous solutions.
The interplay of the disorder potential and the electric field gives rise to the dynamical fluctuation of the electron density matrixĝ (11) k,k+q (ω) that satisfies the following equation obtained from Eq. (12) in the respective expansion, g (11) k,k+q (ω) = ss We take the trace over the spinor indices and sum over k on both sides of Eq. (21) as before, and solve this equation for the density fluctuation ρ (11) (q, ω) ≡ k Trĝ (11) k,k+q (ω), k,k+q (0)P s n k+q s k+q −s k −ω where the sum is taken over s, s = ±1 and the momentum k. Note that the combination of the left hand side and the last term on the right hand side of Eq. (21) ensure that the density fluctuation ρ (11) (q, ω) is proportional to the inverse of the dynamical dielectric function 1/ε(q, ω).
Using the explicit result for ρ (11) (q, ω) we can find the density matrix correctionĝ (11) k,k+q (ω). The full solutionĝ (11) k,k+q (ω) is straightforward to obtain. The result however is rather cumbersome (and not particularly illuminating) and therefore below we will present only the terms that play a key role in the plasmon enhanced scattering rate which dominates the frequency-dependent optical response at low frequencies 1/τ ω E F .

C. Electric current
Using Eq. (12) we can express the linear response current in terms ofĝ (11) k,k+q (ω) which spares us the necessity of calculatingĝ (12) k,k+q (ω). We multiply Eq. (12) by Σ j (where j = x, y) take the trace in the spinor space and sum over k. Comparing the result to the definition in Eq. (13) we arrive at, The term in the second line in Eq. (23) is proportional to the homogeneous component of the correlatorĝ k,k and therefore does not give any contribution to the disorder induced correction, see Sec. IV B for details. The presence of this free electron contribution is an artifact of the perturbation theory which does not correctly describe the broadening of the singular response at ω = 0. In the following we focus on the high frequency correction to the optical conductivity due to disorder scattering which is accurately described by the first line in Eq. (23), and is denoted J j (ω) in the following. Using Eqs. (8,9) and the cyclic property of trace we write the current response due to the disorder scattering accurate in the second order in disorder as, where the sum is taken over s, s = ±1 and the momenta k and q. The numerator in Eq. (24) contains the bracket (ĝ k−q,k (ω) −ĝ k,k+q (ω)). We shift the sum over momentum k → k +q in the first term of this bracket so that the corresponding difference appears in the factor in front of the bracket denoted Ξ j , J j (ω) = ss ,k,q Tr Ξ j ss ĝ k,k+q (ω) where we introduced, This procedure has to be formally justified by an introduction of a band cut off. The details of this formulation are described in Appendix A. Finally, using the expansion (14) and the result (18) the current response Eq. (25) to the external electric field in the second order in the screened disorder potential is written as, k,k+q (ω)

V. PLASMON PROPAGATOR IN GRAPHENE
The excitation spectrum of the interacting Dirac electron gas consists of the electron-hole continuum and the collective plasmon mode. The latter is described by the pole of the inverse dielectric function given within RPA by 51,52 , Here χ(q, ω) is the free electron polarization function Eq. (19). Eq. (27) has solutions only when ω > vq, that correspond to the plasmon band with the dispersion 51,52 , In this frequency range ω > vq, the imaginary part of the dynamical dielectric function is proportional to the plasmon propagator, In the relativistic electron gas the Landau damping is absent, i.e. mχ(q, ω) = 0, in the frequency range vq < ω < 2E F − vq 51,52 . Therefore within this range we can approximate the imaginary part of the plasmon propagator by a delta function neglecting the relatively weak damping by disorder, The plasmon momentum is restricted to vq < ω which therefore allows an expansion in both the frequency ω E F and momentum q k F in the prefactor in front of the plasmon propagator in Eq. (26). At higher frequencies, when the condition ω + vq > 2E F is satisfied, the plasmon is damped by electron-hole excitations, and the frequency-dependent plasmon width (i.e. Landau damping) has to be taken into account. The plasmon damping does not give rise to any qualitatively different behavior. Therefore we consider only the delta-function form of the plasmon propagator which gives quantitatively accurate results at frequencies ω E F . Electronhole excitations determine the behavior of the dielectric function for ω < vq, At low frequencies ω E F the frequency dependence of the dynamical screening due to electron-hole excitations given by Eq. (30) is weak. Therefore the frequency dependence of the current response is dominated by the terms proportional to the plasmon propagator in Eq. (29). This allows us to drop the frequency dependence in all other terms in the current response, where J j pl (ω) is the frequency-dependent part of Eq. (26) that contains the plasmon propagator, where we kept only the last term in Eq. (21) forĝ (11) k,k+q (ω) which is the only one containing the plasmon propagator. The sum in the first line of Eq. (31) is over s, s = ±1 and the momenta k and q, in the second line the sum is over t, t = ±1.
The final analytical expression for the plasmon enhanced optical conductivity is rather cumbersome and therefore is left for the Appendix B. In the next section we calculate the low frequency asymptotic of the solution as well as the frequency-independent response J j (ω = 0).

A. Intraband solution of the equation of motion
At low frequencies and high doping ω/E F 1 both the static scattering rate and the frequency-dependent plasmon induced correction to it are dominated by the intraband processes. Therefore in the right hand side of Eq. (12) we can neglect all interband terms that contain large denominators ≈ 2E F in the low frequency limit, k,k+q (0) and ρ (11) (q, ω) and using Eq. (26) and performing some tedious but straightforward algebraic calculations we arrive at, where, and j = x, y. We introduced above the correlation functions, and Here we introduced a notation, The first term in Eq. (33) reproduces the high frequency expansion of the standard Drude conductivity, with the frequency-independent transport relaxation rate due to impurity scattering. Here ν is the density of states of the Dirac electrons per spin and valley. To demonstrate the relation between Eqs. (34) and (39) explicitly we keep only the main term in the low frequency expansion of Eq. (36), Substituting the above result into Eq. (34) and comparing it to Eq. (39) we reproduce the standard definition 1 of the static elastic scattering rate due to screened Coulomb disorder in graphene, where we restore the Planck's constant. The second term in Eq. (33) describes the frequencydependent contribution to the optical conductivity due to the plasmon anomaly. This part can be simplified by taking the low frequency asymptotic limit of Eq. (37), Taking the integral over the angle ϕ we arrive at, where Ω ≡ vq ω and ζ = 3/2. Therefore the low frequency expansion of Eq. (35) with the use of Eq. (41) gives,  Figure 1. The real part of the conductivity (in units of conductance quantum) approximated by the sum of the Drude peak with the frequency-independent scattering rate and the plasmon peak Eq. (43). The Drude peak corresponds to the very narrow feature at ω = 0. The dashed line shows the limit of applicability of the perturbative expansion ωτ ∼ 1. The impurity density is ρimp = 6 × 10 12 cm −2 located at d = 10 × 10 −9 m, the electron density varies from 10 12 cm −2 to 3 × 10 12 cm −2 top to bottom on the right, the strength of interactions is α = 0.58 corresponding to graphene on hBN.
B. Low frequency, ω EF , asymptotic limit of the full solution We can verify the intraband result (35) by taking the low frequency asymptotic limit of the plasmon anomaly induced correction to the optical conductivity (B1) obtained by solving the full equation of motion valid at frequencies ω < 2E F . We find that the result of this low frequency expansion coincides with Eq. (42).
The result (42) is somewhat analogous to the one derived for the case of the 2D electron gas with the parabolic dispersion 34 , albeit with a different numerical coefficient. However, in the case of graphene Eq. (42) is a low frequency expansion that describes only the frequency-dependent scattering rate due to the plasmon anomaly, the frequency dependence of screening due to the electron-hole excitations is neglected. At higher frequencies ω ∼ E F Eq. (B1) provides a more quantitatively accurate description. At yet higher frequencies ω 2E F the solution of Eq. (21) with (26) has to be used including all interband effects.

VII. RESULTS
A. High frequency regime: ωτ 1 Substituting the plasmon propagator Eq. (29) into the expression for the optical conductivity Eq. (42) and taking the trivial integral over the momentum q which is equivalent to the substitution q → q ω ≡ ω 2 /(2αvE F ), we calculate the frequency-dependent part of the optical conductivity (in units of conductance quantum σ 0 ≡ 2e 2 /h), 40 pl Ω Figure 2. The ratio of the static transport scattering time τ (0) to the frequency-dependent scattering rate due to plasmon anomaly τ pl (ω) for graphene. The parameters are, d = 20nm, the dimensionless interaction strength is α = 0.58, the impurity density ρimp = 3 × 10 12 cm −2 . Different lines correspond to the different electron densities from ρ = 0.5 × 10 12 cm −2 to ρ = 5.5 × 10 12 cm −2 top to bottom on the right. Maximum of the scattering rate ωmax = EF 5α/(2kF d) is dependent on the electron density in the system. (a) and (b) correspond to the imaginary and real part of the complex frequency dependent scattering rate 1/τ pl (ω), respectively.
Here and throughout this section we restore the Planck's constant. Fig. 1 shows Eq. (33) using Eq. (43) and the Lorentzian peak eσ D (ω) ≈ eσ Drude (ω) with the static elastic relaxation rate given by Eq. (40). The effect of the plasmon pole on the optical conductivity is weak at lower frequencies ω cE F , however, at higher frequencies ω cE F there is a sharp increase in the optical conductivity with increasing frequency that is limited by the flattening of the potential of charged impurities as determined by their position d. A consequence of this behavior is a broad peak in the frequency-dependent optical conductivity with the maximum located at, and the width defined by the solution of the equation as shown in Fig. 1.
We define a frequency-dependent relaxation rate associated with the plasmon enhanced dissipation using the Drude form Eq. (39) and comparing it to Eq. (43) we get, In a wide range of parameters the plasmon enhanced scattering rate varies greatly with frequency and can be well over an order of magnitude larger than the static transport scattering rate, see Fig. 2(a). This strong frequency dependence of the imaginary part of the relaxation rate gives rise to a non-zero real part shown in Fig. 2(b). The strong dynamical variation of the scattering rate leads to the spectral weight redistribution in the optical conductivity which will be particularly important in the strongly disordered samples at low densities for which the above perturbation theory is expected to fail. This regime will be discussed in the following section.
In the following we take into account the finite spread of the spatial distribution of impurities along the z-axis perpendicular to the graphene layer. Including the distribution spread requires the additional averaging of Eq. (43) over positions of impurities along the z-axis which results in an overall suppression of the magnitude of the effect of the plasmon enhanced scattering. For comparison, we calculate the plasmon enhanced frequencydependent optical conductivity for the case of a uniform density of impurities along Oz in a layer d 1 ≤ z ≤ d 2 , so that the 3D impurity density equals where c 1 = α/(k F d 1 ). Note that the 3D impurity distribution is reflected in a different power in the frequencydependent optical conductivity and therefore in a different shape of the plasmon induced peak. The maximum of the peak given by Eq. (47) is suppressed with respect to the peak maximum of Eq. (43) which is not necessarily a very small number for typical parameters. Therefore we expect the effect of plasmon enhanced scattering to be observable even in the case of the impurities distributed three-dimensionally throughout the substrate layer in the typical field effect geometry of the graphene based devices.
B. Low frequency regime: ωτ 1 In the following we extend the perturbative results obtained above into the low frequency regime, ωτ 1. This procedure allows us to make a qualitative prediction about the low density behavior of strongly disordered samples. In this low-frequency regime the spectral weight redistribution in the optical conductivity is expected to be substantial and the shape of the Drude peak is expected to be non-Lorentzian.
The frequency-dependent scattering is described by a generalized "self-energy" or memory function M (ω) that reflects the divergent nature of the response of the electron gas at ω = 0 , The real part of the optical conductivity then reads, In the limit eM ω 1 and mM ω 1 we can match the expression for the memory function with the perturbative solution found above, where the right hand side is given by Eq. (33). On the other hand at ω = 0 the scattering rate in the perturbative expression for the optical conductivity in Eq. (33) coincides with the static elastic transport scattering rate Eq. (40). Note that the contribution due to the plasmon anomaly Eq. (43) vanishes with ω → 0. The memory function is therefore a function of frequency that equals mM (0) = 1 τ (0) at ω = 0 and is given by Eq. (48) at ωτ 1. Therefore using the expansion (48) with the right hand side given by Eq. (33) in the whole range of frequencies 0 ≤ ω E F provides a reasonable interpolation for the behavior of the low-frequency optical conductivity. This interpolation would fail in the presence of any divergence in the memory function at low frequencies ωτ 1. We do not anticipate such a divergence at frequencies of interest here ω 2E F (where the interband processes are not expected to play a role).
The result of the low-frequency interpolation procedure is presented in Fig. 3, which shows the strong evolution of the Drude peak shape with the changing density (that corresponds to the different lines on the plot). In particular, the top (blue) line in Fig. 3 demonstrates the importance of the non-zero real part of the memory function which gives rise to the strong redistribution of the spectral weight in the optical conductivity and the resulting non-trivial shape of the tail of the Drude peak.

VIII. DISCUSSION
In summary, we have developed an equation of motion approach to describe the linear response of Dirac electrons in the presence of electron-electron interaction and disorder scattering due to charged impurities. We obtain quantitative predictions for the optical conductivity in the high frequency regime ωτ 1. We show that the presence of a plasmon pole in the dynamical dielectric function leads to a strong enhancement of the dissipation at finite frequency reflected in a wide peak in the real part of the optical conductivity. Characteristics of the predicted peak feature are strongly dependent on the strength of the electron-electron interaction, the location of the impurities with respect to the 2D Dirac electron gas, and can be tuned by tuning the Fermi level by changing the carrier density. Extrapolating to the low frequency regime ωτ 1 we find that the plasmon enhanced scattering may determine the non-Lorentzian shape of the Drude response.
The theory developed here is quantitatively accurate only in the high frequency regime ωτ 1. An accurate description at low frequencies must be based on the solution of the Bethe-Salpeter equation for the vertex function or some analog of it, which is beyond the scope of the current work where we focus on the interplay between disorder and interaction perturbatively (thus necessitating the high-frequency perturbative expansion). We have, however, provided a low-frequency extrapolation of our theory which should have qualitative, if not quantitative, validity, as long as the relevant memory function does not have a singularity at low frequencies. This extrapolation is really an interpolation based on an extension of our high-frequency theory to lower frequencies by using the known zero-frequency Drude result.
We have neglected effects of the many body renormalization of the electronic spectrum and exciton energies by electron-electron interactions. Within the parameter range considered in this paper, high Fermi level and relatively low frequencies 1/τ ω E F we expect the renormalization effects to give only a small quantitative rather than qualitative corrections to the effect of plasmon enhanced scattering described here. Thus, our main qualitative prediction of the broad peak and the associated non-Lorentzian line shape of the optical conductivity should remain valid independent of our neglect of the many-body reconstruction of the Dirac spectrum which is very small at finite carrier densities by virtue of the fine structure constant in graphene (particularly for graphene on substrates) being not too large (of the order of unity or less) 53 . Length scales involved in the plasmon enhanced scattering are short compared to the mean free path and therefore the disorder related quantum interference ef-fects give no contribution to the phenomenon discussed in this paper.
The relatively large magnitude of the finite frequency peak in the optical conductivity may enable an experimental confirmation of our predictions. Moreover, the possibility to tune the position of the plasmon induced peak in the optical conductivity of graphene, controlled by the parameter c = α/(k F d), by tuning the Fermi level in the system allows an unambiguous identification of the plasmon enhanced dissipation effect in the optical response, which can otherwise be masked by the interband transitions, scattering by phonons and other inelastic processes. Also, varying the location of impurities d allows additional tunability of the effect. This can be achieved using hexagonal boron nitride spacers of different thicknesses d between the graphene layer and SiO 2 /Si substrate with the latter containing the majority of charged impurities in the vicinity of its surface. In particular, samples with relatively wide hBN spacers d ∼ 10nm would correspond to c 1 at densities n ∼ 10 12 cm −2 in which case the plasmon induced peak appears at the relatively low frequency ω ∼ cE F unobstructed by the interband transitions. The plasmon enhanced dissipation discussed here is a generic disordered Fermi liquid phenomenon which nevertheless has not been unambiguously identified experimentally so far, and graphene is an ideal system to look for signatures of this phenomenon due to graphene's high electronic quality and tunability of various system parameters, and the strong magnitude of the plasmon enhanced scattering in graphene.
With regards to the existing optical measurements the deviation of the low frequency response from the Lorentzian shape has been reported by some 17,19 but not all 3,12 experiments indicating that the effect is sample dependent which is consistent with the mechanism (i.e. disorder plus plasmon) predicted here. Redistribution of the spectral weight in the optical conductivity and the resulting non-Lorentzian line shape of the Drude peak, if present, have to be accounted for when extracting the Drude weight that have been recently used to analyze the spectrum renormalization due to electron-electron interactions 3,12 .
The theoretical approach developed here is generic and can be straightforwardly generalized to the case of the 2D Dirac fermions observed in other materials such as topological insulators and 2D dichalcongenites. The most important qualitative message of our theory is that the finite-frequency dynamical conductivity appropriate for optical measurements in doped graphene cannot be approximated to have a simple constant disorder broadening given by the corresponding static transport scattering rate because of the substantial interplay between plasmons and impurity scattering. In particular, the finitefrequency scattering time appropriate for the optical conductivity manifests very nontrivial frequency dependence which could produce unexpected qualitative phenomena such as the broad finite-frequency peak in the conductiv-ity with both the peak position and the peak width being determined nontrivially by both impurity scattering and Fermi energy.