A numerical study of the complex propagation constants of a surface plasmon polariton gap waveguide (SPGW) that was nanometric in size is performed by the method of lines (MoL). The validity of the code based on the MoL is examined by comparing the present results with those calculated using a volume integral equation, which is a completely different numerical technique from the MoL. The dependences of the complex propagation constants on the sizes of the SPGWs are investigated in detail and the fundamental propagation characteristics of SPGWs are revealed. Three kinds of SPGW structures (slab-slab, slab-plate and staggered slab-slab) are examined with a view to reducing the attenuation constants and the spot size to nanometric size. It is found that the nanometric field confinement can be controlled by using the staggered slab-slab structure of SPGW without a large change in the propagation constants.
©2009 Optical Society of America
Construction of optical circuits on nanometric scales has attracted the attention of many researchers in the fields of nanophotonics and nanooptics [1-3]. The diffraction limit of light makes it difficult to construct optical devices with dimensions that are much smaller than optical wavelengths and that have much higher integration densities than current optical integrated circuits. Recent theoretical and experimental studies have shown that optical circuits that use surface plasmon polaritons (SPPs) are promising candidates for future optical integrated circuits. Optical waveguides based on SPPs can be miniaturized much further than conventional diffraction-limited optical waveguides [4-7], opening up the possibility of developing nanometric integrated optical circuits. Although SPPs can travel no more than a few micrometers before extinguishing, such distances will be sufficiently long for future nanometric integrated optical circuits. Many interesting experimental and theoretical studies into practical and specific nanometric optical circuits that employ SPPs have been reported [8-12].
The possibility of practical optical circuits of nanometric size using a guiding structure known as a SPP gap waveguide (SPGW) was first proposed in ref. . Numerical simulations have demonstrated that SPGWs have the potential to guide, divide, and bend optical waves with acceptable losses in nanometric optical circuits [11, 12]. Many researchers are interested in SPP waveguides that have similar structures to SPGWs and have reported many interesting results [13-25]. The waveguide mechanism of SPGWs involves controlling the phase velocity in nanometric gap regions by varying the gap size.
In the design of nanometric optical circuits using SPGWs, one of the most fundamental and important parameters is the propagation constant of the guided modes. Since metals are dielectrics with complex permittivities in the optical frequency region, the propagation constants of guided modes in a SPGW inevitably have complex values; i.e., they always consist of both an attenuation constant and a phase constant. Furthermore, the electromagnetic fields penetrate into the metal walls of the SPP waveguide making it difficult to derive a rigorous or analytical solution of guided modes, even for SPGWs having simple cross sections. For SPGWs that have a limited structure where the gap-width difference between wide- and narrow- gap regions (a x and b x in Fig. 1) is not so large, approximate propagation constants can be obtained by the effective index modeling (EIM) . Numerical techniques have been applied to the calculation of the complex propagation constants of SPGWs and similar waveguides. The finite-difference time-domain method (FDTD), finite-different frequency-domain method (FDFD) and finite-element frequency-domain method (FEFD) have been employed [13-24]. The authors have also calculated the complex propagation constants of SPGW by three-dimensional numerical simulations based on the volume integral equation (VIE) method .
The method of lines (MoL) is a well-known and effective numerical technique and has been explored in a number of studies for calculating the propagation constants, including those of SPP waveguides [27-29]. The authors have also successfully used the MoL to calculate the propagation constants of SPPs in a metallic hollow rectangular waveguide with high accuracy . The MoL is expected to give more accurate results than the FDTD, because it partially employs analytical solutions, whereas FDTD is a purely numerical method. However, to the best of our knowledge, the computation of the propagation constants of the SPGWs by the MoL has not been reported. Furthermore, the important and fundamental propagation characteristics of SPGWs have not been sufficiently revealed.
In this paper, the fundamental and unusual propagation characteristics of SPGWs are investigated by the MoL. Three kinds of SPGW structures (slab-slab, slab-plate and staggered slab-slab) are examined with a view to reducing the attenuation constants and the spot size of the field distribution. It is found that the nanometric field confinement can be controlled by using the staggered slab-slab structure of SPGW without a large change in the propagation constants.
2. Geometry of the problem
The geometry of the problem considered in the present study is shown in Fig. 1. Although a lot different structures of SPGWs could be considered, we first consider the I-shaped cross section of the SPGW shown in Fig. 1. A hole structure in which a rectangular narrow-gap region is sandwiched between two rectangular wide-gap regions is created in the metallic material whose size is given by A×B shown in Fig. 1. The narrow- and wide-gap regions have dimensions of ax × ay and b x × b y, respectively, in the x-y plane shown in Fig. 1. In this paper, we refer to a x and a y as the narrow-gap width and depth, respectively, and we refer to b x and b y as the wide-gap width and height, respectively. The permittivities of the metal and space in the gap regions are given by ε 2 and ε 1, respectively. In this paper, we calculate the complex propagation constant of the SPP mode that propagates along the z-direction in Fig. 1 and consider the fundamental mode only. We employ the MoL to solve the complex generalized eigenvalue problem in the two-dimensional cross section of the SPGW shown in Fig. 1. The whole size of the metallic material used in this paper is given by A=3b x and B=2(2b y+a y) shown in Fig. 1.
3. Method of lines
The formulations of the MoL used in this paper are based on the formulation in ref. . The following field equations for the electromagnetic fields E and H by the two potentials ∏ e and ∏ h under the assumption of a harmonic time dependence exp(+/jωt) form the basis for the numerical analysis of the waveguide:
where μ 0 is the permeability in a vacuum, k 0 = ω(ε0 μ 0)1/2 is the free-space wave number, and η 0 = (μ 0/ε 0)1/2. The distribution of the relative permittivity in the x-y plane in Fig. 1 is independent of z and is assumed to be given by εr (x) = εr (x,y) in Eq. (1). The two potentials ∏ e and ∏ h are defined as
where ix is the unit vector in the x direction in Fig. 1 and jkz=(α+jβ)k 0 is the complex propagation constant to be determined numerically, where α and β are normalized attenuation and phase constants, respectively. Two scalar functions, ψ e and ψ h, satisfy the wave equations, and we discretize these equations using the MoL. We enclose the region of metallic material that contains the SPGW (see Fig. 1) with electric and/or magnetic walls in order to give fixed boundary conditions. The guiding structure with a dielectric constant of εr(x) = εr(x,y) is approximated by a set of single layers. These layers show the relative permittivities, which depend on only one transverse coordinate εr (x, const.). The x axis and the function εr(x, const.) are discretized using two shifted line systems parallel to the y axis. Field and wave equations, which are partial differential equations, are set up. These equations are subject to the boundary conditions at the lateral walls. A partial discretization of the potentials and field components is performed. The partial replacement of derivatives with finite-difference formulas yields systems of coupled ordinary differential equations. Applying suitable difference operators simultaneously satisfies the lateral boundary conditions. A transformation of the discretized potentials diagonalizes the systems of equations. The obtained uncoupled equations can be solved. Relationships between the tangential fields at the upper boundary of a layer and those at the lower boundary of a layer are obtained. Matching of the tangential fields, which must be continuous across all interfaces of the layered structure, yields a dispersion relation. The last field matching condition, applied near the center of the structure can be derived as follows :
where the column vector Ē is composed of the transformed tangential components of the electric field strength E at a suitably chosen interface, and Ȳ(kz) is a matrix equation containing kz. The complex propagation constant of the modes kz can be obtained by searching for values that satisfy det[Ȳ(kz)] = 0. This process constitutes the generalized and complex eigenvalue problem. Once the propagation constant of a mode has been determined using Eq. (4), the field distributions of the guided modes can be easily generated
4. Comparison of numerical results by the MoL with those by the VIE method
Throughout this paper, we assume that the wavelength is fixed at λ = 633 nm, the metal that constitutes the SPGWs is gold (Au) with a relative permittivity of ε 2/ε 0= -13.2 – j1.08  and the space inside the waveguide is a vacuum ε 1/ε 0= 1. In order to check the validity of the code based on the MoL used in this paper, we compare the numerical results obtained with those given in ref. . In ref. , the authors first solve the problem of the excitation and propagation of SPPs in the SPGW for the case when b x→∞ and b y→∞ in Fig. 1 by a numerical technique based on the VIE method. Then, the complex propagation constants are calculated by applying a least squares fitting to the resultant field distributions inside the SPGW. Notice that the numerical technique based on the VIE method differs completely from that based on the MoL used in this paper. A comparison between the results obtained using the MoL and those obtained using the VIE method are shown in Figs. 2 and 3. The dependences of the normalized phase constant β and the attenuation constant α obtained by the MoL and by the VIE method on the narrow-gap width ax and depth ay are shown in Fig. 2. The waveguide structures are shown in the insets of these figures. The open and solid circles indicate the results obtained using the MoL and by the VIE method, respectively. The results obtained using the MoL agree well with those obtained by the VIE method. The agreement between the results obtained by the MoL and by the VIE method for the attenuation constant α is not as good as that for the phase constant in Fig. 2. This is due to the values of the normalized attenuation constant being smaller than those of the normalized phase constant. In Fig. 2, the ranges of the ordinates β and α are 1.0-2.2 and 0.0-0.2, respectively.
The electric field distributions calculated by the MoL and those by the VIE method are shown in Fig. 3 for the case of ax × ay = 101 nm×101 nm. The electric field distributions calculated by the MoL also agree with those calculated by the VIE method. The difference in the surrounding free space region is due to the incident waves used in the simulations by the VIE method. The results shown in Figs. 2 and 3 confirm the validity of the present code based on the MoL. Since the simulation based on the VIE method employed a rather complicated process to calculate the propagation constants, we expect that the results obtained by the MoL are more accurate than those obtained by the VIE method. The same code has been also applied to the analysis of a SPP hollow rectangular waveguide and its validity has been also verified by comparing the numerical results with those given in published papers . The computational time by VIE was about 24 hours to calculate one propagation constant, because a large scale 3D scattering problem must be solved. Contrary, since the analysis by MoL is 2D, it was about only 30 minutes.
5. The effect of the wide-gap region dimensions
We first investigate the effect of the dimensions of the wide-gap region on the propagation constants. The dimensions of the narrow-gap region were fixed at ax×ay = 20 nm×20 nm, and the dependence of the complex propagation constant on the wide-gap height by was calculated (see Fig. 4). Hereafter, the propagation length defined as L=λ/(4πα)  is used instead of the normalized attenuation constant α in evaluating the propagation attenuation. The dependences of β and L on the wide-gap height b y are shown respectively in Fig. 4(a) and (b) for five different values of the wide-gap width b x.
It is reasonable that, for large values of b x and b y, the calculated propagation constants are the same as those for the mode which propagate to the z-direction along the narrow-gap region and are independent of the dimensions of the wide-gap regions (see the results for b x=202 nm and 303 nm in Fig. 4). The normalized phase constant β can be less than unity for small wide-gap heights and widths, as shown in Fig. 4(a). This is because the SPP mode in the narrow-gap region will approximate the hollow waveguide mode . For example, it is considered that the waveguide cutoff condition is realized for the case of b x= 40 nm and b y< 60 nm, as shown Fig. 4.
When the b y is large, the propagation length L decreases (α increases) with an increase in the wide-gap width b x (see Fig. 4(b)). In this case, when the width-difference between narrow-and wide-gap regions becomes small, the difference between phase constants in the narrow-and wide-gap regions also becomes small. This causes the decrease of the confinement of the fields into the narrow-gap region i.e., the confined field in the narrow-gap region expand into the wide-gap region. So, these results show the possibility of reducing the attenuation constant by using a small wide-gap width. Figures 5 (a) and (b) show the dependences of β and L on the wide-gap width b x respectively for a large wide-gap height (b y = 600 nm) for three different values of the cross section of the narrow-gap region a x × a y. The case of the smallest value of b x in Fig. 5 corresponds to the case when b x = a x; i.e., a waveguide that consists of two plates. From the results shown in Fig. 5, the optimum value of the wide-gap width is the one that maximizes the propagation length L for a given cross section of the narrow-gap region. As the wide-gap width b x increases, the propagation constants approach those of the mode of the narrow-gap region in Fig. 5 and the results become independent of the wide-gap width and depth.
6. Slab-slab structure of a SPGW
When b x and b y become very large in the I-shaped cross section shown in Fig. 1, the SPP can be considered to propagate only along the narrow-gap region; i.e., the mode of the SPGW that consists of a slab-slab structure shown in the inset of Fig. 6. Hereafter, we investigate only the propagation constants of the propagation mode along the narrow-gap region. In this case, the narrow-gap depth a y is equal to the slab thickness, as shown in the inset of Fig. 6. The basic characteristics of the slab-slab structure have been already reported by the many papers [11-25]. We also investigate this case in detail in this paper.
We first show the dependence of the normalized phase constant β and the propagation length L on the slab thickness a y for three different values of the narrow-gap width a x (see Fig. 6). The straight lines show the results for the case of a y=∞ i.e., the results for the SPP which propagates along the gap between the two metallic plates. The color of each straight line represents the results for the narrow-gap width indicated by the same color. β decreases and L increases monotonically with an increase in the slab thickness a y. This is due to the SPP mode in the narrow-gap region being coupled with the short-range SPP mode along the slabs of the structure shown in the inset. Figures 2(b) and (d) show that the electric field component E y of the mode is asymmetric relative to the two surfaces of the slabs. This field distribution is consistent with the short-range SPP mode along the slab. It is well known that a decreasing the slab thickness increases the attenuation constant in the short-range SPP mode of the slab .
When the slab thickness a y increases, the results approach the case of a y→∞. However, even if the slab thickness ay is very large, the results will not coincide with those of a y→∞ (indicated by the straight lines). This characteristic is shown in the propagation length L shown in Fig. 6(b) and it is because the waveguide structure in the MoL calculation always contains a edge mode, which propagates along the four corners of the slabs.
We next show the dependences of β and L on the narrow-gap width a x for five values of the slab thickness ay in Fig. 7. The straight lines indicate the results for the case when a x→∞; i.e., the results for the SPP mode that propagates only along the edge of the slab (edge mode) [20, 23]. In this case, for very large narrow-gap widths a x, all the results are coincident with those of the edge mode (see Fig. 7). It is interesting that β and L do not monotonically change with an increase in the narrow-gap width a x. The narrow-gap width a x that gives the maximum propagation length L (i.e., minimum attenuation constant) for all slab thicknesses (a y=20, 40, 60, 81 and 101 nm) is evident from Fig. 7(b). Similar characteristics were first reported for the plate-plate structure and it was demonstrated that this is due to the field distribution between the slabs . The authors of this present study have also reported similar characteristics in the slab-slab structure of a SPGW . In this paper, it is found that the maximum propagation length depends on the slab thickness for the slab-slab structure.
Figure 7(a) shows that there is a narrow-gap width a x that gives a minimum phase constant for a narrow-gap depth a y of 20 nm. For the other values of the gap depth (i.e., a y=40, 60, 81 and 101 nm) is no minimum value for the phase constant in Fig. 7(a). Similar characteristics of β have been also reported in ref. . This result shows that β can be smaller than that of the edge-mode in the slab-slab structure for small values of a y and is related to the field distribution of the mode. This aspect requires further investigation.
7. Slab-plate structure of SPGW
It is possible to reduce the attenuation constant by increasing the slab thickness of the SPGW shown in Fig. 6(b). This suggests that the attenuation constant can be reduced by replacing one of the slabs in the slab-slab structure by a plate, as shown in the inset in Fig. 8 [12, 23]. The results for this slab-plate structure for a SPGW are shown in Fig. 8. The dependences of β and L on the narrow-gap width a x of the slab-plate structure are shown in Figs. 8(a) and (b) respectively, for five different slab thicknesses a y. The straight lines having the same colors indicate the results for the edge mode, in which the plate is removed.
The open squares represent the results for the slab-slab structures shown in Fig. 7. The maximum values of L are larger than those of the slab-slab structure and the narrow-gap widths a x that give the maximum values of L are smaller than those of the slab-slab structure in Fig. 7.
Typical distributions of the electric field components of the slab-plate structure are shown in Fig. 9. Since the same scale is used for all three electric field components in Fig. 9, it is clear that the main component of the electric field is the x-component, as is the case for the slab-slab structure . Comparing Re[E z(x,y)] shown in Fig. 9(c) with that of the slab-slab structure shown in Fig. 2(c), it is evident that Re[E z(x,y)] is small close to the plane. In waveguide theory, it is well known that E z(x,y) is related to the attenuation in the propagation in the z-direction. Thus, the small value of E z(x,y) near the plane causes a reduction in the attenuation constant compared with the slab-slab structure whose Re[E z(x,y)] is shown in Fig. 2(c). The field distribution is extended over a larger area than that of the slab-slab structure shown in Fig. 3.
8. Staggered slab-slab structure of SPGW
The results given above for the slab-slab and slab-plate structures of SPGW show that the strong confinement of the field generally results in large phase and attenuation constants in the SPGW. If thick slabs are used in the slab-slab structure of a SPGW, we can reduce the attenuation constant (i.e., increase the propagation length L) shown in Fig. 6(b). However, the field distribution will be extended in the y-direction, because the size of the narrow-gap region is larger in the y-direction in Fig. 1. This suggests a structure in which the narrow-gap region has a staggered slab-slab structure, as shown in the insets of Fig. 10. Since the slab thickness is unchanged in this structure, we expect that it will reduce the spot size of the field distribution without increasing the attenuation constant. The results for the staggered slab-slab structure are shown in Fig. 10. The dependences of the normalized phase constant β and the propagation length L on the narrow-gap width a x for five values of the overlap length S of the staggered slab-slab structure are shown in Figs. 10(a) and (b), respectively. The definition of the overlap length S is shown in the insets of Fig. 10. The slab thickness is fixed to a y =101 nm in Fig. 10. The propagation characteristics shown in Fig. 10 are not strongly dependent on the overlap length S; i.e., the propagation constants do not change much with a change in the overlap length S.
Typical distributions of the electric field component Re[E x(x, y)] are shown in Figs. 11(a), (b), (c) and (d), for S=101 nm, 60 nm, 40 nm and 20 nm respectively, for a x =40 nm and a y=101 nm. We can see that the size of the field distribution can be reduced in the gap by reducing the overlap length S. These results demonstrate that the staggered slab-slab structure of SPGW can realize a small spot size of the confined field without increasing the attenuation constant. For example, the propagation constant of the staggered slab-slab structure is k z=(1.532-j0.0317)k 0 (L= 4277 nm) for a x =40 nm, a y =101 nm and S=20 nm, shown in Fig. 11(d). A similar spot size can be achieved in the slab-slab structure with a x =40 nm and a y =20 nm, but in this case the propagation constant is k z=(2.07-j0.095)k 0 (L= 1666 nm). In this case, the propagation length L of the staggered slab-slab structure is 2.6 times longer than that of the slab-slab structure. The attenuation constant can be effectively reduced by using a staggered slab-slab structure when the narrow-gap width a x is small.
For the larger value of the slab thickness of a y=202 nm, the field distributions for S=202 nm and S=20 nm are respectively shown in Figs. 12(a) and (b), for a x=40 nm. The propagation constants for Figs. 12(a) and (b) are given by k z=(1.62-j0.0314)k 0, (L= 5040 nm) and k z=(1.48-j0.0284)k 0 (L= 5572 nm), respectively. The propagation constants for S=20 nm do not differ much from those for S=202 nm. However, there is a drastic reduction in the spot size of the confined field in the y-direction in the staggered slab-slab structure of SPGW for S=20 nm compared with that for S=202 nm, as Fig. 12 shows. It is interesting that the staggered slab-slab structure of SPGW can control the spot size of confined fields without changing the propagation constant by very much. The staggered slab-slab structure can be applied to produce an optical probe that generates a strong near field and has a large throughput.
A numerical study of the complex propagation constants of a SPGW by the MoL has been performed. A SPGW having a general I-shaped cross section was considered first, followed by SPGWs with slab-slab, slab-plate and staggered slab-slab structures. Fundamental and important propagation characteristics, propagation constants, and field distributions of these SPGWs were investigated numerically in detail. It was demonstrated that the nanometric field confinement can be controlled by using a SPGW with a staggered slab-slab structure without a large change in the propagation constant. The results of the present study have important implications for future theoretical and experimental studies of nanometric optical circuits, as well as for practical applications.
References and links
1. M. Ohtsu, K. Kobayashi, T. Kawazoe, T. Yatsui, and M. Naruse, “Principles of Nanophotonics,” Chapman & Hall (2008).
2. V. M. Shalaev and S. Kawata ed., Nanophotonics with Surface Plasmons (Elsevier Science Ltd. 2007).
3. M. Ohtsu and H. Hori “Near-Field Nano-Optics: From Basic Principles to Nano-Fabrication and Nano-Photonics,” Plenum Pub. Corp. (1999).
5. S. A. Maier, M. L. Brongersma, P. G. Kik, S. Meltzer, A. A. G. Requicha, and H. A. Atwater, “Plasmonics - A route to nanoscale optical devices,” Adv. Mater. 13, 1501–1505 (2001). [CrossRef]
6. W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon sub-wavelength optics,” Nature (London) 424, 824830 (2003). [CrossRef]
7. E. Feigenbaum and M. Orenstein, “Modeling of Complementary (Void) Plasmon Waveguiding,” J. Lightwave Technol. 25, 2547–2562 (2007). [CrossRef]
8. D. F. P. Pile, T. Ogawa, D. K. Gramotnev, T. Okamoto, M. Haraguchi, and M. Fukui, “Theoretical and experimental investigation of strongly localized plasmons on triangular metal wedges for subwavelength waveguiding,” Appl. Phys. Lett. 87, 061106 (2005). [CrossRef]
9. D. K. Gramotnev and D. F. P. Pile, “Single-mode subwavelength waveguide with channel plasmon-polaritons in triangular grooves on a metal surface,” Appl. Phys. Lett. 85, 6323 (2004). [CrossRef]
10. S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, J.-Y. Laluet, and T. W. Ebbesen, “Channeling surface plasmons,” Appl. Phys. A 89, 225–231 (2007). [CrossRef]
11. K. Tanaka and M. Tanaka, “Simulations of nanometric optical circuits based on surface plasmon polariton gap waveguide,” Appl. Phys. Lett. 82, 1158–1160 (2003). [CrossRef]
12. K. Tanaka, M. Tanaka, and T. Sugiyama, “Simulation of practical nanometric optical circuits based on surface plasmon polariton gap waveguides,” Opt. Express 13, 256–266 (2005). [CrossRef] [PubMed]
13. G. Veronis and S. Fan, “Bends and splitters in metal-dielectric-metal subwavelength plasmonic waveguides,” Appl. Phys. Lett. 87, 131102 (2005). [CrossRef]
14. D. F. P. Pile, D. K. Gramotnev, M. Haraguchi, T. Okamoto, and M. Fukui, “Numerical analysis of coupled wedge plasmons in a structure of two metal wedges separated by a gap,” J. Appl. Phys. 100, 013101 (2006). [CrossRef]
17. S.H. Chang, T. C. Chiu, and C. Tai, “Propagation characteristics of the supermode based on two coupled semi-infinite rib plasmonic waveguides,” Opt. Express 15, 1755–1761 (2007). [CrossRef] [PubMed]
19. G. Veronis and S. Fan, “Guided subwavelength plasmonic mode supported by a slot in a thin metal film,” Opt. Lett. 30, 3359–3361 (2005). [CrossRef]
20. D. F. P. Pile, T. Ogawa, D. K. Gramotnev, Y. Matsuzaki, K. C. Vernon, K. Yamaguchi, T. Okamoto, M. Haraguchi, and M. Fukui, “Two-dimensionally localized modes of a nanoscale gap plasmon waveguide,” Appl. Phys. Lett. 87, 261114 (2005). [CrossRef]
23. G. Veronis and S. Fan, “Modes of subwavelength plasmonic slot waveguides,” J. Lightwave Technol. 25, 2511–2521 (2007). [CrossRef]
26. K. Tanaka, M. Tanaka, K. Katayama, and D. Miyahara, “Propagation constants of guided waves in surface plasmon polariton gap waveguides excited through an I-shaped aperture,” Comptes Rendus Physique 9, 16–23 (2008). [CrossRef]
28. R. Pregla and W. Pascher “The Method of Lines,” in Numerical Techniques for Microwave and Millimeter-Wave Passive StructuresT. Itoh, Ed. New York: Wiley, 1989.
29. U. Rogge and R. Pregla, “Method of lines for the analysis of dielectric waveguides,” J. Lightwave Technol. 11, 2015–2020, (1993). [CrossRef]
31. David W. Lynch and W. R. Hunter, Handbook of Optical Constants of SolidsE. D. Palik, ed. (Academic, New York, 1985).
32. J. J. Burke, G. I. Stegeman, and T. Tamir, “Surface-polariton-like waves guided by thin, lossy metal films,” Phys. Rev. B 33, 5186–5201 (1986). [CrossRef]