## Abstract

A method to engineer *complex shapes* into the dispersion curves (DCs) of surface plasmon (SP) in flat multilayer structures composed of doped semiconductors is shown analytically and numerically. It is shown that the shapes of SP DCs in multilayer structures can be engineered by controlling not only the free electron density in each layer, as has been well described in past literature on the subject, but also by controlling the spatial profiles of the electromagnetic fields associated with surface plasmons. It is the ability to impart *complex shapes* in the SP DCs that is new and described in this paper—and not just raising or lowering the surface plasmon energies as a function of free electron concentration. DCs with stair-step or dove-tail shapes are possible. In this work, a method to design structures that support SPs with DCs of arbitrary shapes is described. An analytical description of the design methodology is developed, and the resulting structures are numerically modeled using continuous wave excitation and time-dependent pulsed excitation. Time domain studies of the creation and decay in which the SPs are performed, along with spatial sorting of the SPs, yield results in agreement with the predictions of the analytic model.

© 2015 Optical Society of America

## 1. INTRODUCTION

Surface plasmons (SPs) are electron density waves that propagate along the surface of a conductor [1,2] and have been extensively explored for use in a number of applications [3
–9]. There have been many types of metal/dielectric structures studied that support SPs with dispersion curves (DCs) of complex shapes, including differently shaped particles covering the metal surface [10], combinations of insulator and metal films [11], flat metal/dielectric interfaces for plasmonic slab waveguides [12], multiscale ferromagnetic nanowires metamaterials [13], and layered or stratified dielectric films atop a metal [14]. Such dispersion engineering of electromagnetic (EM) modes is of great use for the study of slow and fast light phenomena [15,16]. However, the electric permittivity ${\u03f5}_{m}$ of the metal in the structures mentioned above is independent of position (i.e., a single metal is used), and this generally limits the shapes of the SP DCs that one can realize with these structures. The exception to this is the work by Karalis *et al.* [17,18], which uses multiple dielectric layers atop a metal film in an approach very similar to the work described in our previous work [19] and the work described in this paper that uses multiple doped semiconductors that act either as a dielectric or a metal depending on the frequency of the incident radiation. In both approaches (i.e., [17
–19]), the in-plane momentum of the SP (${k}_{\mathrm{sp}}$ or ${k}_{x}$) determines the extent to which the SP’s EM fields are distributed in different layers of the structure. This evolving distribution of the EM fields with respect to ${k}_{\mathrm{sp}}$ imparts upon the SPs different energies as a function of ${k}_{\mathrm{sp}}$, yielding a DC whose shape is highly customizable using the techniques and equations given in [17
–19] and in this work.

In this paper, it will be shown that multilayer structures can yield DCs that exhibit customizable stair-step shapes. Our major interest for customizable stair-step DCs is controlling the spatial profiles of the SP fields in the multilayers. It is important to note that SP DCs for a wide variety of structures have been studied in the past, which exhibit negative and/or zero group velocity by controlling the DCs [10
–13], but in ways that are largely *not* customizable, at least when compared to the approaches described in [17
–19] and in this work. It is also important to note that SPs in doped semiconductors have been researched extensively in the last five years and are described in well-cited works by Boltasseva and Atwater [20], Ginn *et al.* [21], and Law *et al.* [22]. Yet this work is different in an important way from these prior works, which used doped semiconductors, in that the electromagnetic field distributions of the SPs are more important than the simple dependence of ${\omega}_{p}$ on free carrier concentration. Shekhar and Jacob [23] also studied the strong coupling in hyperbolic metamaterials consisting of a highly doped semiconductor acting as a metallic layer and an active multiple quantum well dielectric slab. However, in our work, we use only doped semiconductor materials to control dispersion relation and the SP sorting in them. In fact, the use of doped semiconductors to control ${\u03f5}_{m}$ in each layer to control DCs is only an example of the set of materials; another set of materials can be using multiple layers of different metal alloys with different values of ${\u03f5}_{m}$.

This paper builds upon the concept, introduced in [19], of using multiple layers of doped semiconductors. In [19], an equation was derived that predicted the *beginning* of each step in the SP DC, namely, minimum value of ${k}_{x}$ at which the EM fields of the SPs would be concentrated at a particular interface in the multilayer structure. This, in turn, yielded predicable and customizable values for ${\omega}_{\mathrm{sp}}$ at particulars of values of ${k}_{x}$, as described in [19]. This work will first describe an analytic analysis of the SPs in the structure, yielding an equation that predicts the *maximum* value of $\mathfrak{R}({k}_{x})$. Next, pulsed excitation and CW excitation of the structure are performed using finite difference time domain (FDTD) and standard transfer matrix method (TMM), respectively, and the behavior of the generated SPs is studied. This modeling shows the following: the stair-step shape of the SP DC, the existence of the end of each step in the SP DC for *both* the pulsed and CW excitation conditions, and the spatial sorting of the EM fields associated with the SP of a particular ${k}_{x}$ value.

## 2. ANALYTIC METHOD

The multilayer structures considered in this work are composed of simple metal, semiconductor, or dielectric layers. All of these types of materials can be realized with doped semiconductor materials that each has a Drude dielectric response given by

where $\omega $ is the frequency, ${\omega}_{p}$ is the plasma frequency, $\gamma $ is the damping constant, and ${\u03f5}_{\infty}$ is the dielectric permittivity as $\omega \to \infty $, and ${m}^{*}$ is the electron effective mass. For $\omega <{\omega}_{p}$, $\u03f5$ is less than zero, and for $\omega >{\omega}_{p}$, $\u03f5$ is greater than zero; thus, the frequency ${\omega}_{p}$ establishes the frequency ranges of dielectric and metallic behavior of each layer. For $N$ differently doped layers, we have $N$ different ${\omega}_{p}$’s. Arranging these from lowest to highest, we obtain ${\omega}_{p,1},{\omega}_{p,2}\dots {\omega}_{p,N}$ with ${\omega}_{p,M}<{\omega}_{p,M+1}$. The number of metallic and dielectric layers can then be readily determined from the value of the driving frequency $\omega $. If $\omega <{\omega}_{p,1}$, then all layers are metallic; if $\omega >{\omega}_{p,N}$, then all layers are dielectric; and if ${\omega}_{p,M}<\omega <{\omega}_{p,M+1}$, then there are $M$ dielectric layers and $N-M$ metallic layers. These transitions split the frequency band of interest into $N+1$ bands where different EM phenomena are possible. In the first band (i.e., $\omega <{\omega}_{p,1}$), all layers are metallic and two SP modes occur, one at the topmost surface of the structure and the second at bottom-most surface. In the last band (i.e., $\omega >{\omega}_{p,N}$), all layers are dielectric and no surface plasmon modes exist, although dielectric waveguide modes may be present in the system. At most, for $N$ layers, there are $N+1$ (if $N$ is odd) or $N$ (if $N$ is even) SP modes if the layers alternate between metal-like and dielectric-like throughout the structure. In general, the number of possible surface modes at a given frequency depends on the number of interfaces between metallic/dielectric interfaces because an SP can only be sustained at such an interface. The number of SP modes will depend on the material layer configuration and the plasma frequency of every layer in the system.The frequencies of the possible SP modes can be determined in the high wavevector (momentum) limit. In this limit, where ${k}_{\mathrm{sp}}\to \infty $, the SP is bound closely to the metal/dielectric interface with the fields only experiencing the metal layer at the interface and dielectric directly adjacent to the interface. The SP frequency in this case is obtained from the dispersion relation of a surface plasmon on a single metal dielectric interface, given by

As the frequency of the SP (of the type with real ${\omega}_{\mathrm{sp}}$ and complex ${k}_{\mathrm{sp}}$) approaches the value given by Eq. (4), for any real material with loss, even for a single interface metal/dielectric system, as shown in Fig. 1, the SP DC reaches a maximum momentum and then curves back around (i.e., back-bends) to ever-decreasing momentum values as its frequency increases beyond ${\omega}_{\mathrm{sp}}$. The same property occurs for SPs at each interface in the multilayer structure, and Eq. (5) can be used to calculate the *maximum*
${k}_{x}$ possible (i.e., ${k}_{x,\mathrm{max}}$) for each SP associated with a particular metal/dielectric interface of the system. Thus, the steps in the DC of the system, when a step ceases to exist, occur at the ${k}_{x}$ values predicted by

In this work, a particular structure is studied (Fig. 2) to illustrate the concept of SP dispersion engineering; however, this concept can be used in a wide variety of structures with different layers and layer properties. Figure 2 shows a structure with three doped InSb layers, Layers 1–3, which are either metallic or dielectric in a manner depending on the frequency of the EM field. These three layers have the same Drude parameters, ${\u03f5}_{\infty}=15.7$ and $\gamma =1.52\times {10}^{10}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, and layer-specific doping concentrations of ${n}_{1}=1.9\times {10}^{15}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-3}$, ${n}_{2}=2.6\times {10}^{15}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-3}$, and ${n}_{3}=4.3\times {10}^{15}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{-3}$. The thickness of the layers are ${t}_{1}=2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, ${t}_{2}=2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, and ${t}_{3}=20\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, respectively. The three different doping concentrations result in plasma frequencies of ${\omega}_{p,1}=2.01\times {10}^{13}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{s}$, ${\omega}_{p,2}=2.35\times {10}^{13}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{s}$, and ${\omega}_{p,3}=3.02\times {10}^{13}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{s}$. These values are achievable with doped InSb where the mobility of charges ${\mu}_{e}$ of InSb is $1.5\times {10}^{5}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cm}}^{2}/\mathrm{Vs}$ at low temperatures [24]. A semi-infinitely thick dielectric superstrate of dielectric permittivity ${\u03f5}_{\mathrm{sup}}=4$ resides above the semiconductor layers. Because this structure has layers of monotonically increasing doping concentration for layers progressively further from the superstrate, it is expected that only a single metal/dielectric interface is present at each frequency below ${\omega}_{p,3}$ (no metal/dielectric interface is present for $\omega >{\omega}_{p,3}$); thus, only one true SP mode (of real ${\omega}_{\mathrm{sp}}$) exists for each frequency.

Present in this structure (Fig. 2) are three interfaces: Interface (1): Superstrate/Layer 1 with the superstrate supplying ${\u03f5}_{\mathrm{sup}}$ and Layer 1 supplying ${\u03f5}_{m}$ in Eq. (3), and similarly for Interface (2): Layer 1/Layer 2, and Interface (3): Layer 2/Layer 3. Thus, one expects and indeed sees three separate SP DCs [Fig. 3(a)]. Also in Fig. 3(a) are vertical lines calculated using Eq. (5) that denote the expected maximum value of ${k}_{x}$ for each stair-step level of the SP DC; these are the maximum ${k}_{x}$ values observed for each level in the combined SP DC for the entire structure, but only for the case when CW excitations are used. These aspects of the SP DC, namely, the values for ${\omega}_{\mathrm{sp}}$ and ${k}_{x,\mathrm{max}}$ for each level, will be compared to what is observed when calculating the SP DC using the standard FDTD and TMM method. What is built into the analytic model, but not directly modeled or calculated by it, is the spatial sorting of the SP modes that are expected in the structure; verification and detailed analysis of such sorting require the FDTD or TMM, as described below.

## 3. MODELING

For the FDTD simulations, the modeling tool Lumerical FDTD is used to plot the SP DCs and demonstrate the spatial SP sorting. In this model, multiple TM polarized dipole excitation is distributed randomly throughout the structure and with random orientations. Bloch boundary conditions are used for the $x$ direction that determine the ${k}_{x}$ value of the excitation, and perfectly matched layer (PML) boundaries are used in the $y$ direction. The dipoles excitations are initiated at time $t=0$, and a certain amount of time (the apodization time) is allowed to pass so that all nonsurface-confined EM fields escape the structure. The apodization procedure eliminates transient effects by applying a Gaussian-shaped windowing function about a particular time after the excitation; such a process filters out the short-lived transients in $E(t)$. The remaining fields are surface-confined, or at least structure-confined fields, and are time-dependent (i.e., $E(t)$). The resulting time domain data is then Fourier transformed, resulting in the frequency domain field $E(\omega )$. The frequency at which $E(\omega )$ has the highest amplitude for each ${k}_{x}$ value is identified as the frequencies of SPs; these frequencies are plotted in Fig. 3(b), yielding the SP DCs that the structure supports. Since Fig. 3(b) shows the Fourier transform of the $E$ fields from a portion of the time signal, the unit is arbitrary. The fields profiles $\mathfrak{R}({E}_{y})$ are plotted in Figs. 4(a)–4(c) for each value of ${\omega}_{\mathrm{sp}}$ that is associated with different levels in Fig. 3(b). Figures 4(a)–4(c) most clearly show the spatial sorting at the expected interfaces using this SP DC engineering methodology. Figure 4(a) shows that the SP with $f=0.734\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{THz}$ and $k=0.0406\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, the SP with $f=0.878\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{THz}$ and $k=0.0188\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$, the SP with $f=1.041\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{THz}$ and $k=0.0604\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=-2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\mu m}$. It is seen that the SP frequencies predicted by Eq. (4) match those obtained by FDTD.

The TMM is used to plot the SP DCs in Fig. 3(c) for the frequency and wavevector values [25]. In Fig. 3(c), high values of a buildup of EM field (in white) indicate the high field intensities of SPs with an arbitrary unit. In the TMM, the EM fields are expressed as a linear combination of plane waves with expansion coefficients that are determined by the initial conditions (incident beams), and no other approximations or assumptions are made except the boundary conditions, i.e., continuity of the components of $E$ and $H$ tangential to the interfaces in the structure. Specifically, the fields are expressed as ${E}_{y}{e}^{i({k}_{x}x\mp {k}_{y}y)}$ and similarly for ${H}_{z}$, with $-{k}_{y}$ used for the downward $y$ propagating waves or evanescently growing (for increasing $y$) waves, or the $+{k}_{y}$ used for the upward $y$ propagating waves or exponentially decaying (for increasing $y$) waves, and with $\u03f5{k}_{o}^{2}={k}_{x}^{2}+{k}_{y}^{2}$ where $\u03f5$ is the permittivity of the material and ${k}_{o}=\omega /c$. To assess the EM properties of the material both inside and outside the light cone defined by the superstrate, ${k}_{x}$ is allowed to vary from 0 to values large enough to be well outside the light cone (i.e., ${k}_{x}>\sqrt{{\u03f5}_{\mathrm{sup}}}{k}_{o}$). For ${k}_{x}$ values greater than $\sqrt{{\u03f5}_{\mathrm{sup}}}{k}_{o}$, the input excitation beam is not a propagating mode in the superstrate but rather an evanescent beam that is exponentially growing in strength as $y$ increases. This is fine from mathematical and experimental perspectives; in fact near-field scanning optical microscopy can provide such excitation. The resulting output beam is also evanescent in manner, exponentially decaying for increasing $y$, and indicates a buildup of fields at the interface, i.e., indicative of a SP. Alternatively, one can use the TMM to determine the poles in the scattering matrix that are also indicative of SPs.

One thing to note about the difference between the FDTD and TMM approaches is the difference in excitation of the system. In the TMM approach, the structure is only excited from above; thus, potentially only a subset of the complete set of SP modes may be excited in this analysis. This is opposed to the FDTD case, where the random position of the dipole excitation throughout the structure ensures that all possible SP modes are excited and included in the complete DC of the structure.

## 4. DISCUSSION

According to Archambault *et al.* [26], the SPs excited by short pulses (e.g., pulsed dipoles or other transient sources) will be SP modes with complex values for $\omega $ and real values of ${k}_{x}$ and do not exhibit the back-bending that is necessary to produce an *end* to the steps of the SP DC. To a certain extent, this is seen when comparing Figs. 3(b) and 3(c), where the predicted frequencies of the SPs are in agreement but the *end* of the stair-steps for the SP DC generated by FDTD are much less clear when compared to the TMM results of Fig. 3(c). For the exact solution, one might use the analytic approach to calculate the beginning and end of the stair-steps of the DCs. The differences in the generated SP curves could be further explored in future works where one studies the dependence of the SP DC on the nature of the excitation.

The level of well-defined SP sorting can be controlled by how each SP field is localized, as shown in Fig. 4. For the SP with $f=0.734\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{THz}$ and $k=0.0406\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{\mu m}}^{-1}$ [Fig. 4(a)], the field concentration at the first interface (superstrate/Layer 1) is clearly evident. Importantly, though, the numerical simulation shows that the fields penetrate into the lower metal layers (these layers are metallic at $f=0.734\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{THz}$). This can be a significant effect because in Eq. (4) one assumes that the majority of the energy of the SP fields are localized in the layers immediately adjacent to one particular interface, and its validity, therefore, depends on the thicknesses of these layers and the decay lengths of the SP modes in these layers. The structure investigated here was selected with thicknesses that are appropriate, given the predicted SP decay lengths; thus, the SP’s frequencies and momenta are not significantly perturbed from their predicted values. Figure 4(b) shows the fields for the SP that occur in the second frequency band, i.e., when the top layer is a dielectric. In this case, the SP is spatially localized at the second interface (Layer 1/Layer 2). Figure 4(c) shows the situation in the third band, where the SP is strongest at the third interface (Layer 2/Layer 3), but extends slightly into dielectric-acting Layers 1 and the superstrate, thereby slightly perturbing the SP frequency predicted by Eq. (4). It is seen in these plots that the structure has spatially sorted the SPs as predicted.

The behavior of the SPs also depends strongly on the damping constant. For example, the dispersion plots for a structure with an increased damping of $\gamma =1.52\times {10}^{11}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ are shown in Fig. 5 to demonstrate how to control the *complex shapes* of SP DCs. One can see that the SP associated with the first interface (superstrate/Layer 1) has a typical SP DC. In Layer 2, however, i.e., in the second frequency band between the two SP frequencies, there is an SP mode that rapidly overdamps, showing the negative slope of the $\omega -k$ dispersion seen in overdamped systems. In Layers 1 and 3, SPs have a momentum larger than the light line, whereas the SP in Layer 2 has a momentum smaller than the momentum of the light line at this frequency. Thus, the SP in Layer 2 will have strong radiative losses if the fields penetrate through Layer 1 into the superstrate. The choice of material thickness needed to prevent this can be determined from the transverse decay length predicted using Eq. (5). While the thickness of each layer in this work is chosen with the appropriate SP decay lengths to demonstrate the spatial SP sorting in multilayer structures, Blazek *et al.* [27] showed that a thin transition layer (TL) with a linear change in carrier concentration also supports one branch dispersion curve regardless of the TL thickness. Thus, any SP confined layer in a multilayer system contributes to the shape of the DCs. It can be also noted that one might expect a bonding–antibonding type of effect for the surface modes at each interface that influence the form of dispersion relation, which is highly sensitive to the values of those thicknesses [28].

## 5. CONCLUSION

We have demonstrated an analytical approach to predict and engineer the dispersion curves of stacked doped semiconductors and verified the analytical results with FDTD and TTM modeling. Particularly, a system with three spatial and frequency-dependent dispersion profiles was analyzed, but the method is generalizable to a wide range of multilayer dielectric/metal structures. The validity of the approach with regard to material thickness and loss was discussed. Time domain studies were performed to study the creation, decay the SPs, and confirm the assumptions of the analytic model. It was demonstrated that the spatial sorting of plasmons is possible for different frequencies at different layers by controlling the SP DCs. In future works, one can further investigate these effects and how the inclusion of dielectric layers allows for additional modes, including coupled metal-insulator-metal modes and dielectric waveguide modes. These frequencies, momenta, and spatial concentration of the additional modes, in turn, can be engineered using the methods described in this paper. Our model works for most of the real devices where the intrinsic plasma frequency (IPF) is smaller than the scattering frequency [24,29]. The low loss plasmonic materials (PMs), especially in the doped oxide materials in the near-IR and visible ranges, can be also found from recently studied literature [20,30,31]. Applications include tunable broadband perfect absorbers; band-stop filters; solar cells; slow light and buffering in telecommunication systems; stopping, trapping, and fast light effect; and optical analogs of gravitational black and white holes.

## FUNDING INFORMATION

National Science Foundation Industry/University Cooperative Research Center for Metamaterials (IIP-1068028).

## REFERENCES

**1. **R. H. Ritchie, “Plasma losses by fast electrons in thin films,” Phys. Rev. **106**, 874–881 (1957). [CrossRef]

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

**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).

**4. **D. Crouse, “Surface plasmon effects in metal-semiconductor-metal photodetectors,” Proc. SPIE **5594**, 45–56 (2004).

**5. **D. Crouse and Y. Lo, “Nonsteady-state surface plasmons in periodically patterned structures,” J. Appl. Phys. **95**, 4163 (2004). [CrossRef]

**6. **K. Li, X. Li, M. Stockman, and D. Bergman, “Surface plasmon amplification by stimulated emission in nanolenses,” Phys. Rev. B **71**, 115409 (2005).

**7. **D. Crouse and R. Solomon, “Numerical modeling of surface plasmon enhanced silicon on insulator avalanche photodiodes,” Solid-State Electron. **49**, 1697–1701 (2005). [CrossRef]

**8. **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]

**9. **I. Mandel, E. Lansey, J. N. Gollub, C. H. Sarantos, R. Akhmechet, A. B. Golovin, and D. Crouse, “Photon sorting in the near field using subwavelength cavity arrays in the near-infrared,” Appl. Phys. Lett. **103**, 251116 (2013). [CrossRef]

**10. **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]

**11. **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]

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

**13. **C. Caloz, “Metamaterial dispersion engineering concepts and applications,” Proc. IEEE **99**, 1711–1719 (2011). [CrossRef]

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

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

**16. **E. Feigenbaum, N. Kaminski, and M. Orenstein, “Negative dispersion: a backward wave or fast light? Nanoplasmonic examples,” Opt. Express **17**, 18934–18939 (2009).

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

**18. **A. Karalis, J. D. Joannopoulos, and M. Soljačić, “Plasmonic-dielectric systems for high-order dispersionless slow or stopped subwavelength light,” Phys. Rev. Lett. **103**, 043906 (2009). [CrossRef]

**19. **I. M. Mandel, I. Bendoym, Y. U. Jung, A. B. Golovin, and D. T. Crouse, “Dispersion engineering of surface plasmons,” Opt. Express **21**, 31883–31893 (2013). [CrossRef]

**20. **A. Boltasseva and H. A. Atwater, “Low-loss plasmonic metamaterials,” Science **331**, 290–291 (2011). [CrossRef]

**21. **J. Ginn, R. L. Jarecki Jr., E. A. Shaner, and P. S. Davids, “Infrared plasmons on heavily-doped silicon,” J. Appl. Phys. **110**, 043110 (2011). [CrossRef]

**22. **S. Law, L. Yu, and D. Wasserman, “Epitaxial growth of engineered metals for mid-infrared plasmonics,” J. Vac. Sci. Technol. B **31**, 03C121 (2013). [CrossRef]

**23. **P. Shekhar and Z. Jacob, “Strong coupling in hyperbolic metamaterials,” Phys. Rev. B **90**, 045313 (2014).

**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. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic Theory of Propagation, Interference, and Diffraction of Light* (Cambridge University, 1999).

**26. **A. Archambault, T. V. Teperik, F. Marquier, and J. J. Greffet, “Surface plasmon Fourier optics,” Phys. Rev. B **79**, 195414 (2009).

**27. **D. Blazek, M. Cada, and J. Pistora, “Surface plasmon polaritons at linearly graded semiconductor interfaces,” Opt. Express **23**, 6264–6276 (2015). [CrossRef]

**28. **E. N. Economou, “Surface plasmons in thin films,” Phys. Rev. **182**, 539–554 (1969). [CrossRef]

**29. **M. Cada, D. Blazek, J. Pistora, K. Postava, and P. Siroky, “Theoretical and experimental study of plasmonic effects in heavily doped gallium arsenide and indium phosphide,” Opt. Mater. Express **5**, 340–352 (2015). [CrossRef]

**30. **G. Naik, J. Kim, and A. Boltasseva, “Oxides and nitrides as alternative plasmonic materials in the optical range [Invited],” Opt. Mater. Express **1**, 1090–1099 (2011). [CrossRef]

**31. **A. K. Pradhan, R. M. Mundle, K. Santiago, J. R. Skuza, B. Xiao, K. D. Song, M. Bahoura, R. Cheaito, and P. E. Hopkins, “Extreme tunability in aluminum doped Zinc Oxide plasmonic materials for near-infrared applications,” Sci. Rep. **4**, 6415 (2014). [CrossRef]