Tunable surface plasmon instability leading to emission of radiation

Possible sources of terahertz (THz) radiation have been investigated for several years now. For example, by epitaxially growing layers of different III–V semiconductors, THz quantum-cascade lasers, usually with a fixed frequency, have been fabricated.¹ The existence of instability in an electrically driven plasma demonstrates new approaches for efficient energy conversion as well as for designing and fabricating tunable THz emitters. Initially, the instability of a two-dimensional (2D) plasma above a conducting half-space was addressed by Krasheninnikov and Chaplik.² However, we show that when a layer lies above a thick conductor, it is incorrect to simply replace one of the frequencies in a two-layer plasma dispersion equation by the surface plasmon frequency of the underlying substrate, as suggested by Krasheninnikov and Chaplik.² It was later predicted by Kempa et al.³ (see also Ref. 4) that when a current is passed through a stationary two-component 2D electron gas (2DEG), the Doppler shift in response frequency leads to a spontaneous generation of plasmon excitations without a wavenumber cutoff and subsequent Cherenkov radiation⁵ at sufficiently high carrier drift velocities. In the present work, however, the Doppler shift in frequency is found not in general applicable to massless Dirac fermions in graphene at arbitrary wavelengths. For our model, the process is proved irreversible based on the spontaneous breaking of the time reversal symmetry. Similar conclusions are expected for monolayer graphene that is characterized by massless Dirac fermions or a synthesized nanosheet of silicene.⁶ In the same group of the periodic table with graphene, doped silicene is predicted to exhibit similar electronic properties. Additionally, it has the advantage over graphene in its compatibility with Si-based device technologies.

In addition to the work cited above, the role played by plasma excitations in the THz response of 2D microstructures has received considerable attention,^7–16 as well as plasmon properties in fullerenes.^17–19 Plasmon modes of quantum-well transistor structures in the THz range may be excited with the use of far-infrared (FIR) light and other means.²⁰ A split-grating-gate design has been found to significantly enhance FIR response.^21–24 Under this scheme, however, the stimulated electromagnetic (EM) radiation requires either a population inversion²⁵ quantum coherence²⁶ or condensation.²⁷ The EM radiation can also be generated by transferring energy from optical field to another.²⁸ An alternative theoretical model for a THz emitter was proposed in Ref. 29, which is based on a periodic spatial modulation to the velocity of an incident electron beam, similar to Bloch oscillations in a dc-field driven semiconductor superlattice. Initially, the generation of a plasma wave by applying an electron beam based on induced plasma instabilities was addressed in Ref. 3 and the references therein. Consequently, the incident electrons will be scattered by the generated plasmon excitations and may no longer be considered to be ballistic. Such phenomena are expected to be experimentally observable in electron transport measurements. The instability could also be established by converting the plasma-wave energy to electromagnetic radiation by employing a surface grating. This would allow coupling of the near field of a plasma wave with a radiation field, which may then be used to design a tunable THz emitter having variable photon energy and wavelength by varying the driving current.

Graphene has been considered as a promising material for photo-detectors due to its unique electronic properties. A terahertz detector based on an antenna-coupled field effect transistor was demonstrated using the non-linear response to an oscillating radiation field.^30–32 

Here, we explore a new energy conversion approach, that is, from a dc electric field to an optical field based on a current-induced plasma instability with a tunable wavenumber cutoff. The primary objective of this paper is to expand in a careful way the materials platform by exploiting the new functionality of a composite nano-system consisting of a thick conductor (including a heavily doped semiconductor) that is Coulomb coupled to 2D layered materials. This coupling has recently been investigated experimentally and was shown to lead to interesting effects on the plasmon excitations such as nonlinear dispersion and damping.^33–35 The Coulomb coupling of the plasmons in a layer to the surface plasmon on the conductor results in a tailored surface plasmon instability that leads to the emission of radiation (spiler). The instability is found always following a characteristic loop in the dispersion beyond a critical in-plane wave vector for two current-split low-energy plasmon branches. The predicted tunable spiler radiation relies on a current-induced plasmon instability and comes after the self-growth of the plasmon localized field in the time domain at a rate which is determined by the surface-plasmon frequency, the 2D layer separation, the distance of the 2D layers from the conducting surface, and the driving-current strength. This emission tunability provides us with new tools for spatially imaging and spectrally identifying a single molecule adsorbed on the graphene layer.

In our formalism, we consider a nano-scale system consisting of a pair of 2D layers and a thick conducting material. The layer may be monolayer graphene or a 2DEG such as a semiconductor inversion layer or HEMT (high electron mobility transistor). We define a 2D crystalline layer as a layer having a low-frequency plasmon mode dispersion in the form of $Ω≃q$⁠. A conventional 2DEG, graphene, newly discovered buckled honeycomb lattices, such as silicene and germanene, may be suitable components for our hybrid system. The graphene layer may have a gap, thereby extending the flexibility of the composite system that also incorporates a thick layer of conducting material as depicted in Fig. 1. The excitation spectra of allowable modes will be determined from a knowledge of the non-local dielectric function $ϵ(r,r′;ω)$⁠, which depends on position coordinates $r,r′$ and frequency ω or its inverse $K(r,r′;ω)$ satisfying $∫dr′ K(r,r′;ω) ϵ(r′,r″;ω)=δ(r−r″)$⁠. The self-consistent field equation for $K(r,r′;ω)$ is now determined.

In operator notation, the dielectric function for the pair of 2D layers and semi-infinite structure is given by

where $ϵ̂=K̂−1$ with $K̂$ as the inverse dielectric function satisfying

In our notation, $α̂2D(i)$ j and $α̂SI$ are the polarization functions of the composite system and the polarization function of the 2D layer and semi-infinite substrate, respectively. Additionally, $K̂SI$ is the inverse dielectric function for the semi-infinite substrate whose surface lies in the z = 0 plane. In integral form, after Fourier transforming with respect to coordinates parallel to the xy-plane and suppressing the in-plane wave number $q||$ and frequency ω, we obtain for a double layer positioned at z = a₁ and z = a₂ (⁠$0<a1<a2$⁠) interacting with each other as well as the semi-infinite conducting substrate with its surface located at z = 0. Making use of

where $v(z,z′)=(2πe2/ϵs) exp(−q∥|z−z′|)$ with $ϵs=4πϵ0ϵb$ in terms of the permeability ϵ₀ of free space and the backgrounds dielectric constant ϵ_b. Our calculation shows that

By setting z₁ = a₁ and z₁ = a₂ in turn in Eq. (4) and solving the pair of simultaneous equations for $K(a1,z2)$ and $K(a2,z2)$⁠, we obtain

where $Sc(2)(q||,ω)=DetM↔(q||,ω)$ with

Substituting the results for $K(a1,z2)$ and $K(a2,z2)$ into Eq. (4), we obtain the complete inverse dielectric function for a pair of 2D planes interacting with each other and a semi-infinite conducting material. The plasmon excitation frequencies are determined by the zeros of $Sc(2)(q||,ω)$⁠. Furthermore, the effect of the inverse dielectric function for the semi-infinite structure $KSI(z,z′;q||,ω)$ screens coupling within and between two layers. As a matter of fact, our result for the plasmon dispersion relation generalizes that obtained by Das Sarma and Madhukar³⁶ for a bi-plane. We carry out calculations using the local limit for K_SI, which is³⁷ 

where $ϵB(ω)=1−ωp2/ω2$ in terms of the bulk plasma frequency ω_p.

In the present case, we obtain in the local limit

With $Π2D(0)(q||,ω)≈Cq∥2/ω2$ at T = 0 K and in the long-wavelength limit (⁠$q∥≪kF$ with k_F as the Fermi wave vector), we now introduce further notation, $C¯j=2πe2Cj/(ϵsωp2)$ for j = 1, 2. The spectral function yields real frequencies. A plane interacting with the half-space has two resonant modes. Each pair of 2D layers interacting in isolation far from the half-space conducting medium supports a symmetric and an antisymmetric mode.³⁶ In the absence of a driving current, the analytic solutions for the plasmon modes for a pair of 2D layers that are Coulomb coupled to a half-space are given by

where $A≡(C¯1a1−C¯2a2)2+4C¯1C¯2a12$ and the term $4C¯1C¯2a12$ plays the role of “Rabi coupling.”

Clearly, for long wavelengths, only $Ω1(q∥)$ depends on ω_p. However, the excitation spectrum changes dramatically when a current is driven through the configuration. Under a constant electric field, the carrier distribution is modified as may be obtained by employing the relaxation time approximation in the equation of motion for the center-of-mass momentum (see  Appendix). For carriers in a parabolic energy band with effective mass $m∗$ and drift velocity v_d determined by the electron mobility and the external electric field, the electrons in the medium are redistributed. This is determined by a momentum shift in the wave vector $k∥→k∥−m∗vd/ℏ$ in the thermal-equilibrium energy distribution function $f0(ϵk)$⁠. By making a change of variables in the well-known Lindhard polarization function $Π2D(0)(q∥,ω)$⁠, this effect is equivalent to a frequency shift $ω→ω−q∥·vd$⁠. For massless Dirac fermions in graphene with linear energy dispersion, this Doppler shift in frequency is not in general valid for arbitrary wave vector as we prove in our  Appendix. This is our conclusion after we relate the surface current density to the center-of-mass wave vector in a steady state. Our calculation shows that the redistribution of electrons leads to a shift in the wave vector appearing in the Fermi function by the center-of-mass wave vector $K0=(kF/vF)vd$⁠, where k_F and v_F are the graphene Fermi wave vector and velocity, respectively. However, in the long-wavelength limit, $q∥→0$⁠, the Doppler shift in frequency is approximately obeyed. This is discussed in  Appendix. Consequently, regardless of the nature of the 2D layer represented in the dispersion equation, we may replace $ω→ω−q∥·vd$ in the dispersion equation in the presence of an applied electric field at long wavelengths.

In general, we shall treat the solution frequencies $ω±(q∥)$ as complex variables with $Im[ω±(q∥)]≥0$⁠, where $Im[ω±(q∥)]>0$ implies a finite growth rates $γ±(q∥)=Im[ω±(q∥)]$ for two split plasmon modes. Since $ϵ(q∥, ω)$ is a complex function, we ask for $Re[ϵ(q∥, ω)]=Im[ϵ(q∥, ω)]=0$⁠. Therefore, we are left with damping-free plasmon modes in the system but they still face possible instability due to $Im[ω±(q∥)]>0$ (or spontaneous time-reversal symmetry breaking due to self-growth of the plasmon localized field).

We shall treat the solution frequencies $ω±(q∥)$ as complex variables with $Im[ω±(q∥)]≥0$⁠, where $Im[ω±(q∥)]>0$ implies a finite growth rates $γ±(q∥)=Im[ω±(q∥)]$ for two split plasmon modes. Since $ϵ(q∥, ω)$ is a complex function, we ask for $Re[ϵ(q∥, ω)]=Im[ϵ(q∥, ω)]=0$⁠. Therefore, we are left with damping-free plasmon modes in the system, but they still face possible instability due to $Im[ω±(q∥)]>0$⁠. Since the instability only occurs for $q∥vd$⁠, hereafter, we will simply denote the component of wave vector $q$ along the current direction as $q∥$⁠.

Setting $a1=a$ and letting $a2→∞$ in Eq. (6), the off-diagonal matrix elements tend to zero and the element in the first row and first column reduces to unity. Subsequently, the dispersion equation for a single layer interacting with the substrate is given by the zeros of the matrix element in the second row and second column. Using the long-wavelength limit (⁠$q≪kF$⁠), we find $Π2D;j(0)(q∥,ω)≈−Cjq∥2/ω2$⁠. As shown above, we have $C=n2D/m2D∗$ for a 2DEG; $C=(2μ/πℏ2)[1−(Δ2/μ2)]$ for doped graphene. However, for intrinsic graphene whose plasmon excitations are induced by temperature, we obtain³⁸ $C=(2ln2) kBT/πℏ2$⁠. Consequently, using the notation for $C¯$ defined above, we rewrite the plasmon frequency as follows:³⁹ $ω2=K1±K2$ with K₁ and K₂ defined by $K1=C¯q∥ωp2/2+(ωp/2)2$ and $K2=C¯q∥ωp4exp(−2q∥a)/2+(ωp/2)4(1−2C¯q∥)2$⁠. Additionally, within this long-wavelength limit, these expressions yield the plasmon excitation frequencies $ω1/ωp≃q∥2C¯a$ and $ω2/ωp≃1/2+C¯q∥/2$⁠, which are both linear in $q∥$ and unlike the $q∥$-dependence for free-standing graphene or the 2DEG.^40–45 

In Ref. 46, it was demonstrated that the interband and intraband π plasmon excitations in graphene have interesting depencies on the in-plane wave transfer momentum. Specifically, a linear dispersion rather than a square root dependence on the wave vector was reported. This result came as a surprise because theoretical calculations on free-standing graphene clearly do not predict a linear dependence in the long-wavelength limit for the low-frequency plasmon excitations. As a matter of fact, this linear dependence of plasmon frequency on wave vector was attributed to local field corrections to the random-phase approximation.

In order to highlight the effect due to a surface in spiler, as a comparison we first present in Figs. 2(a) and 2(b) the plasmon dispersion for an isolated pair of 2D layers in the absence³⁶ and presence³ of a current, respectively. When a substrate surface plasmon is not contributing [in Figs. 2(a) and 2(b)], the plasmon instability starts at $q∥=0$ and exists over a finite range until a bifurcation point is reached for $q∥$⁠. By going beyond this bifurcation point, the interlayer Coulomb coupling is effectively suppressed, leading to one uncoupled 2D-sheet plasma and two current-split 2D-sheet plasmon modes. However, when a surface plasmon interacts with a single 2D layers [in Figs. 2(c) and 2(d)], the plasmon instability may be moved to shorter wavelengths, as we clearly illustrated in Fig. 2(d), where we show the results when spiler consists of a 2D layer and a semi-infinite conducting medium. Plasmon remains stable in an isolated single current-driven 2D layer (see the discussion of Fig. 3). In the presence of the surface plasmon, as $vd/ωp=0.6C¯1$ in Fig. 2(d), the plasmon instability occurs outside the closed Rabi “loop,” that is, $q∥>qc$ (q_c is the critical wave vector). Here, q_c is pushed significantly above zero to a shorter wavelengths, leaving a finite wavenumber cutoff.

We endeavor to obtain a full understanding of the mechanism for instability shown in Fig. 2 by numerically investigating the effect of a passing current through a layer of 2DEG, graphene, or silicene which is Coulomb coupled to a conductor. Specifically, we consider a pair of 2D layers and a semi-infinite medium such as a heavily doped semiconductor. We present in Fig. 4 both the plasmon dispersion and the growth rate. Each panel shows results for a different drift velocity given by $vd/ωp=0, 0.5C¯1, 0.8C¯1$⁠, and $1.6C¯1$ as labeled. In the absence of a current, panel (a) of Fig. 4 shows that there are three plasmon excitation branches, which are stable as given by Eq. (8). At low v_d, panel (b) of Fig. 4 demonstrates that the plasmons are still remained stable. However, as v_d is increased further, the lowest branches may become unstable as in Figs. 4(c) and 4(d) through the appearance of a positive imaginary part for the complex frequency. There exists a threshold value for v_d beyond which the plasmon excitation becomes unstable. On the other hand, the existence of the surface plasmon greatly screens both the interlayer and intralayer Coulomb couplings as $q∥a1≪1$ in the range of $ω/ωp≪1$⁠. This stabilizes the plasmon excitation for $q∥<qc$ by suppressing the interlayer coupling as shown in Fig. 4(c). As v_d increases to $1.6C¯1ωp$ in Fig. 4(d), q_c reduces almost to zero, but the wide stable region in Fig. 4(c) is squeezed into a narrow belt. The occurrence of such a new unstable region starting from $q∥=0$ is a combined result from both the surface-induced softening of the two 2D-sheet plasma modes to two acoustic-like plasmon modes as well as the strong interlayer coupling for small layer separation. To some extent, this first instable region is similar to the result displayed in Fig. 2(b) for the isolated pair of 2D layers.

Figure 5 illustrates our results for larger 2D layer separations from each other compared to the case in Fig. 4. Panel (a) of Fig. 5 again corresponds to the analytic solutions in Eq. (8) when v_d = 0. The plasmon excitations are all still stable in Fig. 5(b) at $vd/ωp=0.4C¯1$⁠. However, as v_d is increased further in Figs. 5(c) and 5(d), an instability appears at q_c, which is exactly where the characteristic Rabi “loop” for plasmon excitations closes in the two lower panels. Beyond q_c, the lower two current-split plasmon branches become unstable due to switching from the quadruple-like energy-storage mode to dipole-like energy-radiation mode. The loop shape is quite different for current passing through the 2D layer [in Fig. 5(c)] and the semi-infinite medium [in Fig. 5(d)]. We also note that by adjusting the layer separation, we may change q_c value for controlling the onset of the plasmon instability and the wavenumber cutoff. With such a large interlayer separation, the Coulomb interaction between two 2D layers becomes very weak and the system effectively behaves like the single current-driven 2D layer coupled to a conducting surface, similar to that in Fig. 2(d).

Physically, from the point of view of momentum space, electrons may only occupy momentum space within the range of $|k∥|≤kF$ at zero temperature in a state of thermal equilibrium, where k_F is the electron Fermi wave number and $ε(kF)=ε(−kF)=EF$ is the Fermi energy. When a low current is passed through the electron gas, electrons are slightly driven out from this thermal-equilibrium state and their population becomes asymmetrical with respect to $k∥=0$⁠. In this case, the Fermi energy E_F is split into $EF,+=ε(kF+K0)$ and $EF,−=ε(−kF+K0)$ with $EF,+>EF,−$⁠, where $ℏK0$ represents the electron center-of-mass momentum. In this shifted Fermi-Dirac distribution picture, electrons in such a non-equilibrium state are energetically unstable, and the higher-energy electrons in the range $kF≤k≤kF+K0$ tend to decay into lower-energy empty states by emitting EM waves and phonons to ensure the conservations of total momentum and energy.

The current-driven asymmetric electron distribution in $k∥$ space leads to an induced oscillating polarization current or a “dipole radiator.” If two electron gas layers are placed close enough, the in-phase interlayer Coulomb interaction will give rise to a dipole-like plasmon excitation that releases energy by radiation, similar to that of a single layer. On the other hand, the out-of-phase interlayer Coulomb coupling will lead to a quadruple-like excitation that stores energy like a “capacitor.” This quadruple-like plasmon excitation produces a near EM field and can even be effectively converted into a transverse far EM field if a surface grating is employed.

One aim of our investigation is to show how the particle-hole modes change in the presence of finite drift velocity since they determine the regions of undamped plasmons. Our calculations show that both inter-band particle-hole modes' (PHM's) (above the diagonal) and the intra-band modes are clearly modified. We have demonstrated that the PHM regions below the diagonal $ω=vFq$ are suppressed, whereas the inter-band regions above the main diagonal are enhanced.

We now turn our attention to detailed numerical calculations for the plasmon excitations beyond the long wavelength limit, that is, for finite values of the wave vector $q∥$⁠, which are comparable with the Fermi wave vector k_F in graphene. The plasmon frequencies clearly depend on the doping concentration of each graphene layer, as well as the inter-layer separation and their distances from the surface. If these distances are equal to the inverse of k_F, we discover three distinct plasmon modes, two of which start at the origin (corresponding to the graphene layers) and the other one originated at a finite frequency with a depolarization shift from $ωp/2$ and attributed to the surface plasmon of the conductor. All three branches have finite group velocity in the limit $q→0$⁠, which is in agreement with our previous results in the long wavelength limit. If the individual distances are increased, as it is shown in Fig. 6(f), the plasmon branches tend to get close to each other at different values of $q∥$⁠, due to a strong decrease of the Coulomb interaction between each graphene layer and the conductor. Here, we limited our consideration to the real parts of the plasmon frequencies in order to make sure that all the allowed plasmon branches have been accounted for.

The surface-induced instability in spiler may lead to EM radiation, and the cutoff of its radiation wavenumber can be adjusted directly by q_c to give rise to a tunable emission spectrum. Here, we present in Fig. 3 the dependence of controlling parameter q_c as functions of v_d and separation a for a single driven 2D layer coupled to surface plasmon, where q_c is increased with either reducing v_d or increasing a. These results clearly demonstrate significant shifts of q_c within the desired ranges for operations of both photodetectors and EM-wave devices.

In summary, we are proposing a tunable spiler quantum plasmonic device that employs 2D layers in combination with a thick conducting material. We found that the spiler spontaneously emits EM radiation when a current is passed through the 2D layer or the underlying conductor to make the plasmons become unstable at a specific wave number corresponding to a wavenumber cutoff. The plasmon instability is defined through finite, non-zero imaginary part of the plasmon frequency, which indicates its inverse lifetime. One of the important features of this instability is the wave vector and frequency values at which it occurs. The non-zero critical wave vector makes our system (which includes a conductor substrate) distinguished from previous proposals and may be tuned through the coupling of the layer and surface plasmon excitations. It is possible to tune the onset of plasmon instability (or wavenumber cutoff) by selecting the properties of the nanosheet or frequency of the surface plasmon. The surface plasmon introduces a characteristic Rabi “loop” in the dispersion of two current-split low-energy plasmon branches and plays a crucial role in giving rise to the splitting and the concomitant streams of quasiparticles. We would also like to emphasize that the existence of the critical wave vector q_c has introduced a cutoff wavenumber in the dispersion of the emitted terahertz radiation. More importantly, this cutoff wavenumber can be tuned by the 2D-layer separation to the surface of the metal film. This leads to a new prospective device for single-molecule imaging and identification while it is adsorbed on a graphene layer.

The near field arising from the surface-plasmon excitation within our system effectively reduces the so-called diffraction limitation imposed on the spatial resolution in nano-imaging (see Refs. 48 and 49). Consequently, single-molecule imaging is made possible. In addition, the out-coupling of generated THz wave by plasmon instability in our system can be employed for detecting vibration modes of adsorbed molecules on a graphene layer and makes single-molecule identification achievable.

The emitted EM radiation is naturally associated with a near field and may even be coupled outward to free space by a grating on the surface. Experimentally, the emitted EM radiation of a spiler device can be detected by a heterodyne mixing technique using a planar Schottky diode,⁴⁷ where a high precision measurement of the EM radiation frequency can be measured and a transient turn-on behavior can be shown. The generated EM radiation in the THz range is able to penetrate materials that block visible light and have a wide range of possible applications, including chemical analysis, security scanning, medical imaging, and telecommunications. We emphasize that the reported work by Krasheninnikov and Chaplik² is flawed. Physically, one cannot simply replace one of the frequencies in a two-layer plasmon dispersion equation by the surface plasmon frequency. The proper mathematical formalism to include the surface plasmon interacting with a 2D layer was first presented by Horing et al.²⁰ Detailed numerical investigation of the plasmon excitation for such systems demonstrated the extended regions of undamped plasmons.⁵⁰ 

In addition, the work reported by Mikhailov²⁹ on voltage-tunable THz emitter relies on a structure to perform a periodic spatial modulation to the velocity of an incident electron beam. This modulation can be realized by introducing a biased metallic grating adjacent to the electron beam. In the present paper, on the other hand, our proposed THz emitter is based on the instability of electrically driven graphene plasmon modes, which does not require velocity modulation or broken translational symmetry of the system. In fact, our model demonstrates spontaneous breaking of the time-reversal symmetry due to self-growth of the plasmon localized field, whereas both the Smith-Purcell²⁹ and similar Bloch oscillations⁵¹ schemes involve no broken time-reversal symmetry but translational-symmetry breaking.

In this paper, we concentrated our investigation to the long-wavelength regime. If we go beyond this region, we must use the numerical results for the full expression of polarization derived in the  Appendix. At the same time, we must also use the nonlocal inverse dielectric function for the semi-infinite medium. Such a study, however, goes beyond the scope of the present work.

This research was supported by Contract No. FA 9453-13-1-0291 of AFRL. D.H. thank the Air Force Office of Scientific Research (AFOSR) for its support. W.P. was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. We thank Oleksiy Roslyak and Antonios Balassis for helpful discussions.

Let us first consider the case when a current is passed through a 2DEG layer. This current creates a Doppler shift in the response function $Π(q,ω−q·vd)$⁠. The derivation is presented in the paper by Kempa et al.³ and the argument is as follows. The energy dispersion for an electron in a 2DEG is $εk=ℏ2k2/2m∗$⁠, where $m∗$ is the electron effective mass. The current flow leads to the replacement in wave vector $k→k−m∗vd/ℏ$ everywhere in the polarization function

where $f0(ε)$ is the thermal-equilibrium distribution function for electrons, and the drift of electrons can be considered as a perturbation as long as $vd≪vF$ (with v_F as the Fermi velocity of electrons). After this wave vector replacement is carried out and a change of variables is made in the resulting integral, we simply obtain an expression for the polarization function that is exactly the same as that given in Eq. (A1), except with the frequency shifted by $q·vd$⁠.

We now turn to the case of graphene that is characterized by massless Dirac fermions for which the energy dispersion is linear in the wave vector k. For a spatially uniform system, the first-order moment of the Boltzmann equation for a single valley gives us⁵² 

where $Nc=1/(πℏ2vF2), η=μ0(t)/kBT$⁠, μ₀ is the chemical potential of electrons in graphene, the quantity

and

is the electron surface current density, $A$ is the sample area, T is the electron temperature, $ℏk$ is the electron wave vector, $εk=ℏvF k$ is the electron kinetic energy, v_F is the Fermi velocity of graphene, $vk=∇kεk/ℏ$ is the electron group velocity, τ₀ is the average momentum-relaxation time, and $F(t)$ is the external electric field. Additionally, we have the following relation:

where $ρ(t)$ is the electron areal density. Considering a steady state under a constant electric field $F0$⁠, we obtain

At $T≈0$ K, we have

where $EF=ℏvF kF$ is the electron Fermi energy and $kF=2πρ0$ is the Fermi wave number. As a result, this leads to

In the second-quantization picture, the Hamiltonian operator for N electrons in graphene in the presence of an electric field may be written as

where $σ→=(σx, σy)$ is the Pauli-matrix vector, $p̂j=−iℏ∇j$ is the electron momentum operator, and r_j is the electron position vector. For this system, we define the center-of-mass momentum operator as $P̂c=∑j=1N p̂j$⁠. Therefore, the Heisenberg equation gives us

where we have employed the momentum-relaxation time approximation. For a steady state, we have

Finally, we are able to connect the electron surface current density with the center-of-mass wave vector in a steady state simply through

Recalling that we have $j=ρ0 vd$⁠, where v_d is the drift velocity of electrons in the system, we arrive at the relation

Consequently, for drifted electrons, we find from the Lindhard polarization function that

where v_d is determined by the product of the electron mobility and the external electric field $F0$⁠, and

As $v/vd≪1$ for $T≈0$ K, we get

If we further assume $Δ≪ℏvFkF$⁠, this leads to

In a similar way, we get for $q≪k≈kF$

As a result, we find

Finally, the prefactor is given by