Analytical results regarding electrostatic resonances of surface phonon/plasmon polaritons: separation of variables with a twist

TL;DR

This paper uses group-theoretical analysis and boundary integral equations to explicitly determine surface phonon and plasmon polariton resonances for separable shapes, unifying previous results and enabling new applications like cloaking and bound states of light.

Contribution

It introduces a symmetry-based method to calculate eigenvalues and eigenfunctions of the electrostatic operator for separable shapes, extending analysis to full-wave regimes.

Findings

Eigenvalues and eigenfunctions for elliptic and circular cylinders calculated

Unified framework for surface phonon and plasmon polariton resonances

Potential applications in cloaking and bound states of light

Abstract

The boundary integral equation 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 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 the opportunity to expand the study of phenomena like cloaking by anomalous localized resonance. For that reason we calculate the eigenvalues and eigenfunctions of elliptic and circular cylinders. We illustrate the benefits of using the boundary integral equation method to interpret recent experiments involving localized surface phonon polariton resonances and the size scaling of plasmon resonances in graphene nano-disks. Finally, symmetry-based operator analysis can be extended from electrostatic to full-wave regime. Thus, bound states of light in the continuum can be studied for shapes beyond spherical configurations.

Equations

$$\alpha =V\sum _k\frac{w_k}{\frac{1}{2\lambda }-{\chi }_k},$$

$$\hat{M}\left[u\right]=-\underset{{\mathbf{x}}^{\prime }\in \mathbf{\Sigma }}{\int }u\left({\mathbf{x}}^{\prime }\right)\frac{\partial }{\partial {\mathbf{n}}_{\mathbf{x}}}G\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)d{\Sigma }_{{\mathbf{x}}^{\prime }}$$

$$G\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\frac{1}{4\pi }\frac{1}{\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|}$$

$$C_{ext}=\frac{2\pi }{\lambda }Im\left(\alpha \right).$$

$${ϵ}_1={\epsilon }_m-\frac{{\omega }_p^2}{\omega \left(\omega +\text{mbox}i\gamma \right)},$$

$${\alpha }_{LSPR}(\omega )=V\sum _k\frac{w_k({\epsilon }_m-{\epsilon }_d)}{{\epsilon }_k}-\frac{w_k}{1/2-{\chi }_k}\frac{{\epsilon }_d}{{\epsilon }_k}\frac{{\tilde{\omega }}_{pk}^2}{\omega (\omega +i\gamma )-{\tilde{\omega }}_{pk}^2}.$$

$${\tilde{\omega }}_{pk}={\omega }_p\sqrt{\left(1/2-{\chi }_k\right)/{\epsilon }_k}$$

$$ϵ={\epsilon }_{\infty }\left(1+\frac{{\omega }_L^2-{\omega }_T^2}{{\omega }_T^2-{\omega }^2+i\Gamma \omega }\right).$$

$${\alpha }_{LSPhR}(\omega )=V\sum _k\frac{w_k({\epsilon }_{\infty }-{\epsilon }_d)}{{\epsilon }_k}-\frac{w_k}{1/2-{\chi }_k}\frac{{\epsilon }_d}{{\epsilon }_k}\frac{{\tilde{\omega }}_{phk}^2}{\omega (\omega -i\Gamma )-{\omega }_T^2-{\tilde{\omega }}_{phk}^2},$$

$${\tilde{\omega }}_{phk}=\sqrt{{\epsilon }_{\infty }\left(1/2-{\chi }_k\right)\left({\omega }_L^2-{\omega }_T^2\right)/{\epsilon }_k}$$

$${\alpha }_{LSPhR\_v}(\omega )=V\sum _k\frac{w_k({\epsilon }_d-{\epsilon }_{\infty })}{{\epsilon }_{k\_v}}+\frac{w_k}{1/2+{\chi }_k}\frac{{\epsilon }_d}{{\epsilon }_{k\_v}}\frac{{\tilde{\omega }}_{phk\_v}^2}{\omega (\omega -i\Gamma )-{\omega }_T^2-{\tilde{\omega }}_{phk\_v}^2},$$

$$\frac{{ϵ}_1\left(\omega \right)}{{ϵ}_0}=f_{lm}\left({\xi }_0\right).$$

$$\frac{{ϵ}_1\left(\omega \right)}{{ϵ}_0}=-\frac{1/2+{\chi }_{lm}}{1/2-{\chi }_{lm}},$$

$${\tilde{\omega }}_{lm}^2={\omega }_T^2\frac{\left(1/2+{\chi }_{lm}\right){\epsilon }_d+\left(1/2-{\chi }_{lm}\right){\epsilon }_0}{\left(1/2+{\chi }_{lm}\right){\epsilon }_d+\left(1/2-{\chi }_{lm}\right){\epsilon }_{\infty }}.$$

$$\frac{{\epsilon }_0}{{\epsilon }_{\infty }}=\frac{{\omega }_L^2}{{\omega }_T^2}.$$

$${\chi }_{lm}=\frac{1}{2}+\frac{P_{lm}^{\prime }{\left({\xi }_3\right)}_{{\xi }_3={\xi }_0}{Q}_{lm}{\left({\xi }_3\right)}_{{\xi }_3={\xi }_0}}{W{\left[P_{lm}\left({\xi }_3\right),Q_{lm}\left({\xi }_3\right)\right]}_{{\xi }_3={\xi }_0}},$$

$$\Delta G\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=-\delta \left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)$$

$${\chi }_{lm}=\frac{1}{2}\frac{1}{2l+1}.$$

$${\chi }_{lm}=\frac{1}{2}-(-1{)}^m\frac{\left(l-m\right)!}{\left(l+m\right)!}\left({\eta }_0^2-1\right)P_{lm}{\left({\eta }_0\right)}^{\prime }{Q}_{lm}\left({\eta }_0\right)$$

$${\chi }_{lm}=\frac{1}{2}+(-1{)}^m\frac{\left(l-m\right)!}{\left(l+m\right)!}\left({\eta }_0^2+1\right)P_{lm}{\left(i{\eta }_0\right)}^{\prime }{Q}_{lm}\left(i{\eta }_0\right).$$

$$\frac{x^2}{a^2+\theta }+\frac{y^2}{b^2+\theta }+\frac{z^2}{c^2+\theta }=1,$$

$$\begin{array}[]{l}x^{2}=\frac{\left({a^{2}+\rho}\right)\left({a^{2}+\mu}\right)\left({a^{2}+\nu}\right)}{\left({a^{2}-b^{2}}\right)\left({a^{2}-c^{2}}\right)}\\ y^{2}=\frac{\left({b^{2}+\rho}\right)\left({b^{2}+\mu}\right)\left({b^{2}+\nu}\right)}{\left({b^{2}-c^{2}}\right)\left({b^{2}-a^{2}}\right)}\\ z^{2}=\frac{\left({c^{2}+\rho}\right)\left({c^{2}+\mu}\right)\left({c^{2}+\nu}\right)}{\left({c^{2}-a^{2}}\right)\left({c^{2}-b^{2}}\right)}.\\ \end{array}$$

$$4\Delta \left(\rho \right)\frac{d}{d\rho }\left\{\Delta \left(\rho \right)\frac{dL\left(\rho \right)}{d\rho }\right\}=\left(n\left(n+1\right)+B\right)L\left(\rho \right)$$

$$K_n^m\left(\rho \right)=\left(2n+1\right)L_n^m\left(\rho \right)\underset{\rho }{\overset{\infty }{\int }}\frac{d{\rho }^{\prime }}{{\left[L_n^m\left({\rho }^{\prime }\right)\right]}^2\Delta \left({\rho }^{\prime }\right)}.$$

$${\chi }_{nm}=\frac{1}{2}-\frac{\Delta \left({\rho }_0\right)L_n^m{\left({\rho }_0\right)}^{\prime }{K}_n^m\left({\rho }_0\right)}{2n+1}.$$

$$\frac{d^{2}V\left(q,v\right)}{d{v}^2}+\left(\alpha -2q\cos \left(2v\right)\right)V\left(q,v\right)=0,$$

$${\chi }_{qm}=\frac{1}{2}+\frac{C{e}_m{\left(q,u_0\right)}^{\prime }F{e}_m\left(q,u_0\right)}{W_e(q)},$$

$$\frac{d^{2}U\left(q,u\right)}{d{v}^2}-\left(\alpha -2q\cosh \left(2u\right)\right)U\left(q,u\right)=0.$$

$${\chi }_{qm}=\frac{1}{2}+\frac{S{e}_m{\left(q,u_0\right)}^{\prime }G{e}_m\left(q,u_0\right)}{W_o(q)},$$

$$G=\frac{-1}{2{\pi }^2}\underset{-\infty }{\overset{\infty }{\int }}dk{e}^{ik(z-z^{\prime })}\{\sum _{m=0}^{\infty }\frac{1}{W_e(q)}c{e}_m(q,v)c{e}_m(q,v^{\prime })C{e}_m(q,u_{<})F{e}_m(q,u_{>})+$$

Full text

I Introduction

Materials with negative permittivity allow light confinement to sub-diffraction limit and field enhancement at the interface with ordinary dielectrics Barnes et al. (2003). At optical frequencies noble metals exhibit such behaviour leading to numerous applications of the field called plasmonics Ozbay (2006). These applications include photonic circuits Ozbay (2006) and biosensing Brolo (2012). Metals have been promising since, 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) Noguez (2007).

Nowadays, plasmonics has reached the level of maturity such that new plasmonic materials are required for real applications Guler et al. (2015). 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 Agranovich and Mills (1982). Hence, plasmon polaritons are mixed states made from collective excitations of free electrons in metals and photons, while phonon polaritons mix phonons and photons Agranovich and Mills (1982).

When the interfaces of these materials are bounded, the surface polaritons become localized. Localized surface phonon polaritons come along with surface optical phonons which were first studied with dielectric continuum models Englman and Ruppin (1966); Fuchs (1975). Later on, with the progress made in nanofabrication techniques, surface optical phonons were studied due to their influence on the electron-phonon coupling in layered nanostructures Licari and Evrard (1977); Wendler (1985); Trallero-Giner et al. (1992); Nash (1992); Comas et al. (1997). Mutschkel et al. analyzed infrared properties of SiC particles of various polytypes observing size and shape-dependent resonances analogous to plasmon polaritons Mutschke1 et al. (1999). Like plasmons, phonon polaritons can be used for applications regarding sensing Hillenbrand et al. (2002), near-field optics Taubner et al. (2004) or superlensing Taubner et al. (2006), but also for thermal coherent infrared emission Greffet et al. (2002) and enhanced energy transfer Shen et al. (2009). Moreover, in the last few years there is a great interest in phonon polaritons due to their intrinsic low losses in comparison with plasmon polaritons with many experiments probing phonon polariton properties of micro- and nano-patterned arrays of polar semiconductors Chen et al. (2014); Caldwell et al. (2015); Feng et al. (2015); Gubbin et al. (2016).

Polaritons can be successfully described with dielectric continuum models by integrating Maxwell equations or Laplace equation in the quasi-static approximation Tsukerman (2008). 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 localized surface phonon polariton resonances were quasi-static Englman and Ruppin (1966); Fuchs (1975) which, as it will be seen in the following, are simpler by providing general properties for arbitrary shaped particles Fuchs (1975).

Based on potential theory Kellogg (1929); D. Khavison (2007), a boundary integral equation (BIE) method was developed in Ref. Fuchs (1975). The same theory was invoked independently for plasmon resonances Ouyang and Isaacson (1989) and it reached wide recognition when it was used for localized plasmon resonances in metallic nanoparticles Fredkin and Mayergoyz (2003). Needless to say, the formalism was applied also to radiofrequency properties of living cells Vrinceanu and Gheorghiu (1996); Sandu et al. (2010). Even though in most applications only a purely numerical treatment is possible, the BIE method provides direct information about the modes (the resonance strength and frequency) and their relations with particle shape and dielectric permittivity Fuchs (1975); Ouyang and Isaacson (1989); Fredkin and Mayergoyz (2003); Vrinceanu and Gheorghiu (1996); Sandu et al. (2010).

In Ref. Englman and Ruppin (1966) as well as in the following paper Englman and Ruppin (1968) 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 of cylindrical shell. The BIE method in a separation of variables scheme was used for quantum wires of elliptical and circular sections Knipp and Reinecke (1992a) and quantum dots of spheroidal shape Knipp and Reinecke (1992b). Later on, the Laplace equation was solved with the separation of variables method for spheroidal Comas et al. (2002) and ellipsoidal Reese et al. (2004) shapes.

The BIE method is associated with a non-symmetric boundary integral operator that is called the electrostatic operator, while its adjoint is called the Neumann-Poincaré operator D. Khavison (2007). 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 Kellogg (1929); D. Khavison (2007). This operator can be made self-adjoint using Plemelj symmetrization principle which, from practical point of view, relates the eigenfunctions of the electrostatic operator with those of its adjoint Sandu (2013). In the BIE method the eigenvalues of the electrostatic operator are related to plasmon resonances Fredkin and Mayergoyz (2003), with a straightforward eigenvalue-resonance relationship for Drude metals Sandu (2013); Sandu et al. (2011). Mathematical literature offers the whole set of eigenvalues and eigenfunctions for disks, ellipses, and spheres D. Khavison (2007), as well as for spheroids Ahner and Arenstorf (1986); Ahner (1994); Ahner et al. (1994) and ellipsoids Ritter (1995).

In this paper, similar to localized surface plasmon resonances (LSPRs) in metallic nanoparticles Sandu (2013); Sandu et al. (2011) we establish explicit relationships between the localized surface phonon resonances (LSPhRs) and the eigenvalues of the electrostatic operator and dielectric properties of nanocrystals.

While mathematical literature offer analytic forms of the eigenvalues and eigenfunctions for various shapes we will show a way of obtaining the spectral properties of the electrostatic operator from direct resolution of 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 Refs. Englman and Ruppin (1966, 1968); Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004) and in plasmon resonance frequencies found in Refs. Moussiaux et al. (1977); Prodan and Nordlander (2004); Brandl and Nordlander (2007); Morandi (2008); Wan et al. (2013). Thus, we unify all scattered results from literature and we put them in a general setting that can be used to interpret easily and intuitively the experiments concerning localized surface phonon polaritons.

Additionally, the symmetry analysis of Laplace equation Boyer et al. (1976); Miller (1977) has been rarely if ever exploited explicitly to characterize the electrostatic operator and its adjoint which are key elements of the resolution of 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 mathematical literature) in an appropriate eigenfunction system. From Lie group-theoretical analysis it is known that the solutions of Laplace equation for each separable coordinate system are the eigenfunctions of a pair of two commuting operators Boyer et al. (1976); Miller (1977); Makarov et al. (1967). It turns out that these functions are also the eigenfunctions of the adjoint of the electrostatic operator and are proportional to the eigenfunctions of electrostatic operator itself. Therefore, we can also recover the eigenvectors for spheres, prolate and oblate spheroids, and ellipsoids. In addition to that we calculate both the eigenvalues and eigenfunctions for circular and elliptic cylinders.

We further analyze the interaction of IR light with an array of GaN micro-disks studied in Ref. Feng et al. (2015) 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 nano-disks Fang et al. (2013).

Finally, there is a new interest in spectral properties of the electrostatic operator regarding some specific applications like cloaking by anomalous localized resonance Nicorovici et al. (1994); Milton et al. (2005); Milton and Nicorovici (2006); Ando and Kang (2016), inverse problems Ammari et al. (2014, 2016), and materials characterization Ammari et al. (2013). We will discuss how our results fit into these new areas of interest and how symmetry-based analysis may be expanded from electrostatic regime to the full-wave treatment of polaritons with the possibility of studying localized states of light in the continuum Silveirinha (2014); Monticone and A.Alu (2014); Hsu et al. (2016).

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

II Surface plasmon versus surface phonon polariton resonances

Continuum models describe successfully both plasmon and phonons 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 Fig. 1, where a dielectric domain of volume $V$ and relative complex permittivity $\epsilon_{1}$ is bounded by surface $\Sigma$ and is surrounded by a medium of relative complex dielectric permittivity $\epsilon_{0}$. In the quasi-static limit, the applied field is uniform of strength $\bf{E}_{0}$ and its time-dependence is assumed harmonic, i. e., a complex factor $\exp\left({-i\omega t}\right)$ is understood.

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 Sandu et al. (2011); Sandu (2013)

[TABLE]

where $\chi_{k}$ are the eigenvalues of the electrostatic operator

[TABLE]

while $w_{k}$ are subunitary numbers determined by $\bf{E}_{0}$ and the eigenfunctions of $\hat{M}$ and its adjoint, and $\lambda=(\epsilon_{1}-\epsilon_{0})/(\epsilon_{1}+\epsilon_{0})$. In Eq. (2)

[TABLE]

is the free-space Green’s function and $\partial/\partial\bf{n}_{x}$ is the normal derivative. The electrostatic operator and its adjoint are defined on the Hilbert space $L^{2}(\Sigma)$ of square-integrable functions on $\Sigma$. The extinction cross-section is related to nanoparticle polarizability by Maier (2007)

[TABLE]

Also, the near-field enhancement has an eigenmode decomposition similar to Eq. (1) Sandu (2013).

If we assume that the dielectric permittivity $\epsilon_{0}$ of the embedding medium is a real constant $\varepsilon_{d}$ and the metallic nanoparticle is a Drude metal whose permittivity is

[TABLE]

we can obtain compact relations between the polarizability $\alpha_{LSPR}$ characterizing the localized plasmon resonances (LSPRs) and the eigenvalues of $\hat{M}$ Sandu (2013)

[TABLE]

Here $\varepsilon_{k}=(1/2+\chi_{k})\varepsilon_{d}+(1/2-\chi_{k})\varepsilon_{m}$ is an effective dielectric constant and

[TABLE]

is the expression of LSPR frequency. Equations (6)-(7) show clearly the relation between localized surface plasmon resonances and the eigenvalues of $\hat{M}$. For instance, if $\varepsilon_{m}=\varepsilon_{d}=1$ then $\tilde{\omega}_{pk}=\omega_{p}\sqrt{1/2-\chi_{k}}$. On the other hand, the resonance strength is included in the numerical factor $w_{k}$.

A similar expression can be obtained for a void in which $\epsilon_{0}$ takes a Drude form (4) and $\epsilon_{1}$ is a real constant. The expressions for voids are obtained by the swap $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2-\chi_{k}}\right)\leftrightarrow\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2+\chi_{k}}\right)$ and replacing $w_{k}$ with $-w_{k}$. Often in numerical simulation, instead of a Drude metal, experimental values are considered for metals like gold and silver Maier (2007). The main difference between 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 ($\omega_{T})$ and the frequency of transverse optical phonons ($\omega_{L})$ is called the reststrahlen band, where the dielectric permittivity is negative and can be modeled as a Lorentz oscillator Kittel (1996)

[TABLE]

In the reststrahlen band localized surface phonon polaritons of polar nanocrystals have similar behaviour as localized surface plasmons in metallic nanoparticles. Similar to Eqs. (6)-(7) we can derive the phonon polarizability that defines the LSPhRs. The phonon polarizability is obtained by the following algebraic substitutions: $\omega^{2}\to\omega^{2}-\omega_{T}^{2}$, $\omega_{p}^{2}\to\varepsilon_{\infty}\left({\omega_{L}^{2}-\omega_{T}^{2}}\right)$, $\varepsilon_{m}\to\varepsilon_{\infty}$, and $i\gamma\omega\to-i\Gamma\omega$. Thereby, the counterpart of Eq. (6), which is the polarizability $\alpha_{LSPhR}$ characterizing the LSPhRs, takes the form

[TABLE]

where

[TABLE]

is the LSPhR frequency and $\varepsilon_{k}=(1/2-\chi_{k})\varepsilon_{\infty}+(1/2+\chi_{k})\varepsilon_{d}$. Equation (10) illustrates the dependence of LSPhRs on the nanocrystal shape in the same fashion as Eq (7) shows the dependence of LSPR on the metallic nanoparticle shape. For a nanovoid in a polar crystal $\varepsilon_{0}$ is replaced by Eq. (8), $\varepsilon_{1}$ is a real constant denoted by $\varepsilon_{d}$, and

[TABLE]

with $\tilde{\omega}_{phk\_v}=\sqrt{{\varepsilon_{\infty}\left({1/2+\chi_{k}}\right)\left({\omega_{L}^{2}-\omega_{T}^{2}}\right)}/{\varepsilon_{k\_v}}}$ and $\varepsilon_{k\_v}=(1/2-\chi_{k})\varepsilon_{d}+(1/2+\chi_{k})\varepsilon_{\infty}$. We have given explicitly the expressios of LSPhRs for both nanocrystals and nanovoids because often on can encounter both situations like semiconductor nanodots embedded in a crystal host that has its own reststrahlen band. To conclude this section, with the boundary integral equation method we may calculate the eigenvalues and eigenvectors of the electrostatic operator to obtain a set of numbers $\left\{{\chi_{k},w_{k}}\right\}$ which provides the localized surface plasmon/phonon polariton spectrum with the help of Eqs. (6)-(11).

III Eigenvalues and eigenfunctions of the electrostatic operator by separation of variables method

III.1 Retrieval of eigenvalues from the solutions of Laplace equation

As we have already mentioned the whole set of eigenvalues and eigenfunctions of the electrostatic operator is known for disks, ellipses (see also Ahner (1986) for an earlier account), spheres D. Khavison (2007), spheroids Ahner and Arenstorf (1986); Ahner (1994); Ahner et al. (1994) and ellipsoids Ritter (1995). 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 Moon and Spencer (1961). The separation of variables problem dates back to the end of 19th century Stackel (1891); Bocher (1894) when the BIE method was laid out on solid ground (a historic account can be found in Costabel (2007)) with important contributions made by Robertson Robertson (1927) and Eisenhart Eisenhart (1934a, b). The scalar Helmholtz equation and subsequently the Laplace equation are separable in eleven coordinate systems, while the Laplace equation is $R$-separable in other six coordinate systems (see Boyer et al. (1976); Miller (1977) for a detailed discussion on separability and $R$-separability of Laplace equation) These coordinates are called cyclidic coordinates which are quartic surfaces Moon and Spencer (1961); Miller (1977). In Refs. Englman and Ruppin (1968); Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004) the separation of variables method is invoked to solve the Laplace equation associated with a specific physical problem. We consider a separable coordinate system is given by coordinates $\left({\xi_{1},\xi_{2},\xi_{3}}\right)$ and the surface $\Sigma$ is defined by the equation $\xi_{3}=\xi_{0}$, where $\xi_{0}$ is a constant. Different values of $\xi_{3}$ generate confocal surfaces. The boundary conditions are determined by the continuity of the normal component of electric induction D at the surface $\xi_{3}=\xi_{0}$ which reads Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004)

[TABLE]

Here $f_{lm}\left({\xi_{0}}\right)$ is the ratio of logarithmic derivatives of two functions, $P_{lm}\left({\xi_{3}}\right)$ and $Q_{lm}\left({\xi_{3}}\right)$. The indices $l$ and $m$ define these functions which are factors of Laplace equation solutions in the separable coordinate system Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004). We analyze Eq. (12) from BIE point of view. A close inspection of Laplace equation in separable coordinate systems shows that, for each $l$ and $m$, the solution of Laplace equation is also proportional to the eigenfunction $u_{lm}({\xi_{1},\xi_{2}})$ of $\hat{M}$ such that the continuity on the normal component of D reduces to

[TABLE]

where $\chi_{lm}$ is the eigenvalue of $\hat{M}$ corresponding to $u_{lm}\left({\xi_{1},\xi_{2}}\right)$. Relation (13) has a general validity for arbitrary shapes and it may be used not only for surface phonon polaritons but also for surface plasmons. Combining (8) and (13) 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

[TABLE]

Here $\varepsilon_{0}$ is the static dielectric constant of the crystal which is linked to $\varepsilon_{\infty}$ by Lyddane-Sachs-Teller relation Kittel (1996)

[TABLE]

Equation (14) is another form of Eq. (10) via Lyddane-Sachs-Teller relation. From (14) it is easy to check that if there are shapes with $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2-\chi_{lm}}\right)\to 0$ then $\tilde{\omega}_{lm}\to\omega_{T}$. We will see in the next section that infinitely long nanorods fulfill this condition. In the same time, infinitely thin nanodisks have an eigevalue fulfilling the above condition and other satisfying the condition $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2+\chi_{lm}}\right)\to 0$ with the corresponding resonance frequency $\tilde{\omega}_{lm}\to\omega_{L}$.

III.2 The eigenvalues and the eigenfunctions of the electrostatic operator from symmetry analysis

In this subsection we will present a symmetry-based procedure by which the eigenvalues and the eigenfunctions of the electrostatic operator are calculated for various shapes originating from separable coordinate systems admitting. According to the program linking symmetry with separation of variables the solutions for each coordinate system are common eigenfunctions of a pair of commuting operators Miller (1977). 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 (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 $\xi_{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 $\xi_{3}=\xi_{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}\left({\xi_{3}}\right)$ and $Q_{lm}\left({\xi_{3}}\right)$ be, respectively, the regular solution in origin and at infinity for the third differential equation. According to Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004) $f_{lm}\left({\xi_{0}}\right)$ is the ratio of logarithmic derivatives of $P_{lm}\left({\xi_{3}}\right)$ and $Q_{lm}\left({\xi_{3}}\right)$. Then with the help of (12) and (13) we may easily obtain the eigenvalues as

[TABLE]

where the prime means the derivative and $W\left\{{P_{lm},Q_{lm}}\right\}=P_{lm}Q_{lm}^{\prime}-P_{lm}^{\prime}Q_{lm}$ is the wronskian Morse and Feshbach (1953) of $P_{lm}\left({\xi_{3}}\right)$ and $Q_{lm}\left({\xi_{3}}\right)$, all evaluated at $\xi_{3}=\xi_{0}$. Hence direct resolution of Laplace equation performed in Knipp and Reinecke (1992a, b); Comas et al. (2002); Reese et al. (2004) leads to the full set of the eigenvalues of $\hat{M}$ for the shapes considered in those papers.

The explicit expressions of the eigenvalues and the eigenfunctions are given in mathematical literature for spheres, spheroids and ellipsoids Ahner and Arenstorf (1986); Ahner (1994); Ahner et al. (1994); Ritter (1995). These calculations are based on the expression of the Green’s function (3) in the separable coordinate systems. A similar approach was used for disks and ellipses in two-dimensional (2D) space Ahner (1986). Below we will also provide explicit expressions for the eigenvalues and the eigenfunctions of the electrostatic operator for three-dimensional (3D) counterparts of disks and ellipses, i. e., circular and elliptic cylinders, which to our knowledge is new. In principles, we have to solve the equation

[TABLE]

in an appropriate coordinate systems by separation of variables method. The caveat of our analysis and calculations is the fact that the eigenfunctions of the two commuting operators associated with each separable coordinate system form a basis in which the free-space Green’s function may be expressed and consequently the electrostatic operator can become diagonal. In fact the adjoint operator is diagonal in this basis, while the eigenfunction of the electrostatic operator are only proportional to these functions. This is not a surprise since the adjoint of the electrostatic operator is used in potential theory to generate the solution of Laplace equation with Dirichlet boundary conditions. In the previous subsection we have effectively retrieved the eigenvalues of the electrostatic operator from the solution of Laplace equation for surface phonons and plasmons. In the following we present also a systematic way to obtain not only the eigenvalues but also the eigenfunctions of the electrostatic operator and its adjoint for shapes like spheres, prolate and oblate spheroids, ellipsoids, and elliptic and circular cylinders.

(a) sphere

Working in spherical coordinates, $x=\rho\sin\theta\cos\varphi,\;y=\rho\sin\theta\sin\varphi,\;z=\rho\cos\theta$, the sphere is evidently determined by $\rho=\rho_{0}$, while the eigenfunctions of the electrostatic operator are spherical harmonics $Y_{lm}\left({\theta,\varphi}\right)$ and the eigenvalues are

[TABLE]

Here $l$ is a natural number and $m$ is integer with $\left|m\right|\leq l$. In 3D sphere is the only surface for which the electrostatic operator is symmetric D. Khavison (2007). Expression (18) can be also deduced using Eq. (13) from Ref. Comas et al. (2002).

(b) prolate spheroid

In these coordinate system the transformations are: $x=a\sqrt{({\eta^{2}-1})({1-\zeta^{2}})}\cos\varphi$, $y=a\sqrt{({\eta^{2}-1})({1-\zeta^{2}})}\sin\varphi$, $z=a\eta\zeta$ with $\eta\geq 1$ and $-1\leq\zeta\leq 1$. Prolate spheroids characterized by equation $\eta=\eta_{0}$ have the eigenfunctions of the electrostatic operator proportional to $e^{im\varphi}P_{lm}\left(\zeta\right)$ and its eigenvalues equal to

[TABLE]

where $P_{lm}$ and $Q_{lm}$ are the associated Legendre functions, $l$ is a natural number and $m$ is integer with $\left|m\right|\leq l$. We considered the definitions of $P_{lm}$ and $Q_{lm}$ such that their wronskian is $W\left\{{P_{lm}\left(\eta\right),Q_{lm}\left(\eta\right)}\right\}=(-1)^{m}{\left({l+m}\right)!}\mathord{\left/{\vphantom{{\left({l+m}\right)!}{\left({\left({l-m}\right)!\left({1-\eta^{2}}\right)}\right)}}}\right.\kern-1.2pt}{\left({\left({l-m}\right)!\left({1-\eta^{2}}\right)}\right)}$ Lebedev (1965).

(c) oblate spheroid

The oblate spheroidal coordinates are obtained by replacing $a$ with –ia and $\eta$ with $i\eta$. Hence the eigenfunctions of the electrostatic operator are proportional to $e^{im\varphi}P_{lm}\left(\zeta\right)$ and the eigenvalues are

[TABLE]

The eigenvalues of $\hat{M}$ for prolate and oblate spheroids as well as its eigenfunctions are calculated in Ahner and Arenstorf (1986); Ahner (1994); Ahner et al. (1994). These eigenvalues can also be retrieved with Eq. (13) from Refs. Knipp and Reinecke (1992b); Comas et al. (2002).

(d) ellipsoid

The ellipsoidal coordinates, appropriate for an ellipsoid with axes $a>b>c$, are the solutions $\theta$ of the equation

[TABLE]

which are denoted as $\rho\in\left[{\left.{-c^{2},\infty}\right)}\right.$, $\mu\in\left[{-b^{2},-c^{2}}\right]$, and $\nu\in\left[{-a^{2},-b^{2}}\right]$. The Cartesian coordinates are then expressed in ellipsoidal coordinates by the following

[TABLE]

The ellipsoid is governed by simple equation $\rho=\rho_{0}$. Laplace equation separates in this coordinate system such that for each ellipsoidal coordinate a Lamé equation Arscott and Khabaza (1962) of the form

[TABLE]

is obeyed by all three ellipsoidal coordinates. Here $n$ is a natural number, $B$, which is a separation constant, is real Ritter (1995); Dobner and Ritter (1998), and $\Delta(\rho)=\sqrt{(\rho+a^{2})(\rho+b^{2})(\rho+c^{2})}$. For each $n$ we obtain $2n+1$ values of $B$ denoted as $B_{n}^{m}$ with $m=0,1,\ldots\ldots,2n$. Thus, the solution $L_{n}^{m}$ is the first kind Lamé function of order $n$. This solution is regular, the second solution, which is regular at infinity, is obtained by standard procedure Morse and Feshbach (1953)

[TABLE]

The Wronskian considered here reads $W\left\{{L_{n}^{m}\left(\rho\right),K_{n}^{m}\left(\rho\right)}\right\}=-{\left({2n+1}\right)}\mathord{\left/{\vphantom{{\left({2n+1}\right)}{\Delta\left(\rho\right)}}}\right.\kern-1.2pt}{\Delta\left(\rho\right)}$. All the above lead to the following eigenvalues

[TABLE]

The eigenfunctions of the electrostatic operator are proportional to $L_{n}^{m}\left(\mu\right)L_{n}^{m}\left(\nu\right)$. We have adopted the definition of ellipsoidal coordinates used in Ritter (1995) and in the standard textbooks of electrodynamics like Stratton (1941). It provides the straightforward calculation of the wronskian and of the second kind Lamé function. Slightly different definitions provided in Moon and Spencer (1961) and used in Reese et al. (2004) give the same results Feng and Kang (2016). The eigenvalues (25) were first calculated in Ritter (1995) and can be deduced with the help of (13) from Ref. Reese et al. (2004).

(e) elliptic cylinder

Such cylinders are suitably described in elliptic cylindrical coordinates: $x=d\cosh\left(u\right)\cos\left(v\right)$, $\;y=d\sinh\left(u\right)\sin\left(v\right)$, $\;z=z$, with $d$ an arbitrary positive constant. The surface of the cylinder is given by equation $u=u_{0}$. For convenience we take $d=\sqrt{R^{2}-r^{2}}$, where $R$ and $r$ are the lengths of the longer and, respectively, shorter semi-axis of the cross-sectional ellipsis. Accordingly, the constant $u_{0}$ is equal to ${\ln\left[{{\left({R+r}\right)}\mathord{\left/{\vphantom{{\left({R+r}\right)}{\left({R-r}\right)}}}\right.\kern-1.2pt}{\left({R-r}\right)}}\right]}\mathord{\left/{\vphantom{{\ln\left[{{\left({R+r}\right)}\mathord{\left/{\vphantom{{\left({R+r}\right)}{\left({R-r}\right)}}}\right.\kern-1.2pt}{\left({R-r}\right)}}\right]}2}}\right.\kern-1.2pt}2$. The eigenvalues of the electrostatic operator can be calculated with Eq. (13) from Ref. Knipp and Reinecke (1992a). Also we will see below the separation of variables method provides the eigenvectors of the electrostatic operator which are proportional to either $e^{ikz}ce_{m}\left({q,v}\right)$ or $e^{ikz}se_{m}\left({q,v}\right)$. Here $k$ is real number and $ce_{m}\left({q,v}\right)$, $se_{m}\left({q,v}\right)$ are respectively the “even” and the “odd” solutions of the Matheiu equation Knipp and Reinecke (1992a); McLachlan (1947)

[TABLE]

where $q=-k^{2}\left({R^{2}-r^{2}}\right)$ and $\alpha$ is a separation constant which takes a characteristic value $a_{m}$ or $b_{m}$ of the $m^{th}$ solution $ce_{m}\left({q,v}\right)$ or $se_{m}\left({q,v}\right)$, respectively. If the characteristic values of the Mathieu equation are $a_{m}$ the eigenvalues are

[TABLE]

where $Ce_{m}\left({q,u_{0}}\right)$ and $Fe_{m}\left({q,u_{0}}\right)$ are the solutions of the corresponding associated Mathieu equation

[TABLE]

Here $\alpha$ takes the characteristic value $a_{m}$. Similarly, if $\alpha$ is the characteristic value $b_{m}$ the eigenvalues are

[TABLE]

with $Se_{m}\left({q,u_{0}}\right)$ and $Ge_{m}\left({q,u_{0}}\right)$ the solution of (28) when $\alpha$ 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 Arscott (1964). Thus, the wronskians of the solutions of Mathieu equation as well as of associated Mathieu equation depend only on $q$ (see sections 28.5.8, 28.5.9 and 28.20.3-28.20.7 of Olver et al. (2010)). The eigenvalues expressed by (27) and (29) are consistent with the surface phonon spectrum given in Ref. Knipp and Reinecke (1992a).

The eigenvalues as well as the eigenfunctions of the electrostatic operator for elliptic cylinder can be conveniently calculated by expressing the free-space Green’s function $G\left({\bf{x},\bf{x}^{\prime}}\right)$ in separable elliptic cylindrical coordinates. First, we use the identity $\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{e^{ik\left({z-z^{\prime}}\right)}dk}=\delta\left({z-z^{\prime}}\right)$. Then, it can be shown that the Mathieu functions $\left\{{ce_{m}\left({q,v}\right),se_{m}\left({q,v}\right)}\right\}$ are not only orthogonal on the space of square-integrable function $L^{2}\left({-\pi,\pi}\right)$, but also complete since they are the solutions of a Sturm-Liouville problem Higgins (1977). Hence the completeness relation on $L^{2}\left({-\pi,\pi}\right)$ can take the form $\frac{1}{\pi}\sum\limits_{m=0}^{\infty}{ce_{m}\left({q,v}\right)ce_{m}\left({q,v^{\prime}}\right)+\sum\limits_{m=1}^{\infty}{se_{m}\left({q,v}\right)se_{m}\left({q,v^{\prime}}\right)}}=\delta\left({v-v^{\prime}}\right)$. Based on the separation of variables method we then may show that the free-space Green’s function $G\left({\bf{x},\bf{x}^{\prime}}\right)$ depending on $\left({u,v,z}\right)$ and $\left({u^{\prime},v^{\prime},z^{\prime}}\right)$ has the following expression in elliptic cylindrical coordinates:

[TABLE]

Here $u_{<}(u_{>})$ is the smaller (larger) of $u$ and $u^{\prime}$, $k$ is real and $q=-k^{2}\left({R^{2}-r^{2}}\right)$. The normal derivative to elliptic cylindrical surface is

[TABLE]

Then it is easy to check that the eigenvalues are those expressed by Eqs. (27) and (29) and the eigenfunctions of the $\hat{M}$ are

[TABLE]

while the eigenfunctions of its adjoint have the same expression without the factor $1\mathord{\left/{\vphantom{1{\left({\sinh\left(u\right)^{2}+\sin\left(v\right)^{2}}\right)^{1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2}}}}\right.\kern-1.2pt}{\left({\sinh\left(u\right)^{2}+\sin\left(v\right)^{2}}\right)^{1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2}}$. These results are the 3D generalization of 2D ellipses D. Khavison (2007); Ahner (1986); Ando and Kang (2016) which can be recovered if we take $k=q=0$. So, if we consider $k=q=0$, then $a_{m}=b_{m}=m^{2}$, $ce_{m}\left({q,v}\right)=\cos\left({mv}\right)$, $se_{m}\left({q,v}\right)=\sin\left({mv}\right)$, and the eigenvalues are

[TABLE]

During the preparation of the manuscript a recent paper has come to our attention Cohl and Volkmer (2012). In that paper the eigenfunction expansions of the free-space Green’s function are calculated in elliptic and parabolic cylinder coordinates Cohl and Volkmer (2012). The authors used a different proof based on Lipschitz and Lipschitz-Hankel integral identities but their expansion given in theorem 4.3 is similar to our expression (III.2). Also, in theorem 5.3 the authors used other two linearly independent solutions of the associated Mathieu equations whose wronskian is -1.

(f) circular cylinder

In cylindrical coordinates $x=\rho\sin\theta,\;y=\rho\cos\theta,\;z=z$ the cylinder with circular cross-section is described by $\rho=\rho_{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 coordinate that is given, for example, in Jackson’s textbook Jackson (1975). In cylindrical coordinates $\left({\rho,\theta,z}\right)$ and $\left({\rho^{\prime},\theta^{\prime},z^{\prime}}\right)$ the free-space Green’s function has the following form

[TABLE]

where $I_{m},K_{m}$ are the modified Bessel functions whose wronskian is $W\left\{{I_{m}\left(\rho\right),K_{m}\left(\rho\right)}\right\}=-1\mathord{\left/{\vphantom{1\rho}}\right.\kern-1.2pt}\rho$ and $\rho_{<}(\rho_{>})$ is the smaller (larger) of $\rho$ and $\rho$’. It is easy to check that the eigenfunctions of $\hat{M}$ and of its adjoint are the same, namely $e^{ikz}e^{im\theta}$, while the eigenvalues are equal to

[TABLE]

If we consider $k=0$ then

[TABLE]

and

[TABLE]

for $m\neq 0$ which are the eigenvalues of a disk in 2D D. Khavison (2007). $\hat{M}$ and its adjoint having the same eigenfunctions is a reminiscense of the fact that the electrostatic operator is symmetric for disks in 2D Ahner (1986).

Both types of cylinders have a continuum spectrum, which is a signature of accommodating running wave along the z-direction, since cylindrical surfaces are unbounded hence, non-compact. Similar to phonon polaritons and implicitly using the separation of variables method the full spectrum of surface plasmons was calculated for spherical Prodan and Nordlander (2004), for spheroidal Brandl and Nordlander (2007), and for circular cylindrical shapes Morandi (2008); Wan et al. (2013) in the so called hybridization model of plasmon resonances. Analysis of plasmon resonances in Prodan and Nordlander (2004); Brandl and Nordlander (2007); Morandi (2008); Wan et al. (2013) leads to the same conclusion regarding the eigenvalues of the electrostatic operator.

All these shapes presented above generate also 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 Refs. Englman and Ruppin (1966, 1968); Prodan and Nordlander (2004); Brandl and Nordlander (2007); Morandi (2008); Wan et al. (2013).

IV Discussions

The BIE method can be a 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 Ref. Fonoberov and Balandin (2004). The BIE method can be also a convenient tool to comprehend and calibrate results with mixed states made by coupling continuum and localized states like those studied in a recent work Gubbin et al. (2016). Moreover, a boundary integral equation method provides more insight about the localized resonance modes. For instance in Ref. Feng et al. (2015) there are studied the LSPhRs in an array of micro-disks of gallium nitride on silicon carbide substrate. The effect of localized surface phonon polariton resonances were 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-disk 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 disk, while in the TM configuration both types of modes (i. e., parallel and perpendicular to the disk surface) are excited. Although our calculations are qualitative since they are made for a single disk some insight can be obtained since we have obtained all possible resonances of such micro-disk. In Fig. 2 we may see that for the electric parallel to the disk surface there are only two resonances of which one is predominant. On the contrary, for the electric field perpendicular to the disk there are several resonances of which two are predominant. With the BIE method we may see that the small discrepancy noticed in Ref. Feng et al. (2015) 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 has sparked new interest arising from mathematical theory regarding plasmons, cloaking as well as various inverse problems Nicorovici et al. (1994); Ando and Kang (2016); Ammari et al. (2014, 2016, 2013); Milton et al. (2005); Milton and Nicorovici (2006). For bounded and smooth surfaces the electrostatic operator is compact with a countable number of eigenvalues that accumulates at the origin which belongs to the essential spectrum D. Khavison (2007). The eigenvalues lie within (-1/2,1/2] with the eigenvalue 1/2 corresponding to the equilibrium charge distribution D. Khavison (2007); Sandu et al. (2013). For shapes resulted from separable coordinate systems, the capacitance of equilibrium charge can be calculated by simple integration Sandu et al. (2013).

Questions regarding spectral properties like 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 Nicorovici et al. (1994); Ando and Kang (2016); Ammari et al. (2013); Milton et al. (2005); Milton and Nicorovici (2006). Several 2D systems have been considered for cloaking including concentric disks Nicorovici et al. (1994); Milton et al. (2005); Milton and Nicorovici (2006), ellipses Ando and Kang (2016), and confocal ellipses Chung et al. (2014). All these approaches need also the eigenfunctions beside all the eigenvalues. The results of previous section regarding the eigenvalues and the eigenfunctions of the electrostatic operator defined on different surfaces and on their shelled structures counterparts may be used to extend the analysis of cloaking from disks and ellipses to 3D setting.

Since 0 belongs to the essential spectrum a legitimate question is if there are shapes which have 0 as an eigenvalue and more general if there is a shape that may have as an eigenvalue any number within -1/2 and 1/2. The general answer is given in D. Khavison (2007) and in Ahner et al. (1994) for spheroids and more recently in Feng and Kang (2016) for general ellipsoids. In Fig. 3 we present a graphical proof where we plotted the eigenvalues $\chi_{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 $z$-direction and $a_{x}$ is the semi-axis along $x$-direction. An aspect ratio $>1$ defines a prolate spheroid and an aspect ratio $<1$ defines an oblate spheroid. Prolate spheroids reach the value 1/2 very fast at a rate inverse proportional to the square of the aspect ratio (i. e., $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2-\chi_{10}}\right)\sim\left({{a_{x}}\mathord{\left/{\vphantom{{a_{x}}{a_{z}}}}\right.\kern-1.2pt}{a_{z}}}\right)^{2})$ thus, practically, for aspect ratios $>10$, $\chi_{10}$ is pretty close to 1/2. From Eqs. (7) and (10) we can see that this eigenvalue generates a resonance that is redshifted to IR (theoretically it goes to 0) for plasmons or is moved toward $\omega_{T}$ for phonon polaritons.

It can be seen in Fig. 3 that there is an oblate spheroid with an aspect ratio between 0.1 and 1 for which $\chi_{10}=0$. Furthermore, as the aspect ratio tends to 0 $\chi_{10}$ approaches -1/2 as $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2+\chi_{10}}\right)\sim\left({{a_{z}}\mathord{\left/{\vphantom{{a_{z}}{a_{x}}}}\right.\kern-1.2pt}{a_{x}}}\right)^{2}$ and $\chi_{11}$ goes to 1/2 as $\left({1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2-\chi_{11}}\right)\sim\left({{a_{z}}\mathord{\left/{\vphantom{{a_{z}}{a_{x}}}}\right.\kern-1.2pt}{a_{x}}}\right)^{2}$. Similarly, from Eqs. (7) and (10) we see that if the electrostatic operator has an eigenvalues close to -1/2 the plasmon resonance frequency approaches the plasma frequency of the bulk material and the phonon polariton resonances shifts toward $\omega_{L}$.

In other words, asymptotically for $l$=1 a 2D disk embedded in 3D has an eigenvalue -1/2 for fields perpendicular to the surface of the disk and an eigenvalue equals to 1/2 for fields parallel to the surface. The rates of approaching these values are proportional to the square of the aspect ratio. On the other hand, in contrast to oblate disks, for finite and constant-thickness disks, like those discussed at the beginning of this section, there is an eigenvalue that approaches 1/2 at a rate proportional to the aspect ratio, which is defined as the ratio of the thickness and disk diameter. This assertion is shown in the following. This mode defines the graphene plasmon resonance which is excited by an electric field parallel to the graphene plane Fang et al. (2013). For instance, considering graphene as a 2D sheet, the graphene plasmon frequency of nano-disks scales as $1\mathord{\left/{\vphantom{1{\sqrt{d}}}}\right.\kern-1.2pt}{\sqrt{d}}$, where $d$ is the diameter of of graphene nano-disk Fang et al. (2013). In 3D models, graphene nano-disks are simulated as finite-thickness disks with a bulk complex conductivity obtained from graphene sheet-conductivity divided by the disk thickness $h$. Hence the finite thickness nano-disk has the following dielectric relative permittivity Vakil and Engheta (2011)

[TABLE]

Here $\sigma\left(\omega\right)={ie^{2}E_{F}}\mathord{\left/{\vphantom{{ie^{2}E_{F}}{\left[{\left({\pi\hbar^{2}}\right)\left({\omega+i\mathord{\left/{\vphantom{i\tau}}\right.\kern-1.2pt}\tau}\right)}\right]}}}\right.\kern-1.2pt}{\left[{\left({\pi\hbar^{2}}\right)\left({\omega+i\mathord{\left/{\vphantom{i\tau}}\right.\kern-1.2pt}\tau}\right)}\right]}$ is the graphene sheet conductivity that has a Drude-like form and $\varepsilon_{v}$ is the vacuum permittivity. We denoted by $E_{F}$ the Fermi energy in graphene disk, $\tau$ is the relaxation time, and $\hbar$ is the reduced Planck constant. Thus, the plasma frequency for a 3D model of graphene takes the form $\omega_{p}^{2}={e^{2}E_{F}}\mathord{\left/{\vphantom{{e^{2}E_{F}}{\left({\pi\hbar^{2}\varepsilon_{v}}\right)}}}\right.\kern-1.2pt}{\left({\pi\hbar^{2}\varepsilon_{v}h}\right)}$. Bearing in mind that the eigenvalues are scale invariant they are also a function of the aspect ratio $h$/$d$ only. Since the plasmon frequency is $\tilde{\omega}_{pk}=\omega_{p}\sqrt{1/2-\chi_{k}}$ (see Eq. (7)) we may deduce that the graphene plasmon frequency scales as $1\mathord{\left/{\vphantom{1{\sqrt{d}}}}\right.\kern-1.2pt}{\sqrt{d}}$ and that its corrsponding eigenvalue goes to $1/2$ at the rate linear with $h$/$d$ and not quadratically like in the case of oblate disks.

Owing to the success of the BIE in the quasi-static regime there is a renewed interest in the modal approach based on BIE to plasmon problems in the full-wave electromagnetic regime Makitalo et al. (2014). The full-wave electromagnetic treatment discussed in Ref. Makitalo et al. (2014) has at its core the free-space Green’s function of scalar Helmholtz equation that is also separable in those eleven coordinate system in which Laplace eqation separates Boyer et al. (1976); Miller (1977). In general, in the full-wave regime, due to radiation losses, the eigenvalues are no longer real as in the case of quasistatic regime Makitalo et al. (2014). It will be interesting to apply the same symmetry arguments about 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 Hsu et al. (2016) with no radiation losses in structures having other forms than spherical shapes which have been considered recently Silveirinha (2014); Monticone and A.Alu (2014).

V Conclusions

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

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 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 2D to 3D 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 of studying and obtaining bound states in continuum with shapes different from spherical geometry.

Acknowledgements.

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