## Abstract

New results are reported on investigation of dispersion curves for surface plasmon polaritons (SPPs) at an inhomogenously doped semiconductor/dielectric interface whereby the dielectric is represented by the same undoped semiconductor. The doped semiconductor is described by its frequency-dependent permittivity that varies with the depth. It is shown that a transition layer (TL) with a linear change in carrier concentration supports one branch dispersion curve regardless of the TL thickness. The obtained dispersion curves reach a maximum at a finite frequency depending on the TL thickness, and subsequently asymptotically approach the zero frequency in the shortwave limit. Therefore two surface plasmon modes are supported at a given frequency: a long-wave mode with a positive group velocity and a short-wave mode with a negative group velocity. A condition of a zero group velocity can be satisfied by tuning the TL layer. It is shown that the conventional dispersion relation for SPPs at a TL with a zero thickness is an asymptotic solution, and the convergence of real dispersion curves is point-wise instead of an expected uniform convergence.

© 2015 Optical Society of America

## 1. Introduction

SPPs are transverse magnetic (TM) surface electromagnetic waves that propagate typically along a metal–dielectric interface [1–3]. The model of a surface plasmon (SP) was introduced [4] after the discovery of bulk plasmons. SPs have been observed at many interface structures [5] including the semiconductor ones [6], and various applications have been proposed and investigated [7].

The structure of SPPs that propagates along an interface between metal and dielectric materials is well known. The excited SPPs waves consist of associated evanescent fields that penetrate both the dielectric and the metal. The field is concentrated mainly near the surface [8]. Due to its concentration, the SPPs are very sensitive to the surface properties such as surface roughness or inhomogeneity of the permittivity. The properties of a transition layer (TL) between the two media should thus be considered in a model.

The problem of the TL in optics is quite an old one, but the attention of researchers has been concentrated on the problem of the reflection and transmission of light [2]. For the metals, a TL can be introduced by depositing a thin metal layer on a substrate. The dispersion curve is split into two branches with a frequency gap appearing between them if the plasma frequency of the deposited layer is much smaller than that of the plasma frequency of the substrate [9]. The transition layer can give rise to an SP dispersion curve that has a section with negative slopes, i.e. negative group velocities [2,9]. The predicted phenomena have been experimentally observed [10–12]. The existence of negative and zero group velocities have thus been known since at least the seventies.

It is worth noting that today a great attention is devoted to the negative index metamaterials (NIM). The unusual properties of NIM are most prominently revealed at an interface between positive and negative index materials [13]. Only recently it has been realized that inhomogeneous metamaterials enable such fascinating functionalities as cloaking or wave concentrators [14,15]. Moreover, there exists an analytical solution based on confluent hypergeometric functions for gradually changing material properties (permittivity and permeability simultaneously) from the positive to the negative values if the profile is linear or exponential [16]. The phenomenon of localization and resonant enhancement of the field in the vicinity of the point where material properties change signs was confirmed [13,16].

Moreover, it has been shown that the heterostructure waveguide mode is an extraordinary entity where the Poynting vector points towards the + z direction in the material with positive permittivity, while it points towards –z direction in the NIM [17]. The resulting energy velocity (which equals to the group velocity in the lossless structure) may be both parallel or antiparallel to the SPPs wave vector depending on the overall energy flow. This implies the regime with zero or near-zero group velocities when the positive and negative energy flow is balanced [17]. Our investigation of a TL in semiconductors was motivated by the surface effects that lead to the existence of depletion, accumulation and inversion layers [18]. In contrast to the metals, the permittivity of a semiconductor can be made a continuous function of spatial coordinates. An analytical solution does not exist in these situations, and series solutions based on the Frobenius method were used [18–20]. A complicated structure of branches with both positive and negative slopes was discovered [20]. In addition to solutions for complex permittivity profiles, we have found only one article dealing with the geometry where the permittivity is linearly graded from its positive value in the dielectric to the negative bulk value in the substrate [19]. It was concluded there that no bound modes exist in the idealized case of lossless materials, which is a different result than the one reported here.

The studied structure is shown in Fig. 1. A TL with a linearly graded doping concentration is situated between the doped and undoped semi-infinite layers. All three layers are made of the same semiconductor. At a sufficiently low frequency the permitivity of the doped semiconductor becomes negative enough for an SP confined to the TL to exist.

Our investigation is thus devoted to the properties of SPPs at an inhomogeneously doped semiconductor/dielectric interface. We assume that the linear grading of the permittivity is caused by doping profile in the semiconductor ranging from zero free carrier concentration in the superstrate and continuously rising to the bulk concentration in the substrate. The local relation between the displacement vector, *D*, and the electric field vector, *E*, has been employed, where the permittivity depends only on the local plasma frequency. The validity of this assumption is discussed elsewhere [18]. In addition, the dielectric permittivity tensor is assumed, for simplicity, to be isotropic and characterized by a single variable *ε*. This is not generally the case in, for example, III/V semiconductors, however, the assumption here is satisfactory for the purposes of our analysis.

The linear dependence of the charge density on the position has been chosen since it is the most fundamental problem while it is still simple and treatable analytically. More complex profiles, which subject is certainly interesting and worth pursuing, would complicate the study to the point that an analytical asymptotic solution would not be attainable and thus the obtained differential equation would possess a more complex singularity. The analysis would become burdened with a number of new parameters and the insightful and rather elegant solution would be lost

What has transpired while considering a more general situation was that the interface position, i.e. the plane where the permittivity is equal to zero, depends on the frequency. For a zero frequency this point is closest to the dielectric (or even at the edge with the transition layer); it then moves with increasing frequency towards the lower edge of the transition layer, reaching it at the plasma frequency. If the permittivity is not a linear function of the position, the character of the interface depends on frequency, which complicates the analysis significantly and, again, the insightful analytical solution is lost.

We have calculated SPPs dispersion curves for different TL. We have found only one nonradiative branch for each TL thickness. This branch follows the light-line at small wavenumbers, then reaches a maximum frequency that depends on the TL thickness, and finally decreases and approaches the zero frequency in the limit of large wavenumbers. The obtained dispersion curves converge point-wise to the usual dispersion curve as the TL thickness approaches zero. The tangential magnetic field *H*_{y} and the perpendicular displacement field *D*_{x} have a rounded peak, i.e. with a continuous derivative, at a position *x*_{0} where the permittivity crosses the zero value, while the tangential displacement field *D*_{z}*(x _{0})* is zero. The confinement of the SPPs increases with rising wavenumber and finally the whole field is localized within the TL.

One can also note that the nature of our results is, indeed, an opposite to an engineering approach that calls for a model with several free parameters that enable one to tailor the dispersion curve based on an application need. In [21,22], for example, such an engineering problem is addressed by creating a multilayered dielectric structure deposited on the conductive substrate. By modifying properties of individual layers the dispersion curve can be properly adjusted. It should be noted that for a practical use of presented phenomena, especially the frozen light mode, a more complex structure might have to be proposed, which would take into account the real characteristics of the used semiconductors; that would also extend the modal index bandwidth (MIB). A large MIB is neccesary to be able to couple to the slow light mode and also to enable the efficient energy accumulation without spatial spreading [17].

## 2. The mathematical model

The macroscopic permittivity of a substrate semiconductor with free carriers can be described by the well-known Drude model [23]:

^{−12}– 10

^{−13}s at the room temperature irrespective of the electron mobility because the damping is inversely proportional to the product of mobility and effective mass. The second approximation in our model is an assumption of a zero plasma frequency ${\omega}_{p}$ in undoped semiconductors due to the low intrinsic density. The intrinsic density is a product of the density of states and Boltzman exponential factor containing the bandgap energy and temperature. As a result the intrinsic plasma frequency (IPF) at room temperature lays in a wide range of 10

^{8}s

^{−1}– 10

^{13}s

^{−1}. It may be concluded that our model will work for most of the semiconductors where the IPF is smaller than the scattering frequency. Specifically for GaAs and InP, for example, we measured relaxation times to be less than 100 fs [24] (IPF are 10

^{8}s

^{−1}and 10

^{9}s

^{−1}, respectively) while the operation frequencies are almost two orders of magnitude higher. Thus both the damping term and IPF can be neglected and the dispersion equation derived describes the problem adequately. Both the damping term and IPF can be decreased rapidly by cooling the semiconductor to the low temperatures thus substantially improving the reliability of (1) even further.

A variation of the permittivity with the semiconductor depth (*x-*axis) is assumed as:

_{${x}_{N}$}is the TL thickness. Figure 1 shows the permittivity’s depth variation for a given frequency. Note that the permittivity becomes negative at a position

*x*

_{0}, which is required for an SPP to exist, depending on the plasma frequency and the depth. The dielectric, in this case the same yet undoped semiconductor, extends into

*x > 0*.

Since SPPs are TM polarised, only TM modes propagating in the *z-*direction are studied here. The spatial field distribution can be described by complex field amplitudes *H*_{y}(*x*), *D*_{x}(*x*) and *D*_{z}(*x*), each of them multiplied by the time/propagation term exp[j (*ωt* - *βz*)]. Under these assumptions the set of Maxwell’s equations yields three independent equations, written in the frequency domain, as:

It is convenient to rescale spatial dimensions by a factor $\omega /c$ to normalize the expressions for further numerical analysis and interpretations. Let $v$ be the dimensionless rescaled spatial coordinate $v=\left(\omega /c\right)\text{\hspace{0.17em}}x$ and ${n}_{eff}=\left(c/\omega \right)\text{\hspace{0.17em}}\beta $ be the modal index [25]. Equation (4) assumes the form:

Equation (4) needs to be solved in order to find a field distribution in the TL. An analytical solution is not available, however a series solution expanded about the singular point${v}_{0},$ based on the Frobenius method is known [18–20]. The regular singularity appears at ${v}_{0}$ when $\epsilon \left({v}_{0}\right)=0.$ According to Eq. (2), the permittivity around the singularity is a linear function $\epsilon \left(v\right)=A\left(v-{v}_{0}\right),$where, obviously:

*O(.)*represents a mathematical remainder of the order of (.).

It can be concluded that the tangential magnetic field *H _{y}* and the perpendicular displacement field

*D*have a rounded peaks at the position ν

_{x}_{0}while the tangential displacement field ${D}_{Z}\left({\nu}_{0}\right)$is zero. On the other hand, the normal component,

*E*, has a singularity of a ${\left(v-{v}_{0}\right)}^{-1}$type while the tangential component

_{x}*E*exhibits a logarithmic singularity. This singular behaviour of the intensity of the electric field is an artefact that arises from the simultaneous use of the local-permittivity approximation together with neglecting losses. If a correct non-local calculation were made, these infinite fields would not appear [17]. An introduction of small losses i.e. the very small imaginary part of the permittivity, ${\epsilon}_{i},$ modifies the argument of the logarithm function in Eqs. (10a) and (10b) as:

_{z}It can be noted that the singularity in our numerical solution was removed by introducing a small imaginary part to the Drude permittivity, independent of the frequency and the position. The resulting space distribution of the magnetic field is represented by the real part of the solution. The validity of this approach was verified by variatying the introduced imaginary part while the field solution did not display any observable effects. We do not claim this to be a physical fact, although it does support our assumption of neglecting damping; it also proves the stability of our procedure that employed a standard numerical method (built-in in Mathematica).

The presence of singularities in the intensity of the electric field precludes a possibility to calculate the energy velocity v_{e} that is given by the ratio of the time averaged Poynting vector *S*_{z} and the average energy density *W* integrated along the *x*-direction [26]:

*W*cannot be removed in the framework of presented model. On the other hand, the singularity in the Poynting vector:

## 3. Results

Conventional numerical methods have been used to solve Eq. (6) with applying boundary conditions in Eqs. (8a) and (8b). For each desired wavenumber the following process was used. First, a frequency was chosen, the condition in Eq. (8a) was applied, then the field distribution was calculated and finally the condition in Eq. (8b) was checked. The process continued until a frequency was found that satisfied both boundary conditions for a calculated field. Exactly one solution exists for each wavenumber and TL thickness. It should be noted that a numerical method based on the bisection method whereby an algorithm repeatedly bisects an interval and then selects a subinterval in which a root must lie for further processing, would find second solution which is obviously an artefact. This is a solution that does not represent a physical reality and may be termed as a “virtual” non-physical solution. It coincides with initial conditions that cause the magnetic field ${H}_{y}(0)=0$ and a discontinuity on the left-hand side of the condition (8b).

The SPPs dispersion curves were in the focus of our interest. A wide range of the TL thicknesses was investigated. Although the silicon parameters were used in numerical calculations, an effort was made to calculate dispersion curves as general as possible by introducing appropriate dimensionless parameters. The frequency was normalized to the plasma frequency, $\tilde{\omega}=\omega /{\omega}_{p},$ which made the dispersion curves conveniently independent of the bulk doping level. The propagation constant was normalized as $\tilde{\beta}=\beta /{\beta}_{p},$ with

being a value of $\beta $where the asymptotes of the light line and the conventional maximum SP’s frequency intersect. It can also be noted that the well-known dispersion relation for SPPs with zero thickness of TL (conventional surface plasmon with${\beta}_{0}$) becomes independent of the material permittivity and the plasma frequency when the above normalization is used, thus:Dispersion curves of SPPs are plotted in Fig. 2 for different widths of the TL, i.e. ${x}_{N}$. One can see that the whole curve monotonically decreases in frequency with an increase in ${x}_{N}$. The phase velocity of SPPs is decreasing with increasing the TL thickness. A very important discovery is the point-wise character of the convergence of obtained dispersion curves when compared to the conventional ones, as the TL width,${x}_{N},$ approaches zero. There is a maximum on each and every dispersion curve except for the case of ${x}_{N}=0$; the SP’s frequency approaches asymptotically zero in the limit $\beta \to \infty .$ The position of the maximum can be controlled by the thickness of TL, ${x}_{N}$. One observes that the maximum frequency of the SPPs is reduced and a stationary surface plasmon (frozen mode), i.e. a plasmon with a zero group velocity, appears. The long wave plasmons are only slightly affected while their group velocities are determined by the lattice permittivity ${\epsilon}_{\infty}.$ On the other hand, the shortwave SPPs propagate with a negative group velocity. Their frequency and phase velocities are decreasing as the wavenumber increases, and the electromagnetic field of the SPPs become localised in the TL.

Each point $\left(\tilde{\beta},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\tilde{\omega}\right)$in Fig. 2 is uniquely associated with some TL thickness _{${x}_{N}$}. Thanks to the normalization, the dispersion curves are valid for all semiconductors (not considering the dispersion of the lattice permittivity ${\epsilon}_{\infty}$); however, the width of the TL must be renormalized with respect to the lattice permittivity of a specific material chosen, ${\epsilon}_{\infty}^{mat}$, as:

The attractive feature of the SPPs being able to become spatially stationary deserves detailed attention. Conditions necessary for the existence of such a surface plasmon clearly transpire from its dispersion properties. For every selected frequency, there exists a corresponding TL thickness whereby a stationary surface plasmon is supported by the corresponding structure. Table 1. summarizes the required values for such conditions for the case of silicon with a lattice permittivity ${\epsilon}_{\infty}=11.68.$

The values in Table 1 are plotted in Fig. 3 where the required (un-normalized) width of the transition region, ${x}_{N}$, that supports the new spatially stationary (standing) plasmon is shown as a function of the normalized frequency, $\omega /{\omega}_{p}$. The TL thickness decreases monotonically over several orders of magnitude as the normalized frequency increases. The SPPs propagating at a TL thinner than $0.1\text{\hspace{0.17em}}\mu m$ display almost the same properties as the conventional ones propagating along a well-defined interface between two different media. While the second row of Table 1. enables one to adjust the right TL thickness for obtaining stationary SPPs with a desired frequency, the third row allows one to predict the wavenumber of resulting SPPs. It is very interesting that the stationary SPPs have approximately three times shorter wavelength compared to the conventional SPPs with the same frequency.

Our most important finding is that the dispersion curve exists for any thickness of the TL and that its asymptotic frequency is zero. It is thus obvious that it must have a maximum for a finite propagation constant and thus, as in the agreement with other authors’ findings with layered structures [21,22], it possesses a negative as well as a zero group velocity.

The maximum of the magnetic field is at the interface between the TL and the dielectric at the zero frequency, while it shifts towards the doped region with increasing the frequency. In the case of short-wave plasmons, the frequency decreases and the magnetic field maximum shifts back towards the TL. The field penetration depth decreases faster than the speed of the maximum shift, therefore the plasmon energy concentrates more and more in the TLand the confinement increases.

Some typical field distributions are shown in the Figs. 4–7. The general field distribution of all SPPs propagating on a surface with a graded TL is very similar. The magnetic fields decrease rapidly on both sides and contrary to the conventional SPPs, they have rounded peaks. Depending on the SPP’s wavelength and the TL thickness, the central part of the SPPs may be only an interpolation of the field exponentially decreasing on the sides in the case of small wavenumbers and thin TL (Fig. 4), while it becomes the most important part of the short-wave SPPs with a negative group velocity. The field distributions of SPPs are shown with both small (Fig. 4) and large wavenumbers (Fig. 6) propagating at a TL with a thickness of $1/{\beta}_{p}$. If the thickness of the TL is greater than $1/{\beta}_{p}$, the field becomes concentrated near the boundary between the TL and the undoped dielectric (Fig. 7). It does not penetrate the doped substrate in the case of a very thick TL.

An important property of surface plasmons is the confinement of energy. It is generally clear; for all surface plasmons, the confinement is determined by the horizontal distance of their point on the dispersion curve from the light line. The resulting penetration depth decreases as the plasmon’s wavelength decreases and thus the plasmon is more tied to the interface. In our case when the interface is not homogeneous, the dispersion curve lies below the conventional one so the field penetration decreases with the increasing wavenumber faster than in the case of a homogeneous interface. Thus in our case, the plasmon confinement is higher than that of the ones at a homogenous interface.

It has been shown that one stationary SPP exists for each thickness of TL. However, these modes differ in their ability to couple to an external source of light. A large MIB is neccesary to be able to couple practically to the slow light mode and also to enable the efficient energy accumulation without spatial spreading [17]. The second derivative of the dispersion curve at the maximum, which equals the first derivation of the group velocity, can be used as a criterium for MIB. The graph of group velocities versus propagation constant is presented in Fig. 8. The slopes of curves that are reciprocal to the MIB are increasing with rising the TL thickness. It may be concluded that a narrow TL supports stationary SPPs with better MIB, but at very short wavelengths. Some compromise would need to be found for practical applications. As can be seen the SPPs with a negative group velocity possess very low dispersion.

## 4. Conclusions

Details of properties of SPPs propagating in doped semiconductor structures with a continuous change in the permittivity are presented. In contrast to conventional SPPs, the magnetic field and its first spatial derivative are continuous functions. The maximum of the magnetic field is located at a plane with the zero permittivity.

The presented model enables one to determine the electric field of the new surface plasmon. Unlike the continuous magnetic field, there appears a logarithmic singularity in the tangential component and a singularity of the ${\left(x-{x}_{0}\right)}^{-1}$ type in the perpendicular component of the electric field. It is also possible to calculate the volume density of the electric charge in the transition region that has a singularity of the ${\left(x-{x}_{0}\right)}^{-2}$ type. Both singularities can be eliminated by considering material’s damping, i.e. introducing a small imaginary part to the permittivity.

The most important discovery is the effect of the TL with linearly graded permittivity on the dispersion behavior of the SPPs. There is exactly one branch regardless of the thickness of the TL. A maximum appears on the dispersion curve with its position being determined by the width of the TL. The SPPs exhibit a zero group velocity at this maximum, which is obviously the highest frequency the plasmon can be supported at. At lower frequencies there thus exist two surface plasmons: the conventional long-wave forward-propagating one and the new short-wave with a negative group velocity. Finally the SP’s frequency approaches asymptotically zero in the limit $\beta \to \infty $ and its electromagnetic field becomes localised in the TL. It is shown that the conventional dispersion relation for SPPs at a TL with zero thickness is an asymptotic solution, and the convergence of real dispersion curves is of a point-wise type.

As to our knowledge, this is the first time a graded index model with a linear function of doping concentration in the TL has been used to investigate the SPPs properties in a semiconductor structure with continuous permittivity. The discovered dispersion curves are unique by the frequency asymptotically approaching the zero value in the shortwave limit. These results enable one to significantly reduce or even stop and reverse plasmon’s energy-carrying velocity for a chosen frequency by appropriate structuring the transition region or by manipulating the doping profile. Attractive potential applications such as energy transport control or information processing and storage are just a few that come to one’s mind.

The presented model has, in fact, three degrees of freedom, i.e. design/engineering parameters, namely the lattice permittivity, the doping concentration, and the transition layer thickness. We showed that our solution is of a universal validity since it can be fully normalized in all these three free parameters. That yielded the universal Fig. 2. We believe that this represents the most valuable scientific achievement as the number of degrees of freedom has been reduced to one.

Our results are very different from the Kruger’s one [19]. Kruger made an attempt to investigate the SP propagating on metal-dielectric interface that was described by a linearly graded permittivity equivalent to our definition in Eq. (2). He obtained a convex dispersion curves with apparently unlimited grow of the group velocity.

## Acknowledgments

This work was supported by Natural Sciences and Engineering Council of Canada (NSERC) and its Collaborative Research and Training Experience (CREATE) Advanced Science in Photonics and Innovative Research in Engineering (ASPIRE) Program, by IT4 Innovations Centre of Excellence project, reg. no. CZ.1.05/1.1.00/02.0070, by the Czech Republic grant (GACR) P205/11/2137, and by the New Creative Teams in Priorities of Scientific Research grant CZ1.07/2.3.00/30.0055.

## References and links

**1. **A. D. Boardman, *Electromagnetic Surface Modes* (Wiley, 1982).

**2. **V. M. Agranovich, “Effects of the transition layer and spatial dispersion in the spectra of surface polaritons,” in *Surface Polaritons: Electromagnetic Waves at Surfaces and Interfaces, Modern Problems in Condensed Matter Sciences*, V. M. Agranovich and D. L. Mills, eds. (**1**, 1982), pp.187–236.

**3. **S. A. Maier, *Plasmonics: Fundamentals and Application,* (Springer Science and Business Media LLC, 2007).

**4. **H. Raether, *Surface Plasmons on Smooth and Rough Surfaces and on Gratings*, (Springer, 1986).

**5. **P. Berini, “Long-range surface plasmon polaritons,” Adv. Opt. Photon. **1**, 484–588 (2009).

**6. **M. Cada and J. Pistora, “Optical Plasmons in Semiconductors,*” **in Proceedings of International Symposium on Optical and Microwave Technologies*, (ISMOT, 2011), pp. 23–29.

**7. **T. W. Ebbesen, C. Genet, and S. I. Bozhevolnyi, “Surface-plasmon circuitry,” Phys. Today **61**(5), 44–50 (2008). [CrossRef]

**8. **J. Zhang, L. Zhang, and W. Xu, “Surface plasmon polaritons: physics and aplications,” J. Phys. D Appl. Phys. **45**(11), 113001 (2012). [CrossRef]

**9. **V. M. Agranovich, Y. R. Shen, R. H. Baughman, and A. A. Zakhidov, “Optical bulk and surface waves with negative refraction,” J. Lumin. **110**(4), 167–173 (2004). [CrossRef]

**10. **E. A. Vinogradov and T. A. Leskova, “Polaritons in thin films on metal surfaces,” Phys. Rep. **194**(5-6), 273–280 (1990). [CrossRef]

**11. **V. A. Yakovlev, V. G. Nazin, and G. N. Zhizhin, “The surface polariton splitting due to thin surface film LO vibrations,” Opt. Commun. **15**(2), 293–295 (1975). [CrossRef]

**12. **T. Lopez-Rios, F. Abeles, and G. Vuye, “Splitting of the Al surface plasmon dispersion curves by Ag surface layers,” J. Phys. **39**(6), 645–650 (1978). [CrossRef]

**13. **N. M. Litchinitser, A. I. Maimistov, I. R. Gabitov, R. Z. Sagdeev, and V. M. Shalaev, “Metamaterials: electromagnetic enhancement at zero-index transition,” Opt. Lett. **33**(20), 2350–2352 (2008). [CrossRef] [PubMed]

**14. **D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science **314**(5801), 977–980 (2006). [CrossRef] [PubMed]

**15. **J. Yang, M. Huang, C. Yang, Z. Xiao, and J. Peng, “Metamaterial electromagnetic concentrators with arbitrary geometries,” Opt. Express **17**(22), 19656–19661 (2009). [CrossRef] [PubMed]

**16. **P. C. Ingrey, K. I. Hopcraft, E. Jakeman, and O. E. French, “Between right and left handed media,” Opt. Commun. **282**(5), 1020–1027 (2009). [CrossRef]

**17. **S. Foteinopoulou and J. P. Vigneron, “Extended slow-light field enhancement in positive-index/negative-index heterostructures,” Phys. Rev. B **88**(19), 195144 (2013). [CrossRef]

**18. **S. L. Cunningham, A. A. Maradudin, and R. F. Wallis, “Effect of a charge layer on the surface plasmon polariton dispersion curve,” Phys. Rev. B **10**(8), 3342–3355 (1974). [CrossRef]

**19. **B. A. Kruger and J. K. S. Poon, “Optical guided waves at graded metal-dielectric interfaces,” Opt. Lett. **36**(11), 2155–2157 (2011). [CrossRef] [PubMed]

**20. **C. C. Kao and E. M. Conwell, “Surface plasmon dispersion of semiconductors with depletion or accumulation layers,” Phys. Rev. B **14**(6), 2464–2479 (1976). [CrossRef]

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

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

**23. **M. Fox, *Optical Properties of Solids*, (Oxford University Press, 2001, 2010, 2012).

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

**25. **J. C. Knight, “Photonic crystal fibres,” Nature **424**(6950), 847–851 (2003). [CrossRef] [PubMed]

**26. **V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of ε and μ,” Sov. Phys. Usp. **10**(4), 509–514 (1968). [CrossRef]