Analysis of long-range surface plasmon polaritons in nonlinear plasmonic waveguides using pseudospectral method

A full-vectorial pseudospectral method is reported for solving the mode characteristics of nonlinear dielectric and plasmonic waveguides. The coupled equations are formulated in terms of transverse magnetic-field components, and self-consistent solutions are obtained through an iterative procedure. The proposed scheme applies in a saturable medium with biaxial anisotropy of practical interest. The accuracy and efficiency of this scheme are demonstrated by solving the mode bistability of a nonlinear dielectric optical waveguide, analyzed by the well-known finite-element-methodbased imaginary-distance beam propagation method. Furthermore, the relationship between geometry and input power is studied by analyzing the power dispersion curve of the long-range surface plasmon polariton modes of a nonlinear plasmonic waveguide. ©2012 Optical Society of America OCIS codes: (130.2790) Guided waves; (230.7380) Waveguides, channeled; (190.4390) Nonlinear optics, integrated optics; (240.4350) Nonlinear optics at surfaces. References and links 1. T. H. Wood, “Multiple quantum well (MQW) waveguide modulators,” J. Lightwave Technol. 6(6), 743–757 (1988). 2. S. Radic, N. George, and G. P. Agrawal, “Optical switching in λ/4-shifted nonlinear periodic structures,” Opt. Lett. 19(21), 1789–1791 (1994). 3. Y. D. Wu, “New all-optical wavelength auto-router based on spatial solitons,” Opt. Express 12(18), 4172–4177 (2004). 4. T. Fujisawa and M. Koshiba, “All-optical logic gates based on nonlinear slot-waveguide couplers,” J. Opt. Soc. Am. B 23(4), 684–691 (2006). 5. A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar, “Nonlinear plasmonic slot waveguides,” Opt. Express 16(26), 21209–21214 (2008). 6. I. D. Rukhlenko, A. Pannipitiya, and M. Premaratne, “Dispersion relation for surface plasmon polaritons in metal/nonlinear-dielectric/metal slot waveguides,” Opt. Lett. 36(17), 3374–3376 (2011). 7. I. D. Rukhlenko, A. Pannipitiya, M. Premaratne, and G. Agrawal, “Exact dispersion relation for nonlinear plasmonic waveguides,” Phys. Rev. B 84(11), 113409 (2011). 8. A. R. Davoyan, I. V. Shadrivov, A. A. Zharov, D. K. Gramotnev, and Y. S. Kivshar, “Nonlinear nanofocusing in tapered plasmonic waveguides,” Phys. Rev. Lett. 105(11), 116804 (2010). 9. J. R. Salgueiro and Y. S. Kivshar, “Nonlinear plasmonic directional couplers,” Appl. Phys. Lett. 97(8), 081106 (2010). 10. A. Degiron and D. R. Smith, “Nonlinear long-range plasmonic waveguides,” Phys. Rev. A 82(3), 033812 (2010). 11. K. Hayata and M. Koshiba, “Full vectorial analysis of nonlinear-optical waveguides,” J. Opt. Soc. Am. B 5(12), 2494–2501 (1988). 12. R. D. Ettinger, F. A. Fernandez, B. M. A. Rahman, and J. B. Davies, “Vector finite element solutions of saturable nonlinear strip-loaded optical waveguides,” IEEE Photon. Technol. Lett. 3(2), 147–149 (1991). 13. X. H. Wang and G. K. Gambrell, “Full vectorial simulation of bistability phenomena in nonlinear-optical channel waveguides,” J. Opt. Soc. Am. B 10(6), 1090–1095 (1993). 14. S. S. A. Obayya, B. M. A. Rahman, K. T. V. Grattan, and H. A. E. Mikati, “Full vectorial finite–element solution of nonlinear bistable optical waveguides,’,” IEEE J. Quantum Electron. 38(8), 1120–1125 (2002). 15. K. Hayata and M. Koshiba, “Full vectorial analysis of nonlinear-optical waveguides,” J. Lightwave Technol. 20, 1876–1884 (2002). 16. J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springer–Verlag, 2001). #171463 $15.00 USD Received 28 Jun 2012; accepted 23 Jul 2012; published 31 Jul 2012 (C) 2012 OSA 13 August 2012 / Vol. 20, No. 17 / OPTICS EXPRESS 18665 17. C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron. 11(2), 457–465 (2005). 18. P. J. Chiang, C. P. Yu, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008). 19. C. C. Huang, “Simulation of optical waveguides by novel full-vectorial pseudospectral-based imaginary-distance beam propagation method,” Opt. Express 16(22), 17915–17934 (2008). 20. J. B. Xiao and X. H. Sun, “Full-vectorial mode solver for anisotropic optical waveguides using multidomain spectral collocation method,” Opt. Commun. 283(14), 2835–2840 (2010). 21. C. C. Huang, “Numerical investigation of mode characteristics of nanoscale surface plasmon-polaritons using a pseudospectral scheme,” Opt. Express 18(23), 23711–23726 (2010). 22. C. C. Huang, “Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method,” Opt. Express 19(4), 3363–3378 (2011). 23. C. C. Huang, “Pseudospectral mode solver for analyzing nonlinear optical waveguides,” Opt. Express 20(12), 13014–13029 (2012). 24. P. Berini, “Plasmon polariton modes guided by a metal film of finite width,” Opt. Lett. 24(15), 1011–1013 (1999). 25. P. Berini, “Plasmon-polariton modes guided by a metal film of finite width bounded by different dielectrics,” Opt. Express 7(10), 329–335 (2000). 26. P. Berini, “Plasmon-polariton modes guided thin lossy metal films of finite width: bound modes of symmetric structures,” Phys. Rev. B 61(15), 10484–10503 (2000). 27. P. Berini, “Plasmon-polariton modes guided thin lossy metal films of finite width: bound modes of asymmetric structures,” Phys. Rev. B 63(12), 125417 (2001). 28. S. J. Al-Bader, “Optical transmission on metallic wires-fundamental modes,” IEEE J. Quantum Electron. 40(3), 325–329 (2004). 29. A. Degiron and D. R. Smith, “Numerical simulations of long-range plasmons,” Opt. Express 14(4), 1611–1625 (2006). 30. P. Berini, “Long-range surface plasmon polaritons,” Adv. Opt. Photon. 1(3), 484–588 (2009). 31. G. I. Stegeman and C. T. Seaton, “Nonlinear surface plasmons guided by thin metal films,” Opt. Lett. 9(6), 235– 237 (1984). 32. J. Ariyasu, C. T. Seaton, G. I. Stegeman, A. A. Maradudin, and R. F. Wallis, “Nonlinear surface polaritons guided by meal films,” J. Appl. Phys. 58(7), 2460–2466 (1985). 33. R. F. Oulton, V. J. Sorger, D. A. Genov, D. F. P. Pile, and X. Zhang, “A hybrid plasmonic waveguide for subwavelength confinement and long-range propagation,” Nat. Photonics 2(8), 496–500 (2008).


Introduction
The intensity-dependent properties of nonlinear dielectric optical waveguides enrich propagation phenomena and induce effects that are entirely different than those observed for materials with linear permittivity responses.As a result, many functional photonic devices such as waveguide modulators [1], optical switches [2], wavelength auto-routers [3], and logic gates [4] are based on nonlinear mechanisms such as low threshold power, self-focusing, and bistability.In recent years, the dispersion relations of surface plasmon polaritons (SPPs) surrounded by a nonlinear dielectric medium [5][6][7][8][9][10] have been reported to offer additional factors for tailoring the mode characteristics of SPPs.In Ref [5], Davoyan et al. studied a metal-dielectric slot waveguide and predicted the symmetry-breaking bifurcation of symmetry modes with critical power that depends on the slot width.In Refs [6,7], Rukhlenko et al. derived an exact dispersion relation for SPPs using an exact field decomposition of transverse magnetic (TM) modes, and showed that it enables backward-propagating modes.For enhancing nanoscale nonlinear effects, Davoyan et al. [8] proposed nonlinear nanofocusing in a tapered plasmonic waveguide.In Ref [9], the authors proposed a nonlinear plasmonic directional coupler and showed how it substantially differs from a nonlinear dielectric coupler using the finite difference time domain (FDTD) approach.However, these studies [5][6][7][8][9] concerned only planar structures, not practical two-dimensional (2D) channel waveguides.Until now, thin metal strips surrounded by Kerr-type nonlinear materials [10] have been the only devices studied with regard to the mode properties of long-range (LR) SPPs.These studies used the commercial eigenmode solver (COMSOL Multiphysics), which is based on the finite-element method (FEM).
We know that no analytical solutions exist for practical 2D channel nonlinear waveguides.Therefore, it is essential to use an efficient and accurate computational scheme to properly design nonlinear waveguide components and integrated optical circuits.Most studies use the FEM [10][11][12][13] or the FEM-based imaginary-distance beam propagation method (IDBPM) [14,15] to obtain the power dispersion relation of nonlinear channel waveguides.Recently, the pseudospectral method (PSM) [16], which expands the mode fields by Lagrange-type orthogonal basis functions, has shown that fast convergence and high-order accuracy are possible for solving the mode characteristics of the linear permittivity response of many optical and plasmonic waveguides [17][18][19][20][21][22].These studies proved that the PSM requires significantly less computational memory and time than the FEM while achieving the same order of accuracy.In addition, the PSM is also used to study how the geometric parameters and power levels affect the dynamics of nonlinear SPP planar waveguides [23].Thus, the present study extends the previous work [23] to study the mode behaviors of practical 2D nonlinear waveguides.In particular, we investigate in detail how the geometric parameters and power levels of nonlinear plasmonic waveguides affect power dispersion.The remainder of this paper is organized as follows: Section 2 presents the wave equations for a 2D waveguide structure.Section 3 derives the computational schemes including the algorithm of pseudospectral method and the iterative approach.Section 4 demonstrates the accuracy of the proposed scheme by comparing the solution obtained by the proposed scheme with that obtained by other numerical schemes, and analyzes the mode features of nonlinear plasmonic waveguides.Finally, Section 5 presents the conclusions.

Mathematical formulations
The magnetic vector field (H) of a monochromatic light wave with time dependence of exp(jωt) propagating along the z direction obeys the wave equation ( ) where ω is the angular frequency, µ 0 is the vacuum permeability, and [ε] is the nonlinear permittivity tensor.For the biaxial nature of materials, the nonlinear permittivity tensor can be expressed as ( ) 0 0 ( , , ), ( , , ), where ε 0 is the vacuum permittivity, [ε r ] is the relative permittivity tensor, ε i is the linear relative permittivity, c 0 is the speed of light in vacuum, n is the nonlinearity coefficient, and the function f (E x , E y , E z ) expresses the flux-dependent behavior of the nonlinear medium.
Considering the stationary guided modes of the nonlinear medium, the magnetic-field components are assumed to be of the form exp(−jßz), where β = k 0 n e is the propagation constant along the z direction, 1 .
In the proposed scheme, we divide the computational window at interfaces between different materials into several subdomains with uniform or continuous refractive-index profiles, and then assemble these subdomains using physical interfacial conditions.The interfacial conditions used include continuous normal and tangential components of the magnetic fields H x and H y at each intra-element boundary, and the other two interfacial conditions are obtained from the divergence free condition of the magnetic field vector 0 ∇ ⋅ = H and the Ampere's law ∇ × = H jω[ε]E, where E is the electric field vector.To represent the mode fields with only H x and H y components, we express the condition 0 ∇ ⋅ = H as 1 , where H z is the longitudinal component of the magnetic field.For the condition Equations ( 6) and ( 7) are used as interfacial conditions according to the continuities of H z and E z at interfaces.To compute the flux-dependent behavior of the nonlinear medium, the electric-field components that are needed to compute the nonlinear relative permittivities in Eq. ( 3) can be derived by Maxwell's equations: where Z 0 denotes the intrinsic impedance in vacuum (Z 0 = 377 Ω).The input power along the z direction is given by where * denotes the complex conjugate.

Computational schemes
In the PSM, the computational window is divided into several subdomains at interfaces between different materials with uniform or continuous refractive-index profiles.In each subdomain, the transverse magnetic-field components are expanded as follows into a product of suitable basis functions comprising the Lagrange-type interpolations θ(x) in the x direction, ψ(y) in the y direction, and the H x and H y field components denoted H x,pq and H y,pq at the (n x + 1) × (n y + 1) grid (collocation) points:  Substituting Eqs.(10a) and (10b) into Eq.( 4), we find that Eq. ( 4) is satisfied exactly at the specific (n x + 1) × (n y + 1) collocation points determined by the chosen basis functions.The operators P xx , P xy , P yx , and P yy in Eq. ( 4) are (1) (1) (1) where ( ) ( ) indicate the hth-order derivatives of θ p (x) with respect to x and ψ q (y) with respect to y, respectively.After obtaining a matrix eigenvalue equation for each subdomain, a global matrix equation is then obtained by assembling all the matrix equations of all subdomains.Assuming that the computational window is divided into m subdomains, the pattern of the global matrix forms the following block-diagonal matrix: In Eq. ( 12), the linear equations at the interface points shared by more than two subdomains have to be removed by imposing physical interface conditions, including the continuous normal and tangential components of the magnetic fields H x and H y boundary, and the continuities of H z and E z in Eqs. ( 6) and (7), respectively.Accordingly, the number of grid points at each subdomain is reduced to (n x − 1) × (n y − 1).For brevity, we formulate the coupling of H x and H y using only two subdomains (labeled 1 and 2), as shown in Fig. 1(b).At interface boundary points [i.e., (1)   (2) 0 ( , ) ( , ), x n i i x y x y = where i = 0, 1, 2, …, n y , shown by the red circles] located in the vertical red line and shared by two subdomains, the condition of continuous normal and tangential components of the magnetic fields H x and H y is implicitly satisfied, and the continuity of H z in Eq. ( 6) gives where In Eqs. ( 14) and ( 16), (1)   x H and (1)   y H denote the magnetic field components of H x and H y , respectively, at the collocation point in subdomain 1, including the boundary points.Through Eqs. ( 14) and ( 16), the field unknowns at interface boundary points appearing in Eq. ( 12) can be represented by interior collocation points inside subdomains 1 and 2. After considering the physical interface conditions, the final matrix equation is no longer in a block-diagonal form but becomes a matrix with approximately half full elements.
Considering the basis functions, we chose to use Chebyshev polynomials in order to expand the mode fields in interior subdomains because they are mathematically robust with respect to nonperiodic structures, and Laguerre-Gaussian functions (LGFs) with exponential decay characteristics are used to expand the mode fields in semi-infinite exterior subdomains.For Chebyshev polynomials, the explicit form of the Lagrange-type interpolation function is [16] ( ) where T v (x) is the general Chebyshev polynomial of order v and x p denotes the collocation points for the Chebyshev polynomials.For LGFs, the explicit form is [16] where L v (αx) is the Laguerre polynomial of order v and x p denotes the collocation points for LGFs.In Eq. ( 19), the scaling factor α significantly affects the numerical accuracy of the results obtained with LGFs.The definite value of α has been derived in a previous study [17].
The following iterative procedure is used to solve the stationary mode characteristics of nonlinear waveguides: (i) For a given input power P, specify the mode effective index and mode field of the corresponding linear condition as initial values. ( where T denotes the transpose.The actual magnetic field can be expressed as actual 0 2 / ( ) (iii) Once we obtain the new n e and H actual , the electric-field components are calculated by Eqs.8(a)-8(c), and then the updated nonlinear permittivities are obtained through Eq. ( 3).
(iv) Repeat steps (ii) and (iii) until the criterion

Simulation results and discussion
To verify the numerical accuracy and effectiveness of the proposed PSM, we calculate the power dispersion relation of a nonlinear strip dielectric optical waveguide, which has been solved by Obayya et al. using the FEM-based IDBPM [14].Furthermore, for a nonlinear plasmonic waveguide composed of a thin metal film surrounded by a saturable nonlinear substrate and linear cladding, we analyze how the stationary LR-SPP mode depends on the geometrical parameters that vary according to the power level.

A nonlinear strip dielectric optical waveguide
A strip waveguide with width w = 2 µm and thickness t = 1.2 µm is illustrated in Fig. 2(a).The permittivities of the linear core and linear cladding regions are ε co = 1.57 2 and ε cl = 1.55 2 , respectively, at the operating wavelength λ = 0.515 µm.The nonlinear permittivity of the saturable nonlinear substrate is considered as  [14].By the PSM, the geometry structure of the waveguide is divided into 9 subdomains, as shown in Fig. 2(b), and the mode fields in each subdomain are expanded with suitable basis functions.For example, subdomains 1, 3, 7, and 9 are semi-infinite in both x and y directions; thus the mode fields are all expanded with LGFs.As for subdomains 2 and 8, the mode profiles in the x direction are expanded with Chebyshev polynomials and those in the y direction are expanded with LGFs.Considering the fundamental quasi-transverse-electric (quasi-TE) mode 11 y H , the dispersion curve calculated by the present scheme with 15 basis- function terms in each subdomain is shown in Fig. 3 along with that obtained by the FEMbased IDBPM [14].To satisfy the convergent criterion of 10 −6 , the total unknowns needed the present PSM are 2 × (3 × 15 − 2) × (3 × 15 − 2) = 3698.The numerical results calculated by the two schemes agree well and are shown in Fig. 3.This result demonstrates the accuracy and effectiveness of the proposed PSM.For lower input power, the effective index increases slightly as the input power increases.Once the power exceeds the threshold P th = 83 µW, the bistable phenomenon begins, as shown by the sharp jump of the effective index (see Fig. 3).The transition of the mode behavior is due to self-focusing within the nonlinear substrate and leads to power shifting from the core region to the substrate region.To display mode switching, the major (H y ) and minor (H x ) magnetic fields of the stationary mode at a power of P = 80µW on the lower branch of the dispersion curve are shown in Fig. 4(a) and 4(b), respectively.The power of the H y is confined mainly in the linear core region, and the amplitude of the H y is approximately 10 times greater than that of the H x .show the H y and H x , respectively, at P = 80 µW for the upper branch of the dispersion curve.The profile of the H y shifts into the nonlinear substrate region, and the amplitude of the H y is of the order 10 3 greater than that of the H x .The result validates the view that the minor field can be ignored while the mode is operated in the upper branch.The bistability of the nonlinear waveguide is often used to construct all-optical switching devices.

Fundamental long-range surface plasmon polariton mode of a nonlinear strip plasmonic waveguide
The low-loss LR-SPP mode of a thin metal strip bounded by a linear medium has been extensively studied [24][25][26][27][28][29][30], but the configuration in which the nonlinear medium surrounds the metal strip is mostly unexplored.To date, to the best of our knowledge, only a single report exists [10] that investigates the stationary LR-SPPs and that considers Kerr-type nonlinearities.Most reports [5][6][7][8][9] dealing with nonlinear SPPs focus on planar plasmonic waveguide structures.To fill this void, herein, we report the characteristics of the fundamental LR-SPP mode in a plasmonic waveguide composed of a metal strip of width w and thickness t bound between a nonlinear substrate and a linear cladding media.The waveguide structure is shown in Fig. 1(a) except that the core region is a metal strip.When operated at the telecommunication wavelength λ = 1.55 µm, the permittivities of the metal (Au) strip and the linear cladding are ε f = -132 -12.65i and ε cl = 1.75 2 , respectively.While considering the more practical condition in which the nonlinear coefficient of real Kerr-type media is much lower and saturates, the permittivity of the substrate is assumed to be saturable following the form of Eq. ( 21) with ∆ε sat = 0.31 but with a different value ε sub = 1.75 2 .Figures 6(a) and 6(b) show, for width w = 4 µm and thickness t = 20, 50, and 80 nm, the variations of real and imaginary parts, respectively, of the effective index as a function of power using 15 terms of the basis functions in each subdomain.As the input power increases, both the real and imaginary parts of the effective index increase monotonically if the refractive index is not saturated.Figures 7(a)-7(c) show, for t = 80 nm, the mode patterns (|H x | field, the dominate component) for P = 0, 50, and 100 µW, respectively.The monotonic increase in the real and imaginary parts of the effective index with power reveals a better mode confinement and higher loss, respectively.Note that the degree of mode asymmetry increases significantly and the mode maximum is increasingly localized at the metal-cladding interface.The increase of refractive index with power in the nonlinear substrate region, which results from the self-focusing effect, resembles the fundamental LR-SPP modes that are supported in asymmetric structures [25,27].Similar results are also found for an infinitely wide metal film surrounded by the same medium [23,31,32].In contrast, in dielectric waveguides, the highest power is localized in the nonlinear substrate, as shown in the first example.These results validate the view that the nonlinear LR-SPP mode exhibits mode focusing, but not mode switching, as observed in nonlinear dielectric waveguides (at least for a realistic power intensity).Reducing the thickness to t = 20 nm, the variation in the imaginary part of the effective index with respect to power becomes rather small, as shown in Fig. 6(b).The mode patterns at powers P = 0, 50, and 100 µW are shown in Figs.8(a), 8(b), and 8(c), respectively.The degree of mode asymmetry in this case is smaller than that for greater thicknesses.To analyze the impact of a decreased width, we calculate the power dispersion curves for w = 2, 4, and 8 µm and thickness t = 50 nm.The real and imaginary parts of the effective index of refraction are shown in Figs.9(a) and 9(b), respectively, as functions of power.The variations in the real part of the effective index with respect o the input power are analogous for these widths; however, for w = 8 µm, the variations in the imaginary parts exhibit fairly different features.The energy losses for w = 2 and 4 µm increase linearly with the input power, whereas it increases exponentially for w = 8 µm.To gain insight into the results, the mode patterns for w = 2 µm at P = 0, 50, and 100 µW are shown in Figs.10(a)-10(c), respectively.The variation in mode patterns with respect to power also exhibits mode asymmetry and is increasingly localized at the metal-cladding interface, which is similar to the results shown above for w = 4 µm.To understand the sharp variation in energy loss with respect to input power, the mode patterns for w = 8 µm at P = 0, 50, and 100 µW are shown in Figs.11(a)-11(c), respectively.We observe that the mode patterns are different from that for narrower widths.As the input power increases, the mode profiles localized in the substrate region are better confined and the mode intensities localized in the cladding region decrease significantly.Therefore, the net result is an exponential increase in the energy loss.The abovementioned results indicate that no bistable phenomenon occurs in the nonlinear fundamental LR-SPP mode of a twodimensional metal strip waveguide at a reasonable input power intensity of P = 100 µW for the geometry sizes considered in this study.The guided mode found in a nonlinear plasmonic waveguide may be regarded as a hybrid mode composed of a SPP mode and a dielectric waveguide mode.This view is similar to that proposed by Oulton et al. [33] for a hybrid optical waveguide composed of a dielectric nanowire separated from a metal surface by a dielectric gap.The coupling between the LR-SPP mode and the dielectric waveguide mode supported by the sufficiently strong self-focusing effect may lead to mode switching if the input power is increased without any limit.Although mode switching is not observed for P = 100 µW, it is observed in a planar plasmonic structure at a certain thicknesses of the metal film and for reasonable input power [23].

Conclusion
To analyze nonlinear waveguide problems, herein, we develop a full-vectorial nonlinear mode solver based on the pseudospectral method (PSM).The accuracy and efficiency of the proposed PSM are first verified by comparing with results obtained by applying the imaginary-distance beam-propagation method based on the popular finite element method to a nonlinear bistable dielectric optical waveguide.Furthermore, we characterize and discuss the power dispersion relationships of the long-range surface plasmon polariton (LR-SPP) modes of a nonlinear plasmonic waveguide composed of a metal strip surrounded by linear cladding and a saturable nonlinear substrate.The mode maximum is increasingly localized at the metal-cladding interface as the input power increases, and the degree of mode asymmetry increases with the thickness.Moreover, the mode energies of the fundamental LR-SPP mode attenuate linearly with the input power for the narrower metal strips.Beyond a specific width of metal strip, however, the loss of mode energy begins to increase exponentially.In addition, no bistable phenomenon or mode switching is observed at reasonable input powers of P = 100 µW in a two-dimensional nonlinear plasmonic waveguide considered herein unless the input power is increased without a limit.These results provide valuable insights for designing functional photonic devices based on the LR-SPP modes of nonlinear plasmonic waveguides.
δ mp denotes the Kronecker delta.For an arbitrary interior rectangular subdomain, the (n x + 1) × (n y + 1) collocation points are illustrated in Fig.1(a) (where n x = n y = 10 is illustrated).

Fig. 1 .
Fig. 1.(a) Mesh division of an arbitrary interior subdomain and (b) mesh division of a problem with two subdomains (labeled 1 and 2).
where s denotes the s-th iteration.

Fig. 2 .
Fig.2.Schematic of (a) cross section of nonlinear strip waveguide with a nonlinear substrate and (b) division of its computational domain.

Fig. 3 .
Fig. 3. Effective index versus input power for fundamental quasi-TE mode 11 y H .

Fig. 4 .H
Fig. 4. (a) Major magnetic field component (Hy) and (b) minor field component (Hx) of 11 y H mode for an input power of P = 80 µW on the lower branch of the power dispersion curve.

Fig. 5 .
Fig. 5. (a) Major magnetic field component (Hy) and (b) minor field component (Hx) of 11 y H mode for an input power of P = 80 µW on the upper branch of the power dispersion curve.

Fig. 6 .
Fig. 6.Real (a) and imaginary (b) parts of the effective index of the fundamental LR-SPP mode as a function of power for width w = 4 µm and various thicknesses.

Fig. 9 .
Fig. 9. Real (a) and imaginary (b) parts of effective index of fundamental LR-SPP mode as a function of power at the thickness t = 50 nm and different widths.