Pumping electrons in graphene to the M -point in the Brillouin zone: The emergence of anisotropic plasmons
PPumping electrons in graphene to the M -point in the Brillouin zone: The emergence ofanisotropic plasmons A. J. Chaves ∗ and N. M. R. Peres † Department of Physics and Center of Physics, University of Minho, P-4710-057, Braga, Portugal
Tony Low ‡ Department of Electrical and Computer Engineering,University of Minnesota, Minneapolis, Minnesota 55455, USA (Dated: November 7, 2018)We consider the existence of plasmons in a non-equilibrium situation where electrons from thevalence band of graphene are pumped to states in the Brillouin zone around the M -point by a highintensity UV electromagnetic field. The resulting out-of-equilibrium electron gas is later probed by aweak electromagnetic field of different frequency. We show that the optical properties of the systemand the dispersion of the plasmons are strongly anisotropic, depending on the pumping radiationproperties: its intensity, polarization, and frequency. This anisotropy has its roots in the saddle-likenature of the electronic dispersion relation around that particular point in the Brillouin zone. It is foundthat despite the strong anisotropy, the dispersion of the plasmons scales with the square root of thewave number but is characterized an effective Fermi energy, which depends on the properties of thepumping radiation. Our calculations go beyond the usual Dirac cone approximation taking the fullband structure of graphene into account. This is a necessary condition for discussing plasmons at the M -point in the Brillouin zone. I. INTRODUCTION
Light matter interaction at low energies occurs due tothe interaction of electromagnetic radiation with weaklybound electrons in a material. The electron gas in a con-ductor is a well known form of weakly bound electronswhich couple to electromagnetic radiation. Usually, weassume that the electron gas is in equilibrium and thatthe external perturbation is small, in which case the re-sponse of the gas can be computed in terms of its equi-librium properties —linear response theory. The situa-tion is markedly different when a high-intensity electro-magnetic radiation interacts with an electron gas drivingit to an out-of-equilibrium regime. In this case the re-sponse of the system depends on intensity of incomingradiation and on the orientation of its polarization rela-tively to the real space lattice of the crystal. Moreoverthe distribution function of the electron gas occupancy isno longer a Fermi-Dirac distribution. It is this situationthat will be studied in this paper. Here we consider anintense pumping electromagnetic field interacting withthe weakly bound electrons in graphene thus generatingan out-of-equilibrium gas. The pumping is followed bya weak-probe electromagnetic-field, of frequency muchsmaller than that of the pumping field, which probes theout-of-equilibrium plasma created by the pumping. Thephysics of this process is depicted in Fig. 1.Pumping graphene with electromagnetic radiation is apossible tool to study the dynamics of the charge carri-ers in graphene. In the work of George et al. the recom-bination dynamics and carrier relaxation in graphenewas studied with terahertz spectroscopy. When elec-trons in graphene are pumped by an intense light fieldpulse, the response is highly anisotropic as was shownin Ref. [2]. After about 1ps of the initial pump- FIG. 1. (Color on-line.) This figure represents the physical situationwe are considering in this paper: a possibly doped graphene is drivenout-of-equilibrium by an electric field of frequency ω p that creates anelectron gas around the M -point in the Brillouin zone. The createdplasma is then probed by a field of frequency ω (cid:28) ω p . We note thatthe mechanism we are discussing does not require an initially dopedgraphene, that is, we can have E F = 0 . Also note that the electronicspectrum at the M -point has a saddle-like nature (in the drawing onlythe steepest descent direction is shown; see Fig. 17 for a drawing ofthe full band structure). The probe field allows the excitation of anelectron, belonging to the out-of-equilibrium gas, to higher energies. ing pulse, the photogenerated electrons are describedby an isotropic Fermi-Dirac distribution with a hightemperature . The graphene optical properties undersuch conditions were studied by Malic et al. and Sun e t al. . The electron dynamics of photo-excited elec-trons, including the stimulated electron-hole recombina-tion, was studied by Li et al. . The plasmon dispersionrelation under a non-equilibrium hot Fermi-Dirac distri-bution was studied very recently by Page et al. and ex-perimentally measured by Ni el al. . We note that plas- a r X i v : . [ c ond - m a t . m e s - h a ll ] O c t mons in graphene also offer a possible decay channelto cool down the hot electron gas . The optical con-ductivity of doped and gapped graphene taking into ac-count a non-equilibrium distribution and interband pro-cesses was studied theoretically by Singh et al. . Allthese studies where made in a regime where the valid-ity of the Dirac cone approximation holds, that is, thepumping field has a frequency in the IR/Vis region of theelectromagnetic spectrum.For a pulsed laser beam with a pulse duration muchlarger than 1 ps —the case we will consider in thiswork—, the carrier distribution will remain anisotropicfor the duration of the pulse. In this case, electron-phonon and electron-electron interactions, as well as theeffect of disorder can be encoded into a hot carrier relax-ation rate. Under such conditions the surface plasmon-polariton (SPP) was studied in systems described by agapped Dirac equation by Kumar et al. , who discussedthe response of the electron gas to a circular polarizedlight. In this approach, the density matrix equations ofmotion are solved to determine the non-equilibrium elec-tronic distribution. In this paper we study the electronicdistribution and the plasmon spectrum of an opticallypumped graphene. The material is subjected to a lin-early polarized light beam, with a frequency that createsan electron gas beyond the regime where the Dirac ap-proximation holds. This is relevant when graphene issubjected to UV radiation. In this case, the spectrumis no longer Dirac-like and the full band structure of thesystem has to be taken into account. Therefore, thiswork goes beyond that of Anshuman et al. and in-cludes also the regime studied by these authors.The plasmons in graphene were first probed in real-space in the studies of Fei et al. and Chen et al. .In the work [8] the plasmons in graphene are generatedby an infrared beam focused in the metalized tip of anatomic force microscope, after a first pumping pulse ofelectromagnetic radiation. The tip is also used to detectthe plasmons that propagate along the graphene sur-face and after reflection in the sample edges standingwaves are produced. This kind of experiment can beused to detect the plasmons discussed in the presentwork using a pulsed laser in the UV range (pulse du-ration much larger than 1 ps) for pumping the elec-trons in graphene. Due to excitonic effects the posi-tion of the maximum of absorption associated to inter-band transitions at the M − point is reduced from . eV(independent-electron result) to about ∼ . eV ( λ ∼ nm), a wavelength for which there are available lasers(see also Sec. VIII).Under intense and long optical excitation, the carrierdistribution maintains a non-equilibrium state and doesnot follow the Fermi Dirac distribution. The new elec-tronic distribution has to be calculated using the von-Neumann equation of motion, with a phenomenologi-cal relaxation-term that tends to drive the system to-wards thermal equilibrium, characterized by a Fermi-Dirac distribution. Since we want to discuss plasmons in the non-equilibrium electron gas created around the M -point in the Brillouin zone, we need to describe the π − electrons in graphene using a tight-binding Hamilto-nian. For graphene in the tight-binding approximation,we have a two-band (valence and conduction) systemlabeled by the crystal momentum k , which runs over thefull hexagonal Brillouin zone. In our calculations, car-rier scattering is accounted for via a relaxation rate. Assuch, the resulting equations need to be solved for eachpoint in the first Brillouin zone and different momentumvalues are not explicitly coupled.The paper is organized as follows: in Sec. II we de-rive the Bloch equations for the electrons in graphenewithin the tight-binding model, which is valid beyond theDirac cone approximation. We study the transient re-sponse under a pulse laser and obtain analytical expres-sions for the out-of-equilibrium electronic distribution. InSec. III we obtain the equations to calculate the out-of-equilibrium susceptibility from the new electronic distri-bution. In Sec. IV we derive a semi-analytical formula,valid in the long-wavelength regime, for the susceptibil-ity and in Sec. V the optical conductivity is obtained inthe same conditions. In section VI we compute numer-ically the susceptibility of the out-of-equilibrium electrongas and obtain results for the plasmon dispersion andthe loss function. In Sec. VII we use the semi-analyticalequations for the conductivity tensor to numerically cal-culate the relation dispersion of the surface plasmon-polariton and discuss its anisotropic properties. In Sec.VIII we provide a discussion on the feasibility of an ex-periment to observed the predicted effects and the gen-eral conclusions of our work. II. NON-EQUILIBRIUM DENSITY-MATRIX AND BLOCHEQUATION IN THE TIGHT-BINDING APPROXIMATION
We consider the electrons in graphene described by atight-binding Hamiltonian H , with a nearest neighborshopping term only, and subjected to an external electricfield E (see Appendix A): H = H + V = (cid:88) n t TB (cid:16) ˆ a † n ˆ b n + ˆ b † n ˆ a n (cid:17) + e E · R , (1)where the index n extends over all unit cells of the crys-tal, t TB = 2 . eV is the hopping parameter, ˆ a † n ( ˆ b † n ) cre-ates an electron in the site n of the sublattice A(B), and R = R A + R B is the position operator: R A = (cid:88) n R n ˆ a † n ˆ a n , (2a) R B = (cid:88) n ( R n + δ ) ˆ b † n ˆ b n , (2b)where δ = a (1 , is the vector connecting the A and B sub-lattices in the same primitive cell, and R n = n a + n a , n and n are integers numbers, and theprimitive vectors are defined as: a = a (3 / , √ / and a = a (3 / , −√ / , where a ≈ . nm is the carbon-carbon distance in graphene.The eigenvectors of H , satisfying the eigenvalueequation H | k , λ (cid:105) = λE k | k , λ (cid:105) , are: | k , λ (cid:105) = 1 √ (cid:0) | A, k (cid:105) + λe i Θ k | B, k (cid:105) (cid:1) , (3)with λ = ± and | A, k (cid:105) = (cid:80) n e i k · R n ˆ a † n | (cid:105) (thesame holds for the states | B, k (cid:105) upon the replacement { a, A } ↔ { b, B } ). We define the useful auxiliary func-tions φ k and Θ k through the relation: φ k = (cid:88) i e i k · δ i = | φ k | e i Θ k , (4)with the positive eigenvalues given by E k = t TB | φ k | andwhere δ i are the vectors connecting an atom A to all thethree nearest neighbor B atoms [see Eq. (D2)]. Fromnow on, all the momenta are measured in units of theinverse of the lattice parameter a , thus we make thereplacement k → k a .The pumping of the electrons in graphene changesthe electronic density. Thus the system is describedby a non-equilibrium —but stationary— distribution. Tocalculate the system properties in an out-of-equilibriumregime we use the density matrix ρ ( t ) formalism, whosetime-evolution obeys the von-Neumann equation : i (cid:126) ∂ t ρ ( t ) = [ H, ρ ( t )] . (5)We assume a time-dependent and uniform electricfield (that is, with null in-plane wave number). Thus theinteraction term in the Hamiltonian does not couple elec-tronic states from different points in the Brillouin zone.We then project Eq. (5) into the eigenvectors of H givenby Eq. (3) with the same wave number k and differentbands. The diagonal part of the density matrix corre-sponds to the new distribution functions: n c, k ( t ) = (cid:104) k , + | ρ ( t ) | k , + (cid:105) , (6) n v, k ( t ) = (cid:104) k , −| ρ ( t ) | k , −(cid:105) , (7)and the off-diagonal elements correspond to transitionsprobabilities: p cv, k ( t ) = (cid:104) k , + | ρ ( t ) | k , −(cid:105) , (8) p vc, k ( t ) = (cid:104) k , −| ρ ( t ) | k , + (cid:105) . (9)After the calculation of the commutator in Eq. (5) andthe introduction of two phenomenological relaxation-rateterms responsible for relaxing the non-equilibrium elec-tron gas back to its equilibrium Fermi-Dirac distribution,we obtain a set of coupled equations: − ∂ t n c, k = γ ( n c, k − f c, k ) + i Ω k ( t )∆ p k , (10a) − ∂ t n v, k = γ ( n v, k − f v, k ) − i Ω k ( t )∆ p k , (10b) ( ∂ t + iω k + γ p ) p cv, k = − i Ω k ( t )∆ n k , (10c) ( ∂ t − iω k + γ p ) p vc, k = i Ω k ( t )∆ n k , (10d) where the interband relaxation rate is represented by γ p and intraband one by γ ; also we have (cid:126) ω k =2 E k , ∆ n k = n c, k − n v, k , ∆ p k = p cv, k − p vc, k , and f c/v, k is the equilibrium Fermi-Distribution for the con-duction/valence band. The time dependence of n c/v, k , p vc/cv, k , and E has been omitted for simplicity of nota-tion, and finally the Rabi frequency Ω k is given by: Ω k ( t ) = ea E ( t ) · ∇ k Θ k (cid:126) , (11)and couples the diagonal to the off-diagonal elements ofthe density matrix through the external pumping electricfield E ( t ) .Equations (10) are the Bloch equations in graphene ,with only interband contributions included (note that wewant to excite electrons deep in the valence band to highup in the condution band), and describe the evolution ofthe electronic distribution and the rate of interband tran-sitions when an external intense and highly energeticelectric field E is applied. The vector field ∇ k Θ k en-tering in the Rabi frequency does not depend on theexternal parameters and is shown in Fig. 2. The twoinequivalents Dirac points are located at the corners ofthe hexagon in this figure. Near these points the func-tion Θ k becomes the angle between the momentum k and the x -axis.We now comment on the possible values of γ (intra-band scattering rate; note that this controls what is quan-tum optics is called the population ) and γ p (interbandscattering rate; note that this controls what in quantumoptics is called the coherence ). Let us first remark thatthese quantities have been scarcely studied in the UVrange ; probably the most comprehensive study isthat of Oum et al. . Roberts et al. suggest an in-traband electron-electron scattering time τ e − e < ps(the interband electron-electron scattering time has notbeen measured), followed by an electron-phonon intra-band scattering time larger than τ > ps. Theoreti-cally, these scattering times have not yet been studiedat the M − point. On the other hand, the carrier dynam-ics is much better studied when the charge carriers areexcited with IR/Vis radiation. The dynamics after the ini-tial pulse time of duration ∆ t p is the following: intrabandelectron-electron collisions leads the out-of-equilibriumelectron gas to a quasi-thermal and transient distribu-tion after a time τ th . This distribution is characterized bya temperature T el that essentially controls the broaden-ing of the optical conductivity features at high frequen-cies. For longer times the system relaxes towards ther-mal equilibrium via electron-phonon coupling, occurringin a time scale τ C (cooling time scale), finally recombi-nation of electron-hole pairs, occurring in a time scale τ R , takes place. For a pumping-field of photon energy (cid:126) ω p = 1 . eV George et al. identified three time scales:rapid thermalization, τ e − e ∼ − fs (smaller than 60fs, according to Ref. [17], and around 10 fs, accordingto Ref. [20], measured using a Z − scan technique), fol-lowed by carrier cooling on a time scale of τ C ∼ . − FIG. 2. (Color on-line.) Plot of the vector field ∇ k Θ k / || ∇ k Θ k || .This field controls the Rabi frequency and it can be probed by thepolarization of the pumping field. Note that the rotation of the vec-tor field in the two non-equivalent Dirac points has opposite senses.The red hexagon represents the first Brillouin zone and the inten-sity refers to the absolute value of the vector field ∇ k Θ k ; it is moreintense around the Dirac points (brighter spots) and along the direc-tions connecting two Dirac points and passing through the M − point. ps ( (cid:126) /τ C ∼ − meV; τ C ∼ . ps, according toRef. [17]), and finally carrier recombination takes placecharacterized by a time scale τ R ∼ − ps (for smalldoping, theoretical estimations give τ R > ps ). Weassociate the time τ C with /γ , with τ C ∼ . − ps.The assignment of γ p to a given relaxation rate is moredifficult. Intuition dictates that τ p should be determinedby interband electron-electron interactions, a quantitynot easily accessible via pump-probe experiments. Wetherefore concluded that τ R in the IR/Vis/UV range islargely undetermined at present. Luminescence stud-ies of graphene irradiated with Vis/UV electromagneticpulses (30 fs duration) found a characteristic emis-sion time of τ em ∼ . − . ps ( (cid:126) /τ em ∼ − meV). Assuming that the luminescence transition is con-trolled by γ p we can consider the longer time of 0.1 psto estimate (cid:126) γ p ∼ meV (in the figures we shall use (cid:126) γ p = 2 γ = 28 meV). A. Real Time Analysis
We can simplify the set of Eqs. (10), defined in termsof two real ( n c, k , n v, k ) and two complex p vc, k , p cv, k quantities, to three equations involving real quantities only. We will show that under an intense monochromaticwave, the electronic density can reach a steady-state.Using Eqs. (B5) and (B6), we define the deviationfrom the equilibrium density as ρ k ( t ) : n c, k ( t ) = f c, k + ρ k ( t ) , (12a) n v, k ( t ) = f v, k − ρ k ( t ) , (12b)and we split the transition rate into real and imaginaryparts: p vc, k ( t ) = x k ( t ) + iy k ( t ) , (13)where x k ( t ) and y k ( t ) are real and p cv, k ( t ) = x k ( t ) − iy k ( t ) . From Eqs. (10) (see Appendix B), we can write: ˙ x k = − γ p x k − ω k y k , (14a) ˙ y k = ω k x k − γ p y k − Ω k ( t ) (2 ρ k + ∆ f k ) , (14b) ˙ ρ k = − γ ρ k + 2Ω k ( t ) y k , (14c)where the time dependence in x , y and ρ has been omit-ted and ∆ f k = f v, k − f c, k . Also note the different signsin front of ρ k ( t ) in Eqs. (12). FIG. 3. (Color on-line.) Plot of the time evolution of the electronicdistribution ρ k ( t ) for k = 2 π/a (3 , , (cid:126) ω p = 2 t TB , τ = 300 fs(see Ref. [23]), which corresponds to (cid:126) γ = 14 meV (see also Ref.[24]), (cid:126) γ p = 28 meV, E = 0 . GV/m, and E F = 0 . eV. Thepumping field is linearly polarized along the x -axis. The steady stateis attained after 1 ps. The set of coupled Eqs. (14) describe, using threereal functions x , y , and ρ , both the time evolution ofthe transition probability and the electronic density inthe reciprocal space. Since we have included the ef-fect of both electron-electron and electron-phonon inter-actions using only a constant relaxation rate, there isno coupling between excitations from different k . Thus,for each point in the Brillouin zone we can solve Eqs.(14). In Fig. 3 we plot the time evolution of the function ρ k , for k = 2 π/a (3 , (that corresponds to the M -pointin the Brillouin zone) for a monochromatic electric fieldof frequency (cid:126) ω p = 2 t TB (that corresponds to a verticaltransition at the M − points) and intensity E = 0 . GV/m(this is a moderate field intensity), with linear polarizationalong the x-axis.
B. Steady-State Solution
As shown in Fig. 3, the electronic distribution con-verges to a well defined value which we calculate in Ap-pendix C, when the electric field is given by a monochro-matic wave with pumping frequency ω p and intensity E : E = E e iω p t + h.c.. In this case the steady-state solu-tion for the densities can be written as: n c, k = (1 + α k ) f c, k + α k f v, k α k , (15a) n v, k = (1 + α k ) f v, k + α k f c, k α k , (15b)and from Eq. (12) ρ k = α k α k ( f v, k − f c, k ) , (16)where α k is given by Eq. (C5) and ρ k represents the de-viation from the equilibrium Fermi-Dirac distribution. Itshould be noted that the steady state distribution func-tions are not given by a Fermi-Dirac distribution, but canbe written in terms of combinations of f c, k and f v, k .The non-equilibrium distribution depends on the pump-ing frequency ω p and on the complex electric field E .For the linear polarization the electric field can be writ-ten as a real quantity that depends on the intensity ofthe electric field E and the angle of polarization θ .The distribution ρ k is plotted in Fig. 4 for differentpumping frequencies ω p and polarization angles θ of thepumping field. As the pumping frequency increases, theelectronic distribution departs from the Dirac points (thecorners of the blue hexagon). At (cid:126) ω p = 2 t TB , the M -point is populated. For (cid:126) ω p > t TB , the electronic dis-tribution becomes a circle around the Γ -point (center ofthe hexagon). We also see in this figure the polarizationdependence coming from the term ∇ k Θ k · E . For exam-ple, although we have three independent M -points, for (cid:126) ω p = 2 t TB and θ = π/ only two are populated. We canuse Fig. 2 to predict, for a given pumping polarization,what points in the Brillouin zone can be optically popu-lated, noting that the electric field E and the vector field ∇ k Θ k need to be parallel to maximize the electronic oc-cupation. This anisotropy in the population of the M -points is at the heart of other anisotropic effects that wewill discuss ahead. III. INTRABAND TRANSITIONS OF THENON-EQUILIBRIUM GAS DUE TO THE PROBE FIELD
The optical response of graphene is determined byintraband and interband transitions . As shown in theprevious section, in the steady state the pumping fieldchanges the electronic distribution and therefore the op-tical conductivity of the material. This quantity is relatedto the charge-charge correlation function of graphene.The charge-charge correlation function can be calcu-lated using the new electronic distribution obtained inEqs. (15) and (16), instead of the equilibrium Fermi-Dirac distribution, as: χ ( q , ω ) = 2 e (cid:126) a (cid:88) λ,λ (cid:48) = ± (cid:90) ◦ BZ d k (2 π ) n λ k + q − n λ (cid:48) k ω − ω λ,λ (cid:48) k , q + iε N k , q λ (cid:48) λ , (17)where the factor accounts for the spin degeneracy, n λ k is the electronic distribution given by Eq. (15), (cid:126) ω λ (cid:48) ,λ k , q = λ (cid:48) E k − λE k + q is the energy transition, and the overlapof the eigenfunctions is given by: N k , q λ (cid:48) λ = 12 [1 + λ (cid:48) λ cos (Θ k − Θ k + q )] , (18)Using Eq. (12) to split the density into the equilib-rium f c/v k and fluctuation ρ k parts, the susceptibility canbe decomposed into an equilibrium χ ( q , ω ) part, thatis calculated using the equilibrium Fermi-Dirac distribu-tion, and two pumped components, one intraband andthe other interband, as: χ ( q , ω ) = χ ( q , ω ) + χ intrapump ( q , ω ) + χ interpump ( q , ω ) , (19)where ω is the frequency of the probe. The intrabandpumped component of the susceptibility reads χ intrapump ( q , ω ) = 2 e (cid:126) a (cid:88) λ = ± (cid:90) ◦ BZ d k (2 π ) λ ( ρ k + q − ρ k ) ω − λω k , q + iε N k , q λ,λ , (20)where ω k , q = E k − E k + q . Note that χ intrapump is deter-mined by the deviations to the Fermi-Dirac distribution.Since we are only interested in the physics of the elec-tron gas created in the conduction band, we neglect inthe following the contribution from the interband transi-tions χ interpump ( q , ω ) to the total susceptibility, an assump-tion that is valid when ω p (cid:29) ω . This is always the casein this work as we are considering pumping to the M -point, which resides in the UV-part of the electromag-netic spectrum.We now want to introduce the effect of relaxation intothe calculation of the charge-charge susceptibility. Thiscan be done using Mermin’s approach developed for the3D electron gas . Following Mermin’s work and mak-ing the necessary modifications for the graphene case,the total susceptibility, taking into account relaxation pro-cesses, is given by: χ M ( q , ω ) = (cid:0) i ( τ ω ) − (cid:1) χ ( q , ω + iγ )1 + i ( τ ω ) − χ ( q , ω + iγ ) /χ ( q , ω = 0) , (21) FIG. 4. (Color on-line.) Plot of the electronic density ρ k for different values of the pumping frequency ω p and pumping orientation θ withintensity E = 0 . GV/m, E F = 0 , (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. The bright regions in the Brillouin zone depend on the value of ω p and have the orientation dependence of the term ∇ k Θ k · E , which in its turn depends on the polarization angle, θ , of the incident field.Note that for low frequencies (left panels, ω p /t TB = 0 . ) only momentum values near the Dirac points are excited. On the other hand, for ω p /t TB = 2 the brightest spots occur at the M − point. Also note that the M − points are not all excited at the same time, but depend on thepolarization of the pumping field. This result constrasts with the case where the frequency of the pumping field pumps electrons to the Diraccone (left panels). In this case, all the Dirac points are excited simultaneously. where /τ = γ is the relaxation rate. For calculatingMermin’s susceptibility we need the Lindhard suscepti-bility χ ( q , ω + iγ ) , which needs to be computed for acomplex frequency. In addition we also need the staticsusceptibility χ ( q , ω = 0) . The dielectric function can beobtained from the suceptibility as : ε ( q , ω ) = 1 − v q χ M ( q , ω ) , (22)where v q = ea / (2 ε m q ) is the 2D Fourier transform ofthe Couloumb potential and ε m = ε + ε is the effectivedielectric constant of the environment for a grapheneclad between two media of dielectric constants ε and ε . We recall that the term a appears in v q because thewave number q is measured in units of the inverse lat-tice parameter a − . For consistency with the Mermin’sformula, we take γ p = γ = γ in the forthcoming equa-tions. In the all the figures we have kept γ p (cid:54) = γ , whichis in agreement with Mermin’s equation in the long wave-length limit. IV. LONG WAVELENGTH LIMIT: ANISOTROPICPLASMON DISPERSION RELATION
The calculation of the integral in Eq. (20) needs tobe done for every different frequence ω and wavenum-ber q . However, as shown in Fig. 5, as q decreases the conductivity reaches the long wavelength limit andwe can show that in this regime the susceptibility in Eq.(20) behaves as q . In this limit, the static susceptibil-ity appearing in the denominator of Eq. (21) tends toa constant value when q → , and therefore Eq. (21)becomes: χ M ( q , ω ) = (cid:0) i ( τ ω ) − (cid:1) χ ( q , ω + iγ ) . (23)We now split the susceptibility in the right hand side ofthe Eq. (23) in the same way as we did in Eq. (19) —thatis in an equilibrium and an out-of-equilibrium parts. Theequilibrium component χ ( q , ω ) can be approximated bythe Drude term for (cid:126) ω < E F , where E F is the Fermienergy: χ doped ( q , ω ) = 4 eE F π (cid:126) a q ( ω + iγ ) . (24)For undoped graphene we have E F = 0 and the Drudecontribution vanishes. For the out-of-equilibrium compo-nent, we obtain a similar expression in the long wave-length limit for the pumped susceptibility using Eq. (20)(details of the calculations are given in Appendix E) inthe form: χ intrapump ( q , ω ) = (cid:88) i,j C ij q i q j ( ω + iγ ) . (25)The term in the right hand side of Eq. (25) correspondsto the intraband pumped contribution and can also bewritten as a quadratic dependence on the modulus ofthe wavevector q . This is one of the central results ofthis paper with far reaching implications.Comparing Eq. (25) with the susceptibility of dopedgraphene in the long wavelength limit in Eq. (24), wecan define an effective Fermi energy, that depends onthe polarization angle ϕ of the probe field relative to thegraphene lattice, as: E eff F ( ϕ ) = E F + f + f m cos(2 ϕ + φ ) , (26)where f , f m , and φ depend only on the properties of thepumping field — E pump , θ , and ω p — which are defined inAppendix E. Finally the susceptibility in Eq. (23) can bewritten as: χ M ( ϕ, ω ) = 4 eπ (cid:126) a E eff F ( ϕ ) q ω ( ω + iγ ) . (27)The plasmon dispersion is obtained from the condition ε ( q , ω ) = 0 in Eq. (22), leading to: (cid:126) ω ( ϕ, q ) = (cid:114) α (cid:126) ca E eff F ( ϕ ) q − i γ , (28)where α ≈ / is the fine structure constant of atomicphysics. Equation (28) has the same √ q dependence asthat of plasmons in doped graphene without the pump-ing field . The difference lies in the presence of aneffective Fermi energy E eff F ( ϕ ) that depends on the di-rection of the wavevector. Equation (28) and is one ofthe central results of this work. Note that the disper-sion will be anisotropic, as the effective Fermi energydepends on the orientation of the pumping electric fieldrelatively to the graphene lattice. Furthermore, even inthe case of neutral graphene, the system support plas-mons since E eff F ( ϕ ) is finite even for E F = 0 , due to theconstant illumination of the pumping field. V. THE ANISOTROPIC CONDUCTIVITY OF GRAPHENEUNDER PUMPING
In this section we show that in an out-of-equilibriumsituation we can define an anisotropic optical conductiv-ity for graphene. The optical conductivity tensor σ ij ( q , ω ) can be obtained via the continuity equation: q · J − ωρ = 0 , (29)where ρ is the charge density and J the surface densitycurrent. The current is described by J i = (cid:88) j σ ij ( q , ω ) E j = i (cid:88) j σ ij ( q , ω ) q j Φ . (30)The previous result follows from the relation between theelectric potential Φ , with well defined momentum q , and the electric field E via the relation E = − ∇ Φ = i q Φ . Onthe other hand, the charge density is obtained from thecharge-charge susceptibility via ρ = χ M ( q , ω )Φ . Thus,using Eq. (29), the relation between the conductivity ten-sor and the susceptibility is: (cid:88) i,j σ ij ( q , ω ) q i q j = iωχ M ( q , ω ) . (31)The Equation (31) is not enough to determine the con-ductivity tensor from the susceptibility, but in the longwavelength limit, q → , the dependence of each ele-ment of the conductivity tensor on the wavenumber dis-appear, and we can obtain three independent equationsto the four quantities σ ij . These three equations can beobtained changing the direction of the wavevector q or,equivalently, we can compare the Taylor expansion of χ M ( q , ω ) to the left hand side of Eq. (31). This proce-dure would give four equations but one of them wouldnot be independent of the other three. The missingequation can be obtained from the current-current re-sponse, calculated in appendix F, where the intrabandcontributions to the conductivity tensor read σ intra ij ( q , ω ) = 2 ie (cid:126) ωS (cid:88) k ,λ = ± n λ k + q / − n λ k − q / ω − λω intra k , q + iγ v i intra k , q v j intra k , − q , (32)with v i intra k , q defined in Appendix F and S = N c a , where N c is the number of unit cells. In this appendix an ex-pression for the interband term is also provided. FromEq. (32) we can show that in the limit q → , we have: σ intra ij ( q → , ω ) = σ intra ji ( q → , ω ) , (33)thus it follows from Eq. (31) that: σ pumped ij = σ C ij ω + iγ , (34)where σ = e / (4 (cid:126) ) and C ij are the coefficients of theexpansion of χ pumped ( q , ω ) defined in Eq. (25), and thetotal intraband conductivity can be written as function ofan effective Fermi energy tensor E eff ij as: σ ij σ = i (cid:126) π E eff ij ω + iγ , (35)where we have defined the effective Fermi energy tensoras: E eff ij = E F δ ij + π C ij . (36)Although the tensor E eff ij can be reduced to diagonalform by a rotation, doing so we loose the direct connec-tion of the tensor components to the orientation of thegraphene lattice.We can also define the longitudinal conductivity alongthe direction defined by the unit vector u ϕ for a probingelectric field of the form E = E u ϕ as: σ ϕ = J · u ϕ E = 4 i (cid:126) π E eff F ( ϕ ) ω + iγ , (37)where the angle ϕ (the polarization angle of the probingfield) is the same as that of the momentum q , since theelectric field is proportional to q via the gradient of thepotential.It is worth remembering that the effective parameters f , f m , and φ depend solely on the pumping field proper-ties, that is, on the intensity E , the polarization angle θ ,and the frequency ω p . VI. NUMERICAL RESULTS
As shown before (see Fig. 4), graphene under in-tense and energetic light pumping presents a stronganisotropic electronic distribution. This changes theoptical response due to intraband and interband tran-sitions. In doped graphene, without electromagneticpumping, the intraband transitions dominates for pho-ton energy (cid:126) ω < E F , while interband transitions dom-inate for (cid:126) ω > E F . For pumped graphene, we havea similar result, where intraband transitions dominate for ω < ω p , where ω p is the frequency of the pumping radi-ation, and interband transitions dominates for ω ≈ ω p .To show the effects of the pumping in graphene, wesolve numerically the Eq. (20) and compute the pumpedcomponent of the intra-band susceptibility. The χ ( q , ω ) component is calculated with the analytical expressionsderived with the Dirac equation , since it is not nec-essary here to account for the full band structure ofgraphene. This is because for the equilibrium distribu-tion, the probe frequency ω in the range we are con-sidering can only excite electron-hole pairs around theDirac cone. This is not the case for the pumped electrongas around the M -point, that cannot be described bythe Dirac equation. The out-of-equilibrium distribution iscalculated using Eq. (16).We plot in Fig. 5 the imaginary part of the longitudinalconductivity, Eq. (37), as function of the probe incidenceangle ϕ , for different wave numbers q ; this quantity con-trols the dispersion of the surface plasmon-polariton inthe out-of-equilibrium electron gas, as will be discussedin a forthcoming section. We show that the long wave-length limit is reached around q = 10 − , where the con-ductivity have a co-sinusoidal shape as predicted by Eq.(26). Note that q is measured in units of /a . In thesame figure we also depict an example of the longitu-dinal conductivity away from the long wavelength limit( q = 10 − ). It is clear that in this regime the distributionis, for some ϕ values, substantially different form the an-alytical approximation. Since we are interested here inthe long wavelength limit, this results does not interestus and will not be discussed further.Figure 6 shows the numerically computed imaginarypart of the longitudinal conductivity, as function of theprobing polarization angle ϕ , for q = 10 − , as definedby Eq. (37), compared with the semi-analytical result,which depends on the effective Fermi energy definedin Eq. (26). The two approaches show a very good FIG. 5. (Color on-line.) Imaginary part of the longitudinal conduc-tivity as function of the polarization angle of the pump field for differ-ent values of q . The parameters are: E = 0 . GV/m, (cid:126) ω p = 2 t TB , θ = π/ , E F = 0 . eV and (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. For q < − the susceptibility reaches the long wavelength limit. For q = 10 − we can see the strong influence of the static susceptibility(see Fig. 14). agreement, showing that indeed for q = 10 − the sys-tem is already in the long wavelength regime. The os-cillatory variation of the imaginary part of the longitudi-nal conductivity will lead to an anisotropy in the spec-trum of the surface plasmon-polariton, as it is this quan-tity that determines the behavior of the latter. Note that (cid:61) σ ∈ [10 σ , σ ] (see Fig. 6).From here on our analysis is focused on two ways ofparameterizing the effective Fermi energy. In the firstapproach, we discuss Eq. (26), which is suitable for an-alyzing the susceptibility (24) and plasmons modes (28)at long wavelengths. In the second approach, the resultof Eq. (35) is useful to calculate the optical conductiv-ity at long wavelengths and the dispersion of the surfaceplasmon-polariton, defined by Eq. (40).Figure 7 shows that the parameters f and f m have astrong dependency on the intensity of the pumping field.The angle φ , in contrast, changes very little by as muchas ∼ . rad, and tends to saturate for large intensityfields. The parameters f and f m can have a strong im-pact in the optical response of the system, dependingon the initial doping level of graphene, characterized by E F . For large doping, the effect of f and f m is small,except for large pumping field intensities. However, forvanishing small Fermi energies, the effect of these twoparameters have a large impact in the optical propertiesof the system, as the effective Fermi energy is essen-tially controlled by them.The dependence of the effective Fermi energy E Fij on FIG. 6. (Color on-line.) Comparision between the semi-analtical ap-proach and the numerical one for the imaginary part of the longitudi-nal conductivity, showing the validity of the semi-analytical approxi-mation obtained in section IV. The dots correspond to the black solidcurve in Fig. 5 and the solid line is the semi-analytical calculation.Note that (cid:61) σ ∈ [10 σ , σ ] ; the difference between the maximumand the minimum of the conductivity depends on the magnitude of E . The parameters are: q = 10 − , E = 0 . GV/m, (cid:126) ω p = 2 t TB , θ = π/ , E F = 0 . eV, and (cid:126) γ = 14 meV, and γ p = 28 meV. the intensity of the pumping radiation is depicted in Fig.8. A clear anisotropy is seen in this quantity. Particu-larly interesting is the finite value of E Fxy , which leads toa finite off-diagonal term for the non-equilibrium opticalconductivity. The dependence of the parameters f , f m ,and φ on the energy of the pumping photons is depictedin Fig. 9. We see that there is a non-monotonous de-pendence on ω p with a local maximum (for f , f m , and φ )when the photon energy is equal to the electronic tran-sition at the M -point ( (cid:126) ω p /t TB = 2 ). This is, most likely,due to the enhanced density of states associated withthe van-Hove singularity. In Fig. 10 the effective Fermienergy is depicted as function of the frequency of thepumping field. Clearly its behavior is controlled by thevalues of the parameters f , f m , and φ , as can be seenfrom comparing Figs. 9 and 10. Again a local maximumis seen at the value of photon energy given by (cid:126) ω = 2 t TB .It is worthwhile to remark that the absolute maximumof the f and f m parameters, in Fig. 9, takes place for (cid:126) ω p ≈ . t TB ( ∼ (cid:126) ω p ∼ t TB .This energy scale is controlled by the electric field inten-sity. Indeed, the system has an energy scale ∆ , for theparameters of Fig. 9, given by ∆ /t TB ∼ (cid:114) E a t TB ∼ . , (38) FIG. 7. (Color on-line.) Dependence of the parameters f , f m , and φ on the intensity of the pumping radiation. The parameters are: θ = π/ , (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. The importance of the parameters f snd f m grows with theintensity of the pumping field. The minimum value of field intensityconsidered in this figure is 0.1 GV/m. We see that the anisotropy isobservable for this field intensity. We note that the field intensitiesscanned in this figure are experimentally attainable. which is of the same order of magnitude of (cid:126) ω p ≈ . t TB ,the position of the absolute maximum of the parame-ters f and f m . Note that apart from the gradient of thephase Θ k , E a is essentially the Rabi frequency. Wehave verified that by reducing the field intensity by fivetimes, the position of the maximum red-shifts to an en-ergy of about two times smaller the value of (cid:126) ω = 0 . t TB .This effect is represented in the bottom panel of Fig. 10.The scaling of the position of the maximum of f with √E is evident. Note, however, that, for these energyscales, the anisotropy for the plasmon spectrum will bevery small, as E Fxx ≈ E Fyy . Let us also note here thatthe intensity of the density of the states at the van-Hovesingularity is presumably controlled by the value of γ :the larger this parameter is the smaller is the density ofstates at the M − point, which otherwise would be a di-vergence in the absence of relaxation.In Figs. 11 and 12 we show the strong anisotropyin the optical response. The parameter f m , that mea-sures the amplitude of the effective Fermi energy mod-ulation, has maxima where the f presents minima for0 FIG. 8. (Color on-line.) Dependence of the effective Fermi energy E eff ij − E F δ ij on the intensity of the pumping radiation (in GV/m).The anisotropy grows with the increase of E . The parameters are: θ = π/ , (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. The minimum value of field intensity considered in this figureis 0.1 GV/m. some specific angles. This is the origin of the stronganisotropy in the optical response of the system, whichimparts in the anisotropy of the dispersion relation of theplasmons. In the Fig. 11 the parameter φ is also de-picted showing a strong variation with the angle of po-larization of the pumping field. The strong variation of f , f m , and φ on θ controls the dispersion of the plas-mon in this system.We emphasize that the results presented in Figs. 7-12correspond to the contribution from intraband transitionsthat take place near the three independent M -points (inthis case the concept of valley is meaningless). Notethat Fig. 4 shows the effect of the anisotropic electronicdistribution near each M -point and the different occupa-tions of each M -point. For this electronic distribution,the parity symmetry is broken (see Fig. 4), and, asa consequence, we can have a finite off-diagonal con-ductivity. The same symmetry is broken in the in theHamiltonian studied by Kumar et al. . However, in thiscase the parity symmetry is broken by a circular polar-ized pumping field that populates each valley differently(in graphene each valley is connected by the parity sym-metry).One experimental way of accessing the dispersion ofthe plasmons in a given material is to perform a EELSexperiment. This spectroscopic technique is based onthe excitation of plasmons by moving charges. Whenexciting a plasmon wave, the incoming electrons losepart of their kinetic energy. Theoretically, the loss func-tion, which encodes the excitation of the plasmons bythe moving electrons, is defined in terms of the dielec- FIG. 9. (Color on-line.) Dependence of the parameters f , f m , and φ on the energy of the photon of the pumping field. The parametersare: E = 0 . GV/m, θ = π/ , (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. A non-monotonous dependenceon ω p is seen for the three parameters. We must however stress thatfor (cid:126) ω p ≈ the behavior of f , f m , and φ is not accurate, as wehave not included the effect of interband transitions, due to probe offrequency ω , which become relevant for ω p ∼ ω , specially in thecase of neutral graphene. See Fig. 10 for a discussion of the positionof the maximum of f and f m located at low energies. tric function as: L ( q , ω ) = −(cid:61) (cid:26) ε ( q , ω ) (cid:27) . (39)This quantity is depicted in Fig. 13, for different values ofthe probing polarization angle ϕ . The dielectric functionwas calculated using Eq. (22) and the pumping suscep-tibility is given by Eq. (20). In Fig. 13 we can see thecharacteristic plasmon signature in the loss function. Itis clear that the plasmon spectrum depends significantlyon the polarization of the probing field, or, in other terms,on the direction of the momentum in the Brillouin zone.The width of the plasmon spectrum is proportional to therelaxation rate γ . For making apparent the anisotropywe also depict (solid line) the dispersion of the plasmon1 FIG. 10. (Color on-line.) Top panel: Dependence of the componentsof the effective Fermi energy tensor E eff ij − E F δ ij on the pumpingfield frequency. The parameters are: E = 0 . GV/m, θ = π/ , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. Note the localmaximum of the Fermi-energy tensor-elements around the photonenergy (cid:126) ω = 2 t TB . Also note that the largest difference between E Fxx and E Fyy occurs at the M − point which implies the largest anisotropyin the properties of the system, including the plasmon spectrum. Wemust stress that for (cid:126) ω p ≈ the behavior the effective Fermi energycomponents are not accurate, as we have not included the effect ofinterband transitions, due to the probe of frequency ω , which becomerelevant for ω p ∼ ω , specially in the case of neutral graphene. Cen-tral panel: Zoom in of the dependence of the parameter f with E near the absolute maximum. Bottom panel: Scaling of the frequencyof the maximum, ω m , with the √E (right panel); the linear scalingis evident. The values of ω m are extracted from the central panel, andcorrespond to the position of the maximum of the curves for f . Notethat the larger E is the broader is the maximum and more intense is f . after an average of the effective Fermi energy on the po-larization angle ϕ ; the anisotropy is obvious.Let us now discuss the reason why the plasmon char-acterizing pumped graphene out-of-equilibrium is similarto that of doped graphene in equilibrium, in what con-cerns their small energy values. In the latter case, for (cid:126) ω + (cid:126) v F q < E F and ω > v F q , where q is the wavenum- FIG. 11. (Color on-line.) Dependence of the parameters determiningthe effective Fermi energy on the pumping polarization angle. Notethat for some values of θ the magnitudes of f and f m are almostidentical. Also the angle φ varies substantially with θ . The largestanisotropy in the properties of the system occurs for the largest dif-ference between f and f m . The parameters are: E = 0 . GV/m, (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. ber and ω the frequency, interband process are sup-pressed by Pauli-blocking and the susceptibility is dom-inated by intraband processes, where losses are pro-portional to the relaxation rate γ (for γ = 0 the usualplasmons are infinitely long-lived in this momentum-frequency window). In this regime, graphene supportsplasmons with small attenuation with a relation disper-sion proportional to √ q . On the other hand, when weconsider the case of the pumped distribution, the situ-ation is similar, because interband process, that atten-uates the plasmon, only occur for frequencies ω nearthe pumped frequency ω p . Since we are considering theregime ω (cid:28) ω p , the attenuation of the plasmons of thenon-equilibrium electron gas is essentially controlled bythe value of γ (the plasmons cannot decay via particle-hole processes in this regime, as it happens in the caseof an equilibrium plasma). Therefore, the correspon-dent pumped susceptibility is similar in the sense thatthe imaginary part is proportional to the scattering time[see Eq. (24)]. Thus, we can expect for plasmons in theout-of-equilibrium electron gas the same level of attenu-2 FIG. 12. (Color on-line.) Dependence of the effective Fermi-energy E eff ij − E F δ ij on the polarization angle of the pumping field. Weemphasize that our calculations take the three M − points into ac-count simultaneously since we are making a tight-binding calcula-tion. Therefore there is no cancellation of E xy . The parameters are: E = 0 . GV/m, (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. Note the periodic behavior of the different parame-ters. ation of the conventional plasmons in graphene. As con-sequence the former anisotropic plasmons are expectedto be long lived as are their siblings in the equilibriumelectron gas.In Fig. 14 we show that the graphene staticsusceptibility have zeros that renders the term ( τ ω ) − χ ( q, ω ) /χ ( q ) ( τ = 1 /γ ) in Mermin’s susceptibilitylarge, even at small q . In this case the use of the Mer-min’s equation is no longer valid, since the assumptionthat the fluctuations of the local Fermi Energy, which areproportional to /χ ( q ) , are small is no longer true andthe approximation leading to Mermin’s equation breaksdown. VII. SPECTRUM OF THE SURFACEPLASMON-POLARITONS IN THE OUT-OF-EQUILIBRIUMREGIME
As a conductive two dimensional system, graphenesupports surface plasmon-polaritons. We now want toaddress the propagation of the these quasi-particles onthe surface of graphene due to the electron gas cre-ated by the pumping field. We will see that the surfaceplasmon-polariton spectrum in pumped graphene showsa dispersion strongly dependent on the ϕ angle, the po-larization angle of the probing field. A surface plasmon-polariton (SPP) is an hybrid particle that couples elec-tromagnetic radiation to the free oscillations of an elec-tron gas in a conductor. In graphene, the spectrum ofan SPP depends critically on the nature of the optical FIG. 13. (Color on-line.) Loss function for different polarizations ofthe probe field as function of the dimensionless wavenumber (mul-tiplied by ). The parameters are E = 0 . GV/m, θ = π/ , (cid:126) ω p = 2 t TB , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. Thesolid (cyan) curve is the plasmon dispersion for the semi-analyticalresult in Eq. (28) after an average of the effective Fermi energy onthe polarization angle ϕ conductivity of the system (for a discussion about sur-face plasmon-polariton in graphene see Refs. [29 and30]). Indeed, it can be shown that the condition for theexistence of an SPP is given by (cid:20) ε k + ε k + iσ xx ωε (cid:21) (cid:20) k + k ωµ − iσ yy (cid:21) − σ xy σ yx ωε = 0 , (40)for a wave propagating along the x direction and de-caying exponentially along the direction perpendicular tothe graphene plane. When σ xy = σ yx = 0 , the trans-verse electric and the transverse magnetic modes de-couple. In the case we are considering here this is notthe case, since the non-equilibrium nature of the elec-tron gas created by the pumping induces a finite valuefor σ xy . However, since time reversal symmetry is notexplicitly broken in this case, we have the condition that σ xy = σ yx . Using the calculated conductivity tensor inEq. (35) and the coefficients C ij calculated through Eq.(E3), the spectrum of the SPP can be obtained.3 FIG. 14. (Color on-line.) Static susceptibility for q = 2 . − , E =0 . GV/m, and θ = π/ as function of the polarization angle ofthe probe field. Note the existence of points where the susceptibilityis zero. Near and and at these points Mermin’s approach breaksdown.FIG. 15. (Color on-line.) Plasmon dispersion relation (top panel)and plasmon lifetime (bottom panel) as function of the dimensionlesswave number (multiplied by ), for E = 0 . GV/m, (cid:126) ω p =2 t TB , θ = π/ , for two probing polarization angles. Note that γ pl (cid:28) ω pl . The dispersion relation of the surface plasmon-polariton due to the non-equilibrium electron gas de-pends on the orientation of the direction of propaga-tion of the wave with respect to the crystalline lattice.To describe the propagation along another direction, wecan rewrite Eq. (40) in the new reference frame or,alternatively, rotate the conductivity tensor. The latter
FIG. 16. (Color on-line.) Surface plasmon-polaritons dispersion re-lation, of the out-of-equilibrium electron gas, for different values ofthe angle ϕ of the polarization of the probing radiation (shadded ra-gion) as function of the dimensionless wave vector (multiplied by ). Note that all the angles in the interval ϕ ∈ [0 , π ] are con-tained in the shaded region. The boundary values are determined bythe co-sinusoidal form of the effective Fermi energy. Note that thevariation of ω pl as ϕ varies is substantial. Therefore the anisotropy isamenable of experimental verification. The parameters are: E = 0 . GV/m, (cid:126) ω p = 2 t TB , θ = π/ , E F = 0 . eV, (cid:126) γ = 14 meV, and (cid:126) γ p = 28 meV. can be achieved with the usual 2D rotation matrix M ϕ , σ (cid:48) = M ϕ σM − ϕ : M ϕ = (cid:18) cos ϕ − sin ϕ sin ϕ cos ϕ (cid:19) . (41)In Fig. 16 we show, for fixed E , ω p , E F , and θ , thesurface-plasmon polariton in graphene from the solutionof Eq. (40). The shaded region corresponds to differentvalues of the variable ϕ , between those represented bythe black solid lines at the borders of the shaded region.We see again the strong dependence of the optical prop-erties upon the probe angle ϕ , which is measured by theanisotropy in the SPP spectrum. Note that the variationof the spectrum with ϕ is quite substantial and thereforeamenable to experimental verification. VIII. FINAL COMMENTS
In this work we have considered a pump-probe prob-lem, where the pumping field is a relatively intense and4pulsed wave field, with a pulse duration much larger than1 ps. In this situation we can reach a stationary statewhere an out-of-equilibrium electron gas is maintainedin the conduction band in graphene. We have con-sidered the case where the frequency of the pumpingfield lies in the UV-range of the electromagnetic spec-trum. In this case the electrons are pumped to the M − point in the Brillouin zone. In addition to the pump-ing, a probe field of much smaller frequency probes theout-of-equilibrium electron gas. This allows us to ac-cess the collective plasma wave —plasmons— in theout-of-equilibrium electron gas. We have shown that forpumping field of this frequency the excitation of the three M − points in the Brillouin zone is uneven, at odds withthe excitation of an electron gas near the Dirac points.This is a consequence of the strong deviation of theband structure of graphene from the Dirac cone approx-imation. Indeed, near the M − point the band structurehas a saddle point nature being, therefore, very differ-ent from the Dirac cone. Interestingly enough, we havefound that the plasmon in the out-of-equilibrium elec-tron gas still scales with the √ q as in the case of theDirac plasmons. This is a consequence of the formof the charge-charge susceptibility, which scale as q i q j ( i = x, y ) in the long wavelength limit (note that in theDirac cone approximation the charge-charge suscepti-bility scales as q ). This scaling can still be writtenin terms of q if we introduce an effective Fermi en-ergy, depending on the properties of the pumping field.The anisotropy of the plasmon dispersion in the Brillouinzone originates from the scaling q i q j and is encondedin the effective Fermi energy. At the more fundamen-tal level, the fact that the out-of-equilibrium susceptibilityscales with q i q j in the long wavelength limit is a con-sequence of the continuity equation (31) that links thesusceptibility with the conductivity. If in the long wave-length limit the susceptibility scales with a power lowerthan q the conductivity would diverge and if the poweris greater than q , the conductivity would be null.Due to the relation between the charge-charge sus-ceptibility and the optical conductivity, it is possible todefine an out-of-equilibrium optical conductivity. In-terestingly, the non-linear dependence of the out-of-equilibirum distribution function on the pumping field al-lows for a finite value of σ intra xy = σ intra yx (cid:54) = 0 (for theinterband conductivity the situation is identical σ inter xy = σ inter yx (cid:54) = 0 . Onsagar relation requires σ xy ( H ) = − σ yx ( H ) in the presence of magnetic field H , or similarly brokentime reversal symmetry. A Hamiltonian with circular po-larized external light field is not time invariant, in whichcase we would have a different result from above). Thishas an impact on the spectrum of the surface plasmon-polaritons (SPPs) that can be supported by the out-of-equilibirum electron gas, as in this case, the TE and TMpolarization are coupled to each other. We have foundthat the measured values for SPP spectrum depend onthe orientation of the polarization of the probing field.This is a consequence of the anisotropy of the optical conductivity of graphene in the regime considered.What is missing form this work is a detailed study ofthe effect of electron-phonon and electron-electron in-teractions, which has been included only at the level ofa phenomenological scattering rate. Therefore phenom-ena such as carrier multiplication is not included in ourdescription. It would be an interesting to discuss thisproblem in the regime we have considered, a problemthat was not analysed in the literature before, but this isoutside the scope of this paper.Let us comment briefly on the nature of the lightsource needed to deliver the required electric field in-tensity to observe the anisotropy. As noted in the in-troduction, due to excitonic effects the position of theabsorption maximum at the M − point (due to inter-bandtransitions) is shifted from 5.4 eV to about ∼ λ ∼ nm), a wavelength for which there are avail-able lasers. The required wavelength can be obtainedfrom the fourth harmonic (cascade two harmonic gener-ation) of a Q-switch diode-pump solid-state laser. Thisis a very popular wavelength, 266 nm, which corre-sponds to a transition of ∼ P ∼ kW, and a FWHM for the minimum beamwaist of w = 600 µ m (these are figures of commerciallyavailable lasers) it follows that the intensity of the electricfield is about E ∼ (cid:114) P log 2 πw c(cid:15) ∼ .
008 GV / m , (42)which is not yet enough for observing the effects wehave discussed above for weakly doped graphene,which become apparent for E ≥ . GV/m. Howevera lens can be used to increase the value of E (seeahead). We need to comment here on the importanceof the magnitude of γ and γ p on the value of E neededfor observing the effects we have addressed in this pa-per. In this work we have assumed that (cid:126) γ p = 28 meV(larger than (cid:126) γ ) in connection with the characteristictime τ em = 0 . ps implied by luminescence experiments;if we had chosen the value of τ em = 0 . ps, correspond-ing to the lowest figure suggested by experiments, thevalue of E necessary for observing the anisotropy wouldhave decreased by about √ . We have also verifiedthat in the opposite regime γ (cid:29) γ p the value of E in-creases significantly over the value given by Eq. (42)and the observation of the anisotropy reported in thispaper may be out of experimental reach.Let us note here that the anisotropy is not a prop-erty of the M − point alone, so we can choose a largerwavelength than that necessary to excite electrons atthe M − point but still small enough for the Dirac-coneapproximation to hold. In this condition, similar effectsto those described in this paper for the M − point wouldstill be visible. Note that our approach takes the full bandstructure into account and therefore is valid for all ener-gies as long as ω p (cid:29) ω ; the choice for the M − point wasonly a matter of selecting a high symmetry point in the5Brillouin zone.Another important experimental constrain is thethreshold power per unit area above which graphene isdamaged. It was found that for a laser of λ = 248 nmthe threshold power per unit area was of the order of . × W/m (The experiment considered 500 laserpulses of 20 ns duration and a repetition rate of 100 Hz;see however Ref. [32]). This value is an order of mag-nitude smaller than P/w ∼ . × W/m estimatedfrom the numbers given above. This implies a reduc-tion in the power P and therefore a smaller value of E .Indeed, using P/w ∼ . × W/m it follows E ∼ .
004 GV / m , (43)which is only a factor of 2 smaller than that found in Eq.(42). On the other hand, a single pulse of 20 ns durationis much larger than 1 ps needed for attaining the steadystate (see Fig. 3). Therefore, within the duration of asingle 20 ns pulse it is possible to excite the plasmon inthe out-of-equilibrium gas and measure their existence.The use of a single pulse allows to multiply the powervalue of . × W/m by 500 (number of repetitionsneeded to observe clear changes in the D-peak of theRaman spectrum of graphene) thus allowing a value of E , for a single pulse, of the order of E ∼ . / m , (44)and introducing only a small amount of laser ablationon graphene. Therefore, the value used in the paperof E ∼ . / m is acceptable for illustrative purposes,since it is about 7 times higher than the estimation madeabove. Therefore, the used value for E has been cho-sen for making the effects more apparent to the nakedeye. Naturally, focusing the laser spot onto an areasmaller than w , makes the figure of . / m accept-able. Indeed, for the laser indicated above with an exitspot waist of w = 800 µ m we can focus it down toa spot waist of µ m which implies that the value of E ∼ . / m in well within the experimental range. Inconclusion, it is expected that by lowering the repetitionrate, the damage threshold would increase beyond the . × W/m reported in Ref. [31]. This is justified bythe fact that the average power incident on the samplewill decreased and thus, also the damage due to ther-mal effects. Furthermore, focusing the laser beam leadsto values of E more than one order of magnitude largerthat given by Eq. (42), implying that the effects predictedin this paper become experimentally accessible.Let us finally comment on a difference between ourresults and those of Ref. [10]: the Rabi frequency inEq. (11) depends on the gauge choice for the electro-magnetic field. When we add the relaxation time webreak gauge invariance. In the work of Singh et al. the equations of motion were obtained with the minimalcoupling in the Couloumb gauge. The correspondingRabi frequency has a factor ω k /ω in comparision withEq. (11). However, when the relaxation time goes to infinity, γ → , the two approachs give identical results.Indeed, from the equations in Appendix F, that were writ-ten in the same gauge as that used in Ref. [10], we canmake q = 0 (the wave number of the probing field) andproceed in the same way as we did in Sec. II to obtainthe same results as those found in Ref. [10]. ACKNOWLEDGMENTS
A. J. Chaves acknowledges the scholarship fromthe Brazilian agency CNPq (Conselho Nacional deDesenvolvimento Científico e Tecnológico). N.M.R.Peres acknowledges support from the European Com-mission through the project “Graphene-Driven Revolu-tions in ICT and Beyond" (Ref. No. 696656) andthe Portuguese Foundation for Science and Technol-ogy (FCT) in the framework of the Strategic FinancingUID/FIS/04650/2013. The authors acknowledge JoséCarlos Viana Gomes for discussions that led to the esti-mations made in Sec. VIII.
Appendix A: Tight-binding model for graphene subjectedto an external electric field
The graphene tight-binding Hamiltonian in secondquantization reads: H = (cid:88) i,n t TB ˆ a † R n ˆ b R n + δ i + h.c. , (A1)where ˆ a † R n and ˆ b R n + δ i obey anti-commutation relationsand δ i are the nearest neighbors vectors connecting anatom in sub-lattice A to another one in sub-lattice B [seeEq. (D2)]. We can define the Fourier transform and itsinverse as: ˆ a k = 1 √ N c (cid:88) n e − i k · R n ˆ a R n , (A2a) ˆ b k = 1 √ N c (cid:88) n e − i k · ( R n + δ i ) ˆ b R n + δ i , (A2b) ˆ a R n = 1 √ N c (cid:88) k ∈ o B.Z. e i k · R n ˆ a k , (A3a) ˆ b R n + δ i = 1 √ N c (cid:88) k ∈ o B.Z. e i k · ( R n + δ i ) ˆ b k , (A3b)where the sum over n is performed over the entire lat-tice and the sum in k is performed over the first Brillouinzone.After a Bogoliubov transformation the basis that diag-onalize H is: ˆ c k = e iϕ k √ (cid:16) ˆ a k + e i Θ k ˆ b k (cid:17) , (A4a) ˆ d k = e iϕ k √ (cid:16) ˆ a k − e i Θ k ˆ b k (cid:17) , (A4b)6 FIG. 17. (Color on-line.) Electronic band structure of graphenefor π − electrons. The zoom-in shows the band structure around the M − point in the Brillouin zone. The saddle-like nature of the bandstructure is clearly visible around that point. and the inverse transformation reads: ˆ a k = e − iϕ k √ (cid:16) ˆ c k + ˆ d k (cid:17) , (A5a) ˆ b k = e − iϕ k e − i Θ k √ (cid:16) ˆ c k − ˆ d k (cid:17) , (A5b)where ϕ k is a global arbitrary phase. The phase Θ k isthe argument of φ k = (cid:88) i =1 e i k · δ i , (A6)that is, Θ k = arg φ k . In the basis (A4) the Hamiltonian H is written as: H = (cid:88) k E k (cid:16) ˆ c † k ˆ c k − ˆ d † k ˆ d k (cid:17) . (A7)The band structure given by ± E k is depicted in Fig. 17together with a zoom-in around the M − point. For writingthe the interaction term with the electric field we need theposition operator written as: ˆR A = (cid:88) n R n ˆ a † R n ˆ a R n , (A8a) ˆR B = (cid:88) n ( R n + δ ) ˆ b † R n + δ ˆ b R n + δ , (A8b) which in the basis given by Eq. (A4) it reads: R = − (2 π ) i N c a (cid:88) k , q [ ∇ q δ ( q )] e i ( ϕ k + q − ϕ k ) (cid:104)(cid:16) e i (Θ k + q − Θ k ) (cid:17) (cid:16) ˆ c † k + q ˆ c k + ˆ d † k + q ˆ d k (cid:17) ++ (cid:16) − e i (Θ k + q − Θ k ) (cid:17) (cid:16) ˆ c † k + q ˆ d k + ˆ d † k + q ˆ c k (cid:17)(cid:105) , (A9)where δ ( q ) is the Dirac delta-function of zero momen-tum. Therefore, the radiation-electron interaction is fi-nally written as: H I = e E · (cid:88) k (cid:20) i ∇ k (ˆ n c, k + ˆ n v, k ) + ∇ k Θ k p cv, k + ˆ p vc, k ) (cid:21) , (A10)where we have defined: ˆ n c, k = c † k c k , (A11a) ˆ n v, k = d † k d k , (A11b) ˆ p cv, k = c † k d k , (A11c) ˆ p vc, k = d † k c k , (A11d) ϕ k + q = − Θ k + q . (A11e) Appendix B: Derivation of the Bloch equations forgraphene
To obtain the Bloch equations in graphene we calcu-late the commutator (5) with the introduction of two phe-nomenological damping terms and with the density ma-trix written in the basis (A4). After the calculation of theexpectation values we obtain the set of equations: − ∂ t n c, k = γ ( n c, k − f c, k ) + i Ω k ( t )∆ p k , (B1a) − ∂ t n v, k = γ ( n v, k − f v, k ) − i Ω k ( t )∆ p k , (B1b) ( ∂ t + iω k + γ p ) p cv, k = − i Ω k ( t )∆ n k , (B1c) ( ∂ t − iω k + γ p ) p vc, k = i Ω k ( t )∆ n k , (B1d)where (cid:126) ω k = 2 E k , ∆ n k = n c, k − n v, k , ∆ p k = p cv, k − p vc, k , f c/v, k is the Fermi-Distribution for the conduc-tion/valence band, and γ ( γ p ) is a relaxation term. Thetime dependence on n c/v, k , p vc/cv, k , and E is omitted,and we have defined: Ω k ( t ) = ea E ( t ) · ∇ k Θ k (cid:126) . (B2)Summing Eqs. (B1a) and (B1b) we obtain: ∂ t ( n c, k + n v, k ) = − γ ( n c, k + n v, k − ( f c, k + f v, k )) , (B3)which has the exact solution: n c, k ( t ) + n v, k ( t ) = c ( k ) e − γ t + f c, k + f v, k , (B4)where c ( k ) depends on the initial conditions.7For a system that is initially in thermal equilibrium, c ( k ) = 0 and: n c, k ( t ) + n v, k ( t ) = f c, k + f v, k , (B5)thus we introduce the deviation ρ k ( t ) through: n c, k ( t ) = f c, k + ρ k ( t ) , (B6a) n v, k ( t ) = f v, k − ρ k ( t ) . (B6b)We also note that the complex conjugate of (B1d) reads: (cid:18) ∂∂t + iω k + γ p (cid:19) p ∗ vc, k = − i e E · ∇ k Θ k n k , (B7)and using Eq. (B1c) (cid:18) ∂∂t + iω k + γ p (cid:19) (cid:0) p cv, k ( t ) − p ∗ vc, k ( t ) (cid:1) = 0 , (B8)we find the solution: p cv, k ( t ) = p ∗ vc, k ( t ) + c ( k ) e ( − iω k − γ p ) t , (B9)where again c ( k ) depends on the initial conditions. Ifwe assume that the system is initially in thermal equilib-rium it follows that c ( k ) = 0 and: p vc, k ( t ) = p ∗ cv, k ( t ) = x k ( t ) + iy k ( t ) . (B10)The set of four complex equations (10) can be re-duced to a set of three real equations for the functions x k , y k , and ρ k . From Eqs. 10 we have: ∂ t ρ k = − γ ρ k + 4 v k ( t ) y k , (B11a) ( ∂ t − iω k + γ p ) [ x k + iy k ] = i Ω k ( t ) (cid:0) ∆ n k + 2 ρ k (cid:1) , (B11b)from which finally follows that: ˙ x k = − γ p x k − ω k y k , (B12a) ˙ y k = ω k x k − y k γ p − Ω k ( t ) (cid:0) ρ k + ∆ n k (cid:1) , (B12b) ˙ ρ k = − γ ρ k + 2Ω k ( t ) y k . (B12c)These latter set of equations is the one we have solvedin the bulk of the paper. Appendix C: Steady-state equations for the distributionfunctions under continuous pumping
With the assumption that ∂ t n c k = ∂ t n v k = 0 and con-sidering a monochromatic incident field with frequency ω p , we can write Eqs. (10) as: γ (cid:0) n c, k − n c, k (cid:1) + i (cid:104) Ω k ( t )∆ p k (cid:105) t = 0 , (C1a) γ (cid:0) n v, k − n v, k (cid:1) − i (cid:104) Ω k ( t )∆ p k (cid:105) t = 0 , (C1b) ( ∂ t + iω k + γ p ) p cv, k = − i Ω k ( t )∆ n k , (C1c) ( ∂ t − iω k + γ p ) p vc, k = i Ω k ( t )∆ n k , (C1d) where we use (cid:104)(cid:105) t for time average. The solution to thisset of equations is of the form: p cv, k , ( t ) = A ( ω p ) e iω p t + B ( ω p ) e − iω p t , (C2a) p vc, k ( t ) = A ( ω p ) e iω p t + B ( ω p ) e − iω p t , (C2b)where A i , B i can be obtained from Eqs. (C1c) and (C1d)as: A ( ω p ) = n c, k − n v, k ω p + ω k − iγ p ¯Ω k , (C3a) B ( ω p ) = n c, k − n v, k − ω p + ω k − iγ p ¯Ω ∗ k , (C3b) A ( ω p ) = n v, k − n c, k ω p − ω k − iγ p ¯Ω k , (C3c) B ( ω p ) = n v, k − n c, k − ω p − ω k − iγ p ¯Ω ∗ k , (C3d)with: ¯Ω k = ea E · ∇ k Θ k (cid:126) . (C4)If we define: α k = τ τ p | ¯Ω k | τ p (cid:0) ω k + ω p (cid:1) τ p ( ω p − ω k ) + 2 τ p ( ω p + ω k ) + 1 , (C5)with τ = 1 /γ , τ p = 1 /γ p , we have from Eqs. (C1a) and(C1b) that: n c, k − f c, k = α k ( n v, k − n c, k ) , (C6a) n v, k − f v, k = α k ( n c, k − n v, k ) . (C6b)Note that expression for α k is well defined even takingthe collisionless regime ( τ , τ p ) → ∞ . Appendix D: Expressions for the gradient of the phase Θ k In this appendix we present some useful functions thatappear in the rest of the paper. First we recall the func-tion defined in Eq. (4): φ k = (cid:88) j =1 e i k · δ j , (D1)where we have used the following choice of vectors forthe orientation of the nearest neighbor hopping: δ = a (1 , , (D2a) δ = a / (cid:16) − , −√ (cid:17) , (D2b) δ = a / (cid:16) − , √ (cid:17) . (D2c)8With this choice of vectors, the eigenvalues of H arethe solution of: (cid:18) E k t T B (cid:19) = 1+4 cos (cid:18) k x (cid:19) cos (cid:32) √ k y (cid:33) +4 cos (cid:32) √ k y (cid:33) , (D2d)and the Θ k function (the argument of φ k ) is written as: tan Θ k = sin k x − k x cos √ k y cos k x + 2 cos k x cos √ k y . (D3)We can calculate ∇ k Θ k through: ∇ k Θ k = u ∇ k v − v ∇ k uE ( k ) , (D4)where we have split the function φ k in Eq. (D1) into realand imaginary parts φ k = u + iv . It then follows thatwe can obtain the components of the gradient of the Θ k function as: ∂ k x Θ k = 1 − (cid:16) √ k y (cid:17) + cos (cid:0) k x (cid:1) cos (cid:16) √ k y (cid:17) (cid:0) k x (cid:1) cos (cid:16) √ k y (cid:17) + 4 cos (cid:16) √ k y (cid:17) , (D5a) ∂ k y Θ k = √ (cid:0) k x (cid:1) sin (cid:16) √ k y (cid:17) (cid:0) k x (cid:1) cos (cid:16) √ k y (cid:17) + 4 cos (cid:16) √ k y (cid:17) . (D5b) Appendix E: Semi-analytical formula for thecharge-charge correlation function
In the long wavelength limit, the susceptibility for finitefrequency is written in power of q . If we expand ρ k + q and ω k , q until order q , we have: χ intrapump ( q , ω ) = 2 e (cid:126) a (cid:88) λ (cid:90) d k (2 π ) λ ∇ k ρ k · q ω − λ ∇ k E k · q + iγ , (E1)where we used that N ( k , q ) = 1 + O ( q ) . Expandingalso the denominator we find: χ intrapump ( q , ω ) = 4 e (cid:126) a ω + iγ ) (cid:90) d k (2 π ) ∇ k ρ k · q ∇ k E k · q , (E2)and thus we can define: C ij = (cid:90) d k (2 π ) ∂ i ρ k ∂ j E k . (E3)In terms of C ij we rewrite Eq. (E2) as: χ intrapump ( q , ω ) = 4 e (cid:126) a ω + iγ ) (cid:88) ij C ij q i q j . (E4) Making q x = q cos ϕ and q y = q sin ϕ it follows that: (cid:88) ij C ij q i q j = C xx cos ϕ + C yy sin ϕ + ( C xy + C yx ) sin ϕ cos ϕ . (E5)With an integration by parts we can show from Eq. (E3)that C xy = C yx . Using trigonometric identities we canwrite Eq. (E5) as: (cid:88) ij C ij q i q j = C xx + C yy C xx − C yy ϕ + C xy sin 2 ϕ . (E6)Defining: f = 1 π C xx + C yy , (E7a) f m = 1 π (cid:115)(cid:18) C xx − C yy (cid:19) + C xy , (E7b) φ = arctan 2 C xy C xx − C yy , (E7c)the susceptibility (E4) is written as: χ intrapump ( q, ϕ, ω ) = 4 e ( f + f m cos(2 ϕ − φ )) (cid:126) a π q ( ω + iγ ) . (E8)It is then possible to defined an effective Fermi energy E eff F ( ϕ ) = E F + f + f m cos(2 ϕ − φ ) , which allows to writethe total susceptibility as: χ ( q, ϕ, ω ) = 4 eE eff F ( ϕ ) (cid:126) a π q ( ω + iγ ) . (E9)The last result has the same functional form on fre-quency as that of the charge-charge susceptibility in theindependent electron gas model. Appendix F: Current-current response function in theout-of-equilibrium regime
The tight-binding Hamiltonian for a system with twoatoms per unit cell (and thus two sub-lattices A and B),and considering only nearest neighbor hoping, is givenby: H = (cid:88) i,n t i,n ( A )ˆ a † R n ˆ b R n + δ i + h.c. . (F1)where we have Introduced the Peierls substituion as: t i,n ( A ) = t TB exp (cid:18) − i e A ( R n ) · δ i (cid:126) (cid:19) , (F2)9which is valid when the vector potential A ( R n ) changessmoothly with the position in the lattice. The vectors δ i connect two nearest neighbors atoms from different sub-lattices.Expanding the Hamiltonian in the field A up to secondorder we obtain: H = H + V + V , (F3)where H is the usual tight-binding hamiltonian (A1): H = t TB (cid:88) i,n ˆ a † R n ˆ b R n + δ i + h.c. , (F4)and V accounts for the paramagnetic contribution to thecurrent and V to the diamagnetic one: V = − iet TB (cid:126) (cid:88) n,i A n · δ i ˆ a † R n ˆ b R n + δ i + h.c. , (F5) V = − e t TB (cid:126) (cid:88) n,i ( A n · δ i ) ˆ a † R n ˆ b R n + δ i + h.c. . (F6)In the following we will neglect the diamagnetic term V as it only contributes to the Drude conductivity (we willreturn to this point later in the Appendix).Using the momentum basis (A3), we can write V as: V = et TB iN c (cid:126) (cid:88) k , k (cid:48) ,i (cid:32)(cid:88) n A n e i ( k (cid:48) − k ) · R n (cid:33) · δ i e i k (cid:48) · δ i ˆ a † k ˆ b k (cid:48) + h.c.(F7)and defining the Fourier transform of the vector potential: A ( q ) = 1 N c (cid:88) n e − i q · R n A n , (F8)and making the change of variables q = k − k (cid:48) it followsthat: V = − iet TB (cid:126) (cid:88) k , q (cid:88) i A ( q ) · δ i e i ( k − q ) · δ i ˆ a † k + q / ˆ b k − q / + h.c. . (F9) Using the identity A ( q ) = A ∗ ( − q ) (since A ( r ) is real)we obtain: V = − iet TB (cid:126) (cid:88) k , q (cid:88) i A ( q ) · δ i (cid:16) e i ( k − q / · δ i ˆ a † k + q / ˆ b k − q −− e − i ( k + q / · δ i ˆ b † k + q / ˆ a k − q / (cid:17) , (F10)or V = et TB (cid:126) (cid:88) k , q A ( q ) · (cid:16) ∇ k φ k − q / ˆ a † k + q / ˆ b k − q / ++ ∇ k φ ∗ k + q / ˆ b † k + q / ˆ a k − q / (cid:17) , (F11)now we use the basis that diagonalizes the Hamiltonian H (A4) with the choice of global phase ϕ k = 0 ; thus itfollows: V = e (cid:88) k , q A ( q ) · v intra k , q (cid:16) ˆ c † k + q / ˆ c k − q / − ˆ d † k + q / ˆ d k − q / (cid:17) ++ A ( q ) · v inter k , q (cid:16) ˆ c † k + q / ˆ d k − q / − ˆ d † k + q / ˆ c k − q / (cid:17) , (F12)where we have defined: E k = t TB | φ k | , (F13) v inter k , q = a t TB (cid:126) (cid:16) e − i Θ k − q / ∇ k φ k − q / − e i Θ k + q / ∇ k φ ∗ k + q / (cid:17) , (F14) v intra k , q = a t TB (cid:126) (cid:16) e − i Θ k − q / ∇ k φ k − q / + e i Θ k + q / ∇ k φ ∗ k + q / (cid:17) . (F15)Using the Hamiltonian (F3) in Eq. (5), neglecting thecontribution of the term V , and taking the expectationvalue of the resulting equation for finite k and q we arriveat the coupled integral equations (Bloch equations) forthe elements of the density matrix: (cid:0) i∂ t − ω intra k , q − iγ (cid:1) n c k , q = − iγ f c k δ q , + (cid:88) q (cid:48) (cid:110) − Ω intra k − ( q + q (cid:48) ) , q (cid:48) n c k − q (cid:48) , q + q (cid:48) + Ω intra k + ( q + q (cid:48) ) , q (cid:48) n c k + q (cid:48) , q + q (cid:48) −− Ω inter k − ( q + q (cid:48) ) , q (cid:48) p cv k − q (cid:48) , q + q (cid:48) − Ω inter k + ( q + q (cid:48) ) , q (cid:48) p vc k + q (cid:48) , q + q (cid:48) (cid:111) , (F16a) (cid:0) i∂ t + ω intra k , q − iγ (cid:1) n v k , q = − iγ f v k δ q , + (cid:88) q (cid:48) (cid:110) +Ω intra k − ( q + q (cid:48) ) , q (cid:48) n v k − q (cid:48) , q + q (cid:48) − Ω intra k + ( q + q (cid:48) ) , q (cid:48) n v k + q (cid:48) , q + q (cid:48) ++Ω inter k − ( q + q (cid:48) ) , q (cid:48) p vc k − q (cid:48) , q + q (cid:48) + Ω inter k + ( q + q (cid:48) ) , q (cid:48) p cv k + q (cid:48) , q + q (cid:48) (cid:111) , (F16b)0 (cid:0) i∂ t − ω inter k , q − iγ p (cid:1) p cv k , q = (cid:88) q (cid:48) (cid:110) − Ω intra k − ( q + q (cid:48) ) , q (cid:48) p cv k − q (cid:48) , q + q (cid:48) + Ω intra k + ( q + q (cid:48) ) , q (cid:48) p cv k + q (cid:48) , q + q (cid:48) ++Ω inter k − ( q + q (cid:48) ) , q (cid:48) n c k − q (cid:48) , q + q (cid:48) − Ω inter k + ( q + q (cid:48) ) , q (cid:48) n v k + q (cid:48) , q + q (cid:48) (cid:111) , (F16c) (cid:0) i∂ t + ω inter k , q − iγ p (cid:1) p vc k , q = (cid:88) q (cid:48) (cid:110) − Ω intra k − ( q + q (cid:48) ) , q (cid:48) p vc k − q (cid:48) , q + q (cid:48) , +Ω intra k + ( q + q (cid:48) ) , q (cid:48) p vc k + q (cid:48) , q + q (cid:48) −− Ω inter k − ( q + q (cid:48) ) , q (cid:48) n v k − q (cid:48) , q + q (cid:48) + Ω inter k + ( q + q (cid:48) ) , q (cid:48) n c k + q (cid:48) , q + q (cid:48) (cid:111) , (F16d)where two phenomenological relaxation times havebeen introduced and where: n c k , q = (cid:104) c † k + q / c k − q / (cid:105) , (F17a) n v k , q = (cid:104) d † k + q / d k − q / (cid:105) , (F17b) p cv k , q = (cid:104) c † k + q / d k − q / (cid:105) , (F17c) p vc k , q = (cid:104) d † k + q / c k − q / (cid:105) , (F17d) Ω inter/intra k , q = e A ( q ) · v inter/intra k , q / (cid:126) , (F17e) ω inter k , q = E k − q / + E k + q / , (F17f) ω intra k , q = E k − q / − E k + q / . (F17g)We now consider a field given by a pumping compo-nent with frequency ω p and null wavenumber and a prob-ing component with frequency ω and wavenumber q : A n ( R n , t ) = 12 (cid:104) A p e iω p t + A e i ( ωt − q · R n ) + h.c. (cid:105) , (F18)with the Fourier transform: A ( q ) = A p cos( ω p t ) δ q , + A (cid:0) δ q , q e − iωt + δ − q , q e iωt (cid:1) , (F19)Next we simplify the set of Eqs. (F16). First we con-sider the effects of the pumping field on the new elec-tronic distribution n c k , , in a similar way to what has beendone in Sec. II, that is we neglect the effect of A since itis much smaller than A p . Using the result that the elec-tronic distribution converges to a steady-state we find: n c k , = f c k + iτ (cid:10) Ω inter k , (cid:0) p cv k , + p vc k , (cid:1)(cid:11) t , (F20a) n v k , = f v k − iτ (cid:10) Ω inter k , (cid:0) p cv k , + p vc k , (cid:1)(cid:11) t , (F20b) (cid:0) i∂ t − ω inter k , + iγ p (cid:1) p cv k , = Ω inter k , (cid:0) n c k , − n v k , (cid:1) , (F20c) (cid:0) i∂ t + ω inter k , + iγ p (cid:1) p vc k , = − Ω inter k , (cid:0) n v k , − n c k , (cid:1) . (F20d)Making the ansatz that the off-diagonal part of the p vc and p cv tensor oscillates with the same frequency ω p itfollows from Eqs. (F20c) and (F20d) that: p vc k , = A ( ω p ) e iω p t + B ( ω p ) e − iω p t , (F21) p cv k , = A ( ω p ) e iω p t + B ( ω p ) e − iω p t , (F22)and A = 12 (cid:126) e A p · v inter k , − ω p − ω inter k , + iγ p (cid:0) n c k , − n v k , (cid:1) , (F23a) B = 12 (cid:126) e A ∗ p · v inter k , ω p − ω inter k , + iγ p (cid:0) n c k , − n v k , (cid:1) , (F23b) A = 12 (cid:126) e A p · v inter k , − ω p + ω inter k , + iγ p (cid:0) n c k , − n v k , (cid:1) , (F23c) B = 12 (cid:126) e A ∗ p · v inter k , ω p + ω inter k , + iγ p (cid:0) n c k , − n v k , (cid:1) , (F23d)and if we define: β k = τ τ p (cid:12)(cid:12) Ω inter ( ω ) (cid:12)(cid:12) (cid:16) τ p (cid:16) ω + ω inter k , (cid:17)(cid:17) τ p ( ω − ω inter k , ) + 2 τ p ( ω + ω inter k , ) + 1 , (F24)where: ˜Ω inter k ( ω ) = e A p · v inter k , / (cid:126) = ea ω inter k , E · ∇ k Θ k ω (cid:126) (F25)where τ p = 1 /γ p , we have: n c k − f c k = β k ( n v k − n c k ) , (F26) n v k − f v k = β k ( n c k − n v k ) , (F27)1that are the equivalent of Eqs. (C6) in the Couloumbgauge, except that the Rabi frequency, ˜Ω inter k ( ω ) , hasa factor ω k /ω with relation to Eq. (C4): ˜Ω inter ( ω ) = ω k /ω ¯Ω k .Returning to the set of Eqs. (F16) we now considerthe effect of the two fields. Since the pumping fieldhas zero momentum it only contributes to the out-of-equilibrium distribution function calculated above. Inpractical terms this amounts to replace the equilibriumdistribution functions by the out-of-equilibrium ones forthe optical conductivity due to the probe. Following thisprocedures with the new electronic distribution we cancalculate the electronic current. For q = q , we have (cid:0) i∂ t − ω intra k , q + iγ (cid:1) n c k , q = Ω intra k , − q (cid:16) n c k − q , − n c k + q , (cid:17) , (F28a) (cid:0) i∂ t + ω intra k , q + iγ (cid:1) n v k , q = Ω intra k , − q (cid:16) n v k + q , − n v k − q , (cid:17) , (F28b) (cid:0) i∂ t − ω inter k , q + iγ p (cid:1) p cv k , q = Ω inter k , − q (cid:16) n c k + q , − n v k − q , (cid:17) , (F28c) (cid:0) i∂ t + ω inter k , q + iγ p (cid:1) p vc k , q = Ω inter k , − q (cid:16) n c k − q , − n v k + q , (cid:17) . (F28d)In the rotating wave approximation Eq. (F28) becomes: n c k , − q = − e A · v intra k , q (cid:126) n c k − q , − n c k + q , ω − ω intra k , − q + iγ e − iωt , (F29a) n v k , − q = e A · v intra k , q (cid:126) n v k − q , − n v k + q , ω + ω intra k , − q + iγ e iωt , (F29b) p cv k , − q = e A · v inter k , q (cid:126) n c k − q , − n v k + q , ω − ω inter k , − q + iγ p e − iωt , (F29c) p vc k , − q = − e A · v inter k , q (cid:126) n v k − q , − n c k + q , ω + ω inter k , − q + iγ p e iωt , (F29d)We can write the current in the n th cell: ˆ j n = − δV δ A n = et TB N c (cid:88) k , k (cid:48) ,i δ i (cid:104) e − i ( k − k (cid:48) ) · R n e i k (cid:48) · δ i ˆ a † k ˆ b k (cid:48) ++ e i ( k − k (cid:48) ) · R n e − i k (cid:48) · δ i ˆ b † k (cid:48) ˆ a k (cid:105) , (F30)after making the Bogoliubov transformation (A4) in Eq.(F30) and taking the expectation value we have: J n = e (cid:126) N c (cid:88) k , q e i q · R n v intra k , q (cid:0) n c k , q − n v k , q (cid:1) − v inter k , q (cid:0) p cv k , q − p vc k , q (cid:1) . (F31)The current up to first order in the probing field reads: J n ( R n , t ) = 12 (cid:16) J e i ( ωt − q · R n ) + h.c. (cid:17) , (F32)and finally the current J can be written as: J = 2 e (cid:126) N c (cid:88) k (cid:40) − v intra k , − q A · v intra k , q (cid:16) n c k − q , − n c k + q , (cid:17) ω − ω intra k , − q + iγ − v intra k , − q A · v intra k , q (cid:16) n v k − q , − n v k + q , (cid:17) ω + ω intra k , − q + iγ ++ v inter k , − q A · v inter k , q (cid:16) n c k − q , − n v k + q , (cid:17) ω − ω inter k , − q + iγ p + v inter k , − q A · v inter k , q (cid:16) n v k − q , − n c k + q , (cid:17) ω + ω inter k , − q + iγ p (cid:41) , (F33)where we included a factor of 2 to account for the spindegeneracy. The relation between the electric field andthe potential vector in frequency space is: A = E iω . (F34) We can finally write the conductivity tensor as: σ inter ij ( q , ω ) = 2 ie (cid:126) ωS (cid:88) k ,λ = ± λ n v k + λ q / − n c k − λ q / ω − λω inter k , q + iγ p v i inter k , q v j inter k , − q , (F35) σ intra ij ( q , ω ) = 2 ie (cid:126) ωS (cid:88) k ,λ = ± n λ k − q / − n λ k + q / ω − λω intra k , q + iγ v i intra k , q v j intra k , − q , (F36)2where we defined S = N c a with N c the number of unitcells in the crystal and n λ k = n λ k , is expected value of thediagonal element of the density matrix, with λ = + ( − )the conductance (valence) band. Introducing the rela-tion ω (¯ ω − λω k , q ) = − ωλω k , q + 1 λω k , q (¯ ω − λω k , q ) (F37)in Eqs. (F35) and (F36), where ¯ ω = ω + iγ , we find thatthe term proportional to − / ¯ ω has exactly the same form—that is, proportional to / ¯ ω — as the diamagnetic termwe have ignored, and therefore these two terms shouldgrouped. This procedure allows the determination theregular part of the conductivity. Finally, we have for theregular part the conductivity tensor the results: σ R,inter ij ( q , ω ) = 2 ie (cid:126) S (cid:88) k ,λ = ± ( n c k + λ q / − n v k − λ q / ) v i inter k , q v j inter k , − q ω inter k , q ( ω − λω inter k , q + iγ p ) , (F38) σ R,intra ij ( q , ω ) = 2 ie (cid:126) S (cid:88) k ,λ = ± ( n λ k − q / − n λ k + q / ) v i intra k , q v j intra k , − q λω intra k , q ( ω − λω intra k , q + iγ ) , (F39)whose both real and imaginary parts are not divergentwhen ω → . The divergent piece of the conductiv-ity when ω → is associated to the imaginary part ofthe neglected contributions and is nothing but the Drudeconductivity. Note that when q → , σ intra ij ( q , ω ) is minusthe neglected term and therefore they cancel each otherin that limit.Let us now return to the contribution to the total con-ductivity of the term V , termed diamagnetic contribu-tion. Following the same steps that we used to calculatethe paramagnetic term due to V , we find that the dia-magnetic current in the n th unit cell is given by: ˆ j dia n = − δV δ A n = e t TB (cid:126) (cid:88) i ( A n · δ i ) δ i ˆ a † R n ˆ b R n + δ i + h.c. . (F40)After using the Fourier representation of the creation andannihilation operators, and the Bogoliubov transforma-tion that diagonalizes H , we find: ˆ j dia n = e t TB (cid:126) N c (cid:88) k , q (cid:32)(cid:88) i δ i ( A n · δ i ) e i ( k + q / · δ i (cid:33) e − i Θ k − q / e i q · R n (cid:40) ˆ c † k + q / ˆ c † k − q / − ˆ d † k + q / ˆ d † k − q / − ˆ c † k + q / ˆ d † k − q / + ˆ d † k + q / ˆ c † k − q / (cid:41) + h.c. . (F41)Next we take the expectation value of the previous equa-tion, insert a factor to account for the spin degener-acy, truncate at first order in the field A n the expectationvalue, and perform a summation in q , obtaining: J dia n = 2 e t TB (cid:126) N c (cid:88) k (cid:32)(cid:88) i e i k · δ i δ i ( A n · δ i ) e − i Θ k + h.c. (cid:33) × (cid:0) n c k , − n v k , (cid:1) . (F42)The term inside the first pair of braces can be rewrittenwith the aid of the following relation: X ij = 12 t TB (cid:88) l e i k · δ l δ l · u i ( δ l · u j ) e − i Θ k + h.c. , (F43)with u i a versor in the direction i = x, y . The function(F43) can be rewritten as [note that k in Eq. (F43) isdimensionful and therefore the derivatives in order to thecomponents of k are also dimensionful] : X ij ( k ) = t TB a (cid:60) (cid:2) e − i Θ k ∂ i ∂ j φ ( k ) (cid:3) . (F44) Finally the diamagnetic term of the conductivity read as: σ dia ij ( ω ) = − ie (cid:126) S ( ω + iγ p ) (cid:88) k X ij ( k ) (cid:0) n c k , − n v k , (cid:1) . (F45)Putting all together, we have the Drude term for finite q : σ D ij ( q , ω ) = σ dia ij ( ω ) + σ D,intra ij ( q , ω ) + σ D,inter ij ( q , ω ) , (F46)where σ D,intra ij ( q , ω ) = − ie (cid:126) S (cid:88) k ,λ = ± ( n λ k − q / − n λ k + q / ) v i intra k , q v j intra k , − q λω intra k , q ( ω + iγ ) , (F47) σ D,inter ij ( q , ω ) = − ie (cid:126) S (cid:88) k ,λ = ± ( n v k + λ q / − n c k − λ q / ) v i inter k , q v j inter k , − q ω inter k , q ( ω + iγ p ) . (F48)The total conductivity is the sum of all the three terms: σ ij ( q , ω ) = σ D ij ( q , ω ) + σ R,inter ij ( q , ω ) + σ R,intra ij ( q , ω ) . Let usstress again that in the limit q → we have σ D,intra ij ( q → , ω ) = − σ R,intra ij ( q → , ω ) , in which case the Drude con-3ductivity is given solely by σ D ij (0 , ω ) = σ dia ij ( ω ) + σ D,inter ij (0 , ω ) . (F49) This result concludes the discussion of the optical re-sponse of graphene within the tight-binding approxima-tion. ∗ [email protected] † peres@fisica.uminho.pt ‡ [email protected] P. A. George, J. Strait, J. Dawlaty, S. Shivaraman, M. Chan-drashekhar, F. Rana, and M. G. Spencer, Nano Letters ,4248 (2008). M. Mittendorff, T. Winzer, E. Malic, A. Knorr, C. Berger, W. A.de Heer, H. Schneider, M. Helm, and S. Winnerl, NanoLetters , 1504 (2014). I. Gierz, J. C. Petersen, M. Mitrano, C. Cacho, I. E. Turcu,E. Springate, A. Stöhr, A. Köhler, U. Starke, and A. Caval-leri, Nature materials , 1119 (2013). E. Malic, T. Winzer, E. Bobkin, and A. Knorr, Phys. Rev. B , 205406 (2011). B. Y. Sun, Y. Zhou, and M. W. Wu, Phys. Rev. B , 125413(2012). T. Li, L. Luo, M. Hupalo, J. Zhang, M. C. Tringides,J. Schmalian, and J. Wang, Phys. Rev. Lett. , 167401(2012). A. F. Page, F. Ballout, O. Hess, and J. M. Hamm, Phys. Rev.B , 075404 (2015). G. Ni, L. Wang, M. Goldflam, M. Wagner, Z. Fei, A. McLeod,M. Liu, F. Keilmann, B. Özyilmaz, A. C. Neto, F. Fogler, andD. Basov, Nature Photonics , 244 (2016). J. M. Hamm, A. F. Page, J. Bravo-Abad, F. J. Garcia-Vidal,and O. Hess, Phys. Rev. B , 041408 (2016). A. Singh, T. Satpati, K. I. Bolotin, G. Saikat, and A. Agarwal,arXiv:1606.05072 [cond-mat.mes-hall]. A. Kumar, A. Nemilentsau, K. H. Fung, G. Hanson, N. X.Fang, and T. Low, Phys. Rev. B , 041413 (2016). Z. Fei, A. Rodin, G. Andreev, W. Bao, A. McLeod, M. Wag-ner, L. Zhang, Z. Zhao, M. Thiemens, G. Dominguez, et al.,Nature , 82 (2012). J. Chen, M. Badioli, P. Alonso-González, S. Thongrat-tanasiri, F. Huth, J. Osmond, M. Spasenovi´c, A. Centeno,A. Pesquera, P. Godignon, et al., Nature , 77 (2012). V. G. Kravets, A. N. Grigorenko, R. R. Nair, P. Blake, S. Anis-simova, K. S. Novoselov, and A. K. Geim, Phys. Rev. B ,155413 (2010). U. Fano, Rev. Mod. Phys. , 74 (1957). T. Winzer, E. Mali´c, and A. Knorr, “Graphene bloch equa-tions,” in Low-Dimensional Functional Materials, edited by R. Egger, D. Matrasulov, and K. Rakhimov (Springer, Dor-drecht, 2013) pp. 35–61. J. Shang, Z. Luo, C. Cong, J. Lin, T. Yu, and G. G.Gurzadyan, Applied Physics Letters , 163103 (2010). A. T. Roberts, R. Binder, N. H. Kwong, D. Golla, D. Cor-mode, B. J. LeRoy, H. O. Everitt, and A. Sandhu, Phys.Rev. Lett. , 187401 (2014). K. Oum, T. Lenzer, M. Scholz, D. Y. Jung, O. Sul, B. J. Cho,J. Lange, and A. Müller, The Journal of Physical ChemistryC , 6454 (2014). D. Giovanni, G. Yu, G. Xing, M. L. Leek, and T. C. Sum,Opt. Express , 21107 (2015). F. Rana, Phys. Rev. B , 155431 (2007). C. H. Lui, K. F. Mak, J. Shan, and T. F. Heinz, Phys. Rev.Lett. , 127404 (2010). S. Winnerl, M. Orlita, P. Plochocka, P. Kossacki, M. Potem-ski, T. Winzer, E. Malic, A. Knorr, M. Sprinkle, C. Berger,W. A. de Heer, H. Schneider, and M. Helm, Phys. Rev. Lett. , 237401 (2011). L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A.Bechtel, X. Liang, A. Zettl, Y. R. Shen, and F. Wang, Nat.Materials , 630 (2011). N. M. R. Peres, Rev. Mod. Phys. , 2673 (2010). N. D. Mermin, Phys. Rev. B , 2362 (1970). P. A. D. Gonçalves and N. M. R. Peres,An Introduction to Graphene Plasmonics (World Scien-tific, 2016). M. Jablan, H. Buljan, and M. Soljaˇci´c, Phys. Rev. B ,245435 (2009). Y. V. Bludov, A. Ferreira, N. M. R. Peres, and M. I.Vasilevskiy, Int. J. of Mod. Phys. B , 1341001 (2013). S. Xiao, X. Zhu, B.-H. Li, and N. A. Mortensen, Frontiers ofPhysics , 1 (2016). F. Wakaya, T. Teraoka, T. Kisa, T. Manabe, S. Abo, andM. Takai, Microelectronic Engineering , 144 (2012). M. Currie, J. D. Caldwell, F. J. Bezares, J. Robinson, T. An-derson, H. Chun, and M. Tadjer, Applied Physics Letters99