Semiclassical fluid model of nonlinear plasmons in doped graphene

Strathprints is designed to allow users to access the research output of the University of Strathclyde. Unless otherwise explicitly stated on the manuscript, Copyright © and Moral Rights for the papers on this site are retained by the individual authors and/or other copyright owners. Please check the manuscript for details of any other licences that may have been applied. You may not engage in further distribution of the material for any profitmaking activities or any commercial gain. You may freely distribute both the url (https://strathprints.strath.ac.uk/) and the content of this paper for research or private study, educational, or not-for-profit purposes without prior permission or charge.


I. INTRODUCTION
7][8] Two-dimensional (2D) plasmons have been observed in monolayer graphite, 9 single layer graphene on SiC substrates, 10,11 epitaxial graphene, 12 etc.6][17][18][19][20] The wave frequency of low-energy 2D plasmons is typically proportional to the square root of the wavenumber, and observed 2D plasmons have wave-frequencies most commonly in the THz range or below.The properties of massless Dirac fermions also lead to that the frequency of 2D plasmons in graphene scales as the electron density raised to 1/4, 13,14 in contrast to that of Schr€ odinger electrons 21,22 where the frequency of 2D plasmons scales as the square root of the electron density.9][40][41] The aim of this paper is to develop a nonlinear fluid model of plasmons in graphene by using an adiabatic closure of the governing kinetic model.The model captures the main collective and statistical effects (e.g., the equation of state for the pressure) depending on the shape and evolution of the electron distribution function, as well as nonlinear effects leading to harmonic generation.Comparisons are made with fluid models of plasmons using a density of state closure of the fluid hierarchy, 42,43 with the linear dielectric response using random phase approximations, 13,14,21 and with the highly nonlinear response proposed for applications of harmonic generation. 23,24,37

II. SEMI-CLASSICAL KINETIC MODEL
We here derive fluid-like equations for massless Dirac fermions starting from a semi-classical kinetic model.MKS (meter-kilogram-second) units are used throughout.While undoped graphene has the conduction band empty and the valence band filled at zero temperature, with no free carriers and a zero Fermi energy, doped or gated graphene has free carriers that support plasmons and screening of Coulomb interactions. 6,14The latter is the subject of our study, where the free carriers consist of massless Dirac fermions.We adopt the semi-classical one-particle Hamiltonian H¼v F jpjÀe/ðr; tÞ for massless Dirac fermions with constant speed v F ,w h e r e v F % 10 6 m / si st h eF e r m iv e l o c i t y ,p is the momentum, / is the potential, and -e is the electron charge.It is one of the remarkable characteristics that the Fermi velocity is constant for massless Dirac fermions in graphene, in contrast to Schr€ odinger electrons where the Fermi velocity is densitydependent.Hamilton's equations then lead to the collision-less one-particle Boltzmann (or Vlasov) equation The potential is given by 6 /ðr; tÞ¼À e 4p 0 ð ðn e ðr 0 ; tÞÀn 0 Þ jr À r 0 j where n 0 represents the areal number density of the neutralizing ion background, n e ðr; tÞ¼ Ð fd 2 p is the electron areal number density, 0 is the electric vacuum permittivity, and is the mean dielectric constant of the surrounding medium. 6he Vlasov equation (3) describes the evolution of the distribution function of massless Dirac fermions having velocity v ¼ v F p=jpj and constant speed jvj¼v F in the presence of the mean-field potential /.This simplest kinetic model may be expanded by the inclusion of various nonideal effects. 44aking into account the quantum tunneling of the electrons at small scales, the Vlasov equation would be replaced by a Wigner equation for a pseudo-distribution function, 45,46 and taking into account spin effects where the single-particle wave-function is a 2-spinor would give rise to a 2 Â 2 matrix of coupled kinetic equations 47,48 or an extension of the phase space to include the spin polarization variable. 49,50For a fully degenerate 2D electron gas, it has, however, been found that by calculation of the gradient correction to the Fermi-Thomas model the quantum tunneling vanishes. 51,52The present model takes into account the effect of low-dimensionality on the long-range Coulomb interactions 16 and the shape of the electron distribution function in momentum space of a twodimensional electron gas (2DEG), 21 but neglects the shortrange interactions between electrons that are important for plasmons in graphene at short wavelengths. 13,14Linearizing and Fourier analyzing Eqs. ( 3) and ( 4) with f ¼ f 0 þ f 1 and / ¼ / 1 , and f 1 and / 1 proportional to exp ðiq Á r À ixtÞ, yields the dielectric function e e ¼ 1 þ v e , where the electron susceptibility is given by and the second step is obtained by an integration by parts.The equilibrium distribution function f 0 for the massless Dirac fermions can be taken to be a Fermi-Dirac momentum distribution function 23,24 and a a normalization constant, which in the ultra-cold limit bl ) 1 takes the form where H(n) is the Heaviside function which equals one for n > 0 and zero for n < 0, g s ¼ 2 and g v ¼ 2 are the spin and valley degeneracies for graphene, l ¼ hv F k F is the chemical potential (the same as the Fermi energy in the ultra-cold limit), and k F ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4pn 0 =g s g v p is the Fermi wavenumber.The equilibrium distribution function f 0 is normalized such that Ð f 0 d 2 p ¼ n 0 .The ultra-cold limit is assumed here for simplicity, but the obtained results can be extended to finite temperatures in a straightforward manner.Using the distribution function (6), the integral in Eq. ( 5) is evaluated to give with q ¼jqj, and where k s ¼ 0 l/(e 2 n 0 ) is a screening length.The dispersion relation for plasmons, 1 þ v e ¼ 0, then gives The dependence of the wave frequency on the wavenumber is plotted as a solid line in Fig. 1.
For small wavenumbers k s q < 1, we have from Eq. ( 8) the approximate dispersion relation In the last step, we denoted the density-dependent effective mass m Ã ¼ l=v 2 F ¼ hk F =v F , and for the comparisons of models below, a correction parameter ¼ 3/4 is introduced.In experiment, 10 the measured density n 0 ¼ 1.9 Â 10 14 m À2 gives k F ¼ ffiffiffiffiffiffiffi pn 0 p ¼ 0:077 ÅÀ1 and m * ¼ 0.089m e , in reasonable agreement with the obtained best-fit value m * ¼ 0.077m e .I n Eq. ( 9), the first term on the right-hand side represents the space-charge effect, while the second term is due to the shape of the distribution function in momentum space.The typical square-root behavior of the frequency on wavenumber, keeping only the first term on the right-hand side of Eq. ( 9),i s seen in the dashed curve in Fig. 1.The dispersion relation (9), plotted as a dash-dotted curve in Fig. 1, is on the same form as obtained by Stern 21 for a 2D electron gas using a random phase approximation.In the fluid model of Ref. 42 and subsequently used in Ref. 43, ¼ 1/2 (instead of 3/4) was obtained using a density-of-states closure for the pressure.In graphene, it is predicted using a random phase approximation 13,14 that band overlap of the wavefunctions leads to a negative FIG. 1. Dispersion curves for plasmons in graphene obtained from kinetic theory (8) (solid black line), the approximate dispersion relation (9) (dashdotted blue line), the RPA result (11) (solid magenta line), and the longwavelength limit x ¼ v F ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q=ð2k s Þ p (dashed red line), for the case of graphene on a SiC substrate with n 0 ¼ 1.9 Â 10 17 m À2 and ¼ 5.5, for which k F k s % 0.63.The dotted lines show the Fermi velocity cone x ¼ v F q, and the boundary x ¼ 2l= h À v F q between the undamped region I and the Landau damped regions II and III due to interband and intraband transitions, respectively.Starting from the electron susceptibility in the high-frequency (x > v F k) long wavelength limit 14 v e;RPA ¼À qv 2 one obtains from the dispersion relation 1 þ v e,RPA ¼ 0 that As a comparison, the wave frequency obtained from Eq. ( 11) is plotted in Fig. 1 (solid magenta line), showing a frequency down-shift compared with the 2DEG case due to the shortrange electron interactions in graphene.Keeping only the first correction for small q in Eq. ( 11) yields where a ee ¼ 1=ð4k s k F Þ¼e 2 =ð4p 0 hv F Þ is the fine structure constant; 6 hence, the correction parameter is negative, ¼Àa 2 ee .Taking both the positive and negative corrections into account, Ref. 18 obtained ¼ 3=4 À a 2 ee .It was pointed out in Ref. 18 that can be either positive or negative depending on the value of a ee .For example, a suspended graphene sheet with ¼ 1 gives a ee ¼ 2.2 and ¼ -4.0, while for graphene with one side exposed to air ( 1 ¼ 1) and the other to SiC ( 2 ¼ 10) giving ¼ ( 1 þ 2 )/2 ¼ 5.5 yields a ee ¼ 0.40 and ¼ 0.59.The measured plasmon dispersion in Ref. 10 seems, however, to be consistent with the negative correction obtained in Refs.13 and 14 at long wavelengths.Even though x/q > v F for all wavenumbers in our model (and hence no intraband transitions), it should be noted that Landau damping is predicted to take place due to interband transitions 13,14,18 in the region This restricts the validity of our semiclassimodel at large wavenumbers.As an example, interband transitions take place in region II in Fig. 1 for the case of graphene on a SiC substrate 10 with ¼ 5.5 and n 0 ¼ 1.9 Â 10 17 m À2 .The dispersive and dissipative effects due to nonlinearity 27 are more complex, but may include interband transitions due to harmonic generation, etc.

III. NON-LINEAR FLUID MODEL DERIVED FROM KINETIC THEORY WITH ADIABATIC CLOSURE
Next, nonlinear fluid equations are derived using moments of the Vlasov equation.Integrating Eq. (3) over momentum space gives the continuity equation @n e @t þrÁðn e v e Þ¼0; where the fluid density n e and velocity field v e are defined as and respectively.The fluid density n e enters into Eq.( 4) for the scalar potential.To obtain an equation for the fluid velocity, the Vlasov equation is multiplied by p and integrated over momentum space @ðn e p e Þ @t þ v F rÁ ð pp jpj fd 2 p À en e r/ ¼ 0; (16)   where the fluid momentum is defined as and pp denotes a dyadic product.The second term in the left-hand side of Eq. ( 16) represents the momentum flux due to the mean flow and statistical motion of the electrons.Equation ( 15) will be used below to connect v e to p e and n e .Higher moments of the Vlasov equation can be found by multiplying Eq. ( 3) by pp and integrating over momentum space, and so on, leading to a hierarchy of fluid equations.
Here, a closure of the fluid hierarchy is found by assuming f to be on a specific form that approximately solves the Vlasov equation (3).Since the Vlasov equation describes an incompressible flow in phase space, the key idea is that the phase fluid density (the distribution function f)d o e sn o t change in amplitude but only in width and offset in momentum space.Since the phase speed of the plasmon is larger than the Fermi velocity, the electrons are not able to stream to neutralize the wave potential.Similar to a three-dimensional electron gas, 53 the high-frequency plasmons thus lead to a deformation of the Fermi surface.We use a shifted Fermi-Dirac distribution function 54 to model adiabatic compression, keeping the chemical potential constant.Temporarily choosing a coordinate system such that the spatial variation is along the x-axis, the distribution function is assumed to be of the form where p ex (x, t) represents a shift of the mean momentum, and the density variations n e (x, t) are achieved by a widening and contractions of the distribution function along the p x component in momentum space.The particular form of Eq. ( 18) ensures that Eqs. ( 14) and ( 17) are fulfilled exactly.Equation ( 15) is used to derive a nonlinear relation between v e , p e and n e , where p ey ¼ 0 and v ey ¼ 0 due to symmetry.Using that p=jpj¼r p jpj, the x-component of the right-hand side of Eq. ( 15) is readily integrated over p x , with the result and density n e , and when multiplied with the electron charge gives the x-component of the electric current density, j x ¼ -en e v ex .The nonlinear dependence of the current density on momentum has been proposed for the generation of odd harmonic in the case of a monochromatic electromagnetic driving field. 23,24It should be noted that Eq. ( 19) can be writt eno nasi m il arfo r masE q .(6) of Ref. 23: Using the identity and a change of integration variable p y ¼ð l=v F Þ sin u; dp y ¼ð l=v F Þ cos u du in Eq. ( 19) gives, after some straightforward manipulations where Here, P F ¼ v F p ex /l ¼ p ex /p F is the dimensionless fluid momentum (equal to minus the field parameter Q F in Ref. 23) and p F ¼ l/v F the Fermi momentum.Equation ( 20) becomes of the same form as Eq. ( 6) of Ref. 23 for homogeneous density n e ¼ n 0 .Solving Eq. ( 20) for the fluid velocity, a reasonably good approximation is where c takes the role of an effective semi-relativistic gamma factor, and G ¼ 4/(1 þ 3n e /n 0 ) takes into account the compression of the electron density.Figure 2 shows v ex as a function of p ex for different values of n e /n 0 using Eq. ( 20) and the numerical fit in Eq. ( 22).In the small amplitude limit (with n e /n 0 ¼ 1), the linear relationship is v ex ¼ p ex /m * .O n the other hand, in the large amplitude limit p ex v F /l ) n e /n 0 , we have v ex ¼ v F sign(p ex ), consistent with the results in Refs.23 and 24.
With the distribution function given in Eq. ( 18), the Eq. ( 16) reduces to @ @t ðn e p ex Þþ @P @x À en e @/ @x ¼ 0; (23)   where the pressure P ¼ P dyn þ P stat is expressed as the sum of the dynamic pressure mainly due to the mean flow of the electrons and the statistical pressure due to electron degeneracy with x being a unit vector in the x-direction.By a change of integration variable, the fluid momentum p ex is eliminated from the expression for P stat .Then, introducing polar coordinates in momentum space and carrying out the integration along the radial momentum direction, we arrive at the equation of state where, in the last step, an accurate Pad e approximation was used for the integral, giving a small relative error less than 1% for 0.4 < n e /n 0 < 4.
To evaluate the dynamic pressure, we carry out the integral over p x in Eq. ( 24) to obtain where An approximation of Eq. ( 27) is P dyn 'n 0 p 2 ex =M, where M ¼16m Ã =ð9þn e =n 0 Þ depends weakly on the compression of the electron density.The choice of compression along the x-direction was arbitrary, and in 2D, the momentum equation becomes @ðn e p e Þ @t þrÁ n 0 p e p e M þrP stat À en e r/ þrw kin ¼ 0; (28)   with / given by Eq. ( 4) and P stat by Eq. ( 26).The velocity field in the continuity equation ( 13) is approximated as [cf.
Eq. (22)] Equations ( 13) and ( 28), coupled with Eqs. ( 4) and ( 29), constitute the nonlinear fluid Higher order dispersive corrections based on the exact kinetic dispersion relation (8)  have been taken into account by the addition of the rw kin term in the left-hand side of Eq. ( 28), with the Fourier component of w kin given by ŵkin ¼ an e and a ¼ (l/8)k In real space, w is obtained via a convolution product between (n e -n 0 ) and the inverse Fourier transform of a(k); the latter can be expressed in terms of special functions but is omitted here for brevity.
To assess the validity of the model, we study the linearized system with n e ¼ n 0 þ n e1 ; P ¼ P 0 þ P 1 ; p e ¼ p e1 ; v e ¼ v e1 ; / ¼ / 1 .The linearized statistical pressure is P stat;1 ¼ð3=4Þln e1 , while the fluid velocity and momentum are linearly related as v e1 ¼ p e1 =m Ã ¼ðv 2 F =lÞp e1 .Then, assuming n e1 , p e1 and / 1 to be proportional to exp ðik Á r À ixtÞ, and using that the Fourier component / of the potential is related to the electron density ne as / ¼Àen e =ð2 0 jkjÞ, the Fourier components are eliminated from the linearized system.Neglecting w kin in Eq. ( 28) leads to a dispersion relation identical to Eq. ( 9), obtained from the Vlasov-Poisson system in the long wavelength limit and previously derived using a random phase approximation. 21The inclusion of the w kin term yields the exact kinetic dispersion relation Eq. (8).
Using a simplified nonlinear fluid model of a 2D plasma layer 41 it was found that large-amplitude waves give rise to a non-linear frequency up-shift (similar to Stokes waves), and a modulational instability due to self-amplifying sidebands of the carrier wave.Such an investigation of the present fluid model, including the effects of the nonlinear relation (29)  between v e , p e , and n e , is planned for future research.

IV. CONCLUSIONS
We have derived a non-linear fluid model for highfrequency plasmons in graphene, by taking moments of the collision-less Boltzmann (Vlasov) equation, and finding a closure corresponding to adiabatic compression and rarefaction of the electron density in the wave-field.Higher order linear dispersive effects are added to the fluid model by using the exact linear dispersion relation obtained from semiclassical kinetic theory.We emphasize that the nonlinearity in Eq. ( 29) and in the second and third terms of Eq. ( 28) give rise to nonlinear currents that can be used for high-order harmonic generation and frequency multiplication in the response of harmonic large-amplitude electromagnetic waves, 23,24 and may lead to the localization of plasmon wave energy via modulational instabilities 41 and to wave localization and the formation of solitons. 39,40A positive nonlinear frequency shift through self-interaction of waves was predicted in Ref. 41 for a 2D electron gas, but for graphene, the increase of the effective mass through the semi-relativistic gamma factor in Eq. ( 29) may lead to a negative frequency shift, in analogy with large-amplitude Langmuir waves 55,56 and electromagnetic waves 55,57,58 in 3D classical plasmas where the relativistic mass increase of the electrons leads to an effective decrease in the plasma frequency.Thus, a few applications of the nonlinear fluid model can be to study the nonlinear wave-wave interaction between plasmons and the nonlinear coupling and scattering between plasmons and electromagnetic waves in graphene.

FIG. 2 .
FIG. 2. The normalized fluid velocity v ex /v F as a function of the normalized fluid momentum p ex /(m * v F ) for n e /n 0 ¼ 0.5, 1 and 1.5 as indicated in the figure.The solid lines show the velocity using Eq.(20) while the dashed lines show the approximation in Eq. (22).