A A planar graphene sheet

# Localized surface plasmons in a continuous and flat graphene sheet

## Abstract

We derive an integral equation describing surface-plasmon polaritons in graphene deposited on a substrate with a planar surface and a dielectric protrusion in the opposite surface of the dielectric slab. We show that the problem is mathematically equivalent to the solution of a Fredholm equation, which we solve exactly. In addition, we show that the dispersion relation of the localized surface plasmons is determined by the geometric parameters of the protrusion alone. We also show that such system supports both even and odd modes. We give the electrostatic potential and the stream plot of the electrostatic field, which clearly show the localized nature of the surface plasmons in a continuous and flat graphene sheet.

## I Introduction

Light-matter interaction at the nanoscale is the realm of plasmonics. In the visible and near-IR spectral range one generally relies on the plasmonic properties of noble metals, such as Gold, Silver, and Copper. On the other hand, in the mid-infrared (IR) and in the THz the use of noble metals is excluded due to poor confinement of the surface plasmons. It is in this context that graphene emerges as a platform for plasmonics in the mid-IR and in the THz spectral range, since this material supports strongly confined surface plasmons in this frequency region Ju et al. (2011); Yan et al. (2012).

In the field of plasmonics one can distinguish between two different types of surface plasmons. For graphene, as for noble metals, the types are: (i) surface-plasmon polaritons and (ii) localized surface plasmons. The former are propagating surface waves on the graphene surface, whereas the latter are localized excitations in graphene nanostructures, such as ribbons Ju et al. (2011) and disksYan et al. (2012). Thus, in general, for obtaining localized plasmons in graphene one has to pattern the graphene sheet, which hinders the quality factor of these excitations. It would be, therefore, convenient, to provide a method to induce localized plasmons in a continuous graphene sheet. To investigate this possibility is the purpose of this paper.

The coupled mode of an electromagnetic field with charge density oscillations of a conductor is called surface plasmon-polaritons (SPP). The synthesis of graphene and other new two dimensional (2D) materials opened the door to natural candidates to support this type of surface modes, therefore creating a new field inside plasmonics Gonçalves and Peres (2016); Basov et al. (2016); Low et al. (2017). The most interesting properties of the SPP are the confinement of light below the diffraction limit Atwater (2007). The list of applications for SPP is extensive, including biochemical sensing Brolo (2012); Acimovic et al. (2014); Rodrigo et al. (2015), solar cells Green and Pillai (2012), optical tweezers Reece (2008); Juan et al. (2011), and transformation optics Huidobro et al. (2010); Pendry et al. (2012). In planar dielectric-metal interfaces, the SPP is confined along the direction transverse to the surface. However, at dielectric gaps, such as V-shaped groove and nanogaps, such plasmons can be also be confined along the non translational-invariant direction, being classified as channel plasmon-polaritons (CPP) Smith et al. (2015). Those systems can be used to steer SPP, thus forming plasmonic waveguides Nielsen et al. (2008); Gramotnev and Pile (2004).

Inside the sub-field of 2D photonics, graphene attracted much attention due to its unique properties, including the fact that its optical properties can be controlled externaly by electrostatic gating, originating long-lived plasmons, with large field confinement in the THz and mid-IR spectral range Gonçalves and Peres (2016); Xiao et al. (2016); García de Abajo (2014). On the other hand, there are few studies on CPP in graphene. The existing ones discuss 2D nano-slitsGonçalves et al. (2017a) and covered grooves and wedges Liu et al. (2013); Smirnova et al. (2016); Gonçalves et al. (2017b). In all these approaches graphene is deposited in a deformed substrate, thus assuming the same shape of the substrate. This approaches reduce the quality factor of the SPP. Therefore it would be ideal to find a method where a continuous graphene sheet is deposited on a flat substrate but would still support localized (or channel) plasmons. Here we propose a new approach: we take a susbtrate that is flat in one of its surfaces and patterned in the other surface, as can be seen in Fig. 1. The graphene sheet is then deposited on the flat region of the substrate, therefore keeping its natural flatness.

Studies of plasmon resonances in a surface with a protuberance or depression were first performed in Ref. Maradudin and Visscher, 1985, for the case of noble metals. In Ref. Sturman et al., 2014 localized plasmons in nanoscale pertubations were studied using the integral equation eigenvalue method, within a quasi-static approximation Mayergoyz et al. (2005). The scattering of SPP by a localized defect in a dielectric with axial symmetry was performed in Ref. Arias and Maradudin, 2013 using the reduced Rayleigh equations Zayats et al. (2005). There are numerous works about plasmonic resonances in rough or periodic surfaces, using similar procedures as the one employed for the study of a single protuberance Raether (1988); Pereira et al. (2003); Chubchev et al. (2017). All these works refer to noble-metal plasmonics and similar geometries have not been yet considered for graphene.

In this paper we calculate, using an electrostatic approximation, the plasmonic localized modes of a graphene sheet deposited on a flat substrate which has either a protuberance or an indentation in the opposite face. In Section II we discuss the procedure to solve a generic 2D dielectric protrusion in the electrostatic approximation and in Section III we derive an integral equation for the Fourier coefficient of the potential field for a generic 1D deformation. In Section IV we discuss the classification of the integral equation, the condition for the existence of localized plasmons, and the numerical procedure. In Section V we discuss our results for a Gaussian profile.

## Ii A planar graphene sheet on a dielectric defect: 2D protrusion

Let us consider the geometry of Fig. 1. There are three regions in this system: the region , with dielectric function , the region between , with dielectric function , and the region , with dielectric function . We have to define the electrostatic potential in these three regions.

In all regions we write the electrostatic potential as a Fourier integral. In the first region we write

 ϕ1(ρ,z,t)=∫dk∥(2π)2A(k∥)eik∥⋅ρ−k∥ze−iωt, (1)

where and . In the central region both real exponentials have to be present, that is,

 ϕc(ρ,z,t)=∫dk∥(2π)2[B(k∥)ek∥z+C(k∥)e−k∥z]eik∥⋅ρ−iωt. (2)

Finally, in the third region we have

 ϕ2(ρ,z,t)=∫dk∥(2π)2D(k∥)eik∥⋅ρ+k∥ze−iωt. (3)

Next we assume that the above expressions for and hold in the selvage region. The boundary conditions are imposed at and at , where is some even function, for example an inverted Gaussian:

 ζ(x,y)=−ζ0e−ρ2/s2, (4)

where . At the boundary conditions are the same we have used for solving the flat graphene case discussed in Appendix A, whereas for , the boundary conditions are adapted from those at considering that (the optical conductivity), and therefore the normal component of the electric displacement field is continuous through the interface. Although we have formulated the problem for a defect with cylindrical symmetry, we can also consider one-dimensional profiles, such as , which will be the case considered next.

## Iii A planar graphene sheet on a dielectric defect: 1D protrusion

From here on we consider a 1D defect. In this case the Fourier representation of the field is one-dimensional, reading

 ϕ1(ρ,z,t)=∫dkx2πA(kx)ei(kxx+kyy)−k∥ze−iωt, (5)

and similar equations for and . Next we want to obtain an eigenvalue equation, thus allowing us to determine the eigen-frequencies, in terms of a single coefficient. In particular, we want that coefficient to be . For implementing the boundary conditions we need an expression for the normal derivative along the surface of the defect. This is given by

 ∂∂n=^n⋅∇=[1+(∂xζ(x))2]−1/2(−∂xζ(x)∂x+∂z). (6)

At the interface we have simply

 ∂∂n=∂∂z. (7)

Thus the boundary conditions are

 ϕ1(ρ,0,t) =ϕc(ρ,0,t), (8a) ϵ1∂ϕ1(ρ,0,t)∂z −ϵ(ω)∂ϕc(ρ,0,t)∂z=−iσϵ0ω∇22Dϕ(ρ,0,t), (8b)

which turn into

 A(kx) =B(kx)+C(kx), (9a) −ϵ1A(kx) −ϵ(ω)[B(kx)−C(kx)]=κA(kx), (9b)

where we defined:

 κ≡iσk∥ϵ0ω, (10)

and the solution for and reads

 B(kx) =A(kx)ϵ(ω)−ϵ1−κ2ϵ(ω)≡A(kx)f+(ω,k∥), (11a) C(kx) =A(kx)ϵ(ω)+ϵ1+κ2ϵ(ω)≡A(kx)f−(ω,k∥). (11b)

This allows us to write the field in the central region as a function of the alone. The next step is the implementation of the boundary conditions at the interface . These are

 ϕc(ρ,z2c(x),t) =ϕ2(ρ,z2c(x),t), (12a) ϵ(ω)∂ϕc(ρ,z2c(x),t)∂n =ϵ2∂ϕ2(ρ,z2c(x),t)∂n. (12b)

The boundary condition (12a) is simply given by

 Missing or unrecognized delimiter for \left +f−(ω,k∥)e−k∥z2c(x)]=∫dkx2πD(kx)ei(kxx+kyy)ek∥z2c(x) (13)

The second boundary condition, Eq. (12b), reads

 ϵ2∫dkx2πD(kx)[−∂ζ(x)∂xikx+k∥]ei(kxx+kyy)ek∥z2c(x)= ϵ(ω)∫dkx2πA(kx)ei(kxx+kyy)[f+(ω,k∥)ek∥z2c(x)× ×(−∂ζ(x)∂xikx+k∥)+f−(ω,k∥)e−k∥z2c(x)× ×(−∂ζ(x)∂xikx−k∥)] (14)

Now we need to combine Eqs. (13) and (14) for obtaining a single integral equation for the coefficient . This is a more difficult task since we have the function in the exponent together with derivatives of . For circumventing this difficulty we introduce the Fourier representation of the exponential as

 eαζ(x)=1+α∫dQ2πJ(α;Q)eiQx, (15)

where

 J(α;Q)=∫dxe−iQxeαζ(x)−1α. (16)

Equation (15) also implies that

 ∂ζ(x)∂xeαζ(x)=1α∂eαζ(x)∂x=∫dQ2πiQJ(α;Q)eiQx. (17)

Eqs. (15) and (17) allow the simplification of Eqs. (13) and (14). For eliminating the coefficient, we multiply Eq. (13) by

 (iqx∂ζ(x)∂x+q∥)e−i(qxx+qyy)eq∥ζ(x), (18)

where , and multiply Eq. (14) by . We then use Eqs. (15) and (17), and integrate over . After lengthy calculations we obtain a single equation involving the coefficient only:

 Γ(q,ω)A(qx)=∫∞−∞dPKqy,ω(qx,P)A(P), (19)

where , , , , and

 Γ(q,ω)=[ϵ2−ϵ(ω)]−1{[ϵ2−ϵ(ω)]e−qdf+(ω,q)+ +[ϵ2+ϵ(ω)]eqdf−(ω,q)}, (20)

with the Kernel:

 Kqy,ω(qx,P)=pJ(q−p;qx−P)f−(ω,p)epd(1− −^p⋅^q)−pJ(q+p;qx−P)f+(ω,p)e−pd(1+^p⋅^q). (21)

## Iv Numerical solution

Eq. (19) is a homogeneous integral equation. In the following we discuss the classification of the integral equation as a Fredholm equation of second or third kind depending on the values of and . If the equation:

 Γ(√q2y+q2x,ω)=0, (22)

has real solutions for a real , the integral equation (19) is a Fredholm equation of the third kind. If not, is a Fredholm equation of the second kind. The curve separates the two regimes (see Fig. 2): the continuous solution and the localized plasmon states (see next Section for details).

### iv.1 Continuous solution

We first assume that Eq. (22) has solutions for real . In this condition Eq. (19) has not regular solutions, however considering a generalized functional space, the solution has the form Bart and Warnock (1973):

 A(qx)=α1δ(qx−q0x)+α2δ(qx+q0x)+Areg(qx), (23)

where is the solution of Eq. (22) and is the regular part of the . The coefficients and are determined putting Eq. (23) back to Eq. (19) and making . We can separate the odd and even solutions by making . In the case that we set the protrusion to zero (), the regular part , and we simply recover the solution for the protrusion-free problem. Eq. (23) can be interpreted as the sum of the propagating waves from the protrusion-free problem plus a term that comes from the geometric effect of the protrusion.

The integral equation satisfied by the regular part of the solution is obtained substituting Eq. (23) back to Eq. (19):

 α1Kqy,ω(qx,q0x)+α2Kqy,ω(qx,−q0x)+ +∫∞−∞dPKqy,ω(qx,P)Areg(P)=Γ(q,ω)Areg(qx), (24)

where we used the property obtained by the construction of the function (23):

 Γ(√q2y+q2x,ω)A(qx)=Γ(√q2y+q2x,ω)Areg(qx). (25)

Now, making , the coefficients satisfy the system of equations:

 α1Kqy,ω(q0x,q0x)+α2Kqy,ω(q0x,−q0x)+ +∫∞−∞dPKqy,ω(q0x,P)Areg(P)=0, (26a) α1Kqy,ω(−q0x,q0x)+α2Kqy,ω(−q0x,−q0x)+ +∫∞−∞dPKqy,ω(−q0x,P)Areg(P)=0, (26b)

using the following kernel property and that we obtain:

 ∫∞−∞dP[Kqy,ω(qx,P)−Kqy,ω(q0x,P)]Areg(P)= =Γ(q,ω)Areg(qx), (27)

so, for a given and , one has to solve (27) to obtain the field in the presence of the protrusion. Note that we have a continuous set of frequencies in this case.

### iv.2 Localized plasmons

In the case where Eq. (22) has no real solutions for , Eq. (19) is a homogeneous Fredholm equation of the second kind. This equation, for a given , has solutions for some particular values of . In the following, we consider the case where , that is, the dielectric function is independent of the frequency. In this case, Eq. (19) can be rewritten as:

 λ(ω)D1(qx,qy)A=D2(qx,qy)A, (28)

where we have the following integral operators:

 D1(qx,qy)A=q2[ε2εsinh(qd)+cosh(qd)]A(qx)+ +ε2−ε2ε∫∞−∞dP2πp[J(q+p,qx−P)e−pd(q2y+pq+ +Pqx)+J(q−p,qx−P)epd(q2y−pq+Pqx)]A(P), (29a) D2(qx,qy)A=q2ε[(ε2−ε)(ε−ε1)e−qd+(ε2+ε)(ε+ +ε1)e+qd]A(qx)+ε2−ε2ε∫∞−∞dP2π[(ε−ε1)× ×J(q+p,qx−P)e−pd(q2y+pq+Pqx)+ +(ε+ε1)J(q−p,qx−P)epd(q2y−pq+Pqx)]A(P), (29b)

and we defined:

 λ(ω)=−iσ(ω)ω, (30)

that has the dimension of length. To obtain those results we used:

 f+(ω,p)=ε−ε12ε+λ(ω)p2ε, (31a) f−(ω,p)=ε+ε12ε−λ(ω)p2ε. (31b)

The integral operators and do not depend on the frequency . To proceed, we discretize the integrals in Eqs. (29a) and (29b). First we apply a cutoff in the momentum : , and is chosen to be large enough such that the solution converges (we checked that all boundary conditions are obeyed by the numerical solution). The integral can be discretized by applying Gauss-Legendre quadrature. This will reduce the integral [Eq. (28)] to a generalized eigenvalue problem:

 λ(ω)D1a=D2a, (32)

where are matrix, with the number of Gauss points, and is a vector with dimension that represents the discretized version of the function . Solving Eq. (32) we have the spectrum of eigenvalues . The plasmon frequency then is given by the solution of

 λ(ω)=λn(qy). (33)

If we assume that the conductivity of graphene is given by the Drude formula:

 σ(ω)=σ04πEFℏω+iℏγ, (34)

with , the Fermi energy, and the relaxation rate, the plasmon dispersion is given by:

 ωn(qy)=√4cαEFℏ1λn(qy)−iγ2, (35)

with the fine structure constant and the speed of light. Recalling the condition for the existence of localized plasmons, the solution of equation is:

 λ(ω)=b(qy)qy, (36)

with:

 b(qy)=ε(ε2+ε1)+(ε2+ε1ε2)tanh(qyd)ε+ε2tanh(qyd), (37)

and using the definition of , Eq. (30):

 ωspp(qy)=√4cαEFℏqyb(qy)−iγ2. (38)

Using the condition that, for the localized plasmons, , we arrive at

 λn(qy)qy>b(qy). (39)

The latter relation defines the region of existence of SPP and does not depend on properties of the graphene sheet, that is, it is a purely geometric condition.

Once we compute the coefficient , the coefficients and are calculated using Eqs. (11a) and (11b):

 B(kx) =A(kx)ϵ(ω)−ϵ1−κ2ϵ(ω), (40a) C(kx) =A(kx)ϵ(ω)+ϵ1+κ2ϵ(ω). (40b)

The equation for can be obtained from the boundary condition (13), using the same procedure that was used to obtain the Eq. (19):

 qe−qdD(qx)+∫∞−∞dP2πJ(q+p,qx−P)e−pd(q2y+pq+qxP)D(P)=q(ε−ε12εe−qd+ε+ε12εeqd− −qλ(ω)sinh(qd)ε)A(qx)+∫∞−∞dP2π[(q2y+pq+Pqx)f+(ω,p)J(q+p,qx−P)e−pd+ +(q2y−pq+Pqx)f−(ω,p)J(q−p,qx−P)epd]A(P), (41)

which can be written in a matrix from as:

 G1d=G2a, (42)

where and are the discrete versions of the operators appearing in Eq. (41) and is the corresponding discretized vector. From the previous equation, have the elementary solution:

 d=G−11G2a, (43)

which, since has been previously obtained, readily gives the values for by matrix multiplication.

## V Results

From here on, we consider that the protrusion/indentation is described by a Gaussian profile:

 ζ(x)=−ζ0e−4x2/R2, (44)

with the sign of defining the two different cases schematically illustrated in Figs. 1(A) and 1(B). In Appendix B we calculate the function [Eq. (16)] for the Gaussian profile.

From here on, we will consider the localized plasmons case only, that is, for a given , the maximum frequency that we consider is given by Eq. (38). Otherwise specified, we use the following parameters: , , , eV, , m, m, and m. The numerical parameters are Gauss numbers and the cutoff . These parameters illustrate the implications of the method, but choosing other values amounts to quantitative changes only.

### v.1 Parity

The Kernel of the integral Eq. (19) obeys the identity:

 K(qx,P)=K(−qx,−P), (45)

and the function is even in the variable. From this condition the solutions can be classified in odd and even. Therefore, the limits of integration can be changed as: , which simplifies the numerical solution.

### v.2 Scale invariance

Here we consider how the spectrum changes upon a scale transformation. Making the scale transformation: , , and makes the Kernels of Eq. (28) transform Eq. (29a) to and Eq. (29b) . Therefore the eigenvalue of the matrix equation (28) transform to . From Eq. (35) the frequency of the localized plasmon scale as . This simple transformation of the eigen-frequencies upon a scale transformation is due to the electrostatic limit we have considered from the outset.

From this discussion, only the ratios and matters for the calculation of the dispersion relation. In Fig. 3 we show the dependence on the plasmon frequency for fixed and as function of , where we can see clearly two regimes: for , we have a fast change in the localization frequency starting from the continuous solution (maximum localized frequency), and for the system reaches a plateau and the change in the plasmon frequency is negligible.

### v.3 Discussion

First we show in Fig. 4 the solution for the generalized eigenvalue problem (28) for the first even and odd solutions and , where we can see that the functions approach zero for .

Using Eqs. (11a), (11b), and (41) we can compute all the other functions , , and . The potential field can be calculated now from [see Eqs. (1), (2),and (3)]:

 ϕ1(x,0,z,0) =∫∞−∞dP2πA(P)eiPxe−pz, (46a) ϕc(x,0,z,0) =∫∞−∞dP2πeiPx(B(P)e−pz+C(P)epz), (46b) ϕ2(x,0,z,0) =∫∞−∞dP2πD(P)eiPxepz, (46c)

where for simplicity we are only interested for the results to and , because of time and translation invariance. The electric field can be obtained from , with , where labeling the three regions. The other two components are:

 Ex,1(x,y=0,z,t=0) =−i∫∞−∞dP2πPA(P)eiPxe−pz, (47a) Ez,1(x,y=0,z,t=0) =∫∞−∞dP2πpA(P)eiPxe−pz, (47b)

and similar expressions for the regions and c. First we note that from the parity symmetry of the system, the normalization of the field can be choose such that will be always a real quantity. With a real , the electric fields and will also be real and will be a pure imaginary quantity, i.e., it will always be out-of-phase by with the other electric field components.

From Eqs. (46a)–(47b) we calculate the potential and electrical field in Fig. 5, for the first even solution and in Fig. 6 for the first odd solution in a Gaussian 1D protrusion. For those solutions the plasmon frequencies are THz and THz, respectively. The even solution has a node at , as it should be, and the field strength, as can be seen in the panel D, is concentrated in the axis for and in the axis around , far below the wavenumber m for the light in air. We have verified that the fields obtained by Eqs. (46a)–(47b) satisfy all the boundary conditions. In Figs. 7 and 8 we show the localized plasmons in a groove (Gaussian indentation), where we can see that the field is less localized in the axis in comparison with the protrusion case. However, the indentation “squeezes” the plasmon in the central region. A remarkable characteristic of the electrostatic approximation is that all the fields profile are only a geometric solution of the integral equation, that is, they do not depend on the properties of the graphene sheet. However, they can only exist if Eq. (33) has solution. Therefore, without the graphene sheet there are no localized plasmons.

We also note that the Fermi-energy can be used to tune the frequency of the localized plasmons as per Eq. (35). Another important result of our study is that the localized plasmon dispersion is always below the usual SPP dispersion (see Fig. 2). Therefore, for a given frequency, the wavelength of the localized plasmons is always smaller than its value for a continuous graphene sheet on a homogeneous dielectric, implying a higher degree of confinement of the plasmons in dielectric with a protrusion/indentation.

The charge density can be calculated using the equivalent of Eq. (52):

 n2D(x)=−λ(ω)∫∞−∞dP2π(P2+k2y)A(P)eiPx, (48)

and we show the charge density for the first two modes in the groove in Fig. 9, where we can see again that the charge is localized around the center of the wedge. Finally, we note that we have used a Gaussian profile but our approach can be used to any differentiable profile.

## Vi Conclusions

In this paper we have developed an approach of creating transversely localized plasmons in a flat graphene sheet. This is possible in a configuration where graphene rests on a flat substrate with the opposite surface of the latter showing a protrusion or and indentation (a defect). The localized plasmons dispersion relation appears below the dispersion relation of the propagating plasmons when graphene rests on a flat dielectric of thickness . Above this latter dispersion relation, we have found a continuum of states, which would be needed for describing scattering by the defect. Therefore, we have shown that a defect (in this case with even symmetry) can trap localized surface plasmons. Since the defect is 1D, the wave number along the axis of symmetry of the defect is well defined and, therefore, this defect can also act as channel for propagation of the transversely localized surface plasmons. This geometry has the advantage of being unnecessary to pattern the graphene sheet, therefore it works without deteriorating the electronic mobility of graphene. The generalization of the problem dealt in this paper to a 2D defect is straightforward, involving only extra computer power. This is no impediment as our codes are fast enough and run in a laptop in only few minutes.

###### Acknowledgements.
A.J.C. acknowledges for a scholarship from the Brazilian agency CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico). N.M.R.P. acknowledges the European Commission through the project “Graphene-Driven Revolutions in ICT and Beyond” (Ref. No. 696656) and the Portuguese Foundation for Science and Technology (FCT) in the framework of the Strategic Financing UID/FIS/04650/2013. D. R. C. and G. A. F. acknowledges CNPq under the PRONEX/FUNCAP grants and the CAPES Foundation.

## Appendix A A planar graphene sheet

Let us assume a graphene sheet located at in the plane. The graphene is capped by two dielectrics with dielectric functions , for , and , for . We want to find the spectrum of graphene SPP. The solution of Laplace’s equation for reads

 ϕ1(ρ,z,t)=A1eik∥⋅ρe−k∥ze−iωt≡ϕ1(ρ,z,ω)e−iωt, (49)

where and , and for it is given by

 ϕ2(ρ,z,t)=A2eik∥⋅ρek∥ze−iωt≡ϕ2(ρ,z,ω)e−iωt. (50)

The boundary conditions are

 ϕ1(ρ,0,t) =ϕ2(ρ,0,t) (51a) ϵ1∂ϕ1(ρ,0,t)∂n −ϵ2∂ϕ2(ρ,0,t)∂n=−n2D(ρ,0,t)ϵ0, (51b)

where is the charge density in graphene, whose time dependence can be explicitly made as . The first boundary condition expresses the continuity of the electrostatic potential and the second one the discontinuity of the normal component of the displacement vector. In addition, the electronic density obeys the continuity equation in frequency space: , where . Since the electric current density obeys Ohm’s law, , it follows that

 iωn2D(ρ,0,ω) =−σ∇22Dϕ(ρ,0,ω)=σk2∥ϕ(ρ,0,ω). (52)

Finally, we have for the 2D electronic density the result

 n2D(ρ,0,ω)=−iσωk2∥ϕ(ρ,0,ω) (53)

The first boundary condition implies and the second boundary condition gives

 −ϵ1k∥−ϵ2k∥=iσωϵ0k2∥ (54)

or

 ϵ1k∥+ϵ2k∥+iσωϵ0=0, (55)

which is the condition giving the dispersion relation of the SPP in graphene, Note that for we recover the condition giving the dispersion of SPP at the interface between two dielectrics. In general, we should have written the electrostatic potential as

 ϕ1(ρ,z,t) =∫dk∥(2π)2A1(k∥)eik∥⋅ρe−k∥ze−iωt, (56)

and an identical expression for , except for the dependence which should by replaced by . This way of writing the electrostatic potential is appropriate for discussing rough surface and defects.

## Appendix B The case of a Gaussian profile

In this appendix we give the evaluation of the function for the Gaussian profile. This is accomplished by expanding the exponential of in the integrand, that is,

 J(α;Q) Extra open brace or missing close brace Extra open brace or missing close brace (57)

Therefore we need to compute the integral (the Fourier transform of a Gaussian)

 I(n;Q) =(−ζ0)n∫∞−∞dxe−iQxe−4nx2/R2 =(−ζ0)n