## Abstract

In this paper, we present a full-vector finite difference method to solve for optical modes in one and two dimensional subwavelength plasmonic waveguides. We have used the Implicitly Restarted Arnoldi method to directly calculate the propagation constants of the dominant modes. The method has low computational complexity and can be applied to accurately model complex geometries and structures with fast-varying field profiles. When applied to solve for purely bounded modes, our method automatically separates evanescent and low-loss guided modes.

© 2006 Optical Society of America

## 1. Introduction

The interface between two materials with dielectric constants of opposite signs can support propagation of surface electromagnetic modes known as surface plasmon-polaritons (SPPs). Unlike the conventional dielectric waveguides that confine electromagnetic waves to an optically dense core, SPPs are localized at the interfaces between dielectric materials and metals or ionic solids that support charge density oscillation. Energy and information can be transmitted in through plasmonic waveguides beyond the diffraction limit in highly integrated devices [1, 2]. The most commonly used method to excite SPPs has been prism coupling through attenuated total reflection which in fact matches the SPP wave vector parallel to the metal-dielectric surface to that of the incident radiation.

Prism coupling to the desired modal wave vector parallel to the metal-dielectric surface is the most commonly used method to excite SPPs. This method has been widely utilized in SPP-based chemical senors. However, it is inefficient and topologically inappropriate for waveguid-ing application because the energy of the incident radiation weakly couples to the waveguide [3]. On the other hand, end-fire excitation has been experimentally verified to be an efficient way to excite the symmetric mode of metallic slabs [4]. The efficiency of an end-fire excitation directly depends on matching the field profile of the desired propagating mode and the field profile of the raditation that excites it. Therefore, the numerical modeling of the field profile is crucial for the realization of SPP-based applications.

Two important considerations for utilizing plasmonic waveguides as optical nanowires in future optical integrated circuits are the light confinement and the propagation loss. A significant limitation of SPP waveguides is the propagation loss, which originates from resistive heating losses in the metal. As a result, SPPs cannot travel more than a few micrometers before experiencing a high degree of attenuation. However, these distances are sufficient for many applications nanoscale optical integrated circuits. Although the bounded modes of plasmonic waveguides are mainly concentrated at the interfaces between the metal and dielectric materials, the electromagnetic field spreads into both media. Therefore, the light is less confined to the waveguide structure, which may produce undesired coupling between adjacent waveguides.

Investigating the transmission properties of plasmonic waveguides requires accurate and efficient modeling techniques. The method of lines (MoL) and a full-vector finite difference analysis have been the two primary numerical approaches used to model 2D plasmonic waveguides [6, 7]. Both methods have been applied to characterize single wide and thin metal strips. The computational complexity of these methods is exasperated by the eigenvalue problem associated with modal solutions of waveguides. MOL also requires a complex zero-searching algorithm to find the zeros of the determinant of a large matrix [8]. In addition to the bound modes, the leaky modes of surface plasmon waveguides can be also determined using vectorial finite difference analysis with transparent or absorbing boundary conditions [9].

In this paper, we present a fast and accurate full-vector finite difference analysis method to simulate 2D plasmonic waveguides by efficiently solving the eigenvalue problem associated with the formulation. The method does not generate spurious modes. Furthermore, the evanescent modes are well-separated from the propagating modes. We show that this method can be easily applied to 1D structures. We validated the method against one of the most established and accurate numerical approaches for modeling 1D planar waveguides [10]. The method calculates the field profile directly and can be applied to model any 1D and 2D structure. The simulation results for metallic strips embedded in dielectric medium are presented. These results will be a step toward generalizing one-dimensional plasmonic structures [5, 11, 12] to practical two-dimensional cases.

## 2. Numerical solution approach

#### 2.1. Theory

In [13], an accurate and numerically stable finite difference method was proposed for the vectorial analysis of optical waveguides. The method is based on transverse magnetic field components to prevent spurious modes and was applied to solve for the optical modes of subwave-length 2D plasmonic waveguides [6]. The absence of spurious modes is beneficial in the simulation of realistic plasmonic structures where the large imaginary part of the dielectric constant of metals at optical frequencies results in numerous unbounded and dissipative modes.

In this method, the homogeneous Helmholtz equations for transverse H-field, Eqs. (1) and (2), and the continuity of the longitudinal field components, Eqs. (3) and (4), are

where *β* is the propagation constant assuming a functional form of *exp*(*i*(*βz* - *ωt*)) for the steady state modes, *k*
_{0} is the wave number in the free space, *ε* is the dielectric constant, and *n*_{eff}
= *β*/*k*
_{0} is the effective index. If the derivatives in these equations are approximated with their equivalent finite differences, all the discrete equations can be arranged into a linear system

where *H̅* is a vector containing field values at all the points on the mesh [13]. All material interfaces fall on the grid points. The computational window must be large enough to ensure zero field values at its boundaries. We have chosen even meshing to discretize the main structure and uneven meshing for the rest of the computational window. If the eigenvalue problem in Eq. (5) is solved using the algorithm proposed in [13], its solution requires an intractable period of time. This is also a limiting factor for the simulation of complex structures, where a large number of points are present in the computational mesh. Furthermore, when solving Eq. (5) using direct eigenvalue decomposition techniques, the number of the calculated eigenvalues equals twice the number of the total points in the computational window. However, most of the eigenvalues have large imaginary parts, and therefore, correspond to evanescent modes. Therefore, the small number of propagating modes coupled with the sparse nature of *L* motivate the use of the Arnoldi method to iteratively calculate the eigenvalues. Alternative iterative eigenvalue algorithms such as [14] that are based on the subspace iteration method are significantly less efficient since they only calculate one eigenvalue at a time and require the direct computation of the system’s eigenvectors [15]. In contrast, using techniques based on the Arnoldi iteration, we can simultaneously solve for all the propagating modes.

#### 2.2. Eigenvalue computation

Even for modest problems sizes with uneven 2-D meshing, the order of *L* in Eq. (5) can exceed 10000. Direct methods for solving the eigenvalue decomposition problem typically require *O*(*n*
^{2}) storage and *O*(*n*
^{3}) operations. Therefore, solving eigenvalue problems for matrices with order *n* larger than several thousand elements is intractable. Since we are interested in the dominant resonant modes of the simulated structure, we only need to calculate a few eigenvalues (*β*
^{2}) that meet the following criteria: (a) *Re*(*β*
^{2}) > 0, (b) *Im*(*β*
^{2}) > 0, and (c) *Re*(*β*
^{2}) >> *Im*(*β*
^{2}).

Iterative eigenvalue decomposition methods based on the Arnoldi iteration have the capability to calculate a small number of eigenvalues near a desired value. To motivate the use of the Arnoldi iteration, we first discuss the power method, which is a traditional technique for iterative eigenvalue computation. In its simplest implementation, the power method successively multiplies a random vector *b* by powers of a given matrix *L* in order to find the largest eigenvalue of *L*. The power method explicitly forms the Krylov Subspace (*K*_{n}
), which is defined as

The power method has several drawbacks. First, it cannot locate multiple eigenvalues without employing block decomposition methods, which have significant computational overhead [16]. Furthermore, as higher powers of *L* are explicitly generated, the resulting basis vectors that span *K*_{n}
become numerically ill-conditioned, which can ultimately lead to inaccurate results. In contrast to the power method, the Arnoldi iteration projects- the matrix *L* onto the Krylov Subspace (*K*_{n}
) iteratively forming a set of orthogonal basis vectors (*V*_{m}
)

where ${V}_{m}^{*}$
*V*_{m}
= *I* and *f*_{m}${e}_{m}^{T}$
is the residual error at each step of the iteration. The set of basis vectors generated by the Arnoldi iteration forms an upper Hessenberg matrix (*H*_{m}
) of order *m* where *m* << *n*,

The eigenvalues of *H*_{m}
are estimates for several of the eigenvalues of *L*. Note that the Arnoldi iteration is closely related to the Lanczos iteration since Lanczos is essentially Arnoldi applied to symmetric matrices. In this case, *H*_{m}
is reduced to a tri-diagonal matrix, which provides additional computational savings over the Arnoldi iteration. However, for the 2D finite difference solver developed in this paper, *L* is not symmetric. The standard Arnoldi iteration locates the eigenvalues with the largest magnitude. However, by applying shifts of *μI* to *H*_{m}
during the Arnoldi iteration, we transform the eigenspace of *L* so that the Arnoldi iteration locates eigenvalues near *μ*.

To efficiently calculate the eigenvalues of interest and their associated eigenvectors, we utilize the Implicitly Restarted Arnoldi method, which is based on the Arnoldi iteration described above [18, 17]. Since the size of *V*_{m}
and *H*_{m}
grows as the number of iterations increases, using Eq. (7) directly to compute the eigenvalues can become expensive. Therefore, restarting the Arnoldi iteration based on the information obtained in the previous iterations can lead to a substantial performance improvement since the size of *H*_{m}
does not exceed *k* where *k* = *m*+*p*, *m* is the number of desired eigenvalues, and *p* is the number of extra iterations performed. The algorithm is implicitly restarted by using the *p* extra eigenvalue estimates as shifts in a manner similar to the implicitly shifted QR algorithm. The shifted *QR* factorization effectively reduces the subdiagonal elements of *H*_{m}
to 0 to form the matrix *Ĥ*. Therefore, an approximate partial Schur decomposition is obtained, Q *Ĥ**Q*
^{*}, where *Q* is an orthogonal vector and *Ĥ* is an upper triangular matrix. Based on the properties of Schur decomposition, the diagonal entries of *Ĥ* are the eigenvalue estimates of *L* at the end of each iteration of the algorithm. More information on the Implicitly Restarted Arnoldi algorithm and critical implementation considerations can be found in [17, 19].

Figure 1 displays the eigenvalues of a 1680×1680 *L* matrix and the 13 eigenvalues of largest real part calculated using Implicitly Restarted Arnoldi. It is apparent that only a few of the calculated eigenvalues correspond to bounded propagation modes. Therefore, we signifcantly reduce the computational complexity by calculating only the eigenvalues that satisfy the criteria. For the 1680×1680 matrix, the solution to the eigenvalue problem using direct methods requires approximately 86.8 seconds of CPU time while the iterative eigenvalue computation only needs 1.1 seconds. As the order of *L* is increased, the computational savings provided by the iterative method becomes even more significant. For a 2100×2100 *L* matrix, the direct eigenvalue solution requires 677.8 seconds while the iterative method only needs 2.25 seconds to calculate 20 eigenvalues. Given the efficiency of iterative eigenvalue computation using Implicitly Restarted Arnoldi, we are able to calculate the desired eigenvalues and their associated eigenvectors for matrices with *n* greater than 100,000, which are required for complex structures.

## 3. Results and discussion

#### 3.1. Planar multilayered waveguide

We have validated our numerical technique and implementation against a well-established and accurate numerical method for computing the electromagnetic modes supported by multilayered planar optical waveguides constructed from lossy media [10]. The results are obtained for the wavelength *λ* = 1.55 *μm*. The metal is silver with a dielectric constant of *ε*_{m}
= -125.72 + *i*3.23. A dielectric constant of *ε*_{d}
= 12.25 is assumed for the dielectric material. In order to solve for the bounded electromagnetic modes, zero boundary conditions are applied to the points lying on the boundary of the computational window. Also, as described in the left inset of Fig. 2, if periodic boundary conditions are imposed, a 1D waveguide with infinite width can be modeled. The dispersion curves for the slab fundamental (symmetric bound *S*_{b}
) mode of a 1D metal-insulator-metal (MIM) structure with a dielectric region thickness of 100 *nm* is calculated using our method and the method proposed in [10], which are compared in Fig. 2. Our results are in excellent agreement the standard 1d method from [10]. The maximum error of the calculated *n*_{eff}
values using our modeling technique with a moderate resolution of 20 is less than 0.3% with respect to the values obtained using [10]. We define the resolution as the number of points in the computational mesh on the shortest dimension of th waveguide.

Furthermore, when applied to 1D structures, the method still provides an efficient modeling solution because the computational complexity of solving the eigenvalue problem has been greatly reduced. In addition, for 1D structures, typically fewer eigenvalues are needed. Therefore, in these cases, the simulation is significantly faster than 2D ones. Table 1 summarizes the size of *L*, the required CPU time, the *n*_{eff}
values obtained from our simulation method for an MIM waveguide of 100 *nm* thickness with different resolutions. With a resolution of 10, the error is less than 0.5%, and the waveguide is simulated in 26 seconds.

#### 3.2. Subwavelength metallic strip in a dielectric medium

We model a metallic strip embedded in dielectric medium and compare the results with the method in [6]. The dispersion relation curves of this waveguide structure are shown in Fig. 3. For *λ* = 1.55 *μm*, *ε*_{m}
= -125.72+*i*3.23, *ε*_{d}
= 12.25, we calculate the *n*_{eff}
values of the first three guided modes, *M*
_{00} (Fig. 4), *M*
_{01} (Fig. 5), and *M*
_{10} (Fig. 6), for a single metallic strip of fixed width *w* = 1 *μm* and varying thickness. Two almost degenerate corner modes are represented by curve *M*_{C}
. This degeneracy originates from the fact that the absolute value of the field profiles of these two modes are almost the same. These modes are concentrated at the four corners and are highly lossy [6]. The symmetric (*S*_{b}
) and antisymmetric (*A*_{b}
) slab modes shown in Fig. 3 correspond to the guided modes of an insulator-metal-insulator (IMI) planar waveguide.

In Fig. 3, the *M*
_{00} and *M*
_{10} dispersion diagrams are the same as those obtained in [6]. In contrast, the y-antisymmetric *M*
_{01} mode and the corner modes curves are slightly different. However, *M*
_{01} shown in this figure is well predicted by the dielectric model for guided surface polaritons that use the effective index of a planar IMI waveguide to predict the effective indices of the guided modes for a 2D metallic strip waveguide of the same thickness [20]. This discrepancy may be due to the discretization error introduced by modeling the fast varying fields of the y-antisymmetric modes with a finite-difference scheme [20]. Using the iterative eigenvalue solver introduced in Section 2.2, we can reduce the grid size in order to decrease the discretization error while simulating the structure in a reasonable time. Also, more complex structures can be characterized. For a waveguide consisting of two parallel metallic strips of *t* = 100 *nm* and *w* = 1 *μm* with a separation of 100 *nm*, we found that the fundamental mode has an effective index of *n*_{eff}
= 4.24+*i*0.012. The field profile associated with this mode is shown in Fig. 7. Interestingly, the field profile and effective index of this mode are similar to those of the fundamental mode of a single dielectric strip in a metallic medium [5]. This indicates that multiple-strip metallic waveguides can be used to effectively confine the light to subwavelength area.

## 4. Conclusion

In this paper, we presented a full-vector finite difference method to solve for propagating modes supported by plasmonic waveguides of subwavelength size. We have used the Implicitly Restarted Arnoldi method to calculate the propagation constants for guided SPP modes. By directly solving for the dominant modes, the computational cost of the simulations has been dramatically reduced. The method can be applied to accurately model complex geometries and structures with fast varying field profiles. When applied to solve for purely bounded modes of plasmodic waveguide, our method automatically separates evanescent and low-loss propagating modes.

## References and links

**1. **M. Hochberg, T. Baehr-Jones, C. Walker, and A. Scherer, “Integrated plasmon and dielectric waveguides,” Opt. Express **12**, 5481–5486, (2004). [CrossRef]

**2. **K. Tanaka, M. Tanaka, and T. Sugiyama, “Simulation of practical nanometric circuits based on surface plasmon polariton gap waveguide,” Opt. Express **13**, 256–266, (2005). [CrossRef]

**3. **G. I. Stegeman, R. F. Wallias, and A. Maradudin, “Excitation of surface polaritons by end-fire coupling,” Opt. Lett. **8**, 386-(1983). [CrossRef]

**4. **R. Charbonneau, P. Berini, E. Berolo, and E. Lisicka-Shrzek, “Experimental observation of plasmon-polariton waves supported by a thin metal film of finite width,” Opt. Lett. **25**, 844–846, (2000). [CrossRef]

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

**6. **S. J. Al-Bader, “Optical transmission on metallic wires-fundamental modes,” IEEE J. Quantum Electron **40**, 325–329, (2004). [CrossRef]

**7. **P. Berini, “Plasmon-polariton waves guided by thin lossy metal films of finite width: bound modes of symmetric strucutres,” J. Phys. Rev. B **61**, 10484–10503, (2000). [CrossRef]

**8. **P. Berini, A. Stohr, K. Wu, and D. Jager, “Normal mode analysis and characterization of an InGaAs/GaAs MQW field-induced optical waveguide including electrode effects,” J. Lightwave Technol. **14**, 2422–2435, (1996). [CrossRef]

**9. **R. Zia, M. D. Selker, and M. L. Brongersma, “Leaky and Bound Modes of Surface Plasmon Waveguide,” Phys. Rev. B **71**, 165431, (2005). [CrossRef]

**10. **C. Chen, P. Berini, D. Feng, S. Tanev, and V. Tzolov, “Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media,” Opt. Express **7**, 260–272, (2000). [CrossRef]

**11. **J. A. Dionne, L. A. Sweatlock, and H. A. Atwater, “Planar metal plasmon waveguides: frequency-dependent dispersion, propagation, localization, and loss beyond the free electron model,” J. Phys. Rev. B **72**, 075405, (2005). [CrossRef]

**12. **I. El-Kady, M. M. Sigalas, R. Biswas, K. M. Ho, and C. M. Soukoulis, “Metallic photonic crystals at optical wavelength,” J. Phys. Rev. B **62**, 15299–15302, (2000). [CrossRef]

**13. **P. Lusse, P. Stuwe, J. Schule, and H.-G. Unger, “ Analysis of vectorial mode fields in optical waveguides by new finite difference method,” J. Lightwave Technol. **12**, 487–494, (1994). [CrossRef]

**14. **K. Ramm, P. Lusse, and H.-G. Unger, “Multigrid Eigenvalue Solver for Mode Calculation of Planar Optical Waveguides,” IEEE Photonics Technol. Lett. **9**, 967–969, (1997). [CrossRef]

**15. **V. Hernandez, J. E. Roman, A. Tomas, and V. Vidal, “A Survey of Software for Sparse Eigenvalue Problems,” Technical report, Universidad Politecnica de Valencia, (2005).

**16. **W. J. Stewart and A. Jennings, “Algorithm 570: LOPSI: A Simultaneous Iteration Method for Real Matrices [F2],” ACM Trans. Math. Softw. **7**, 230–232, (1981). [CrossRef]

**17. **R. B. Lehoucq and D. C. Sorensen, “Deflation techniques within an implicitly restarted iteration,” SIAM Journal on Matrix Analysis and Applications **17**, 789–821, (1996). [CrossRef]

**18. **D. C. Sorensen, “Implicit application of polynomial filters in a K-step Arnoldi method,” SIAM Journal on Matrix Analysis and Applications **13**, 357–385, (1992). [CrossRef]

**19. **R. Radke, “A MATLAB implementation of the implicitly restarted Arnoldi method for solving large scale eigenvalue problems,” Technical report, Dept. of Applied and Computational Mathematics, Rice University, Houston, TX, (1996).

**20. **R. Zia, A. Chandran, and M. L. Brongersma, “ Dielectric waveguide model for guided surface polaritons,” Opt. Lett. **30**, 1473–1475, (2005). [CrossRef]