## Abstract

Photonic functions are programmed by designing the interference of local waves in inhomogeneous resonant guided wave networks composed of power-splitting elements arranged at the nodes of a nonuniform waveguide network. Using a compact, yet comprehensive, scattering matrix representation of the network, the desired photonic function is designed by fitting structural parameters according to an optimization procedure. This design scheme is demonstrated for plasmonic dichroic and trichroic routers in the infrared frequency range.

©2010 Optical Society of America

## 1. Introduction

The design and development of new photonic components is at the heart of the electrooptics field. Artificial photonic materials allow for the possibility of device designs based on new underlying physics. One way of classifying artificial photonic designs is by their basic resonating elements with which light interacts to give the desired artificial dispersion. Under this classification scheme, we can think of materials that operate based on the nonlocal interference of periodic waves (i.e. photonic crystals) and of materials in which the waves interact locally with subwavelength resonating elements (i.e. metamaterials). The discovery of photonic crystals [1,2] led to extremely high quality factor cavities using all dielectric materials based on the idea of Bragg wave physics [3,4]. The field of metamaterials [5,6] enabled the exciting research area of transformation optics [7] based on the ability to locally tune the material electromagnetic properties. Different in principle from these existing concepts, resonant guided wave networks (RGWNs) were recently reported as a new class of artificial material [8] where localized waves resonate in closed paths throughout a network of isolated waveguides connected by wave splitting elements. The resulting multiple resonances within the network give rise to wave dispersion that is sensitive, and thus tunable, according to the network layout. Here we utilize these distinctive RGWNs properties to formulate a new method for designing photonic components.

A RGWN is comprised of power splitting elements connected by isolated waveguides. The function of the splitting element is to distribute a wave entering any of its terminals between all of its terminals, as illustrated in Fig. 1a . The waves are then propagated in isolated waveguides between the splitting elements, where the local waves from different waveguides are coupled. Four splitting elements arranged in a rectangular network layout form a 2x2 RGWN (see Fig. 1b). When one of the terminals is excited, the multi splitting occurrences of the incident wave within the network form closed path resonances that reshape the dispersion of the emerging waves according to the network layout and is different from the dispersion of the individual waveguides. Furthermore, the interference of the local waves, for a properly designed network layout, considerably enhances the resonance quality factor due to the contribution of coherence effects [8]. The 2x2 RGWN supports one closed loop resonance, however larger two- or three-dimensional networks can support multiple resonances that enrich the network dynamics.

One possible implementation of a RGWN is by using plasmonics via a mesh of intersecting sub-wavelength air gaps in a metal matrix as illustrated in Fig. 1c. Surface plasmon polaritons (referred here as plasmons for brevity) are slow surface waves propagating on metal-dielectric interfaces [9,10]. Metal-insulator-metal (MIM) waveguides support a highly confined plasmon wave (TM_{0} mode) that does not get structural cutoff as the dielectric gap between the metal layers becomes vanishingly small [11–14], allowing for such plasmonic components as high transmission sharp waveguide bends and power splitters [15]. In this implementation, the intersection of two sub-wavelength MIM waveguides forms an X-junction that functions as the power splitting elements in the network [8,16] (see illustration I in Fig. 1c), and the MIM segments between the intersections serve as the isolated waveguides connecting the X-junctions. As the air gap thickness is varied the power splitting between the terminals can be tuned both in terms of amplitude and phase. This, in addition to determining the phase accumulation in the waveguide segments, set independent controls over the design of the interference pattern in the RGWN. The utility of these degrees of freedom was demonstrated in the design of artificial dispersion in homogeneous RGWNs designed to have a frequency band gap at infrared frequencies, where the metal dispersion is negligible [8].

Because the underlying physics of RGWNs is based on the interference of local waves, it allows for layouts that are inhomogeneous and non-periodic across the network. Unlike photonic crystals, which are restricted to Bragg wave effects in periodic structures, the flexibility of RGWNs open up possibilities for designs where the wave properties are varied across the structure. With respect to metamaterials, which could inherently be non-homogeneous due to the local nature of the interaction between light and the meta-atoms, RGWNs have the advantage of having interference effects within the network, which allows for frequency spectrum reshaping designs through interference effects.

An additional unique feature of RGWNs relates to the constraints on wave propagation within the structure. Unlike other photonic designs, RGWNs have a limited number of modes that are allowed to propagate within the structure (e.g., TM_{0} mode). Furthermore, the waves can propagate only inside the waveguides connecting the splitting elements. The different waveguides are coupled only in the X-junctions, which each have only a limited number of terminals. This level of control is beneficial for several reasons. First, the interference pattern in the network can be controlled more directly. Schematically, there is complete control over the combination of the local waves that interfere as well as the phase and the amplitude that they gain in their different trajectories, as illustrated in Fig. 1c (II). Second, it allows for a comprehensive mathematical representation of the RGWN by scattering matrix (S-matrix) formalism which reduces the computational complexity of programming the network, as demonstrated below. Third, since the waveguides are isolated from each other, their only contribution to the network is to serve as phase retardation elements between the splitting elements. As a result, the waveguide length is the only effective parameter in its contour, as long as the bending is not too severe, as is the case for MIM waveguides. This feature of the waveguides allows for the network to maintain its engineered function even when distorted, as schematically illustrated in Fig. 1d. The ability to utilize curved or bent waveguides to accommodate long contours is useful when designing the interference pattern of RGWNs.

These distinctive RGWN characteristics open up new opportunities in designing photonic devices by programming the entire network rather than by assembling interconnected discrete components with traceable functions. The usual way of designing photonic devices is to target the desired subsystem functions, map them logically into sub-functions, and then assembling components that carry out these sub-functions into the desired system. For example, a wavelength router could be designed using add/drop ports where the input and output waveguides are coupled by wavelength sensitive ring resonators [17] or by defects in a photonic crystal [18]. Similarly, in free space optics, this function could be achieved through the use of collinear beam splitters, each designed to deflect a desired wavelengths band. In these schemes, the couplers and waveguides are discrete components that are associated with a specific function, and are combined in a logical way to carry out the overall system function. An alternative approach is to use a network of components that carries out the desired function but, unlike traditional designs, there is no specific logical sub-function associated with any individual component. While the inner connectivity of the device will be less intuitive, it has the potential to result in efficient designs of complex and compact devices. In this Letter we exemplify this design paradigm by programming a RGWN to function as a color router.

## 2. Programming RGWN using S-matrix formalism

One possible way of representing system function in a RGWN is through the use of a scattering operator that maps the set of local waves entering the device terminals to the set of the waves exiting from the same terminals [19]. Since a RGWN is composed of a discrete set of components (waveguides and X-junctions) and terminals, the system function is represented by a scattering matrix (S-matrix) connecting the vectors of the waves entering and emerging from the RGWN via the external ports (see Fig. 2a ). Designing the system function of the RGWN is then mathematically equivalent to designing the S-matrix to yield a desired output, given a set of inputs.

To elucidate the design principle, programming will first be demonstrated for a 2x2 RGWN. This network contains four waveguides intersecting each other in four X-junctions arranged in a 2x2 array. The waveguides are implemented here as air voids in a Au matrix and could differ in width, length, and contour. The device has eight terminals, numbered from ‘1’ to ‘8’, as illustrated in Fig. 2a. The input vector lists the complex amplitudes of the magnetic fields (H-fields) entering the network in the eight terminals, and similarly the output vector describes the complex H-field amplitudes of the waves exiting the network through these same terminals.

The network S-matrix is assembled from the mathematical representation of its components according to the network layout. As a first step, a ‘function library’ of mathematical representations is generated for all the possible network components (i.e., waveguides and X-junctions) using finite difference time domain (FDTD) full wave electromagnetic simulations. Once this library is established, the RGWN S-matrix can be assembled according to its layout. Once the component library is established, the S-matrix calculation scheme is almost always much faster than resolving the RGWN behavior from full wave electromagnetic simulations, yet reproduces the same information about the network. This becomes significant for optimization tasks and especially as the network size increases.

The RGWN contains two types of components, waveguides and X-junctions, which need to be represented mathematically. The waveguides are mathematically represented by their complex phase retardation, determined by the complex propagation constant of the wave and the length of the waveguide. The propagation constants are extracted from FDTD simulations of waveguides with various widths at different frequencies, excited by a continuous wave (CW) modal source. Figures 3a
and 3b depict the real and imaginary parts of the modal effective index, respectively, where n_{eff} = β/k_{0} (the ratio of the wave propagation coefficient to the free space wavenumber). Figure 3a illustrates that as the waveguide size is reduced, the effective index increases, with no significant difference across the spectral range of interest. Figure 3b shows that Ohmic loss in the waveguides will cause modal attenuation with decay length of tens of microns in the near infrared. However, as long as the attenuation length is substantially larger than the closed path of the dominant resonances, the RGWN function will not be significantly altered as demonstrated here and shown in [8]. Throughout this Letter, the Au is material modeled according to tabulated data of its optical coefficients [20].

The other component of the RGWN is the X-junction, comprised of two intersecting waveguides with four terminals, which is mathematically represented by a (4x4) S-matrix. For a given frequency, since the intersection of two waveguides is defined by a small number of structural parameters, namely the waveguide widths and intersecting angle, the parameter space for defining the complete library of X-junction components is greatly reduced. To further simplify our analysis, we set here the two waveguides to be normally intersecting (i.e. setting the intersection angle to 90 degrees). For a given set of waveguide widths, the complex transmission coefficients of the X-junction ports are extracted from FDTD simulations by measuring the amplitude and phase of the wave transmitted to the different ports when excited from one of the terminals at a given wavelength. The complex transmission coefficients into the forward, sideways, and backwards arms are denoted by *t _{F}*,

*t*,

_{S}*t*respectively. Because of the symmetry of the two normally intersecting waveguides, the transmission coefficients of the two side ports are the same. Figures 3c and 3d plot the amplitude and the phase of the transmission coefficients as a function of the two waveguides widths at wavelength of 1.5 μm. When the waveguides widths both equal 200 nm, the incident wave power is split equally (relative field amplitude equals 0.5) and the reflected wave is π phase-shifted from the forward and sideways waves. Increasing the two waveguides widths together (i.e., d

_{B,}_{1}=d

_{2}) deviates from the equal split as more power is transmitted forward; however, having different widths for the two waveguides expands the scope of possible mathematical representation of components as illustrated in Fig. 3.

The S-matrix of the 2x2 RGWN is assembled from the mathematical representation of components according to the network layout. The explicit S-matrix representation is given in appendix B. We validated the complex field predictions of the S-matrix representation with FDTD simulations, finding less than 5% difference for various test cases. The two major contributions to the small deviation result from the interpolation between the parameter space points, where the library components were calculated, and from the error added when the waveguides are bent. For cases where no interpolation or waveguide bending occurs, the FDTD results differ only by 1% from the S-matrix predictions. The ability to accurately predict the RGWN interference using S-matrix representation allows designing it by an efficient optimization scheme.

Programming the RGWN is then performed by minimizing the difference between the network output and the desired one (for a given input) by sweeping the network parameter space, which includes the waveguide widths and lengths. The optimization process will result with a set of network parameters that can be translated to a network layout and then validated with FDTD simulations.

## 3. Designing a color router

We exemplify the S-matrix programming method by designing a 2x2 RGWN to function as a dichroic router (Fig. 4a
). Although simple in concept, setting a passive device to have different functions at different wavelengths is quite instructive. Explicitly, the required function is to route two different wavelengths (*λ _{1}* and

*λ*) to a different set of ports (‘1’ and ‘6’ for λ

_{2}_{1}and ‘2’ and ‘5′ for

*λ*) when the two bottom ports (‘7’ and ‘8’) are simultaneously excited with equal power. Mathematically, we can represent the device as an 8x8 S-matrix

_{2}**S**(λ

_{1,}λ

_{2}) connecting the input and the output vectors. For both wavelengths, the input vector is nonzero for the bottom ports (i.e.

**In**= (0,0,0,0,0,0,1,1)) and the desired output vectors would be

**Out**(λ

_{1}) = (1,0,0,0,0,1,0,0) for

*λ*and

_{1}**Out**(λ

_{2}) = (0,1,0,0,1,0,0,0) for

*λ*. Because we do not have enough degrees of freedom in this small 2x2 network to exactly attain the desired outputs, we optimize the ratio of power going to the two sets of ports at the different wavelengths. Full details of the optimization procedure can be found in Appendix C.

_{2}After defining the optimization function we constrain the parameter space based on practical considerations. The parameter space includes the width and length of the upper, lower, and side waveguides as well as the two wavelengths of operation (*λ _{1}* and

*λ*). We decrease the number of parameters to optimize by restricting the device to have left-right symmetry based on the desired operation. We restrain the design to operate in the infrared frequency range (

_{2}*λ*= 1.2 μm – 2 μm) where the material dispersion and loss are less pronounced than in the visible. Furthermore, the waveguide thickness is constrained to be small enough to only support the lowest order plasmonic mode (

_{0}*d*= 100 nm - 500 nm).

_{1,2}The optimization procedure yields the network parameters, given in Table 1
, which reveal that the required RGWN for color routing is distributed inhomogeneously. When translating the optimization parameters into layout, we learn that the upper waveguide is longer than the lower one, and therefore needs to be bent. Importing the resulting layout into FDTD, we obtain the steady state H-field distribution shown in Fig. 4b and 4c which show time snap shots for the two operation wavelengths. A buildup of local resonance, which results in the filtering out of the desired outputs, is evident from the temporal dynamics of the fields, as is illustrated in Media 1 and Media 2. We note that the transmission (‘3′ and ‘4’) and reflection (‘7’ and ‘8’) ports from the device are not identically zero since the device does not have enough degrees of freedom and therefore are not included in the optimization function. The FDTD simulation results validate the S-matrix design, with *λ _{1}* and

*λ*clearly routed to a different set of sideways ports as illustrated in Fig. 4b and 4c, respectively.

_{2}The matrix representation can also be used to understand the conditions for which the RGWN fulfills its desired function. From the known input vector and the network S-matrix, the wave complex amplitudes can be identified at any point in the network. For each wavelength, we resolve the excitation conditions of the X-junctions that have the ports that are to be filtered out. For example, for *λ _{1}* to be filtered out from terminals ‘2’ and ‘5′, we examine the excitation conditions in X-junction ‘3,’ which has four terminals. Two of the terminals are external device ports (‘4’ and ‘5′) and the other two are internal network terminals. There is no input signal incident on the two external ports, so it is the excitation conditions of the remaining two junction terminals that null the output in terminal ‘5.’ Indeed, the excitation amplitudes of junction ‘3′ obtained from the S-matrix representation are 0.23exp{-j0.21π} and 0.34exp{j0.64π}, which are close in amplitude and ~π phase-shifted. This is consistent with the results from [8], which show that when an ‘X-junction’ is simultaneously excited π phase-shifted from two adjacent terminals, the two other terminals will be filtered out. The excitations are not simply the same in amplitude and π phase-shifted, which is attributed to the additional constraints the design has on the other wavelength as well as the limitations imposed on the parameter space.

Similarly, the excitation conditions necessary for filtering out terminals ‘1’ and ‘6’ at *λ _{2}* (Fig. 4b) are examined by focusing on the S-matrix amplitudes of X-junction ‘4.’ In this case there are three terminals being excited: the lower terminal of the X-junction (port ‘7’) is given by the network excitation, so the excitation of the other two internal ports will determine the filtering out of port ‘6.’ Intuitively, the condition to filter out terminal ‘6’ will be simply a π phase-shifted excitation of the upper and lower terminals of junction ‘4’, with zero excitation from the side port. From the case of

*λ*we also know that additional constraints might cause a residual wave emerging from terminal ‘6’, which could be compensated by a small amplitude excitation at the other side terminal of the junction ‘4.’ Indeed, the excitation amplitudes of junction ‘4’ in the S-matrix representation are 1 in lower terminal, 0.9exp{j0.82π} in upper terminal, and 0.3exp{-j0.32π} in the side terminal, which is consistent with the intuitive picture and is evident from Media 1 and Media 2.

_{1}To further exemplify the programmability of RGWNs via S-matrix formalism, we consider a 3x3 RGWN programmed to function as a trichroic router. In order to implement the more complex task of routing three wavelengths we allow for more degrees of freedom in the network by increasing the number of components, effectively increasing the amount of data contained. The function is defined as an extension of the dichroic router, but here when the three bottom terminals are simultaneously excited at three different frequencies, the frequencies are filtered out to three different sets of side terminals as illustrated in Fig. 5 . The analysis results in the optimal RGWN parameters shown in Table 2 . FDTD simulations validate the resulting layout design, as shown in Fig. 5.

## 4. Conclusion

In this Letter we have implemented a scheme for programming resonant guided wave networks by exploiting their unique interference properties. The constraints on the propagation and coupling of the local waves in RGWNs allow for the device operation to reduce to a simple mathematical representation using S-matrix formalism. This allows for the network programming to take the form of an optimization procedure over a relatively small parameter space. The RGWN S-matrix representation was demonstrated here where the inputs are given and the S-matrix of the device is designed to give the desired outputs (e.g., routing, modal converters). However, this formalism could be extended to different type of functions, such as sensing, in which the inputs are given and the output changes are monitored.

The specific RGWN function examined here is color routing, which is challenging in terms of having to program a different function for different wavelengths all within the same device. Even though it was meant as a concept demonstration, incorporating more components into the RGWN and therefore increasing the possible degrees of freedom, may evolve to a more complex device or alternatively a device with enhanced performance.

## Appendix A. Tunabillity of the wave properties in RGWN components.

The three-dimensional data of Fig. 3 is plotted in Fig. A1 in two-dimensional contours to allow for easier read in a quantitative manner.

## Appendix B. S-matrix representation of 2x2 RGWN

This appendix describes the assembly of the 2x2 RGWN S-matrix from the properties of its components according to the network layout.

The phasor representation of the local wave H-fields in the network is represented by three column vectors (transposed for brevity):

_{out}and A

_{in}hold the values of the local input and output waves of the RGWN at its ports, and A

_{net}represents the input wave on the X-junctions from the internal terminals of the RGWN. The

*i/o*denotes input/output waves in respect to the X-junction, the junctions are numbered as in Fig. 2a, and the junction ports are numbered in brackets, in accordance to Fig. 2c.

The coupling of the H-field vectors by the network connectivity is represented with the system:

_{FS}and M

_{RS}are diagonal 8-by-8 matrices that originate from the splitting relations in the X-junctions and K is a sparse 8-by-8 matrix that stands for the wave propagation in the waveguides. These matrices are defined as:

Algebra of equation set B-2 gives the matrix representation of the 2x2 RGWN S-matrix:

## Appendix C. Optimization procedure

The optimization procedure is implemented in Matlab using the pre-calculated mathematical representation data set of the RGWN components obtained from full-field electromagnetic FDTD simulations excited with continuous wave sources. The procedure is illustrated with the following flow chart in Fig. C1.

FDTD is used to calculate the complex effective index of Au/Air/Au MIM waveguides with air gap thickness ranging from 100 nm to 520 nm in steps of 10 nm. The effective index is calculated from the mode field attenuation and its phase retardation relative to the speed of light. The S-matrix data is calculated for X-junctions of two MIM intersecting waveguides with air gap thickness ranging from 100 nm to 520 nm in steps of 10 nm (with the two waveguides not necessarily having the same width). The transmission coefficients are calculated by monitoring the relative field transmitted to each of the ports relative to the incident field through any one of its ports.

The dichroic router network is defined by eight parameters: the length and width of the upper, lower and side waveguides and the two wavelengths. The waveguide widths determine the effective index in the waveguides as well as the transmission coefficients of the X-junctions. The optimization procedure is conducted using Matlab *fmicon* function using ‘interior point’ settings and a random set of 200 initial points over the parameter space to reduce local minimum point effects. The main script starts at an initial point within the parameter space and use the *fmicon* function to move in the parameter space towards minimization of the optimization function. The optimization function for the dichroic router is defined as follows:

*O*and

_{1}*O*represent the two terminals output ratios that need to be maximized at the two different wavelengths,

_{2}*λ*and

_{1}*λ*. The

_{2}*f*function is used to merge the two ratios together into one weighted optimization parameter, where the first term in squared brackets maximizes the total power routed into the selected terminals, the second term maximizes the two ratios, and the third term is a weighting factor that ensures that the two ratios are maximized equally. The target function is defined as the inverse of these three terms in multiplication since the

*fmicon*function is a minimization function. At each point in the parameter space, the network output vector is calculated as the multiplication of the 2x2 RGWN S-matrix evaluated at the parameter values times the input vector representing excitation only from terminal ‘8’. The optimized network parameters are given in Table 1 and the corresponding output field amplitude values are given in Table C1.

For the Trichroic function the optimization function is:

The logic of constructing the function is similar to Table C1 but for weighting the three ratios to be maximized. The target function in (C-2) is optimized yielding the parameter values given in Table 2. The output field amplitude values for the trichroic router are given in Table C2.

## Acknowledgments

We wish to acknowledge the support from the Office of Basic Energy Sciences under Contract No. DE-FG02-07ER46405.

## References and links

**1. **E. Yablonovitch, “Photonic crystals: semiconductors of light,” Sci. Am. **285**(6), 46–51, 54–55 (2001). [CrossRef]

**2. **J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of light, 2’nd Ed. (New Jersey, Princeton, 2008).

**3. **M. Fujita, S. Takahashi, Y. Tanaka, T. Asano, and S. Noda, “Simultaneous inhibition and redistribution of spontaneous light emission in photonic crystals,” Science **308**(5726), 1296–1298 (2005). [CrossRef] [PubMed]

**4. **J. Scheuer, W. M. Green, G. A. DeRose, and A. Yariv, “Lasing from a circular Bragg nanocavity with an ultra small modal volume,” Appl. Phys. Lett. **86**(25), 251101 (2005). [CrossRef]

**5. **D. R. Smith, J. B. Pendry, and M. C. K. Wiltshire, “Metamaterials and negative refractive index,” Science **305**(5685), 788–792 (2004). [CrossRef] [PubMed]

**6. **V. M. Shalaev, “Optical negative-index metamaterials,” Nat. Photonics **1**(1), 41–48 (2007). [CrossRef]

**7. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science **312**(5781), 1780–1782 (2006). [CrossRef] [PubMed]

**8. **E. Feigenbaum and H. A. Atwater, “Resonant guided wave networks,” Phys. Rev. Lett. **104**(14), 147402 (2010). [CrossRef] [PubMed]

**9. **W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon subwavelength optics,” Nature **424**(6950), 824–830 (2003). [CrossRef] [PubMed]

**10. **E. N. Economou, “Surface Plasmons in Thin Films,” Phys. Rev. **182**(2), 539–554 (1969). [CrossRef]

**11. **B. Prade, J. Y. Vinet, and A. Mysyrowicz, “Guided optical waves in planar heterostructures with negative dielectric constant,” Phys. Rev. B Condens. Matter **44**(24), 13556–13572 (1991). [CrossRef] [PubMed]

**12. **R. Zia, M. D. Selker, P. B. Catrysse, and M. L. Brongersma, “Geometries and materials for subwavelength surface plasmon modes,” J. Opt. Soc. Am. A **21**(12), 2442–2446 (2004). [CrossRef]

**13. **J. A. Dionne, H. J. Lezec, and H. A. Atwater, “Highly confined photon transport in subwavelength metallic slot waveguides,” Nano Lett. **6**(9), 1928–1932 (2006). [CrossRef] [PubMed]

**14. **E. Feigenbaum and M. Orenstein, “Modeling of Complementary (Void) Plasmon Waveguiding,” J. Lightwave Technol. **25**(9), 2547–2562 (2007). [CrossRef]

**15. **G. Veronis and S. Fan, “Bends and splitters in metal-dielectric-metal subwavelength plasmonic waveguides,” Appl. Phys. Lett. **87**(13), 131102 (2005). [CrossRef]

**16. **E. Feigenbaum and M. Orenstein, “Perfect 4-way splitting in nano plasmonic X-junctions,” Opt. Express **15**(26), 17948–17953 (2007). [CrossRef] [PubMed]

**17. **B. E. Little, S. T. Chu, H. A. Haus, J. Foresi, and J. P. Laine, “Microring resonator channel dropping filters,” J. Lightwave Technol. **15**(6), 998–1005 (1997). [CrossRef]

**18. **S. Fan, P. Villeneuve, J. Joannopoulos, and H. Haus, “Channel drop filters in photonic crystals,” Opt. Express **3**(1), 4–11 (1998), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-3-1-4. [CrossRef] [PubMed]

**19. **S. E. Kocabas, G. Veronis, D. Miller, and S. Fan, “Transmission Line and Equivalent Circuit Models for Plasmonic Waveguide Components,” IEEE J. Sel. Top. Quantum Electron. **14**(6), 1462–1472 (2008). [CrossRef]

**20. **E. D. Palik, Handbook of optical constants of solids, 2'nd Ed. (San-Diego, Academic, 1998).