Probing Nonlocal Effects in Metals with Graphene Plasmons
Eduardo J. C. Dias, David Alcaraz Iranzo, P. A. D. Gonçalves, Yaser Hajati, Yuliy V. Bludov, Antti-Pekka Jauho, N. Asger Mortensen, Frank H. L. Koppens, N. M. R. Peres
PProbing Nonlocal Effects in Metals with Graphene Plasmons
Eduardo J. C. Dias, ∗ David Alcaraz Iranzo, P. A. D. Gon¸calves,
3, 4, 5
Yaser Hajati, Yuliy V. Bludov, Antti-Pekka Jauho,
7, 5
N. Asger Mortensen,
3, 8, 5
Frank H. L. Koppens,
2, 9 and N. M. R. Peres † Department of Physics and Center of Physics, and QuantaLab,University of Minho, PT–4710–057, Braga, Portugal ICFO - The Institute of Photonic Sciences, The Barcelona Instituteof Science and Technology, 08860 Castelldefels (Barcelona), Spain Center for Nano Optics, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark Department for Photonics Engineering, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark Center for Nanostructured Graphene, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark Department of Physics, Faculty of Science, Shahid Chamran University of Ahvaz, Ahvaz, Iran Department of Micro- and Nanotechnology, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark Danish Institute for Advanced Study, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark ICREA - Instituci´o Catalana de Recera i Estudis Avan¸cats, Barcelona, Spain (Dated: April 20, 2018)In this paper we analyze the effects of nonlocality on the optical properties of a system consistingof a thin metallic film separated from a graphene sheet by a hexagonal boron nitride (hBN) layer.We show that nonlocal effects in the metal have a strong impact on the spectrum of the surfaceplasmon-polaritons on graphene. If the graphene sheet is shaped into a grating, we show that theextinction curves can be used to shed light on the importance of nonlocal effects in metals. Therefore,graphene surface plasmons emerge as a tool for probing nonlocal effects in metallic nanostructures,including thin metallic films. As a byproduct of our study, we show that nonlocal effects lead tosmaller losses for the graphene plasmons than what is predicted by a local calculation. We showthat these effects can be very well mimicked using a local theory with an effective spacer thicknesslarger than its actual value.
I. INTRODUCTION
Nanoplasmonics is a field of optics dealing with theinteraction of electromagnetic radiation with metallicnanostructures and nanoparticles . Over the lastdecade, the characteristic size of plasmonic structures hasbeen steadily approaching the few-nanometer scale , withconcomitantly ultra-confined plasmonic modes . Themost common description of the electrodynamics govern-ing plasmonic systems typically ignores that the responseof a metallic nanostructure is controlled by a nonlocaldielectric tensor . Indeed, the general linear-responseexpression for the electric displacement vector reads D ( r , ω ) = ε (cid:90) d r (cid:48) ¯ ε ( r , r (cid:48) , ω ) E ( r (cid:48) , ω ) , (1)where D ( r , ω ) is the electric displacement vector, E ( r , ω )is the electric field, and ¯ ε ( r , r (cid:48) , ω ) is the nonlocal dielec-tric tensor. For a translational invariant system, we have¯ ε ( r , r (cid:48) , ω ) = ¯ ε ( r − r (cid:48) , ω ). Equation (1) embodies the state-ment that the electric displacement field at point r de-pends on the electric field at all points r (cid:48) .In general, and in particular for systems with transla-tional invariance, it is convenient to transform Eq. (1) tomomentum space, obtaining D ( k , ω ) = ε ¯ ε ( k , ω ) E ( k , ω ),where translational invariance has been assumed and k denotes the wave vector. The local response approxi-mation (LRA) is equivalent to neglecting the wave vec-tor dependence of the dielectric tensor by taking thelong-wavelength limit ( k → − k × [ k × E ( k , ω )] = − k [ k · E ( k , ω )] + k E ( k , ω )= ω c ¯ ε ( k , ω ) E ( k , ω ) . (2)Equation (2) has two types of solutions. The divergence-free solution ( k · E = 0) corresponds to the usual waveequation k E = ¯ ε ( ω/c ) E ; in this context, this solution isusually dubbed as the transverse mode. However, withinthe nonlocal framework, there is an additional curl-freesolution ( k × E = ) which, for a given frequency ω , isobtained when the condition ε ( k , ω ) = 0 is satisfied. Thissolution is traditionally referred to as the longitudinalmode, although one should bear in mind that this modeis also perpendicular to the direction of propagation ofthe field.The longitudinal solution is generally overlookedwithin the LRA, which may reveal itself as an inaccurateapproximation when either the size of the metallic nanos-tructures or the separation between two metallic surfacesfall below a couple of tens of nanometers. . Here,we consider a geometry similar to the latter, namely aconfiguration in which a graphene sheet is placed par-allel to and extremely close to a metal surface—seeFig. 1. The spectrum of the aforementioned structurescan be computed either analytically or numerically—depending of the geometry—and such calculations showthat the resonances associated with the excitation of sur- a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r grapheneairhBNmetalair ε ε ε ε IIIIIIIV shz xy
FIG. 1. Schematic representation of one type of layered het-erostructure investigated in this work. We assume the systemto be invariant in the x – and z –directions. face plasmon-polaritons (either localized or propagating)exhibit strong deviations from the ones obtained withinthe LRA . Additionally, the local-response calcula-tion predicts field-enhancements which are larger than inthe nonlocal case . This is particularity true whenwe have two nanoparticles in close proximity .Experimentally, one can investigate nonlocal effects us-ing diverse techniques, such as electron energy-loss spec-troscopy and far-field spectroscopy . Therefore, aproper theoretical description of nonlocal plasmonics re-quires the calculation of a suitable nonlocal dielectricfunction (or, equivalently, the nonlocal conductivity). Apossible approach is using ab-initio methods . How-ever, typical (albeit small) plasmonic structures involvea large number of atoms which render those methods dif-ficult to use in a routinely fashion, since they quickly be-come time-consuming and computationally demanding.Hence, it is natural to seek for an alternative approachto the calculation of the nonlocal dielectric function. Asearly as in the 1970s, it was observed that the plasmonicproperties of thin metallic films did not follow the pre-diction of the local-response theory . The need toreconcile theory and experiment required a description ofan electron gas subjected to external fields introduced byBloch and Jensen . This model of an electronic fluidcame to be known as the hydrodynamic model .This model combines Maxwell’s equations with Newton’ssecond law of motion for a charged particle supplementedwith a term taking into account the statistical pressurein the electron gas due to charge inhomogeneity .The revival of plasmonics unburied the hydrodynamicmodel and applied it to the optical response of metal-lic nanoparticles with great success . As alreadymentioned, both the position and width of the plasmonresonance, and the field enhancement measured experi-mentally do not agree with the approach using the LRA,whereas the hydrodynamic model is able to explain theexperimental results . For this reason, the hydrodynamicmodel has now become a popular approach for nonlocaloptics .More recently, the birth of 2D crystals introducedanother possibility for studying nonlocal effects in theoptical response of materials. In particular, the emer- gence of graphene plasmonics constitutes a new play-ground for studying these effects down to the one atomthick limit . The combination of graphene with two-dimensional (2D) insulators—such as hexagonal BoronNitride (hBN)—has allowed an unprecedented control ofthe distance that a 2D material can be placed in thevicinity of metallic films or nanostructures. In fact, theplacement of a graphene sheet at a distance of a fewnanometers away from a metal surface has recently beenexperimentally demonstrated , in a similar setting asdepicted in Fig. 1. This geometry provides a suitableway to obtain deep subwavelength confinements downto the atomic limit. This is possible due to the inter-play of graphene plasmons and the screening exerted bythe nearby metal film. This approach has revealed thatthe nonlocal properties of graphene’s conductivity haveto be included in a proper account of the experimentaldata . In graphene, it is more useful to describe thenonlocal response via the material’s nonlocal dynamicalconductivity, which within linear-response theory relatesthe surface current, K s ( r , ω ), and the in-plane electricfield, E (cid:107) ( r , ω ), via K s ( r , ω ) = (cid:90) d r (cid:48) ¯ σ ( r − r (cid:48) , ω ) E (cid:107) ( r (cid:48) , ω ) . (3)Notice that this statement is merely a reformulation ofEq. (1). Here, ¯ σ ( r , ω ) refers to the nonlocal (surface)conductivity tensor of graphene. This quantity can beprobed experimentally using graphene plasmons whenthe material is in the proximity of a metallic surface,and varying the graphene-metal distance, thereby re-trieving the Fourier transform of ¯ σ ( r , ω ), ¯ σ ( q , ω ) . Inthe absence of strain, graphene is isotropic and thus¯ σ ( q , ω ) = ¯ σ ( q, ω ), where q labels the in-plane wave vec-tor.If the graphene sheet is transformed into a grating (oran external grating is deposited on the graphene sheet),the properties of the plasmons depend sensitively on thelenght of the imposed period . Therefore, adjust-ing the grating period can be used to tune the grapheneplasmons’ frequency to the desired spectral range. Eithergraphene itself or a grating may be patterned . Whengraphene is placed in the vicinity of a metallic surface,e.g., as illustrated in Fig. 1, nonlocal effects arising fromthe metal’s response may influence indirectly the behav-ior of the fields in the region between graphene and themetal, namely through a larger penetration of the fieldinside the metal due to a less effective screening inducedby nonlocality. The quantification of the importance ofsuch nonlocal effects constitutes the goal of this work.Here, we investigate the influence of the nonlocal re-sponse of a metal on the optical response of systemsbased on graphene lying in close proximity to a metal-lic film or surface. In what follows, we consider boththe case in which the metal is constituted of Gold andof Titanium. We shall focus on the latter, as this metalexhibits a very strong nonlocal behavior (compared, forexample, with Gold). The nonlocal effects arising fromthe metal’s response are evaluated both on the plas-monic (near-field physics) and on the optical properties(far-field spectroscopy) of these structures. Through-out the manuscript, graphene’s conductivity is taken asbeing nonlocal . The aim of this work is therefore touse graphene plasmonics to probe nonlocal effects withinthe metal, in opposition to the study of nonlocality ingraphene performed in a recent publication .Specifically, we study the impact of nonlocal effects dueto the metal as a function of the graphene-metal separa-tion, and discuss its implications on the field distributionand plasmonic losses. Finally, we consider a prospec-tive system in which the graphene sheet is replaced by agraphene diffraction grating made of a periodic array ofgraphene ribbons. In this case, we compute the system’sresponse in the far-field and determine the influence ofthe metal’s nonlocality on the measured spectra. Wetake realistic parameters for the dielectric functions ofthe materials that constitute each system, which allowsus to compare our results directly to the experiments. II. MATHEMATICAL DETAILS
Throughout this work, we consider nonmagnetic ( µ =1) media. The system under study consists in a multi-layered heterostructure, with either dielectric or metal-lic layers stacked along the y –direction, as portrayed inFig. 1. The structure is further assumed to be infinite inthe xz –plane. We look for TM-modes with a harmonictime dependence in the form of e − i ωt , and thus the cor-responding magnetic field may be written as H ( r , t ) = H z ( x, y )e − i ωt ˆz . (4)The form of H z ( x, y ) follows from Maxwell’s equations,and, in general, is different of the regions defined in Fig. 1.As we consider nonlocal effects only in graphene and inthe metal, it is useful to distinguish the description of thefields in the dielectric(s) and metal regions. A. Dielectric Regions
In the dielectric media (source-free)the magnetic fieldobeys Helmholtz’s equation, ∇ H + εk H = 0, where k = ω/c . Assuming a magnetic field in the form ofEq. (4), the wave equation admits the following (trans-verse) solution: H z T ( x, y ) = (cid:0) C + e i k T y + C − e − i k T y (cid:1) e i qx , (5)where the in-plane wave vector, q (assumed to be alongthe x –direction without loss of generality), and the per-pendicular wave vector, k T , are related by k T = (cid:113) εk − q . (6) The ‘T’ subscript was introduced to make explicit thetransverse nature of these solutions. For uniaxial me-dia (characterized by different permittivities in the xz –plane, ε x , and in the y –direction, ε y ), such as hBN, q and k T are connect through an alternative condition,i.e., k T = (cid:112) ε x ω /c − q ε x /ε y . Maxwell’s equation E T = ( − i ωε ε ) − ∇ × H T enables us to write the cor-responding electric field components as E x T ( x, y ) = − k T ωε ε x (cid:0) C + e i k T y − C − e − i k T y (cid:1) e i qx , (7) E y T ( x, y ) = qωε ε y (cid:0) C + e i k T y + C − e − i k T y (cid:1) e i qx . (8)Naturally, if the dielectric is isotropic, then ε ≡ ε x = ε y . B. Metal Region
Within the metal, we assume a Drude dielectric func-tion ε m of the form ε m ( ω ) = ε ∞ − ω ω + i γω , (9)where ω p and γ are the plasma and damping frequencies,respectively, and ε ∞ a background permittivity to ac-count for interband polarization effects. We take nonlocaleffects in the metal within the framework of the hydrody-namic model—see Refs. 31,33 for details. When takingnonlocality into account, Amp`ere’s law becomes ∇ × H = − i ωε ε m [ E − ξ ∇ ( ∇ · E )] , (10)with the parameter ξ defined as ξ = β (cid:113) ω /ε ∞ − ω − i γω, (11)where β is a nonlocal parameter proportional to theFermi velocity v F of electrons in the metal. In this work,we will take the most usual definition β = (cid:112) / v F .Naturally, the local limit is recovered upon setting β = 0.As discussed above, Maxwell’s equations can be shownto admit two kinds of solutions : divergence-free andcurl-free fields. The former are the usual transversewaves that exist in the local regime.Within the LRA, theelectromagnetic fields in the metal are then describedsimilarly to the fields in the dielectric case, upon re-placing ε → ε m ( ω )—see Eqs. (5)–(8). On the otherhand, curl-free waves do not have an associated mag-netic field, as imposed by Faraday’s law. However, unlikethe local case, the electric field has a non-trivial solutiongiven by the vanishing of the term in square parenthe-sis figuring in Eq. (10), equivalent to the wave equation ∇ E − (1 /ξ ) E = 0. This equation describes longitudinalsolutions of the form E L = ( E x L ˆx + E y L ˆy )e − i ωt , with E x L ( x, y ) = (cid:0) D + e i k L y + D − e − i k L y (cid:1) e i qx , (12) E y L ( x, y ) = k L q (cid:0) D + e i k L y − D − e − i k L y (cid:1) e i qx , (13)with (longitudinal) wave vector k L = (cid:112) − ξ − − q = β − (cid:113) ω + i γω − ω /ε ∞ − q . (14)Furthermore, E x L and E y L are related by the curl-free con-dition, ∂E y L /∂x = ∂E x L /∂y = i qE x L .The general solution for the fields inside the metal istherefore H z = H z T and E x/y = E x/y T + E x/y L . C. Boundary Conditions
The coefficients which characterize the fields in eachlayer are determined by the boundary conditions (BCs).For bound modes (like surface plasmons) we have q > √ εω/c . This renders k imaginary, and the signal must bechosen judiciously to ensure that the fields satisfy Som-merfeld’s radiation condition. For an interface betweentwo dielectrics, the usual BCs apply, that is, the con-tinuity of the tangential component of the electric field( E x ) and the (dis)continuity of the magnetic field ( H z )in the (presence) absence of a surface current density atthe interface. The presence of a finite surface currentdensity is needed in order to describe a graphene sheet(or any other 2D material) placed at an interface be-tween the two media , and enters through Ohm’s law K s ( q, ω ) = σ ( q, ω ) E x ( q, ω ) ˆx ; see also Eq. (3).Note thathere σ ( q, ω ) entails both frequency and momentum de-pendencies in order to account for nonlocal effects ingraphene. We employ graphene’s nonlocal conductivityusing Mermin’s particle-conserving prescription, which isdetailed in Appendix C.Finally, for a dielectric/metal interface (withoutgraphene ), although the BCs described above remainvalid, these need to be augmented by an additional BC(due to the existence of a longitudinal mode within themetal). Typically, this additional BC dictates that thenormal component of the polarization vector vanishes atmetal’s surface, given by P = (i /ω ) ∇ × H − ε ε ∞ E . (15)Note that this polarization field refers to the one associ-ated with the free electrons in the surface of the metal,responsible for the transport of electric currents. There-fore, from a physical perspective, this condition merelyimposes that cannot exist currents flowing from the metalto the neighboring dielectric regions (e.g., the interfaceis a hard wall). In possession of this additional BC theamplitudes of the fields may now be determined; see Ap-pendices A and B for a more detailed account. III. NONLOCAL EFFECTS IN THEPLASMONIC PROPERTIES
We consider a system composed by a thin metallic layer(with thickness h ) with a graphene sheet lying at a dis-tance s from its surface. In the spacer region, i.e., be-tween graphene and the metal, we assume to have slabof hBN; however, our formalism is general and thus al-lows for the consideration of different dielectric media.Below the metal and above the graphene, we assume tohave air. For the sake of clarity, the system has been di-vided into four regions, I–IV, as depicted in Fig. 1. Thecalculation of the nonlocal plasmonic properties of thesystem follows the guidelines discussed in the previoussection, and can be consulted with a greater degree ofdetail in the Appendix A. In order to assess the influenceof the metal’s nonlocal effects in the plasmonic propertiesof the system, we compare our nonlocal results with thecorresponding predictions of the local-response theory. A. Nonlocal Effects in the Plasmon Dispersion
For the study of the metal’s nonlocal effects in thedispersion relation of graphene plasmons, we consid-ered here two different metals—Gold (Au) and Titanium(Ti)—described by the parameters presented in Table I.The dielectric function of hBN, on the other hand, isadopted from Refs. 36 (out-of-plane direction) and 41(in-plane direction).
TABLE I. Drude model parameters used for Titanium a (Ti)and Gold (Au). The superscripts indicate the correspondingreferences. Ti Au ω p [eV] 2.80 γ m [meV] 82.0 ε ∞ v F /c . . Note that the authors of Ref. 42 described Titanium using bothDrude and Lorentz terms, but, for the purpose of this work, weonly considered the Drude contribution, what nonethelessprovides a very good approximation.
The influence of the nonlocal effects arising from themetal’s response is investigated by computing the disper-sion relation, ω ( q ), of graphene plasmons and compar-ing the ensuing spectra obtained assuming a local- anda nonlocal-response. The outcome is presented in Fig. 2(white curves) for the allowed plasmonic modes in thestructure pictured in Fig. 1. Since there is dissipationin the system (both in the metal, hBN and graphene),either q or ω needs to be regarded as a complex quan-tity in order to fulfill the boundary conditions. In whatfollows, we have chosen q to be a real number and hence ω is complex. For the time being, we focus on the realpart of the frequency, Re { ω } , and denoting it by ω forsimplicity. The corresponding imaginary part, Im { ω } ,associated with the plasmonic losses, will be discussed ata later stage. FIG. 2. Dispersion relation (DR) of SPPs in anair/metal/hBN/graphene/air configuration, both for a local(dashed white) and nonlocal (solid white) metal, in the caseof Gold (top) and Titanium (bottom). The black dot-dashedline corresponds to the light dispersion in vacuum, whereasthe green dot-dashed line corresponds to the graphene’selectron-hole continuum boundary, ω = v F q , with v F ≈ c/ s = 1 nm, h = 10 nm, E F = 0 . It is apparent from Fig. 2 that the dispersion curvesare not continuous, but rather present two asymptotes around 750 and 1350 cm − ; this feature is shared in bothlocal and nonlocal frameworks. These frequencies corre-spond to the beginning of the hBN’s Restrahlen bands,where this material is hyperbolic , and are associatedwith the excitation of surface phonon-polaritons (whichin this case hybridize with graphene plasmons) .Moreover, Fig. 2 clearly demonstrates that nonlocaleffects (associated with the metal) impacts the plasmondispersion very differently depending on the metal: forAu, the nonlocal and local curves lay very close, whereasin the case of Ti, there is a significant blueshift of thedispersion of graphene plasmons due to the nonlocal ef-fects of the Ti (we stress that graphene is treated asbeing nonlocal in both cases). Quantitatively, the re-spective blueshifts are around 2% and 20%. We havedetermined that this difference in the magnitude of thenonlocal effects originates from the cumulative influenceof two main factors: the Ti’s larger nonlocal parameterwith respect to Au’s ( β Ti /β Au (cid:39) . ω p , Ti /ω p , Au (cid:39) . s of the spacer) is determinant, as Fig. 3(a)plainly shows. The figure shows the difference betweenthe local and nonlocal calculations increases as the thick-ness of the spacer region is decreased. In particular, whilefor s = 10 nm the difference is negligible, nonlocal ef-fects become clearly perceptible for s = 5 nm, and for s = 1 nm nonlocal effects become of paramount impor-tance for an accurate description of the system’s plas-monic response. This behavior is due to a greater spatialconfinement of the fields inside narrower spacers, whichis enhanced owing to the screening of graphene plasmonsby the metal, and in turn gives rise to an acoustic-likegraphene plasmons with very high momenta and thusmore susceptible to nonlocality.Figure 3(a) also shows that an increment of the spacerthickness induces a blueshift in the dispersion relationof the plasmonic modes.For that reason, a na¨ıve ap-proach to account for nonlocal effects while carryingout a local calculation is to consider an effective spacerthickness, s eff , larger than the actual value, s , in asimilar fashion to what has been proposed in earlierworks . Although this method can indeed be re-garded as somewhat na¨ıve version of quantum-correctedboundary conditions , it can mimic the proper nonlo-cal calculation as illustrated in Fig. 3(b). In particular,the nonlocal dispersion is can be reproduced using a localformalism with an effective parameter α = s eff /s between1 .
75 and 1 .
80. However, it should be stressed that the (a) (b)
FIG. 3. (a) Uppermost branch of the dispersion rela-tion ( ω > − ) of the plasmonic modes in anair/Ti/hBN/graphene/air configuration, for different valuesof the hBN’s slab thickness, s , calculated for a local (dashed)and nonlocal (solid) metal. (b) Dispersion relation in thesame configuration, calculated for local and nonlocal metalsfor different values of the effective parameter α = s eff /s , with s = 1 nm. The remaining parameters, for both panels, are: h = 10 nm, E F = 0 . value of the effective parameter highly depends both onthe considered materials (particularly the metal’s prop-erties) and on the configuration itself. For this reason,this procedure is generally hard to implement because theparticular value of α = s eff /s that reproduces the nonlo-cal effects is difficult to predict theoretically, which ren-ders this approach unsuitable. Instead, it can be merelyregarded as a fitting parameter when describing an ex-periment through local calculations. Before concluding the present section, we investigatehow nonlocal effects vary with the thickness of the metalfilm. To that end, in Fig. 4(a) we have plotted the plas-mon dispersion relation for an air/Ti/hBN/graphene/airstructure using three different values for the Titaniumthickness: 1, 10 and 100 nm. The plasmonic spectrumwas obtained both nonlocally (left panel) and locally(right panel). In both cases, the graphene is treatedas a nonlocal medium. The figure demonstrates thatthe 10 and 100 nm cases produce the same dispersion,whereas the 1 nm curve lies toward smaller frequencies(for the same q ). This results suggest that, above a cer-tain threshold thickness, the plasmonic properties of asystem with a finite-thickness metal are equivalent tothose of a configuration with a semi-infinite metal. Thisbehavior becomes particularly evident upon inspection ofFig. 4(b), where the plasmon wavevector, for a given fre-quency, is shown as a function of the metal’s thickness.Clearly, for h (cid:38) ∼
1. Comparison to the semi-infinite metal case
We have seen that for metal thicknesses h (cid:38) . This convenient be-cause it not only simplifies the analysis, but it also al-lows us to write a closed-form expression for the plas-mon dispersion in the heterostructure—now a dielec-tric/graphene/dielectric/metal configuration—, reading (cid:32) ε x κ (4)T + ε x κ (3)T + iσωε (cid:33) (cid:32) ε m κ (3)T ε x κ mT + δ nl (cid:33) == (cid:32) ε x κ (4)T − ε x κ (3)T + iσωε (cid:33) (cid:32) − ε m κ (3)T ε x κ mT − δ nl (cid:33) e − κ (3)T s , (16)where κ ( ν )T = (cid:112) q ε xν /ε yν − ε xν k for ν = { , } , κ mT = (cid:112) q − ε m k , and δ nl is a nonlocal correction term givenby δ nl = q κ mL κ mT ε m − ε ∞ ε ∞ , (17)with κ mL = (cid:113) q − (cid:0) ω + i γω − ω /ε ∞ (cid:1) /β . Indeed, ourcalculations demonstrate that Eq. (16) is able to repro-duce extremely well the plasmon dispersion presented,for instance, in Figs. 2 and 3. Furthermore, it is instruc-tive to note that by taking the ε m → ∞ limit in Eq. (16), (a) (b) FIG. 4. Effect of the variation of the dispersion relation ofthe plasmonic modes in a air/Ti/hBN/graphene/air structurewith the metal thickness, h . (a) Nonlocal (left, solid) andlocal (right, dashed) dispersions for three distinct values ofthe metal thickness (the orange curve is behind the greenone). (b) Plasmon wave vector, for a fixed frequency ω =2000 cm − (dot-dashed on the top panel) as a function of themetal thickness. In both panels, the remaining parametersare: s = 1 nm, E F = 0 . which corresponds to the case where the metal becomesa perfect conductor, one obtains (neglecting nonlocal ef-fects) ε x κ (4)T coth (cid:104) κ (3)T s (cid:105) + ε x κ (3)T + iσωε = 0 . (18)Equation (18) coincides with the dispersion relation ofthe acoustic plasmon branch in double-layer graphene in a symmetric dielectric environment, where the individ- ual graphene layers are separated by a distance 2 s . Thisresult reflects the scenario in which the screening exertedby the perfect conductor mirrors exactly the charges in-duced in the graphene sheet. Lastly, notice that for largeseparations, i.e., s → ∞ , the plasmon dispersion (16) re-duces to that of an isolated graphene sheet between twodielectric media ε x and ε x . B. Nonlocal Effects in the Field Distributions
Having discussed the influence of nonlocal effects inthe plasmon dispersion, we now study their impact in thespatial distribution of the fields associated with the plas-monic modes. The fields amplitudes follow from the BCs,and we determine the y –dependence of the fields at agiven ( q, ω )–point which satisfies the plasmon dispersionshown in the previous sectio. In this spirit, Fig. 5 depictsthe variation of the absolute value of the magnetic andelectric field (in logarithmic units), along the heterostruc-ture, for a (real part of the) frequency ω = 2000 cm − .It is apparent from Fig. 5(a) that magnetic field dis-tribution remains essentially unchanged under nonlocalcorrections. This is a natural consequence of the hy-drodynamic model, in which nonlocality only enters inthe longitudinal components. Since the magnetic fieldis purely transverse, nonlocal effects do not influence it,and the only differences between the local and nonlocalcases arise from small differences between the coefficients(and the value of q for the same frequency) that describethe magnetic field on each case.In contrast to the magnetic field, the inclusion of non-local effects in the metal’s response renders significantchanges in the spatial distribution of the electric field,as can be seen from Fig. 5(b)–(c). Naturally, the differ-ences between the local and nonlocal calculations arisemostly in fields within the metal region, as could be an-ticipated, due to the fact that we have allowed the ex-istence of a longitudinal mode inside the metal, in thenonlocal case. The figure demonstrates that the nonlo-cality introduced by this additional longitudinal wave hasprofound implications in the eletric field’s spatial distri-bution. Specifically, notice that the electric field insidethe metal in the local case is significantly smaller than thefield in the spacer region (as expected for a good metal).However, when nonlocal effects are taken into account,the electric field is increased (when compared with theLRA) by several orders of magnitude in the vicinity ofthe metal’s surface. This feature is consequence of thesmearing of the electron density introduced by nonlocal-ity, which translates into a larger penetration of the fieldsinside the metal. This penetration depth is practicallynegligible when taking the local approximation, but be-comes important when taking nonlocality into account,as Fig. 5 shows. In the latter, the penetration lengthreaches a couple nanometers, and can become compara-ble to the metal’s thickness for ultra-thin films.Another feature visible in Fig. 5(b)–(c) is the pres- - - - - - - - - - (a) - - - - - - - (b) FIG. 5. Spacial distribution of (a) the magnetic field and (b)the electric field components, for a air (blue)/Ti (gray)/hBN(green)/graphene (white)/air (blue) structure, calculated us-ing local and nonlocal metal models, for plasmonic modeswith the same real part of the frequency ω = 2000 cm − . H and E are respectively the nonlocal magnetic and totalelectric fields calculated at x = 0. The parameters used are: h = 10 nm, s = 1 nm, E F = 0 . ence, in the nonlocal case, of a sharp dip in the electricfield magnitude around y = − . y –component, which is caused by a destructive inter-ference between the transversal and longitudinal modesinside the metal. Note that both the transversal andlongitudinal modes have x - and y –components, with thetotal field being given by E x/y = E x/y T + E x/y L ; the nodeoccurs when E y T + E y L = 0. For that reason, across thenode there is a change of the sign of the y -componentof the electric field, dividing the regions where the lon-gitudinal mode amplitude is higher (to the right of thenode) or lower (to the left) than the transversal modeamplitude. C. Nonlocal Effects on the Plasmonic Losses
We conclude the study of the system portrayed inFig. 1 with an analysis on the effect of the nonlocalityin the losses associated with the plasmonic modes. Asmentioned above, in the presence of losses and for a real-valued q , the ensuing condition for the plasmon disper-sion requires a complex-valued ω . So far we have limitedour discussing to its real part, but here we focus on itsimaginary part, Im { ω } , since this quantity is intrinsicallyrelated with the plasmonic losses and plasmon life-time τ p , in particular, τ − = − Im { ω } /
2. In Fig. 6 we haveplotted − Im { ω } as function of the real part of the polari-ton frequency, Re { ω } , both within the local and nonlocalresponse formalism (for the metal; graphene is modeledas nonlocal in all cases).We have considered two cases, one with a relativelyhigh (for graphene) scattering rate of Γ = 16 meV, andanother illustrating high-quality, low-loss graphene, inwhich Γ = 1 meV. We note that the change of graphene’selectronic scattering rate does not significantly alter thedispersion relation presented in Fig. 2 (in fact, there is novisible difference to the eye). However, it significantly im-pacts the losses affecting graphene plasmons sustained inthe heterostructure. As Fig. 6 shows, the differences be-tween the calculation using a local and a nonlocal metal ismodest. This suggests that the losses that the plasmonsincur originate from the graphene, which is natural sincethe metal only participates indirectly—via screening—and the mode’s spectral weight is essentially attributed tographene’s.Furthermore, Fig. 6(b) shows that this differ-ence becomes more evident when the spacer width is re-duced, since the losses, when calculated locally, stronglyincrease for small spacers, whereas nonlocally their in-crease is very small. This means that, by reducing thespacer size, an increasingly higher confinement can beachieve without a strong increasing in the losses, what isa very important result.Let us now understand the sudden decrease of thelosses seen, in Fig. 6 [panel (a)] around the hBN phononfrequencies, when large damping in graphene (16 meV)is considered. To that end we recall that the damping ofthe phonon modes in hBN are 2.4 meV for the in-planephonon and 1.9 meV for the out-of-plane phonon .At the frequency of the phonons the polariton spectrumhas, essentially, a phononic nature. Therefore, the lossesare essentially controlled by those due to phonons. Sincethese are much smaller than 16 meV, we see a sudden de-crease in the losses. On the other hand, for a damping ingraphene of 1 meV, the phonon damping is comparableto the damping in graphene. As consequence the curvesof the losses in Fig. 6 [panel (a)] do not present the sud-den drop at the phonon frequencies. In conclusion, awayfrom the phonon frequencies, the losses are essentiallydue to graphene, whereas near the phonon frequencies ofhBN the losses are essentially controlled by the behav-ior of the imaginary part of the dielectric function of thehBN spacer.
500 1000 1500 2000020406080 500 1000 1500 2000 (a) (b)
FIG. 6. (a) Imaginary part of the frequency of the SPPs inan air/Ti/spacer/graphene/air configuration, with the spacerbeing hBN (left panel) or air (right panel), calculated usinglocal (dashed) and nonlocal (continuous) models, for two dif-ferent values of the graphene’s electric relaxation energy Γ,and for s = 1 nm (dot-dashed on the bottom panel). (b)Variation of the imaginary part of the SPPs, in a configura-tion with a hBN spacer, with the spacer width s , for a fixedfrequency ω = 2000 cm − (dot-dashed on the top panel). Theremaining parameters are: h = 10 nm and E F = 0 . IV. PROBING NONLOCALITY IN METALSUSING A GRAPHENE NANORIBBON GRATING
A common approach to experimentally access the plas-monic properties of a system consists on performing re-flectance and/or transmittance measurements of the far-field spectra, upon illuminating the sample.In this con-text, plasmon excitations appear as resonances in thespectra (as peaks or dips). Therefore, one may inves-tigate the influence of nonlocal effects by studying theresulting spectra and compare it to experimental data.In the extended, continuous layered system studied in the previous section, we have seen that the metal’s nonlo-cal response affects the plasmon properties of the modessupported in the system. However, the wave vector mis-match between the graphene plasmons in Fig. 1 structureand the one of a photon in free-space differs by more thantwo orders of magnitude, and thus it has been primarilyinvestigated by near-field techniques which are able toovercome this kinematic limitation .In what follows we consider a different configuration.It is similar to the one considered in the previous sec-tion, but we now assume that the graphene monolayerhas been patterned into a periodic array of graphenenanoribbons—see Fig. 7. These effectively act as adiffraction grating, whose Fourier components providethe necessary in-plane momentum to excite plasmons inthe system . Indeed, this mechanism not only allowsthe excitation of graphene plasmons by free-space pho-tons but it also serves as a platform susceptible to en-hance nonlocal effects, since the higher diffraction ordercan carry substantial momentum and thus promote non-local effects. Hence, in order to illustrate a case wherenonlocality plays a major role on the optical properties ofthe system, we will replace the graphene sheet consideredin the previous section by a graphene diffraction gratingwith a period d , and ribbon width w (see Fig. 7). We wd graphene hBNmetalSiO Siair ε ε ε ε ε IIIIIIIVV bsh kBE θ FIG. 7. Schematic representation of the layered system con-sidered in this section. It is assumed to be periodic in the x -direction and uniform in the z -direction. All the differentmaterials, dimensions and regions considered are marked inthe figure. Graphene is considered as being nonlocal through-out and the metal is assumed to be either local or nonlocal. calculate the optical properties (namely, the reflectance R , transmittance T and absorbance A ) of such struc-ture when it is illuminated by a p –polarized plane-wavecoming from region V, with monochromatic frequency ω and incident angle θ , as depicted in Fig. 7 Apart fromthe introduction of the diffraction grating, we will alsobe considering henceforth a system where the ribbon ar-ray is encapsulated between the spacer (hBN, as before)and a thick (285 nm) layer of silicon dioxide (SiO ). Thelatter’s permittivity is taken from Ref. 53. On top ofthe SiO , we further consider a semi-infinite layer of sil-icon (Si) described by a isotropic dielectric constant of ε Si = 11 . . Our theoretical calculations for the grating0system follow from the considerations outlined in the pre-vious sections and it is based on a Fourier modal expan-sion of the electromagnetic fields described elsewhere .In the interest of selfcontainedness, we provide the math-ematical details in Appendix B.Figure 8(a) shows the outcome of our computations forthe reflectance, R , transmittance, T , and absorbance, A ,spectrum for a representative graphene ribbon diffractiongrating placed near a metal, with period d = 25 nm andribbon width of d/
2. Normal incidence ( θ = 0) is as-sumed hereafter, but our formalism is valid for arbitraryincident angles (though the dependence on the imping-ing angle is weak since sin θ (cid:28) π/d ). Furthermore, asbefore, the impact of the nonlocal effects arising from thenearby metal’s response are evaluated by comparing thelocal with the nonlocal case (while maintaining graphenetreated at the nonlocal level). The spectral features vis-ible in Fig. 8(a) are rich, but these can be split into twocategories, detailed below. The first three peaks/dips—located approximately at 750, 1100, and 1350 cm − —do not present any variation due to the nonlocal ef-fects. These peaks correspond to the excitation of opticalphonons in either the hBN or the SiO slabs (indicated inFig. 8(a)). Note that both the hBN and the SiO supportphononic modes around 750 cm − . On the other hand,the position of the fourth resonance shifts considerablyto higher frequencies upon inclusion of nonlocal effectsin the metal. Moreover, in Fig. 8(b) we have plotted theextinction spectra, defined as 1 − T / T CNP (where T isthe transmittance at some finite Fermi level and T CNP is the transmittance at the graphene’s charge neutralitypoint), for several values of the Fermi energy. We fo-cus on the spectral window where the graphene plasmonlies. Clearly, the peak disperses towards higher frequen-cies with increasing the Fermi level. This unarguablydemonstrates that this resonance corresponds indeed tothe excitation of a graphene plasmons in the array. Theadvantage of exciting these plasmons using a grating isthat their frequency is highly dependent on the valuechosen for the period and ribbon width of the graphenearray, with smaller values of the period yielding concomi-tant larger wave vectors and thus plasmon resonances athigher frequencies. For a Fermi level of 0 . − , which in turncorresponds to a relative variation of more than 10%.This plainly shows that, in the case of a graphene sheetplaced very close a metallic film, the proper account ofnonlocal effects both in the metal and in graphene arekey to properly model and interpret experimental resultsin such a system.A solution for this problem may be to consider thespacer width s to be an effective fitting parameter, s eff ,higher than the actual experimental value, as shown inFig. 8(c).In Fig. 8(c) is depicted the nonlocal curve for s = 1 nmoverlaid by local curves for different values of the effectivespacer thickness; for s eff = 1 .
75 nm, the correspondence (a) (b) (c)
FIG. 8. (a) Reflectance ( R ), transmittance ( T ), and ab-sorbance ( A ) spectra of an air/Ti/hBN/graphene/SiO /Sisystem. Graphene ( E F = 0 . α = s eff /s , compared to the nonlo-cal counterpart. In all cases, the remaining parameters are: h = 10 nm, s = 1 nm, b = 285 nm, d = 25 nm, w = d/ θ = 0 ◦ . V. CONCLUSIONS
Nonlocality in metals is strongly influenced by sev-eral material-dependent parameters, namely the plasmafrequency, the Fermi velocity, the momentum relaxationrate, and the background dielectric function of the coreelectrons in the metal. In bulk metals, nonlocal effectscan be safely neglected since in that case the wavenumberof the radiation interacting with the metal is, in general,small ( q/k F (cid:28) q/k F ∼ ACKNOWLEDGEMENTS
The authors thank S´ebastien Nanot and Itai Ep-stein for valuable discussions and comments. E.J.C.D.,Yu.V.B. and N.M.R.P. acknowledge support from theEuropean Commission through the project Graphene-Driven Revolutions in ICT and Beyond (Ref. No.785219), and from the Portuguese Foundation for Scienceand Technology (FCT) in the framework of the Strate-gic Financing UID/FIS/04650/2013. E.J.C.D. acknowl-edges FCT for the grant CFUM-BI-14/2016. D.A.I. ac-knowledges the FPI grant BES-2014-068504. F.H.L.K.acknowledges financial support from the Government ofCatalonia trough the SGR grant (2014-SGR-1535), andfrom the Spanish Ministry of Economy and Competi-tiveness, through the Severo Ochoa Programme for Cen-tres of Excellence in R&D (SEV-2015-0522), support byFundacio Cellex Barcelona, CERCA Programme / Gen-eralitat de Catalunya and the Mineco grants Ramn yCajal (RYC-2012-12281) and Plan Nacional (FIS2013-47161-P and FIS2014-59639-JIN). Furthermore, the re-search leading to these results has received funding fromthe European Union Seventh Framework Programme un-der grant agreement no.696656 Graphene Flagship, theERC starting grant (307806, CarbonLight), and projectGRASP (FP7-ICT-2013-613024-GRASP). N. A. M. is aVILLUM Investigator supported by VILLUM FONDEN(grant No. 16498). Center for Nano Optics is financiallysupported by the University of Southern Denmark (SDU2020 funding). Center for Nanostructured Graphene issupported by the Danish National Research Foundation(DNRF103).
Appendix A: Calculation of the PlasmonicProperties
In this appendix, we describe the procedure we em-ployed to calculate the dispersion relation and all subse-quent analysis of the system presented in Fig. 1. We willassume, with full generality, that all the dielectric mediais described by two dielectric functions along the in-plane( ε x ) and out-of-plane ( ε y ) directions.Following the guidelines discussed in Section II, wewrite the magnetic field in either of the 4 regions thatcompose the system as H z I ( x, y ) = η e κ (1)T y e i qx , (A1) H z II ( x, y ) = (cid:16) γ + e κ (2)T y + γ − e − κ (2)T y (cid:17) e i qx , (A2) H z III ( x, y ) = (cid:16) α + e κ (2)T y + α − e − κ (2)T y (cid:17) e i qx , (A3) H z IV ( x, y ) = ζ e − κ (4)T y e i qx , (A4)where η , γ ± , α ± and ζ are undetermined coefficients, q is the in-plane momentum of the modes and κ ( ν )T = (cid:112) q ε xν /ε yν − ε xν ω /c is the out-of-plane momentum ofthe transversal modes.2The electric fields, on the other hand, are given by theexpressions E x I ( x, y ) = i κ (1)T ωε ε x η e κ (1)T y e i qx , (A5) E x III ( x, y ) = i κ (3)T ωε ε x (cid:16) α + e κ (3)T y − α − e − κ (3)T y (cid:17) e i qx , (A6) E x IV ( x, y ) = − i κ (4)T ωε ε x ζ e − κ (4)T y e i qx , (A7)for the dielectric regions, and by E x II ( x, y ) = (cid:34) i κ (2)T ωε ε x (cid:16) γ + e κ (2)T y − γ − e − κ (2)T y (cid:17) ++ (cid:16) δ + e κ (2)L y + δ − e − κ (2)L y (cid:17)(cid:105) e i qx , (A8)for the metallic region, with δ ± being two additional un-known coefficients that arise from the nonlocality, and κ (2)L = (cid:113) q − (cid:0) ω + i γω − ω /ε ∞ (cid:1) /β is the verticalmomentum of the longitudinal modes. Furthermore, inregion II, the normal component of the electric field isgiven by E y II ( x, y ) = (cid:20) qωε ε x (cid:16) γ + e κ (2)T y + γ − e − κ (2)T y (cid:17) ++ κ (2)L i q (cid:16) δ + e κ (2)L y − δ − e − κ (2)L y (cid:17)(cid:35) e i qx , (A9)which fully determines the polarization in the y -direction, though Eq. (15), with the form P y II ( x, y ) =( q/ω ) H z II ( x, y ) − ε ε ∞ E y II ( x, y ).There are then 8 undetermined coefficients in the prob-lem: η , γ ± , δ ± , α ± and ζ , which are determined throughthe application of the boundary conditions of the prob-lem: H z I ( x, − h ) = H z II ( x, − h ) , (A10) E x I ( x, − h ) = E x II ( x, − h ) , (A11) P y II ( x, − h ) = 0 , (A12) H z II ( x,
0) = H z III ( x, , (A13) E x II ( x,
0) = E x III ( x, , (A14) P y II ( x,
0) = 0 , (A15) H z III ( x, s ) − H z IV ( x, s ) = − σE x IV ( x, s ) , (A16) E x III ( x, s ) = E x IV ( x, s ) , (A17)where σ ( q, ω ) is the conductivity of the graphene sheet,calculated using Mermin’s formula in Appendix C.It is straightforward to check that the 8 equations(A10)–(A17) form an undetermined system of equations,meaning that it cannot be solved for all the unknown co-efficients. In order to calculate the dispersion relation ofthe allowed plasmonic modes, we need therefore to define ζ ≡ H as a free parameter of the problem, and find the remaining coefficients as a function of this parameter.Using arbitrarily the boundary equations (A10)–(A16),we find easily the coefficients η , γ ± , δ ± and α ± normal-ized to H ; we do not show the actual expressions becausethey are rather ugly and very unintuitive.The remaining boundary condition (A17), on the otherhand, can be solved in order to find the relation between q and ω which ensures the solubility of the system ofequations; it corresponds to the dispersion relation of theallowed plasmonic solution. This equation is given by κ (3)T ε x (cid:18) α + H e κ (3)T s − α − H e − κ (3)T s (cid:19) = − κ (4)T ε x e − κ (4)T s , (A18)where we need to bear in mind that α ± /H are, at thispoint, totally determined coefficients (function of both q and ω , in general).The process described above takes explicitly in con-sideration the nonlocal effects. The counterpart localdispersion relation, on its turn, can be calculated whenexplicitly considering δ ± = 0 and ignoring the bound-ary conditions (A12) and (A15), proceeding analogouslyotherwise. Appendix B: Calculation of the Optical Properties
In this appendix, we describe the procedure we em-ployed to calculate the optical properties of the systempresented in Fig. 7. We will be considering that the me-dia in regions III and IV may be axial, and we will thusdescribe it by two dielectric functions along the in-plane( ε x ) and out-of-plane ( ε y ) directions. We limit this gen-eralizations to these regions, because if media I and/or Vwere axial, it would have consequences in the definitionsof the reflectance and the transmittance that are out ofthe scope of this work.Let us consider that the impinging light carries a mo-mentum k = k x ˆx − k y ˆy with k x = k sin( θ ), k y = k cos( θ ) and k = √ ε ω/c . Due to the grating, thisfield will be diffracted, so the scattered fields need to bedescribed by a combination of several different diffractionmodes n , written in each region as H z I ( x, y ) = H (cid:88) n τ n e − i k (1)T ,n y e i ρ n x , (B1) H z II ( x, y ) = H (cid:88) n (cid:16) γ + n e i k (2)T ,n y + γ − n e − i k (2)T ,n y (cid:17) e i ρ n x , (B2) H z III ( x, y ) = H (cid:88) n (cid:16) α + n e i k (3)T ,n y + α − n e − i k (3)T ,n y (cid:17) e i ρ n x , (B3) H z IV ( x, y ) = H (cid:88) n (cid:16) φ + n e i k (4)T ,n y + φ − n e − i k (4)T ,n y (cid:17) e i ρ n x , (B4) H z V ( x, y ) = H e i k x x e − i k y y + H (cid:88) n r n e i k (5)T ,n y e i ρ n x , (B5)3where ρ n = k x + 2 nπ/d is the in-plane momentum of the n th mode, as settled by the Bloch Theorem, and k ( ν )T ,n = (cid:112) ε xν ω /c − ρ n ε xν /ε yν is the out-of-plane momentum inthe region ν for the transversal waves.The electric field, in its turn, is written as E x I ( x, y ) = H (cid:18) ωε ε x (cid:19) (cid:88) n τ n k (1)T ,n e − i k (1)T ,n y e i ρ n x , (B6)and analougously for regions III, IV and V, whereas, forregion II, it is written as, E x II ( x, y ) = (cid:34) − k (2)T ,n ωε ε x (cid:16) γ + e i k (2)T ,n y − γ − e − i k (2)T ,n y (cid:17) ++ (cid:16) δ + e i k (2)L ,n y + δ − e − i k (2)L ,n y (cid:17)(cid:105) e i ρ n x , (B7)where k L ,n = (cid:113)(cid:0) ω + i γω − ω /ε ∞ (cid:1) /β − ρ n is the out-of-plane momentum for the longitudinal waves. The y-component E y II is defined analogously. Note that k L ,n and k T ,n were defined appropriately to describe propagatingfields.There is a total of 10 undetermined coefficients in theproblem for each mode n ( τ n , γ ± n , δ ± n , α ± n , φ ± n and r n ),which must be determined through the boundary condi-tions of the system. Boundary conditions (A10)–(A15)and (A17) still hold, to which we need to add two moreconditions, H z IV ( x, s + b ) = H z V ( x, s + b ) , (B8) E x IV ( x, s + b ) = E x V ( x, s + b ) . (B9)To illustrate the next step, let us consider boundarycondition (A10), with the form H (cid:88) n τ n e i k (1)T ,n h e i ρ n x = H (cid:88) n (cid:16) γ + n e − i k (2)T ,n h + γ − n e i k (2)T ,n h (cid:17) e i ρ n x . (B10)Multiplying the previous equation, on both sides, bye − i ρ (cid:96) x and integrating the resulting equation in one unitcell of the system ( | x | < d/ (cid:82) d/ − d/ dx e i( ρ n − ρ (cid:96) ) x = dδ n(cid:96) , which then yield τ (cid:96) e i k (1)T ,(cid:96) h = γ + (cid:96) e − i k (2)T ,(cid:96) h + γ − (cid:96) e i k (2)T ,(cid:96) h . (B11)Notice that the previous equation, for each mode (cid:96) ,relates only the coefficients of the same index (cid:96) —or, inother words, the boundary condition is obeyed by each in-dividual mode that composes the total field, and not onlyby the total field itself, what strongly simplifies the prob-lem. This property is equally obeyed by the remaining(A11)–(A15), (A17) and (B8)–(B9) conditions, meaningthat, proceeding as illustrated above, we arrive at 9 equa-tions which relate all the coefficients for some particular mode (cid:96) . These form a determined system of equationswhich allows for the calculation of 9 out of the 10 differ-ent coefficients, in function of the remaining one. We willhenceforth take the undetermined coefficient to be the r (cid:96) ;in that case, one finds that all the other coefficients arerelated to r (cid:96) by a linear relation of the form α ± (cid:96) = θ ± δ (cid:96) + r (cid:96) ψ ± (cid:96) , φ ± (cid:96) = λ ± δ (cid:96) + r (cid:96) χ ± (cid:96) , (B12)and analogously for τ (cid:96) , γ ± (cid:96) and δ ± (cid:96) .The remaining unknown coefficient, r (cid:96) , must nowbe determined using the remaining boundary condition,which corresponds to the discontinuity of the magneticfield across the graphene grating, due to the presence ofsurface currents in the graphene. This condition, anal-ogous to the condition expressed in Eq. (A16), has theexplicit form (cid:88) n (cid:16) α + n e i k (3)T ,n s + α − n e − i k (3)T ,n s (cid:17) e i ρ n x −− (cid:88) n (cid:16) φ + n e i k (4)T ,n s + φ − n e − i k (4)T ,n s (cid:17) e i ρ n x = (cid:88) n (cid:32) σ n k (4)T ,n ωε ε x (cid:33) (cid:16) φ + n e i k (4)T ,n s − φ − n e − i k (4)T ,n s (cid:17) e i ρ n x , (B13)where σ n = σ ( ρ n , ω ) is the conductivity of the graphene.The big difference between this condition and all the oth-ers is that graphene’s conductivity is only non-zero for | x | < w/
2; this means that, when multiplying both sidesof the equation by e − i ρ (cid:96) x and integrating it on the unitcell of the system, we obtain˜ α + (cid:96) + ˜ α − (cid:96) − ˜ φ + (cid:96) − ˜ φ − (cid:96) = (cid:88) n (cid:32) σ n k (4)T ,n ωε ε x (cid:33) (cid:16) ˜ φ + n − ˜ φ − n (cid:17) S (cid:96)n , (B14)with ˜ α ± n ≡ α ± n e ± i k (3)T ,n s (and equivalently for ˜ φ ± n ), and S (cid:96)n = (cid:90) w/ − w/ dx e i( ρ n − ρ (cid:96) ) x = sin [ πw ( n − (cid:96) ) /d ] π ( n − (cid:96) ) . (B15)This equation, unlike all the others, relate coefficientsof all the allowed diffraction orders n . To solve it, weemploy equations (B12) and, upon some mathematicalmanipulation, we arrive at an expression of the form (cid:80) n M (cid:96)n r n = − F (cid:96) , where F (cid:96) = (˜ θ +0 + ˜ θ − − ˜ λ +0 − ˜ λ − ) δ (cid:96) − (cid:34) σ k (4)T , ωε ε x (cid:35) (˜ λ +0 − ˜ λ − ) S (cid:96) , (B16) M (cid:96)n = ( ˜ ψ + n + ˜ ψ − n − ˜ χ + n − ˜ χ − n ) δ (cid:96)n − (cid:34) σ n k (4)T ,n ωε ε x (cid:35) ( ˜ χ + n − ˜ χ − n ) S (cid:96)n (B17)This equation may be written in the matrix form as M · R = − F , with M being a square matrix with elements4[ M ] (cid:96)n = M (cid:96)n , and R and F being columns with elements[ R ] (cid:96) = r (cid:96) and [ F ] (cid:96) = F (cid:96) . The r (cid:96) are hence determinedthrough the solution of that matrix equation, written fora high-enough matrix dimension to ensure the conver-gence of the solution. The remaining coefficients are thendetermined by equations (B12) and the analogous for τ n , γ ± n and δ ± n .After this process, we can find the reflectance andtransmittance of the system by the equations R = (cid:88) n ∈ PM Re (cid:34) k (5)T ,n ε (cid:35) Re (cid:20) ε k y (cid:21) | r n | , (B18) T = (cid:88) n ∈ PM Re (cid:34) k (1)T ,n ε x (cid:35) Re (cid:20) ε k y (cid:21) | τ n | , (B19)where the summations are performed over the propagat-ing modes (PM). The absorbance, on the other hand, isdefined as A = 1 − R − T .Once again, the optical properties for the local caseare calculating when repeating the same procedure whilstsetting δ ± n = 0 and ignoring boundary conditions (A12)and (A15).Apart from the optical properties, it is interesting tonote that, from the reflectance amplitudes r n , one cancalculate the loss function L ( ω, q ) = − (cid:80) n Im[ r n ] (where k x is substituted by an arbitrarily changeable momentum q ), which is an alternative way to calculate the dispersionrelation of the allowed bound modes of the problem. Thisis shown in Fig. 2, where the loss function was overlaidby the dispersion relation calculated as detailed in Ap-pendix A, and the correspondence between the two plotsis excellent. Appendix C: Graphene’s Nonlocal Conductivity
For the conductivity of the graphene sheet, we haveused Mermin’s formula , which includes nonlocal ef-fects. Let x ≡ q/k F and y ≡ (cid:126) ω/E F be dimensional vari-ables constructed from q and ω , respectively. E F refersto the graphene’s Fermi energy, k F = E F / ( (cid:126) v F ) is theFermi momentum ( v F ≈ c/
300 is the Fermi speed) and Γis the material’s relaxation energy. The formula we haveused was retrieved from Gon¸calves and Peres , σ ( q, ω ) = 4i σ (cid:126) ωq χ τ (cid:18) qk F , (cid:126) ωE F (cid:19) , (C1)with σ ≡ e / (4 (cid:126) ) and χ τ ( x, y ) = (cid:16) Γ yE F (cid:17) χ g (cid:16) x, y + i Γ E F (cid:17) Γ yE F χ g (cid:16) x, y + i Γ E F (cid:17) /χ g ( x, . (C2)The function χ g ( x, y ) is calculated differently accord-ing to the region where it is calculated in the xy -phase FIG. 9. Regions for the calculation of the graphene suscepti-bility in the xy -phase space. space represented in Fig. 9. It can be written as χ g ( x, y ) = χ (1)B ( x, y ) , Re[ y ] > x ∧ Re[ y ] < − x,χ (2)B ( x, y ) , Re[ y ] > x ∧ Re[ y ] > − x,χ (3)B ( x, y ) , Re[ y ] > x + 2 ,χ (1)A ( x, y ) , Re[ y ] < x ∧ Re[ y ] < − x,χ (2)A ( x, y ) , Re[ y ] < x ∧ Re[ y ] > − x,χ (3)A ( x, y ) , Re[ y ] < x − , (C3)where Re[ y ] stands for the real part of y , and each of thefunctions in the previous expression are given by χ (1)B ( x, y ) = − π E F ( (cid:126) v F ) + 14 π E F ( (cid:126) v F ) x (cid:112) y − x ·· (cid:20) F (cid:18) y + 2 x (cid:19) − F (cid:18) − yx (cid:19)(cid:21) , (C4) χ (2)B ( x, y ) = − π E F ( (cid:126) v F ) + 14 π E F ( (cid:126) v F ) x (cid:112) y − x ·· (cid:20) F (cid:18) y + 2 x (cid:19) + i G (cid:18) − yx (cid:19)(cid:21) , (C5) χ (3)B ( x, y ) = − π E F ( (cid:126) v F ) + 14 π E F ( (cid:126) v F ) x (cid:112) y − x ·· (cid:20) − i π + F (cid:18) y + 2 x (cid:19) − F (cid:18) y − x (cid:19)(cid:21) , (C6)5 χ (1)A ( x, y ) = − π E F ( (cid:126) v F ) − i4 π E F ( (cid:126) v F ) x (cid:112) x − y ·· (cid:20) F (cid:18) y + 2 x (cid:19) − F (cid:18) − yx (cid:19)(cid:21) , (C7) χ (2)A ( x, y ) = − π E F ( (cid:126) v F ) + i4 π E F ( (cid:126) v F ) x (cid:112) x − y ·· (cid:20) i π − F (cid:18) y + 2 x (cid:19) + i G (cid:18) − yx (cid:19)(cid:21) , (C8) χ (3)A ( x, y ) = − π E F ( (cid:126) v F ) + 14 π E F ( (cid:126) v F ) x (cid:112) x − y ·· (cid:20) − π + G (cid:18) y + 2 x (cid:19) − G (cid:18) y − x (cid:19)(cid:21) . (C9)with the functions F ( x ) ≡ x √ x − − arccosh( x ) and G ( x ) ≡ x √ − x − arccosh( x ). ∗ eduardo.dias@fisica.uminho.pt † peres@fisica.uminho.pt S. A. Maier,
Plasmonics: Fundamentals and Applications (Springer, 2007). M. Pelton and G. W. Bryant,
Introduction to Metal-Nanoparticle Plasmonics (Wiley, 2013). P. A. D. Gon¸calves and N. M. R. Peres,
An Introduction toGraphene Plasmonics , 1st ed. (World Scientific, Singapore,2016). V. Klimov,
Nanoplasmonics (CRC Press, 2014). A. I. Fern´andez-Dom´ınguez, F. J. Garc´ıa-Vidal, andL. Mart´ın-Moreno, Nature Photonics , 8 (2017). F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy,Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi,B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg,Science , 726 (2016). C. Cirac`ı, R. Hill, J. J. Mock, Y. Urzhumov, A. I.Fern´andez-Dom´ınguez, S. A. Maier, J. B. Pendry,A. Chilkoti, and D. R. Smith, Science , 1072 (2012). S. Raza, S. Kadkhodazadeh, T. Christensen, M. Di Vece,M. Wubs, N. A. Mortensen, and N. Stenger, Nature Com-munications , 8788 (2015). M. B. Lundeberg, Y. Gao, R. Asgari, C. Tan, B. Van Dup-pen, M. Autore, P. Alonso-Gonz´alez, A. Woessner,K. Watanabe, T. Taniguchi, R. Hillenbrand, J. Hone,M. Polini, and F. H. L. Koppens, Science , 187 (2017). S. Raza, G. Toscano, A.-P. Jauho, M. Wubs, and N. A.Mortensen, Physical Review B , 121412(R) (2011). N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, andS. I. Bozhevolnyi, Nature Communications , 3809 (2014). J. D. Jackson,
Classical Electrodynamics , 3rd ed. (Wiley,New York, 1998). S. Raza, T. Christensen, M. Wubs, S. I. Bozhevolnyi, andN. A. Mortensen, Physical Review B , 115401 (2013). C. David and F. J. Garc´ıa de Abajo, The Journal of Phys-ical Chemistry C , 19470 (2011). A. I. Fern´andez-Dom´ınguez, A. Wiener, F. J. Garc´ıa-Vidal,S. A. Maier, and J. B. Pendry, Physical Review Letters , 106802 (2012). J. A. Scholl, A. L. Koh, and J. A. Dionne, Nature ,421 (2012). J. Khurgin, W.-Y. Tsai, D. P. Tsai, and G. Sun, ACSPhotonics , 2871 (2017). F. J. Garc´ıa de Abajo, The Journal of Physical ChemistryC , 17983 (2008). C. David and F. J. Garc´ıa de Abajo, AIP Conference Pro-ceedings , 43 (2010). A. Wiener, H. Duan, M. Bosman, A. P. Horsfield, J. B.Pendry, J. K. W. Yang, S. A. Maier, and A. I. Fern´andez-Dom´ınguez, ACS Nano , 6287 (2013). O. Schnitzer, V. Giannini, R. V. Craster, and S. A. Maier,Physical Review B , 041409 (2016). G. Toscano, S. Raza, A.-P. Jauho, N. A. Mortensen, andM. Wubs, Optics Express , 4176 (2012). H. Shen, L. Chen, L. Ferrari, M.-H. Lin, N. A. Mortensen,S. Gwo, and Z. Liu, Nano Letters , 2234 (2017). L. Stella, P. Zhang, F. J. Garc´ıa-Vidal, A. Rubio, andP. Garc´ıa-Gonzlez, The Journal of Physical Chemistry C , 8941 (2013). A. Varas, P. Garc´ıa-Gonz´alez, J. Feist, F. J. Garc´ıa-Vidal,and A. Rubio, Nanophotonics , 409 (2016). A. D. Boardman, B. V. Paranjape, and R. Teshima, Sur-face Science , 275 (1975). A. D. Boardman,
Electromagnetic surface modes (JohnWiley & Sons, 1982). F. Bloch, Annalen der Physik , 285 (1933). F. Bloch, Zeitschrift f¨ur Physik A Hadrons and Nuclei ,363 (1933). H. Jensen, Zeitschrift f¨ur Physik , 620 (1937). A. Moreau, C. Cirac`ı, and D. R. Smith, Physical ReviewB , 045401 (2013). C. Cirac`ı, J. B. Pendry, and D. R. Smith, ChemPhysChem , 1109 (2013). S. Raza, S. I. Bozhevolnyi, M. Wubs, and N. A. Mortensen,Journal of Physics: Condensed Matter , 183204 (2015). G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl,F. Evers, H. Xu, N. A. Mortensen, and M. Wubs, NatureCommunications , 7132 (2015). K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V.Khotkevich, S. V. Morozov, and A. K. Geim, Proceedingsof the National Academy of Sciences of the United Statesof America , 10451 (2005). A. Woessner, M. Lundeberg, Y. Gao, A. Prin-cipi, P. Alonso-Gonz´alez, M. Carrega, K. Watanabe,T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillen-brand, and F. H. L. Koppens, Nature Materials , 421(2015). A. Principi, M. Carrega, M. B. Lundeberg, A. Woessner,F. H. L. Koppens, G. Vignale, and M. Polini, PhysicalReview B , 165408 (2014). P. A. D. Gon¸calves, E. J. C. Dias, Y. V. Bludov, andN. M. R. Peres, Physical Review B , 195421 (2016). E. J. C. Dias and N. M. R. Peres, ACS Photonics , 3071(2017). We do not consider in this work the case of a dielec-tric/metal interface separated by a graphene sheet. Themain reason is that, under those conditions, doping thegraphene in order to adjust its Fermi level would be im-practical, since the metal would drain the free electrons inthe graphene sheet. V. W. Brar, M. S. Jang, M. Sherrott, S. Kim, J. J. Lopez,L. B. Kim, M. Choi, and H. Atwater, Nano letters ,3876 (2014). A. D. Raki´c, A. B. Djuriˇsi´c, J. M. Elazar, and M. L.Majewski, Applied Optics , 5271 (1998). T. Brandt, M. H¨ovel, B. Gompf, and M. Dressel, PhysicalReview B , 205409 (2008). D. Barchiesi and T. Grosges, Journal of Nanophotonics ,083097 (2014). A. Derkachova, K. Kolwas, and I. Demchenko, Plasmonics , 941 (2016). J. S. Lehtinen, T. Sajavaara, K. Y. Arutyunov, M. Y. Pres-njakov, and A. L. Vasiliev, Physical Review B , 094508(2012). N. W. Ashcroft and N. D. Mermin,
Solid State Physics ,HRW International Editions (Holt, Rinehart and Winston,1976). Y. Luo, A. I. Fern´andez-Dom´ınguez, A. Wiener, S. A. Maier, and J. B. Pendry, Physical Review Letters ,093901 (2013). Y. Luo, R. Zhao, and J. B. Pendry, Proceedings of theNational Academy of Sciences , 18422 (2014). W. Yan, M. Wubs, and N. A. Mortensen, Physical ReviewLetters , 137403 (2015). T. Christensen, W. Yan, A.-P. Jauho, M. Soljaˇci´c, andN. A. Mortensen, Physical Review Letters , 157402(2017). P. Alonso-Gonz´alez, A. Y. Nikitin, Y. Gao, A. Woess-ner, M. B. Lundeberg, A. Principi, N. Forcellini, W. Yan,S. V´elez, A. J. Huber, K. Watanabe, T. Taniguchi,F. Casanova, L. E. Hueso, M. Polini, J. Hone, F. H. L.Koppens, and R. Hillenbrand, Nature nanotechnology ,31 (2017). E. D. Palik,
Handbook of optical constants of solids , Vol. 3(Academic Press, 1998). A. Dargys and J. Kundrotas,
Handbook on physical prop-erties of Ge, Si, GaAs and InP (Science and EncyclopediaPubl., 1994). N. D. Mermin, Physical Review B , 2362 (1970). A. K. Das, Journal of Physics F: Metal Physics5