# 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) |

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

(2) |

Finally, in the third region we have

(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:

(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

(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

(6) |

At the interface we have simply

(7) |

Thus the boundary conditions are

(8a) | ||||

(8b) |

which turn into

(9a) | ||||

(9b) |

where we defined:

(10) |

and the solution for and reads

(11a) | ||||

(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

(12a) | ||||

(12b) |

The boundary condition (12a) is simply given by

(13) |

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

(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

(15) |

where

(16) |

Equation (15) also implies that

(17) |

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

(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:

(19) |

where , , , , and

(20) |

with the Kernel:

(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:

(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):

(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):

(24) |

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

(25) |

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

(26a) | |||

(26b) |

using the following kernel property and that we obtain:

(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:

(28) |

where we have the following integral operators:

(29a) | |||

(29b) |

and we defined:

(30) |

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

(31a) | |||

(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:

(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

(33) |

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

(34) |

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

(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:

(36) |

with:

(37) |

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

(38) |

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

(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):

(40a) | ||||

(40b) |

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

(41) |

which can be written in a matrix from as:

(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:

(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:

(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:

(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)]:

(46a) | ||||

(46b) | ||||

(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:

(47a) | ||||

(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):

(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

(49) |

where and , and for it is given by

(50) |

The boundary conditions are

(51a) | ||||

(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

(52) |

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

(53) |

The first boundary condition implies and the second boundary condition gives

(54) |

or

(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

(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,

(57) |

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