## Abstract

The boundary integral equation (BIE) method ascertains explicit relations between localized surface phonon and plasmon polariton resonances and the eigenvalues of its associated electrostatic operator. We show that group-theoretical analysis of the Laplace equation can be used to calculate the full set of eigenvalues and eigenfunctions of the electrostatic operator for shapes and shells described by separable coordinate systems. These results not only unify and generalize many existing studies, but also offer us the opportunity to expand the study of phenomena such as cloaking by anomalous localized resonance. Hence, we calculate the eigenvalues and eigenfunctions of elliptic and circular cylinders. We illustrate the benefits of using the BIE method to interpret recent experiments involving localized surface phonon polariton resonances and the size scaling of plasmon resonances in graphene nanodiscs. Finally, symmetry-based operator analysis can be extended from the electrostatic to the full-wave regime. Thus, bound states of light in the continuum can be studied for shapes beyond spherical configurations.

## 1. Introduction

Materials with negative permittivity allow light confinement to the subdiffraction limit and field enhancement at the interface with ordinary dielectrics [1].

At optical frequencies, noble metals exhibit such behaviour, leading to numerous applications of the field known as plasmonics [2]. These applications include photonic circuits [2] and biosensing [3]. Metals have been considered promising, because, with the advent of nanotechnology, a large variety of structures, shapes and sizes can be tailored to obtain large tunability of plasmon resonances from ultraviolet (UV) to visible and near infrared (IR) [4].

Nowadays, plasmonics has reached a level of maturity such that new plasmonic materials are required for real applications [5]. Plasmon excitations belong to the class of polaritons which are mixed states made from an elementary excitation (i.e. dipole-active such as a phonon, plasmon, magnon, exciton) coupled to a photon [6]. Hence, plasmon polaritons are mixed states made from collective excitations of free electrons in metals and photons, whereas phonon polaritons mix phonons and photons [6].

When the interfaces of these materials are bounded, the surface polaritons become localized. Localized surface phonon (LSPh) polaritons are found with surface optical phonons, which were first studied with dielectric continuum models [7,8]. Later, with the progress made in nanofabrication techniques, surface optical phonons were studied owing to their influence on the electron–phonon coupling in layered nanostructures [9–13]. Mutschkel *et al*. [14] analysed the infrared properties of SiC particles of various polytypes, observing size- and shape-dependent resonances analogous to plasmon polaritons. Like plasmons, phonon polaritons can be used for applications regarding sensing [15], near-field optics [16] or superlensing [17], and also for thermal coherent infrared emission [18] and enhanced energy transfer [19]. Moreover, in the last few years, there has been great interest in phonon polaritons owing to their intrinsically low losses in comparison with plasmon polaritons with many experiments probing the phonon polariton properties of micro- and nanopatterned arrays of polar semiconductors [20–23].

Polaritons can be successfully described with dielectric continuum models by integrating Maxwell equations or the Laplace equation in the quasi-static approximation [24]. The Laplace equation suffices to describe polaritons if the size of the particle is much smaller than the wavelength of the incident light such that the electric field associated with the incident light is spatially uniform and retardation effects can be neglected. In fact, the first treatments of surface phonons and LSPh polariton resonances were quasi-static [7,8] which, as will be seen in the following, are simpler because they provide general properties for arbitrary-shaped particles [8].

Based on potential theory [25,26], a boundary integral equation (BIE) method was developed in [8]. The same theory was invoked independently for plasmon resonances [27] and it reached wide recognition when it was used for localized plasmon resonances in metallic nanoparticles [28]. Needless to say, the formalism was applied also to the radiofrequency properties of living cells [29,30]. Even though in most applications only a purely numerical treatment is possible, the BIE method provides us with direct information about the modes (the resonance strength and frequency) and their relations with the particle shape and dielectric permittivity [8,27–30].

In [7] as well as in [31], Englman and Ruppin used the separation of variables method to solve the Laplace equation for spherical and cylindrical shapes together with their corresponding shells, as well as the slab geometry as a limiting case for a cylindrical shell. The BIE method in a separation of variables scheme was used for quantum wires of elliptical and circular sections [32] and quantum dots of spheroidal shape [33]. Later, the Laplace equation was solved with the separation of variables method for spheroidal [34] and ellipsoidal [35] shapes.

The BIE method is associated with a non-symmetric boundary integral operator that is called the electrostatic operator, whereas its adjoint is called the Neumann–Poincaré operator [26]. The electrostatic operator is well behaved in the sense that, if the boundary of the particle/domain is smooth, the operator is compact and its spectrum has discrete real eigenvalues [25,26]. This operator can be made self-adjoint using the Plemelj symmetrization principle which, from a practical point of view, relates the eigenfunctions of the electrostatic operator with those of its adjoint [36]. In the BIE method, the eigenvalues of the electrostatic operator are related to plasmon resonances [28], with a straightforward eigenvalue–resonance relationship for Drude metals [36,37]. The mathematical literature offers us the whole set of eigenvalues and eigenfunctions for discs, ellipses and spheres [26], as well as for spheroids [38–40] and ellipsoids [41].

In this paper, similar to localized surface plasmon resonances (LSPRs) in metallic nanoparticles [36,37], we establish explicit relationships between the localized surface phonon resonances (LSPhRs) and the eigenvalues of the electrostatic operator and dielectric properties of nanocrystals.

While the mathematical literature offers analytic forms of the eigenvalues and eigenfunctions for various shapes, we show a way of obtaining the spectral properties of the electrostatic operator from direct resolution of the Laplace equation performed on specific physical problems regarding surface phonon or plasmon polaritons. Combining the BIE method with the separation of variables, we recover all the eigenvalues embodied in the energies of surface phonons given in [7,31–35] and in plasmon resonance frequencies found in [42–46]. Thus, we unify all the scattered results from the literature, and we put them in a general setting that can be used to interpret easily and intuitively the experiments concerning LSPh polaritons.

Additionally, the symmetry analysis of the Laplace equation [47,48] has been rarely, if ever, exploited explicitly to characterize the electrostatic operator and its adjoint, which are key elements of the resolution of the Laplace equation. Essentially, we use the same separation of variables method to expand the free-space Green’s function (also called the fundamental solution in the mathematical literature) in an appropriate eigenfunction system. From Lie group-theoretical analysis, it is known that the solutions of the Laplace equation for each separable coordinate system are the eigenfunctions of a pair of two commuting operators [47–49]. It turns out that these functions are also the eigenfunctions of the adjoint of the electrostatic operator and are proportional to the eigenfunctions of the electrostatic operator itself. Therefore, we can also recover the eigenvectors for spheres, prolate and oblate spheroids, and ellipsoids. In addition,we calculate both the eigenvalues and eigenfunctions for circular and elliptic cylinders.

We further analyse the interaction of IR light with an array of GaN micro-discs studied in [22], and we show that the results of the BIE method may offer more information than the calculations obtained with finite-element methods. With the BIE method, we also provide another interpretation of the scaling law of the plasmon resonance frequency with the size of graphene nanodiscs [50].

Finally, there is a new interest in spectral properties of the electrostatic operator regarding some specific applications such as cloaking by anomalous localized resonance [51–54], inverse problems [55,56] and materials characterization [57]. We discuss how our results fit into these new areas of interest and how symmetry-based analysis may be expanded from the electrostatic regime to the full-wave treatment of polaritons with the possibility of studying localized states of light in the continuum [58–60].

The paper is structured as follows. In the following section, we show the way by which the surface phonon and plasmon polariton resonances are calculated from the eigenvalues of the electrostatic operator. In §3, we present two procedures for calculating the eigenvalues and the eigenfunctions of the electrostatic operator for shapes generated by separable coordinate systems. One procedure relies on the solutions of the Laplace equation in separable coordinate systems, and the other is based on the symmetry properties of the Laplace equation applied to the electrostatic operator. In §4, we discuss our results in the context of current research interests and future developments, and in §5 we conclude our work.

## 2. Surface plasmon versus surface phonon polariton resonances

Continuum models describe successfully both plasmon and phonon polaritons. They are based on the same boundary conditions expressing the behaviour of optical modes at the interfaces. The boundary conditions are the continuity of the electric potential and of the normal component of electric induction. The model is schematically shown in figure 1, where a dielectric domain of volume *V* and relative complex permittivity *ϵ*_{1} is bounded by surface *Σ* and is surrounded by a medium of relative complex dielectric permittivity *ϵ*_{0}. In the quasi-static limit, the applied field is uniform of strength **E**_{0} and its time dependence is assumed harmonic, i.e. a complex factor

In the local approximation, the response of the system is given by a surface polarization charge density inducing a nanoparticle polarizability which has an eigenmode decomposition given by [36,37]
*χ*_{k} are the eigenvalues of the electrostatic operator
*w*_{k} are subunitary numbers determined by *E*_{0} and the eigenfunctions of *λ*=(*ϵ*_{1}−*ϵ*_{0})/(*ϵ*_{1}+*ϵ*_{0}). In equation (2.2),
*n*_{x} is the normal derivative. The electrostatic operator and its adjoint are defined on the Hilbert space *L*^{2}(*Σ*) of square-integrable functions on *Σ*. The extinction cross section is related to nanoparticle polarizability by [61]

If we assume that the dielectric permittivity *ϵ*_{0} of the embedding medium is a real constant *ε*_{d}, and the metallic nanoparticle is a Drude metal whose permittivity is
*α*_{LSPR} characterizing the LSPRs and the eigenvalues of *ε*_{m}=*ε*_{d}=1 then *w*_{k}.

A similar expression can be obtained for a void in which *ϵ*_{0} takes a Drude form (2.4) and *ϵ*_{1} is a real constant. The expressions for voids are obtained by the swap *w*_{k} with −*w*_{k}. Often in numerical simulation, instead of a Drude metal, experimental values are considered for metals such as gold and silver [61]. The main difference between the experimental values and the Drude model of dielectric permittivity is the interband contributions that become significant for wavelengths below 500 nm.

In polar crystals, the frequency band between the frequency of transverse optical phonons (*ω*_{T}) and the frequency of transverse optical phonons (*ω*_{L}) is called the reststrahlen band, where the dielectric permittivity is negative and can be modelled as a Lorentz oscillator [62]
*ω*^{2}→*ω*^{2}−*ω*^{2}_{T}, *iγω*→−*iΓω*. Thereby, the counterpart of equation (2.6), which is the polarizability *α*_{LSPhR} characterizing the LSPhRs, takes the form
*ε*_{0} is replaced by equation (2.8), *ε*_{1} is a real constant denoted by *ε*_{d}, and
*χ*_{k},*w*_{k}} which provides the localized surface plasmon/phonon polariton spectrum with the help of equations (2.6)–(2.11).

## 3. Eigenvalues and eigenfunctions of the electrostatic operator by the separation of variables method

### (a) Retrieval of eigenvalues from the solutions of the Laplace equation

As we have already mentioned the whole set of eigenvalues and eigenfunctions of the electrostatic operator is known for discs, ellipses (see also [63] for an earlier account), spheres [26], spheroids [38–40] and ellipsoids [41]. These calculations are based on the expansion of the free-space Green’s function in a separated form of the corresponding coordinate systems, i.e. spherical, spheroidal and ellipsoidal coordinate systems [64]. The separation of variables problem dates back to the end of the nineteenth century [65,66] when the BIE method was laid out on solid ground (a historic account can be found in [67]) with important contributions made by Robertson [68] and Eisenhart [69,70]. The scalar Helmholtz equation and subsequently the Laplace equation are separable in 11 coordinate systems, while the Laplace equation is *R*-separable in six other coordinate systems (see [47,48] for a detailed discussion on the separability and *R*-separability of the Laplace equation). These coordinates are called cyclidic coordinates, which are quartic surfaces [48,64]. In [31–35], the separation of variables method is invoked to solve the Laplace equation associated with a specific physical problem. We consider a separable coordinate system given by coordinates (*ξ*_{1},*ξ*_{2},*ξ*_{3}) and the surface *Σ* is defined by the equation *ξ*_{3}=*ξ*_{0}, where *ξ*_{0} is a constant. Different values of *ξ*_{3} generate confocal surfaces. The boundary conditions are determined by the continuity of the normal component of electric induction ** D** at the surface

*ξ*

_{3}=

*ξ*

_{0}, which reads [32–35]

*f*

_{lm}(

*ξ*

_{0}) is the ratio of the logarithmic derivatives of two functions,

*P*

_{lm}(

*ξ*

_{3}) and

*Q*

_{lm}(

*ξ*

_{3}). The indices

*l*and

*m*define these functions, which are factors of Laplace equation solutions in the separable coordinate system [32–35]. We analyse equation (3.1) from a BIE point of view. A close inspection of the Laplace equation in separable coordinate systems shows that, for each

*l*and

*m*, the solution of the Laplace equation is also proportional to the eigenfunction

*u*

_{lm}(

*ξ*

_{1},

*ξ*

_{2}) of

**reduces to**

*D**χ*

_{lm}is the eigenvalue of

*u*

_{lm}(

*ξ*

_{1},

*ξ*

_{2}). Relation (3.2) has a general validity for arbitrary shapes and it may be used not only for surface phonon polaritons but also for surface plasmons. Combining (2.8) and (3.2), we may obtain the resonance frequency of surface phonon polaritons not only when the surface is described by a separable coordinate system, but also in general for arbitrary shape

*ε*

_{0}is the static dielectric constant of the crystal which is linked to

Equation (3.3) is another form of equation (2.10) via the Lyddane–Sachs–Teller relation. From (3.3), it is easy to check that if there are shapes with

### (b) The eigenvalues and the eigenfunctions of the electrostatic operator from symmetry analysis

In this section, we present a symmetry-based procedure by which the eigenvalues and the eigenfunctions of the electrostatic operator are calculated for various shapes originating from admitting separable coordinate systems. According to the programme linking symmetry with separation of variables, the solutions for each coordinate system are common eigenfunctions of a pair of commuting operators [48]. It turns out that the eigenfunctions of the electrostatic operator are generated by the eigenfunctions of these two commuting operators. This fact can be seen by expressing the free-space Green’s function (2.3) in a separated form, which is determined by three second-order differential equations and two separation constants. The eigenvalues of the electrostatic operator are given by the third differential equation in the third coordinate *ξ*_{3}. The third differential equation, an ordinary second-order equation, has two solutions, of which one is regular in origin (inside the surface) and the other is regular at infinity or outside the surface. The two solutions for *ξ*_{3}=*ξ*_{0} will provide the eigenvalues. This recipe is applied in the following. We designate by *l* and *m* the indices defining the eigenvalues of the two commuting operators. These indices also define the eigenvalues of the electrostatic operator. Let *P*_{lm}(*ξ*_{3}) and *Q*_{lm}(*ξ*_{3}) be, respectively, the regular solution in origin and at infinity for the third differential equation. According to [32–35], *f*_{lm}(*ξ*_{0}) is the ratio of logarithmic derivatives of *P*_{lm}(*ξ*_{3}) and *Q*_{lm}(*ξ*_{3}). Then with the help of (3.1) and (3.2), we may easily obtain the eigenvalues as
*W*{*P*_{lm},*Q*_{lm}}=*P*_{lm}*Q*_{lm}′−*P*_{lm}′*Q*_{lm} is the Wronskian [71] of *P*_{lm}(*ξ*_{3}) and *Q*_{lm}(*ξ*_{3}), all evaluated at *ξ*_{3}=*ξ*_{0}. Hence direct resolution of the Laplace equation performed in [32–35] leads to the full set of eigenvalues of

The explicit expressions of the eigenvalues and the eigenfunctions are given in the mathematical literature for spheres, spheroids and ellipsoids [38–41]. These calculations are based on the expression of the Green’s function (2.3) in the separable coordinate systems. A similar approach was used for discs and ellipses in two-dimensional space [63]. In the following, we also provide explicit expressions for the eigenvalues and the eigenfunctions of the electrostatic operator for three-dimensional counterparts of discs and ellipses, i.e. circular and elliptic cylinders, which to our knowledge is new. In principle, we have to solve the equation

#### (i) Sphere

Working in spherical coordinates, *ρ*=*ρ*_{0}, whereas the eigenfunctions of the electrostatic operator are spherical harmonics *Y* _{lm}(*θ*,*φ*) and the eigenvalues are
*l* is a natural number and *m* is an integer with |*m*|≤*l*. The only surface for which the electrostatic operator is symmetric is the sphere in a three-dimensional space [26]. Expression (3.7) can also be deduced using equation (3.2) from [34].

#### (ii) Prolate spheroid

In these coordinate systems, the transformations are: *z*=*aηζ* with *η*≥1 and −1≤*ζ*≤1. Prolate spheroids characterized by equation *η*=*η*_{0} have the eigenfunctions of the electrostatic operator proportional to *e*^{imφ}*P*_{lm}(*ζ*) and its eigenvalues equal to
*P*_{lm} and *Q*_{lm} are the associated Legendre functions, *l* is a natural number and *m* is an integer with |*m*|≤*l*. We considered the definitions of *P*_{lm} and *Q*_{lm} such that their Wronskian is *W*{*P*_{lm}(*η*),*Q*_{lm}(*η*)}=(−1)^{m}(*l*+*m*)!/((*l*−*m*)!(1−*η*^{2})) [72].

#### (iii) Oblate spheroid

The oblate spheroidal coordinates are obtained by replacing *a* with −*ia* and *η* with *iη*. Hence the eigenfunctions of the electrostatic operator are proportional to *e*^{imφ}*P*_{lm}(*ζ*) and the eigenvalues are

#### (iv) Ellipsoid

The ellipsoidal coordinates, appropriate for an ellipsoid with axes *a*>*b*>*c*, are the solutions *θ* of the equation
*μ*∈[−*b*^{2},−*c*^{2}] and *ν*∈[−*a*^{2},−*b*^{2}]. The Cartesian coordinates are then expressed in ellipsoidal coordinates by the following:
*ρ*=*ρ*_{0}. The Laplace equation separates in this coordinate system such that for each ellipsoidal coordinate a Lamé equation [73] of the form
*n* is a natural number, *B*, which is a separation constant, is real [41,74], and *n*, we obtain 2*n*+1 values of *B* denoted as *m*=0,1,……,2*n*. Thus, the solution *n*. This solution is regular; the second solution, which is regular at infinity, is obtained by the standard procedure [71]

The eigenfunctions of the electrostatic operator are proportional to

#### (v) Elliptic cylinder

Such cylinders are suitably described in elliptic cylindrical coordinates: *z*=*z*, with *d* an arbitrary positive constant. The surface of the cylinder is given by equation *u*=*u*_{0}. For convenience, we take *R* and *r* are the lengths of the longer and, respectively, shorter semi-axes of the cross-sectional ellipsis. Accordingly, the constant *u*_{0} is equal to *e*^{ikz}*ce*_{m}(*q*,*v*) or *e*^{ikz}*se*_{m}(*q*,*v*). Here *k* is a real number and *ce*_{m}(*q*,*v*), *se*_{m}(*q*,*v*) are, respectively, the ‘even’ and the ‘odd’ solutions of the Matheiu equation [32,77]
*q*=−*k*^{2}(*R*^{2}−*r*^{2}) and *α* is a separation constant which takes a characteristic value *a*_{m} or *b*_{m} of the *m*^{th} solution *ce*_{m}(*q*,*v*) or *se*_{m}(*q*,*v*), respectively. If the characteristic values of the Mathieu equation are *a*_{m} the eigenvalues are
*Ce*_{m}(*q*,*u*_{0}) and *Fe*_{m}(*q*,*u*_{0}) are the solutions of the corresponding associated Mathieu equation

Here *α* takes the characteristic value *a*_{m}. Similarly, if *α* is the characteristic value *b*_{m} the eigenvalues are
*Se*_{m}(*q*,*u*_{0}) and *Ge*_{m}(*q*,*u*_{0}) the solution of (3.17) when *α* takes the characteristic value *b*_{m}. We denoted by *W*_{e}(*q*) the Wronskian of *Ce*_{m}(*q*,*u*) and *Fe*_{m}(*q*,*u*) and by *W*_{o}(*q*) the Wronskian of *Se*_{m}(*q*,*u*) and *Ge*_{m}(*q*,*u*). The Wronskians of the solutions of both the Mathieu equation and the associated Mathieu equation are constant. On the other hand, the Mathieu equation has a second solution which is non-periodic. We can take the second solution to be odd or even if the first (periodic) solution is even or odd [78]. Thus, the Wronskians of the solutions of the Mathieu equation as well as of the associated Mathieu equation depend only on *q* (see §§28.5.8, 28.5.9 and 28.20.3–28.20.7 of [79]). The eigenvalues expressed by (3.16) and (3.18) are consistent with the surface phonon spectrum given in [32].

The eigenvalues as well as the eigenfunctions of the electrostatic operator for an elliptic cylinder can be conveniently calculated by expressing the free-space Green’s function *G*(**x**,**x**′) in separable elliptic cylindrical coordinates. First, we use the identity *ce*_{m}(*q*,*v*),*se*_{m}(*q*,*v*)} are not only orthogonal on the space of the square-integrable function *L*^{2}(−*π*,*π*), but also complete because they are the solutions of a Sturm–Liouville problem [80]. Hence, the completeness relation on *L*^{2}(−*π*,*π*) can take the form *G*(**x**,**x**′) depending on (*u*,*v*,*z*) and (*u*′,*v*′,*z*′) has the following expression in elliptic cylindrical coordinates:
*u*_{<}(*u*_{>}) is the smaller (larger) of *u* and *u*′, *k* is real and *q*=−*k*^{2}(*R*^{2}−*r*^{2}). The normal derivative for an elliptic cylindrical surface is

Then, it is easy to check that the eigenvalues are those expressed by equations (3.16) and (3.18) and the eigenfunctions of the *k*=*q*=0. So, if we consider *k*=*q*=0, then *a*_{m}=*b*_{m}=*m*^{2},

#### (vi) Circular cylinder

In cylindrical coordinates *ρ*=*ρ*_{0}. The eigenvalues and eigenfunctions of the electrostatic operator can be easily calculated if we express the free-space Green’s function in cylindrical coordinates. The Green’s function has a well-known expression in cylindrical coordinates that is given, for example, in Jackson [82]. In cylindrical coordinates (*ρ*,*θ*,*z*) and (*ρ*′,*θ*′,*z*′), the free-space Green’s function has the following form:
*I*_{m},*K*_{m} are the modified Bessel functions whose Wronskian is *W*{*I*_{m}(*ρ*),*K*_{m}(*ρ*)}=−1/*ρ* and *ρ*_{<}(*ρ*_{>}) is the smaller (larger) of *ρ* and *ρ*’. It is easy to check that the eigenfunctions of *e*^{ikz}*e*^{imθ}, while the eigenvalues are equal to
*k*=0, then
*m*≠0, which are the eigenvalues of a disc in two dimensions [26].

Both types of cylinder have a continuum spectrum, which is a signature of accommodating a running wave along the *z*-direction; because cylindrical surfaces are unbounded, they are, non-compact. Similar to phonon polaritons and implicitly using the separation of variables method, the full spectrum of surface plasmons was calculated for spherical [43], spheroidal [44] and circular cylindrical shapes [45,46] in the so-called hybridization model of plasmon resonances. Analysis of plasmon resonances in [43–46] leads to the same conclusion regarding the eigenvalues of the electrostatic operator.

All these shapes presented above also generate shelled configurations made of confocal surfaces determined by the above coordinate systems. These shelled systems support cavity and particle polaritons that interact among themselves. The interaction between cavity and particle modes however takes place in the subspace determined by each eigenvalue, thus it is possible to calculate all eigenmodes in shelled configurations. Calculations of shelled structures were performed in [7,31,43–46].

## 4. Discussions

The BIE method can be the method of choice in many situations, even when nanoparticles are fabricated from crystals with an anisotropic dielectric constant. In this case, the BIE method can be perfectly adapted by adopting the approach used in [83]. The BIE method can also be a convenient tool to comprehend and calibrate results with mixed states made by coupling continuum and localized states such as those studied in a recent work [23]. Moreover, a BIE method provides more insight about the localized resonance modes. For instance, Feng *et al.* [22] studied the LSPhRs in an array of micro-discs of gallium nitride on a silicon carbide substrate. The effect of LSPh polariton resonances was studied by reflectance measurements in two polarization configurations, transverse electric (TE) and transverse magnetic (TM). The authors observed a little discrepancy between the measurements, which show only one resonance, and the finite-element simulations, which predict two close resonances. We calculated the eigenvalues and their weights for a smooth but also a similar micro-disc that in turn has an aspect ratio of 1/4.7. In the reflectance measurements in the TE configuration the electric field is parallel to the surface of the disc, whereas in the TM configuration both types of modes (i.e. parallel and perpendicular to the disc surface) are excited. Although our calculations are qualitative because they are made for a single disc some insight can be obtained because we have obtained all possible resonances of such a micro-disc. In figure 2, we may see that, for the electric field parallel to the disc surface, there are only two resonances, of which one is predominant. On the contrary, for the electric field perpendicular to the disc, there are several resonances, of which two are predominant. With the BIE method, we may see that the small discrepancy noted in [22] is from two small resonances that come from different field polarizations. These two resonances are highlighted by a dotted rectangle.

Spectral studies of the electrostatic operator have sparked new interest arising from mathematical theory regarding plasmons, cloaking as well as various inverse problems [51–57]. For bounded and smooth surfaces, the electrostatic operator is compact with a countable number of eigenvalues that accumulate at the origin that belongs to the essential spectrum [26]. The eigenvalues lie within

Questions regarding spectral properties such as the location of the essential spectrum of the electrostatic operator and the asymptotic form of its eigenfunctions turn out to be of great interest. Some form of cloaking is related to anomalous localized resonance that takes place at the accumulation point of the eigenvalues [51–54,57]. Several two-dimensional systems have been considered for cloaking, including concentric discs [51–53], ellipses [54] and confocal ellipses [85]. All these approaches need the eigenfunctions in addition to all the eigenvalues. The results of the previous section regarding the eigenvalues and the eigenfunctions of the electrostatic operator defined on different surfaces and on their shelled structure counterparts may be used to extend the analysis of cloaking from discs and ellipses to a three-dimensional setting.

Because 0 belongs to the essential spectrum a legitimate question is whether there are shapes which have 0 as an eigenvalue and more generally whether there is a shape that may have as an eigenvalue any number between *χ*_{lm} with *l*=1, *m*=0,1 for prolate and oblate spheroids. We denote the ratio *a*_{z}/*a*_{x} as the aspect ratio, where *a*_{z} is the semi-axis along the *z*-direction and *a*_{x} is the semi-axis along the *x*-direction. An aspect ratio >1 defines a prolate spheroid and an aspect ratio <1 defines an oblate spheroid. Prolate spheroids reach the value *χ*_{10} is fairly close to *ω*_{T} for phonon polaritons.

It can be seen in figure 3 that there is an oblate spheroid with an aspect ratio between 0.1 and 1 for which *χ*_{10}=0. Furthermore, as the aspect ratio tends to 0, *χ*_{10} approaches *χ*_{11} goes to *ω*_{L}.

In other words, asymptotically for *l*=1 a two-dimensional disc embedded in three dimensions has an eigenvalue *d* is the diameter of the graphene nanodisc [50]. In three-dimensional models, graphene nanodiscs are simulated as finite-thickness discs with a bulk complex conductivity obtained from graphene sheet conductivity divided by the disc thickness *h*. Hence the finite thickness nanodisc has the following dielectric relative permittivity [86]:
*ε*_{v} is the vacuum permittivity. We denoted by *E*_{F} the Fermi energy in a graphene disc, *τ* is the relaxation time and *h*/*d* only. Because the plasmon frequency is *h*/*d* and not quadratically as in the case of oblate discs.

Owing to the success of the BIE in the quasi-static regime, there is a renewed interest in the modal approach based on the BIE to plasmon problems in the full-wave electromagnetic regime [87]. The full-wave electromagnetic treatment discussed in [87] has at its core the free-space Green’s function of the scalar Helmholtz equation that is also separable in those 11 coordinate systems in which the Laplace equation separates [47,48]. In general, in the full-wave regime, owing to radiation losses, the eigenvalues are no longer real as in the case of the quasi-static regime [87]. It will be interesting to apply the same symmetry arguments about the Helmholtz equation to the spectrum of the boundary integral operators used in the full-wave treatment for the shapes and their corresponding shells discussed in this work. In this way, one can study optical bound states [60] with no radiation losses in structures having other forms than spherical shapes which have been considered recently [58,59].

## 5. Conclusion

In this work, we systematically apply the BIE method to obtain explicit relations between the eigenvalues of the electrostatic operator associated with the method and the LSPh/plasmon polariton resonances. The BIE method shows real benefits by explicitly assigning modes, spectral features and resonances to experimental data. Active modes of thin discs are analysed asymptotically with respect to disc thickness. Two examples were considered: the LSPh resonances in GaN micro-discs and the scaling of localized surface plasmon resonance frequency with the size of graphene nanodiscs.

We have shown that, regarding the spectra of both surface phonons and plasmons, all the calculations performed so far in separable coordinates contain not only all the eigenvalues but also all the eigenfunctions of the electrostatic operators.

Group-theoretical analysis applied to the Laplace equation may be exploited in obtaining the eigenvalues and eigenvectors of the electrostatic operator for shapes and their corresponding shell structures described by separable coordinate systems. We calculated explicitly the eigenvalues and the eigenfunctions of elliptic and circular cylinders, which with the known results about other shapes may be used to extend the study of cloaking by anomalous localized resonance from two- to three-dimensional structures. The results of Lie group analysis may be further employed to study the spectrum of boundary operators associated with the full-wave regime. The extension to the full-wave regime creates the opportunity to study and obtain bound states in continuum with shapes different from spherical geometry.

## Data accessibility

All data accompanying this publication are available within the publication and are generated using standard routines in Mathematica^{®} (see electronic supplementary material).

## Authors' contributions

The authors contributed equally to the manuscript.

## Competing interests

We have no competing interests.

## Funding

The work was supported by a grant from the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project no. PNII-ID-PCCE-2011-2-0069.

## Acknowledgements

The authors thank the National Funding Authority for their financial support.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3702211.

- Received October 24, 2016.
- Accepted February 13, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.