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 ${ϵ}_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 ${ϵ}_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 ${ϵ}_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

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,\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, ${ϵ}_{\infty }=15.7$ and $\gamma =1.52\times {10}^{10}{\mathrm{s}}^{-1}$, and layer-specific doping concentrations of $n_1=1.9\times {10}^{15}{\mathrm{cm}}^{-3}$, $n_2=2.6\times {10}^{15}{\mathrm{cm}}^{-3}$, and $n_3=4.3\times {10}^{15}{\mathrm{cm}}^{-3}$. The thickness of the layers are $t_1=2\mathrm{\mu m}$, $t_2=2\mathrm{\mu m}$, and $t_3=20\mathrm{\mu m}$, respectively. The three different doping concentrations result in plasma frequencies of ${\omega }_{p,1}=2.01\times {10}^{13}\mathrm{rad}/\mathrm{s}$, ${\omega }_{p,2}=2.35\times {10}^{13}\mathrm{rad}/\mathrm{s}$, and ${\omega }_{p,3}=3.02\times {10}^{13}\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{\mathrm{cm}}^2/\mathrm{Vs}$ at low temperatures [24]. A semi-infinitely thick dielectric superstrate of dielectric permittivity ${ϵ}_{\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.

 

Fig. 1. SP dispersion curve for a thin real (i.e., includes loss) silver film bounded by vacuum. SP modes with real values for $\omega$ and complex values for $k_x$ are used. Back-bending occurs with these modes in the vicinity of $\omega ={\omega }_p/\sqrt{2}$ (shown by the black horizontal line at 261 THz). The real momenta are shown as a solid line and the imaginary components are shown as a dashed line. It is seen that, as $\omega$ approaches ${\omega }_p/\sqrt{2}$, the DC curves back upon themselves, yielding a maximum $\mathfrak{R}(k_x)$ that the SP can possess. This same property exists for SPs in the multilayer system studied in this work (Fig. 2) when a time-harmonic wave is used to excite the SPs. This bending-back of the SP DC is used to predict the maximum $\mathfrak{R}(k_x)$ values of the steps in the SP DC for the multilayer structure.

Download Full Size | PDF

 

Fig. 2. Schematic of the multilayer structure. The three layers have different doping densities leading to different plasma frequencies and properties, in terms of them being metallic or dielectric in a manner at a particular frequency.

Download Full Size | PDF

Present in this structure (Fig. 2) are three interfaces: Interface (1): Superstrate/Layer 1 with the superstrate supplying ${ϵ}_{\sup }$ and Layer 1 supplying ${ϵ}_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,\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.

 

Fig. 3. (a) Analytic dispersion curves of the frequency and $k$ vector dependency for the SP modes of Layer 1 (red), Layer 2 (blue), and Layer 3 (black). SP modes with real values of $\omega$ and complex values of $k_x$ are used for this plot; thus, back-bending occurs in the vicinity of the $\omega ={\omega }_{\mathrm{sp}}$, with ${\omega }_{\mathrm{sp}}$ given by Eq. (4). The real momenta are shown as solid lines, and the imaginary components are shown as dashed lines. The horizontal lines are the frequencies by Eq. (1) at which the materials in Layers 1, 2, and 3 change from being a dielectric at lower frequencies to a metal at higher frequencies. The light line associated with the superstrate ${ϵ}_{\mathrm{sub}}$ as well as ${ϵ}_{\inf }$ is also shown. (b) SP DCs calculated using the FDTD method with an arbitrary unit. The asymptotic behavior of the DCs is shown, and the amplitude decay of SPs is observed due to Ohmic losses in the metal. (c) DCs of SPs calculated through the TTM with an arbitrary unit. Vertical lines indicate the expected limiting momentum (maximum $k_x$) from Eq. (5), i.e., blue dashed line = Interface 2; red dashed line = Interface 1; green dashed line = Interface 3.

Download Full Size | PDF

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\mathrm{THz}$ and $k=0.0406{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=2\mathrm{\mu m}$, the SP with $f=0.878\mathrm{THz}$ and $k=0.0188{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=0\mathrm{\mu m}$, the SP with $f=1.041\mathrm{THz}$ and $k=0.0604{\mathrm{\mu m}}^{-1}$ concentrates field intensities at the interface located at $x=-2\mathrm{\mu m}$. It is seen that the SP frequencies predicted by Eq. (4) match those obtained by FDTD.

 

Fig. 4. (a) $\mathfrak{R}(E_y)$ calculated from the discrete Fourier transform of the time domain study for $f=0.734\mathrm{THz}$ and $k=0.0406{\mathrm{\mu m}}^{-1}$ for the first SP mode, occurring on the interface at $y=2\mathrm{\mu m}$. Likewise, we have (b) $\mathfrak{R}(E_y)$ for $f=0.878\mathrm{THz}$ and $k=0.0188{\mathrm{\mu m}}^{-1}$ for the second SP mode, occurring on the interface at $y=0\mathrm{\mu m}$. (c) $\mathfrak{R}(E_y)$ for $f=1.041\mathrm{THz}$ and $k=0.0604{\mathrm{\mu m}}^{-1}$ for the third SP mode, occurring on the interface at $y=-2\mathrm{\mu m}$. The spatial sorting of the SP modes matches the SP frequencies predicted by Eq. (4). The initial dipoles excitation located on the interface at $y=2\mathrm{\mu m}$ is still shown. The unit of $\mathfrak{R}(E_y)$ is arbitrary.

Download Full Size | PDF

 

Fig. 5. Analytic dispersion plot of three different SP modes for $\gamma =1.52\times {10}^{11}{\mathrm{s}}^{-1}$. In Layers 1 and 3, plasmons have a momentum larger than the superstrate light line, whereas the plasmon in Layer 2 has a lower momentum. Thus, the plasmon in Layer 2 would have strong radiative losses if the fields penetrate through Layer 1 into the superstrate, whereas the first and third SP modes are still bounded.

Download Full Size | PDF

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 $ϵk_o^2=k_x^2+k_y^2$ where $ϵ$ 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{{ϵ}_{\sup }}{k}_o$). For $k_x$ values greater than $\sqrt{{ϵ}_{\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\mathrm{THz}$ and $k=0.0406{\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\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}{\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.