Solving dielectric and plasmonic waveguide dispersion relations with a pocket calculator

We present a robust iterative technique for solving complex transcendental dispersion equations routinely encountered in integrated optics. Our method especially befits the multilayer dielectric and plasmonic waveguides forming the basis structures for a host of contemporary nanophotonic devices. The solution algorithm ports seamlessly from the real to the complex domain--i.e., no extra complexity results when dealing with leaky structures or those with material/metal loss. Unlike several existing numerical approaches, our algorithm exhibits markedly-reduced sensitivity to the initial guess and allows for straightforward implementation on a pocket calculator.


Introduction
With an ever-improving ability to fabricate miniature functional components, the field of nanophotonics is in a phase of explosive growth. Not only are novel phenomena being routinely observed in wavelength-and sub-wavelength-scale components, but these concepts are being rapidly translated into exotic devices for applications ranging from solar energy and information technology, to biology [1,2,3]. At the heart of this development is the ability of nanostructures to confine, guide, and scatter light in ways that are not achievable with bulk materials. It follows that understanding and improving the performance of these nanostructures requires a detailed knowledge of the electromagnetic modes they support.
The allowed electromagnetic modes in microphotonic structures are determined by solving Maxwell's equations for the given geometry. The last step in this procedure is the application of boundary conditions at the interfaces, yielding the dispersion equations which must be solved to obtain the allowed modes. These dispersion equations are in general transcendental and except in a few simple cases, their solutions cannot be expressed in terms of elementary mathematical functions.
When the solutions of the dispersion equations are known to be real, they are usually determined using a graphical search algorithm such as the bisection method [4]. However, in lossy material systems or low index-contrast asymmetric waveguides, the dispersion equations have complex solutions [5]. This calls for a search in two dimensions (i.e., in the complex plane) and severely limits the effectiveness of graphical algorithms.
Material loss and waveguide leakage are even more commonly encountered in plasmonic waveguides [1]. Due to the growing technological and scientific clout of plasmonic systems [6,7,8,9,10,11,12,13]various techniques for solving plasmon waveguide dispersion equations are in current use. Examples include the reflection-pole method (RPM) [14,15], Newton's method [4], and the argument-principle method [16]. Practical implementation of these methods requires careful programming customized for solving the problem at hand. In the methods that rely on curve-fitting, increasing the accuracy beyond a few decimals is often a challenging task, often requiring one to iterate the procedure manually based on previous results. Additionally, the root-finding algorithms of several commercial mathematical software suites (e.g., Mathematica™) can be very sensitive to initial guesses provided by the user; estimation of this initial guess is especially difficult when the solution is complex.
In this paper, we present an easy-to-implement iterative procedure for solving complex transcendental dispersion equations that is relatively insensitive to initial guess. Our method applies to rectangular multilayer dielectric and plasmonic waveguides which may have either material or leakage loss. We first use a simple numerical example to illustrate the use of this technique Fig. 1. Graphs of the left-(red) and the right-(blue) hand sides of Eq. (4). The red dot indicates the approximate location of the solution near x 0.5. and its convergence behavior. We then successively apply the procedure to find modes of dielectric slab waveguides, photonic wire waveguides, and plasmonic waveguides.
Papers in the past have discussed the use of iterative techniques for solving transcendental equations in general [17]. Our aim in this paper is to demonstrate how this method can be applied to automate the design and analysis of waveguide structures of significant contemporary interest.

Iterative technique: A simple example
Although the general philosophy of iterative methods is well known, we will use a simple example to introduce our terminology, illustrate the procedure use, and build an intuitive understanding of convergence issues. Suppose we need to solve the equation: where F(x) is a combination of elementary mathematical functions. We reformulate this as an iterative problem by converting Eq. (1) to: where f is obtained by manipulating the parent function F. Next, we start with an initial guess x 1 and obtain a sequence according to: If the sequence defined by Eq. (3) converges, then the limiting value is the solution to Eq. (2). We will illustrate the method by a simple example. Consider the transcendental equation Fig. 1 shows the graphs of the left-hand side (LHS) and right-hand side (RHS) of Eq. (4), which intersect around x 0.5. To find the root more accurately by the iterative method, we recast the equation in a form similar to Eq. (2) by choosing f (x) = 1 − sin x, setting up an iteration scheme following Eq. (3): the procedure converges to the intersection point of the two curves, starting from the initial guess. Making an initial guess corresponds to choosing a point on the curve y = x shown by the red curve in Fig. 2(a). This point is indicated in Fig. 2(a) by the red dot. The vertical lines in the spiral represent computation of f (x) for the chosen x and the horizontal lines represent re-substitution of x by f (x) for the next iteration. Before we proceed to apply the iterative technique to practical waveguide problems, it is useful to gain an intuitive understanding of how the above scheme converges to the solution. This is important since, for any given F(x), there are usually multiple ways to choose f (x) which differ in their convergence behavior. For example, an equally-legitimate way of setting up an iteration scheme for Eq. (4) would be to choose f (x) = sin −1 (1 − x) and obtain the sequence {x n } according to: The behavior of successive iterates is shown in Fig. 2(d-f). Regardless of how close the initial guess is to the actual solution, the iterative scheme diverges (spirals away) from the intersection point. Although the real part of the solution appears to converge after about 25 iterations, the overall complex solution does not converge, as seen from the undamped oscillations in the imaginary part of the solution. Thus, this iteration scheme, although derived from the same parent equation, does not converge. This behavior is typical for transcendental equations involving trigonometric functions and is therefore encountered for waveguide problems, as we will show in section 4.2.
On the other hand, convergence of Eq. (4) can be improved considerably by choosing an alternative iterative form as follows. Squaring both sides of Eq. (4) yields: from which we can obtain yet another sequence of {x n } of iterations according to: The preceding examples illustrate the crucial importance of choosing appropriate iterative forms. Arriving at such a form necessitates an understanding of the convergence of the iterative scheme. Figure 3 pictorially shows the criterion for convergence of the iterative scheme in Eq. (2). Whether the iteration spirals toward the intersection point (the solution) or away from it is governed by the relative slopes of the intersecting curves. If f (x) is steeper than x (i.e., if | f (x)| > 1), the successive computations and re-substitutions spiral away from the intersection and the solution diverges. On the other hand, if f (x) rises gently compared to x (i.e., if | f (x)| < 1), then the scheme spirals toward the intersection and the solution converges. For convergent functions, the rate of convergence is governed by the magnitude of | f (x)|. The RHS of Eq. (8) is much flatter (| f (x)| 1) near the solution compared to the RHS of Eq. (5), which results in the latter's considerably slower convergence. While the preceding justification is not a rigorous proof, it provides a basis for choosing the manipulations required to obtain convergent forms of the practically-useful equations we consider in the forthcoming sections [18].

Dispersion equation of a general asymmetric three-layer structure
The design of many integrated photonics systems of practical importance-including both novel plasmonic waveguide architectures as well as conventional dielectric waveguides-can often be reduced to solving for the effective mode index and electromagnetic field distributions in a three-layer planar structure. As such, we focus our analysis on a generalized three-layer slab waveguide, as shown in Fig. 4(a). The central layer is called the core and has a complex relative permittivity ε f . The bottom and top cladding layers are called the substrate (with permittivity ε s ) and cover (with permittivity ε c ), respectively.
Depending on which materials comprise the three layers, it is possible to classify the most commonly-encountered waveguides into three basic types. The dielectric waveguides, shown in Figures 4(e-i) schematically show the possible modes of a general three-layer slab waveguide. These modes are distinguished according to several criteria: mode number (number of zero-crossings in the field), polarization (transverse electric or transverse magnetic), field symmetry (even or odd), and field profile in the core (sinusoidal or hyperbolic). Of the possible field profiles shown in Fig. 4, the mode shown in (e), having sinusoidal core-field variation with no zero-crossings, is exclusive to dielectric waveguides and is known as the fundamental mode. Modes (h) and (i), having hyperbolic core-field variation, are exclusive to plasmonic waveguides and are known as the gap-plasmon modes. Modes (f) and (g), having sinusoidal core-field variation with ≥ 1 zero-crossings, are common to both dielectric and plasmonic waveguides.
This rich variety of modes arises out of the solutions of the dispersion equation written for an asymmetric three-layer slab waveguide. In principle, a single dispersion equation can completely describe all the modes of both dielectric and plasmonic waveguides. However, because the nominally-assumed core-field distributions for the two types are different (sinusoidal for dielectric waveguides and hyperbolic for plasmonic waveguides), we choose to write two separate equations for the two types [19, 20]: For convenience, we have collected the symbol definitions and their expressions in Table 1.
Using these definitions, Eqs. (9) and (10) can be cast explicitly as transcendental equations in a single complex variable k or κ. Furthermore, it is evident that Eq. (9) transforms to Eq. (10) for purely imaginary values of k (i.e., k → iκ, κ ∈ R). Having solved for k or κ, the mode's complex effective index n eff may be calculated using relations in Table 1.  We will now begin the description of the actual iterative solution procedure. We saw in section 1 that the convergence of this technique depends crucially on the iteration function. As a result, the iterative form for determining modes with sinusoidal core-fields [ Fig. 4(e-g)] is different from the one needed for determining modes with hyperbolic core-fields [ Fig. 4(hi)]. Therefore, we split our description into two broad categories: dielectric waveguides and plasmonic waveguides. These two categories are further divided depending upon the degree of confinement (strong or weak) for dielectric waveguides, and core/cladding type (MDM or DMD) for plasmonic waveguides. Figure 5 depicts the division of the waveguide modes that we have made for the purposes of describing our iterative solution process. In the forthcoming sections, we will obtain and test rapidly-convergent iterative forms for calculating the mode indices of these four sub-categories of three-layer asymmetric slab waveguides. These useful forms will be boxed for the convenience of the reader.

Strong-confinement dielectric waveguides
Using the half-angle identity tan(z/2) = − cot z ± √ 1 + cot 2 z, Eq. (9) can be transformed to: Fig. 6. Iterative solution of a strong-confinement waveguide using Eq. (12). (a) The convergence plot depicting the normalized LHS (magenta) and RHS (gray) of Eq. (12); f (k) on the y-axis of of (a) refers to the function on the right hand side. The convergence of the real and imaginary parts of the effective index are shown in (b) and (c), respectively.
Inverting the tangent function in Eq. (11) gives us the convergent iterative form for strongconfinement waveguides as: Here M is an integer and the quantities γ cn , γ sn , G cn , and G sn are related to the n th iterate k n through the relationships in Table 1. Although the form of Eq. (12) may seem daunting, we note that it is a general equation which can calculate odd and even modes of both transverse electric (TE) and transverse magnetic (TM) polarizations for asymmetric slab waveguides. All of the functions involved in Eq. (12) are found on scientific calculators and in mathematical software tools. The iterative procedure obviates any need for a graphical search and has an inherent ability to give arbitrary-precision solutions-a feat difficult to achieve using graphical interpolation. Before embarking on actual computation, we need to specify the parity (+ or − sign), order (M), and the polarization (p and q) of the desired mode. For TE polarization (p = q = 1) and even parity (+ sign), the effective indices of successive even TE modes may be simply calculated by setting M = 0, 1, . . . m and iterating according to Eq. (12); this yields the TE 0 , TE 2 , . . . TE 2m modes. For TE polarization with odd parity (− sign), setting M = 1, 2, . . . , m yields the TE 1 , TE 3 , . . . TE 2m−1 modes. Note that for even modes, the mode order begins with M = 0, whereas for odd modes, it begins with M = 1. The procedure for obtaining TM mode indices is identical, with proper input of (p, q) as shown in Table 1.
To illustrate the use of the iterative method, we choose a 1µm-thick silicon-on-insulator (SOI) slab waveguide operating at 1550 nm as a test structure. The refractive indices of the silicon film, oxide substrate, and air cover are assumed to be √ ε f = 3.50, √ ε s = 1.45, and √ ε c = 1.00 respectively. Using the TE 0 mode as an example, Fig. 6 depicts how Eq. (12) converges to the solution. Fig. 6(b) and (c) show the real and imaginary values of the first twenty iterates of Eq. (12). Although the iteration is carried out in terms of k, we have chosen to depict the convergence in terms of the effective index n eff = β /k 0 since this is a more familiar quantity in waveguide analysis. Fig. 6(a) shows the LHS and RHS of Eq. (12) and how the iterative scheme converges to the solution. Notice that the slope condition mentioned in section 1 is satisfied for this particular case.
To examine the iterative technique in contrast with standard numerical techniques, we solved for the modes of the same structure using Mathematica™. Table 2 compares the two solutions. Mathematica™ calculations were performed by feeding Eq. (11) to the FindRoot function with an initial guess dependent on the mode order, parity, and polarization. To obtain the initial guesses, we plotted the left-and right-hand sides of Eq. (11) Table 1 were initiated with the same initial guess k 1 = K s . Moreover, the convergence of the iterative technique is completely insensitive to the exact value of this initial guess: we obtained the same effective indices for initial guesses ranging from K s to 10 5 K s (including complex values). While a precise mathematical characterization of the convergence behavior is outside the scope of our paper, this exercise does suggest the robustness of the iterative technique with respect to the accuracy of the initial guess.

Symmetric strong-confinement waveguides
If the relative permittivities of the cover and the substrate are equal, then the structure is called a symmetric waveguide. Since several important photonic device structures employ symmetric waveguides, we now give explicit iterative forms for calculating their mode indices. These iterative forms arise out of the simplifications in Eq. (12) due to the equality of ε s and ε c : We also note that, like Eq. (12), M assumes values starting from 0 for even modes and 1 for odd modes.

Weak-confinement dielectric waveguides
Eq. (12) converges for a very wide range of refractive indices, wavelengths, and waveguide thicknesses which are commonly used in practical integrated photonics designs. It performs poorly, however, when applied to the weak-confinement single-mode waveguide structures routinely fabricated out of III-V (e.g., GaAs/AlGaAs) or organic materials. In the context of multimode waveguides, the strong-confinement formula encounters convergence problems when used to calculate indices of modes near the waveguide cutoff. We highlight this limitation and its remedy through the following example. Consider a 1µm-thick air-clad GaAs waveguide on an Al 0.1 Ga 0.9 As substrate operating at 1550 nm. The refractive indices of GaAs and Al 0.1 Ga 0.9 As are assumed to be √ ε f = 3.300 and √ ε s = 3.256 respectively. We attempt to find the fundamental TE mode iteratively by inputting the corresponding parameters to Eq. (12) (M = 0, p = q = 1, and + sign). Fig. 7(b) and (c) show the behavior of the real and imaginary parts of the calculated effective index for first 50 iterates. The real part of the effective index initially diverges to reach a maximum at iteration number 30 and begins to converge thereafter. Starting from the thirtieth iteration, the imaginary part oscillates between ±0.01. These oscillations do not damp out even if we increase the number of iterations to 10 4 . Mathematica™, however, indicates that this structure supports TE 0 and TM 0 modes with indices of 3.266 and 3.263, respectively. We conclude that the iterative scheme in Eq. (12) fails to converge for this example. This failure to converge can be understood by examining the behavior of the LHS and RHS of Eq. (12) near the intersection point, as shown in Fig. 7(a). Near the solution, | f 1 (κ)| > 1. As a result, the iterates begin to diverge away from the intersection point. In fact, | f 1 (κ)| continues to increase away from the intersection, causing a rapidly-increasing divergence. However, at κ 0.537k 0 , | f 1 (κ)| abruptly becomes < 1, causing the real part of the iteration to converge. This convergence of the real part is misleading, however since the imaginary part shows undamped oscillations. The behavior of the strong confinement formula in this case is similar to the example problem in section 1 illustrated in Fig. 2(d-f).
In summary, the basic cause of the failure of convergence is the breakdown of the slope condition mentioned in section 1. This is remedied by recasting Eq. (9) into a form which satisfies the convergence condition. To this end, we use the identity tan 2 z = 1/ cos 2 z − 1, rewriting Eq. (9) as: The convergent form suitable for weak confinement is obtained by solving for k contained in the γ c γ s term on the right hand side of Eq. (14): where G cn and G sn are related to k n through the relationships in Table 1. We emphasize that Mode (p, q) Iterative scheme Mathematica™ TE 0 (1,1) 3.26599646645606654 3.26599646645606654 TM 0 (ε f /ε c , ε f /ε s ) 3.26338400537407312 3.26338400537407312 Table 3. Comparison of mode indices for an Al 0.1 Ga 0.9 As/GaAs/air slab waveguide computed via the iterative method and the FindRoot function in Mathematica™ Eq. (15) is a general equation capable of solving for near-cutoff TE and TM modes of odd and even parity in three-layer asymmetric slab waveguides. The choice of + or − sign in Eq. (15) depends on whether the cutoff occurs at the odd (+) or the even (−) mode. When analyzing a given structure, it is difficult to determine a priori the parity of the cutoff mode. The prescription is therefore to use the strong-confinement formula until it fails to converge, continuing with the weak-confinement formula thereafter. The domain of convergence of the strong-and weakconfinement formulas may be determined rigorously by setting the derivatives of the RHS of Eqs. (12) and (15) equal to 1 and solving for the corresponding value of k. Solving for the mode indices of the previously-considered GaAs/Al 0.1 Ga 0.9 As waveguide using Eq. (15) yields the effective indices of the TE 0 and TM 0 modes in excellent agreement with the Mathematica™ results, as seen from Table 3. The convergence behavior of the effective index is plotted in Fig. 7(d-f). The convergence of Eq. (15) improves as the core-cladding index contrast decreases; similarly, the convergence of Eq. (12) improves for increasing corecladding index contrast.

Symmetric weak-confinement dielectric waveguides
Similarly to strong-confinement waveguides, having identical cover and substrate indices considerably simplifies the iterative form for weakly-confined slab waveguides, which can be written as:

Extension to photonic wire waveguides
Waveguides in which light is confined in two dimensions and propagates in the third dimension are known as photonic wire (PW) waveguides; optical fibers, as well as ridge, rib, and channel waveguides are examples. Because they offer stronger light confinement compared to slab waveguides, PW waveguides are commonly employed in on-chip optical devices where a small size is essential for achieving several critical device specifications ( e.g., speed, low power, integration density, etc.). Although the exact field distribution in PW waveguides is best determined by numerical calculations, the effective index method (EIM) is remarkably successful at quickly calculating the mode indices for several common PW waveguide configurations [19]. A knowledge of the effective index suffices for many important design problems. Figure 8 illustrates the two steps in the implementation of EIM. Figure 8(a) shows an idealized silica-clad SOI photonic wire waveguide having a width w and thickness h. In the first step, we disregard the confinement in the x-direction and solve for the effective index n . Depending on the polarization of the desired mode, we will need to solve either the TE or TM equation. In the second step, we consider a silica-clad slab waveguide with a core index n confined in the x-direction. Note that the orientation of the effective waveguide is orthogonal to that in the first step and hence the equation in this step is for the polarization orthogonal to that in the previous step. That is, if the first step uses the TE dispersion equation, then the second step uses the TM equation(and vice versa). The effective index n determined in this second step is our final answer for the effective mode index of the original PW waveguide.
As an example, suppose that we need to determine the effective index of the fundamental E y -polarized mode of a sub-micron PW waveguide shown in Fig. 9 (a) with h = 450 nm and w = 300 nm. In the first step, we solve for the TE mode equation for the structure shown in Fig. 8(b). Since SOI waveguide constitutes a strong-confinement system, the first step is readily accomplished by iteration of Eq. (12) for the TE condition (namely, M = 0, p = q = 1). This yields the first effective index as n = 3.073930677459340. In the second step, we consider an xconfined slab with core-index n and iterate Eq. (12) with the fundamental TM mode condition. This yields the effective index of the PW waveguide as n = 2.652766507502340.
To test the relative accuracy of the iterative method, we use the finite element method (FEM) to compute the mode index. For the waveguide with the above dimensions, FEM gives n = 2.612594. It is evident that even for sub-mircron dimensions, the relative error defined as   n iter − n FEM   /n FEM is only 1.5%. Figure 9(b) shows the variation of the relative error with waveguide dimensions for a fixed width-to-height aspect ratio of 1.5, from which we note the exponential drop in the relative error with increasing waveguide size. The reason for the decrease of accuracy with decreasing size lies in the increased interaction of the electromagnetic fields with the corner regions of the waveguide. By decomposing a two-dimensional mode-solving problem into two one-dimensional problems, we effectively chose to neglect the interaction of the fields with the waveguide corners. This assumption becomes decreasingly accurate with decreasing waveguide size. Figure 9(c) and (d) show the y-component of the electric field, calculated using FEM and the iterative method, for the sub-micron waveguide considered above. Note the strong electric field at the waveguide corners in Fig. 9(c) calculated using FEM, and its absence in Fig. 9(d) calculated using the iterative method. In addition to highlighting this difference for small waveguides, this exercise illustrates the well-known fact that while EIM successfully calculates the effective indices, it does not guarantee satisfaction of the boundary conditions. It should therefore not be used to calculate fields for wavelength-sized structures.
In spite of this limitation, the EIM remains an intuitive way of accomplishing designs which rely primarily on the knowledge of mode indices. Examples of such designs include several contemporary device structures such as SOI ring resonators, Bragg gratings, Mach-Zehnder interferometers, and directional couplers. By offering a simple means of solving the slab waveguide dispersion equation, the iterative method greatly enhances the efficacy and implementation speed of EIM.

Modes of plasmonic waveguides
The examples in the preceding sections indicate the relative simplicity of the iterative method when applied to dielectric waveguides. However, the true advantage of this technique emerges when dealing with systems having material loss (complex permittivity) or leakage loss (complex propagation constant). The modal indices are complex in both these cases and providing graphical root-finders with a complex initial guess becomes difficult. As such, the selfconverging behavior of the iterative algorithm becomes an invaluable asset.
Propagation loss and waveguide leakage arise routinely in sub-wavelength-scale metallic waveguides. Intense research efforts are currently underway in the field of metal-based optics, also known as plasmonics (see, e.g., [1] and references contained therein). Extended metal structures support surface plasmon-polariton (SPP) modes that are electromagnetic waves strongly coupled to collective electron oscillations in the metal. To exploit the strong light localization achievable in plasmonic structures, a variety of waveguide configurations have recently been proposed and demonstrated [1]. However, because of its simplicity, the MDM geometry remains a canonical structure for achieving and studying strong light confinement. As such, a convenient technique for determining the optical modes of MDM waveguides is highly valuable.
We begin the description of the iterative solution technique for plasmonic waveguides by referring to the division of slab waveguide types illustrated in Fig. 4(c) and (d). We first consider the modes of an MDM-type plasmonic waveguide followed by a treatment of the DMD-type.

Metal-dielectric-metal waveguides
A MDM waveguide supports gap-plasmon and TM-like waveguide modes. Even and odd gapplasmon modes are practically important, as they offer the strongest sub-wavelength field confinement. However, TM-like modes are also of theoretical value and are necessary for gaining a complete understanding of reflection and transmission phenomena in MDM waveguides and antennas [16]. Effective indices for both these types of modes can be conveniently determined using the iterative technique, as we show in the following.

Gap-plasmon modes
We start with our master plasmonic Eq. (10) and rearrange it as: Here, α c ,α s , and S are related to κ as specified in Table 1. Treating this as a quadratic equation in κ and solving yields the convergent iterative form for the MDM gap-plasmon modes as: This simple-looking equation can calculate even (+ sign) and odd (− sign) gap-plasmon mode indices of a wide variety of deep sub-wavelength asymmetric plasmonic waveguides. We illustrate the use of Eq. (18) through two examples. Our first structure is a 50-nm-thick goldsilica-silver slab waveguide operating at 1550 nm. The relative permittivities of gold, silica, and silver are assumed to be −95.92 − i10.97, 2.1025, and −143.49 − i9.52, respectively. The effective index of the even gap-plasmon mode, calculated using the + sign in Eq. (18), is given in Table 4. Fig. 10 shows plots of the LHS and RHS of Eq. (18) and the convergence of the real and imaginary parts of the effective index. The thin 50-nm slab considered above does not support an odd (antisymmetric) gap-plasmon mode. However, increasing the silica thickness to 3 µm allows the structure to support gapplasmon modes of both symmetries. We calculate the effective indices iteratively using Eq. (18) with + and − signs for odd and even gap-plasmon modes respectively. Once again, the iterative method agrees with the calculations of the FindRoot function in Mathematica™ (Table 4).

TM-like waveguide modes
The other class of modes supported by a MDM waveguide are the TM-like waveguide modes with profiles as shown in Fig. 4(f) and (g). For these modes, κ is purely imaginary and plasmon Eq. (10) transforms to dielectric Eq. (9). The iterative form for obtaining the indices of these modes is identical to Eq. (12), with only a few slight differences as regards its implementation. For the case of dielectric waveguides, the mode index M assumed integer values starting from 0 for even modes and 1 for odd modes. For MDM waveguides, this is reversed: M assumes integer values starting from 1 for even modes and 0 for odd modes. This is a consequence of the signs of p and q being reversed for MDM waveguides due to the negative permittivity of the metal "claddings." Additionally, because of the large index difference between metals and dielectrics, the TM-like modes of MDM waveguides are almost always calculated using the strong-confinement formula in Eq. (12).
To show how TM-like modes are calculated for MDM waveguides, we choose a 300-nm thick silica layer sandwiched between semi-infinite gold and silver layers. The permittivities of   each material are the same as assumed in the previous subsection. The indices of the first five TM-like modes, calculated using the strong-confinement formula Eq. (12), are listed in Table 4 alongside the corresponding solutions from Mathematica™.

Symmetric MDM waveguides
Many contemporary high-confinement architectures employ the symmetric MDM waveguide as their skeleton structure. Equality of the relative permittivities of the substrate and cover further simplify Eq. (18) for the gap-plasmon modes. We can express the resulting iterative forms as: . . . for even gap plasmon. (19a) κ n+1 = −p κ 2 n + K 2 c coth (κ n h/2) . . . for odd gap plasmon. (19b)

Dielectric-metal-dielectric waveguides
The second type of plasmonic waveguide structure commonly encountered in practice is an DMD waveguide, as shown in Fig. 4 where A and B are defined in Table 1. Considering this as a quadratic equation in A and B leads us to the desired iterative forms. We start by identifying the lower-index dielectric as the substrate and using initial guesses of A 0 = k 0 , B 0 = 0, and κ 0 = k 0 . We then iterate to obtain five different quantities successively, using the following equations in the order shown: a * n = −κ n coth κ n h ± B 2 n + κ 2 n csch 2 (κ n h), b * n = a * 2 n + κ 2 n + 2a * κ n coth 2 (κ n h), (21b) A n+1 = (pξ c,n+1 + qξ s,n+1 ) /2, (21d) B n+1 = (pξ c,n+1 − qξ s,n+1 ) /2.
Here, a * n , b * n are intermediate dummy parameters and ξ c,n+1 and ξ s,n+1 are related to κ n+1 through Table 1. The indices of the low-and high-energy plasmon modes are obtained by using the + and the − signs, respectively, in Eq. (21a). The large asymmetry in the decay constants in the metal core (κ) and dielectric claddings (ξ c , ξ s ) makes the determination of the effective indices of DMD modes a numerically-challenging problem. For the iterative procedure, this translates into a difficulty in obtaining convergent iteration expressions. Therefore, unlike the previous cases, the iterative procedure for DMD waveguides involves multiple steps instead of a single closed-form iteration function. Admittedly, this multi-step algorithm goes against the "pocket calculator" philosophy, but nevertheless is convenient and efficient compared to direct numerical solutions.
We will now illustrate the use of Eq. (21) through specific examples. Consider the operation of a 50-nm-thick silver film on a silica substrate at 1550 nm. The relative permittivities of silver and silica are the same as assumed previously; the top cover is air (ε c = 1). We start  Table 5. Effective indices of various modes of an DMD waveguide obtained using the iterative method and their comparison with direct solutions using Mathematica™.
by specifying A 0 , B 0 , and κ 0 , and obtain their successive values according to Eq. (21). The convergence of A n , B n , and κ n is shown in Fig. 11; Table 5 compares the solution obtained using the iterative scheme with the direct solution computed by the Mathematica™ FindRoot function. In our computations, the FindRoot function was unable to return the complex mode index unless we gave a very precise initial guess for both the real and the imaginary parts. On the other hand, the iterative method computed the complex mode index regardless of the initial guess. Our next example is a symmetric structure with a 100-nm silver film sandwiched between two silica layers. Although 100 nm is at the boundary of the film behaving like an DMD waveguide versus two separate interfaces, we choose these dimensions to highlight the robustness of the iterative technique even for the most extreme cases of root-finding. For this structure, both high-and low-energy modes exist, whose indices are conveniently calculated by using the − and + signs, respectively, in Eq. (21). The indices calculated using the iterative method and Mathematica™ for this example are included in Table 5 and show good agreement. Calculating these indices was especially difficult using FindRoot because of their close numerical proximity. The iterative technique, on the other hand, specifies separate functions (the + and− signs) which are guaranteed to converge to these different modes, irrespective of their numerical proximity.

Symmetric DMD waveguides
Symmetric DMD waveguides appear in the form of idealized waveguide geometries such as metal films in homogenous dielectric media (including air). Fortunately, for the symmetric case, the iterative scheme consists of a single equation which can be written as: . . . for odd plasmon mode.
These equations do not follow automatically from Eq. (21) but instead have to be derived separately by considering the dispersion equation for the symmetric DMD waveguide. Eq. (22) are preferable to Eq. (21) for analyzing symmetric DMD structures owing to better convergence behavior and an obvious ease in programming.

Conclusion
We have presented a robust and an easy-to-implement iterative technique for determining complex propagation constants of asymmetric dielectric and plasmonic waveguides. At the heart of our procedure are the closed-form iteration functions, namely the boxed Eqs. (12), (15), (18), and (21). In addition to the programming ease, the iterative technique has an inherent ability to give arbitrary-precision answers-a feat difficult to achieve using graphical or curve-fitting algorithms such as the reflection-pole method. Because of its insensitivity to initial guess, it is our hope that this technique will help facilitate design tasks that require rapid and automated calculation of mode indices for a variety of nanophotonic structures having a rectangular geometry.