Optics Communications 282 (2009) 135–140
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Surface plasmon modes in multi-layer thin-films T.J. Davis 1 CSIRO Materials Science and Engineering Private Bag 33, Clayton Victoria 3169, Australia
a r t i c l e
i n f o
Article history: Received 17 June 2008 Received in revised form 12 September 2008 Accepted 16 September 2008
a b s t r a c t The surface plasmon modes in multi-layer thin-film structures of the metal–insulator–metal (MIM) type are modeled using a matrix method for the propagation of electromagnetic waves in stratified media. It is shown that the dispersion relation for the allowed surface plasmon modes is obtained from one of the matrix elements. The fundamental electromagnetic modes are the eigenfunctions of the differential equation for the magnetic field distribution and the eigenvalues are obtained from the dispersion relation. The expansion of an arbitrary wave profile in terms of the eigenfunctions is discussed and applied to the problem of surface plasmons propagating in a structure consisting of seven layers of alternating metal films and dielectrics. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Surface plasmon polaritons (often referred to as surface plasmons) provide a means for propagating optical energy in nano-scale structures. The surface plasmons arise from the coupling between the electromagnetic field and the free charges near the surface of a metal film. The wavelengths of the surface plasmons may be much smaller than those of freely propagating electromagnetic waves, making them ideal candidates for high frequency and compact electronic-photonic systems. The guided transmission of signals is obtained using surface plasmon waveguides based on a variety of geometries, including wedges [1], slots [2], strip-lines [3,4], rectangular waveguides [5] and metal–insulator–metal guides. Feigenbaum and Orenstein [6] reviewed the work of others in this field and also examined a variety of these structures using theoretical and numerical methods. An extensive review of the optics of surface plasmons on metal films and in waveguide structures was given by Zayats et al. [7]. The focus of this paper is the propagation of surface plasmons in multi-layered thin-film structures, of the metal–insulator–metal (MIM) type. The surface plasmons are associated with the natural electromagnetic modes of the structure. They are confined to the thin-film layers but are free to propagate over the surfaces or interfaces. The natural modes are obtained from the eigenvalue solutions to Maxwell’s equations subject to the boundary conditions imposed by the structure. The modes depend on the electric permittivities of the metal layers and are represented by complex numbers to properly account for losses on propagation. Dionne et al. [8] used complex permittivities in their analytical models of the MIM three-layer structure and described some of the associ-
E-mail address:
[email protected] Honorary Senior Researcher, School of Physics, Monash University, Victoria 3800, Australia. 1
0030-4018/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2008.09.043
ated electromagnetic modes, while Han and He [9] used a numerical method to analyze the splitting of modes at the junction between two MIM structures. The modes in a MIM nano-cavity were studied theoretically and experimentally by Kurokawa and Miyazaki [10,11] who showed that it exhibits Fabry–Perot-like resonance. A more general eigenfunction treatment was used by Sturman et al. [12] to analyze slits, one-dimensional photonic crystals and holes in metal films. Oulton et al. [13] use a semi-analytical mode matching technique to model the scattering of surface plasmons at abrupt interfaces for surface plasmons propagating along metal-dielectric interfaces (open waveguide structure). This involves the eigenmodes for the open waveguide structure that consist of bounded and unbounded (radiative) modes. The analysis becomes progressively more difficult as the number of layers is increased (e.g. MIMIM structures) leading to complex expressions for the dispersion of the surface plasmons. This paper addresses the problem of evaluating the eigenfunctions associated with multi-layer thin-film structures to model the behaviour of surface plasmons. The paper is divided into two parts. Firstly, the theory of eigenmodes in multi-layer thin-film structures is developed similar to the matrix method for stratified media [14] and it is shown how this leads to the dispersion relation for the allowed surface plasmons modes. Secondly, the method is demonstrated by modeling the propagation of electromagnetic waves in a seven layer structure using an expansion in terms of eigenfunctions. A natural consequence of the expansion is the alternating transfer of energy between two guiding layers within the structure. 2. Eigenmodes in multi-layered films Maxwell’s equations linking the electric and magnetic fields (E and H respectively) appropriate for the description of the surface plasmons are (in Gaussian units)
136
T.J. Davis / Optics Communications 282 (2009) 135–140
r Hðr; xÞ þ ikeðr; xÞEðr; xÞ ¼ 0 r Eðr; xÞ ikHðr; xÞ ¼ 0
ð1Þ ð2Þ
where the fields have a common time dependence expðixtÞ and k ¼ x=c is the free space wave number. The electric permittivity is e and the magnetic permeability is l = 1, since most magnetization effects are negligible at optical frequencies. The thin-film structure contains dielectric layers aligned with the x–y plane but stacked in the z-direction. Each layer has a constant electric permittivity leading to piece-wise continuous solutions. Since surface plasmons are excited only when there is a component of the electric field perpendicular to the surface [15], the solution is simplified by taking the electric field in the x–z plane and the magnetic field directed along the y-axis. The usual boundary conditions require that the tangential components of the electromagnetic fields are continuous between the layers. Only the magnetic field satisfies this condition completely (since the z component of the electric field is not tangential to the boundary), so that the eigenfunctions are found from the equations for the magnetic field. The eigenvalue equation is obtained by combining (1) and (2) and taking the x dependence as expðian xÞ; to represent a wave travelling in the positive x-direction. Then the z dependence of the magnetic field is described by
eðzÞ
d dz
1 dun ðzÞ 2 þ eðzÞk un ðzÞ ¼ a2n un ðzÞ eðzÞ dz
ð3Þ
In this equation, un ðzÞ is the nth eigenfunction and a2n is the nth eigenvalue. Since the structure is piecewise continuous, the solution to (3) for each layer l is conveniently written as
un ðzÞ ¼
X
Hþn;l expðicn;l ðz z0l ÞÞ þ Hn;l expðicn;l ðz z0l ÞÞ
ð4Þ
l
where the sum is over layers and z0l is the location of the lower boundary of layer l. The wavenumber cn;l is obtained from
c2n;l ¼ el k2 a2n :
ð5Þ
for each layer of permittivity el . When taking the square root of (5), there are two possible signs for cn;l . Since the coefficients Hþ n;l and H n;l are the amplitudes of upward and downward propagating waves respectively, the sign of cn;l is chosen such that the imaginary part is positive, to ensure that the waves decay as they propagate. Note that the coefficients have different values in different layers but are constant within each layer. They are set to zero when z lies outside their respective layer. The magnetic and electric fields are then given by
^ Hðr; xÞ ¼ y
X þ ðHn;l expðicn;l ðz z0l ÞÞ þ Hn;l expðicn;l ðz z0l ÞÞÞ l
expðian xÞ
ð6Þ
and
X 1 þ ^ an ^zÞ expðicn;l ðz z0l ÞÞ H ðc x kel n;l n;l l X 1 ^ þ an ^zÞ expðicn;l ðz z0l ÞÞ ð7Þ Hn;l ðcn;l x expðian xÞ k e l l
Eðr; xÞ ¼ expðian xÞ
Because the tangential components of the fields are continuous across the layers, an is constant for the entire structure and the coef ficients Hþ n;l and Hn;l in layer l are related to those in the adjacent layers. Matching the magnetic fields and the x components of the electric fields between layers l and l + 1 leads to the matrix equation
"
Hþn;lþ1 Hn;lþ1
#
"
¼
#"
1 1 1 1 elþ1 =cn;lþ1 cn;l =el cn;l =el 2 1 elþ1 =cn;lþ1 " #" þ # expðicn;l tl Þ 0 Hn;l 0 expðicn;l t l Þ Hn;l
#
ð8Þ
where layer l is of thickness of tl. Eq. (8) propagates the field amplitudes from one layer to the next and the magnetic fields in the entire structure are obtained by multiplying together the matrices for all the layers. If the lowest layer is labeled l = 0, then the field amplitudes in layer L are given by a matrix equation of the form
"
Hþn;L
#
¼
Hn;L
M11
M12
M21
M22
"
Hþn;0
#
Hn;0
ð9Þ
where the matrix elements are obtained from the set of matrix multiplications like in (8) for all layers. These depend on the electric permittivities for each layer and the eigenvalues a2n ; through cn;l and Eq. (5). Matrix formulations are commonly used for modeling wave propagation in stratified media (see e.g. [14]) but can take different forms depending on whether the field amplitudes or the field components form the basis set. In this paper the matrix is based on the amplitudes of the upward and downward propagating waves. The important point is that the free surface plasmon modes are those for which there are no incident fields [15], only outward directed fields (which are usually evanescent). If L is the top-most layer then the condition for no incident fields is Hþ n;0 ¼ H n;L ¼ 0. Expanding out (9) then gives a relation between the two outward directed fields, Hþ n;L ¼ M 12 H n;0 , and a constraint on one of the matrix . For there to be non-trivial solutions, it is elements, 0 ¼ M22 H n;0 necessary that M22 ¼ 0. This condition constrains the allowed values of a2n and gives rise to a discrete set of surface plasmon modes in the structure. This is the general dispersion relation linking an to the vacuum wavenumber k for the multi-layer film and represents the central result of this paper. There are usually an infinite number of modes because of the oscillatory nature of the factors in (8). The values a2n that solve the dispersion relation cannot be found in closed form for an arbitrary structure because of the complexity of the equation. Instead, numerical methods must be used and an example of one method is given in Section 3. However, it is straightforward to verify that M22 = 0 leads to the dispersion relation for surface plasmons in the cases of a single metal surface and a thin metal film. 2.1. Plasmon dispersion for a metal surface For a metal surface adjacent to a dielectric, L = 1 and expanding Eq. (8) gives M 22 ¼ ð1 þ e1 cn;0 =e0 cn;1 Þ=2. Equating this to zero and writing each cn;l in terms of an using (5) gives
a ¼ k
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
e0 e1 e0 þ e1
ð10Þ
There is only one eigenvalue and two possible values for a, corresponding to left and right propagating surface plasmons. This is the well-known expression for the wave number of a surface plasmon propagating along a metal-dielectric interface. 2.2. Plasmon dispersion for a metal film For a thin metal film (l = 1) between two dielectric layers (l = 0, 2), the matrix component can be expanded in the form
M22 ¼
expðicn;1 t 1 Þ fðe1 cn;0 þ e0 cn;1 Þðe2 cn;1 þ e1 cn;2 Þ 4e0 e1 cn;1 cn;2 þ ðe1 cn;0 e0 cn;1 Þðe2 cn;1 e1 cn;2 Þ expð2icn;1 t1 Þg
ð11Þ
The term in brackets is the surface plasmon dispersion relation for the metal film given by Raether [15]. This must be solved numerically and it yields an infinite number of eigenvalues. However, most of them have large imaginary components and the modes decay rapidly with distance.
T.J. Davis / Optics Communications 282 (2009) 135–140
2.3. Eigenfunctions and orthonormal expansions Once the eigenvalues are found (5), (8), and (9) can be employed to calculate the eigenfunctions (4). These are useful for representing an arbitrary surface plasmon in the structure in terms of the natural modes. However, the eigenfunctions are not orthogonal because the electric permittivity is a complex number and the differential equation (3) is not self-adjoint. This problem was investigated thoroughly by Botten et al. [16] in their description of the finitely conducting lamellar diffraction grating. Even though the thin-film structure is not periodic, Eq. (3) applies to both periodic and non-periodic structures. The main difference is that, for the grating, the boundaries are closed and periodic, whereas for the multi-layer structure, the boundaries are open and the fields are either evanescent or are matched to propagating waves. Botten et al. [16] showed that a function un ðzÞ defined by (4) is associated with an adjoint function equal to the complex conjugate, uþ n ðzÞ ¼ un ðzÞ. In this case, the orthogonality condition is found with respect to an inner product defined by
huþj ; un i ¼
Z
þ1
1 þ1
Z
uþj ðzÞ uk ðzÞ=eðzÞ dz
uj ðzÞuk ðzÞ=eðzÞ dz ¼ djk
ð12Þ
1
which has the inverse of the complex electric permittivity as a weighting function. The functions are normalized to unity. This expression is necessary for the proper expansion of functions in terms of the natural modes and has been used for analyzing a number of plasmonic structures [12,13]. With condition (12), any function f that satisfies the boundary conditions can be written as a sum of eigenfunctions
X
f ðzÞ ¼
an un ðzÞ
ð13Þ
n
with coefficients given by
an ¼
Z
þ1
ðf ðzÞun ðzÞ=eðzÞÞdz
ð14Þ
1
The expansion (13) with coefficients (14) applies to functions confined within the bounds of the multi-layer structure. The eigenfunctions do not form a complete set over the entire space because the eigenvalues are not continuous. In general, propagating waves outside the multilayer structure are not represented by these expansions and need to be included separately, as discussed by Oulton et al. [13]. In this paper the application of (12)–(14) is re-
137
stricted to surface plasmon fields confined to the multi-layer structure. In the following section, an example of the surface plasmon modes for a complex multi-layer film is given. 3. Surface plasmons in a multi-layer structure The methods of Section 2 are used here to model surface plasmon modes in a complex multi-layer structure consisting of three thin gold films separated by insulating layers Fig. 1a forming an IMIMIMI structure. The lower and upper gold films are adjacent to glassy materials and the central gold film is separated by low refractive index layers (in this case MgF2). These low refractive index layers act as waveguides and will be examined in more detail below. 3.1. Dispersion relation The zeroes of the complex function M22 in (9) are found by systematically calculating its magnitude and looking for minima while stepping through a range of complex a values. In the vicinity of each minimum, the eigenvalue is found using a complex version of Newton’s method [17] to locate the point where M 22 ¼ 0. In practice, the method is stopped when the magnitude of the residual d in M22 ¼ 0 þ d is reduced below 1010. The method works by adjusting the complex value of a to shift the solution along the gradient in solution space towards the zero. It was found that the method can miss one of a pair of eigenvalues that are degenerate or are almost degenerate. However, once an eigenvalue is found, the degenerate eigenvalues can be located by exploring the region about each known eigenvalue. The numerical procedure was modified so that the value of a was stepped by a small amount Reda and Imda away from the known eigenvalue and Newton’s method restarted. In effect, this forces the search to approach the eigenvalue position from different directions in the complex a plane, greatly improving the chance of locating additional eigenvalues. Finding a robust and reliable numerical method still remains a problem. The value of jM22 ðaÞj can be represented by a colour in the complex a plane, producing an image to visualize the distribution of eigenvalues Fig. 1b. The eigenvalues depend on frequency, which is represented in terms of a vacuum wavelength k ¼ 2p=k. In this example, the vacuum wavelength is 780 nm. The permittivity of gold at this wavelength is 21.19 + 0.7361i [18]. The permittivities of glass (2.310) and of the MgF2 layers (1.891) were taken to be real (that is, lossless). There are an infinite number of eigenvalues that
Fig. 1. (a) The geometry of the multi-layer film. It consists of three gold films (layers 1, 3 and 5), two MgF2 films (layers 2 and 4) sandwiched between two glass layers, forming an IMIMIMI structure. (b) An ‘‘image” of the value of the matrix element |M22| in the complex a plane corresponding to a free space wavelength of 780 nm (k = 2p/k). The eigenvalues are located in the dark red regions (dark grey in print version).
138
T.J. Davis / Optics Communications 282 (2009) 135–140
Fig. 2. (a) Six eigenfunctions associated with the multi-layer structure. Only the real parts of modes 1, 2, 3 and 4 are shown, since the imaginary components are small. Both real and imaginary parts of modes 5 and 6 are shown – the real parts are almost identical. The eigenfunctions are normalized according to Eq. (12) and the eigenvalues are given in Table 1. (b) The representation of a Gaussian function as a sum of 30 eigenmodes. The dominant modes are 3 and 4 shown in (a) (see text for details).
form a repetitive sequence for this geometry. Most of the eigenvalues have large imaginary terms – these are highly lossy and do not propagate very far. The set that lie close to the real axis are the most important ones and represent the propagating waveguide modes. The first six eigenfunctions ordered in terms of the imaginary component of a are shown in Fig. 2a and the eigenvalues are given in Table 1. Included in the table are the corresponding wavelengths (kn ¼ 2p=Rean ) and a ‘‘decay length” representing the distance over which the field amplitude falls to 1/e. The disper-
Table 1 The real and imaginary parts of the eigenvalues for the first six surface plasmon modes associated with a vacuum wavelength of 780 nm Mode 1 2 3 4 5 6
Re an (nm1) 2
1.2969 10 1.2971 102 3.0454 102 3.2794 102 2.1254 104 1.2634 103
Im an (nm1)
kn (nm)
Decay length (nm)
2.7301 105 2.7644 105 3.7872 104 4.6749 104 5.4538 102 5.4604 102
484.5 484.4 206.3 191.6 29,560 4973
36,630 36,170 2641 2139 18.34 18.31
The table also shows the wavelength of each mode and the distance over which the mode amplitudes decay.
Fig. 3. The dispersion relation for the first four eigenmodes shown in Fig. 2a. All modes lie below the light line and are not able to propagate in vacuum.
sion relations for the propagating modes (1–4) are shown in Fig. 3. These modes all lie below the light line which means they are confined to the multi-layer structure and are not able to propagate in free space. Modes 1 and 2 are long range and are associated with the glassgold interfaces. These are not tightly confined but extend quite far from the interfaces. These modes are a symmetric–antisymmetric pair that have almost identical wavenumbers. This symmetry is broken if the top surface has a different refractive index from the bottom surface. Because the dielectric layers between the metal films are very thin, modes 5 and 6 are evanescent and their propagation is highly damped along the planes of the layers for all wavelengths. If layers 2 and 4 Fig. 1a were thick, these would be propagating modes associated with surface plasmons coupled across the inner surfaces of the metal layers. Instead, these modes are resonances associated with charge oscillations within the metal regions. The magnetic fields and the z-component of the electric fields associated with these modes oscillate 90° out of phase, which is a consequence of the large imaginary component of a (compare (6) with (7)). This phase difference is characteristic of a coupled capacitor–inductor configuration [19] in which a resonance occurs due to the phase shift between the electric force acting on the charge (capacitance) and the magnetic force acting on the electric current (inductance). The magnetic field profile is given by the eigenfunction in Fig. 2a which is seen to peak within the metal films. There are Ohmic losses associated with the flow of charge through the metal film leading to strong damping of these resonances. It is likely that all of the higher order modes will have this property since they are associated with eigenvalues a that have large imaginary components. Note that the real part of a for mode 5 is negative because the imaginary part, which is associated with loss on propagation, must be chosen positive. This means that the excitation of this mode is associated with a small backward flow of energy in the x-direction, as well as outward flows in ±z-directions. The component of the Poynting vector in the plane of the structure, propor (where the bar denotes a complex conjugate), tional to E H depends on a and is mainly imaginary. This represents the alter-
T.J. Davis / Optics Communications 282 (2009) 135–140
139
nating flow of reactive or stored energy [20] as would be found in a coupled capacitor–inductor circuit. While all these higher-order modes are necessary to properly represent an arbitrary field profile within the multi-layer structure, they are not necessary for the description of the propagating modes. 3.2. Waveguide modes The two modes of interest here are the ‘‘waveguide” modes, 3 and 4, corresponding to surface plasmons travelling along the low refractive index regions, supported by the gold surfaces. These have short wavelengths but are also quite lossy, having decay distances of a few microns. They are distributed across both of the waveguide regions and form a symmetric/antisymmetric pair. Consequently, any wave propagating in one of the waveguide sections must be a combination of both modes. For example, consider a field with a Gaussian profile at the left of the lower dielectric region (layer l = 2 in Fig. 1a) chosen with a unit amplitude and a full width at half maximum of 47 nm, centred 80 nm above layer l = 0 to predominantly excite the waveguide modes. This field might be derived experimentally using the evanescent field from a nanoscale aperture placed just to the left of the multi-layer structure. The Gaussian and its representation as a sum of the first 30 modes are shown in Fig. 2b. The decomposition and reconstruction are calculated using (13) and (14). Because only a finite number of modes are used, there are some small errors in the representation of the Gaussian. In addition, the Gaussian will have some small values beyond the boundaries of the multi-layer structure that lie outside the domain where (13) and (14) are applicable (see the discussion in Section 2.3). Table 2 gives the mode amplitudes for the first 6 modes showing that only modes 3 and 4 are significant. The propagation of the modes is shown in Fig. 4 for the magnetic field and the z component of the electric field. The colour indicates the relative amplitude which decays over distance. Note the periodic transfer of energy from the lower guide to the upper guide and back. This oscillatory behaviour is a consequence of the wavelength difference between the two eigenmodes that are excited by the Gaussian. When the modes are in-phase, the energy is confined to the lower guide. When the modes are out of phase, the energy is distributed in the upper guide. The oscillation appears with a period given by Reða4 a3 Þs ¼ 2p, or s ¼ k3 k4 =ðk3 k4 Þ. For the modes in Table 1, the period is 2.689 lm, so the energy is transferred to the upper guide at 1.344 lm, as seen in Fig. 4. This behaviour is well-known from the theory of coupled optical fibres and is one mechanism for coupling optical energy between the waveguides [21]. In this example, the thickness of the gold layer l = 3 was chosen specifically to highlight this effect for the surface plasmons. By increasing the thickness of the layer, the coupling is reduced and the oscillation period increased. This is related to the ‘‘splitting” of the two waveguide modes 3 and 4 that become independent for very large layer thickness. It may be possible to use this effect to create a plasmon interferometer distributed within the depth of film. That is, one ‘‘arm” of Table 2 The first six expansion coefficients of the Gaussian wave in terms of surface plasmon eigenmodes Mode 1 2 3 4 5 6
Rean
Iman 2
1.839 10 1.402 102 1.404 1.426 2.965 103 4.908 103
Only modes 3 and 4 are significant.
6.603 104 3.804 104 8.152 103 6.467 103 6.942 103 7.491 103
Fig. 4. The magnetic and electric fields within the multilayer structure arising from an incident field associated with a free space wavelength 780 nm having a Gaussian profile (Fig. 2b). The incident field excites two modes which cancel in the upper guide and reinforce in the lower guide. Because of their different wavelengths, the modes eventually shift in phase which redistributes the energy into the upper region away from the lower region. The colour indicates the relative amplitude of the waves.
the interferometer lies in one plane (say layer l = 2) and the other arm in another plane (l = 4). By increasing the thickness of the middle layer (l = 3) with distance the two modes can be decoupled or by reducing the thickness again the two modes can be re-coupled. Perturbations to the phase of one of these modes, as caused by a change in the refractive index of the top layer, will change the interference between the modes when they are re-coupled. This will change the intensity distribution across the layers. While an analysis of this geometry is beyond the scope of this paper, the matrix method for determining the surface plasmon modes in the multi-layer thin-film structure provides a complete set of information about the modes, including the mode propagation constants, their losses and the distribution of magnetic and electric fields. 4. Conclusions In this paper, a method for finding the surface plasmon modes in a multi-layer thin-film structure is described. The central result is that the dispersion relation for the surface plasmons is found using a matrix method by setting one of the diagonal terms to zero. The dispersion relation yields a set of eigenvalues that determine the allowed surface plasmon modes. An analytical expression for the surface plasmon eigenfunctions is given and it is described how these can be used to represent an arbitrary electromagnetic wave in the multi-layer structure. It was shown that the use of eigenmodes is important for modeling the behaviour of surface plasmons in such structures and an example was given demonstrating the periodic transfer of energy between two coupled waveguiding regions. While the emphasis in this paper has been on structures of the MIM type, the method is applicable to any multi-layer structure that contains a metallic layer capable of supporting surface plasmons. The advantage of this numerical method over finite difference methods such as FDTD is that the eigenmodes of the structure are found independently of the nature of the driving the field.
140
T.J. Davis / Optics Communications 282 (2009) 135–140
References [1] D.K. Gramotnev, D.F.P. Pile, Appl. Phys. Lett. 85 (26) (2004) 6323. [2] N-N. Feng, M.L. Brongersma, L. Dal Negro, IEEE J. Quantum Electron. 43 (6) (2007) 479. [3] J.J. Ju, S. Park, M.-S. Kim, J.T. Kim, S.K. Park, Y.J. Park, M.-H. Lee, Appl. Phys. Lett. 91 (2007) 171117. [4] W. Barnes, A. Dereux, T.W. Ebbesen, Nature 424 (2003) 824. [5] F. Kong, B.-I. Wu, H. Chen, J.A. Kong, Opt. Express 15 (19) (2007) 12331. [6] E. Feigenbaum, M. Orenstein, J. Lightwave Technol. 25 (9) (2007) 2547. [7] A.V. Zayats, I.I. Smolyaninov, A.A. Maradudin, Phys. Rep. 408 (2005) 131. [8] J.A. Dionne, L.A. Sweatlock, H.A. Atwater, A. Polman, Phys. Rev. B 73 (2006) 035407. [9] Z. Han, S. He, Opt. Commun. 278 (2007) 199. [10] Y. Kurokawa, H.T. Miyazaki, Phys. Rev. B75 (2007) 035411. [11] T. Miyazaki, Y. Kurokawa, Phys. Rev. Lett. 96 (2006) 097401. [12] B. Sturman, E. Podivilov, M. Gorkunov, Phys. Rev. B 76 (2007) 125104.
[13] R.F. Oulton, D.F.P. Pile, Y. Liu, X. Zhang, Phys. Rev. B 76 (2007) 035408. [14] M. Born, E. Wolf, Principles of Optics, seventh ed., Cambridge University Press, Cambridge, 1999. [15] H. Raether, Surface plasma oscillations and their applications, in: G. Hass, M.H. Francombe, R.W. Hoffman (Eds.), Physics of Thin Films, vol. 9, Academic Press, New York, 1977, p. 145. [16] L.C. Botten, M.S. Craig, R.C. McPhedran, J.L. Adams, J.R. Andrewartha, Opt. Acta 28 (1981) 1087. [17] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes – The Art of Scientific Computing, Cambridge University Press, Cambridge, 1986. [18] CRC Handbook of Chemistry and Physics, 87th edition, 2006–2007. [19] T.J. Davis, Modeling and fabrication of tuned circuits for optical metamaterials, in: D. Abbott, Y.S. Kivshar, H. Rubinsztein-Dunlop, S. Fan (Eds.), Proc. SPIE, 6038, Photonics: Design, Technology, and Packaging II, 2006, p. 191. [20] J.D. Jackson, Classical Electrodynamics, second ed., Wiley, New York, 1975. [21] D. Marcuse, Light Transmission Optics, Van Nostrand Reinhold, New York, 1973.