An eigenvalue approach to quantum plasmonics based on a self-consistent hydrodynamics method
Plasmonics has attracted much attention not only because it has useful properties such as strong field enhancement, but also because it reveals the quantum nature of matter. To handle quantum plasmonics effects, ab initio packages or empirical Feibelman d-parameters have been used to explore the quantum correction of plasmonic resonances. However, most of these methods are formulated within the quasi-static framework. The self-consistent hydrodynamics model offers a reliable approach to study quantum plasmonics because it can incorporate the quantum effect of the electron gas into classical electrodynamics in a consistent manner. Instead of the standard scattering method, we formulate the self-consistent hydrodynamics method as an eigenvalue problem to study quantum plasmonics with electrons and photons treated on the same footing. We find that the eigenvalue approach must involve a global operator, which originates from the energy functional of the electron gas. This manifests the intrinsic nonlocality of the response of quantum plasmonic resonances. Our model gives the analytical forms of quantum corrections to plasmonic modes, incorporating quantum electron spill-out effects and electrodynamical retardation. We apply our method to study the quantum surface plasmon polariton for a single flat interface.
Keywords: Eigenvalue approach, Quantum plasmonics, Self-consistent hydrodynamics method, Surface plasmon-polariton dispersion
Plasmonic resonances are the intrinsic modes of metallic systems, and these collective excitations have found widespread applications in sensing [1, 2, 3, 4, 5], photonics [6, 7, 8, 9, 10, 11, 12], the lightning rod effect  and the realization of the perfect lens [11, 12]. Recently, with the development of sophisticated and controllable fabrication techniques for metallic nanoparticles, the quantum nature of electrons has become more obvious and the classical description of local electrodynamics has become inadequate [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. Quantum effects, such as the nonlocal response [13, 14, 15, 16, 17, 18, 19, 20] and electron spill-out [21, 22, 23, 24], must be taken into account. The most popular ab initio numerical approach to studying the plasmonic quantum effect is time-dependent density functional theory (TD-DFT) [25, 26, 27, 28, 29, 30]. In addition, the hydrodynamics method has also been employed to study quantum plasmonics [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. While TD-DFT is very useful within the quasi-static framework, the method becomes very computationally demanding if it is coupled with the full Maxwell equations to include electromagnetic (EM) wave characteristics such as retardation [42, 43]. It is also difficult to extract the physics from fully ab initio numerical simulations, as the sheer volume of information produced can be overwhelming. In order to understand the underlying physics, semi-empirical approaches such as the Feibelman d-parameters  have been introduced to model the electron spill-out effect for a single flat interface, and the predicted quantum corrections to surface plasmons have been verified experimentally [45, 46, 47, 48, 49, 50]. The d-parameter framework has subsequently been extended to more complicated shapes through the boundary element method (BEM) . However, the method of Feibelman d-parameters is still formulated within the quasi-static approximation, meaning that no retardation is incorporated. From an optics and photonics point of view, the Hamiltonian (eigenvalue) approach is used less frequently than scattering formalisms, for the reason that many optics problems are formulated with open boundary conditions. But the eigenvalue approach has proved indispensable in situations where periodic boundary conditions are imposed [52, 53], and it is also useful in describing novel properties that are defined using eigenfunctions (such as topological invariants) . Up to now, essentially all classical Hamiltonian approaches assume that the system must have well-defined macroscopic permittivity and permeability . How to include the quantum effects of electron gases in the Hamiltonian is still an open question.
Plasmon-polariton excitations, as bound states outside the light cone in a phase space, can be fully described using an eigenvalue formulation. In this article, we attempt to establish the eigenvalue approach based on the self-consistent hydrodynamics model (SC-HDM) to investigate the quantum correction of plasmonic modes. We find that a global operator must exist in the eigenvalue approach in order to make it Hermitian (in the limit of no dissipation). The existence of this global operator means that the response of the system is nonlocal, which is consistent with our physical understanding that classical local electrodynamics is an approximation that will break down when quantum effects kick in. By employing first-order perturbation theory, we give the analytical form of quantum corrections to the plasmonic modes. As an example, we apply our method to study the surface plasmon polaritons (SPPs) of a single flat surface and find that our method agrees with that of the Feibelman d-parameters in the intermediate region of wavenumber space. The paper is organized as follows. In section 2, we formulate an eigenvalue approach based on SC-HDM. The related issues on the global operator and the inner product are discussed in section 2.2 and section 2.3, respectively. The quantum corrections to a general plasmonic mode and the dispersion of quantum SPPs are given in section 3. Conclusions are drawn in section 4.
To formulate the eigenvalue approach, we employ SC-HDM because it treats electrons and photons on the same footing [55, 56, 57, 58, 59, 60]. However, SC-HDM uses the electron density and not the electronic wavefunction as the basic variable, it is not as accurate as local density functionals. Yet previous studies have shown that SC-HDM gives reasonable results compared with ab initio calculations [55, 56, 57, 58, 59, 60] but at a much lower computational cost. In the following, we will show the procedure for establishing the eigenvalue approach based on SC-HDM.
2.1 Eigenvalue approach
The first step of SC-HDM involves determining the electronic ground state density that minimizes a density functional subject to constraints (chemical potential and electron number). This step requires the numerical calculation of the equilibrium electron density and the effective single electron potential . Once the ground state density is obtained, the excited state calculations can be performed numerically by coupling the linearized equations of motion for the electron gas with Maxwell equations as 
where is the angular frequency, is the speed of light in vacuum, () is the mass (charge) of an electron, is the microscopic electric field, is the induced current, is the induced charge density, and is the loss parameter. The first term in (3) gives the nonlocal and spill-out effect of the electron gas, an important effect at the nano-scale that is absent in the classical electromagnetic approaches. In principle, the density functional has many contributions, including kinetic energy, exchange-correlation, and Coulomb interactions. The exchange-correlation part is explicitly quantum in nature. All other terms depend on the density, which itself is determined variationally by minimizing with respect to certain constraints. We can view those physical effects that originate from as quantum in nature. Nevertheless, some of the terms can be derived through a semi-classical approach. We will discuss the explicit forms of later. Using (2), we see that the currents can be used to replace as the variable in (3), so we introduce an operator defined as
We assume that the dispersion of electrons near the Fermi energy is parabolic, namely
where , , and are the Fermi velocity, the Fermi wave vector and the ion density, respectively. The ion density is defined via a dimensionless quantity as
where is the Bohr radius. Furthermore, the corresponding Thomas-Fermi screening wave vector and plasma frequency are both defined through the ion density as 
where is a dimensionless quantity. Note that we have added the term which stands for interband resonances .
In the absence of losses (), the matrix on the left-hand side of (13) would be a Hermitian matrix, so that
Note that in (14) we have assumed the existence of the inverse of . If we define the operator as
then the matrix could be written as
With the help of , equation (13) gives the operator for a photonic system
If we set , then , as it must be if there is no dissipation. Meanwhile the auxiliary quantity is evaluated as
It is easy to see that the above procedure automatically introduces magnetic fields into the formulations. Let us further define and . Then is simplified as
We can now introduce the four-component vector: . The eigenvalue formulation is then
We have now completely derived the eigenvalue approach for a photonic system expressed in terms of microscopic electron densities rather than macroscopic permittivity. In the photonic crystal literature, in equation (20) is often called the Hamiltonian because it can be transformed to the wavevector space and hence forms a Hilbert space. In the following, we discuss the operator and the related definition of the inner product within the current framework.
2.2 Energy functional and nonlocal operator
In the previous section, we showed that the eigenvalue approach for photonic systems based on microscopic electron densities could be formulated using the operator , which comes from the internal energy of the electron gas. The explicit form of should be specified in order to obtain an expression for . Throughout this paper, we use the form of used commonly in the literature [56, 60]. can then be obtained as
where , and the coefficients are
comes from the kinetic energy of the inhomogeneous electron gas, and describes effective masses. It is worth mentioning that there are different choices of the exchange-correlation energy functional in , but they all give essentially the same result .
In (15), the first term represents the interband resonance and the second term the quantum correction. The exact form of under the energy functional operator is not easy to obtain because it involves fractional derivative operators . The leading differential order in (21) is the first-order derivative, and hence there must exist terms in . Mathematically, the fractional derivative of to order is often defined by means of Fourier or Mellin integral transforms . at a point is a local property only when is an integer, but if is not an integer, then at not only depends on the values of near , but also relates to in the whole domain. This indicates that must be a global operator, illustrating that the response of the electron gas is nonlocal, which is consistent with our physical understanding as shown in (3).
However, if we assume that the second term in (15) is much smaller than the first term, then the operator could be approximately written as
Using (27), the fourth component of can be explicitly written as
2.3 Inner product
For a completeness method, we define the inner product as
The bracket in the integral is the energy density of a particular photonic state (for ) with the internal energy of the electron gas taken into account. Let denote the bracket. Then the orthogonal condition (29) used within this method is
Note once again that our method is an eigenvalue approach, while others are based on the scattering approach . Computationally, the incident EM wave must be specified for the scattering approach, and the focus is on the scattering of the incident wave by the plasmonic object. Our method solves an eigenvalue problem, similar to band dispersion calculations, and the output consists of the eigenmodes (e.g. surface plasmons) supported by the system. Thus the eigenvalue approach shown in this section is quite different from the scattering approach although the partial differential equations are the same.
3 Results and discussions
We have established the eigenvalue approach based on SD-HDM. In this section, we give an explicit expression for quantum corrections to certain plasmonic modes using first-order perturbation theory, and discuss the dispersion of quantum SPPs for a single flat interface.
3.1 Quantum correction of plasmonic modes
To study the quantum correction of particular plasmonic modes, we can treat the terms involving quantum corrections and loss in (17) as the perturbation potential, namely in which
where , and describes the ion density in the jellium model . Mathematically, is the indicator function and stands for the union of all metallic domains. Suppose the eigenmodes of are ( is the plasmonic mode index). Then perturbation theory gives the first-order corrections to the plasmon resonance frequency as
where denotes the classical energy density of the photonic state, i.e. setting in (30). For simplicity, we have omitted the superscript in the fields and will use this notation in the following. Before proceeding, let us discuss the conditions under which the perturbation approximation holds. The condition indicates that the variations in ground state densities and the surface-to-volume ratio of the nanoparticles should be small.
In (35) there are three terms: the first one is due to loss, and the other two terms are quantum corrections. The first term is a purely imaginary number, which could be evaluated as
This is consistent with classical EM results in . In order to relate the other two terms to the Feibelman -parameters ( and ) defined in quasi-static approximations which are used to explain quantum corrections to plasmonic modes , we introduce the quantum correction factor as follows:
Within the quasi-static limit, BEM gives the quantum correction factor as and
where and are the surface charge density and the electrostatic potential in the classical model . Within the quasi-static framework, the quantum correction factor can be factorized into the shape factor and the electron spill-out factor . However, while the d-parameters are well defined for a single flat interface, they are not applicable to the description of nanoparticles with complex geometries, especially those with sharp corners.
Next, we give the explicit form of the quantum correction factor within our method. Using the properties in the classical model
we could explicitly obtain
The first term in the integral depends on the electrostatic surface dipole and electric field intensities, while the second term originates from the nonlocal operators, indicating that the second term depends more on details of the mode than the first one. To obtain physical insights into this correction, we keep the leading derivative term in , namely
Then (42) becomes
where we set because of our focus on plasmonic systems. Firstly, according to our method, the contributions from the electron spill-out effect and shape effects to the quantum corrections are mixed. Moreover, the first term in the bracket does not have any feature length explicitly, but the second term carries the Thomas-Fermi screening length scale . For example, the typical value of is for sodium. Secondly, the major contribution to the integration comes from the surface region of the metallic particles, because is only nonzero near the surface region, is zero in the bulky vacuum domain, and is zero in the bulky metal domain.
To seek a clearer understanding, we use the following relations:
where , , and and are bulk and surface charge densities, and / stands for tangential/normal directions of the surface . Then we could simplify the quantum correction factor in (44) as
where the feature length is determined from classical results and material properties.
There are four terms in (47). The first two are volume integrals in all of the space, but the contributions mainly come from the surface region of the metallic particles. Meanwhile, the last two terms are integrals on the metal surface. Therefore the quantum correction factor is nonnegligible only for particles with a large surface-to-volume ratio.
3.2 Dispersion of quantum surface plasmon polaritons
In this section, we aim to apply our method to a single plasmonic interface as shown in figure 1(a). To see the coefficients in (47), we plot , and of this single interface case in figure 1(b). It is clear that all of the coefficients are nonzero near the surface as expected.
The next task is to explicitly write down the classical results of SPPs for the system in figure 1(a). The magnetic fields are given as
where , is the parallel wavenumber, and . The related electric fields are
The boundary condition of gives the dispersion relation of SPPs as
where the dimensionless quantities and are introduced to simplify the expressions.
with the following notations:
The second term in (47) is
The third term of the quantum correction factor in (47) for the single interface case is zero due to the fact that is perpendicular to the surface. However, for particles have corners, must have a parallel component to the surface, so this part of the quantum correction factor can become prominent. The fourth term of the quantum correction factor in (47) for the single interface case is
The full numerical results of the SPP dispersion with quantum correction factor are shown by the solid red line in figure 2(a). For comparison, we also plot the SPP dispersion of the classical model and the d-parameter method using green and blue lines, respectively [44, 51, 60, 64]. A magnified version of the SPP dispersion in the small wavenumber limit is shown in figure 2(b). We see that the d-parameter approach works well in the intermediate region of space.
Before ending this section, let us discuss the asymptotic behaviors of . Firstly, in the limit of a small , namely , because , as shown by the red and green lines in figure 2(b). Secondly, in the large limit, namely ,
It is not difficult to see that , and if we are only concerned with the behavior near the plasmon resonance (), then the quantum correction factor could be further simplified as
where and . These two integrals are infinite series in , and the leading order in is a constant, indicating that in (67) is linear in . This is numerically shown by the red and blue lines in figure 2(b).
We formulated an eigenvalue approach for plasmonic resonances using the self-consistent hydrodynamics model. We showed that the Hamiltonian carries a global operator, indicating that the response of quantum plasmonic resonances is highly nonlocal at the nano-scale. We derived the analytical forms of quantum corrections to a general plasmonic mode. The calculated dispersions of quantum surface plasmon polaritons for a single interface show that in the intermediate region, our results agree well with the Feibelman d-parameter method.
-  Anker J N, Hall W P, Lyandres O, Shah N C, Zhao J and Van Duyne R P 2008 Nature Materials 7 442–453
-  Li M, Cushing S K and Wu N 2014 Analyst 140 386–406
-  Wong C L and Olivo M 2014 Plasmonics 9 809–824
-  Liang F, Guo Y, Hou S and Quan Q 2017 Science Advances 3 e1602991
-  Zhang S, Bao K, Halas N J, Xu H and Nordlander P 2011 Nano Letters 11 1657–1663
-  Nicoletti O, de la Pena F, Leary R K, Holland D J, Ducati C and Midgley P A 2013 Nature 502 80–84
-  Schuller J A, Barnard E S, Cai W, Jun Y C, White J S and Brongersma M L 2010 Nature Materials 9 193–204
-  Ebbesen T W, Lezec H J, Ghaemi H F, Thio T and Wolff P A 1998 Nature 391 667–669
-  Lezec H J, Dionne J A and Atwater H A 2007 Science 316 430
-  Liao P F and Wokaun A 1982 The Journal of Chemical Physics 76 751–752
-  Pendry J B 2000 Physical Review Letters 85 3966
-  Fang N, Lee H, Sun C and Zhang X 2005 Science 308 534
-  Ciraci C, Hill R T, Mock J J, Urzhumov Y, Fernandez-Dominguez A I, Maier S A, Pendry J B, Chilkoti A and Smith D R 2012 Science 337 1072–1074
-  Scholl J A, Koh A L and Dionne J A 2012 Nature 483 421–427
-  Teperik T V, Nordlander P, Aizpurua J and Borisov A G 2013 Physical Review Letters 110 263901
-  Stella L, Zhang P, Garcia-Vidal F J, Rubio A and Garcia-Gonzalez P 2013 The Journal of Physical Chemistry C 117 8941–8949
-  Wiener A, Fernandez-Dominguez A I, Horsfield A P, Pendry J B and Maier S A 2012 Nano Letters 12 3308–3314
-  Raza S, Toscano G, Jauho A P, Wubs M and Mortensen N A 2011 Physical Review B 84 121412
-  Fernandez-Dominguez A I, Wiener A, Garcia-Vidal F J, Maier S A and Pendry J B 2012 Physical Review Letters 108 106802
-  Luo Y, Fernandez-Dominguez A I, Wiener A, Maier S A and Pendry J B 2013 Physical Review Letters 111 093901
-  Marinica D, Kazansky A, Nordlander P, Aizpurua J and Borisov A G 2012 Nano Letters 12 1333–1339
-  Jin D, Hu Q, Neuhauser D, von Cube F, Yang Y, Sachan R, Luk T S, Bell D C and Fang N X 2015 Physical Review Letters 115 193901
-  Zhu W, Esteban R, Borisov A G, Baumberg J J, Nordlander P, Lezec H J, Aizpurua J and Crozier K B 2016 Nature Communications 7 11495
-  Esteban R, Borisov A G, Nordlander P and Aizpurua J 2012 Nature Communications 3 825
-  Yan W, Wubs M and Asger Mortensen N 2015 Phys. Rev. Lett. 115(13) 137403
-  Zhang P, Feist J, Rubio A, Garcia-Gonzalez P and Garcia-Vidal F J 2014 Phys. Rev. B 90(16) 161407
-  Li J H, Hayashi M and Guo G Y 2013 Phys. Rev. B 88(15) 155437
-  Yan J, Yuan Z and Gao S 2007 Phys. Rev. Lett. 98(21) 216602
-  Marinica D C, Zapata M, Nordlander P, Kazansky A K, Echenique P M, Aizpurua J and Borisov A G 2015 Science Advances 1 e1501095
-  Wang B J, Xu Y and Ke S H 2012 The Journal of Chemical Physics 137 054101
-  CiracÃ¬ C, Pendry J B and Smith D R 2013 ChemPhysChem 14 1109–1116
-  Li X, Fang H, Weng X, Zhang L, Dou X, Yang A and Yuan X 2015 Optics Express 23 29738
-  Esteban R, Zugarramurdi A, Zhang P, Nordlander P, Garcia-Vidal F J, Borisov A G and Aizpurua J 2015 Faraday Discussions 178 151–183
-  Ford G W and Weber W H 1984 Physics Reports 113 195–287
-  Intravaia F and Busch K 2015 Phys. Rev. A 91(5) 053836
-  Toscano G, Raza S, Jauho A P, Mortensen N A and Wubs M 2012 Optics Express 20 4176
-  Toscano G, Raza S, Yan W, Jeppesen C, Xiao S, Wubs M, Jauho A P, Bozhevolnyi S I and Mortensen N A 2013 Nanophotonics 2 161–166
-  Sobhani A, Manjavacas A, Cao Y, McClain M J, Garcia de Abajo F J, Nordlander P and Halas N J 2015 Nano Letters 15 6946–6951
-  Varas A, GarcÃa-GonzÃ¡lez P, Feist J, GarcÃa-Vidal F and Rubio A 2016 Nanophotonics 5 409–426
-  Trugler A, Hohenester U and Garcia de Abajo F J 2017 International Journal of Modern Physics B 31 1740007
-  Moradi A 2015 Physics of Plasmas 22 032112
-  Yabana K, Sugiyama T, Shinohara Y, Otobe T and Bertsch G F 2012 Phys. Rev. B 85(4) 045134
-  Sato S A, Yabana K, Shinohara Y, Otobe T and Bertsch G F 2014 Phys. Rev. B 89(6) 064304
-  Feibelman P J 1982 Progress in Surface Science 12 287–407
-  Liebsch A 1997 Electronic excitations at metal surfaces (Plenum Press, New York)
-  Tsuei K D, Plummer E W and Feibelman P J 1989 Phys. Rev. Lett. 63(20) 2256–2259
-  Tsuei K D, Plummer E W, Liebsch A, Kempa K and Bakshi P 1990 Phys. Rev. Lett. 64(1) 44–47
-  Tsuei K D, Plummer E W, Liebsch A, Pehlke E, Kempa K and Bakshi P 1991 Surface Science 247 302–326
-  Sprunger P T, Watson G M and Plummer E W 1992 Surface Science 269-270 551–555
-  Zacharias P and Kliewer K L 1976 Solid State Communications 18 23–26
-  Christensen T, Yan W, Jauho A P, Soljacic M and Mortensen N A 2017 Physical Review Letters 118 157402
-  Raman A and Fan S 2010 Phys. Rev. Lett. 104(8) 087401
-  Sakoda K 2005 Optical Properties of Photonic Crystals (Springer-Verlag Berlin Heidelberg)
-  Lu L, Joannopoulos J D and Soljacic M 2014 Nature Photonics 8 821
-  Zaremba E and Tso H C 1994 Phys. Rev. B 49(12) 8147–8162
-  Yan W 2015 Phys. Rev. B 91(11) 115416
-  Toscano G, Straubel J, Kwiatkowski A, Rockstuhl C, Evers F, Xu H, Asger Mortensen N and Wubs M 2015 Nature Communications 6 7132
-  Ciraci C and Della Sala F 2016 Phys. Rev. B 93(20) 205405
-  Ciraci C 2017 Phys. Rev. B 95(24) 245434
-  Ding K and Chan C T 2017 Phys. Rev. B 96(12) 125134
-  Kittel C 2005 Introduction to Solid State Physics (John Wiley and Sons, Inc, New York)
-  Herrmann R 2014 Fractional Calculus — An Introduction for Physicists (World Scientific, Singapore)
-  Borisov A G and Shabanov S V 2006 Journal of Computational Physics 216 391–402
-  Mahan G D 1981 Many-Particle Physics (Plenum Press, New York)