Copyright © 2009 American Scientific Publishers All rights reserved Printed in the United States of America Journal of Computational and Theoretical Nanoscience Vol. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks Georgios Veronis1
∗ , Şükrü Ekin Kocabaş2 , David A. B. Miller2 , and Shanhui Fan2 1 Department of Electrical and Computer Engineering and Center for Computation and Technology, Louisiana State University, Baton Rouge, LA 70803, USA 2 Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA We review some of the recent advances in the simulation of plasmonic devices, drawing examples from our own work in metal-insulator-metal (MIM) plasmonic waveguide components and networks. We introduce the mode-matching technique for modeling of MIM waveguide devices. We derive the complete set of orthogonal modes that theby MIM waveguide Delivered Ingenta to: supports and use it to apply the mode-matching technique to the analysis of plasmonic waveguide networks. We also introduce University of Houston several different equivalent models for plasmonic waveguide components, such as the characteristic IP : 129.7.158.43 impedance model for deep subwavelength MIM waveguides, the scattering matrix description of Mon, 05 Oct 2009 14:12:14 MIM waveguide junctions, and equivalent circuit models. The model abstraction provided by these equivalent models is important for the analysis and synthesis of device functions, as illustrated with the design of a waveguide mode converter. Keywords: Plasmonics, Surface Plasmons, Mode-Matching, Equivalent Circuits, Scattering Matrix, Characteristic Impedance. CONTENTS 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. Challenges Involved in Modeling of Plasmonic Devices . . 3. General Purpose Simulation Methods for Plasmonic Waveguide Devices . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Finite-Difference Frequency-Domain Method . . . . . . 3.2. Finite-Difference Time-Domain Method . . . . . . . . . 4. Mode-Matching Method . . . . . . . . . . . . . . . . . . . . . . 4.1. Spectrum of MIM Waveguides . . . . . . . . . . . . . . . 4.2. Mode Orthogonality . . . . . . . . . . . . . . . . . . . . . 4.3. Discretization of the Continuous Spectrum . . . . . . . 4.4. Mode Completeness . . . . . . . . . . . . . . . . . . . . . 4.5. Convergence of the Mode-Matching Method . . . . . . 5. Equivalent Models for Plasmonic Waveguide Components 5.1. Characteristic Impedance Model . . . . . . . . . . . . . . 5.2. Scattering Matrix Description of Junctions . . . . . . . 5.3. Circuit Model for the Waveguide Junction . . . . . . . 5.4. Cascade Connection of Junctions . . . . . . . . . . . . . 6. Summary and Conclusions . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1808 . . . 1810 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1812 1812 1813 1814 1814 1815 1817 1817 1817 1818 1819 1819 1821 1822 1824 1824 1. INTRODUCTION Surface plasmons are electromagnetic waves that propagate along the interface of a metal and a dielectric. In surface plasmons light interacts with the free electrons of ∗ Author to whom correspondence should be addressed. 1808 J. Comput. Theor. Nanosci. 2009, Vol. 6, No. 8 the metal which oscillate collectively in response to the applied field.1 Recently, nanometer scale metallic devices have shown the potential to manipulate light at subwavelength scales using surface plasmons. This could lead to photonic circuits of nanoscale dimensions. The use of nano-metallic structures could also bridge the size mismatch between modern electronic components with critical dimensions on the order of tens of nanometers and the micrometer scaled optical devices.2 Plasmonic waveguides have shown the potential to guide subwavelength optical modes at metal–dielectric interfaces. Several different plasmonic waveguiding structures have been proposed,3–8 such as metallic nanowires4
5 and metallic nanoparticle arrays.6
7 Most of these structures support a highly-confined mode only near the surface plasmon frequency. In this regime, the optical mode typically has low group velocity and short propagation length. It has been shown however that a metal-insulator-metal (MIM) structure with a dielectric region thickness of ∼100 nm supports a propagating mode with a nanoscale modal size at a wavelength range extending from DC to visible.9 Thus, such a waveguide could be potentially important in providing an interface between conventional optics and subwavelength electronic and optoelectronic devices. There have been several theoretical studies of MIM waveguides in the literature.9–18 Because of the 1546-1955/2009/6/1808/019 doi:10.1166/jctn.2009.1244 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. Georgios Veronis received the B.S. degree in electrical engineering from the National Technical University of Athens, Greece, in 1997, and the M.S. and Ph.D. degrees in electrical engineering from Stanford University, Stanford, CA, in 1999 and 2002 respectively. He was an engineering research associate at the Ginzton Laboratory at Stanford University. He is currently an assistant professor jointly at the Department of Electrical & Computer Engineering and the Center for Computation & Technology at Louisiana State University, Baton Rouge, LA. He has published more than 25 refereed journal papers, and is the holder of 1 U.S. patent. His research interests include the theoretical analysis of nanophotonic and plasmonic devices, and computational electromagnetics. Dr. Veronis is a member of IEEE and OSA. Şükrü Ekin Kocabaş received the B.S. and M.S. degrees in electrical engineering from Bilkent University, Ankara, Turkey and Stanford University, Stanford, CA, USA in 2002 and 2004 respectively. Currently, he is a Ph.D. student majoring in electrical engineering with a minor in physics. He has been a member of IEEE since 2000. His main research interest is modeling and fabrication of nanophotonic optoelectronic devices. Delivered by Ingenta to: University of Houston IP : 129.7.158.43 Mon, 05 Oct 2009 14:12:14 David A. B. Miller received the B.Sc. degree from St Andrews University, and the Ph.D. degree from Heriot-Watt University, in 1979, both in Physics. He was with Bell Laboratories, from 1981 to 1996, as a department head from 1987, latterly of the Advanced Photonics Research Department. He is currently the W. M. Keck Professor of Electrical Engineering at Stanford University, Stanford, CA, the Director of the Solid State and Photonics Laboratory at Stanford and a Co-Director of the Stanford Photonics Research Center. He also served as Director of the Ginzton Laboratory at Stanford from 1997–2006. His research interests include nanophotonic and quantum-confined optoelectronic physics and devices, and fundamentals and applications of optics in information sensing, switching, and processing. He has published more than 200 scientific papers, delivered more than 100 conference invited talks, and holds 62 patents. Dr. Miller has served as a Board member for both the Optical Society of America (OSA) and the IEEE Lasers and Electro-Optics Society (LEOS), and in various other society and conference committees. He was President of the IEEE Lasers and Electro-Optics Society in 1995. He has also served on Boards for various photonics companies. He was awarded the Adolph Lomb Medal and the R. W. Wood Prize from the OSA, the International Prize in Optics from the International Commission for Optics, and the IEEE Third Millennium Medal. He is a Fellow of the Royal Societies of London and Edinburgh, OSA, APS, and IEEE, and holds honorary degrees from the Vrije Universiteit Brussel and Heriot-Watt University. Shanhui Fan is an associate professor of electrical engineering at Stanford University. He received his Ph.D. in 1997 in theoretical condensed matter physics from the Massachusetts Institute of Technology (MIT), and was a research scientist at the Research Laboratory of Electronics at MIT prior to his appointment at Stanford. His research interests are in computational and theoretical studies of solid state and photonic structures and devices, especially photonic crystals, microcavities, and nanophotonic circuits and elements. He has published over 140 refereed journal articles, has given over 100 invited talks, and was granted 28 US patents. Professor Fan received a National Science Foundation Career Award (2002), a David and Lucile Packard Fellowship in Science and Engineering (2003), the National Academy of Sciences Award for Initiative in Research (2007), and the Adolph Lomb medal from the Optical Society of America (2007). Dr. Fan is fellow of OSA and APS, a senior member of IEEE, and a member of SPIE. J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 1809 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. predicted attractive properties of MIM waveguides, people 2. CHALLENGES INVOLVED IN MODELING have also started to explore such structures experimentally. OF PLASMONIC DEVICES In particular, Dionne et al.19 have recently demonstrated Surface plasmons can be described by macroscopic waveguiding in a quasi-two-dimensional MIM geometry electromagnetic theory, i.e., Maxwell’s equations, if the experimentally, showing clear evidence of a subwavelength electron mean free path in the metal is much shorter than guided mode with substantial propagation distances. the plasmon wavelength. This condition is usually fulfilled Most of the theoretical studies of nanoscale plasmonic at optical frequencies.24 Assuming an expit harmonic devices involve the use of general purpose electromagtime dependence of all field quantities, Maxwell’s curl netic simulation techniques. These techniques are able equations in the frequency domain take the form to simulate devices of arbitrary geometry and material composition. On the other hand, as described in the  × Er = −i0 Hr (1) next section, these techniques typically require numerical  × Hr = i 0 r rEr (2) grids with resolution of tens to thousands of grid points per wavelength, depending on the application. They are In macroscopic electromagnetic theory, bulk material proptherefore not suitable for the design and optimization of erties, such as the relative dielectric constant r r, are multi-component optical circuits. One approach to address used to describe objects irrespective of their size. However, this problem is the development of special-purpose simulafor particles of nanoscale dimensions a more fundamental tion techniques with superior computational efficiency for description of their optical and electronic properties may a specific class of problems. Previously, for example, in the by Ingenta to:25 Delivered be required. modeling of photonic crystal devices, it has been shown of Houston University Analytical methods, such as Mie theory,26 can only be that the use of Wannier functions results in severalIPorders : 129.7.158.43 applied to planar geometries or to objects with very high of magnitude speedup in the optimizationMon, of photonic 05 Oct 2009 14:12:14 symmetry (spheres, infinite cylinders) and have therefore 20 crystal circuits. A second approach is the development limited importance in the analysis of plasmonic devices of equivalent models with only a few dynamic variables, and structures. Thus, the analysis of plasmonic devices is which are nevertheless capable of describing complex optimostly based on numerical simulation techniques. cal processes in photonic devices in detail. For example, Numerical modeling of plasmonic devices involves sevin the modeling of photonic crystal devices, it has been eral challenges specific to plasmonics which need to be shown that such models can be derived with the use of addressed. The dielectric constant of metals at optical coupled mode theory,21
22 and perturbation theory.22
23 wavelengths is complex, i.e., r  = Re  + i Im  In this context, here, we provide a review of some of and is a complicated function of frequency.27 Thus, sevour own recent research activities aiming to advance the eral simulation techniques which are limited to lossless, theory and simulation of plasmonic waveguide devices non-dispersive materials are not applicable to plasmonic through the development of efficient special-purpose simdevices. In addition, in time-domain methods the disperulation techniques and equivalent models. We first give sion properties of metals have to be approximated by an overview of the challenges involved in modeling of suitable analytical expressions.28 In most cases the Drude plasmonic devices, and we briefly examine two general model is invoked to characterize the frequency dependence purpose simulation techniques which are widely used for of the metallic dielectric function29 modeling of plasmonic waveguide devices. We then intro2p duce the mode-matching technique for modeling of MIM (3) r
Drude  = 1 −  + i  waveguide devices. We derive the complete set of orthogonal modes that the MIM waveguide supports and use where p , are frequency-independent parameters. Howit to apply the mode-matching technique to the analyever, the Drude model approximation is valid over a limsis of plasmonic waveguide networks. We show that this ited wavelength range.29 The range of validity of the Drude special purpose simulation technique is far more effimodel can be extended by adding Lorentzian terms to cient for this class of problems than general purpose elecEq. (3) to obtain the Lorentz–Drude model29 tromagnetic simulation techniques. Finally, we introduce k several different equivalent models for plasmonic waveg
fj 2j (4)  = + r
LD r
Drude 2 uide components, such as the characteristic impedance 2 j=1 j −   − i j model for deep subwavelength MIM waveguides, the scattering matrix description of MIM waveguide juncwhere j and j stand for the oscillator resonant tions, and equivalent circuit models. We show that the frequencies and bandwidths respectively, and fj are model abstraction provided by these equivalent models is weighting factors. Physically, the Drude and Lorentzian important for the analysis and synthesis of device functerms are related to intraband (free-electron) and interband tions, and we illustrate this with the design of a mode (bound-electron) transitions respectively.29 Even though converter. the Lorentz–Drude model extends the range of validity of 1810 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. grid resolution is required at the metal-dielectric interanalytical approximations to metallic dielectric constants, face to adequately resolve the local fields. In addition, it is not suitable for description of sharp absorption edges several plasmonic devices are based on components of observed in some metals, unless a very large number of subwavelength dimensions.1 In fact, most of the potential terms is used.29 In particular, the Lorenz-Drude model applications of surface plasmons are related to subwavecannot approximate well the onset of interband absorplength optics. The nanoscale feature sizes of plasmonic tion in noble metals (Ag, Au, Cu) even if five Lorentzian devices pose an extra challenge to numerical simulation terms are used.29 In Figure 1 we compare the Drude and techniques. Lorentz–Drude models with experimental data for silver. We illustrate the challenges involved in modeling plasWe observe that even a five-term Lorentz–Drude model monic devices using a simple example. We consider an with optimal parameters results in a factor of two error infinite periodic array of silver cylinders illuminated by a at certain frequencies. An alternative approach is the use plane wave at normal incidence (inset of Fig. 2(a)). We use of an analytical expression based on multiple complexthe finite-difference frequency-domain method, described conjugate pole-residue pairs.30 It has been shown that, in more detail below, to calculate the transmission of when such an approach is used with time-domain meththe periodic array. This method allows us to directly use ods, it can lead to substantial savings in both memory and experimental data for the frequency-dependent dielectric computation time.30 constant of metals, including both the real and imagiIn addition, in surface plasmons propagating along nary parts, with no further approximation. The fields are the interface of a metal and a dielectric, the field is discretized concentrated at the interface, and decays exponentially Delivered by Ingenta to:on a uniform two-dimensional grid with grid size x = away from the interface in both the metal and dielecUniversity of Houston y = l. In Figure 2(a) we show the calculated transmission as a function of frequency. We also tric regions.1 Thus, for numerical methods based IP on: dis129.7.158.43 show the transmission of the structure calculated with cretization of the fields on a numerical grid, a very fine Mon, 05 Oct 2009 14:12:14 (a) (a) 0.9 0.8 102 k0 0.7 |εRe| Transmission 2a 101 |εIm| 100 0.6 0.5 a 0.4 0.3 0.2 Drude 0.1 10−1 0 200 400 600 800 200 300 1000 400 500 600 700 800 900 Frequency (THz) Frequency (THz) (b) (b) 0.85 102 |εRe| Transmission 0.8 1 10 |εIm| 100 0.75 0.7 0.65 0.6 −1 10 0.55 0.5 200 400 600 800 1000 1 2 5 10 ∆l (nm) Frequency (THz) Fig. 1. Real and imaginary part of the dielectric constant of silver at optical frequencies. The solid lines show experimental data.27 The dashed lines show values calculated using (a) the Drude model, (b) the Lorentz– Drude model with five Lorentzian terms. The parameters of the models are obtained through an optimization procedure.29 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Fig. 2. (a) Calculated transmission spectrum of an infinite array of silver cylinders (shown in the inset) for normal incidence and TM polarization. Results are shown for a = 100 nm. The dashed line shows the transmission spectrum calculated using the Drude model with parameters p = 137 × 1016 sec−1 , = 729 × 1013 sec−1 . (b) Calculated transmission at 855 THz as a function of the spatial grid size l. 1811 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. the Drude model of Eq. (3). We observe that the use with node coordinates rijk = xi
yj
zk  is the simplest and most commonly-used. A field quantity at nodal location of the Drude model results in substantial error. In genrijk is denoted for convenience as fijk = f rijk . Based on eral, the Drude model parameters are chosen to minimize Eq. (5), the first derivative can be approximated by the the error in the dielectric function in a given frequency following central-difference formula range.31 However, this approach gives accurate results in  a limited wavelength range, as illustrated in this example. df  fi+1 − fi−1 In general, the complicated dispersion properties of met≃ (6) dx i 2x als at optical frequencies pose a challenge in modeling of plasmonic devices not encountered in modeling of low- or which is second-order accurate, based on the discussion high-index-contrast dielectric devices. above, if the rectangular grid is uniform, i.e., xi = ix. In Figure 2(b) we show the calculated transmission at Similarly, the second derivative can be approximated by a specific wavelength of 0 = 351 nm as a function of the formula  the spatial grid size l. We observe that a grid size of d 2 f  fi+1 − 2fi + fi−1 l ≃ 1 nm is required in this case to yield reasonably accu≃ (7) dx2 i x2 rate results. The required grid size is directly related to the decay length of the fields at the metal-dielectric interwhich is also second-order accurate on a uniform grid.28 face. In general, modeling of plasmonic devices requires By replacing derivatives in differential equations with much finer grid resolution than modeling of low- or hightheir finite-difference approximations, we obtain algebraic index-contrast dielectric devices, due to the high localiza- by Ingenta Delivered equationsto: which relate the value of the field at a spetion of the field at metal-dielectric interfaces ofUniversity plasmonic ofcific node to the values at neighboring nodes. To solve Houston devices. The required grid size depends on the shape Maxwell’s equations with the FDFD method, we discretize IP : and 129.7.158.43 feature size of the modeled plasmonic device, the the 14:12:14 system of the three coupled scalar partial differential Mon,metallic 05 Oct 2009 material used and the operating frequency. equations obtained from the wave equation for the electric field 3. GENERAL PURPOSE SIMULATION METHODS FOR PLASMONIC WAVEGUIDE DEVICES Here we briefly examine two general purpose electromagnetic simulation techniques which are widely used for modeling of plasmonic waveguide devices, the finite-difference frequency-domain (FDFD) and the finite-difference timedomain (FDTD) methods. We also examine how they address the challenges mentioned above. Other general purpose simulation techniques for modeling of plasmonic waveguide devices include the Green dyadic method,4 the finite-element method,32 and the method of lines.33 3.1. Finite-Difference Frequency-Domain Method In finite-difference methods, derivatives in differential equations are approximated by finite differences. To approximate the derivative df /dxx0 we consider Taylor series expansions of f x about the point x0 at the points x0 + x and x0 − x and obtain28  df  f x0 + x − f x0 − x = + Ox2  (5)  dx x0 2x Equation (5) shows that a central-difference approximation of the first derivative is second-order accurate, meaning that the remainder term in Eq. (5) approaches zero as the square of x. In finite-difference methods a continuous problem is approximated by a discrete one. Field quantities are defined on a discrete grid of nodes. The rectangular grid 1812  ×  × Er − r r 2 Er = −i0 Jr c2 (8) For simplicity we consider here two-dimensional problems with TE polarization. For TE polarization we have E = Ez ẑ and the wave equation for the electric field becomes32
34  2   2 2 + + k x
y Ez x
y = −i0 Jz x
y (9) 0 r x2 y 2 For simplicity we consider a uniform rectangular grid with xi = ix, yj = jy, and replace the derivatives in Eq. (9) with their finite-difference approximations of Eq. (7) to obtain fi+1
j −2fi
j +fi−1
j fi
j+1 −2fi
j +fi
j−1 + +k02 ri
j fi
j = Ai
j x2 y2 (10) where f = Ez and A = −i0 Jz . Thus, application of finite-difference approximations at the node location rij = xi
yj  results in a linear algebraic equation which relates the field fi
j to the fields at the four adjacent nodes fi+1
j , fi−1
j , fi
j+1 , fi
j−1 . By applying the finite-difference approximation to all nodes of the grid we obtain a system of linear equations of the form Ax = b, where b is determined by the source current J. Since the equation for the field at each point involves only the fields at the four (six in three dimensions, two in one dimension) adjacent points, the resulting system matrix is extremely sparse.34 FDFD can be used to model plasmonic devices with arbitrary geometries. In addition, FDFD is a frequencydomain technique and can thus treat arbitrary material J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. dispersion. Nonuniform and/or nonorthogonal grids are required in FDFD for efficient treatment of curved surfaces and rapid field decays at metal-dielectric interfaces. In FDFD, as in all other methods which are based on discretization of the differential form of Maxwell’s equations in a finite volume, absorbing boundary conditions (ABCs) are required so that waves are not artificially reflected at the boundaries of the computational domain.28
32 Very efficient and accurate ABCs such as the perfectly matched layer (PML) have been demonstrated.35 As mentioned above, FDFD results in extremely sparse systems of linear equations. Such problems can be solved efficiently if direct or iterative sparse matrix techniques are used.32
34 3.2. Finite-Difference Time-Domain Method The FDTD method28 solves directly Maxwell’s timedependent curl equations We observe that, using the Yee lattice, all spatial finitedifference expressions are central and therefore secondorder accurate. Similar finite-difference equations are obtained by discretizing the other components of Eqs. (11) and (12). In summary, FDTD is an explicit numerical scheme which is second-order accurate both in time and in space (in uniform media). 3.2.1. Treatment of Dispersive Media in FDTD One of the major challenges in FDTD modeling of metals at optical frequencies is the treatment of the metallic dispersion properties. As mentioned above, in time-domain methods the dielectric constants of dispersive media have to be approximated by suitable analytical expressions. The most common algorithm for modeling dispersive materials with FDTD is the auxiliary differential equation (ADE) method.28
36 In dispersive materials  relates E and D H Delivered D = E (15) (11) by Ingenta to: t University of Houston ADE is based on integrating an ordinary differential equaIP : 129.7.158.43 E  ×H = 0 r (12) tion14:12:14 in time that relates Dt to Et, concurrently with Mon, 05 Oct 2009 t Maxwell’s equations. This equation is derived by taking so that both space and time have to be discretized. the inverse Fourier transform of Eq. (15). The standard FDTD is based on the Yee algorithm.28 We consider here a simple example where the dielectric As mentioned in the previous section, central-difference constant r  consists of a single Lorentzian term, i.e., approximations are second-order accurate. To achieve second-order accuracy in time, the Yee algorithm uses a 20  = (16) 28 r leapfrog arrangement. E fields are calculated at t = nt 20 − 2  − i 0 using previously calculated and stored H fields. Then H If we substitute Eq. (16) into Eq. (15) and take the inverse fields are calculated at t = n + 1/2t using the preFourier transform, we obtain a second-order differential viously calculated and stored E fields, and the process equation relating D and E continues until time-stepping is concluded. Applying this scheme to Eq. (11) we obtain D  2 D + 2 = 20 0 E (17) 20 D + 0 t t t n+1/2 n−1/2 n H = H −  × E (13) 0 Equation (17) is discretized using a second-order accurate central-difference scheme similar to those described above. We observe that the leapfrog scheme yields centralWe note that, if the ADE method is used, E is obtained difference in time and therefore second-order accurate from H in two steps. First, D is obtained from H by solvapproximations. In addition, since E(H) fields are obtained ing the finite-difference approximation of from previously calculated and stored H(E) fields, the time-stepping is fully explicit, meaning that we do not D  ×H = (18) have to solve a system of simultaneous equations.28 t To achieve second-order accuracy in space, FDTD Second, E is obtained from D by solving the finiteuses a special grid, known as the Yee lattice, where difference approximation of Eq. (17). Calculation of the every E component is surrounded by four H compofinite-difference expressions of the first and second time nents and every H component is surrounded by four E 28 derivatives of D in Eq. (17) requires the storage of 2 precomponents. Based on this arrangement, discretization of vious values of D. In other words not only Dn+1 but also the x-component of Eq. (11) gives n n−1 D and D are required to obtain E from D. Hx n+1/2 Another approach to model dispersive materials with i
j
k  FDTD is the recursive convolution (RC) method.28
36 t Ey ni
j
k+1/2 − Ey ni
j
k−1/2 n−1/2 = Hx i
j
k + FDTD is a finite-difference method, so its performance 0 z in modeling plasmonic devices is similar to the perfor Ez ni
j+1/2
k − Ez ni
j−1/2
k mance of FDFD. However, there are some major differ− (14) ences. First, as already mentioned above, in time-domain y  × E = −0 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 1813 Modeling of Plasmonic Waveguide Components and Networks methods the dispersion properties of metals have to be approximated by suitable analytical expressions which introduce substantial error in broadband calculations. In addition, the implementation of the ADE or RC methods requires additional computational cost and extra memory storage.28
36 On the other hand, in FDTD it is possible to obtain the entire frequency response with a single simulation by exciting a broadband pulse and calculating the Fourier transform of both the excitation and the response.28 Veronis et al. εm x z a εi a′ x=0 (a) x=a+h = a′ + h′ x 4. MODE-MATCHING METHOD h εm z h′ In this section, we introduce the mode-matching technique for MIM waveguide devices. The mode-matching technique37 is commonly used in the microwave and optia cal domains.38–41 We derive the complete set of orthogonal εi a′ modes that the MIM waveguide supports and use it to x=0 apply the mode-matching technique to the calculation of (b) scattering at the junction between two guides Delivered with differ- by Ingenta to: 3. (a) Geometry for the even modes of the MIM waveguide. The University ofFig. Houston ent cross sections. 4.1. Spectrum of MIM Waveguides x = 0 plane contains a fictitious perfect electric conductor (PEC) to simIP : 129.7.158.43 plify the problem when dealing only with even TM modes of the guide. Mon, 05 Oct 2009 14:12:14 This fictitious MIM waveguide is equivalent to an actual guide with an The mode-matching method is based on expanding the fields in terms of the modes of the waveguides. In this context here we derive the modal structure (spectrum) of the MIM waveguide.42 We will specifically focus on the even modes of the waveguide, for which the transverse magnetic (TM) field component is an even function of the transverse coordinate, x. The reason why we focus on even modes is that we will be analyzing the scattering of the fundamental, even mode of the MIM waveguide—which is also a TM mode10
11
15 —off of a symmetric junction with a different sized MIM waveguide. Due to the symmetry of the problem at hand, even modes will be sufficient. Evenness of the function is achieved by putting a fictitious perfect electric conductor (PEC) at the x = 0 plane of the waveguide, which forces the tangential electric field Ez to be an odd function, and the magnetic field Hy to be an even function of x. In other words, the modes of this fictitious waveguide with the PEC at x = 0 are mathematically the same as the even modes of the actual waveguide of interest, and so we will work with this hypothetical waveguide. The geometry is as shown in Figure 3(a). m refers to the permittivity of the metal region and i of the insulator region. At infrared frequencies, m (Fig. 1) is a complex number with a large, negative real part and a relatively small imaginary part (the sign of which is determined by the time convention used, being negative for an expit time dependence). We begin with Maxwell’s equations (Eqs. (1), (2)). The MIM waveguide is a two dimensional structure which does not have any variation in the y direction. Therefore, we can eliminate all the derivatives with respect to y in Maxwell’s equations. Furthermore, our study will be based on the TM 1814 insulator thickness of 2a. The inset shows the equivalent symmetric junction of two MIM waveguides. Dashed line in the inset is the plane of symmetry, which is where the fictitious PEC layer is introduced. (b) Geometry for mode matching. x = a + h plane of Figure 3(a) is terminated by a PEC which leads to a discretization of the continuous spectrum. modes which only have the Hy , Ex and Ez field components. Also, the uniformity of the waveguide in the z direction leads to exp−ikz z as the space dependence in z by using the separation of variables technique for differential equations (kz may, however, be a complex number). After simplifying the curl equations in (1), (2), we have the following relationships between the different field components d E x dx z ikz Hy x = i xEx x ixHy x = ikz Ex x + (19) d H x = i xEz x dx y Using these equations we get the following differential equation for Hy   d 1 d + 2 x x Hy = kz2 Hy (20) x dx x dx and since Ez 0 = 0 by the PEC wall at x = 0, the boundary condition for Hy under (19) becomes dHy x/dxx=0 = 0. 4.1.1. Point Spectrum The dispersion equation that should be solved in order to find the kz values for the modes of the MIM waveguide is J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. derived by satisfying the continuity of tangential electric and magnetic fields at material boundaries and applying the boundary conditions. We refer the reader to Refs. [43 and 9, 11, 44–46] for the details. The eigenvectors &n  and the dispersion equation for the corresponding eigen2  of (20) for the even TM modes of the MIM values kz
n waveguide are    cosh'i
n x 0 < x < a  &n x = H0 cosh'i
n a (21)   e−'m
n x−a a<x< tanh'i
n a = − 'm
n / 'i
n / 2 2 2 kz
n = 'm
n +   m m (22) Re'2m
,  < 0 and Im'2m
,  = 0. These conditions can be written in terms of kz
, by using (26) as 2 2 Rekz
, −  m <0 and 2 2 Imkz
, −  m =0 Note that when (22) holds true, we have - = −1 in (25) which makes (24) and (21) equivalent. 4.2. Mode Orthogonality Orthogonality and completeness are two very valuable properties of modes, which make the mode matching technique possible. Here we use the pseudo-inner product, ··, defined as50 i = '2i
n + 2  i f g = (23) f xgxdx 0 It can be shown that two different eigenfunctions, &1 x where Re'm
n  > 0 so that &n x does not diverge and is 2 and &2 x, corresponding to two different eigenvalues kz
1 integrable. Here n is a discrete index for the eigenvalues −1 50
51 2 and kz
2 are Delivered by Ingenta to: pseudo-orthogonal with x weight and the eigenfunctions. University of Houston  −1 &1  &2  = 0 (27) IP : 129.7.158.43 Mon, 05 Oct 2009 14:12:14 From (19) it can be seen that −1 &1 is proportional to 4.1.2. Continuous Spectrum In this section, we show that a continuous spectrum exists in the MIM waveguide. The utility of the continuous spectrum will be evident in the mode matching analysis. The condition of square integrability of the modes can be replaced by the weaker condition of finiteness of the modes in their domain of definition.47 For the MIM waveguide, this would imply a non-zero, yet finite electromagnetic field at infinity. These infinite-extent and, therefore, infinite energy, continuum modes (which can be normalized through the use of the Dirac delta distributions48
49 ) are integrated to realize any physically possible finite energy field configuration. In this respect, such an approach is similar to the well-known Fourier transform methods, where finite energy functions are expanded in terms of the infinite energy exponentials. Constraining fields to be finite, instead of zero, at infinity leads to the following field profile  cosh'i
, x   0<x<a   cosh'  i
, a  +, x = H0 (24) cosh'm
, x − a       + - sinh'm
, x − a a < x < -= 'i
, / 'm
, / i tanh'i
, a (25) m 2 kz
, = '2m
, + 2  m = '2i
, + 2  i (26) which is calculated very similarly to the case of a dielectric slab.43 Here , is a continuous index for different functions in the continuous spectrum. For finite +, , the arguments inside the hyperbolic functions for x > a in (24), 'm
, , should be purely imaginary which implies that J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 the transverse electric field component Ex of the mode. Therefore, the orthogonality condition can also be written as 0 Ex1 xHy2 xdx = A E1 r × H2 r · dA = 0 which is the well known modal orthogonality condition proved by the Lorentz reciprocity theorem,52 where A denotes the cross section of the waveguide. One can directly verify (27) by integration and using '2m
1 − '2m
2 = '2i
1 − '2i
2 which is a result of (23). The following orthogonality conditions between the elements of the point (&n ) and the continuous (+, ) spectrum can similarly be proved  −1  −1 &n +,  = 0 for all n and , + +,  = 0 for ,= In the following sections, we will be working with fields at the junction of two different waveguides. For notational abbreviation we will use the following convention i /L
R2 i /L
R2 e/L
R2 = Ex
i h/L
R2 = Hy
i where /L
R2 is used to denote the modes of the left and right side of the junction, which leads to the following orthogonality condition i j e/L
R2  h/L
R2  = 3ij 4/L
R2 (28) where 3ij is the Kronecker delta function and 4 is the overlap integral of the electric and magnetic transverse fields. 1815 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. Table I. Adjectives. fields have an exp−ikz z dependence. The argument principle method gives us the 'm value for the modes. By Signifier Signified using (23) we get the kz2 value. We then calculate kz2 1/2 Leaky Re'm  < 0 and choose the root which satisfies Imkz  < 0. Proper Re'm  > 0 In Figure 4 the spectrum of an idealized lossless silverImproper Re'm  = 0 like MIM waveguide is shown on the plane of '2m for Forward Rekz  > 0 Backward Rekz  < 0 m = −143497 which is the real part of the permittivity of silver at a wavelength  of 1550 nm.27
59 There are four real modes for 2a = /4—TM0 , TM2 , TM4 , TM6 — After the classification and analysis of the MIM wavegindexed according to the number of zero crossings in Hy . uide modes, we will now visualize different parts of its There is also an infinite number of complex modes, which spectrum by finding the zeros of the respective disperare those with eight and more zero crossings in the insusion equations through the use of the argument principle lator region. These modes have a 'm with a positive real method.53–58 We will use the adjectives in Table I to further part that is rather small compared to the imaginary part— differentiate between the modes. this can also be deduced from the scale of the imaginary Leaky modes are not normalizable and are not part of the axis of Figure 4. The continuous spectrum is illustrated spectrum. Proper modes can be normalized by the usual by the thick line which corresponds to Re'2m  < 0 and integration and they form the point spectrum. Improper Im'2m  = 0. This line is also the branch cut of the square modes can be normalized by using the Dirac delta funcroot function that is used to get 'm from '2m . The field proDelivered by Ingenta to: tions, 3x. They form the continuous spectrum. Forward files of the modes in the insulator region, as shown in the University of Houston modes have a positive phase velocity, whereas the backinsets of Figure 4, look quite similar to the field profile of : 129.7.158.43 ward modes have a negative phase velocity. WeIPdecide the even modes of a parallel plate waveguide with a plate 05 Oct 14:12:14 on the sign of Rekz  based on Imkz : By Mon, definition, all 2009 separation of 2a. 2 modes are propagating in the +z direction. Therefore, in 'm = 0 is the bifurcation point for the point specthe limit z → , the fields should go to zero. Such a trum when m is purely real. For positive '2m , the point spectrum has real modes, whereas for negative '2m , the behavior is possible only if Imkz  is negative, since the TM8 TM6 0.4 TMC1 2 Im(km )/k02 0.2 TM8 TMC2 TM4 TMC1 TM2 0.0 TM6TM4TM0 TMC2 TM2 –0.2 –0.4 –800 TM0 –600 –400 –200 0 200 2 Re(km )/k02 Fig. 4. Spectrum of the MIM waveguide for m = −143497 and 2a = /4 where  = 1550 nm is the wavelength of operation. There are four real modes and an infinite number of complex modes, all denoted with the • symbol. The thick line denotes the continuous spectrum. Due to the fact that m is real, complex modes come in complex conjugate pairs. Insets show the Hy mode shapes in the x direction for the discrete spectrum (TM0 through TM8 ) and the continuous spectrum (TMC1 and TMC2 )—solid lines in the insets are the real part of the mode, dashed lines are the imaginary part. The locations of the drawn continuous modes are shown by the • symbol. Modes in the continuous spectrum are purely oscillatory in the x direction. Complex modes have a small decay, which is not visually apparent in the inset for TM8 . 1816 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Veronis et al. point spectrum splits into two branches that are complex conjugates of one another. '2m = 0 corresponds to kz2 = m k02 which then implies kz = −i  m k0 —bounded modes should have Imkz  < 0. When loss is introduced to the metal, the spectrum moves on the complex plane.42 The forward, proper, complex modes of the lossless case turn into leaky modes by migrating into the third quadrant of the complex 'm plane. Modeling of Plasmonic Waveguide Components and Networks waveguides convergence of the fields on both sides of the junction is obtained only when the continuous spectrum is also taken into consideration. Thus, it is clear that the point spectrum on its own is not sufficient to describe the behavior of the waveguide junctions. Inclusion of the continuous spectrum is essential. This is further illustrated in the next section. 4.5. Convergence of the Mode-Matching Method 4.3. Discretization of the Continuous Spectrum Now that we know how to treat the continuous spectrum and are confident that the collection of the point and the continuous spectrum results in a complete basis set, we can proceed with the mode-matching formalism. We will begin by assuming that the pth mode of the left waveguide propagates toward the right, scatters and creates the following set of fields at the right and left sides of the junction, which by the continuity of the tangential magnetic and electric fields, are set equal The presence of a continuous spectrum leads to the formation of integral equations when the mode-matching method is applied.48 The integral equation is then expanded using an orthogonal basis set—not necessarily that of the modes—to solve the scattering problem. Another way to approach the scattering problem is to limit the transverse coordinates by a PEC wall. This approach has the effect of discretizing the continuum part Delivered by Ingenta to: of the spectrum47
60 —turning it into a discrete spectrum. University of Houston

m k To limit parasitic reflections from the PEC walls, absorb3mp + Rmp hL x = Tkp hR x (31) IP : 129.7.158.43 38 ing layers can be positioned before the PEC termination. m=1 k=1 Mon, 05 Oct 2009 14:12:14 In Ref. 61 detailed analysis of how the continuous spec

m k 3mp − Rmp eL x = Tkp eR x (32) trum appears from a discrete collection can be found. We m=1 k=1 will use a PEC wall to discretize the continuous spectrum. Also, we will not use any perfectly matched layers to limit Here Rmp is the reflection coefficient of the mth mode of parasitic reflections since the metallic sections with perthe left waveguide in response to an incoming field in the mittivity m effectively absorb the fields away from the pth mode. Similarly, Tkp is the transmission coefficient of junction. the kth mode of the right waveguide. Note that we chose The geometry is as shown in Figure 3(b). For the left Rmp to denote the reflection coefficient for the transverse waveguide the dispersion equation for modes becomes magnetic fields, which automatically results in −Rmp as the reflection coefficient for the transverse electric fields. 'm
n / m tanh'i
n a = − tanh'm
n h (29) In Ref. [62], it is shown that the testing of the above 'i
n / i equations (i.e., the discretization of the equation using integration of both sides by a given function) should be done which asymptotes to (22) as h → . The transverse magby the magnetic field of the larger waveguide for enforcing netic field shape is electric field continuity (32) and by the electric field of the  smaller waveguide to enforce the magnetic field continu cosh'i
n x  0<x<a   cosh'i
n a ity (31). Although that analysis was specifically done for &n x = H0 (30) waveguides with perfect metals  m  → , we still use  cosh'm
n x − a − h    that strategy so that the formulation limits to the correct a < x < a+h cosh'm
n h one should metals be made perfect. For those cases where a < a′ , we will take the pseudon n inner product of (31) with eL and of (32) with hR . Fur4.4. Mode Completeness thermore, assuming there are L modes on the left and R modes on the right, we get We tested the completeness of the modes by expanding the R L fundamental mode of an MIM waveguide of a given thick

n k m Tkp eL  hR  3mp + Rmp 4L 3mn = ness in terms of the modes of the MIM waveguide with a m=1 k=1 different thickness. Without the continuous spectrum, the R L

field expansion converges, but to a field profile which is k m n  = Tkp 4R 3kn 3 − R e  h 42 mp mp R L not the same as the desired profile of the right junction. m=1 k=1 On the other hand, inclusion of the continuous spectrum through the discretization of the continuum by a PEC wall with the help of (28). These are linear matrix equations leads to the correct field profile. Similarly, we found that with Rmp and Tkp as the unknowns. After calculating the inner products, the set of equations can be inverted to for the magnetic field profile at the junction of two MIM J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 1817 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. the mode-matching calculations to converge. Without the continuous spectrum, the mode matching results converge to the wrong result. Inclusion of the continuous spectrum 2a′ 2a decreases the error to around 2%, which is probably due 0.3 to the space discretization of FDFD simulations as well as the method used in the de-embedding of the scattering 2a = 0.9λ 0.2 coefficients from fields. As is also evident from Figure 5(a) 2a = 0.5λ the utility of the single mode L = R = 1 mode-matching 0.1 calculations increases as the dimensions of the waveg2a = 0.1λ uides decreases. The single mode approximation is closely 0.0 related to the characteristic impedance model described in 5 10 15 20 25 30 Number of modes Section 5.1 which is a good approximation for deep subwavelength structures. In Figure 5(a) we also show the 90º (b) 120º effect of neglecting the backward modes in the mode60º matching calculations for the 2a = 01 case. Backward modes are important in this subwavelength geometry; however, for the wider geometries of the 2a = 05 and 150º 30º 2a = 09 cases we did not observe any increase in the error when Delivered by Ingenta to:backward modes were neglected in the modematching University of Houstoncalculations. 1 In Figure 5(b) we visualize the scattering coefficient IP : 0º129.7.158.43 180º S22 of main mode of the MIM waveguide. (The extracMon, 05 Oct 2009the 14:12:14 tion of the scattering parameters from full field simula0.2 1 S11 tions is described in detail in the section below.) We do 0.4 the calculations in two different ways, one using FDFD, 330º 210º and the other using the mode-matching technique with 0.6 the point and the continuous spectrum. When applying 0.8 mode-matching, we use the a > a′ formulation for S11 300º 240º calculations and the a < a′ one for S22 . There is a very 1.0 good match between the results of the two techniques, 270º verifying the applicability of the mode-matching method. Fig. 5. (a) Convergence study of the reflection coefficient, S11 , of the In addition, since the convergence of the mode-matching main mode of the left waveguide traveling toward the right waveguide for method requires relatively few modes, the use of the mode2a = 09, 2a = 05 and 2a = 01. a′ /a = 04 and m = −143497 − i9517 for all cases. Dashed lines are for calculations including the point matching technique results in orders of magnitude speedup spectrum only. Solid lines are those with both the point and the conin the calculation of the properties of the junction. (a) Error 0.4 tinuous spectrum. Empty circles denote the calculations done with the forward point spectrum and the continuous spectrum for the 2a = 01 MM FDFD FDFD −S11 /S11  where MM stands for case only. Error is defined as S11 mode-matching and FDFD for finite-difference frequency-domain calculations. The inset shows the junction geometry. (b) Reflection coefficient of the main mode at the junction between two MIM waveguides of different insulator thicknesses, plotted on the complex plane within the unit circle. Both waveguides have m = −143497−i9517 and i = 10. Filled circles, •, are FDFD results, empty circles, , are mode-matching results. The left waveguide thickness is fixed at 2a = 09. The right waveguide thickness varies from 2a′ = /002
004
  
092. The origin is the zero reflection point that corresponds to a = a′ . As a′ decreases progressively toward zero, we move progressively along the curves away from the origin. The first set of curves, S11 , are for the case when the mode of the left waveguide, traveling from left to right, is scattered by the junction. The second set of curves, S22 , are for the case when the main mode of the right waveguide, traveling from right to left, is scattered by the junction. Insets illustrate the respective cases. give the reflection and transmission coefficients for the modes. In Figure 5(a), we compare the mode-matching method with the FDFD technique. It takes relatively few modes for 1818 5. EQUIVALENT MODELS FOR PLASMONIC WAVEGUIDE COMPONENTS Modeling electromagnetic wave propagation using transmission lines has been one of the most important achievements of microwave network theory.63 The concept of impedance64 and understanding the effects of waveguide discontinuities in terms of lumped circuit elements were crucial in this respect. Even though the properties of metals are quite different at optical wavelengths compared to the microwave, designs that are qualitatively similar to their low frequency counterparts have been demonstrated at optical frequencies.65 It has also been shown that the fundamental TM mode of an MIM waveguide continuously changes to the TEM mode of a parallel-plate waveguide with PEC boundaries as the frequency of operation is decreased.66 It is intriguing to ask whether methods of microwave can be applied to this new generation of nanometallic structures to come up with concise descriptions of J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks components that can lead to a simplified approach to the design of functional systems composed of many interacting parts. Here we first introduce the concept of the characteristic impedance for MIM plasmonic waveguides and we show its validity and usefulness for subwavelength guides. We then characterize the modal reflection and transmission from MIM junctions using the scattering matrix approach, a commonly used method in microwave network theory. We also represent the scattering matrix of MIM junctions in terms of an equivalent lumped circuit model. Finally, to test our characterization, we design a cascade connection of MIM junctions to couple the mode of a wavelengthsized MIM waveguide to that of a subwavelength one with zero reflection. Throughout our analysis, we will compare MIM waveguides to PEC parallel plate waveguides and comment on the similarities and the differences between the two. 0.12 0.1 0.08 Reflection Veronis et al. din 0.06 0.04 dout 0.02 0 1 1.5 2 2.5 3 din/dout Fig. 6. Reflection coefficient R of a MIM T-shaped splitter (shown in the inset) as a function of din /dout at 0 = 155 m calculated using FDFD. We also show with dashed line the reflection coefficient R̄ calculated based on the characteristic impedance ZMIM and transmission-line theory. Results are shown for dout = 50 nm. Delivered by Ingenta to: University ofcompared Houstonto the wavelength and the quasistatic approxIP : 129.7.158.43 imation therefore breaks down. We found that similar The characteristic impedance of the fundamental Mon, 05TEM Oct 2009 14:12:14 deviations are observed for PEC parallel-plate waveguides. 5.1. Characteristic Impedance Model mode in a PEC parallel-plate waveguide is uniquely defined as the ratio of voltage V to surface current density I and is equal to67 ZTEM ≡ V E d : = x = TEM d = I Hy  0 0 d (33) 0 where Ex , Hy are the transverse components of the electric and magnetic field respectively, and we assumed a unitlength waveguide in the y direction. For non-TEM modes, such as the fundamental MIM mode, voltage and current are not uniquely defined. However, metals like silver satisfy the condition  metal  ≫ diel at the optical communication wavelength of 1.55 m.27 Thus, Ex metal  ≪ Ex diel  so that the integral of the electric field in the transverse direction can be approximated by Ex diel d and we may therefore define the characteristic impedance of the fundamental MIM mode as15 ZMIM d ≡ Ex diel d :MIM d = d Hy diel  0 (34) where :MIM d = 2;/g d, and g is the guide wavelength. In Figure 6 we show the reflection coefficient of an MIM T-shaped splitter calculated based on ZMIM as      Z − Z0 2  2ZMIM dout  − ZMIM din  2   = R̄ =  L (35)  2Z d  + Z d   ZL + Z0  MIM out MIM in We observe that there is very good agreement between R̄ and the exact reflection coefficient R calculated using FDFD. This agreement suggests that the concept of characteristic impedance for MIM waveguides is indeed valid and useful. The deviation between R̄ and R at large values of din /dout is due to the fact that din is not very small J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Such deviations decrease at longer wavelengths in both the PEC and MIM waveguide cases. 5.2. Scattering Matrix Description of Junctions In this section, we characterize the modal reflection and transmission from MIM junctions using the scattering matrix approach, a commonly used method in microwave network theory.68 We consider the geometry shown in Figure 7(a). The insulating region is free space with a permittivity i = 1, and the metal is silver with a permittivity of m = −143497 − i9517.27
59 The wavelength of operation is fixed at  = 1550 nm, in the L band of optical telecommunications. Using the dispersion equation for even modes of the MIM waveguide it can be shown that only a single even propagating mode can exist for 2a < 097 for our choice of m , i and . The condition for the PEC parallel plate waveguide is similar, where only a single even propagation mode exists for 2a < 10. When there is only one propagating mode, far away from the waveguide junction the fields can be written in terms of that main mode of the system since all higher order modes will have an exponential decay much faster compared to the main propagating mode. Under such circumstances, the effects of the waveguide junction on the propagating modes can be described using the single mode scattering matrix (S) formalism.69 In the terminology of the scattering matrix, the forward and backward mode amplitudes are considered to scatter from one “port" to another. Here we can think of the ports as being the left and right port planes shown in Figure 7(a). These ports are sufficiently far to the left and right of the junction that the fields have settled down again to being 1819 Modeling of Plasmonic Waveguide Components and Networks εm 5λ 5λ Veronis et al. 1 1 S11 x εi 2a 2a′ z (b) (c) 1 + HL Left port Right port 1 S22 R + HR S – (d) – HR HL (e) (a) Fig. 7. (a) Description of the modeling geometry. Dashed lines represent the location of the left and right ports of the overall scattering matrix S that describes this junction (schematically shown in the bottom part of the figure). (b) Description of the steps taken in extracting S from fields. Calibration simulations with uniform insulator widths of 2a and 2a′ , which give the wave vector k and the values of the incoming fields at the left, HL+ &L x, Delivered by Ingenta to: and right ports, HR+ &R x. (c) Field impinging from the left side, which leads to S11 . (d) Field impinging from the right side, which leads to S22 . University of Houston (e) Simulation domain is terminated by a perfect electric conductor at the right input port plane. S12 is extracted from the reflection coefficient R using IP : 129.7.158.43 the previously calculated S11 and S22 . Mon, 05 Oct 2009 14:12:14 the propagating modes of the guides (nearer to the interface, there will in general be other field behavior, including various near-field components that decay rapidly with distance). If we can deduce the scattering matrix for such a junction, then we can have a very simple way of modeling the behavior of structures containing such junctions, as is already well known in the modeling of microwave guides. The elements of the scattering matrix, S11 , S12 , S21 , S22 , are complex numbers which describe the phase and magnitude of the reflection and transmission of the main modes. Thus, in general there are 8 independent real numbers in S. However, under certain conditions the number of independent parameters can be reduced. First of all, if the system is composed of reciprocal media (i.e., symmetric permittivity and permeability tensors) then using the Lorentz reciprocity theorem it can be shown that S12 = S21 . Note that this equality implies a certain normalization of the modes,69 specifically A EL × HL · dA = A ER × HR · dA = 1 (36) where E/L
R2 and H/L
R2 denote the electric and magnetic components of the main propagating modes on the left (L) and the right (R) of the waveguide junction. Also note that for lossless systems S is a unitary matrix69 (though in general here we will be considering systems with loss). As a result, using reciprocity it is possible to describe a lossy junction using six real numbers, two for each of S11 , S12 , and S22 . When there is no loss we only need three real numbers due to the unitarity of S. We now describe the method we used to extract the elements of S from the electromagnetic fields in such waveguide junctions. We solved Maxwell’s equations using the 1820 FDFD method. In the vicinity of the waveguide junction, higher order modes will be excited. We chose the left and right ports of our junction sufficiently (5) away from the physical junction where the amplitudes of the higher order modes are negligible. In the following we will formulate S in terms of the transverse magnetic field component Hy . We will use &L x to denote the main mode of the left waveguide and &R x for the right waveguide. The scattering matrix relates the amplitudes and phases of the modes that arrive at the left and right ports, HL+ , HR+ to the amplitudes and phases of the modes that propagate away from the ports, HL− and HR− . Formally we can write  HL− HR−  = S 11 S21 S12  H +  L S22 HR+ (37) In order to extract S we need to know the fields that arrive at the left and right ports from our numerical sources in the simulation domain. To do that we do two calibration simulations (one for the left waveguide, another for the right waveguide) without any discontinuities, as shown in Figure 7(b), and record the fields. This gives us + the required H/L
R2 &/L
R2 x in addition to the propagation vectors, k/L
R2 , of the two main modes for guides of insulator thicknesses 2a and 2a′ respectively. Then we do two more simulations where we send the mode from the left and from the right waveguide to the discontinuity as schematically shown in Figures 7(c–e). From the results of the simulation in Figure 7(c), for the fields to the left of the left port, HL x
z, we get HL x
z = HL+ e−ikL z + HL− e+ikL z &L x = HL+ e−ikL z + S11 e+ikL z &L x J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Veronis et al. Modeling of Plasmonic Waveguide Components and Networks where the location of the left port determines the origin for z and in (37) we used the fact that HR+ = 0 for the simulation depicted in Figure 7(c). Simple algebra gives Taking these circuit and transmission line approaches together, we can then model a broad range of MIM systems with circuit models. There is no unique way to describe S using lumped cirH x
z −ikL z cuit elements.71 To choose one circuit out of the infinite − e−2ikL z (38) S11 = +L e HL &L x possible set that could correspond to the same S, we will first look at the well studied PEC case. The solution to Very similarly, we also extract S22 from the results of the the scattering problem for the junction of two PEC parsimulation of Figure 7(d). allel plate waveguides was developed and experimentally In order to extract S12 , we terminate our simulation verified.72 It consists of a capacitor with susceptance B and domain at the plane of the right port with a perfect electric a transformer with a turns ratio of n:1. The susceptance conductor. Such a termination results in zero tangential and the turns ratio are described in terms of the geometry electric fields and therefore gives −1 for the reflection of the junction. The square of the turns ratio of the transcoefficient of the transverse electric field, Ex , and +1 for former is equal to the magnetic field, Hy . Thus, at the right port we get HR− = a′ HR+ . Using this equality in (37), we get the reflection coefn2 = (41) a ficient, R = HL− /HL+ , in Figure 7(e) It is worthwhile remembering that the primary-secondary 2 S12 S12 S21 turns ratio of the transformer, n:1, is also the ratio of = S11 + (39) R = S11 + the voltages 1 − S22 1 − S22 Delivered by Ingenta to: at its terminals. From the conservation of power, currents have the inverse ratio and as a result Houston where in the last equality we used the fact thatUniversity S21 = S12 . ofthe impedance ratio at the transformer terminals is n2 :1. We extract R using the same method as we usedIPin: 129.7.158.43 the The derivation of the circuit elements for the PEC case Mon, 05 Oct 14:12:14 S11 and S22 2009 extraction of S11 . From the knowledge of R, can be found in Refs. [73–75]. Note that for the PEC one can easily invert (39) to calculate S12 . case, only two parameters, B and n, are sufficient to After we calculate S for the ports defined in Figure 7(a), describe the junction even though in general three paramewe shift both the left and right reference planes back to ters are required for a lossless reciprocal system. The nonthe exact location of the junction using dispersive nature of the main mode of PEC parallel plate    ik ℓ  ik ℓ waveguides leads to a further symmetry in the junction L L L L 0 e 0 e (40) SJ = ikR ℓR ikR ℓR S which reduces the number of circuit parameters required. 0 e 0 e At optical frequencies where the modes are strongly diswhere ℓL = ℓR = 5 as defined in Figure 7(a) and SJ is persive, a third circuit element is needed in order to be the effective scattering matrix for the case where the left able to fit the elements of S exactly. For that reason we and right ports are projected back to coincide with the have an inductor term with a reactance X. A schematic of junction plane.69 For the sake of notational abbreviation, the circuit diagram is shown in Figure 8(a). The PEC parfrom this point on we will use S to imply SJ . Note that allel plate waveguide circuit is the same, with X = 0. The this effective scattering matrix is defined just for the algenormalization that we defined in (36) leads to transmission braic convenience of having a scattering matrix associated lines with a unit characteristic impedance on both sides directly with the position of the interface. The fields near of the junction. From transmission line theory we get the to the interface are not in fact just describable by these following equalities in terms of the equivalent impedance single main modes because of various near field effects of looking from the left side of the circuit, Z1 , and the equivhigher order modes. alent admittance looking from the right side, Y2 Note that the polar plot of the reflection coefficients Z −1 1 − Y2 is nothing other than the Smith chart of the microwave −S11 = 1 and − S22 = (42) 70 Z + 1 1 + Y2 theory, which we will use in Section 5.4. 1 where 5.3. Circuit Model for the Waveguide Junction Another important approach in microwave waveguide modeling is the use of equivalent circuit models, which can give an intuitive picture of the system as well as allowing the use of circuit simulators for design. Here we relate the scattering matrix to simplified circuit models that can characterize the MIM waveguide interfaces. Since we only have single propagating modes in the guides we consider, we can also use equivalent transmission lines to describe the propagation between interfaces. J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 1 iB + 1/n2    1 + iB n2 Y2 = 1 + iX Z1 = iX + (43) The reason why we have negative signs in front of S11 and S22 in (42) is because we defined S based on the transverse magnetic component of the main mode, Hy . However, the norm in circuit parametrization is to use the voltage reflection and transmission coefficients, which correspond to 1821 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. 1 0.8 iX Y2 0.6 n2 iB 0.4 Z1 n:1 0.2 (a) (b) 0 0 0.2 0.4 0.6 0.8 1 a′/a 0.07 2 (c) (d) 0.06 1.6 0.05 1.2 X B 0.04 Delivered by Ingenta to: University of Houston 0.03 IP : 129.7.158.43 Mon, 05 Oct 20090.0214:12:14 0.8 0.4 0.01 0 0 0 0.2 0.4 0.6 0.8 1 a′/a 0 0.2 0.4 0.6 0.8 1 a′/a Fig. 8. (a) Circuit description of the lossless ( m = −143497) waveguide junction for 2a = 01 (
), 2a = 05 (), 2a = 09 (•). (b) Square of the turns ratio, n2 , which is equivalent to the impedance ratio at the terminals of the transformer. Dashed line is the PEC result described by (41). (c) Susceptance, B, for the MIM (
, , •) and the PEC (△, , ) case. (d) Reactance, X, of the MIM waveguide. X = 0 for the PEC parallel plate waveguide. a scattering matrix description for the transverse electric component, Ex . Just as in transmission line theory where the reflection coefficient for voltage is the negative of that of the current, the same relationship also holds exactly between the reflection coefficients of Ex and Hy . Let the real and imaginary parts of Z1 and Y2 be denoted as ZR = ReZ1 , ZI = ImZ1 , YR = ReY2  and YI = ImY2 . Using (43) we get Bn2 = Y I + YR ZI 1 − YR ZR and X = Z I + Z R YI 1 − YR ZR (44) Once we know Bn2 and X, we can calculate n2 using (43) as   n2 = YR 1 + X 2  = ZR 1 + Bn2 2 (45) Using Eqs. (42), (44) and (45) one can calculate the circuit parameters from S11 and S22 . In Figures 8(b–d) we plotted n2 , B and X as a function of a′ /a for the three different fixed 2a values of 01, 05 and 09. It can be seen that the PEC circuit description and the MIM circuit description lead to parameters which qualitatively have similar behaviors. 1822 5.4. Cascade Connection of Junctions We now test the utility of the scattering matrix description by numerically simulating mode propagation through a cascade connection of junctions and comparing the results with the predictions of the scattering matrix formalism. When different scattering matrices are cascaded, the transfer matrix, T, leads to a much simpler formulation76  +   −  HL T11 T12 HR (46) = HR+ T21 T22 HL− 5.4.1. Conditions for Zero Reflection In order to have HL− = 0, one should have T21 HR− + T22 HR+ = 0, which can be cast in terms of the scattering parameters as S11 HR− = S11 S22 − S12 S21 HR+ (47) a. Lossy Case. Let us investigate the case when two junctions characterized by two different scattering matrices, L S and R S, are separated by a center waveguide of length ℓ as shown in Figure 9(a). Suppose that we adjust J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. – + HC HL (a) RS 1 ℓ LS LS LS 11 21 LS RS 12 RS 22 – + HL RS HC 11 21 RS RS 21 12 22 11 90º (b) (c) 0.7 PL 0.6 22 PR 180º 0º RS 0.2 0.4 0.6 0.8 1. 270º 11 Power reflection LS 0.5 0.4 0.3 0.2 Delivered by Ingenta to: 0.1 University of Houston IP : 129.7.158.43 0 ℓ 0 0.2 300Mon, 05 Oct 2009 14:12:14 0.9λ 0.5λ w 0.4 0.6 0.8 1 1.2 Center waveguide length (ℓ/λ) Fig. 9. (a) Schematic diagram of modal propagation. The left and right junctions are shown as boxes with an S matrix description. The center waveguide is shown as a transmission line of length ℓ. The source that creates the fields is normalized such that the mode that propagates to the right has a unit magnitude at the input of the right junction. (b) Graphical solution of (50) and (51) on the complex plane. Point PL is the location of the left junction on the L S22 curve where 2a = 09 and 2a′ = 05. Point PR is the location of the right junction on the R S11 curve where 2a = 05 and 2a′ = 016. (c) Test of the scattering matrix description. Horizontal axis is the length of the center waveguide normalized to . Vertical axis is the power reflection coefficient. FDFD simulation results •, transfer matrix calculations using lossy junctions (solid line) and lossless junctions (dashed line) are also plotted. Transfer matrix calculations do take into account the loss in the center waveguide for both cases. As the junctions get very close to each other <01 transfer matrix model begins to break down due to higher order modal interactions. our excitation amplitude such that the mode that propagates toward the right junction at its input plane, which is the junction plane, has unit strength. That choice of normalization leads to HR− = eikC ℓ and HR+ = R S11 e−ikC ℓ where kC is the wave vector of the center waveguide. With these definitions, the condition for zero reflection, (47), for the left junction can be written as L e−2ikC ℓ = LS 11 S11 /R S11 L L 22 − S12 S21 LS For reciprocal media, S12 = S21 , we can write L e−2ikC ℓ = S11 /R S11 LS LS − LS2 11 22 21 (48) b. Lossless Case. If the system is lossless, then the scattering matrix should be unitary (SS† = 1) which implies the following three conditions S11 2 = S22 2 = 1 − S12 2 S S12 = − 22 ∗ ∗ S21 S11 (49) Using (48) and (49), the zero reflection condition becomes L S22  = R S11  (magnitude condition) J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 (50) ∡L S22 + ∡R S11 = 2kC ℓ + 2;n (phase condition) (51) where n is any integer value, superscripts R and L denote right and left respectively. “∡” is used to represent the argument of a complex number. What this means is that, to match a left waveguide to a right waveguide, one should choose a center waveguide width which satisfies the magnitude condition, and decide on the length of the center waveguide based on the phase condition. Once a matching left, center and right waveguide triplet is found, the procedure can be recursively repeated to cascade more junctions without getting any reflection at the leftmost waveguide. 5.4.2. Mode Converter Design Now that we have the conditions (48) and (50–51) for zero reflection, we can test their validity. Condition (48) is more general and is applicable to the lossy case. We did a series of simulations in which we extracted S for a hypothetical lossless metal with a real, negative permittivity m = −143497. The results were very similar to the case 1823 Modeling of Plasmonic Waveguide Components and Networks Veronis et al. where the loss was included. That led us to suspect that 6. SUMMARY AND CONCLUSIONS the conditions for the lossless reciprocal junctions, (50) In this paper, we first gave an overview of the challenges and (51), would be essentially sufficient in the design of involved in modeling of plasmonic devices, and briefly a mode converter that converts the mode of a wavelength examined two general purpose simulation techniques which sized MIM waveguide 2a = 09 to that of a subwaveare widely used for modeling of plasmonic waveguide length one with no reflection. devices. In our design we choose the left waveguide width to be We then introduced the mode-matching technique for 09 and the center waveguide width to be 05 as shown modeling of MIM waveguide devices. We derived the in the inset in Figure 9(b). The parameters that we need complete set of orthogonal modes that the MIM wavegare the insulator width of the right waveguide, w, and the uide supports. We showed that in addition to the previlength of the center waveguide, ℓ. ously unreported complex set of discrete modes, there also The width of the right waveguide can be chosen by satexists a continuous set of modes. We then used the comisfying (50). In Figure 9(b) PL is the location of the 2a = plete set of modes to apply the mode-matching technique 09 to 2a′ = 05 junction on the L S22 curve. To satisfy in the investigation of modal scattering at the symmet(50) we need to have L S22  = R S11 . The solution can be ric junction of two MIM waveguides with different cross graphically found by drawing a circle in the complex plane sections. with a radius PL  centered at the origin and finding its This special purpose simulation technique is far more intersection with the R S11 curve. The intersection point is efficient Delivered by Ingenta for to: this class of problems than general purpose denoted by PR . PR corresponds to a right waveguide thickelectromagnetic simulation techniques. The knowledge of University of Houston ness of 016. The phase condition (51) is then easily calthe set of orthogonal modes which form a complete basis IP : 129.7.158.43 culated from the phases of the scattering coefficients, ∡P for a given geometry leads to a much more simplified L Mon, 05 Oct 2009 14:12:14 algebra and speeds up calculations. These results are valuand ∡PR . After some simplification through the use of the able for electromagnetic scattering calculations involving numerical value for kC one gets ℓ/ = 01377 + 04861n, the MIM geometry. The analysis can also be generalized where n is any positive integer. for other related geometries involving metals at optical To check our design, we numerically simulated the frequencies.77–81 structure shown in the inset of Figure 9(b) using FDFD and We also introduced several different equivalent modlooked at the amount of power reflected back as a function els for plasmonic waveguide components. We first showed of the center waveguide length ℓ. We also calculated the that the characteristic impedance model can accurately power reflection coefficient through the use of the transmodel the behavior of deep subwavelength MIM wavegfer matrix formalism in which we multiplied the transfer uides. We also extracted the scattering matrices of juncmatrices for the right junction, TR , a center waveguide of tions of different geometries from full field solutions, and length ℓ, TC , and the left junction, TL , to get the overparametrized the scattering matrix of the MIM junction in 2 all transfer matrix T = TL TC TR , and plotted T21 /T11  of terms of lumped circuit elements. We validated our charT as a function of ℓ/. We did the calculations for two acterization by designing a mode converter that concendifferent sets of /TR
TL 2: one in which we used the scattrates light from an MIM waveguide of wavelength-sized tering matrices for the lossy junctions and another for the dimension to one of subwavelength dimension with zero lossless junctions. The center waveguide of length ℓ had reflection. loss in both cases i.e., kC = 103 − i945 × 10−4 2;/. The model abstraction provided by these equivalent Figure 9(c) verifies that lossless junction models are models is important for the analysis and synthesis of quite effective at modeling the waveguide discontinuities device functions.82–85 The circuit representation of the and the prediction of the length of the center guide for zero junction helps us associate the effects of geometry, matereflection reached by their use, ℓ/ = 01377 + 04861n rial properties and wave propagation in terms of a simis very accurate. The lossy junction model on the other ple network of a capacitor, inductor and a transformer. hand gives results essentially indistinguishable from the The scattering matrix description of junctions can be simulation results as long as the two junctions are not used to design optical circuitry with complex functionvery close to each other (< 01). When the junctions get ality using tools of circuit analysis.86
87 It is conceivvery close, the coupling of higher order non-propagating able to build a library of junction geometries associated modes becomes important and the single mode modeling with their scattering matrices for different waveguides we employed in the construction of scattering matrices including three dimensional nano-metallic ones.88 Such a breaks down. For such closely spaced junctions, the whole library, indexed according to modal scattering and propstructure should be treated as a single unit and its characagation properties, would be invaluable in the design of teristics should be extracted by the techniques described integrated optical circuits composed of many interacting components. in Section 5.2. 1824 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 Veronis et al. Modeling of Plasmonic Waveguide Components and Networks References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. I. Breukelaar, R. Charbonneau, and P. Berini, J. Appl. Phys. 100, 043104 (2006). W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature 424, 824 41. R. F. Oulton, D. F. P. Pile, Y. Liu, and X. Zhang, Phys. Rev. B 76, (2003). 035408 (2007). E. Ozbay, Science 311, 189 (2006). 42. Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, Phys. Rev. J. Takahara, S. Yamagishi, H. Taki, A. Morimoto, and T. Kobayashi, B 79, 035120 (2009). Opt. Lett. 22, 475 (1997). 43. G. W. Hanson and A. B. Yakovlev, Operator Theory for ElectromagJ. C. Weeber, A. Dereux, C. Girard, J. R. Krenn, and J. P. Goudonnet, netics, Springer-Verlag, New York (2002). Phys. Rev. B 60, 9061 (1999). 44. B. Prade, J. Y. Vinet, and A. Mysyrowicz, Phys. Rev. B 44, 13556 J. R. Krenn, B. Lamprecht, H. Ditlbacher, G. Schider, M. Salerno, (1991). A. Leitner, and F. R. Aussenegg, Europhys. Lett. 60, 663 (2002). 45. E. Feigenbaum and M. Orenstein, Journal of Lightwave Technology M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Phys. Rev. B 25, 2547 (2007b). 62, R16356 (2000). 46. A. R. Zakharian, J. V. Moloney, and M. Mansuripur, Opt. Express S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E. 15, 183 (2007). Koel, and A. A. G. Requicha, Nat. Mater. 2, 229 (2003). 47. V. V. Shevchenko, Continuous Transitions in Open Waveguides, S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, and T. W. Ebbesen, Golem Press, Boulder, Colorado (1971). Phys. Rev. Lett. 95, 046802 (2005). 48. T. Rozzi and M. Mongiardo, Open Electromagnetic Waveguides, The E. N. Economou, Phys. Rev. 182, 539 (1969). Institution of Electrical Engineers, London (1997). R. Zia, M. D. Selker, P. B. Catrysse, and M. L. Brongersma, J. Opt. 49. D. A. B. Miller, Quantum Mechanics for Scientists and Engineers, Soc. Am. A 21, 2442 (2004). Cambridge University Press, Cambridge (2008). J. A. Dionne, L. A. Sweatlock, H. A. Atwater, and A. Polman, Phys. 50. M. Mrozowski, Guided Electromagnetic Waves, Research Studies Rev. B 73, 035407 (2006a). Press, UK (1997). E. Feigenbaum and M. Orenstein Journal of Lightwave Technology Delivered by Ingenta to: 51. W. C. Chew, Waves and Fields in Inhomogenous Media, Van Nos25, 2547 (2007a). University of Houston trand Reinhold, New York (1990). R. Gordon, Phys. Rev. B 73, 153405 (2006). IP : 129.7.158.43 52. R. E. Collin, Field Theory of Guided Waves, 2nd edn., WileyZ. Han and S. He, Opt. Commun. 278, 199 (2007). Mon, 05 Oct 2009 Interscience, 14:12:14New York (1991a). G. Veronis and S. Fan, Appl. Phys. Lett. 87, 131102 (2005). 53. E. Anemogiannis and E. N. Glytsis, Journal of Lightwave TechnolK. Tanaka and M. Tanaka, Appl. Phys. Lett. 82, 1158 (2003). ogy 10, 1344 (1992). G. Veronis and S. Fan, Opt. Express 15, 1211 (2007a). 54. A. Bakhtazad, H. Abiri, and R. Ghayour, Journal of Lightwave TechP. Ginzburg and M. Orenstein, Opt. Express 15, 6762 (2007). J. A. Dionne, H. J. Lezec, and H. A. Atwater, Nano Lett. 6, 1928 nology 15, 383 (1997). (2006b). 55. M. S. Kwon and S. Y. Shin, Opt. Commun. 233, 119 (2004). Y. Jiao, S. Fan, and D. Miller, IEEE Journal of Quantum Electronics 56. R. Rodriguez-Berral, F. Mesa, and F. Medina, IEEE Transactions on 42, 266 (2006). Microwave Theory and Techniques 53, 1613 (2005). W. Suh, Z. Wang, and S. Fan, IEEE Journal of Quantum Electronics 57. R. Rodriguez-Berral, F. Mesa, and F. Medina, International Journal 40, 1511 (2004). of RF and Microwave Computer-Aided Engineering 14, 73 (2004). J. Bravo-Abad, S. Fan, S. G. Johnson, J. D. Joannopoulos, and 58. R. E. Smith, S. N. Houde-Walter, and G. W. Forbes, IEEE Journal M. Soljacic, Journal of Lightwave Technology 25, 2539 (2007). of Quantum Electronics 28, 1520 (1992). S. Fan, M. F. Yanik, Z. Wang, S. Sandhu, and M. L. Povinelli, 59. H. J. Hagemann, W. Gudat, and C. Kunz, Journal of the Optical Journal of Lightwave Technology 24, 4493 (2006). Society of America 65, 742 (1975). L. Novotny, B. Hecht, and D. Pohl J. Appl. Phys. 81, 1798 (1997). 60. A. Zettl, Sturm-Liouville Theory, American Mathematical Society E. Prodan, P. Nordlander, and N. Halas, Chem. Phys. Lett. 368, 94 (2005). (2003). 61. L. P. Felsen and N. Marcuvitz, Radiation and Scattering of Waves, C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light IEEE Press, New York (1994). by Small Particles, Wiley Science, New York (1983). 62. G. V. Eleftheriades, A. S. Omar, L. P. B. Katehi, and G. M. Rebeiz, E. D. Palik, Handbook of Optical Constants of Solids, Academic, IEEE Transactions on Microwave Theory and Techniques 42, 1896 New York (1985). (1994). A. Taflove and S. C. Hagness, Computational Electrodynamics, 63. A. A. Oliner, IEEE Transactions on Microwave Theory and TechArtech House, Boston (2005). niques MTT/32, 1022 (1984). A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, Appl. 64. S. A. Schelkunoff, IEEE Transactions on Antennas and Propagation Opt. 37, 5271 (1998). AP/18, 309 (1970). M. Han, R. Dutton, and S. Fan, IEEE Microwave and Wireless Com65. P. J. Schuck, D. P. Fromm, A. Sundaramurthy, G. S. Kino, and W. E. ponents Letters 16, 119 (2006). Moerner, Phys. Rev. Lett. 94, 017402 (2005). A. Vial, A. Grimault, D. Macias, D. Barchiesi, and M. de la 66. T. Takano and J. Hamasaki, IEEE Journal of Quantum Electronics Chapelle, Phys. Rev. B 71, 85416 (2005). QE/8, 206 (1972). J. Jin, The Finite Element Method in Electromagnetics, John Wiley 67. S. Ramo, J. R. Whinnery, and T. V. Duzer, Fields and Waves in & Sons, New York (2002). Communication Electronics, Wiley, New York (1994). P. Berini, Phys. Rev. B 61, 10484 (2000). 68. Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, IEEE Journal G. Veronis, R. W. Dutton, and S. Fan, Opt. Lett. 29, 2288 (2004). of Selected Topics in Quantum Electronics 14, 1462 (2008b). J. Berenger, J. Comput. Phys. 114, 185 (1994). 69. D. M. Pozar, Microwave Engineering, Wiley, New York (1998). J. Young and R. Nelson, IEEE Antennas and Propagation Magazine 70. P. H. Smith, Electronic Applications of the Smith Chart in Waveg43, 61 (2001). uide, Circuit and Component Analysis, Robert E. Krieger Publishing P. J. B. Clarricoats and K. R. Slinn, Proceedings of the Institution Company, Malabar, FL (1983). of Electrical Engineers (1967), Vol. 114, p. 878. 71. E. L. Ginzton, Microwave Measurements, McGraw-Hill Book ComP. Bienstman, Ph.D. thesis, Ghent University, Belgium (2001). pany Inc. (1957); Representation and Measurement of Microwave I. Breukelaar and P. Berini, Journal of the Optical Society of America A 23, 1971 (2006). Circuits, pp. 313–345. J. Comput. Theor. Nanosci. 6, 1808–1826, 2009 1825 Modeling of Plasmonic Waveguide Components and Networks 72. N. Marcuvitz, Waveguide Handbook, Peter Peregrinus Ltd., London, UK (1986), Chap. 5, pp. 307–309. 73. J. Schwinger and D. S. Saxon, Discontinuities in Waveguides, Gordon and Breach Science Publishers Inc., New York (1968), Chap. 5, pp. 99–124. 74. L. Lewin, Advanced Theory of Waveguides, Iliffe & Sons Ltd. (1951), Chap. 5, pp. 98–106. 75. R. E. Collin, Field Theory of Guided Waves, 2nd edn., WileyInterscience (1991b), Chap. 8, pp. 581–588. 76. R. E. Collin, Foundations for Microwave Engineering, McGraw-Hill Book Company (1966). 77. S. Collin, F. Pardo, and J. L. Pelouard, Opt. Express 15, 4310 (2007). 78. F. M. Kong, K. Li, B. I. Wu, H. Huang, H. S. Chen, and J. A. Kong, Progress in Electromagnetics Research 76, 449 (2007a). 79. F. M. Kong, B. I. Wu, H. S. Chen, and J. A. Kong, Opt. Express 15, 12331 (2007b). Veronis et al. 80. F. Lopez-Tejeira, S. G. Rodrigo, L. Martin-Moreno, F. J. GarciaVidal, E. Devaux, T. W. Ebbesen, J. R. Krenn, I. P. Radko, S. I. Bozhevolnyi, M. U. Gonzalez, et al., Nat. Phys. 3, 324 (2007). 81. M. W. Vogel and D. K. Gramotnev, Phys. Lett. A 363, 507 (2007). 82. A. Alu and N. Engheta, Journal of the Optical Society of America B 23, 571 (2006). 83. N. Engheta, Science 317, 1698 (2007). 84. A. Hosseini, H. Nejati, and Y. Massoud, Opt. Express 15, 15280 (2007). 85. A. Hosseini, H. Nejati, and Y. Massoud, Opt. Express 16, 1475 (2008). 86. SPICE http://en.wikipedia.org/wiki/SPICE. 87. HSPICE® Signal Integrity User Guide Synopsys, z-2007.03 edn., (2007), ch. S-parameter Modeling Using the S-element. 88. G. Veronis and S. Fan, Journal of Lightwave Technology 25, 2511 (2007b). Received: 8 October 2008. Accepted: 22 October 2008. Delivered by Ingenta to: University of Houston IP : 129.7.158.43 Mon, 05 Oct 2009 14:12:14 1826 J. Comput. Theor. Nanosci. 6, 1808–1826, 2009