## Abstract

In this work, it is shown how the shapes of surface plasmon dispersion curves can be engineered by manipulating the distribution of the electromagnetic fields in multilayer structures, which themselves are controlled by the free electron density in metal-like materials, such as doped semiconductors in the THz spectral range. By having a nonuniform free electron density profile, reduced relative to that in typical bulk metals, the electromagnetic fields of surface plasmons are distributed in different metallic materials that have different complex dielectric permittivities. As the in-plane component of surface plasmon’s wave-vector increases, they become more confined to a particular layer of the multilayer structure and have energies that are predictable by considering the permittivity of the layer in which the fields are most concentrated. Unusual and arbitrary shapes of surface plasmon dispersion curves can be designed, including stair steps and dovetails shapes.

© 2013 Optical Society of America

## 1. Introduction

Surface plasmons (SPs) have been extensively researched for many years and have been shown to have interesting electromagnetic properties that have been used in a growing number of applications [1]. Properties such as light concentrating [2], directing and filtering of light [3], surface plasmon lasers or spasers [4], and many other properties and applications of SPs have been previously discussed [1]. There are several types of metal/dielectric structures that support SPs that have differently shaped dispersion curves (DCs), including flat metal/dielectric interfaces [5], structures with different combinations of insulator and metal films [6], and differently shaped metallic particles [7]. SPs in layered or stratified media, have also been studied in great detail with SP Bloch modes and other modes being possible in such structures [8]. Yet for all of these structures, the metals are assumed to have electric permittivities *ε _{m}* that are independent of position within the metal, thus limiting the shapes of the SP DCs that one can achieve with these structures. The ability to tailor SP DCs through structural and compositional design is attractive for the study of slow and fast light phenomena [9–12], all-optical analogs of gravitational black and white holes [13,14], and many other light controlling phenomena and applications that are receiving increasing interest [15,16]. These phenomena and applications require structures that support electromagnetic (EM) modes that have group velocities that transition from positive to negative (or vice versa) as the in-plane wave vector component (i.e.,

*k*) increases.

_{x}For a flat metal/dielectric interface, the SP DC is given by the equation [5]:

*k*=

_{o}*ω/c*, and

*ε*and

_{m}*ε*being the electric permittivities of the metal and dielectric materials, respectively. For metals,

_{d}*ε*is dependent on the frequency, with this frequency dependence often described by the Drude model:

_{m}*ω*is the plasma frequency and Γ =

_{p}*e*/(

*m*) is the damping rate with

_{eff}μ_{e}*e*as the electron charge, and

*m*and

_{eff}*μ*being the effective mass and mobility of the charges respectively. The plasma frequency

_{e}*ω*= 2

_{p}*πν*is given by ${\omega}_{p}^{2}=n{e}^{2}/{\epsilon}_{o}{m}_{\mathit{eff}}$, where

_{p}*n*is the number density of free charge carriers and

*ε*is the free-space permittivity.

_{o}In this work, we introduce a method to engineer *complex shapes* in the SP DCs by using flat multilayer structures, composed of metal-like materials, in which one controls how the electromagnetic (EM) fields of SPs are distributed within the structure.

The EM field distribution is dependent on the decay length *l* of the SP from the dielectric-metal interface, which in turn is dependent on the SP momentum *k _{x}*. As

*k*increases (with ${k}_{x}>\sqrt{\epsilon}{k}_{o}$), the value of

_{x}*k*is imaginary with increasing magnitude and the decay length

_{y}*l*decreases (with

*l*defined as

*k*=

_{y}*i/l*), as predicted by the equation ${k}_{y}=\sqrt{\epsilon {k}_{o}^{2}-{k}_{x}^{2}}$ in a medium with a permittivity

*ε*. Thus, for SPs on a flat metal/dielectric interface, as

*l*decreases the fields are more heavily concentrated in regions of the structure that are closer to the dielectric/metal interface, and the effects on the SP field distribution due to the layers farther than

*l*from the interface will become less significant, and the SP frequency will depend on an effective

*ε*of the layers within a distance

*l*of the interface. When

*l*<

*t*with

_{d}*t*the thickness of the dielectric layer right above the metal, the SPs will assume (i.e., take on) the frequency appropriate for

_{d}*ε*of that layer. Alternatively, one can consider how the fields behave as

*k*is decreased from a large to progressively smaller values. In this case the fields are necessarily spread out in the

_{x}*y*direction for decreasing

*k*and experience the effects of more layers in the structure; and for small values of

_{x}*k*, the SPs modes will assume frequencies appropriate for the

_{x}*ε*of the thickest layer. Therefore, the EM field control and SP DC engineering can be performed by controlling

_{m}*ε*as a function of position; and this can be done using structures with either multiple dielectric layers atop a metal, or multiple metal layers atop a dielectric, or some combination thereof. We find in this work that the most effective method to achieve to a large degree of control of the shape of SP DCs is to use doped semiconductors for which

*ε*at a specific position is determined by the free electron density

*n*at that position [17].

In this work, a rigorous couple wave algorithm (RCWA) and scattering matrix method was used to calculate the DCs of SPs. In the RCWA, the reflection and transmission coefficients are determined by propagating the incident field through the various material layers using transfer matrices as outlined in Moharam *et al.* [18]. In the scattering matrix method, a matrix *S* is constructed using transfer matrices that relate the tangential field components at the interface between two adjacent layers as presented in Sambles *et al.* [19], with *S* providing relations between the amplitudes of the input waves and output waves of the structure:

*A*and and

_{t}*B*are the amplitudes of the input beams in the top-most and bottom-most layers respectively, and

_{b}*B*and

_{t}*A*are the output beams in the top-most and bottom-most layers respectively. Ideally (i.e., without absorption or scattering), SP modes and other surface confined resonant modes are self-sustaining and have nonzero values of

_{b}*B*and

_{t}*A*without being excited by the input beams. Equation (3) can be rewritten to stress this as:

_{b}In order for Eq. (4) to allow this solution, the determinant of *S*^{−1} (i.e., |*S*^{−1}|) has to be zero which indicates a pole in the scattering matrix *S*. Thus, to identify the frequencies and *k _{x}* values of SPs, one determines the locations in

*ω*−

*k*space where |

_{x}*S*

^{−1}| = 0 [20]. Additionally, sharp dips or peaks in the reflectance |

*R*|, with

*R*=

*B*=

_{t}*S*

_{11}

*A*, also indicate the excitation of SPs. This is due to the poles in

_{t}*S*that occur when |

*S*

^{−1}| = 0. The solution frequencies and momentum for the roots of the equation |

*S*

^{−1}| = 0 may be complex, however, if the imaginary parts are small, the associated poles in

*S*would produce dips and peaks with bandwidths and lifetimes directly related to the imaginary parts of the solution [20, 21].

## 2. Methods

#### 2.1. Multiple dielectric layers

The most obvious approach to engineer the DCs of SPs is to use multiple dielectric layers below (or above) a single metal layer. With this approach, as *k _{x}* increases (with
${k}_{x}>\sqrt{{\epsilon}_{d}}{k}_{o}$), the EM fields of the SPs go from being distributed largely within the thicker bottom-most dielectric layer to being confined to the thinner dielectric layer immediately below the metal (see Fig. 1). An example of this is shown in Fig. 1, which shows a structure consisting of a thin 2 nm thick dielectric material with

*ε*= 13 placed immediately below a 10 nm thick silver (Ag) film, with the superstrate and substrate as vacuum. The Ag film is assumed to have a Drude-like permittivity with

_{d}*ω*= 1.32 × 10

_{p}^{16}rad/s and

*τ*= 1.45 × 10

^{−14}s. The DCs of the odd and even SP modes in this insulator-metal-insulator (IMI) structure are well known [22, 23], with the odd SP modes being less confined to the metal/dielectric interface and having higher energies compared to the even SP modes. Thus, the thin

*ε*= 13 dielectric layer is expected to affect (in this case lower) the frequencies of the even SP modes at high

_{d}*k*values to a much larger extent than the frequencies of odd SP modes. With this IMI structure, the turn-over of the DC of the odd SP modes is typical, where the term ”turn-over” refers to the cresting of the SP DC followed by a gradual decrease in energy as

_{x}*k*increases. In this work however, the motivation for putting the thin

_{x}*ε*= 13 dielectric film below and in contact with the Ag is to turn over the DC of the even SP modes. There are two parameters of the dielectric film that one can adjust to achieve this effect, the thin film’s dielectric constant and its thickness. The dielectric constant of the film determines the extent of the raising (if

_{d}*ε*<

_{d}*ε*) or lowering (if

_{superstrate}*ε*>

_{d}*ε*) of the energies of the SP modes. The thickness of the film determines at what

_{superstrate}*k*value this raising/lowering will start to occur (the thinner the film, the higher the

_{x}*k*value needed before the EM field of the even SP mode is primarily concentrated in the thin dielectric film). While it is possible to slightly lower the frequency in the asymptotic limit of the of the even SP mode (i.e,

_{x}*ω*as

_{sp}*k*→ ∞), there is only a small turn-over of the even SP mode (relative to what is possible with the second approach described in the rest of this work). And this relatively small turn-over is achieved only after using an unrealistically small thickness for the thin film and very large dielectric constant differences between the semi-infinite superstrate and the thin dielectric layer. Thus another approach is needed to perform complex and interesting DC engineering.

_{x}#### 2.2. Multiple conductive layers

The other approach is to use multiple metal layers or metal-like layers. For real metals, such as gold, copper, or silver, the values for *ε _{m}* are so large from the visible to microwave spectral ranges that the EM fields of the SPs are compressed into a very thin layer in the metals. This greatly limits the ability of the EM field to extend to lower-lying metal layers and incorporate the properties of these layers (i.e., their values of

*ε*) into the SP DC. However, if one could control

*ε*to have smaller real negative values, then this approach may prove fruitful. Such control is possible with semiconductors in the THz frequency range. Much work has been done on using doped semiconductors to control

_{m}*ε*and SPs using semiconductors, e.g., silicon, gallium arsenide, and indium antimonide [17]. In these semiconductors, one can precisely control the free charge carrier density

*n*by doping them with varying concentrations of dopants (assuming high enough temperatures such that all the dopants are ionized). A range for

*n*from 0 (for low temperature, intrinsic semiconductors) to 2 × 10

^{20}

*cm*

^{−3}is possible with phosphorus doped silicon. These materials and concepts can be extended to allow for complex shapes in the SP DCs, as describe in this work.

In order for the semiconductor to support well-defined, low-loss SP modes, *ω _{p}* should be significantly greater than Γ. Of the typical semiconductors, InSb is the most promising. It has a high room temperature value of

*μ*of 70,000

_{e}*cm*

^{2}/

*Vs*that can rise to 200,000

*cm*

^{2}/

*Vs*at low temperatures and doping concentrations [24]. Yet tempering this enthusiasm for InSb is the fact that it has a very small value for

*m*of 0.015

_{e}*m*that undesirably increases the value of Γ (with Γ = 1/

_{o}*μ*), such that Γ = 1.52 × 10

_{e}m_{eff}^{12}rad/s at room temperature and Γ = 5.86 × 10

^{11}rad/s for low temperature. Gallium arsenide and silicon are less promising due to their smaller values of

*μ*, and doping these materials to raise the value of

_{e}*ω*relative to Γ may not improve the matter because higher doping levels lowers

_{p}*μ*and therefore increases Γ.

_{e}Figures 2(a) and 2(b) show the structure of a multilayer InSb structure under consideration and the permittivities of the three InSb layers. The structure consists of three doped InSb layers with Drude parameters *ε*_{∞} = 15.7 and layer-specific values consisting of: *n*_{2} = 1.9 × 10^{15}*cm*^{−3} and *t*_{2} = 15 *nm* (*ω _{p}* = 2.01 × 10

^{13}

*rad/s*) for layer 2;

*n*

_{3}= 2.6 × 10

^{15}

*cm*

^{−3}and

*t*

_{3}= 20

*nm*(

*ω*= 2.35 × 10

_{p}^{13}

*rad/s*) for layer 3;

*n*

_{4}= 4.3 × 10

^{15}

*cm*

^{−3}and

*t*

_{4}= 1

*μm*(

*ω*= 3.02 × 10

_{p}^{13}

*rad/s*) for layer 4. Layers 1 and 5 are lossless, dispersionless dielectrics with

*ε*

_{1}=

*ε*

_{5}= 4. To achieve this low free carrier concentration, the InSb would have to be at a low temperature (i.e., 77 K) because the intrinsic carrier concentration at room temperature is 2 × 10

^{16}

*cm*

^{−3}. Thus, for low-doped, cooled InSb, a value for the mobility of 1.5 × 10

^{5}

*cm*

^{2}/

*Vs*is appropriate [24].

Different doping densities can be chosen to adjust the frequencies that the SP modes asymptotically approach for large *k _{x}* values; these would be the frequencies of each ”stair-step” in the DC (i.e.,

*ν*) (see Fig. 3). Letting

_{step}*k*→ ∞ in Eq. (1), one finds that

_{x}*ε*= −

_{m}*ε*, and then using Eq. (2) for

_{d}*ε*, the asymptotic frequencies of the

_{m}*i*th layer

*ν*of the SP DCs as a function of

_{step,i}*ω*,

_{p,i}*ε*

_{∞}and

*ε*can be obtained:

_{d}The values for *ε′ _{m}* and −

*ε′*as a function of frequency and the values of

_{m}/ε″_{m}*ν*for the three semiconductor layers are shown in Fig. 2(b). It is seen that for a frequency of

_{step}*ν*

_{step}_{,4}, the above-lying Layers 2 and 3 both have

*ε′*> 0 and act as dielectrics that better allow the fields to penetrate to Layer 4, which they will for small

_{m}*k*values, thus forming Step 4. For a frequency of

_{x}*ν*

_{step}_{,3}Layer 3 is now metallic-like with

*ε′*< 0, and Layer 2 is still dielectric-like; and for high

_{m}*k*values, the fields cannot penetrate deeply enough to reside in Layer 4 to any significant degree, thus Step 3 is formed. Finally, for even higher

_{x}*k*values, the fields are even more tightly constrained to the top interface and only reside in and experience Layer 2, with the frequencies of the SPs assuming the value of

_{x}*ν*

_{step}_{,2}. Also, it should be noted that for a structure to support well-defined SPs modes that exists for any appreciable time, −

*ε′*should be greater than 1 [25]. In addition, for the fields to be able to differentiate between layers of such thin thickness instead of seeing an effective dielectric material, the distance over which the fields change in space must be of the same order of the thickness, and therefore will occur only at high

_{m}/ε″_{m}*k*values where the decay constant associated with the SP mode is very small.

_{x}A rough approximation can be obtained for the values of *k _{x}* at which the step changes occur (i.e.,

*k*). Because the SPs have a

_{step,m→n}*k*(with ${k}_{y}=i\kappa =\sqrt{\epsilon {k}_{o}^{2}-{k}_{x}^{2}}$) that is imaginary in each layer, the fields have a

_{y}*y*dependence that is exponentially decreasing with

*y*with a decay constant of

*κ*. The properties of a particular lower-lying layer is not appreciably experienced by the fields when the field intensity at the top surface of this layer is at a decayed value of

*e*

^{−1}relative to its incident intensity. Thus the equation that determines

*k*

_{step}_{,4→3}is ${t}_{2}\sqrt{{k}_{x}^{2}-{\epsilon}_{2}{k}_{o}^{2}}+{t}_{3}\sqrt{{k}_{x}^{2}-{\epsilon}_{3}{k}_{o}^{2}}=1$ where

*ε*

_{2},

*ε*

_{3}and

*k*are calculated assuming a frequency

_{o}*ν*

_{4}given by Eq. (5). The equation that determines

*k*

_{step}_{,3→2}is ${t}_{2}\sqrt{{k}_{x}^{2}-{\epsilon}_{2}{k}_{o}^{2}}=1$ where

*ε*

_{2}and

*k*are calculated assuming a frequency

_{o}*ν*

_{3}given by Eq. (5). In the InSb discussed in this paper,

*k*is typically much larger than

_{x}*εk*in any layer, thus the conditions for

_{o}*k*

_{step,m}_{→}

*can be simplified to ${k}_{\mathit{step},4\to 3}=\frac{1}{{t}_{2}+{t}_{3}}$ and ${k}_{\mathit{step},3\to 2}=\frac{1}{{t}_{2}}$ which predict the transitions to occur for*

_{n}*k*= 28.5

_{x}*μm*

^{−1}and

*k*= 66.7

_{x}*μm*

^{−1}respectively.

## 3. Results

Figure 3(a) shows the frequency and *k _{x}* dependence of

*R*, and shows the stair-step shaped SP DCs. Figure 3(b) shows |

*S*| as a function of

*ν*and

*k*; the black curve indicates the locations where |

_{x}*S*| = 0 that indicate SP modes. SPs produce maxima in the reflectance that largely match the SP DC given by |

*S*|.

The magnetic fields *H _{z}* of the SP modes for three different SP modes with (

*ν*= 1.02 THz,

*k*= 2.5

_{x}*μm*

^{−1}) for Step 4, (

*ν*= 0.91 THz,

*k*= 18

_{x}*μm*

^{−1}) for Step 3, and (

*ν*= 0.75 THz,

*k*= 50

_{x}*μm*

^{−1}) for Step 2 are shown in Figs. 4 and 5. It is seen that the fields behave as expected, with SPs of higher

*k*values being increasingly confined within layers closer to the top-most dielectric/semiconductor interface and assuming (i.e., acquiring) frequencies appropriate for the doping concentrations of these layers.

_{x}In this example structure, the layers lying above the layer primarily responsible for a step have values of *ε′ _{m}* that are greater than zero, and thus are dielectric in character and therefore allow more of the incident field to penetrate to the step-producing metal-like layer. This is not a necessary condition however, and more complicated DCs can be obtained by having the doping density vary in arbitrary ways. For example, if the doping level of Layer 2 was made to be larger than that of Layer 3, the DC would curve back up to higher frequencies at

*k*=

_{x}*k*

_{step}_{,3→2}and produce a dispersion curve with a ”dovetail” shape as can be seen in the case of Fig. 6.

## 4. Discussion

It is clearly possible to use a structure composed of different metallic layers to produce two different frequencies for SP modes with high and low *k _{x}* values. Yet, what is less clear is what happens at the transition region between the one stair step and the next. In Fig. 3(b), it appears that the SP DC overshoots the underlying stair step such that for certain

*k*values, there are multiple possible SP modes with different frequencies and group velocities. The transitions between modes shown occur at lower values of

_{x}*k*(approximately

_{x}*k*= 15

_{x}*μm*

^{−1}and

*k*= 30

_{x}*μm*

^{−1}) than predicted using the rough approximation which are off by a factor of 2, the approximation would be more accurate with higher order terms. The occurrence of three SP modes at particular

*k*values is possible due to the possibility of having distinct, largely independent SP modes excited within each metallic layer. Both of these details of this SP-SP transition coupling and the SP modes at the ends of each stair step need to be studied in greater detail, but are outside the scope of this paper. One possibility is that the SP modes at the ends of each stair step simply fade away, are poorly defined, or cannot be excited, thus making this troubling but interesting issue moot.

_{x}The real material system that was used in this work involved InSb at low doping levels and low temperature. However other temperatures, doping conditions and even different materials can potentially be used, as long as *ω _{p}* is significantly larger than Γ and that

*ε*is relatively slowly varying. Quantifying these statements is not within the scope of this work, but it is expected that other structures should satisfy these necessary conditions that allow for SP DC engineering. For example, one can use InSb at room temperature instead of low temperature and try to compensate the lowering of

*μ*from 200, 000

_{e}*cm*

^{2}/

*Vs*at low temperature to 77, 000

*cm*

^{2}/

*Vs*at room temperature by increasing the doping. There appears to be an optimal doping level that yields a maximum in −

*ε′*. Above this optimal doping level, −

_{m}/ε″_{m}*ε′*decreases while

_{m}/ε″_{m}*ε*continues to become more rapidly varying; both of these aspects are undesirable for dispersion engineering of SPs described in this work. But if materials and conditions can be found to realize this approach, it should provide a powerful way to dial-in a desired shape for the SP DC.

## 5. Conclusion

In this work, it was shown that structures with multiple metallic layers, each with different electric permittivity values can produce SP DCs that have unusual stair-step shapes. By tuning the free electron density within multi-layer doped semiconductors, the DCs can be designed to have a wide variety of complex shapes with multiple excitable energy levels at particular wavevector intervals.

## Acknowledgments

This work is supported by the Air Force Office of Scientific Research project "Hybrid Metamaterials for Solar Biofuel Generation" (grant number FA9550-10-1-0350) and the National Science Foundation Industry University Cooperative Research Center, Center for Metamaterials (grant number IIP-1068028).

## References and links

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

**2. **J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, “Plasmonics for extreme light concentration and manipulation,” Nat. Mater. **9**, 193–204 (2010). [CrossRef] [PubMed]

**3. **F. García-Vidal and L. Martín-Moreno, “Transmission and focusing of light in one-dimensional periodically nanostructured metals,” Phys. Rev. B **66**, 155412 (2002). [CrossRef]

**4. **K. Li, X. Li, M. Stockman, and D. Bergman, “Surface plasmon amplification by stimulated emission in nanolenses,” Phys. Rev. B **71**, 1–5 (2005). [CrossRef]

**5. **H. Raether, *Surface Plasmons on Smooth and Rough Surfaces and on Gratings*, Vol. 11 of Springer Tracts in Modern Physics (Springer-Verlag, 1988).

**6. **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**, 2442–2446 (2004). [CrossRef]

**7. **V. Chegel, Y. Demidenko, V. Lozovski, and A. Tsykhonya, “Influence of the shape of the particles covering the metal surface on the dispersion relations of surface plasmons,” Surf. Sci. **602**, 1540–1546 (2008). [CrossRef]

**8. **P. Yeh, *Optical Waves in Layered Media*, Vol. 95 of Wiley Series in Pure and Applied Optics(Wiley, 1988).

**9. **J. Wang, C. Fan, P. Ding, J. He, Y. Cheng, W. Hu, G. Cai, E. Liang, and Q. Xue, “Tunable broad-band perfect absorber by exciting of multiple plasmon resonances at optical frequency,” Opt. Express **20**, 14871–14878 (2012). [CrossRef] [PubMed]

**10. **Y. A. Vlasov and S. J. McNab, “Coupling into the slow light mode in slab-type photonic crystal waveguides,” Opt. Lett. **31**, 50–52 (2006). [CrossRef] [PubMed]

**11. **M. S. Bigelow, N. N. Lepeshkin, and R. W. Boyd, “Superluminal and slow light propagation in a room-temperature solid,” Science **301**, 200–202 (2003). [CrossRef] [PubMed]

**12. **Y. Okawachi, M. Bigelow, J. Sharping, Z. Zhu, A. Schweinsberg, D. Gauthier, R. Boyd, and A. Gaeta, “Tunable all-optical delays via Brillouin slow light in an optical fiber,” Phys. Rev. Lett. **94**, 1–4 (2005). [CrossRef]

**13. **T. G. Philbin, C. Kuklewicz, S. Robertson, S. Hill, F. König, and U. Leonhardt, “Fiber-optical analog of the event horizon,” Science **319**, 1367–1370 (2008). [CrossRef] [PubMed]

**14. **K. L. Tsakmakidis, A. D. Boardman, and O. Hess, “‘Trapped rainbow’ storage of light in metamaterials,” Nature **450**, 397–401 (2007). [CrossRef] [PubMed]

**15. **R. W. Boyd, “Slow and fast light: fundamentals and applications,” J. Mod. Opt. **56**, 1908–1915 (2009). [CrossRef]

**16. **P. W. Milonni, “Controlling the speed of light pulses,” J. Phys. B At. Mol. Opt. Phys. **35**, R31 (2002). [CrossRef]

**17. **P. West, S. Ishii, G. Naik, N. Emani, V. Shalaev, and A. Boltasseva, “Searching for better plasmonic materials,” Laser Photonics Rev. **4**, 795–808 (2010). [CrossRef]

**18. **M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A **12**, 1077–1086 (1995). [CrossRef]

**19. **N. P. K. Cotter, T. W. Preist, and J. R. Sambles, “Scattering-matrix approach to multilayer diffraction,” J. Opt. Soc. Am. A **12**, 1097–1103 (1995). [CrossRef]

**20. **M. Nevière, E. Popov, and R. Reinisch, “Electromagnetic resonances in linear and nonlinear optics: phenomenological study of grating behavior through the poles and zeros of the scattering operator,” J. Opt. Soc. Am. A **12**, 513–523 (1995). [CrossRef]

**21. **D. A. Bykov and L. L. Doskolovich, “Numerical methods for calculating poles of the scattering matrix with applications in grating theory,” J. Lightwave Technol. **31**, 793–801 (2013). [CrossRef]

**22. **S. Maier, *Plasmonics - Fundamentals and Applications*, Vol. 76 of Chemistry and Materials Science (Springer, 2010).

**23. **A. Karalis, E. Lidorikis, M. Ibanescu, J. Joannopoulos, and M. Soljačić, “Surface-plasmon-assisted guiding of broadband slow and subwavelength light in air,” Phys. Rev. Lett. **95**, 1–4 (2005). [CrossRef]

**24. **E. Litwin-Staszewska, W. Szymanska, and P. Piotrzkowski, “The electron mobility and Thermoelectric Power in InSb at Atmospheric and Hydrostatic Pressures,” Phys. Status Solidi B **106**, 551–559 (1981). [CrossRef]

**25. **A. K. Azad, Y. Zhao, and W. Zhang, “Transmission properties of terahertz pulses through an ultrathin subwavelength silicon hole array,” Appl. Phys. Lett. **86**, 141102 (2005). [CrossRef]