In this paper, we theoretically investigate plasmon propagation characteristics in AB and AA stacked bilayer graphene (BLG) in the presence of energy asymmetry due to an electrostatic field oriented perpendicularly to the plane of the graphene sheet. We first derive the optical conductivity of BLG using the Kubo formalism incorporating energy asymmetry and finite electron scattering. All results are obtained for room temperature (300 K) operation. By solving Maxwell's equations in a dual gate device setup, we obtain the wavevector of propagating plasmon modes in the transverse electric (TE) and transverse magnetic (TM) directions at terahertz frequencies. The plasmon wavevector allows us to compare the compression factor, propagation length, and the mode confinement of TE and TM plasmon modes in bilayer and monolayer graphene sheets and also to study the impact of material parameters on plasmon characteristics. Our results show that the energy asymmetry can be harnessed to increase the propagation length of TM plasmons in BLG. AA stacked BLG shows a larger increase in the propagation length than AB stacked BLG; conversely, it is very insensitive to the Fermi level variations. Additionally, the dual gate structure allows independent modulation of the energy asymmetry and the Fermi level in BLG, which is advantageous for reconfiguring plasmon characteristics post device fabrication.

The two-dimensional (2D) material graphene offers strong light-matter interactions over a broad electromagnetic spectrum, ranging from several hundreds of GHz to optical frequencies.1 An interesting and useful attribute of graphene is its ability to support the propagation of surface plasmon polaritons (SPPs). SPPs (referred to as plasmons, hereafter) are collective oscillations of electrons, typically occurring at the interface between a metal and a dielectric, that can be coupled with electromagnetic waves and enable the manipulation of light below its diffraction limit.2 Plasmon oscillations can occur in various metallic and semiconducting material systems including 2D electron gases.3–5 Typically, plasmons have a smaller wavelength compared with the wavelength of light at the same frequency allowing significant mode compression.

Plasmons supported on a graphene sheet have a tighter confinement and a higher propagation length relative to those of surface plasmons in metals.6,7 Moreover, graphene can also support the propagation of transverse electric (TE) modes,8 which are absent in the case of a metal/dielectric system. The characteristics of graphene plasmons can also be tuned via electrostatic gating or chemical doping over a wide frequency range. In particular, the electrostatic tuning of plasmons enables reconfigurable plasmonic devices with graphene—a feature that is unique to graphene and not afforded by metal plasmonics.

Recently, plasmon oscillations in monolayer and multilayer graphene heterostructures were experimentally imaged.9–13 These experiments pave the way to realize tunable plasmonic devices that are currently being researched and developed including waveguides,14 nanoantennas,15 photodetectors,16 ocsillators,17 and modulators.18 From a theoretical standpoint, plasmons have been studied extensively in gated19 and ungated20,21 monolayer graphene (MLG) and also in ungated22 bilayer graphene (BLG). However, not much attention has been given to plasmons in BLG structures in which the bandstructure and, therefore, optical properties can be readily tuned via electrostatic gating.

In this paper, we develop theoretical models to study the propagation of TE and TM (transverse magnetic) plasmons in gated BLG structures for various material parameters of graphene and the dielectrics surrounding the graphene sheet. We also highlight the distinguishing features of plasmon propagation in BLG and MLG. While both MLG and BLG exhibit tunable optical properties, the tunability is significantly enhanced in BLG. This enhancement is a result of the electric-field-induced energy asymmetry in BLG structures. Essentially, a vertical electric field in a metal-oxide-graphene geometry creates an energy asymmetry between the two parallel layers of BLG.23–25 The energy asymmetry changes the bandstructure of BLG, which in turn alters its optical properties significantly.26,27 In addition to the energy asymmetry, we also consider other phenomena, such as the electron scattering rate and the stacking forms, that affect the optical properties of BLG.

The remainder of this paper is organized as follows: In Sec. II, we derive the plasmon wavevector for both TE and TM polarized modes in a graphene sheet using the solutions of Maxwell's equations. We use a device setup in which the graphene sheet is modeled as an impedance boundary condition with a finite dynamical surface conductivity. To obtain the dynamical conductivity, we first derive the bandstructure of BLG using tight binding model in Sec. III. The bandstructure is obtained in terms of the energy asymmetry between the layers and for two common stacking orders, which are AB and AA stacked BLG. The bandstructure is used to solve for the dynamical conductivity of BLG within the Kubo formalism. The effects of gate voltage and electrostatic screening on the carrier concentration and energy asymmetry are included using a self-consistent Hartree approximation. We present the results of our analysis in Sec. IV in which we quantify the impact of physical parameters of BLG on the plasmon characteristics such as compression factor, propagation length, and mode confinement. The key findings of the paper are summarized in Sec. V.

To obtain plasmon wavevector in graphene heterostructures, we consider the device setup shown in Fig. 1 in which an infinite 2D sheet of graphene is immersed in a uniform, linear dielectric medium with a relative dielectric constant of εr. We note that graphene in this device setup could be either a monolayer or a bilayer sheet. The propagating electric and magnetic fields in the z-direction for TE-polarized plasmon mode are

{E=(Ey+ŷ)ej(kzz+kxx)H=(Hx+x̂+Hz+ẑ)ej(kzz+kxx)x>0,{E=(Eyŷ)ej(kzzkxx)H=(Hxx̂+Hzẑ)ej(kzzkxx)x<0,
(1)

and for TM-polarized plasmon mode are

{H=(Hy+ŷ)ej(kzz+kxx),E=(Ex+x̂+Ez+ẑ)ej(kzz+kxx)x>0,{H=(Hyŷ)ej(kzzkxx),E=(Exx̂+Ezẑ)ej(kzzkxx)x<0,
(2)

where time dependence term ejωt is dropped and kz(kx) is the wavevector in z–(x-) direction as illustrated in Fig. 1. If we plug Eqs. (1) and (2) into Helmholtz equation (2E+klight2E=0), we get kx2+kz2=klight2, where klight=ωεrε0μ0 is the wavevector of light in the medium. Here, ω is the radial frequency, and μ0 and ε0 are the permeability and permittivity, respectively, of free space.

FIG. 1.

Device setup to study the propagation of plasmons in a graphene sheet with a surface conductivity of σ(ω) immersed in a uniform dielectric medium (ϵr). The components of electric and magnetic fields for TE and TM Plasmon modes are labeled in the figure. Plasmons decay by ekx|x| in the direction perpendicular to the graphene sheet.

FIG. 1.

Device setup to study the propagation of plasmons in a graphene sheet with a surface conductivity of σ(ω) immersed in a uniform dielectric medium (ϵr). The components of electric and magnetic fields for TE and TM Plasmon modes are labeled in the figure. Plasmons decay by ekx|x| in the direction perpendicular to the graphene sheet.

Close modal

By solving Maxwell's equations in the device setup of Fig. 1 and applying appropriate boundary conditions (impedance boundary with a finite surface conductivity for graphene), we obtain the following dispersion relations:8,21

1+jμ0ωσ(ω)2kz2εrω2/c2=0,TEMode,
(3a)
1jσ(ω)kz2εrω2/c22εrε0ω=0,TMMode,
(3b)

where σ(ω) is a complex quantity that represents surface conductivity of the graphene sheet and c=1/ϵ0μ0 is the velocity of light in free space. By replacing σ(ω)=σ(ω)jσ(ω) and η=μ0/εrε0 in Eq. (3) and rearranging the equations, we obtain the plasmon wavevectors

kz=klight(1+η24(σ(ω)2σ(ω)2)+jη22σ(ω)σ(ω))1/2,
(4)

for TE mode and

kz=klight(1+4η2σ(ω)2σ(ω)2|σ(ω)|4j8η2σ(ω)σ(ω)|σ(ω)|4)1/2,
(5)

for TM mode. From Eqs. (4) and (5), it is evident that the propagation of TE (TM) mode requires σ(ω) to be negative (positive). This is due to the fact that for propagating plasmon modes, the imaginary part of kz, which represents the propagation loss, must be positive.

While plasmons propagate in both z-direction (longitudinal) and x-direction (transverse), they are expected to decay more rapidly in the x-direction. The decay of plasmons can be represented using a length scale called the propagation length. The propagation length is the distance at which the intensity of the wave degrades to 1/e of its initial value, where e = 2.718 is the base of natural logarithm. For given wavevectors kz=kzjkz and kx=kxjkx, the propagation length is obtained as 1/2kz and 1/2kx for longitudinal and transverse directions, respectively. The measure of confinement of the plasmon wave to the graphene sheet can be represented by its propagation length in the transverse direction. A lower transverse propagation length means higher mode confinement on the graphene surface. Another important characteristic of a plasmon wave is its compression factor, which determines how compressed the plasmon wavelength is relative to the wavelength of light calculated in the absence of the graphene sheet. The compression factor is mathematically given as the ratio of the longitudinal phase constant to the phase constant of light kz/klight.

An important ingredient in obtaining the plasmon dispersion is the dynamical conductivity, σ(ω), which depends on the bandstructure of graphene, electron scattering rate, Fermi level, and the operating frequency. In Sec. III, we present the model for σ(ω) for AB and AA stacked BLG. For brevity, we do not include the model for σ(ω) in MLG. Interested readers are referred to Refs. 21 and 28–31 for a detailed discussion pertaining to MLG conductivity.

To study the optical properties of graphene, in particular, its dynamical conductivity, we start with the calculation of bandstructure of BLG. Here, we investigate two of the most common stacking forms, namely, AB (Bernal stacked) and AA stacked BLG. The natural stacking form of BLG is AB, which is also found in graphite. Although AA stacked BLG is likely to be metastable,32 it is still interesting to study it from a theoretical point of view. Figure 2 illustrates the configuration of both the stacking forms. Each graphene layer consists of a triangular Bravais lattice with a basis of two atoms, which we denote as “A” and “B.” In an AB staked BLG, “A” atoms of one layer align with “B” atoms of the other layer, and all other atoms lie at the center of the other layer's hexagons. In an AA stacked BLG, all the atoms in one layer align exactly with other layer's atoms correspondingly.

FIG. 2.

Honeycomb structure of (a) AB stacked BLG and (b) AA stacked BLG. Top layer (dashed) and bottom layer (solid) are slightly displaced in (b) for clarity. The graphene lattice constant is denoted as b=2.46Å.

FIG. 2.

Honeycomb structure of (a) AB stacked BLG and (b) AA stacked BLG. Top layer (dashed) and bottom layer (solid) are slightly displaced in (b) for clarity. The graphene lattice constant is denoted as b=2.46Å.

Close modal
FIG. 3.

Bandstructure of (a) AB stacked BLG and (b) AA stacked BLG. The bandstructure is plotted for two different values of energy asymmetry Δ. In this figure, γ1 represents the interlayer hopping parameter. When both Δ and γ1 are equal to zero, the BLG bandstructures of Eqs. (7) and (9) reduce the energy dispersion of MLG. The energy bandgap Eg occurs at k=0.

FIG. 3.

Bandstructure of (a) AB stacked BLG and (b) AA stacked BLG. The bandstructure is plotted for two different values of energy asymmetry Δ. In this figure, γ1 represents the interlayer hopping parameter. When both Δ and γ1 are equal to zero, the BLG bandstructures of Eqs. (7) and (9) reduce the energy dispersion of MLG. The energy bandgap Eg occurs at k=0.

Close modal

By utilizing tight binding model, we can write the Hamiltonian of an AB stacked BLG that is exposed to a vertical electrostatic field. The electric field introduces an energy asymmetry between the two layers, which is denoted as Δ. It lowers the energy of the first layer by Δ/2 and raises the energy of the second layer by the same amount. The Hamiltonian for an AB stacked BLG is written in the following form:33 

Ĥ=[Δ/200ϕ*(k)0Δ/2ϕ(k)00ϕ*(k)Δ/2γ1ϕ(k)0γ1Δ/2]·
(6)

Here, γ1=0.4 eV is the interlayer hopping energy between B1 and A2 sites34 and |ϕ(k)|=|vFkeiθ| is the energy dispersion of MLG in the vicinity of the K point of the Brillouin Zone within a radius of 3 eV. Also, vF=3γ0b/2 denotes the Fermi velocity, γ0=3 eV is nearest neighbor hopping energy within each layer, and b = 0.246 nm is graphene's lattice constant. The eigenvalues of the Hamiltonian describe the bandstructure, which in the case of AB stacked BLG is given as

ϵα(k)=±γ122+Δ24+|ϕ(k)|2+(1)αγ144+|ϕ(k)|2(γ12+Δ2)·
(7)

Here, α=1,2 is the subband index. We can see that for γ1=0 and Δ = 0, this energy dispersion reduces to that of MLG. Figure 3(a) shows the bandstructure of AB stacked BLG for two different values of Δ. Contrary to MLG that exhibits a gapless linear energy dispersion relationship, AB stacked BLG shows a nonlinear energy dispersion with a finite energy bandgap. We note that the energy bandgap (Eg=Δγ1/Δ2+γ12) occurs at k=0 and is smaller than Δ. Similarly, we can write the Hamiltonian of AA stacked BLG for a finite electrostatic bias as follows:35 

Ĥ=[Δ/20γ1ϕ(k)0Δ/2ϕ*(k)γ1γ1ϕ(k)Δ/20ϕ*(k)γ10Δ/2]·
(8)

The eigenvalues of this matrix describe the bandstructure of AA stacked BLG, which is given as

ϵα(k)=±(γ12+Δ24+(1)α|ϕ(k)|).
(9)

Expression (9) shows that the energy asymmetry does not open a bandgap in the case of AA stacked BLG, and Δ can be treated as a correcting factor to γ1γ12+Δ2/4. Figure 3(b) shows the bandstructure of AA stacked BLG for two different values of Δ.

We can write the real part of frequency dependent conductivity σ(ω) using Kubo formalism. The following procedure is described in Ref. 33, and for an AB stacked BLG, we have

σAB(ω)=2e2ωdE2π[f(E)f(E+ω)]d2k(2π)2|vk|2[A11(E,Δ)A44(E+ω,Δ)+A44(E,Δ)A11(E+ω,Δ)+A11(E,Δ)A44(E+ω,Δ)+A44(E,Δ)A11(E+ω,Δ)+2(A13(E,Δ)A13*(E+ω,Δ)+A13*(E,Δ)A13(E+ω,Δ))]·
(10)

Here, μ is the Fermi level, f(E)=1/(1+exp((Eμ)/kBT)) is the Fermi-Dirac distribution function with kB as Boltzmann's constant and T as the temperature, and |vk| is the velocity of electrons and is approximated by vF. The inner integral is being evaluated over the K point in the Brillouin zone. The Aij functions are the spectral representations of Green's functions. They are related to Green's functions through the following equation:

Gij(z)=dw2πAij(w)zw,
(11)

and can be written in the form of

Aij(w,Δ)=α=1,2[aij(α,Δ)δ(wϵα)+aij(α,Δ)δ(w+ϵα)],
(12)

in which α denotes the αth subband and ϵα is the subband energy. The coefficients aij(α,Δ) are derived from Eqs. (11) and (12), and the Hamiltonian through Ĝ1=zÎĤ.

a11(α,Δ)=(1)α+1πϵ22ϵ12([(Δ2)2+γ12ϵα2](1Δ2ϵα)+|ϕ(k)|2(1+Δ2ϵα)),a44(α,Δ)=(1)α+1πϵ22ϵ12([(Δ2)2ϵα2](1Δ2ϵα)+|ϕ(k)|2(1+Δ2ϵα)),a13(α,Δ)=(1)απϵ22ϵ12(1Δ2ϵα)γ1ϕ(k)·
(13)

In a similar fashion, we can write the conductivity of AA stacked BLG, which is described in Ref. 35. The real part of conductivity is written as follows:

σAA(ω)=2e2ωdE2π[f(E)f(E+ω)]d2k(2π)24|vk|2[A11(E+ω)A11(E)+A13(E+ω)A13(E)],
(14)

in which the spectral functions are

A11(w,Δ)=α=1,2π2[δ(wϵα)+δ(w+ϵα)],
(15a)
A13(w,Δ)=α=1,2π2[δ(wϵα)δ(w+ϵα)]·
(15b)

To incorporate the effect of electron scatterings in graphene, the delta functions in Eqs. (12) and (15) are broadened using the Lorentzian defined as

δ(w)=1πΓw2+Γ2·
(16)

The electron scattering rate is represented by Γ, which is generally a function of energy and the bandstructure. In the present work, we use a constant value that represents the energy-dependent scattering rate averaged over energy. The impact of bandstructure and energy asymmetry on the scattering rate is discussed in  Appendix A. For all simulations in this work, Γ is chosen as 20 meV, unless otherwise stated, which is a typical experimental value of scattering rate in a graphene sheet.36 The imaginary part of conductivity is derived from its real part using Kramers-Kronig relation

σ(ω)=2ωπP0σr(ω)ω2ω2dω,
(17)

where P denotes the principal value of the integral.

Figure 4 shows the optical conductivity of AB and AA stacked BLG for different values of μ and Δ using an energy-independent value of electron scattering rate Γ = 20 meV at room temperature T = 300 K. The plots are normalized to the universal conductivity of graphene σ0=q2/4. In Figs. 4(a) and 4(c), the real parts of the optical conductivity (σ(ω)) of AB and AA stacked BLG are illustrated for different values of the Fermi level, μ at Δ = 0 eV. We can see that for a fixed energy asymmetry, the optical conductivity can be tuned by the Fermi level. For instance, in the case of AA stacked BLG, the interband transition is shifted in frequency as the Fermi level changes. In the case of AB stacked BLG, there is a strong absorption peak at frequencies around γ1 due to interband transitions. Figures 4(e) and 4(g) illustrate σ(ω) for different values of Δ at a fixed Fermi level μ=0.3 eV. From these plots, we can see that the impact of energy asymmetry on the optical conductivity of BLG is quite significant specially in the case of AB stacked BLG where the absorption peak shifts toward higher frequencies as Δ increases. This phenomenon can be seen as an additional degree of freedom for tuning the optical conductivity of BLG that is absent in the case of MLG. At frequencies close to the visible region, the real part of conductivity converges to 2σ0 as expected.37 The imaginary parts of the optical conductivity σ(ω) are plotted in Figs. 4(b), 4(d), 4(f), and 4(h). These figures show the existence of different frequency bands where TE plasmons (σ(ω)<0) and TM plasmons (σ(ω)>0) can be supported and tuned in AA and AB stacked BLG structures.

FIG. 4.

Optical conductivity of BLG, for different values of Fermi level μ and energy asymmetry Δ, numerically evaluated by using Eqs. (10), (14), and (17), normalized to σ0=q2/4. The figures illustrate the impact of μ on the optical conductivity of (a) and (b) AB stacked BLG, (c) and (d) AA stacked BLG for a fixed Δ and the impact of Δ on the optical conductivity of (e) and (f) AB stacked BLG, and (g) and (h) AA stacked BLG for a fixed μ. All simulations are conducted at 300 K with a constant scattering rate of Γ = 20 meV.

FIG. 4.

Optical conductivity of BLG, for different values of Fermi level μ and energy asymmetry Δ, numerically evaluated by using Eqs. (10), (14), and (17), normalized to σ0=q2/4. The figures illustrate the impact of μ on the optical conductivity of (a) and (b) AB stacked BLG, (c) and (d) AA stacked BLG for a fixed Δ and the impact of Δ on the optical conductivity of (e) and (f) AB stacked BLG, and (g) and (h) AA stacked BLG for a fixed μ. All simulations are conducted at 300 K with a constant scattering rate of Γ = 20 meV.

Close modal
FIG. 5.

Gated BLG within a uniform dielectric medium with relative permittivity of εr. Gate oxide thickness is denoted as tt for the top oxide and tb for the bottom oxide. The interlayer distance in BLG is denoted as d = 0.35 nm. Dual gate structure is used to control Fermi level μ and energy asymmetry Δ independently. The carrier density on graphene layers are represented by n1 and n2. The difference between oxide fields δE=EbEt determines the Fermi level and their average E¯=(Eb+Et)/2 introduces the energy asymmetry. We assume that the relative permittivity of BLG is equal to εr,g=1.

FIG. 5.

Gated BLG within a uniform dielectric medium with relative permittivity of εr. Gate oxide thickness is denoted as tt for the top oxide and tb for the bottom oxide. The interlayer distance in BLG is denoted as d = 0.35 nm. Dual gate structure is used to control Fermi level μ and energy asymmetry Δ independently. The carrier density on graphene layers are represented by n1 and n2. The difference between oxide fields δE=EbEt determines the Fermi level and their average E¯=(Eb+Et)/2 introduces the energy asymmetry. We assume that the relative permittivity of BLG is equal to εr,g=1.

Close modal

Experimental observations,26,27,38,39 as well as theoretical studies,23–25,40–42 have demonstrated that applying an electrostatic field perpendicular to the BLG plane forms an energy asymmetry between the two layers of graphene and consequently alters its bandstructure. In Secs. III A and III B, we discussed that the presence of energy asymmetry creates a bandgap in the bandstructure of AB stacked BLG, while it modifies the interlayer hopping energy γ1 of AA stacked BLG. These phenomena introduce an additional degree of freedom (along with the Fermi level) for tuning the optical conductivity of BLG. This electric-field-induced energy asymmetry can be realized by gating the graphene layers as shown in Fig. 5. We use the dual gate structure shown in this figure to obtain the energy asymmetry, Δ, in terms of top and bottom gate voltages (Vt and Vb, respectively). In this metal-oxide-graphene geometry, the gate voltage applied at the metal electrode modulates the energy asymmetry and the carrier density in BLG via the electrostatic field in the oxide. The physical thickness of the oxide is denoted as tt for the top oxide and tb for the bottom oxide. The parameter d = 0.35 nm is the interlayer distance of BLG and is equal to the thickness of each layer. The use of dual gates allows us to independently control the Fermi level and energy asymmetry in the BLG. The fields of top and bottom gates are Et and Eb, respectively. It has been shown that the difference between these two fields δE=EbEt determines the Fermi level, and their average E¯=(Eb+Et)/2 introduces the energy asymmetry.27 For given gate voltages, Vt>0 and Vb<0, by utilizing Gauss' law, we derive the total carrier density induced on the graphene layers as

n=nt+nb=εrε0qδE=εrε0qttVt+εrε0qtbVb·
(18)

Here, nt and nb are the carrier densities induced by the top and bottom gates, respectively, and q is the magnitude of the electron charge. We note that for the negative charged carriers, n is positive. The Fermi level is related to the total carrier density by the following equation:

n=αnelectronsnholes=α0ρα(E)f(E)dE0ρα(E)[1f(E)]dE,
(19)

where ρα(E) is the density of states of the αth subband and is defined as follows:43 

ρα(E)=4d2k(2π)2δ(Eϵα(k))·
(20)

We have included the valley and spin degeneracies using a factor of four in the above equation. Using Eqs. (18) and (19), we find the Fermi level as a function of the top and bottom gate voltages. Figure 6 shows how the Fermi level of BLG varies with respect to Vt+Vb=ttδE for a fixed Δ. Here, we have assumed that top and bottom oxides are made of SiO2 with permittivity εr=3.9 and that they have equal thicknesses tt=tb=300 nm. As we can see from the plots, the Fermi level of monolayer and AB stacked BLG show a higher sensitivity to the gate voltage than the Fermi level of AA stacked BLG. Also, at high gate voltages, different energy asymmetries lead to different curves. The inset illustrates the curves for gate voltages close to zero. The average electric field induced by the gate contacts (E¯=(VtVb)/2tt) polarizes the graphene layers, which leads to an internal electric field given as Eint=qδn/2εr,gε0 with δn as the charge imbalance between the two layers and εr,g=1 as the relative permittivity of BLG. The internal electric field acts to cancel out the external electric field. This effect, known as electrostatic screening, plays a major role in realistic multilayer graphene devices and is carefully modeled in our analysis. Due to the finite electrostatic screening in BLG, the electron density on the top layer n1=(n+δn)/2 is not equal to the electron density on the bottom layer n2=(nδn)/2. We note that the charge conservation is satisfied, n1+n2=nt+nb. The charge imbalance δn=n1n2 is given as the relative weights of electron wavefunctions as follows:25 

δn=n1n2=2πα(|ψA1α(k)|2+|ψB1α(k)|2|ψA2α(k)|2|ψB1α(k)|2)kdk,
(21)

where (ψA1α(k),ψB2α(k),ψA2α(k),ψB1α(k)) represents the αth eigenvector of the Hamiltonian in Eqs. (6) and (8), and the integral is evaluated over all occupied states. Therefore, we can write the energy asymmetry as

Δ=q(E¯+Eint)d=(VtVb)qd/2ttδn(Δ)q2d2εrε0,
(22)

which is the same as the self-consistent Hartree approximation in Ref. 23. By setting Vt+Vb to a constant value and applying a differential voltage δVt,b to the gate contacts (Vt becomes Vt+δVt,b/2 and Vb becomes VbδVt,b/2), we fix the Fermi level and can independently study the impact of gate voltage on the energy asymmetry. Figure 6 illustrates the energy asymmetry as a function of differential gate voltage, δVt,b, for a variety of fixed Fermi levels. As we can see the energy asymmetry has a linear dependence on δVt,b. Moreover, for the given gate voltages, AA stacked BLG shows a higher energy asymmetry than AB stacked BLG. We note that the curves are also valid for larger differential gate voltages as long as the oxide field is below the dielectric breakdown strength of SiO2, which is larger than 1 V/nm.

FIG. 6.

Impact of gate on the Fermi level μ (top) and the energy asymmetry Δ (bottom) in BLG. The total gate voltage, Vt+Vb, determines μ and the differential gate voltage δVt,b determines Δ. The inset shows the curves near zero gate voltage. The energy asymmetry is calculated in the presence of imperfect electrostatic screening. In the bottom figure, the initial top gate voltage is denoted by Vt0. The initial bottom gate voltage, Vb0=10 V, is the same for all the cases.

FIG. 6.

Impact of gate on the Fermi level μ (top) and the energy asymmetry Δ (bottom) in BLG. The total gate voltage, Vt+Vb, determines μ and the differential gate voltage δVt,b determines Δ. The inset shows the curves near zero gate voltage. The energy asymmetry is calculated in the presence of imperfect electrostatic screening. In the bottom figure, the initial top gate voltage is denoted by Vt0. The initial bottom gate voltage, Vb0=10 V, is the same for all the cases.

Close modal

In this section, we discuss the features of propagating TM and TE plasmons, such as compression factor, mode confinement, and propagation length, in various graphene structures. From the perspective of graphene plasmonic devices (antennas and waveguides), it is desirable to simultaneously achieve confined (low transverse propagation length 1/2kx) waves with a low loss (high longitudinal propagation length 1/2kz) and a high compression factor kz/klight. We use the complex wavevector derived in Sec. II and combine it with the conductivity model from Sec. III B in the presence of gate effects (Sec. III C) for analysis in this section.

In Fig. 7, we plot the following plasmon characteristics versus the operating frequency: (i) relative longitudinal propagation length 1/2λkz normalized to the wavelength of light λ=2π/klight, (ii) relative phase constant kz/klight normalized to klight, and (iii) relative transverse propagation constant 1/2λkx normalized to λ. The results in this figure are obtained for AB and AA stacked BLG, and MLG in the terahertz and infrared regimes for μ=0.3 eV and Δ = 0 eV. The conductivity of MLG is derived using Green's functions approach described for BLG in Sec. III B. It can also be derived using the procedure outlined in Ref. 29. In the terahertz frequency band (left column of Fig. 7), all the graphene structures support TM mode propagation with similar propagation characteristics. As the frequency increases the plasmon waves become more compressed and more confined to the graphene sheet. While MLG has a relatively higher compression factor, its propagation length is smaller than that of BLG. In the infrared regime (right column of Fig. 7), there exists separate frequency bands that allow the propagation of both TM and TE modes in graphene structures. Generally, TE mode in both BLG and MLG has a phase constant (kz) very close to the phase constant of light (klight) and shows a very high propagation length in longitudinal direction but they are very loosely bounded to the graphene sheet. Conversely, TM mode mostly shows a high compression factor, but very low propagation length which is two orders of magnitude lower than the wavelength of light. However, they are more confined in the transverse direction than TE modes and so are better guided by the graphene structure. The trade-off between the highly confined TM mode and the highly propagating TE mode is a general characteristic of plasmons and has been reported previously.20,21,44 The reason is that high confinement to the graphene sheet implies that a large amount of plasmon wave energy is confined to the sheet which has a higher loss than the surrounding dielectric medium.

FIG. 7.

Relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx as a function of frequency for AB stacked BLG (a) and (b), AA stacked BLG (c) and (d), and MLG (e) and (f) in the terahertz (left column) and infrared (right column) regimes. Here, we assume that the Fermi level is μ=0.3 eV, energy asymmetry is Δ = 0 eV. The wavevectors are calculated using Eqs. (4) and (5).

FIG. 7.

Relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx as a function of frequency for AB stacked BLG (a) and (b), AA stacked BLG (c) and (d), and MLG (e) and (f) in the terahertz (left column) and infrared (right column) regimes. Here, we assume that the Fermi level is μ=0.3 eV, energy asymmetry is Δ = 0 eV. The wavevectors are calculated using Eqs. (4) and (5).

Close modal

To study the impact of chemical potential, we plot the various plasmon characteristics in graphene structures as a function of the chemical potential μ for a constant frequency in the terahertz band and a constant frequency in the infrared band in Fig. 8. The gate voltages corresponding to the Fermi levels are labeled on the top axis. We note that very high Fermi levels are not achievable in practice by mere gating since it would require electric fields higher than the break down field of SiO2. In our device setup, oxide voltages cannot exceed 300 V; however, to obtain high Fermi levels, one can utilize chemical doping which would basically act as a gate with a constant induced carrier density. For the 1 THz frequency (left column of Fig. 8), it is evident that as the Fermi level increases, AB stacked BLG and MLG show an increase in both longitudinal and transverse propagation lengths and a decrease in compression factor. However, plasmons in AA stacked BLG show little sensitivity to changes in the Fermi level. Moreover, similar to the trade-off between TM and TE modes, we observe a trade-off between compression factor and propagation length of TM mode plasmons as we change the Fermi level. At infrared frequencies (right column of Fig. 8), by changing the Fermi level, we can switch the mode of propagation between the highly confined TM mode and the loosely bounded TE mode.

FIG. 8.

Impact of Fermi level μ on the relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx for AB stacked BLG (a) and (b), AA stacked BLG (c) and (d), and MLG (e) and (f) for a given frequency in the terahertz band (left column) and infrared band (right column). The total gate voltages corresponding to each Fermi level are labeled on the top axis. The wavevectors are calculated using Eqs. (4) and (5).

FIG. 8.

Impact of Fermi level μ on the relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx for AB stacked BLG (a) and (b), AA stacked BLG (c) and (d), and MLG (e) and (f) for a given frequency in the terahertz band (left column) and infrared band (right column). The total gate voltages corresponding to each Fermi level are labeled on the top axis. The wavevectors are calculated using Eqs. (4) and (5).

Close modal

In Sec. III B, we discussed that in addition to the Fermi level, the energy asymmetry in BLG can be utilized as an extra degree of freedom to tune its optical conductivity. To study the impact of energy asymmetry, we plot the propagation characteristics of plasmons in BLG versus frequency for two values of energy asymmetry (Δ = 0 and 0.8 eV) and a fixed Fermi level in Fig. 9. The Fermi level needs to be sufficiently far from the charge neutrality point so that there are enough carriers to enable plasmon oscillations. Specifically, due to the induced bandgap in AB stacked BLG, there exists a plasmon cut-off region (μ<0.1 eV) where the carrier concentration is not sufficient for plasmon oscillations.9,11,12 Here, the Fermi level is fixed at 0.3 eV. Each curve in Figs. 9(a) and 9(b) corresponds to a different value of scattering rate for both Δ = 0 eV (solid) and Δ=0.8 eV (dashed). As we see from the figure, scattering rate has an enormous impact on plasmon waves decay. Hence, the ability to produce high quality samples with low impurity and defect concentrations becomes an essential challenge of building high performance graphene based plasmonic devices.

FIG. 9.

Impact of Energy asymmetry Δ and scattering rate Γ on the longitudinal propagation length 1/2λzkz (a) and (b) and relative phase constant kz/klight (c) and (d) for AB and AA stacked BLG. The propagation length can be improved by increasing Δ but at the expense of compression factor. The differential gate voltage required to create Δ is denoted as δVt,b. Here, we assume that the Fermi level is μ=0.3 eV. The wavevectors are calculated using Eqs. (4) and (5).

FIG. 9.

Impact of Energy asymmetry Δ and scattering rate Γ on the longitudinal propagation length 1/2λzkz (a) and (b) and relative phase constant kz/klight (c) and (d) for AB and AA stacked BLG. The propagation length can be improved by increasing Δ but at the expense of compression factor. The differential gate voltage required to create Δ is denoted as δVt,b. Here, we assume that the Fermi level is μ=0.3 eV. The wavevectors are calculated using Eqs. (4) and (5).

Close modal

Interestingly, the propagation length of plasmons in AB and AA stacked BLG can be improved by increasing Δ. Plasmons in AB stacked BLG can propagate 2 to 5 times longer, depending on the frequency band and the scattering rate. Similarly, plasmons in AA stacked BLG can propagate 1.5 to 100 times longer. The improvement in propagation length is mainly due to the reduction in absorption (σ(ω)) caused by suppression of interband transitions such as second valance subband to second conduction subband (ϵ,2ϵ+,2) and first conduction subband to second conduction subband (ϵ+,1ϵ+,2). This a significant benefit of BLG as it will enable low energy and low leakage plasmonic antennas and waveguides for THz and low infrared bands communication. However, these improvements come at the expense of smaller compression factor [Figs. 9(c) and 9(d)] as mentioned earlier.

Finally, we illustrate how the plasmon characteristics are affected by the medium dielectric. Figure 10 illustrates the plasmon characteristics as a function of the relative permittivity of dielectric medium for TM and TE modes of propagation. As we can see, while increasing dielectric constant lowers the propagation length of the TM mode, it increases the propagation length of the TE mode. This effect can also be seen directly from Eqs. (4) and (5) where the dependence on dielectric constant appears in η. Moreover, higher compression factors can be achieved for TM mode by increasing dielectric constant.

FIG. 10.

Relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx as a function of dielectric constant εr for AB stacked BLG for TM (1 THz—solid) and TE (80 THz—dashed) modes. Here, we assume that the Fermi level is μ=0.3 eV, Γ = 20 meV, and T = 300 K. The wavevectors are calculated using Eqs. (4) and (5).

FIG. 10.

Relative phase constant kz/klight, longitudinal propagation length 1/2λzkz, and transverse propagation length 1/2λxkx as a function of dielectric constant εr for AB stacked BLG for TM (1 THz—solid) and TE (80 THz—dashed) modes. Here, we assume that the Fermi level is μ=0.3 eV, Γ = 20 meV, and T = 300 K. The wavevectors are calculated using Eqs. (4) and (5).

Close modal

The focus of this paper is to quantify the impact of material parameters on the tunability of plasmon propagation characteristics in bilayer graphene (BLG) structures. Toward this end, we consider a dual gate metal-oxide-graphene geometry to independently tweak the Fermi level and energy asymmetry in BLG through the gate terminal voltages. We numerically obtain the relationship between energy asymmetry and the applied gate voltage using the self-consistent Hartree approximation. The optical conductivity of BLG is modeled as a function of Fermi level and energy asymmetry using the Kubo formalism. The propagation of surface plasmons is studied numerically to quantify their propagation length, mode confinement, and compression factor in various graphene structures including AB stacked and AA stacked BLG, and monolayer graphene (MLG). Important trade-offs between propagation length and mode confinement that must be considered for device design are discussed. The impact of scattering rate on the propagation length of plasmons is also quantified. Finally, we demonstrate that the propagation length of TM plasmons in BLG can be improved by increasing the energy asymmetry. Such an advantage is absent in MLG. Overall, the results presented in this work provide insight into the design and optimization of plasmonic structures utilizing graphene as their platform.

The authors acknowledge the funding support of National Science Foundation through the Grant No. CCF-1565656.

In general, the scattering rate depends on the bandstructure and the Fermi velocity. A larger energy asymmetry would reduce the Fermi velocity of electrons. However, Fermi velocity is not the only factor that affects scattering rate and correspondingly the various figures of merit of plasmons. Factors including polarizability and screening, which depend on Δ, are important to consider, particularly in the case of charged impurity scatterers. These factors could outweigh the reduction in Fermi velocity and even decrease the scattering rate. To study the nontrivial dependence of scattering rate on Δ, we numerically calculate the average scattering rate, Γ=1/2τ, where τ is the average carrier relaxation time and is given as45 

τ=s,αρs,α(E)τs,α(E)(df(E)dE)dEs,αρs,α(E)(df(E)dE)dE·
(A1)

Here, τs,α(E) is the energy-dependent relaxation time corresponding to the band s and subband α. Density of states is denoted as ρs,α(E), and f(E) is the Fermi-Dirac distribution function. The k-dependent relaxation time τs,α(k)=τs,α(E)|E=ϵs,α(k) is written using Boltzmann approximation as follows:46 

1τs,α(k)=2πn0d2k(2π)2|Vs,α(k,k)|2×(1cos(θk,k))δ(ϵ(k)ϵ(k)),
(A2)

where n0 is the density of scatterers, θk,k is the angle between the incident wavevector k and the scattered wavevector k, and Vs,α(k,k) is the matrix element of the disorder potential. Vs,α(k,k) depends on the scattering mechanism and the electron wavefunctions ψs,α(k) and is written as

Vs,α(k,k)=ψs,α(k)ψs,α(k)V(q).
(A3)

Here, q=kk and V(q)=V(r)eiq·rdr is the Fourier transform of V(r), which is the scattering potential. We consider long-range Coulomb scattering and short-range defect scattering mechanisms. In the presence of screening, for charged impurities located at distance d from the BLG sheet we have V(q)=v(q)exp(qd)/ε(q), where v(q)=e2/2εrε0q is the bare coulomb potential. The dielectric function is ε(q)=1+v(q)Π(q), and Π(q) represents the polarizability of BLG. The details of calculation of polarizability are provided in Refs. 47 and 48. For defect scattering, V(r) is a delta function in space; therefore, its Fourier transform is constant and independent of the dielectric function and polarizability, i.e., V(q)=V0. Here, we assume that V0=10eV nm2 (Ref. 49). Figure 11 illustrates k-dependent scattering rates, evaluated numerically using Eq. (A2), for both coulomb and defect scattering mechanisms, and for two values of energy asymmetry (i.e., Δ = 0 and Δ=0.8 eV). Material-specific parameters used to obtain the results reported in Fig. 11 are provided in the figure caption. The figure highlights the opposite behaviors of Coulomb and defect scattering in response to introducing a large Δ. For Coulomb scattering, scattering rate decreases with an increase in momentum36 because of the increase in screening. Due to the absence of screening in the case of defect scattering, scattering rate increases with wavevector.

FIG. 11.

Coulomb and defect scattering rates as functions of wavevector for two values of energy asymmetry, i.e., Δ = 0 and Δ=0.8 eV evaluated using Eq. (A2). The concentration of Coulomb scatterers is ni=1×1016 m2 and the concentration of defects is nd=0.01×1016 m2. The dielectric media is SiO2 with relative permittivity of 3.9.

FIG. 11.

Coulomb and defect scattering rates as functions of wavevector for two values of energy asymmetry, i.e., Δ = 0 and Δ=0.8 eV evaluated using Eq. (A2). The concentration of Coulomb scatterers is ni=1×1016 m2 and the concentration of defects is nd=0.01×1016 m2. The dielectric media is SiO2 with relative permittivity of 3.9.

Close modal

To study the impact of energy asymmetry on the average scattering rate, we use the results from Fig. 11 in Eq. (A1). The average scattering rate for both Coulomb and defect scattering mechanisms versus Δ is plotted in Fig. 12. We see that as Δ increases from 0 to 0.8 eV, the average scattering rate increases by 70% for the case of defect scatterers, whereas it decreases by 30% for the case of Coulomb scatterers. The reason for the unexpected behavior of Coulomb scatterers is the nontrivial behavior of polarizability. Our calculations show that for μ=0.3 eV, introducing a large energy asymmetry will increase the polarizability of BLG and consequently its dielectric function. This is the effect that outweighs the reduction in Fermi velocity and leads to a smaller average scattering rate overall in a sample dominated by Coulomb scatterers.

FIG. 12.

Average scattering rate Γ versus energy asymmetry Δ, for Coulomb and defect scattering mechanisms, evaluated using Eq. (A1). Coulomb and defect scatterers show opposite behavior as Δ changes. Material parameters used for this plot are the same as those in Fig. 11.

FIG. 12.

Average scattering rate Γ versus energy asymmetry Δ, for Coulomb and defect scattering mechanisms, evaluated using Eq. (A1). Coulomb and defect scatterers show opposite behavior as Δ changes. Material parameters used for this plot are the same as those in Fig. 11.

Close modal

Average scattering rate determines the intraband optical absorption and, therefore, the propagation length of plasmons in BLG. Depending on the dominant scattering mechanism, the propagation length could be improved or degraded. In some experimental works, defect scattering is believed to be dominant,36,46 and in some other works, it is Coulomb scattering that limits carrier transport.50,51 For Coulomb-dominated samples, given the specific Fermi level of 0.3 eV, introducing a large Δ will reduce the average scattering rate and consequently increase the plasmon propagation length, while for defect-dominated samples the opposite effect will happen. To elucidate the situation, we plot the plasmon propagation length in Fig. 13 in the absence and presence of energy asymmetry for (i) constant scattering rate, (ii) Coulomb-dominated scattering, and (iii) defect-dominated scattering. For comparison, we choose disorder concentration in a way that for the gapless BLG, all the three cases have the same average scattering rate, i.e., Γ=20 meV (or equivalently 4.84×1012 1/s). The disorder concentrations are ni=9.57×1016 m2 for Coulomb-dominated scattering and nd=0.012×1016 m2 for defect-dominated scattering. Introducing an energy asymmetry of Δ=0.8 eV will increase the average scattering rate of defect-dominated case to Γd=34 meV, which will degrade the propagation length relative to the constant scattering case. On the other hands, it will decrease the average scattering rate of Coulomb-dominated case to Γd=14 meV, which will improve the propagation length relative to the constant scattering case.

FIG. 13.

Propagation length (Left) and compression factor (Right) versus frequency for both zero and finite energy asymmetry, i.e., Δ = 0 and Δ=0.8 eV. Three cases of (i) constant scattering rate, (ii) Coulomb-dominated scattering, and (iii) defect-dominated scattering are considered. For the gapless BLG (Δ = 0), Γ=20 meV is the same for all the three cases.

FIG. 13.

Propagation length (Left) and compression factor (Right) versus frequency for both zero and finite energy asymmetry, i.e., Δ = 0 and Δ=0.8 eV. Three cases of (i) constant scattering rate, (ii) Coulomb-dominated scattering, and (iii) defect-dominated scattering are considered. For the gapless BLG (Δ = 0), Γ=20 meV is the same for all the three cases.

Close modal
1.
A.
Grigorenko
,
M.
Polini
, and
K.
Novoselov
,
Nat. photonics
6
,
749
(
2012
).
2.
M.
Jablan
,
M.
Soljačić
, and
H.
Buljan
,
Proc. IEEE
101
,
1689
(
2013
).
3.
S.
Allen
, Jr.
,
D.
Tsui
, and
R.
Logan
,
Phys. Rev. Lett.
38
,
980
(
1977
).
4.
Y.
Liu
,
R.
Willis
,
K.
Emtsev
, and
T.
Seyller
,
Phys. Rev. B
78
,
201403
(
2008
).
6.
F. H.
Koppens
,
D. E.
Chang
, and
F. J.
Garcia de Abajo
,
Nano Lett.
11
,
3370
(
2011
).
7.
W.
Zhou
,
J.
Lee
,
J.
Nanda
,
S. T.
Pantelides
,
S. J.
Pennycook
, and
J.-C.
Idrobo
,
Nat. Nanotechnol.
7
,
161
(
2012
).
8.
S.
Mikhailov
and
K.
Ziegler
,
Phys. Rev. Lett.
99
,
016803
(
2007
).
9.
Z.
Fei
,
E.
Iwinski
,
G.
Ni
,
L.
Zhang
,
W.
Bao
,
A.
Rodin
,
Y.
Lee
,
M.
Wagner
,
M.
Liu
,
S.
Dai
 et al,
Nano Lett.
15
,
4973
(
2015
).
10.
S. G.
Menabde
,
D. R.
Mason
,
E. E.
Kornev
,
C.
Lee
, and
N.
Park
,
Sci. Rep.
6
,
21523
(
2016
).
11.
J.
Chen
,
M.
Badioli
,
P.
Alonso-González
,
S.
Thongrattanasiri
,
F.
Huth
,
J.
Osmond
,
M.
Spasenović
,
A.
Centeno
,
A.
Pesquera
,
P.
Godignon
 et al,
Nature
487
,
77
(
2012
).
12.
A.
Woessner
,
M. B.
Lundeberg
,
Y.
Gao
,
A.
Principi
,
P.
Alonso-González
,
M.
Carrega
,
K.
Watanabe
,
T.
Taniguchi
,
G.
Vignale
,
M.
Polini
 et al,
Nat. Mater.
14
,
421
(
2015
).
13.
L.
Ju
,
B.
Geng
,
J.
Horng
,
C.
Girit
,
M.
Martin
,
Z.
Hao
,
H. A.
Bechtel
,
X.
Liang
,
A.
Zettl
,
Y. R.
Shen
 et al,
Nat. Nanotechnol.
6
,
630
(
2011
).
14.
A.
Vakil
and
N.
Engheta
,
Science
332
,
1291
(
2011
).
15.
J. M.
Jornet
and
I. F.
Akyildiz
,
IEEE J. Sel. Areas Commun.
31
,
685
(
2013
).
16.
L.
Vicarelli
,
M.
Vitiello
,
D.
Coquillat
,
A.
Lombardo
,
A.
Ferrari
,
W.
Knap
,
M.
Polini
,
V.
Pellegrini
, and
A.
Tredicucci
,
Nat. Mater.
11
,
865
(
2012
).
17.
F.
Rana
,
IEEE Trans. Nanotechnol.
7
,
91
(
2008
).
18.
M.
Liu
,
X.
Yin
,
E.
Ulin-Avila
,
B.
Geng
,
T.
Zentgraf
,
L.
Ju
,
F.
Wang
, and
X.
Zhang
,
Nature
474
,
64
(
2011
).
19.
S.
Rakheja
and
P.
Sengupta
,
IEEE Trans. Nanotechnol.
15
,
113
(
2016
).
20.
M.
Jablan
,
H.
Buljan
, and
M.
Soljačić
,
Phys. Rev. B
80
,
245435
(
2009
).
21.
G. W.
Hanson
,
J. Appl. Phys.
103
,
064302
(
2008
).
22.
M.
Jablan
,
H.
Buljan
, and
M.
Soljačić
,
Opt. Express
19
,
11236
(
2011
).
23.
E. V.
Castro
,
N.
Peres
,
J. L.
Dos Santos
,
F.
Guinea
, and
A. C.
Neto
,
J. Phys.: Conf. Ser.
129
,
012002
(
2008
).
24.
25.
E.
McCann
,
D. S.
Abergel
, and
V. I.
Fal'ko
,
Eur. Phys. J.: Spec. Top.
148
,
91
(
2007
).
26.
T.
Ohta
,
A.
Bostwick
,
T.
Seyller
,
K.
Horn
, and
E.
Rotenberg
,
Science
313
,
951
(
2006
).
27.
Y.
Zhang
,
T.-T.
Tang
,
C.
Girit
,
Z.
Hao
,
M. C.
Martin
,
A.
Zettl
,
M. F.
Crommie
,
Y. R.
Shen
, and
F.
Wang
,
Nature
459
,
820
(
2009
).
28.
L.
Falkovsky
and
A.
Varlamov
,
Eur. Phys. J. B
56
,
281
(
2007
).
29.
L.
Falkovsky
and
S.
Pershoguba
,
Phys. Rev. B
76
,
153410
(
2007
).
30.
L.
Falkovsky
,
J. Phys.: Conf. Ser.
129
,
012004
(
2008
).
31.
J.
Horng
,
C.-F.
Chen
,
B.
Geng
,
C.
Girit
,
Y.
Zhang
,
Z.
Hao
,
H. A.
Bechtel
,
M.
Martin
,
A.
Zettl
,
M. F.
Crommie
 et al,
Phys. Rev. B
83
,
165113
(
2011
).
32.
A.
Rozhkov
,
A.
Sboychakov
,
A.
Rakhmanov
, and
F.
Nori
,
Phys. Rep.
648
,
1
(
2016
).
33.
E.
Nicol
and
J.
Carbotte
,
Phys. Rev. B
77
,
155409
(
2008
).
34.
L.
Zhang
,
Z.
Li
,
D. N.
Basov
,
M.
Fogler
,
Z.
Hao
, and
M. C.
Martin
,
Phys. Rev. B
78
,
235408
(
2008
).
35.
C. J.
Tabert
and
E. J.
Nicol
,
Phys. Rev. B
86
,
075439
(
2012
).
36.
M.
Monteverde
,
C.
Ojeda-Aristizabal
,
R.
Weil
,
K.
Bennaceur
,
M.
Ferrier
,
S.
Guéron
,
C.
Glattli
,
H.
Bouchiat
,
J.
Fuchs
, and
D.
Maslov
,
Phys. Rev. Lett.
104
,
126801
(
2010
).
37.
R. R.
Nair
,
P.
Blake
,
A. N.
Grigorenko
,
K. S.
Novoselov
,
T. J.
Booth
,
T.
Stauber
,
N. M.
Peres
, and
A. K.
Geim
,
Science
320
,
1308
(
2008
).
38.
J. B.
Oostinga
,
H. B.
Heersche
,
X.
Liu
,
A. F.
Morpurgo
, and
L. M.
Vandersypen
,
Nat. Mater.
7
,
151
(
2008
).
39.
E. V.
Castro
,
K.
Novoselov
,
S.
Morozov
,
N.
Peres
,
J. L.
Dos Santos
,
J.
Nilsson
,
F.
Guinea
,
A.
Geim
, and
A. C.
Neto
,
Phys. Rev. Lett.
99
,
216802
(
2007
).
40.
J.
Nilsson
,
A. C.
Neto
,
F.
Guinea
, and
N.
Peres
,
Phys. Rev. B
76
,
165416
(
2007
).
41.
H.
Min
,
B.
Sahu
,
S. K.
Banerjee
, and
A.
MacDonald
,
Phys. Rev. B
75
,
155115
(
2007
).
42.
A.
Avetisyan
,
B.
Partoens
, and
F.
Peeters
,
Phys. Rev. B
80
,
195401
(
2009
).
43.
N. W.
Ashcroft
and
D.
Mermin
,
Solid State Physics
(
Saunders
,
Philadelphia
,
1976
).
44.
S. A.
Maier
,
Plasmonics: Fundamentals and Applications
(
Springer Science & Business Media
,
2007
).
45.
E.
Hwang
and
S. D.
Sarma
,
Phys. Rev. B
77
,
115449
(
2008
).
46.
S. D.
Sarma
,
E.
Hwang
, and
E.
Rossi
,
Phys. Rev. B
81
,
161407
(
2010
).
47.
E.
Hwang
and
S. D.
Sarma
,
Phys. Rev. Lett.
101
,
156802
(
2008
).
48.
M.
Lv
and
S.
Wan
,
Phys. Rev. B
81
,
195409
(
2010
).
49.
S. D.
Sarma
,
S.
Adam
,
E.
Hwang
, and
E.
Rossi
,
Rev. Mod. Phys.
83
,
407
(
2011
).
50.
W.
Zhu
,
V.
Perebeinos
,
M.
Freitag
, and
P.
Avouris
,
Phys. Rev. B
80
,
235402
(
2009
).
51.
S.
Xiao
,
J.-H.
Chen
,
S.
Adam
,
E. D.
Williams
, and
M. S.
Fuhrer
,
Phys. Rev. B
82
,
041406
(
2010
).