Complex band structure and plasmon lattice Green's function of a periodic metal-nanoparticle chain
11 Complex band structure and plasmon lattice Green’s function of a periodic metal-nanoparticle chain
Kin Hung Fung, Ross Chin Hang Tang, and C. T. Chan Department of Physics, The Hong Kong University of Science and Technology, Hong Kong, China When the surface plasmon resonance in a metal-nanoparticle chain is excited at one point, the response signal will generally decay down the chain due to absorption and radiation losses. The decay length is a key parameter in such plasmonic systems. By studying the plasmon lattice Green’s function, we found that the decay length is generally governed by two exponential decay constants with phase factors corresponding to guided Bloch modes and one power-law decay with a phase factor corresponding to that of free space photons. The results show a high level of similarity between the absorptive and radiative decay channels. By analyzing the poles (and the corresponding residues) of the Green’s function in a transformed complex reciprocal space, the dominant decay channel of the real-space Green’s function is understood. PACS number(s): 78.67.Bf, 73.22.Lp, 73.20.Mf, 78.70.-g
I. INTRODUCTION
Band theoretic techniques are widely applied in describing the eigenmodes of periodic systems. With the help of the Bloch’s theorem, a band structure for propagating modes in a lossless periodic system is usually displayed as a plot of the eigenmode as a function of real frequency ( ω ) versus real Bloch wave vector ( k ) in the reduced-zone scheme. Such a “real- k ” band-structure description has been generalized for including evanescent waves within bandgap by allowing the Bloch wave vectors to be complex quantities. This generalized “complex- k ” description is sometimes called the complex-band-structure description. There are several advantages of considering the complex band structure. For example, the magnitude of the imaginary part of the k -vector gives the strength and the physical nature of the bandgaps. In addition, it can describe the decaying waves in lossy systems. It is generally believed that a non-zero imaginary part of k implies an exponential decay (or growth) of the wave with a decay constant given by k = (cid:65) . Although the actual decay profile may depend on the incident field of a particular problem, such kind of exponential decay characteristics are indeed embedded in the Green’s function. For example, the frequency-domain Green’s function for a homogeneous one-dimensional (1D) wave equation takes the form, Re( )| | Im( )| | ( , ) ~ i k x k x
G x e e ω − , where ω is the excitation frequency, k is the corresponding complex wave number, and x is the spatial coordinate. Since the Green’s function contains all the essential information on wave propagation, we can have a useful interpretation of the complex wave properties by studying the decay profile of the Green’s function and its relation to the complex band structure. Although there were some studies in the literature that is concerned with the Green’s function for inhomogeneous lossy media, only a few of them focus on surface plasmons, which have been shown to exhibit many interesting properties in recent years. In plasmonic systems, which are intrinsically lossy due to absorption and radiation, the primary concern is the decay length of the plasmonic modes that can be excited in a certain point. In this paper, we study the Green’s function for the plasmonic waves that propagate along a linear array of metal nanoparticles (MNPs). We will first discuss the numerical results on the complex band structure (Sec. II) and the lattice Green’s function (Sec. III) for such MNP array. Then, we will propose an ansatz which gives a closed-form expression for the lattice Green’s function and we will show that the ansatz is indeed the correct solution by considering a generating function for the lattice Green’s function (Sec. IV).
II. COMPLEX BAND STRUCTURE
We consider a linear array of identical MNPs (with lattice constant a and particle radius r ). For a r ≥ and frequencies close to the Fröhlich frequency, the electromagnetic modes of the MNPs can be accurately described by the coupled-dipole equations, ⎥⎦⎤⎢⎣⎡ −+= ∑ ≠ nmn nmextmm pRRWEp )( (cid:73) α , (1) where m p is the electric dipole moment of the m th particle, α is the single-particle dynamic polarizability, and extm E is the external electric field acting on the m th particle. The dynamic propagator, )( rW (cid:73) , for a dipole in a homogeneous host medium can be written as )( r uv W ]/)()([ rrrrkBrkAk vuhuvhh += δ , where hh ck / ω= , h c is the speed of light in the host medium, )( xA ( ) ix exixx −−− −+= , )( xB ( ) ix exixx −−− +−−= , and = vu are component indices in Cartesian coordinates. Equation (1) can be written as extmuvn nvmunv EpM = ∑ , , (2) where mu p and extmu E are the u th components of the m p and extm E , respectively. For a linear array of MNPs, Eq. (2) can be further split into decoupled equations in the form, extmn nmn EpM = ∑ , for either the transverse modes (with m p perpendicular to the chain axis) or longitudinal modes ( m p parallel to the chain axis). Since the transverse modes have richer properties in its band structure (such as negative slope and strong coupling to free photon ), we will focus on the transverse mode. The band structure can be calculated by setting 0 = extm E and finding pairs of ω and k so that Eq. (2) has non-trivial Bloch solutions of the form ikmam e pp = . More explicitly, we have to find the solutions of the equation, =− ka ωκωα , (3) where ),( k ωκ Σ−Σ+Σ= aikak hh for the transverse mode, n Σ )( )( akkin h eLi − = )( )( akkin h eLi + + , and )( xLi n is the polylogarithm function of order n . For lossless systems, ω and k are real numbers. Since the system we consider here has both radiation loss and absorption loss, real values of ω and k cannot satisfy Eq. (3). There are several ways to define the band structure of lossy systems but here we only focus on the complex band structure calculated by finding the complex k that corresponds to each real ω . A typical result calculated by numerical root searching in the complex- k plane is shown in Fig. 1. In the same figure, we compare the complex band structures for lossy ( γ (cid:61) = 0.02, results display as grey-colored lines) and lossless ( γ (cid:61) = 0, results display as blue-colored lines) metal with electric permittivity, ( ) 1 / [ ( )] p i ε ω ω ωω γ= − + , where γ is the electron scattering rate. Chosen for an easy comparison with the results in the literature, the plasma frequency for the Drude metal is p ω (cid:61) = 6.18 eV and the geometrical parameters are chosen to be
75 nm a = and
25 nm r = . Such parameters will reproduce the band structure in Ref. 13 outside the light cone (i.e., / h k c ω> ) in the case of γ (cid:61) = 0. The band structure has a negative slope for a range of k close to the zone boundary. It should be noted that such a band of negative slope always exists when the particles are close enough. The slope of the band will turn to positive when the inter-particle distance becomes larger. When there is no absorption, there is a critical frequency c ω ≡ (cid:61) k and k , where Re( ) Re( ) k k < ]. In this regime, the guided modes define a pass-band region because Im( ) Im( ) 0 k k = = , indicating that the plasmon excitation is truly guided along the chain without radiation loss. As frequency increases to c ω ω > , we have Re( ) Re( ) k k = and Im( ) Im( ) k k = − . (This is due to the time reversal symmetry.) This region is a band-gap region in which the plasmon modes are leaky and can radiate to the surroundings. When we have 0 γ ≠ (cid:61) , all solutions have Im( ) 0 k ≠ . (The two solutions are distinct because the time reversal symmetry is broken by absorptive dissipation.) For c ω ω < , the propagating plasmon experience a loss that is dominated by absorption. As frequency increases to c ω ω > , both absorption and radiation losses can suppress the plasmon propagation. Details on the features of the band structure are available elsewhere and we will not repeat the discussions. The aim of showing Fig. 1 is to show that the ratio between absorption and radiation losses can be adjusted by considering different frequencies and absorption rates. We can, therefore, find out the similarities and differences between the two kinds of losses by a suitable choice of parameters. III. LATTICE GREEN’S FUNCTION
We now consider the spatial decays of plasmonic excitations along the 1D array. The frequency-domain transverse lattice Green’s function, ( , , ) ( , )
G m n G m n ω ω ≡ − , for the coupled dipoles satisfies mnn mn nmGM δωω =− ∑ ),()( . (4) The response dipoles for a given external driving field are thus given by ( , ) extm nn p G m n E ω = − ∑ . ( , ) G m ω has the usual meaning of the response of the system when one particle is excited by a localized source ( ext extn n E E δ = ) and can be written as a Fourier integral: ∫ − = aa ikmaeig dkekamG // ),(2),( ππ ωαπω , (5) where ),( k eig ωα ]),()(/1[ −− −= ak ωκωα is the Green’s function in the reciprocal space. It should be noted that ),(),( mGmG −= ωω because of the mirror symmetry. While the complex band structure does not display the complete wave nature of the system (such as the exact response to the external field), the Green’s function keeps all the information about the waves in the system. Therefore, all properties shown in the complex band structure can be fully explained by investigating ( , )
G m ω and we can see the actual meaning of the complex k by comparing with ( , ) G m ω . There are at least two ways to evaluate ( , ) G m ω . We can substitute the analytical formula of ),( k eig ωα into Eq. (5) and evaluate the integral. However, a closed-form solution to the integral is not available and numerical evaluation of the integral may run into stability problems. Therefore, we choose to calculate, numerically, the inverse of M by truncating M to a matrix of finite dimension ( N N × ), which is the same as considering a finite chain of N particles. In order to preserve the mirror symmetry, we consider N as an odd number. Then, the Green’s function is ( , ) G m ω ( ) nn mn M δ − = ∑ ( ) m M − = , where m = N − , N − + , …, N − , N with ( 1) / 2 N N = − . Our method is similar to those in Refs. 7 and 8 except the single-site driving electric field is located at the center of the chain so that the mirror symmetry between the two ends of the chain is kept [see Fig. 2(a)]. In our calculations, we take N = 3001. Here, we first briefly highlight some features of ( , ) G m ω . When there is no absorption and radiation loss (i.e., γ = 0 and c ω ω < ), ( , ) G m ω does not decay. As shown in Ref. 8, the results are intuitively obvious and thus the results are not repeated here. However, when the wave experience mainly absorption loss (i.e., γ = 0.02 eV and c ω ω < ), the Green’s function, ( , ) G m ω , has two short-range exponential decays (which dominates for m < ) and a long-range power-law decay ( m > ) [see Fig. 2(b)]. Due to the beating between modes of different wave numbers, we can observe oscillations at the crossovers between different decays. In the case of pure radiation loss (i.e., γ = 0 and c ω ω > ), ( , ) G m ω has only one short-range exponential decay and a long-range power-law decay [see Fig. 2(c)]. In addition to the beating at the crossover points, there are strong oscillations for the whole exponential decay curve in Fig. 2(c). Such oscillations are associated with the time reversal symmetry and will be explained later in this paper. When radiation loss co-exists with absorption loss (i.e., γ (cid:61) = 0.02 eV and c ω ω > ), the result [see Fig. 2(d)] shows no qualitative difference from that of Fig. 2(b) except that the exponential decay lengths are shortened. These interesting features (single/double exponential decays with power-law decay) have been largely ignored (though they indeed exist if we examine the details) in previous published works. In the following, we will understand such features by making an ansatz for the form of ( , )
G m ω and the accurate meaning of the complex Bloch wave vector will emerge accordingly. IV. ANSATZ AND JUSTIFICATIONS
We will use an ansatz of the lattice Green’s function, ),( mG ω that takes the form, mmmans zmfzczcmG )(),( ++≡ ω , (6) where c , c , z , z , z are m -independent complex numbers satisfying =≤≤ zzz and f is a bounded even function satisfying = ∞→ mfm pm for some real number p . The first two terms in Eq. (6) will give short-range exponential decays while the third term represents the long-range power-law decay. We will justify the ansatz in the following, by comparing ),( mG ω and ),( mG ans ω and their generating functions in the complex number plane. By considering the generating functions, we will get extremely accurate values of the complex quantities (including c , c , z , z , and z ) without using any function-fitting method. To verify that ),( mG ω has to take the form as shown in Eq. (6), we consider the generating function of ),( mG ω , defined as ∑ ∞−∞= ≡ m m zmGzg ),(),( ωω , (7) where z is a complex number. By multiplying m z to Eq. (4) and summing the equation over the index m , one can show that the exact ),( zg ω reads ),(),( keg eigika ωαω = , (8) for any complex number k by allowing analytic continuation (as in the usual definition of the polylogorithm functions). Thus, we have a closed-from expression for ),( zg ω and it has a one-to-one correspondence to ),( k eig ωα , which is also the Green’s function in the reciprocal space. To present a concrete picture, we plot ),( zg ω in the complex- z plane (see Fig. 3) for various frequencies. From Fig. 3, we see that the exact generating function, ),( zg ω , has 4 poles (at z z = , z − , z , and z − ) and two branch cuts for each frequency. The poles are the solutions of Eq. (3) and the branch cuts, which show up as straight line segments, are defined by, Arg( ) z k a = for z < < and Arg( ) z k a = − for z > . These branch cuts are the principle branch cuts of the polylogarithm functions. The advantage of considering ),( zg ω as a function of z ( ) ika e = instead of ),( k eig ωα as 0 a function of k is that all the information are compressed within a finite unit circle, z ≤ , because of the mirror symmetry of the MNP chain [i.e., ),(),( − = zgzg ωω ]. In Figs. 3(a) and 3(b), we see that all four poles lie on the unit circle when there is no loss, i.e. c ω ω< and γ = . As we increase ω (keeping γ = ), z moves towards z − while z moves towards z − , and all four poles remain on the unit circle. After the poles meet at c ω ω= (not shown), they leave the unit circle separately [see Fig. 3(c) and 3(d)], keeping z z = < (due to radiation loss). It should be noted that we have z z ∗ = (due to time reversal symmetry) when γ = . When γ ≠ , the absorption breaks the time reversal symmetry, leading to z z ∗ ≠ and smooth changes in the positions of the poles as frequency changes [see Figs. 3(e)-3(h)]. All these results can be directly mapped into the complex band structure and since each z can be mapped into a wave number, k , the poles of the generating functions also give all the information required to produce the complex band structure in Fig. 1. Let us now consider the corresponding generating function for our ansatz Green’s function in Eq. (6), ( , ) ( , ) mans ansm g z G m z ω ω ∞=−∞ ≡ ∑ . Using Eq. (6), we can write ( , ) ( ) m mans m z z z zg z c c f m z zz z z z z z z z ω − − ∞− − =−∞ ⎛ ⎞ ⎛ ⎞ = − + − + ⎜ ⎟ ⎜ ⎟ − − − − ⎝ ⎠ ⎝ ⎠ ∑ . (9) It should be noted that, in writing Eq. (9), we have again used the analytic continuation as in the polylogarithm functions. From Eq. (9), we see that the two short-range exponential decays in ( , ) ans G m ω are associated with the four poles located at z z = , z − , z , and z − , which is consistent with ( , ) g z ω (as shown in Fig. 3). Here, we consider the exact 1 values of the poles and the corresponding residues of ( , ) g z ω . These quantities can be found accurately as long as we have rough estimates for the poles (for example, finding the poles from Fig. 3 by naked eye). With the estimated (and thus approximate) values of the poles, we can calculate a simple contour integral ( , ) / 2 C R g z dz i ω π= ∫ (cid:118) numerically around a closed curve, C , that encloses the pole but not crossing the branch cuts. We thus obtain the four residues (denoted by R , R , R , and R ) of ( , ) g z ω accordingly. It is expected that these four residues should be consistent with that of ( , ) ans g z ω . By equating the residues of ( , ) g z ω and ( , ) ans g z ω (i.e. R c z = ,
12 1 1
R c z − = − , R c z = , and
14 2 2
R c z − = − ), we can solve for c , c , z , and z . The values of z and z obtained in this way are surprisingly accurate in the sense that they pin the poles of ( , ) g z ω with a relative error within − even when the relative error of our initial estimated values is about 0.1. This verifies that the form of the ansatz is essentially exact. To further make sure that the complex quantities, c , c , z , and z , can describe the initial exponential decays of ( , ) G m ω accurately, we also compare m c z and m c z with ( , ) G m ω [see Figs. 2(b), 2(c), and 2(d)]. [We have c c ∗ = and z z ∗ = in the case without absorption loss and thus m c z and m c z overlap in Fig. 2(c).] The results deduced from the ansatz show an exact agreement to ( , ) G m ω . We note that in addition to the exponential decay(s), there are oscillations in the magnitude of ( , ) G m ω . We mentioned that those features near the crossover points are due to the beatings between different modes. Furthermore, there is an additional 2 oscillatory behavior for the whole range of exponential decay in the case of no absorption because we have m m c z c z + m c z = , which oscillates in its magnitude, when 0 γ= . To verify these statements, we subtract both m c z and m c z from ( , ) G m ω and plot the remaining function, ( , ) ( , ) m mrem G m G m c z c z ω ω≡ − − , (dashed lines) in Fig. 2. We can see that the oscillatory features are immediately removed after the subtraction and the remaining function is smooth in its magnitude. This is another evidence to show that the ansatz gives the correct functional form of the Green’s function. Finally, we have to check whether ( , ) rem
G m ω is of the form, ( ) m f m z , with z = . The verification of m z can be done by making a Fourier transform of the normalized function, ( , ) / ( , ) rem rem G m G m ω ω . We take a finite list of ( , ) / ( , ) rem rem
G m G m ω ω for m ≤ ≤ to calculate the Fourier transform numerically and find that the Fourier transform has only one single sharp peak located at h k k ≈ (free-space wave number), which means we have ( , ) / ( , ) h ik marem rem G m G m e ω ω ∼ . This is consistent with fixing h ik a z e = in our ansatz . Furthermore, if ( , ) ( ) mrem G m f m z ω = , we will have ( , ) ( ) rem
G m f m ω = . We can see in the insets of Figs. 2(b)-2(d) that ( , ) rem
G m ω is a smoothly decreasing function that behaves like 1 / p m for some constant p , manifesting as straight lines in the log-log plots (except for m < ). Again, this verifies the form of the ansatz . The long-range power-law decay has a wave number very close to that of the free-space photons, suggesting that such portion of decaying wave is not confined within the chain and is less related to the guided modes. The power factor, p , is found to be 3 close to 1.35, which is larger than the power factor for the Green’s function of photons in free space. V. DISCUSSION AND CONCLUSION
In summary, we have shown semi-analytically that the (transverse) lattice Green’s function for a periodic metal nanoparticle chain is of the form ),( mG ω m zc = m zc + amik h emf )( + , which contains both exponential and power-law decays. This is independent of whether the decay is due to radiation or absorption loss. Our results show that a non-zero imaginary part of the Bloch wave number k in the complex band structure only means an initial short-range exponential decay in the Green’s function for such lossy system. There exists an additional power-law long-range decay that has no apparent relation to the values of k of the guided mode. The existence of such power-law region suggests that extra care must be taken when one try to obtain the decay length by fitting the Green’s function with exponential function. Our results can be helpful in understanding the wave propagation in lossy open systems. In addition, one interesting and subtle feature is that there are two exponential decay constants in the transverse mode, and it is true for frequencies that can excite guided modes or the leaky modes at the band gap above guided mode frequencies. The existence of two (instead of one) decay constants is due to the breaking of time reversal symmetry by absorptive dissipation. If there is no absorption, the two decay lengths become equal, as required by time reversal symmetry. ACKNOWLEDGMENTS
This work was supported by the Central Allocation Grant from the Hong Kong RGC through HKUST3/06C. Computation resources were supported by the Shun Hing Education and Charity Fund. We thank Dr. Dezhuan Han, Dr. Yun Lai, and Prof. Zhaoqing Zhang for comments and useful discussions. 5
REFERENCES M. J. O. Strutt, Ann. d. Physik , 129 (1928); , 319 (1929); F. Bloch, Zeits. f. Physik , 555 (1928). L. P. Bouckaert, R. Smoluchowski, and E. Wigner, Phys. Rev. , 58 (1936). For photonic bands, see N. Stefanou, V. Karathanos, and A. Modinos, J. Phys.: Condens. Matter , 7389 (1992) and T. Suzuki and K. L. Yu, J. Opt. Soc. Am. B , 804 (1995). For electronic bands, see H. J. Choi and J. Ihm, Phys. Rev. B , 2267 (1999) and J. K. Tomfohr and O. F. Sankey, Phys. Rev. B , 245105 (2002). E. N. Economou, Green's Functions in Quantum Physics, 3rd ed. (Springer, Berlin, 2006). P. Sheng, Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena, 2nd ed. (Springer, Heidelberger, 2006). W. H. Weber and G. W. Ford, Phys. Rev. B , 125429 (2004). V. A. Markel and A. K. Sarychev, Phys. Rev. B , 085426 (2007). See, e.g., W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature , 824 (2003). The propagation of surface plasmon in a linear array of metal nanoparticles has been studied experimentally in Ref. 11 and theoretically in Refs. 7, 8, 13, 14, and 15. J.R. Krenn, A. Dereux, J.C. Weeber, E. Bourillot, Y. Lacroute, J.P. Goudonnet, G. Schider, W. Gotschy, A. Leitner, F.R. Aussenegg, and C. Girard, Phys. Rev. Lett. , 2590 (1999); R. de Waele, A. F. Koenderink, and A. Polman, Nano Lett. , 2004 (2007); A. F. Koenderink, R. de Waele, J. C. Prangsma, and A. Polman, Phys. Rev. B , 6 201403(R) (2007); K. B. Crozier, E. Togan, E. Simsek, and T. Yang, Opt. Express , 17482 (2007). M. Quinten, A. Leitner, J. R. Krenn, and F. R. Aussenegg, Opt. Lett. , 1331 (1998); M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Phys. Rev. B , R16356 (2000); Stefan A. Maier, Mark L. Brongersma, and Harry A. Atwater, Appl. Phys. Lett. , 16 (2001); S. Y. Park and D. Stroud, Phys. Rev. B , 125418 (2004); D. S. Citrin, Nano Lett. , 1561 (2004); R. A. Shore and A. D. Yaghjian, Electron. Lett. , 578 (2005); J. J. Xiao, K. Yakubo, and K. W. Yu, Appl. Phys. Lett. , 241111 (2006); A. F. Koenderink and A. Polman, Phys. Rev. B , 033402 (2006); A. A. Govyadinov and V. A. Markel, Phys. Rev. B , 035403 (2008); H. X. Zhang, Y. Gu, and Q. H. Gong, Chinese Phys. B , 2567 (2008). C. R. Simovski, A. J. Viitanen, and S. A. Tretyakov, Phys. Rev. E , 066606 (2005). A. Alù and N. Engheta, Phys. Rev. B , 205436 (2006). K. H. Fung and C. T. Chan, Opt. Lett. , 973 (2007). C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, New York, 1983). V. A. Markel, J. Opt. Soc. Am. B , 1783 (1995). The ),( k eig ωα defined here is the same as the eigenpolarizability defined in Ref. 15. The frequency dependence of every quantity (except m ) on the right-hand side of Eq. (6) is not notated for the sake of a tidy expression. 7 FIG. 1: (Color online) Complex band structure for a 1D array of MNPs. Left (right) panel show the real (imaginary) part of the wave number. The numbers inside the square boxes indicate different non-degenerate modes when absorption is included. FIG. 2: (Color online) Lattice Green’s function ( , )
G m ω . Panel (a) shows a schematic diagram with a dipole source located at the center of the MNP chain. Other panels show the values of ( , ) G m ω and the comparisons with the functions making up the ansatz for the case when energy loss is (b) dominated by absorption ( ω (cid:61) = 3.260 eV and γ (cid:61) = 0.02 eV), (c) purely due to radiation ( ω (cid:61) = 3.275 eV and γ (cid:61) = 0), and (d) contributed fairly by absorption and radiation ( ω (cid:61) = 3.275 eV and γ (cid:61) = 0.02 eV), respectively. The insets show the same graphs but in log-log scale, highlighting the long-range power law decay. Numbers with square boxes indicate the corresponding plasmon modes shown in Fig. 1. FIG. 3: (Color online) Contour plot of ( , ) / g z r ω at different frequencies. The poles at z z = , z , z − , and z − of the function are labeled. Upper (lower) panels show the results for the lossless (lossy) metal with γ (cid:61) = 0 ( γ (cid:61) = 0.02 eV). (a) and (e) ω (cid:61) = 3.260 eV. (b) and (f) ω (cid:61) = 3.270 eV. (c) and (g) ω (cid:61) = 3.275 eV. (d) and (h) ω (cid:61) = 3.280 eV.= 3.280 eV.