## Abstract

For the first time, we demonstrate the application of the time domain transmission line method (TLM) to accurate modeling of surface plasmon polariton (SPP) structures. The constructed TLM node allows for modeling of dispersive materials through simple time-difference equations. Using such node, an ultra-wide band excitation can be applied to obtain the response over the band of interest. Bérenger’s perfectly matched layer (PML) boundary condition can readily be implemented using the same node. We illustrate our TLM approach through the modeling of different challenging structures including SPPs filters and focusing structures.

©2010 Optical Society of America

## 1. Introduction

Surface Plasmon Polaritons (SPPs) have been recently proposed for the design of miniaturized photonic devices [1–4]. It has the unique property of guiding, directing, and focusing light energy below the diffraction limit. The ability of the SPPs to guide the light in subwavelength dimensions provides a promising solution for realizing optical circuits with a small physical size and fast operation speed. These requirements are essential in optical interconnects applications. On the other hand, the ability of subwavelength focusing can be exploited in different applications including near and far field scanning nanoscale optical microscopes [5–13]. The nanolens structure based on SPPs has been proposed for integration with CMOS image sensors and other optoelectronic devices [12]. SPP structures have also been utilized in many biosensing applications [14, 15].

Due to its compatibility in dimensions with electronic circuits, SPPs have been proposed for different circuit-based applications [3]. Accurate modeling for the whole circuit is highly desirable. SPPs are usually modeled with electromagnetic solvers both in frequency and time domain. These techniques are based on direct solution of Maxwell’s equations. They include Finite Difference Time Domain (FDTD) [16] and Finite Difference Frequency Domain (FDFD) [12] methods. However, these techniques are still under development for direct efficient modeling of electronic circuits [16].

The Transmission line method (TLM) provides a universal approach for modeling electromagnetic-based structures [17] and integrated circuits [18]. It has also been utilized for the modeling of circuit components [19–21], Micro-Electro-Mechanical Systems (MEMS) [22], and Electro-Magnetic Compatibility (EMC) [23]. The drift and diffusion processes have also been modeled with TLM [24]. Mechanical and heat systems are also modeled by TLM [25]. Time domain TLM is unconditionally stable [26–28] which allows for any excitation ranging from time-domain impulses to Gaussian-modulated sinusoids. However, the application of the TLM in photonics area has been so far limited to few basic structures using simple boundary conditions [29].

In this paper, we demonstrate, for the first time, the modeling of plasmonic photonic structures using a dispersive TLM algorithm. Unlike the conventional TLM techniques for nondispersive media [26–28], we utilize a TLM technique which allows for efficient and accurate modeling of these media [30–33]. We integrate the dispersive TLM approach with the efficient and accurate Perfectly Matched Layer (PML) proposed in [34]. The technique is exploited not only for modeling SPPs but also for studying their design parameter dependence. To demonstrate the universality of the proposed TLM technique and its capabilities, this study is applied to various structures for in-depth understanding of the physical effects of varying parameters.

The paper starts with a brief review of the TLM technique in Section 2. The application of TLM technique to a number of promising surface plasmon polariton structures is dedicated to Section 3. Our work is concluded in Section 4.

## 2. Transmission line modeling

#### 2.1 Formulation of dispersive TLM

We utilize a formulation for the 3-D TLM method based on the symmetrical condensed node (SCN) [28]. This node is shown in Fig. 1
, where (*V*
_{0}
*… V*
_{11}) are the voltage impulses at different ports of the node. The six total field quantities (*E _{x}, E_{y}, E_{z}, H_{x}, H_{y},* and

*H*) at the node centre are related to the incident voltage impulses through linear mappings. We utilize a revised port numbering in order to preserve the symmetry of the scattering matrix of the TLM node [34].

_{z}At each node, we define the following quantities

where *V** ^{i}* and

*V**are the incident and reflected port voltages, respectively. Note that in Eqs. (1) and (2) we drop the time index (*

^{r}*n*) and the node index (

*j*) to simplify the notations. The total circuit voltages and currents (

**) and the node excitations (**

*F*

*F**) are defined by:*

^{ex}

*V**= [*

_{f}*i*]

_{fx}i_{fy}i_{fz}V_{fx}V_{fy}V_{fz}*.*

^{T}At each time step, the TLM carries out the procedures of scattering and connection. The scattering process is carried out at each node to obtain the vector of reflected impulses *V** ^{r}* using the vector of incident impulses and the free sources. First, the cell excitation is calculated using [33]:

*R**is the cell excitation input matrix whose elements belong to the set {0, −1, 1}. Using the calculated cell excitation, the different currents and voltage values are obtained through:where*

_{ex}*t*and

_{ej}*t*,

_{mj}*j*∈{

*x*,

*y*,

*z*} depend on the computational domain normalized constitutive parameters. For example, we haveandwhere

*χ*,

_{e}*χ*,

_{m}*g*, and

_{e}*r*are the electrical susceptibility, magnetic susceptibility, normalized conductance, and normalized resistance, respectively. $\overline{s}$ is the normalized frequency parameter $\overline{s}=s/\Delta t$, where

_{m}*s*= $\sqrt{-1}$

*ω*and Δ

*t*is the time step. The transmission parameters

*t*and

_{ej}*t*, ∀

_{mj}*j*are transformed to the Z-domain utilizing the bilinear transformation:

Finally, the reflected impulses for each cell are given by:

where**is the TLM reflection matrix, and**

*R***is the TLM input / output matrix. The elements of both**

*P***and**

*R***are from the set {0, −1, 1}.**

*P*The second important step in the TLM time stepping is the connection step. In this step, the reflected voltage impulses at each port given by Eq. (11) are incident on the corresponding ports of the neighboring nodes in the next time step. Figure 2
illustrates the connection process for the 2-D case. The incident voltage waves at the *n* + 1 time step (${V}_{n+1}^{i}$) are calculated from the reflected voltage waves at time step *n* (${V}_{n}^{r}$).

#### 2.2 Modeling of plasmonic devices

The dispersive behavior of SPPs can be modeled by Eqs. (8) and (9) through two different approaches. In the first approach, we supply the electric susceptibility in Eq. (8) for the metal in an SPP structure with the dispersion described by the Drude model [35]:

where*χ*is the frequency independent part of the model,

_{∞}*ω*is the bulk plasma frequency,

_{p}*ω*is the angular frequency, and

*γ*is the collision frequency. The expression (12) is substituted in Eq. (8) to result in a time difference equation that obtains the vector

**using the excitation vector**

*F*

*F*^{ex}.

An alternative approach for modeling the metal dispersion in SPP structures is to introduce a dispersive conductivity using the relationship *χ*(*ω*) = $\sqrt{-1}$
*σ*(*ω*)/*ε _{0}ω* [35], where

*χ*(

*ω*) is the frequency-dependent part of the metal susceptibility with

*χ*excluded. Therefore, the metal material is characterized as a non-magnetized plasma with constant susceptibility

_{∞}*χ*and a frequency dependent conductivity

_{∞}*σ*(

*ω*). The frequency dependent conductivity is given by:

*τ*is the collision time constant. The DC electric conductivity is

_{c}*σ*= ${\omega}_{p}^{2}{\epsilon}_{0}{\tau}_{c}$with

_{0}*ε*denoting the air permittivity. Our TLM realization is based on dispersive equivalent conductivity approach.

_{0}For modeling the frequency dependent conductivity, we utilize a technique based on the *Z* transform [33]. This discrete *Z* transform is compatible with the numerical simulation where the fields are calculated at discrete time steps. The normalized conductance is modified to:

*g*

_{e}_{0}=

*g*(1-

_{ec}*β*

_{0}),

*g*

_{e}_{1}= α /(1-

*z*

^{-}^{1}

*β*

_{0}) and$\alpha ={g}_{ec}(1-{\beta}_{0}^{2})$. Equation (15) is then substituted into Eq. (8) to result in a time domain difference equation. For example, for the calculation of

*V*, the transmission parameter

_{x}*t*in Eq. (8) is calculated utilizing the plasmonic material model in Eq. (15) as follows:

_{ex}*χ*= 0,

_{∞}*g*= 0, and

_{e0}*g*= 0). No memory is required in this case to store past field values. For a material of finite frequency independent permittivity and conductivity, the previous field value needs to be stored. For the plasmonic materials, the field history for the previous two time steps is accumulated. This is due to the

_{e1}*z*

^{−1}and

*z*

^{−2}terms in the transmission coefficients. The application of TLM to other dispersive materials is straight forward and can be achieved with a number of extra storages dependent on the order of the material models.

#### 2.3 Perfectly matched layer

The perfectly matched layer (PML) developed by Bérenger is widely utilized with FDTD codes [36]. For the implementation of a PML at the boundary in TLM, we utilize an adapted version of this approach. This approach is based on the manipulation of Maxwell’s equations to allow for nonphysical contributions that lead to a finite computational domain.

The implementation of PML to the TLM computational domain can be achieved through different approaches. The first approach is to implement the PML through an attached FDTD layer [37]. This method requires an interpolation at the domain interface which may be a source of artificial spurious modes. It also adds the complexity of using two different numerical techniques for the same problem.

A second approach utilizes a modified TLM node suitable for the implementation of the PML for TLM [38]. The overhead of this approach is an increased memory storage, computational effort and complexity. An extra 8 links are attached to each node for implementing PML [38]. This approach, however, leads to a unified algorithm which is more robust and accurate and does not require any interpolation at the interface between different computational domains.

In this work, we developed an optimized version of the PML which utilizes the same node for the computational domain and the PML. This PML realization does not require any extra links [34]. It is simple to implement with only a memory overhead due to the presence of extra field values. The calculation of the reflected waves (scattering) and new incident waves (connection) proceed as explained earlier.

According to the basic idea introduced by Bérenger, the TLM implementation [34] uses a modified set of voltage quantities ** V** = [

*V*]

_{xz}V_{xy}V_{yx}V_{yz}V_{zy}V_{zx}i_{xz}i_{xy}i_{yx}i_{yz}i_{zy}i_{zx}^{T}. The vector of the cell excitation fields

*F**is modified to ${F}^{ex}={[{V}_{xz}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{V}_{xy}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{V}_{yx}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{V}_{yz}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{V}_{zy}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{V}_{zx}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{xz}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{xy}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{yx}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{yz}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{zy}^{ex}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{i}_{zx}^{ex}]}^{T}$. The calculations of both*

^{ex}

*F**and*

^{ex}**follow the same procedure as in Eqs. (5) and (6) with modified**

*F*

*R**and*

_{ex}**matrices. The real**

*T*six field components can be calculated by simply adding the corresponding two field quantities as follows:

andThe scattering procedure utilizes Eq. (11) while the connection step is similar to that explained in Fig. 2.## 3. Numerical results and discussions

In this Section, we demonstrate the application of TLM to the modeling of plasmonic structures. We perform a wideband time-domain simulation of a number of structures using our in-house MATLAB-based TLM simulator. Our examples include surface plasmon polariton Bragg gratings, subwavelength beam focusing using surface dielectric gratings, improved focusing using surface metallic grooves, and focusing using nanoscale slit array. Our results are compared to the results available in the literature. We also utilize the TLM simulation to predict the dependence of the response on some device parameters. Our examples are tested using both Gaussian-modulated sinusoidal excitations and impulse excitations in order to illustrate the stability of the TLM approach.

To test our in-house simulator, we start by simulating a simple slab waveguide (see inset of Fig. 3
). This waveguide is made of a cladding and a core layer with refractive index of 1.5 and 1.7, respectively. The width of the waveguide is *d* = 0.8 μm. This structure is excited with the fundamental mode. The excitation signal is a Gaussian-modulated sinusoid with a wavelength between 1.0 μm to 3.0 μm. The transmitted signal inside the waveguide is calculated using both the line matching boundary condition [29] and the PML boundary condition. The results of this simple waveguide simulation are shown in Figs. 3 and 4
. Figure 3 shows a large time domain reflection for the line matching boundary condition. The transmission spectrum using the same boundary is highly distorted as shown in Fig. 4. The PML boundary provides far more accurate results especially for wide band excitations as shown in Figs. 3 and 4. Our PML realization is utilized in all subsequent examples. These examples are all conducted on an Intel® Xeon® Processor 5160 (3.0 GHz) platform.

#### 3.1 Surface plasmon Bragg grating

In this example, we model a metal-insulator-metal (MIM)-type waveguide with metallic grating (see Figure 5
). Using this grating, we can achieve Bragg filtering. This device is known as Surface Plasmon Polaritons Waveguide Bragg Grating (SPP-WBG) [39]. For wideband modeling of these devices, time domain simulators are preferred. The parameters are chosen to be *W* = 500.0 nm, *W*
_{1} = 100.0 nm, Λ = 660.0 nm, and *L* = 292.0 nm [40]. The values of the last two parameters are chosen to allow the Bragg condition to occur at a wavelength of 1550.0 nm. The number of grating pitches is set to 14. The thickness of the grating *t* is included as a parameter in order to evaluate its impact on the response. We set the minimum value of *t* to 25.0 nm. The metal is assumed to be silver with its permittivity given by the Drude model. The parameters of the model are *ε _{∞}* = 3.7,

*ħω*= 9.1eV, and

_{p}*ħ*γ = 0.018eV [41–43], where

*ħ*is the reduced Planck’s constant. The dielectric layer is assumed to be air.

The SPP Bragg grating is simulated using a TLM space step of 10.0 nm for 30,000 time steps. The simulated structure has dimensions 0.4 × 13 μm which requires 40 × 1300 TLM cells. For this large structure with a large number of time steps, the overall simulation time is about 90 minutes.

A time-domain impulse excitation is utilized in this example. The transmission spectrum is estimated for the wavelength range 1.0 μm to 6.0 μm as shown in Fig. 6
. The predicted center wavelength inside the band gap is 1550.0 nm. Both the center wavelength and the quality factor have good agreement to those reported in [40, 44]. The magnetic field distribution (|*H _{y}*|) inside the SPP-WBG is shown in Figs. 7
and 8
for incident wavelengths 1000.0 nm and 1550.0 nm, respectively. It is clear that the incident wave with a wavelength of 1000.0 nm, which lies outside the band gap, is transmitted. However, at 1550.0 nm, the incident wave is completely reflected. The asymmetry in the wavelength response in Fig. 6 is consistent with the dispersive permittivity dependence of silver.

For robust validation of our model, we estimate the impulse response due to different values of the grating thickness *t*. This is the first reported simulated impulse response of a photonic structure. This is the direct advantage of the TLM which is a physics based numerical technique. The transmission response is calculated for thicknesses ranging from 25.0 nm to 155.0 nm as shown in Fig. 9
. The grating thickness affects both the band gap and the performance of the device. As the thickness increases, the band gap increases with a significant shift in the mid wavelength and the transmitted energy decreases. The thickness increase also allows for a large propagation loss inside the grating which leads to a degradation in the filter response as shown in Fig. 9.

#### 3.2 Sub-wavelength beam focusing using dielectric surface grating

In this example, we model a single subwavelength slit surrounded by chirped dielectric surface grating as shown in Fig. 10
. By introducing a chirped dielectric grating on the metallic surface, the light is focused at a certain focal length that is dependent on the dielectric material and the chirped grating period [45, 46]. This device can be used for focusing the light in nanoscale scanning optical microscopy. It consists of silver slab of thickness *d* with a slit of width *a.*

For the considered structure, we set *d* = 1.0 μm, *a* = 100.0 nm, and *t* = 120.0 nm. The structure consists of 12 surface grating pitches at each side of the slit with a 0.5 fill factor. The periods of these grating are given by Λ = [327.0 294.0 272.0 258.0 249.0 242.0 238.0 234.0 231.0 229.0 228.0 226.0] nm which allows the focal point of an incident wave at 532.0 nm to be at 1.5 μm. The same silver dispersion model of Section 3.1 is utilized with dispersion parameters consistent with [47]. The grating is made of the dielectric material with a refractive index of 1.72.

The structure is simulated with a space step of 5.0 nm for 8000 time steps utilizing different types of excitation. The space step is fairly small to match the contrast in the design dimensional parameters. The simulation is conducted for a fairly large domain size to ensure that the focal point is located inside the computational domain. For a 7 × 10 μm structure, the computational domain includes 1400 × 2000 TLM cells. The simulation time is approximately 300 minutes.

The field distribution (|*H _{y}*|) at an operating wavelength of 532.0 nm is shown in Fig. 11
. The achieved focal point is 1.5 μm which is identical to that reported in [46] with the same full width at half maximum (FWHM = 409.0 nm) as shown in Fig. 12
.

In Fig. 13
, the longitudinal intensity distribution (|*Hy*|)^{2} for a selected set of wavelengths demonstrates the variation of the focal length. As can be seen, focusing can only be achieved at wavelengths 500.0 nm, 550.0 nm, and 600 nm according to the designed structure. For wavelengths ranging from 500.0 nm to 600.0 nm, the structure produces the spectral focal length dependence shown in Fig. 14
. As the light wavelength increases the focal point becomes closer to the slit.

#### 3.3 Sub-wavelength focusing using metallic grooves

We consider the promising focusing structure shown in Fig. 15
. It consists of a metallic slit surrounded with metallic grooves of constant depths [5–8]. These grooves allow for highly focused beam without any chirping. The main parameters are the slit width and groove period, depth, and width. The structure utilizes a gold film with a thickness of 400.0 nm. The film dielectric constant is characterized by the Drude model with parameters given as *ε _{∞}* = 9.0685,

*ħω*= 8.9269eV, and

_{p}*ħ*γ = 0.076eV [41, 42]. The dielectric layer is made of silicon dioxide (SiO

_{2}) with a refractive index of 1.46. The slit width is 200.0 nm. The groove period, width, and depth are 380.0 nm, 60.0 nm and 40.0 nm, respectively. This structure provides focusing at wavelength of 650.0 nm with a focal point of 2.0 μm [6].

This structure is modeled using a computational domain with a space step of 10.0 nm for 8000 time steps. The computational effort for this simulation is less than the design in Section 3.2 because the grating is not chirped. For a 5 × 12 μm structure, the computational domain includes 500 × 1200 TLM cells. The simulation time is approximately 90 minutes.

The structure is excited with a wideband Gaussian-modulated sinusoidal signal (from 400.0 nm to 700.0 nm). For this example, we demonstrate the focusing by showing the change in the field distribution by introducing the grooves. The field distribution (|*H _{y}*|) for the structure without metallic grooves is shown in Fig. 16
. The field distribution of the focusing structure is shown in Fig. 17
. Here, light is focused at a focal length of 2.0 μm. In Fig. 18
, the intensity distribution along the transverse direction at the focal point is compared to that obtained in [6]. The TLM simulation predicts a FWHM of 350.0 nm which is in good agreement with the reported data [6].

The TLM time domain simulation is further exploited to predict the focusing behavior for different wavelengths. The lateral intensity distribution (|*H _{y}*|)

^{2}for a selected set of wavelengths is shown in Fig. 19 . Focusing is achieved in the neighborhood of the wavelength 650.0 nm. For the wavelength range between 600.0 nm to 700.0 nm, the continuous focal-length wavelength relation is shown in Fig. 20 .

#### 3.4 Slit Array in a Metallic Film

In this example, we model the nanoscale slit array shown in Fig. 21 using time-domain TLM. This structure has been utilized for angle compensation which enhances the efficiency of the off-axis pixels in solid state image sensors [11]. It also allows for beaming and focusing effects which can be utilized in optical microscopy [9–15].

The considered structure has a metal thickness *h* = 400.0 nm. The slit widths are ranging from 80.0 nm in the middle to 150.0 nm at edges [12]. The dielectric constant of SiO_{2} is 2.13. The metallic film is made of gold which has the same parameters as in the previous example. These parameters lead to a dielectric constant of −11.04 + 0.78i for the gold at a wavelength of 637.0 nm [12].

We set the cell dimension of the TLM domain to 10.0 nm for 8000 time steps. For 5 × 16 μm structure, the computational domain includes 500 × 1600 TLM cells. The simulation time is approximately 100 minutes. The excitation exploits a wide band Gaussian-modulated sinusoid with wavelength ranging from 400.0 nm to 800.0 nm. In Fig. 22 , the magnetic field distribution is shown for an incidence wavelength of 637.0 nm. The achieved focal length is 5.3 μm which is identical to that reported experimentally [12]. The field distribution around the focal spot in the transverse direction is compared to the reported experimental result in Fig. 23 . Both results have a good agreement.

From the conducted time domain simulation, the focusing effect for different incidence wavelengths is explored. Figure 24 shows the longitudinal field distribution at the structure axis of the symmetry. The maximum field indicates the position of the focal point. It is clear that focusing is achieved for all wavelengths at different focal lengths. The focusing effect for this structure is based on the interference pattern constructed from the different waves coming out of the slits. The optical size of the lens determines the focal point. For an increased wavelength, the optical length is decreased leading to a scaled-down focal length. As the wavelength increases, the focal point becomes closer to the lens as shown in Fig. 25 .

## 4. Conclusion

We successfully demonstrated, for the first time, the application of the TLM method to the time-domain modeling of Surface Plasmon Polaritons (SPPs) devices. A dispersive TLM node is utilized which allows for modeling any dispersion behaviour in the same way through simple difference equations. Wideband responses of various focusing structures have been accurately modeled using impulse excitations or Gaussian modulated sinusoids. Our results match well the reported results in the literature.

## References and links

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

**2. **S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, and T. W. Ebbesen, “Channel plasmon-polariton guiding by subwavelength metal grooves,” Phys. Rev. Lett. **95**(4), 046802 (2005). [CrossRef] [PubMed]

**3. **E. Ozbay, “Plasmonics: merging photonics and electronics at nanoscale dimensions,” Science **311**(5758), 189–193 (2006). [CrossRef] [PubMed]

**4. **C. Genet and T. W. Ebbesen, “Light in tiny holes,” Nature **445**(7123), 39–46 (2007). [CrossRef] [PubMed]

**5. **F. J. García-Vidal, L. Martín-Moreno, H. J. Lezec, and T. W. Ebbesen, “Focusing light with a single subwavelength aperture flanked by surface corrugations,” Appl. Phys. Lett. **83**(22), 4500–4502 (2003). [CrossRef]

**6. **B. Guo, Q. Gan, G. Song, J. Gao, and L. Chen, “Numerical study of a high-resolution far-field scanning optical microscope via a surface plasmon modulated light source,” J. Lightwave Technol. **25**(3), 830–833 (2007). [CrossRef]

**7. **H. Shi, C. Du, and X. Luo, “Focal length modulation based on a metallic slit surrounded with grooves in curved depths,” Appl. Phys. Lett. **91**(9), 093111 (2007). [CrossRef]

**8. **B. Lee, S. Kim, H. Kim, and Y. Lim, “The use of plasmonics in light beaming and focusing,” Prog. Quantum Electron. **34**(2), 47–87 (2010). [CrossRef]

**9. **Z. Sun and H. K. Kim, “Refractive transmission of light and beam shaping with metallic nano-optic lenses,” Appl. Phys. Lett. **85**(4), 642–644 (2004). [CrossRef]

**10. **H. Shi, C. Wang, C. Du, X. Luo, X. Dong, and H. Gao, “Beam manipulating by metallic nano-slits with variant widths,” Opt. Express **13**(18), 6815–6820 (2005). [CrossRef] [PubMed]

**11. **L. Verslegers, P. B. Catrysse, Z. Yu, and S. Fan, “Planar metallic nanoscale slit lenses for angle compensation,” Appl. Phys. Lett. **95**(7), 071112 (2009). [CrossRef]

**12. **L. Verslegers, P. B. Catrysse, Z. Yu, J. S. White, E. S. Barnard, M. L. Brongersma, and S. Fan, “Planar lenses based on nanoscale slit arrays in a metallic film,” Nano Lett. **9**(1), 235–238 (2009). [CrossRef]

**13. **D. Choi, Y. Lim, S. Roh, I. M. Lee, J. Jung, and B. Lee, “Optical beam focusing with a metal slit array arranged along a semicircular surface and its optimization with a genetic algorithm,” Appl. Opt. **49**(7), A30–A35 (2010). [CrossRef] [PubMed]

**14. **J. Homola, S. S. Yee, and G. Gauglitz, “Surface plasmon resonance sensors: review,” Sens. Actuators B Chem. **54**(1-2), 3–15 (1999). [CrossRef]

**15. **R. D. Harris and J. S. Wilkinson, “Waveguide surface plasmon resonance sensors,” Sens. Actuators B Chem. **29**(1-3), 261–267 (1995). [CrossRef]

**16. **A. Taflove, and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method, Third Edition*, 3rd ed. (Artech House Publishers, 2005).

**17. **Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, “Transmission Line and Equivalent Circuit Models for Plasmonic Waveguide Components,” IEEE J. Sel. Topics in Quantum Mechanics **14**(6), 1462–1472 (2008). [CrossRef]

**18. **J. R. Brews, “Transmission line models for lossy waveguide interconnect in VLSI,” IEEE Trans. Electron. Dev. **33**(9), 1356–1365 (1986). [CrossRef]

**19. **P. B. Johns and M. O’Brien, “Use of transmission-line modelling method to solve non-linear lumped networks,” Radio Electron. Eng. **50**, 59–70 (1980). [CrossRef]

**20. **S. Hui, K. FUNG, and C. Christopoulos, “Decoupled simulation of multistage power electronics systems using transmission-line links”. Power electronics specialist conference, (Toledo, Spain, 1992), pp. 1324–1330.

**21. **K. Fung and S. Hui, “Improved TLM link model for reactive circuit components,” IEE Proc. Sci. Meas. Technol. **143**(5), 341–344 (1996). [CrossRef]

**22. **L. Vietzorreck, F. Coccetti, V. Chtchekatourov, and P. Russer, “Modeling of MEMS capacitive switches by TLM,” *in 2000 IEEE MTT-S International Microwave Symposium Digest**,* (2000), pp. 1221–1224.

**23. **C. Christopoulos and J. Herring, “The application of transmission-line modeling (TLM) to electromagnetic compatibility problems,” IEEE Trans. Electromagnetic Compatibility **35**(Part 2), 185–191 (1993). [CrossRef]

**24. **D. Pasalic and R. Vahldieck, “A Hybrid Drift-Diffusion-TLM Analysis of Traveling-Wave Photodetectors,” IEEE Trans. Microw. Theory Tech. **53**(9), 2700–2706 (2005). [CrossRef]

**25. **D. D. Cogan, *Transmission line matrix (TLM) techniques for diffusion applications*, (CRC Press, 1998).

**26. **P. B. Johns and R. L. Beurle, “Numerical solution of 2-dimensional scattering problems using a transmission-line matrix,” Proc. Inst. Electr. Eng. **118**(9), 1203–1208 (1971). [CrossRef]

**27. **C. Christopoulos, *The Transmission-Line Modeling Method: TLM*, (Wiley-IEEE Press, 1995).

**28. **P. B. Johns, “A symmetrical condensed node for the TLM method,” IEEE Trans. Microw. Theory Tech. **35**(4), 370–377 (1987). [CrossRef]

**29. **O. Jacquin, F. Ndagijimana, A. Cachard, and P. Benech, “Application of the TLM technique to integrated optic component modelling,” Int. J. of Numer. Model: Electronic Networks Devices and Fields **14**, 95–105 (2001). [CrossRef]

**30. **L. de Menezes, and W. J. R. Hoefer, “Modeling frequency dependent dielectrics in TLM,” in IEEE Antennas Propagat.-S Symp. Dig., (1994), pp. 1140–1143.

**31. **L. de Menezes, and W. J. R. Hoefer, “Modeling nonlinear dispersive media in 2-D-TLM,” in 24th European Microwave Conference Proc., (1994), pp. 1739–1744.

**32. **L. de Menezes and W. J. R. Hoefer, “Modeling of general constitutive relationships in SCN TLM,” IEEE Trans. Microw. Theory Tech. **44**(6), 854–861 (1996). [CrossRef]

**33. **J. Paul, C. Christopoulos, and D. W. P. Thomas, “Generalized material models in TLM—Part I: Materials with frequency-dependent properties,” IEEE Trans. Antenn. Propag. **47**(10), 1528–1534 (1999). [CrossRef]

**34. **J. Paul, C. Christopoulos, and D. W. P. Thomas, “Perfectly matched layer for transmission line modelling (TLM) method,” Electron. Lett. **33**(9), 729–730 (1997). [CrossRef]

**35. **S. A. Maier, *Plasmonics: Fundamentals and Applications*, (Springer, 2007).

**36. **J. P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**(2), 185–200 (1994). [CrossRef]

**37. **C. Eswarappa and W. Hoefer, “Implementation of Berenger absorbing boundary conditions in TLM by interfacing FDTD perfectly matched layers,” Electron. Lett. **31**(15), 1264 (1995). [CrossRef]

**38. **S. Le Maguer, N. Peňa, and M. Ney, “Matched absorbing medium techniques for full-wave TLM simulation of microwave and millimeter-wave components,” Ann. Telecommun. **53**, 115–129 (1998).

**39. **B. Wang and G. P. Wang, “Plasmon bragg reflectors and nanocavities on flat metallic surfaces,” Appl. Phys. Lett. **87**(1), 013107 (2005). [CrossRef]

**40. **Z. Han, E. Forsberg, and S. He, “Surface plasmon Bragg gratings formed in metal-insulator-metal waveguides,” IEEE Photon. Technol. Lett. **19**(2), 91–93 (2007). [CrossRef]

**41. **C. Sönnichsen, *Plasmons in metal nanostructures*, Ph.D. dissertation, Ludwig-Maximilians-Universtät München, (München, 2001).

**42. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**(12), 4370–4379 (1972). [CrossRef]

**43. **J. Park, H. Kim, and B. Lee, “High order plasmonic Bragg reflection in the metal-insulator-metal waveguide Bragg grating,” Opt. Express **16**(1), 413–425 (2008). [CrossRef] [PubMed]

**44. **A. Hosseini, H. Nejati, and Y. Massoud, “Modeling and design methodology for metal-insulator-metal plasmonic Bragg reflectors,” Opt. Express **16**(3), 1475–1480 (2008). [CrossRef] [PubMed]

**45. **D. Z. Lin, C. K. Chang, Y. C. Chen, D. L. Yang, M. W. Lin, J. T. Yeh, J. M. Liu, C. H. Kuan, C. S. Yeh, and C. K. Lee, “Beaming light from a subwavelength metal slit surrounded by dielectric surface gratings,” Opt. Express **14**(8), 3503–3511 (2006). [CrossRef] [PubMed]

**46. **S. Kim, Y. Lim, H. Kim, J. Park, and B. Lee, “Optical beam focusing by a single subwavelength metal slit surrounded by chirped dielectric surface gratings,” Appl. Phys. Lett. **92**(1), 013103 (2008). [CrossRef]

**47. **E. D. Palik, *Handbook of optical constants of solids*. (Academic Press, 1998).