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-method-based 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
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 , optical switches , wavelength auto-routers , and logic gates  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–10] have been reported to offer additional factors for tailoring the mode characteristics of SPPs. In Ref , 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.  proposed nonlinear nanofocusing in a tapered plasmonic waveguide. In Ref , 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–9] concerned only planar structures, not practical two-dimensional (2D) channel waveguides. Until now, thin metal strips surrounded by Kerr-type nonlinear materials  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–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) , 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–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 . Thus, the present study extends the previous work  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.
2. 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 equationEq. (1) can be represented by two transverse magnetic components Hx and Hy to obtain the following full vector eigenvalue equation:Equations (6) and (7) are used as interfacial conditions according to the continuities of Hz and Ez 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:
3. 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 Hx and Hy field components denoted Hx,pq and Hy,pq at the (nx + 1) (ny + 1) grid (collocation) points:Fig. 1(a) (where nx = ny = 10 is illustrated).
Substituting Eqs. (10a) and (10b) into Eq. (4), we find that Eq. (4) is satisfied exactly at the specific (nx + 1) (ny + 1) collocation points determined by the chosen basis functions. The operators Pxx, Pxy, Pyx, and Pyy in Eq. (4) areEq. (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 Hx and Hy boundary, and the continuities of Hz and Ez in Eqs. (6) and (7), respectively. Accordingly, the number of grid points at each subdomain is reduced to (nx − 1) (ny − 1). For brevity, we formulate the coupling of Hx and Hy using only two subdomains (labeled 1 and 2), as shown in Fig. 1(b). At interface boundary points [i.e., where i = 0, 1, 2, …, ny, 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 Hx and Hy is implicitly satisfied, and the continuity of Hz in Eq. (6) givesEq. (7) givesEqs. (14) and (16), and denote the magnetic field components of Hx and Hy, 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]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 . 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.
- (ii) Solve the final matrix equation to find the new ne and H, and then normalize the numerically calculated field to beby setting
where T denotes the transpose. The actual magnetic field can be expressed as .
- (iv) Repeat steps (ii) and (iii) until the criterion is satisfied, where s denotes the s-th iteration.
4. 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 . 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.
4.1 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.572 and εcl = 1.552, respectively, at the operating wavelength λ = 0.515 μm.
The nonlinear permittivity of the saturable nonlinear substrate is considered as14]. 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, 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 FEM-based IDBPM . To satisfy the convergent criterion of 10−6, the total unknowns needed by the present PSM are 2 (315 − 2) (315 − 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 Pth = 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 (Hy) and minor (Hx) 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 Hy is confined mainly in the linear core region, and the amplitude of the Hy is approximately 10 times greater than that of the Hx.
In addition, Figs. 5(a) and 5(b) show the Hy and Hx, respectively, at P = 80 μW for the upper branch of the dispersion curve. The profile of the Hy shifts into the nonlinear substrate region, and the amplitude of the Hy is of the order 103 greater than that of the Hx. 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.
4.2 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–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  that investigates the stationary LR–SPPs and that considers Kerr-type nonlinearities. Most reports [5–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.752, 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.752. 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 (|Hx| 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 two-dimensional 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.  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 .
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.
The author would like to thank the National Science Council of the Republic of China, Taiwan for financially supporting this research under Contract No. NSC 99-2112-M-005-005-MY3.
References and links
1. T. H. Wood, “Multiple quantum well (MQW) waveguide modulators,” J. Lightwave Technol. 6(6), 743–757 (1988). [CrossRef]
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). [CrossRef]
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). [CrossRef] [PubMed]
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). [CrossRef]
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). [CrossRef] [PubMed]
9. J. R. Salgueiro and Y. S. Kivshar, “Nonlinear plasmonic directional couplers,” Appl. Phys. Lett. 97(8), 081106 (2010). [CrossRef]
10. A. Degiron and D. R. Smith, “Nonlinear long-range plasmonic waveguides,” Phys. Rev. A 82(3), 033812 (2010). [CrossRef]
11. K. Hayata and M. Koshiba, “Full vectorial analysis of nonlinear-optical waveguides,” J. Opt. Soc. Am. B 5(12), 2494–2501 (1988). [CrossRef]
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). [CrossRef]
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). [CrossRef]
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). [CrossRef]
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, 2nd edition, 2001).
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). [CrossRef]
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). [CrossRef]
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). [CrossRef] [PubMed]
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). [CrossRef]
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). [CrossRef]
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). [CrossRef]
28. S. J. Al-Bader, “Optical transmission on metallic wires-fundamental modes,” IEEE J. Quantum Electron. 40(3), 325–329 (2004). [CrossRef]
30. P. Berini, “Long-range surface plasmon polaritons,” Adv. Opt. Photon. 1(3), 484–588 (2009). [CrossRef]
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). [CrossRef]
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). [CrossRef]