Scattering of surface plasmon-polaritons in a graphene multilayer photonic crystal with inhomogeneous doping
Yu. V. Bludov, N. M. R. Peres, G. Smirnov, M. I. Vasilevskiy
SScattering of surface plasmon-polaritons in a graphene multilayer photonic crystalwith inhomogeneous doping
Yu. V. Bludov, N. M. R. Peres, G. Smirnov, M. I. Vasilevskiy
Department of Physics and Center of Physics, University of Minho, P-4710-057, Braga, Portugal
The propagation of a surface plasmon-polariton along a stack of doped graphene sheets is con-sidered. This auxiliary problem is used to discuss: (i) the scattering of such a mode at an interfacebetween the stack and the vacuum; (ii) the scattering at an interface where there is a sudden changeof the electronic doping. The formalism is then extended to the barrier problem . In this system richphysics is found for the plasmonic mode, showing: total reflection, total transmission, Fabry-P´erotoscillations, and coupling to photonic modes.
PACS numbers: 81.05.ue,72.80.Vp,78.67.Wj
I. INTRODUCTION
Plasmonics deals with the excitation, manipulation,and utilization of surface plasmon-polaritons (SPPs),where the latter are hybridized excitations of radia-tion with the collective charge oscillations of an electrongas . In traditional noble-metal plasmonics the elec-tron gas is provided by the free electrons in the metal.Furthermore, SPPs are excited at the interface between ametal and a dielectric and propagate along the interfacewith exponential localization in the direction perpendic-ular to that of their motion.One central idea in plasmonics is to explore thesub-wavelength confinement of light to build plasmonicwaveguides that would propagate, at the same time,an electric signal and a highly confined electromagneticwave . A plasmonic circuitry would involve lenses, mir-rors, beam splitters, and the like. Therefore, the studyof scattering of plasmons by such structures arises.Clearly, the problem of scattering of plasmons is a sci-entific and technological one. The deep understandingof the scattering of plasmons is instrumental for build-ing new technologies. In traditional noble-metal plas-monics, the range of wavelengths where SPPs show sub-wavelength confinement is restricted to the interval span-ning the near infrared (near-IR) to the ultraviolet. In themid-infrared (mid-IR) to the terahertz (THz) spectralrange SPPs in structures with noble metals are essen-tially free radiation, therefore lacking the key advantageof sub-wavelength confinement. This makes them unsuit-able for plasmonics circuitry and sensing .From what has been said above it follows that newplasmonic materials able of showing sub-wavelength con-finement and spanning the frequency interval rangingfrom the THz to the mid-IR are necessary. This is partic-ularly relevant as important biomolecules exhibit uniquespectral signatures in this frequency range. Thus, thesensing capability arising with noble-metal plasmonics inthe near-IR to the ultraviolet could be extended to aregion that the traditional systems cannot cover. Suchpossibility would increase the application of plasmonicsfor sensing and security applications, such as detectionof pollutants, diagnosis of diseases, food control quality, and detection of plastic explosives.It is in the above context that graphene emerges asa promising plasmonic material . SPPs in grapheneexist in the THz to mid-IR range and show a high degreeof sub-wavelength localization, therefore circumventingthe mentioned limitations of noble-metal plasmonics. In-deed, it can be shown that the degree of localization ofplasmons in graphene is given by ζ G ∝ α (cid:126) c E F ( (cid:126) ω ) , (1)where α is the fine structure constant of atomic physics, c is the speed of light, E F is the Fermi energy of graphene,and ω is the frequency of the surface plasmon-polariton.Taking, as an example, a frequency of 150 THz (equiv-alent to the wavelength of λ =2 µ m for the radiationin vacuum, which corresponds to the edge of the mid-IRregion), and considering a typical Fermi energy of 0.5 eV(a value easily attainable by the electrostatic gating) weobtain for ζ G ∼ . µ m, that is λ ζ G ∼ , (2)which is a rather high degree of localization. Thisvalue yields highly intense and localized electromagneticfields. The above estimation highlights the potentialof graphene plasmonics in the THz to mid-IR spectralrange.Being a two-dimensional membrane, graphene isamenable for stacking. The idea is to build a pho-tonic crystal composed by several stacked sheets ofgraphene separated by dielectric layers; this struc-ture has been investigated both theoretically andexperimentally . In such structures charge carriers indifferent graphene layers are able to interact by meansof electromagnetic waves, which can be either propagat-ing or evanescent inside the dielectric. The latter caserefers to the area of plasmonics, where interaction be-tween the SPPs supported by each of the graphene lay-ers results in the formation of polaritonics bands .This fact allows for the existence of a number of inter-esting phenomena such as Bloch and Rabi oscilla-tions of SPPs as well as the formation of nonlinear self-localized wavepackets —lattice solitons . Moreover, a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r it was predicted that a plasmonic biosensor based on agraphene multilayer system shows a much higher sensi-tivity than its counterpart operating using a gold film .The propagation of bulk waves in a graphene stack ischaracterized by several phenomena typical for periodicstructures, like the presence of the omnidirectional low-frequency gap in the spectrum (which is not presentin the photonic crystals without graphene), extraordi-nary absorption decrease , and light pulse delay . Inpractice, these graphene multilayers can be used as ter-ahertz modulators , broadband polarizers , tunableBragg reflectors , and polarization splitters . Also it isinteresting that the graphene stacks exhibit the proper-ties of hyperbolic metamaterials . Finally, the studyof propagation of radiation in a disordered graphene stackwhen light impinges perpendicularly to the graphenesurface has recently been considered , showing thatgraphene can control Anderson localization of radiation.As mentioned above, eventually it will be necessary tobuild some kind of plasmonic circuitry where the problemof scattering of SPPs arises. In graphene, it is possible tocontrol the percentages of reflection and transmittanceof a surface plasmon-polariton by controlling the localvalue of the electronic density. Consider the simplestcase of a graphene sheet on a split gate. Each of the twoparts of the gate is subjected to different gate potentials,and this creates two zones in the material presenting twodifferent electronic concentrations. Assuming now that asurface plasmon-polariton is impinging on the border linedefined by the split gate, the amount of power reflectedand transmitted will be controlled by the difference inthe local electronic densities. The problem just describedhas already been discussed in the literature . Anotherinteresting question is the coupling of SPPs to photonicmodes. The idea can work in two ways: either a photonicmode will excite a surface plasmon-polariton in grapheneor a surface plasmon-polariton propagating on graphenewill radiate to free space as a photonic mode. Both casesfind relevant technological applications.As we will show in the bulk of the article, it is easier toachieve the interaction between SPP and photonic modesin the graphene stack, than in single graphene layer. Dueto the periodicity of the system, this interaction can bedirect (i.e. without using prisms). Also the scatteringof SPPs in a stack of graphene sheets has rich physics,including total reflection, total transmission, and Fabry-P´erot oscillations.The paper is organized as follows. In Sec.II we considerthe eigenvalues and the eigenfunctions of a graphene mul-tilayer photonic crystal (PC), which will serve as a basisfor the following sections. Sec.III is devoted to two prob-lems: (i) scattering of an incident polaritonic mode atthe interface between the graphene multilayer PC and ahomogeneous dielectric; (ii) scattering of an incident po-laritonic mode at the interface between two PCs, charac-terized by different Fermi energies of the graphene sheets.In Sec.IV we consider the scattering of a polaritonic modeon a double interface between PCs. II. EIGENMODES OF THE GRAPHENEMULTILAYER PHOTONIC CRYSTAL
In order to calculate the reflection of SPPs from the in-terface between two graphene multilayer PCs, it is neces-sary to know the spectrum of the electromagnetic wavesin PC. In the present section we consider an auxiliaryproblem of eigenmodes in a multilayer graphene stackcomposed of an infinite number of single graphene layerswith equal Fermi energy E F [see Fig.1(g)]. We supposethat graphene layers are arranged at equal distances d from each other at planes z = md , m ∈ ( −∞ , ∞ ) andare embedded into a uniform dielectric medium with adielectric constant ε . If the electromagnetic field is uni-form along the direction y ( ∂/∂y ≡ p -polarized waves, whose magnetic field is perpendicularto the plane of incidence ( xz ). Such a wave possessesthe electromagnetic field components (cid:126)E = { E x , , E z } ,(cid:126)H = { , H y , } , and is described by the Maxwell equa-tions ∂E x ∂z − ∂E z ∂x = iκH y , (3) ∂H y ∂z = iκεE x , ∂H y ∂x = − iκεE z . (4)Here we assumed the temporal dependence of the elec-tromagnetic field (cid:126)E, (cid:126)H in the form of exp( − iωt ), where ω is the cyclic frequency, κ = ω/c and c is the speed oflight in vacuum. The electromagnetic properties of thegraphene layer are determined by its dynamical conduc-tivity σ g ( ω ), whose form can be found, e.g., in Ref. .In order to find the dispersion relation of the graphenemultilayer PC, the electromagnetic fields should be con-sidered separately in each layer between the adjacentgraphene sheets at planes z = md and z = ( m + 1) d .The solutions of the Maxwell equations (3)–(4) at thespatial domain md ≤ z ≤ ( m + 1) d can be represented as H y ( x, z ) = { H m, + exp [ ik z ( z − md )] + (5) H m, − exp [ − ik z ( z − md )] } exp( ik x x ) ,E x ( x, z ) = k z κε { H m, + exp [ ik z ( z − md )] − (6) H m, − exp [ − ik z ( z − md )] } exp( ik x x ) E z ( x, z ) = − k x κε { H m, + exp [ ik z ( z − md )] + (7) H m, − exp [ − ik z ( z − md )] } exp( ik x x ) . Here k x is the in-plane component of the wavevector(parallel to the graphene sheets), k z = (cid:0) κ ε − k x (cid:1) / , H m, ± are the amplitudes of the forward (sign ”+”) orbackward (sign ”–”) propagating waves. By matchingboundary conditions at z = md [continuity of the tangen-tial component of the electric field across the graphene E x ( x, md + 0) = E x ( x, md − H y ( x, md +0) − H y ( x, md −
0) = − (4 π/c ) j x = − (4 π/c ) σ g E x ( x, md )], one can findthat amplitudes H m, ± can be related to H m − , ± as (cid:18) H m, + H m, − (cid:19) = ˆ M (cid:18) H m − , + H m − , − (cid:19) , (8)where the matrix ˆ M reads asˆ M = (cid:18) exp ( ik z d ) [1 − i Λ k z ] i Λ k z exp ( − ik z d ) − i Λ k z exp ( ik z d ) exp ( − ik z d ) [1 + i Λ k z ] (cid:19) with Λ = 2 πσ g / ( iωε ). Since the considered structure isperiodic, it is possible to use the Bloch theorem, whichdetermines the proportionality between field amplitudesin the adjacent periods through the Bloch wavevector q : H m − , ± = exp ( − iqd ) H m, ± , (9)After substitution of this relation into Eqs.(8), the solv-ability condition of the resulting linear equations requires(here ˆ I is the unit matrix)Det (cid:12)(cid:12)(cid:12) ˆ M − exp ( iqd ) ˆ I (cid:12)(cid:12)(cid:12) = 0 , (10)which results into the dispersion relation for the p -wavein graphene multilayer PCcos ( qd ) − cos ( k z d ) − Λ k z sin ( k z d ) = 0 . (11)Similar expressions for the dispersion relation were ob-tained in Ref. . The dispersion relation Eq.(11) forfixed ω and q possesses an infinite number of solutionsfor k x . Further in the paper we will prescribe an index n ≥ H ( n ) m, + H ( n ) m, − = − exp (cid:16) − ik ( n ) z d (cid:17) Λ k ( n ) z exp (cid:16) ik ( n ) z d (cid:17) (cid:104) − Λ k ( n ) z (cid:105) − exp ( iqd ) == exp ( iqd ) − exp (cid:16) − ik ( n ) z d (cid:17) exp ( iqd ) − exp (cid:16) ik ( n ) z d (cid:17) , which allows to represent these amplitudes as H ( n ) m, + = H ( n ) ( q, ω ) exp ( iqd ) − exp (cid:16) − ik ( n ) z d (cid:17) √ A ( n ) exp ( iqmd ) , (12) H ( n ) m, − = H ( n ) ( q, ω ) exp ( iqd ) − exp (cid:16) ik ( n ) z d (cid:17) √ A ( n ) exp ( iqmd ) , (13)where H ( n ) ( q, ω ) is the magnetic field amplitude and A ( n ) is a normalization factor. Notice that amplitudes being represented in this form also satisfy Bloch con-dition (9). Substituting Eqs.(12) and (13) into (5), weobtain the expression for the component of the electro-magnetic field at spatial domain md ≤ z ≤ ( m + 1) d inthe form H ( n ) y, ± ( x, z || q, ω ) = H ( n ) ± ( q, ω ) × ψ ( n ) ( z || q, ω ) exp( ± ik ( n ) x x ) , (14) E ( n ) x, ± ( x, z || q, ω ) = H ( n ) ± ( q, ω ) iκε × ∂ψ ( n ) ( z || q, ω ) ∂z exp( ± ik ( n ) x x ) , (15) E ( n ) z, ± ( x, z || q, ω ) = ∓H ( n ) ± ( q, ω ) k ( n ) x κε × ψ ( n ) ( z || q, ω ) exp( ± ik ( n ) x x ) , (16)where ψ ( n ) ( z || q, ω ) = (cid:110) exp [ iqd ] cos (cid:104) k ( n ) z ( z − md ) (cid:105) (17) − cos (cid:104) k ( n ) z ( md + d − z ) (cid:105)(cid:111) exp [ iqmd ] √ A ( n ) , is a dimensionless spatial profile function. Here the nor-malization factor A ( n ) = 1 − cos ( qd ) cos (cid:16) k ( n ) z d (cid:17) +cos (cid:16) k ( n ) z d (cid:17) − cos ( qd ) k ( n ) z d sin (cid:16) k ( n ) z d (cid:17) is chosen to satisfy the condition1 d ˆ md + dmd (cid:12)(cid:12)(cid:12) ψ ( n ) ( z || q, ω ) (cid:12)(cid:12)(cid:12) dz = 1 . Also we take into account that all eigenmodes of the PCcan be either forward- or backward propagating: this factis stressed in Eqs.(14)–(16) by adding the signs “+” or“-” before the x -component of wavevector k ( n ) x as well asby the subscript in the amplitude H ( n ) ± .Before considering the dispersion properties in de-tail, it should be noticed that further in the paper weneglect collisional losses and take the relaxation ratein graphene Γ = 0, i.e., the graphene conductivity issupposed to be purely imaginary. This approximationsimplifies the analysis without changing quantitativelythe results. The spatial periodicity of the multilayergraphene PC gives rise to the band-gap structure [seeFig.1(a)]: the spectrum is composed of an infinite num-ber of bands [ n ≥
0, colored domains in Fig.1(a)], whoseboundaries are determined by the Bloch wavevector atthe center q = 0, or at edge q = π/d of Brillouinzone [bold solid and bold dashed lines in Fig.1(a), re-spectively]. One of the bands [ n = 0, depicted by redcolor in Fig.1(a)] is polaritonic, where electromagneticwaves are evanescent in z − direction (with purely imagi-nary k (0) z ), and propagating in x − direction (with purely FIG. 1. (a) Dispersion curves (squared x − component of wavevector k x versus squared frequency ω ) of graphene multilayerPC: colored domains correspond to the allowed bands, which boundaries are determined by Bloch wavevector at center ofBrillouin zone q = 0 (bold solid lines, like B–F), or at its edge q = π/d (bold dashed lines, like B’–F’); (b)–(f) Spatial profilesof the multilayer graphene PC’s eigenfunctions at ω = 5 meV. The parameters, which correspond to each of the eigenfunctionsin panels (b)–(f), are depicted by respective points B–F, B’–F’ in panel (a). In all panels other parameters are: d = 40 µ m, E F = 0 .
157 eV (which correspond to the gate voltage 25 V, applied to graphene on top of the 300 nm thickness SiO substrate), ε = 3 .
9, Γ = 0; (g) Schematic view of the graphene-based PC. real k (0) x ). The other bands ( n ≥
1) are photonic ones,four of which with n = 1 , ..., z − direction (withpurely real k ( n ) z ), while in x -direction they can be eitherpropagating [ (cid:16) k ( n ) x (cid:17) >
0, shaded by lighter colors inFig.1(a)], or evanescent [ (cid:16) k ( n ) x (cid:17) <
0, shaded by darkercolors in Fig.1(a)]. The latter are physically meaningfulonly in confined photonic crystals because they divergein either plus or minus infinity. Notice that, as followsfrom the inset in Fig.1(a), the polaritonic and the firstphotonic band touch each other at a cutoff frequency ω ∗ ≈ (4 αE F c/ (cid:126) εd ) / in the center of the Brillouin zone, q = 0 (see Appendix A for details). Below this frequency,for ω < ω ∗ , the edge of the polaritonic band at q = 0coincides with the light line k x = ω ε/c , while aboveit, for ω > ω ∗ , the boundary of the polaritonic band q = 0 detachs from the light line [which coincides withthe edge of the first photonic band, solid green line inFig.1(a)]. At high frequencies, when the localization ofthe polaritonic modes near each graphene layer is strong, and the SPPs, sustained by neighboring graphene lay-ers, are almost noninteracting. This fact gives rise to thesituation when polaritonic dispersion curves for differentBloch wavevectors q merge together [see e.g. points Band B’ in Fig.1(a)], although their spatial profiles remaindifferent, as it is evident from Fig.1(b).Some of the edges of the photonic bands are charac-terized by an interesting property: at q = π/d the edgesof the photonic bands with odd numbers, n = 2 l − ≤ l < ∞ ) possesses z -components of the wavevector k (2 l − z = (cid:26) ωc ε − (cid:16) k (2 l − x (cid:17) (cid:27) / = 2 l − d π. (18)Examples of such modes are C’ and E’ in Figs.1(c) and1(e). Corresponding expressions for the eigenfunctions(17) can be represented as (details can be found in Ap-pendix A) ψ (2 l − (cid:16) z || πd , ω (cid:17) = √ (cid:20) l − d πz (cid:21) . (19)When q = 0, the edges of the bands with even numbers n = 2 l (1 ≤ l < ∞ ) possess the same property, i.e. k (2 l ) z = (cid:26) ωc ε − (cid:16) k (2 l ) x (cid:17) (cid:27) / = 2 ld π, (20) ψ (2 l ) ( z || , ω ) = √ (cid:20) ld πz (cid:21) . (21)Examples of such modes are D and F in Figs.1(d) and1(f). It is interesting that expressions for solutions(18), (20) do not contain the graphene conductivity σ g (in other words, their spectrum does not dependupon the graphene’s Fermi energy E F ): z -components k (2 l ) z and k (2 l − z , specified by Eqs.(18) and (20), cancelconductivity-dependent term in the dispersion relation(11). This happens because the derivatives of their spa-tial profile functions vanish, ∂ψ ( n ) /∂z =0 at z = md [see(19), (21)]. According to (15), it implies zero tangentialcomponents of the electric field E ( n ) x, ± ( x, md || q ) = 0 atgraphene layers, clearly seen from the spatial profiles ofthe modes C’, D, E’, F in Figs.1(c)–1(f). The light line k x = ω √ ε/c mentioned above is also an eigenmode ofthe multilayer graphene PC with q = 0 [see mode C inFig.1(c)]. It is characterized by absence of the wavevec-tor’s z -component k ( n ) z = 0 and the spatial profile func-tion ψ ( n ) ( z || , ω ) = 1 (further in the paper this modewill be referred to as the light-line mode ). The bandindex for this mode is n = 0 in the frequency domain ω < ω ∗ (where this mode is the edge of polaritonic band)and n = 1 above the critical frequency ω > ω ∗ (wherethis mode is the edge of the first photonic band). As amatter of fact, this mode is a bulk electromagnetic wavepropagating in the x -direction and with zero longitudinalcomponent of the electric field E ( n ) x ( x, z || q, ω ) ≡
0, i.e. apurely transverse wave.
III. SCATTTERING OF POLARITONIC MODEFROM A SINGLE INTERFACE BETWEEN TWOGRAPHENE MULTILAYER PCS
Let us consider now an interface between two graphenemultilayer PCs (described in Sec.II). We suppose thatboth PCs possess the same period d , are embedded inthe same homogeneous dielectric medium with the di-electric permittivity ε , and graphene layers are arrangedin the same planes z = md along axis z . The only differ-ence between these graphene multilayer PCs is the Fermienergy, which is equal to E F in the first PC (occupy-ing the half-space x <
0) and to E F in the second one(occupying the half-space x > x -axis (thus coming from x = −∞ ) and impinges onthe aforementioned interface. The main objective of thepresent section is to describe the scattering of the polari-tonic mode on this interface, that is, we want to find the FIG. 2. A single [panel (a)] or double [panel (b)] interfacebetween two graphene multilayer PCs with different Fermienergies of graphene layers. transmission and reflection coefficients of the polaritonicmode as well as to determine which part of its energy istransferred to each of the photonic modes.In order to do this, we expand the electromagneticfield in series with respect to the graphene multilayerPC eigenmodes (14)–(16). So, the magnetic field and the z -component of the electric field in PC1 ( x <
0) can bewritten as H ,y ( x, z ) = H (0)1 , + ψ (0)1 ( z ) exp( ik (0)1 ,x x ) + ∞ (cid:88) n =0 H ( n )1 , − ψ ( n )1 ( z ) exp( − ik ( n )1 ,x x ) , (22) E ,z ( x, z ) = −H (0)1 , + k (0)1 ,x κε ψ (0)1 ( z ) exp( ik (0)1 ,x x ) + ∞ (cid:88) n =0 H ( n )1 , − k ( n )1 ,x κε ψ ( n )1 ( z ) exp( − ik ( n )1 ,x x ) . (23)In the same manner, the fields in PC2 can be expressedas H ,y ( x, z ) = ∞ (cid:88) l =0 H ( l )2 , + ψ ( l )2 ( z ) exp( ik ( l )2 ,x x ) , (24) E ,z ( x, z ) = − ∞ (cid:88) l =0 H ( l )2 , + k ( l )2 ,x κε ψ ( l )2 ( z ) exp( ik ( l )2 ,x x ) . (25)In Eqs.(22)–(25) we prescribe the PC index j = 1 , ψ ( n ) j , the x -component ofwavevector k ( n ) j,x , and the amplitude H ( n ) j, ± . Also we useband indices n and l , referring to PC1 and PC2, re-spectively, and drop the arguments q and ω (which areequal for PC1 and PC2) in functions and amplitudes forbrevity. Notice that Eqs.(22) and (23) contain only onemode (polaritonic one with index n = 0) propagatingtowards the interface in PC1 (referring to the above-mentioned incident wave with amplitude H (0)1 , + ) and afull set of modes propagating backward from the inter-face (corresponding to the reflected harmonics with am-plitudes H ( n )1 , − ). At the same time, Eqs.(24) and (25)contain only modes propagating in PC2 in the positivedirection of x -axis, which correspond to the transmittedmodes with amplitudes H ( m )2 , + .The next step is to apply the boundary conditions atthe interface x = 0 [continuity of tangential componentsof magnetic field H ,y (0 , z ) = H ,y (0 , z ), and electric field E ,z (0 , z ) = E ,z (0 , z )] and use the orthogonality of thespatial profile functions, ˆ md + dmd ψ ( n (cid:48) ) j ( z ) ψ ( n ) j ( z ) dz = δ n,n (cid:48) d, where the overbar denotes complex conjugation. Afterapplying this orthogonality conditions to Eqs.(22)–(25)we obtain the following equations for the amplitudes: δ n, H (0)1 , + + H ( n )1 , − = ∞ (cid:88) l =0 H ( l )2 , + Ψ l,n , (26) (cid:104) δ n, H (0)1 , + − H ( n )1 , − (cid:105) k ( n )1 ,x = ∞ (cid:88) l =0 k ( l )2 ,x H ( l )2 , + Ψ l,n , (27)where Ψ l,n = 1 d ˆ md + dmd ψ ( l )2 ( z ) ψ ( n )1 ( z ) dz. (28)It should be noticed that the obtained Eqs.(26) and (27)can also be applied to the case of an interface betweenthe PC and a homogeneous medium (formally in this case E F = 0). In this case the spatial profile functions willbe as follows: ψ ( l )2 ( z ) = exp (cid:16) ik ( l ) z z (cid:17) , (29) k ( l ) z = q + 2 ld π = (cid:114) κ ε − (cid:16) k ( l ) x (cid:17) . (30)In order to express the reflection and transmission co-efficients in terms of energy fluxes, we notice thatthe component of the Poynting vector along the direc-tion of propagation ( x -axis) for the n -th mode S ( n ) j, ± = − ( c/ π ) Re (cid:16) E ( n ) j,z, ± H ( n ) j,y, ± (cid:17) , after substututing the ex-plicit forms of the electromagnetic fields (22)–(25), canbe written as S ( n ) j, ± = ± c πκε Re (cid:16) k ( n ) j,x (cid:17) (cid:12)(cid:12)(cid:12) H ( n ) j, ± (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ψ ( n ) j ( z ) (cid:12)(cid:12)(cid:12) . In other words, if the mode is propagating (with purelyreal k ( n ) x ), it carries energy either in positive (sign “+”) or in negative (sign “-”) direction of x -axis. In contrast,evanescent modes (with purely imaginary k ( n ) x ) do notcarry any energy. Thus, we define the coefficients R n , T l as the integral characteristics R n = − ´ md + dmd S ( n )1 , − dz ´ md + dmd S (0)1 , + dz = Re (cid:16) k ( n )1 ,x (cid:17) (cid:12)(cid:12)(cid:12) H ( n )1 , − (cid:12)(cid:12)(cid:12) Re (cid:16) k (0)1 ,x (cid:17) (cid:12)(cid:12)(cid:12) H (0)1 , + (cid:12)(cid:12)(cid:12) , (31) T l = ´ md + dmd S ( l )2 , + dz ´ md + dmd S (0)1 , + dz = Re (cid:16) k ( l )2 ,x (cid:17) (cid:12)(cid:12)(cid:12) H ( l )2 , + (cid:12)(cid:12)(cid:12) Re (cid:16) k (0)1 ,x (cid:17) (cid:12)(cid:12)(cid:12) H (0)1 , + (cid:12)(cid:12)(cid:12) . (32)The coefficients R , T are the reflectance and the trans-mittance of the polaritonic mode, respectively, while theothers ( n (cid:54) =0) are normalized intensities of higher diffrac-tion orders in PC1 and PC2 (coefficients R n and T l , re-spectively).First, we will consider the scattering of an incidentpolaritonic mode on the interface between the graphenemultilayer PC and a homogeneous dielectric, schemati-cally depicted in Fig.3(i). As it was mentioned, in thecase q = 0 [left column in Fig.3] the polaritonic modeexists in the frequency range ω > ω ∗ ≈ . ω < ω ∗ is impossible [dark gray domains inin Figs.3(a)–3(d)]. In the frequency range ω ∗ < ω < πc/ √ εd ≈ . n = 0 [depicted by redline in Fig.3(a)] and the light-line mode n = 1 [depictedby green line in Fig.3(a)]. In the homogeneous dielec-tric in this frequency range there is only one propagatingmode with l = 0 [see Eq.(29)]. Notice the coincidence ofboth the shape and dispersion properties of modes with l = 0 in the homogeneous dielectric and the light-linemode in PC1. Nevertheless, due to the opposite parity of the spatial profile functions integral (28)Ψ , (0) ≡ R = 1) ofthe polaritonic mode in this frequency range, shown inFig.3(b).The narrow frequency range 15 . (cid:46) ω (cid:46) . n = 2 [depictedby blue line in Fig.3(a)] as well as two more modes in thehomogeneous dielectric, with l = − , n = 0, although not being ableto couple to the mode n = 2 in PC due to the oppo-site parity [see modes B in Fig.1(b) and D in Fig.1(d)],still can couple to the modes l = − , T − and T , shownin Fig.3(d), and a decrease of the polaritonic mode re- FIG. 3. (a,e) Dispersion curves [the same as in Fig.1(a)] of graphene multilayer PC for q = 0 [panel (a)], or for q = π/d [panel (e)]; (b)–(d), (f)–(h) Frequency dependence of the reflectance R n ( ω ) of polaritonic mode n = 0 [panels (b) and (f)],diffraction order intensities of PC modes n = 3 or n = 2 [correspondingly panels (c) and (g)] and those T l ( ω ) of homogeneousdielectric modes l = − l = − l = − l = 0 [panel (h)] for the case when the polaritonic mode withBloch wavevector q = 0 [left column, panels (b)–(d)] or q = π/d [central column, panels (f)–(h)] is scattered from the interfacebetween the graphene multilayer PC and homogeneous medium. The parameters of the structure under consideration are ε = 3 . E F = 0 .
157 eV, d = 40 µ m. Notice that in panel (d) transmission coefficients for modes l = − l = 1 are equal toeach other [as well as in panel (h) transmission coefficients for modes for modes l = − l = 0 are also equal]; (i) Schematicview of the interface between the graphene-based PC and a homogeneous medium. flectance [Fig.3(b)]. At higher frequencies, ω (cid:38) . n = 3 in the PC1 with the same par-ity as the incident polaritonic mode [compare modes B inFig.1(b) and E in Fig.1(e)], thus allowing their couplingand nonzero intensity of diffraction order R [Fig.3(c)].Nevertheless, the third diffraction order intensity R is adecreasing function of frequency. At the same time, thediffraction intensities T − , T and the reflectance R [seeFigs.3(d) and 3(b)] possess a maximum and a minimum,correspondingly, at the boundary of the white and theyellow domains, ω ≈ . q = π/d [central column in Fig.3] is characterized by the follow-ing interesting features. In this case the polaritonic mode n = 0 can exist at any frequency [red line in Figs.3(e)],thus total reflectance ( R ≡
1) takes place [Fig.3(f)] inthe frequency range ω (cid:46) . . (cid:46) ω (cid:46) . l = − , T − and T [Fig.3(h)]. At the same time, in the fre-quency range ω (cid:38) . n = 2 becomes possible [Fig.3(g)].The characteristics of the diffraction of the polaritonicmode at the interface between two graphene multilayerPCs are shown in Fig.4. As in the previous case, po-laritonic modes for q = 0 exist above the cutoff fre-quencies ω ∗ and ω ∗ in the PC1 and PC2, respectively[nonexistence domain ω < ω ∗ is depicted in Figs.4(a)and 4(b) by dark gray color]. The frequency range ω ∗ < ω (cid:46) . R = 1 [see Fig.4(b)] of the polaritonicmode from the interface between two PCs. In the fre-quency range ω (cid:38) . n = 1) below the cutoff frequency ω ∗ and polaritonic ( n = 0) above it. The last fact gives riseto a gradual decrease of the reflectance and an increaseof the transmittance [Fig.4(b)] of the incident polaritonicmode. At the same time, when q = π/d there is no cut-off frequency for the existence of polaritonic mode [seeFig.4(c)], which results in the nonzero transmittance atany frequency [Fig.4(d)]. It is interesting that for q = π/d (contrary to the case of q = 0) the reflectance (transmit- FIG. 4. (a), (c) Dispersion curves of graphene multilayer PCwith E F (solid lines) or E F (dashed lines) for q = 0 [panel(a)], or for q = π/d [panel (c)]; (b),(d) Frequency dependenceof the reflectance R ( ω ) and the transmittance T ( ω ) of po-laritonic mode with Bloch wavevector q = 0 [panel (b)] or q = π/d [panel (d)] when it is diffracted on the interface be-tween two graphene multilayer PCs with Fermi energies E F and E F . The parameters of the structure under considera-tion are: ε = 3 . E F = 0 .
157 eV, E F = 0 . d = 40 µ m. tance) is an increasing (decreasing) function of frequency[compare Figs.4(b) and 4(d)]. IV. SCATTERING OF POLARITONIC MODEFROM THE DOUBLE INTERFACE BETWEENTWO GRAPHENE MULTILAYER PCS
We now consider the situation [Fig.2(b)] when agraphene multilayer PC of finite width D (along the x -axis) and with graphene layer’s Fermi energy E F (fur-ther referred to as PC2) is cladded by two semi-infinitePCs which graphene layer are characterized by the Fermienergy E F (these two PCs will be referred to as PC1 andPC3). As in Sec.II, the incident wave is the polaritonicmode of PC1, propagating in the positive direction of the x -axis. Notice, that the only difference between Fig.2(a)and 2(b) is the presence of two interfaces in the last case.We note that the electromagnetic field in the region x < < x < D ,occupied by the PC2, the field components will have theform: H ,y ( x, z ) = (cid:88) l ψ ( l )2 ( z ) × (cid:104) H ( l )2 , + exp( ik ( l )2 ,x x ) + H ( l )2 , − exp( − ik ( l )2 ,x x ) (cid:105) , (33) E ,z ( x, z ) = − (cid:88) l k ( l )2 ,x κε ψ ( l )2 ( z ) × (cid:104) H ( l )2 , + exp( ik ( l )2 ,x x ) − H ( l )2 , − exp( − ik ( l )2 ,x x ) (cid:105) . (34) FIG. 5. Reflectance [panels (a) and (c)] and transmittance[panels (b) and (d)] versus frequency ω for the double interfacebetween two graphene multilayer PCs with parameters ε =3 . E F = 0 .
157 eV, E F = 0 . d = 40 µ m, q = 0 [panels(a) and (b)], or q = π/d [panels (c) and (d)], D = 45 µ m (bluelines), D = 15 µ m (green lines), D = 3 µ m (red lines). Due to the finite width D of the PC2, Eqs.(33) and (34)contain both forward- and backward propagating waves.Finally, in the PC3 (region x > D ) the electromagneticfield components are: H ,y ( x, z ) = (cid:88) n H ( n )3 , + ψ ( n )1 ( z ) × exp (cid:104) ik ( n )1 ,x ( x − D ) (cid:105) , (35) E ,z ( x, z ) = − (cid:88) n H ( n )3 , + k ( n )1 ,x κε ψ ( n )1 ( z ) × exp (cid:104) ik ( n )1 ,x ( x − D ) (cid:105) . (36)Applying the same boundary condition as in Sec.III, wehave: δ n, H (0)1 , + + H ( n )1 , − = (cid:88) l Ψ l,n (cid:104) H ( l )2 , + + H ( l )2 , − (cid:105) , (37) (cid:104) δ n, H (0)1 , + − H ( n )1 , − (cid:105) k ( n )1 ,x = (38) (cid:88) l k ( l )2 ,x Ψ l,n (cid:104) H ( l )2 , + − H ( l )2 , − (cid:105) , H ( n )3 , + = (cid:88) l Ψ l,n × (39) (cid:104) H ( l )2 , + exp( ik ( l )2 ,x D ) + H ( l )2 , − exp( − ik ( l )2 ,x D ) (cid:105) , H ( n )3 , + k ( n )1 ,x = (cid:88) l k ( l )2 ,x Ψ l,n × (40) (cid:104) H ( l )2 , + exp( ik ( l )2 ,x D ) − H ( l )2 , − exp( − ik ( l )2 ,x D ) (cid:105) . For this structure the coefficients T n are expressed as T n = ´ md + dmd S ( n )2 , + dz ´ md + dmd S (0)1 , + dz = Re (cid:16) k ( n )1 ,x (cid:17) (cid:12)(cid:12)(cid:12) H ( n )3 , + ( q ) (cid:12)(cid:12)(cid:12) Re (cid:16) k (0)1 ,x (cid:17) (cid:12)(cid:12)(cid:12) H (0)1 , + ( q ) (cid:12)(cid:12)(cid:12) , while the coefficients R n are the same as before [Eq.(31)].The frequency dependence of the reflectance and trans-mittance for this structure are shown in Fig.5. For q = 0 and large width of the PC2 [light gray domain inFigs.5(a) and 5(b)] the reflectance [blue line in Fig.5(a)]of the structure is slightly less than unity in the fre-quency range ω ∗ < ω (cid:46) . ω (compare withFig.4(b), where in the case of single interface this rangecorresponds to the total reflectance R = 1, T = 0).The reason for this is the tunneling through the region0 < x < D because there is no propagating mode in PC2.The tunneling rate (and, as a consequence, the transmit-tance T ) gradually grows when the PC2 width, D isdecreased, compare blue, green and red lines in Fig.5(b).Notice that for small D = 3 µ m [red line in Fig.5(b)] al-most the whole energy of the incident wave is transmit-ted via tunneling, thus giving T (cid:46)
1. In the frequencyrange ω (cid:39) . T = 1 (respectively, R = 0). The latter takes place atfrequencies, for which the width D matches an integernumber of polaritonic mode half-wavelengths, that is k (0)2 ,x D = lπ, l ≥ . (41)The same phenomenon also occurs when q = π/d [seeFigs.5(c) and 5(d)]. The frequencies at which full trans-mission takes place can be approximated using the Drudemodel for graphene’s conductivity. Correspondingly, theparameter Λ in the dispersion relation (11) can be ap-proximated by Λ = 2 αcE F / (cid:126) ω ε , while we also use theso-called non-retarded approximation for k (0)2 ,x = − ik (0) z .Under these assumptions the dispersion relation (11) canbe rewritten ascos ( qd ) − cosh (cid:16) k (0)2 ,x d (cid:17) + Λ k (0)2 ,x sinh (cid:16) k (0)2 ,x d (cid:17) = 0 , which along with Eq.(41) gives ω = 2 αcE F (cid:126) ε lπD sinh (cid:0) lπD d (cid:1) cosh (cid:0) lπD d (cid:1) − cos ( qd ) . (42)For the case D = 15 µ m Eq.(42) gives ω ≈ .
82 meV for q = 0 and ω ≈ .
81 meV for q = π/d , which qualitativelyagrees with Figs.5(b) and 5(d) [green line].The oscillations of the transmittance can be seen asan optical analogue of the well known effect of electronresonant tunneling in a double-barrier heterostructure(DBHS). The modulation depth of the electromagneticwave transmission is smaller than in the case of electronsbecause of the weaker confinement of the former (so thatdirect tunneling involving only evanescent waves in thePC2 is possible) and also because of the simultaneouspresence of several scattering channels.
V. CONCLUSIONS
We have analysed the scattering of surface plasmon-polaritons generated in a graphene multilayer photoniccrystal and propagating across its lateral surface or inter-face with another graphene multilayer PC with a differentband structure (controlled by the graphene’s Fermi en-ergy). In particular, for the low-frequency region whereonly the polaritonic mode (with imaginary z componentof the wavevector, k z ) is propagating in the directionalong the graphene sheets (while the modes with real k z have imaginary k x , i.e. are evanescent in this sense), wehave shown that this mode is totally reflected from the in-terface between the graphene multilayer PC and a homo-geneous dielectric. Nevertheless, in the higher-frequencyregion where photonic (i.e. propagating with real k z )modes are allowed, the partial transformation of the inci-dent polaritonic mode’s energy into that of other diffrac-tion orders (both in PC and in homogeneous dielectric)becomes possible, thus reducing the reflectance of the in-cident wave. Moreover, by virtue of the reciprocity prin-ciple this gives rise to the possibility to excite the PC’spolaritonic eigenmode by an external wave impinging onits edge. In-phase polaritonic mode (with q = 0) canalso be totally reflected from the interface between twoPCs, while for out-of-phase oscillations (with q (cid:54) = 0) thisscenario is impossible. It is also shown that the trans-mittance and the reflectance of a structure consisting ofthree photonic crystals and two interfaces between them(we asumed PC1=PC3 (cid:54) =PC2) exhibit Fabry-P´erot oscil-lations with a discrete set of frequencies for which totaltransmission of the polaritonic mode takes place. Thiseffect is similar to the resonant tunneling of electrons ina DBHS where the electric current shows a sharp peakas a function of bias. Here, in addition to the differencebetween the Fermi levels of the crystals PC1 and PC2,the electromagnetic wave transmission depends also onthe PC wavevector q . ACKNOWLEDGMENTS
We acknowledge support from the EC under theGraphene Flagship (Contract No. CNECT-ICT-649953).
Appendix A: Approximations for the eigenvaluesand eigenfunctions.
Let us consider first the asymptotical behaviour of thepolaritonic modes n = 0 at high frequencies, where themodes with different Bloch wavevector q merge. SinceSPPs are evanescent waves, it is natural to introduce adecay parameter in z -direction, p = − ik (0) z . In this casethe dispersion relation (11) can be rewritten as p = cosh ( pd ) − cos ( qd )Λ sinh ( pd ) . (A1)0Notice that in the limiting case d → ∞ the dispersion re-lation (A1) transforms into p = Λ − , i.e. into the disper-sion relation for a single graphene sheet. If we considerthe situation of d being large but finite, the single sheetdispersion relation can be used as zeroth approximationfor p . This approximation after being sibstituted into theright hand side of Eq.(A1), allows for an approximate re-spresentation of Eq.(A1) as p = (cid:114)(cid:16) k (0) x (cid:17) − κ ε = 1 − cos ( qd ) exp ( − d/ Λ)Λ . (A2)The above-mentioned zeroth approximation for p also canbe used for approximating the spatial profile function(17), namely: ψ (0) ( z || q, ω ) = (cid:8) exp [ iqd ] cosh (cid:2) Λ − ( z − md ) (cid:3) − cosh (cid:2) Λ − ( md + d − z ) (cid:3)(cid:9) exp [ iqmd ] √ A (0) ,A (0) = 1 − cos ( qd ) cosh (cid:0) Λ − d (cid:1) +Λ cosh (cid:0) Λ − d (cid:1) − cos ( qd ) d sinh (cid:0) Λ − d (cid:1) . Even more, since all the phenomena we are interestedin take place in the frequency range (cid:126) ω (cid:28) E F , wecan take into account only the Drude contribution tographene’s conductivity. Correspondingly, the parame-ter Λ can be approximated by Λ = 2 αcE F / (cid:126) ω ε .This fact can be used in the approximation of thetouching point (cutoff frequency) between the zeroth andfirst band at q = 0, mentioned in Sec.II. Taking into ac-ccount the smallness of k z in the touch points, one canexpand the dispersion relation (11) ascos ( qd ) − k z d ) − ( k z d ) − Λ k z (cid:34) k z d − ( k z d ) (cid:35) = 0 . At q = 0 we have k z = 12 d d − d − . Taking into account the above-mentioned approximationfor Λ we find that k z > k (0) z being purely real)or k z < k (1) z being purely imaginary) in the fre-quency ranges ω < ω ∗ = (4 αE F c/ (cid:126) dε ) / and ω > ω ∗ ,respectively. The last fact well corroborates with thedispersion curve behavior in the vicinity of the cutoff fre-quency ω ∗ , depicted in Fig.1(a).From Fig.1(a) it is evident that in the high-frequencyregion ω (cid:29) ω ∗ the gaps becomes narrow. We can use thisfact in order to approximate the eigenfunctions and theeigenvalues in this region. At the edge of the Brillouinzone ( q = π/d ) we represent z -component of the wavevec-tor in the band with even number n = 2 l (1 ≤ l < ∞ )as k (2 l ) z = k (2 l − z + ∆ (2 l ) with ∆ (2 l ) being small value( (cid:12)(cid:12) ∆ (2 l ) (cid:12)(cid:12) (cid:28) (cid:12)(cid:12)(cid:12) k (2 l − z (cid:12)(cid:12)(cid:12) ) and k (2 l − z is determined in Eq.(18).In this case the dispersion relation after the expansion ofsine and cosine with respect to ∆ (2 l ) can be representedas − (cid:0) ∆ (2 l ) d (cid:1) (cid:20) l − d π + ∆ (2 l ) (cid:21) ∆ (2 l ) d = 0 , from which it is possible to determine∆ (2 l ) = 2 l − d π d − k (2 l ) z = π l − d − . In the same manner it is possible to obtain an approxi-mate expression for the spatial profile function (17): ψ (2 l ) (cid:16) z || πd , ω (cid:17) = ( − m (cid:26) sin (cid:20) l − d πz (cid:21) − ∆ (2 l ) (cid:20) l − d πz (cid:21) [(2 m + 1) d − z ] − (A3) (cid:0) ∆ (2 l ) (cid:1) (cid:20) l − d πz (cid:21) × (cid:2) d − z − md ) ( z − md − d ) (cid:3)(cid:9) (cid:114) A (2 l ) , with A (2 l ) = 1 − ∆ (2 l ) d (2 l − π − (cid:32) − l − π (cid:33) (cid:16) ∆ (2 l ) d (cid:17) . The same formalism can be applied to the boundaries ofthe bands with an odd numbers n = 2 l + 1 (1 ≤ l < ∞ )at the center of Brillouin zone q = 0. Thus, we obtain∆ (2 l +1) = 2 ld π d − k (2 l +1) z = k (2 l ) z + ∆ (2 l +1) = π ld − . An approximate expression for the spatial profile function(17) can be written ψ (2 l +1) ( z ||
0) = (cid:26) sin (cid:20) ld πz (cid:21) − ∆ (2 l +1) (cid:20) ld πz (cid:21) [(2 m + 1) d − z ] − (A4) (cid:0) ∆ (2 l +1) (cid:1) (cid:20) ld πz (cid:21) × (cid:2) d − z − md ) ( z − md − d ) (cid:3)(cid:9) (cid:114) A (2 l +1) , with A (2 l +1) = 1 − ∆ (2 l +1) d lπ − (cid:18) − l π (cid:19) (cid:16) ∆ (2 l +1) d (cid:17) . S. A. Maier,
Plasmonics: Fundamentals and Applications (Springer, New York, 2007). T. Low and P. Avouris, ACS Nano , 1086 (2014). J. A. Dionne and H. A. Atwater, MRS Bulletin , 717(2012). M. I. Stockman, Phys. Today , 39 (2011). J. B. Khurgin and A. Boltasseva, MRS Bulletin , 768(2012). F. J. G. de Abajo, ACS Photonics (2014). Y. V. Bludov, N. M. R. Peres, and M. I. Vasilevskiy, Jour-nal of Optics , 114004 (2013). T. Stauber, J. Phys.: Condens. Matter , 123201 (2014). J. Chen, M. Badioli, P. Alonso-Gonz´alez, S. Thongrat-tanasiri, F. Huth, J. Osmond, M. Spasenovi´c, A. Centeno,A. Pesquera, P. Godignon, et al., Nature , 77 (2012). Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S.McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens,G. Dominguez, et al., Nature , 82 (2012). P. A. D. Gon¸calves and N. M. R. Peres,
An Introduction toGraphene Plasmonics (World Scientific, Singapore, 2016). F. H. L. Koppens, D. E. Chang, and F. J. G. de Abajo,Nano Lett. , 3370 (2011). O. L. Berman, V. S. Boyko, R. Y. Kezerashvili, A. A.Kolesnikov, and Y. E. Lozovik, Physics Letters A ,4784 (2010). O. L. Berman and R. Y. Kezerashvili, Journal of Physics:Condensed Matter , 015305 (2012). Z. Arefinia and A. Asgari, Physica E: Low-DimensionalSystems and Nanostructures , 34 (2013), ISSN13869477. C. Qin, B. Wang, H. Huang, H. Long, K. Wang, and P. Lu,Optics Express , 25324 (2014). A. J. Chaves, N. M. R. Peres, and F. A. Pinheiro, Phys.Rev. B , 195425 (2015). H. Hajian, A. Soltani-Vala, and M. Kalafi, Optics Com-munications , 149 (2013). S. A. El-Naggar, Optical and Quantum Electronics ,1627 (2015). F. Al-sheqefi and W. Belhadj, Superlattices and Mi-crostructures , 127 (2015). A. Madani and S. R. Entezar, Physica B: Condensed Mat-ter , 1 (2013). Y. H. Cheng, C. Chen, K. Y. Yu, and W. J. Hsueh, OpticsExpress , 28755 (2015). C. S. R. Kaipa, A. B. Yakovlev, G. W. Hanson, Y. R.Padooru, F. Medina, and F. Mesa, Physical Review B -Condensed Matter and Materials Physics , 4 (2012). Z. Xu, C. Chen, S. Q. Y. Wu, B. Wang, J. Teng, C. Zhang,and Q. Bao,
Graphene-polymer multilayer heterostructurefor terahertz metamaterials (2013). K. V. Sreekanth, S. Zeng, J. Shang, K.-T. Yong, and T. Yu,Scientific reports , 737 (2012). Y. Fan, Z. Wei, H. Li, H. Chen, and C. M. Soukoulis,Physical Review B , 1 (2013), 1311.7037. B. Wang, X. Zhang, F. J. Garc´ıa-Vidal, X. Yuan, andJ. Teng, Physical Review Letters , 1 (2012). Y. Fan, B. Wang, H. Huang, K. Wang, H. Long, and P. Lu,Optics Letters , 6827 (2014). F. Wang, C. Qin, B. Wang, S. Ke, H. Long, K. Wang, andP. Lu, Optics Express , 31136 (2015). Y. V. Bludov, D. A. Smirnova, Y. S. Kivshar, N. M. R.Peres, and M. I. Vasilevskiy, Physical Review B , 045424(2015). Z. Wang, B. Wang, H. Long, K. Wang, and P. Lu, OpticsExpress , 32679 (2015). K. V. Sreekanth, S. Zeng, K. T. Yong, and T. Yu, Sensorsand Actuators, B: Chemical , 424 (2013). J.-T. Liu, N.-H. Liu, H. Wang, T.-B. Wang, and X.-J. Li,Physica B: Condensed Matter , 66 (2014). B. Sensale-Rodriguez, R. Yan, M. M. Kelly, T. Fang,K. Tahy, W. S. Hwang, D. Jena, L. Liu, and H. G. Xing,Nature communications , 780 (2012). H. Yan, X. Li, B. Chandra, G. Tulevski, Y. Wu, M. Freitag,W. Zhu, P. Avouris, and F. Xia, Nature Nanotechnology , 330 (2012). I. V. Iorsh, I. S. Mukhin, I. V. Shadrivov, P. A. Belov, andY. S. Kivshar, Physical Review B , 075416 (2013). S. V. Zhukovsky, A. Andryieuski, J. E. Sipe, and A. V.Lavrinenko, Physical Review B , 155429 (2014). M. Cheng, P. Fu, M. Weng, X. Chen, X. Zeng, S. Feng,and R. Chen, Journal of Physics D: Applied Physics ,285105 (2015). A. Vakil and N. Engheta, Science (New York, N.Y.) ,1291 (2011). B. Rejaei and A. Khavasi, J. of Opt. , 075002 (2015). Y. V. Bludov, A. Ferreira, N. M. R. Peres, and M. I.Vasilevskiy, International Journal of Modern Physics B ,1341001 (2013). Here the parity is considered with respect to the middleof the dielectric slab between the graphene sheets, i.e. theplane z = ( m + 1 / d . L. L. Chang, L. Esaki, and R. Tsu, Appl. Phys. Lett.24