Pseudospectral mode solver for analyzing nonlinear optical waveguides

Chia-Chien Huang

doi:10.1364/OE.20.013014

1. Introduction

The use of materials with nonlinear permittivity responses enriches optical guided wave propagation, and it induces different effects than those observed for materials with linear permittivity responses. Many experimental [1,2] and theoretical [3–26] studies have analyzed nonlinear optical waveguides characterized by intensity-dependent permittivities, because such waveguides have unique properties such as low threshold power, self-focusing or defocusing, and bistability. Various all-optical signal-processing devices such as waveguide modulators [27], optical switches [28], wavelength auto-routers [29], and logic gates [30] have been implemented based on these nonlinear mechanisms. Recently, combining metallic materials with nonlinear dielectric media [31–37] has been reported to offer increased tailoring of mode characteristics of surface plasmon polaritons (SPPs) in plasmonic waveguides. Davoyan et al. presented dispersion diagrams of nonlinear-guided modes of a nonlinear plasmonic slot waveguide [31] and analyzed phase matching in metal-dielectric nonlinear structures [32]. Rukhlenko et al. [33,34] derived a dispersion relation for SPPs in a nonlinear plasmonic waveguide with arbitrary power-law nonlinearity, and they predicted the existence of backward-propagating modes. The power diagram of nonlinear plasmonic directional couplers [35] is reported to show substantial differences from that of nonlinear dielectric couplers. Nonlinear nanofocusing in tapered plasmonic waveguides [36] is also proposed to enhance nonlinear effects. In [37], a plasmonic waveguide made of a thin strip surrounded by Kerr-type nonlinear materials was reported to analyze the mode properties of long-range (LR) and short-range (SR) SPPs.

The research mentioned above suggests that developing an efficient and accurate computational scheme is essential for understanding the propagation characteristics of nonlinear waveguide components and integrated optical circuits. Unfortunately, only the transverse electric (TE) wave in a slab waveguide has an exact analytical solution [3,5,6,8,11]. No analytical solution exists for the transverse magnetic (TM) wave in a slab waveguide with a biaxial nonlinear refractive index. Until now, the reported analytical solutions [4,9,10], have been limited to the uniaxial approximation, which considers the intensity-dependent permittivities resulting from only the longitudinal component of the electric field. However, the transverse component of the electric field is often stronger than the longitudinal one; therefore, the application of this approximation is limited in practical devices [7]. In addition, no accurate analytical solutions exist for practical 2D channel nonlinear waveguides. Therefore, numerical methods have been proposed to analyze nonlinear optical waveguides [12–26]. Except for these studies [13,17,18,20,24] that solve 2D waveguide structures by using the finite element method (FEM), most numerical schemes only solve one-dimensional slab waveguides. Recently, a pseudospectral method [38] successfully demonstrated fast convergence and high-order accuracy for solving the propagation characteristics of both guided and leaky modes of various linear optical and plasmonic waveguides [39–44]. A pseudospectral scheme requires remarkably less memory storage and achieves the same order of accuracy [39,43,44], as that of the finite difference method (FDM) and FEM. The superior computational performance of the pseudospectral scheme led to the aim of this study: to develop a new mode solver based on the pseudospectral scheme to study the mode characteristics of one- and two-dimensional nonlinear optical and plasmonic waveguides.

In most numerical mode solvers [12,13,16,17,19,24,26], a nonlinear eigenvalue problem is solved iteratively by assuming the effective index and mode field of the corresponding linear mode as the initial eigenvalue and eigenvector, respectively, to achieve the final convergent solution. This approach converges to stable modes, but it fails to find unstable solutions. The effective index may be viewed mathematically as a multi-valued function of the power for certain conditions; therefore, an alternative approach [18] that uses the power as an eigenvalue has been proposed. This new approach can obtain a complete dispersion relation that involves both the stable and unstable stationary solutions of nonlinear waveguides. In this study, two approaches to the pseudospectral scheme are implemented.

The remainder of this study is organized as follows. Section 2 presents the mathematical formulations of wave equations for planar and 2D waveguides. The computational schemes including the pseudospectral scheme and two iterative approaches are derived in Section 3. Section 4 demonstrates the accuracy of the proposed scheme by comparing its solution with the exact analytical solution and solutions from other numerical schemes; it also 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 a time dependence of exp(jωt) propagating along the z direction obeys the following wave equation:

\nabla \times ({[ε]}^{- 1} \nabla \times H) - ω^{2} μ_{0} H = 0,

where ω is the angular frequency, μ₀ is the permeability in vacuum, and [ε] is the nonlinear permittivity tensor. For practical purposes, this study considers the complete biaxial nature of materials; thus, the nonlinear permittivity tensor can be expressed as follows:

[ε] = ε_{0} [ε_{r}] = ε_{0} [\begin{matrix} {\tilde{ε}}_{x} & 0 & 0 \\ 0 & {\tilde{ε}}_{y} & 0 \\ 0 & 0 & {\tilde{ε}}_{z} \end{matrix}],

where ε₀ is the permittivity in vacuum, and [ε_r] is the relative permittivity tensor. Considering the stationary guided modes of nonlinear media, the magnetic field components can be assumed to be of the form exp(−jßz), where β = k₀n_e is the propagation constant along the z direction,

k_{0} = ω \sqrt{μ_{0} ε_{0}}

is the wave number in vacuum, and n_e is the effective refractive index of guided modes.

For slab waveguides with uniform structure in the x direction, $\partial / \partial x = 0$ , the wave equation of the TE polarization is obtained as follows:

\frac{d^{2} H_{y}}{d y^{2}} + k_{0}^{2} ({\tilde{ε}}_{x} - n_{e}^{2}) H_{y} = 0,

where the nonlinear relative permittivity is given by

{\tilde{ε}}_{x} = ε_{x} + a f (E_{x}),

a = c_{0} ε_{0} ε_{x} \bar{n},

where c₀ is the light velocity in vacuum, ε_x is the linear relative permittivity,

\bar{n}

is the nonlinearity coefficient, and f(E_x) is a function of the electric field component in the x direction. For example, f(E_x) = |E_x|² for materials that have Kerr-type nonlinearity. The input power per unit length along the x direction for TE polarization is given by

P = \frac{1}{2} \int_{- \infty}^{\infty} E_{x} H_{y}^{*} d y = \frac{Z_{0} β}{2 k_{0}} \int_{- \infty}^{\infty} \frac{1}{{\tilde{ε}}_{x}} {| H_{y} |}^{2} d y,

where * denotes the complex conjugate and Z₀ denotes the intrinsic impedance in vacuum (Z₀ = 377 Ω). The actual electric field component E_x can be obtained from Maxwell’s equations as follows:

E_{x} = \frac{Z_{0} β}{k_{0} {\tilde{ε}}_{x}} H_{y} .

For TM polarization, the wave equation is shown below:

\frac{d}{d y} (\frac{1}{{\tilde{ε}}_{z}} \frac{d H_{x}}{d y}) + (k_{0}^{2} - \frac{β^{2}}{{\tilde{ε}}_{y}}) H_{x} = 0,

and the nonlinear relative permittivities are given by

{\tilde{ε}}_{y} = ε_{y} + a g (E_{y}) + b h (E_{z}),

{\tilde{ε}}_{z} = ε_{z} + b g (E_{y}) + a h (E_{z}),

where ε_y and ε_z denote the linear relative permittivities along the y and z directions, respectively, and g(E_y) and h(E_z) denote functions of the electric field components E_y and E_z, respectively. The value of b determines the particular nonlinear mechanism. For example, b = a represents electrostrictive nonlinearity, and b = a/3 represents electronic nonlinearity [12]. The input power per unit length along the x direction for TM polarization is given by

P = - \frac{1}{2} \int_{- \infty}^{\infty} E_{y} H_{x}^{*} d y = \frac{Z_{0} β}{2 k_{0}} \int_{- \infty}^{\infty} \frac{1}{{\tilde{ε}}_{y}} {| H_{x} |}^{2} d y .

The actual electric field components E_y and E_z can be obtained from Maxwell’s equations as follows:

E_{y} = - \frac{Z_{0} β}{k_{0} {\tilde{ε}}_{y}} H_{x},

E_{z} = \frac{j Z_{0}}{{\tilde{ε}}_{z}} \frac{d H_{x}}{d y} .

For 2D structures, only the weakly guiding approximation is considered in this study. By simplifying Eq. (1) and considering isotropic media with ${\tilde{ε}}_{x} = {\tilde{ε}}_{y} = {\tilde{ε}}_{z} = \tilde{ε},$ the scalar wave equation is obtained as follows:

\frac{\partial^{2} ξ}{\partial x^{2}} + \frac{\partial^{2} ξ}{\partial y^{2}} + k_{0}^{2} (\tilde{ε} - n_{e}^{2}) ξ = 0,

where

ξ (x, y)

is regarded as E_x here. Similarly, the nonlinear relative permittivity is given by

\tilde{ε} = ε + c f (ξ),

and the input power is obtained as

P = \frac{β}{2 Z_{0} k_{0}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} {| ξ (x, y) |}^{2} d x d y .

3. Computational schemes

The idea of the pseudospectral scheme [38] is to divide the computational window, which generally consists of several different materials in waveguide problems, into several subdomains with homogeneous or continuous refractive index profiles. In each subdomain, the unknown optical field ψ(y) is expanded by a set of appropriate interpolation functions θ_j (y) and unknown grid point values ψ_j as follows:

ψ (y) = \sum_{j = 0}^{n} θ_{j} (y) ψ_{j},

where

θ_{j} (y) = \frac{ρ_{n + 1} (y)}{ρ_{n + 1}^{'} (y) (y - y_{j})}, 0 \leq i \leq n .

In Eq. (15), ρ_{n + 1}(y) denotes a particular basis function of order n + 1, the prime denotes the first derivative of ρ_{n + 1}(y) with respect to y, and y_j denotes the corresponding collocation point fulfilling the condition of θ_j(y_i) = δ_ij, where δ_ij denotes the Kronecker delta. The explicit form of θ_j(y) is written in terms of Chebyshev polynomials as follows [38]:

θ_{j} (y) = \frac{{(- 1)}^{j + 1} (1 - y^{2}) T_{n}^{'} (y)}{c_{j} n^{2} (y - y_{j})}, y \neq y_{j} .

where T_n(y) denotes the Chebyshev polynomial of order n, c₀ = c_n = 2, and c_j = 1 (1 ≦ j ≦ n − 1). For LG functions, namely ψ_{n + 1}(αy) = (αy)exp(−αy/2)L_n(αy), where L_n(αy) denotes the Laguerre polynomial of order n, the explicit form of θ_j(y) is given as follows [38]:

θ_{j} (α y) = \frac{e^{- α y / 2}}{e^{- α y_{j} / 2}} \frac{(α y) L_{n} (α y)}{{(α y L_{n})}^{'} (α y_{j}) (α y - α y_{j})}, y \neq y_{j} .

Here, the parameter α in Eq. (17) is called the scaling factor, and it influences the accuracy of a given number of terms of the basis functions. The definite determination of α has been derived in our previous work [39].

Once the basis functions are determined, Eq. (14) is substituted into Eq. (3) or Eq. (7), which must be satisfied at n + 1 collocation points; therefore,

D {\bar{ψ}} = {(k_{0} n_{e})}^{2} {\bar{ψ}},

where D denotes an (n + 1) × (n + 1) matrix, and

{\bar{ψ}}

denotes the normalized grid point vector (eigenvector) that satisfies the relation

\int_{- \infty}^{\infty} {| \bar{ψ} (y) |}^{2} d y = 1,

where

ψ = γ \bar{ψ}

. For the TE modes, we have

γ_{T E} = \sqrt{(2 P {\tilde{ε}}_{x}^{}) / (Z_{0} n_{e})}

and

\bar{ψ} = {\bar{H}}_{y}

, and the elements of matrix D are

D_{i j} = θ_{j}^{(2)} (y_{i}) + k_{0}^{2} {\tilde{ε}}_{x} (y_{i}) δ_{i j},

where

θ_{j}^{(k)} (y)

denotes the k-th order derivative of a basis function of degree j with respect to y. For the TM modes, the relations are

γ_{T M} = \sqrt{(2 P {\tilde{ε}}_{y}^{}) / (Z_{0} n_{e})}

,

\bar{ψ} = {\bar{H}}_{x}

, and

D_{i j} = \frac{{\tilde{ε}}_{y} (y_{i})}{{\tilde{ε}}_{z} (y_{i})} θ_{j}^{(2)} (y_{i}) - \frac{{\tilde{ε}}_{y} (y_{i})}{{\tilde{ε}}_{z}^{2} (y_{i})} (\frac{d {\tilde{ε}}_{z} (y)}{d y} |_{y = y_{i}}) θ_{j}^{(1)} (y_{i}) + k_{0}^{2} {\tilde{ε}}_{y} (y_{i}) δ_{i j} .

Assembling the contributions from subdomains 1 to t, the global matrix eigenvalue equation is formed as follows:

[\begin{matrix} D_{1} & 0 & 0 & 0 \\ 0 & D_{2} & 0 & 0 \\ 0 & 0 & ⋱ & 0 \\ 0 & 0 & 0 & D_{t} \end{matrix}] [\begin{matrix} {\bar{ψ}}_{1} \\ {\bar{ψ}}_{2} \\ ⋮ \\ {\bar{ψ}}_{t} \end{matrix}] = {(k_{0} n_{e})}^{2} [\begin{matrix} {\bar{ψ}}_{1} \\ {\bar{ψ}}_{2} \\ ⋮ \\ {\bar{ψ}}_{t} \end{matrix}] .

In Eq. (22), the rows including all the interface points between different materials must be replaced by the following interface conditions for TE and TM modes, respectively:

{\bar{ψ}}_{j} (y_{r}^{+}) = {\bar{ψ}}_{j} (y_{r}^{-}), {\bar{ψ}}_{j}^{(1)} (y_{r}^{+}) = {\bar{ψ}}_{j}^{(1)} (y_{r}^{-}),

and

{\bar{ψ}}_{j} (y_{r}^{+}) = {\bar{ψ}}_{j} (y_{r}^{-}), {\tilde{ε}}_{z} (y_{r}^{-}) {\bar{ψ}}_{j}^{(1)} (y_{r}^{+}) = {\tilde{ε}}_{z} (y_{r}^{+}) {\bar{ψ}}_{j}^{(1)} (y_{r}^{-}),

where

{\bar{ψ}}_{j}^{(1)} (y)

denotes the first-order derivative of

\bar{ψ} (y)

of degree j with respect to y, and

y_{r}^{+}

and

y_{r}^{-}

denote positions infinitesimally close to the interfaces y_r (r = 1, 2, 3, …, t − 1) from the top and bottom sides, respectively.

Mathematically, Eq. (22) with Eq. (23) or Eq. (24) leads to a nonlinear eigenvalue problem with eigenvalue n_e and eigenvector H_x or H_y, and the iterative steps for solving the TE modes are described as follows:

(i) For a given input power P, specify the mode effective index and mode field of the corresponding linear condition (i. e., $\bar{n} = 0$ ) as the initial value $n_{e}^{0}$ and initial mode field ${{\bar{H}}_{y}}^{0}$ , respectively, to compute the first relative permittivity ${\tilde{ε}}_{x}^{1}$ via the initial actual value of ${E_{x}}^{0}$ in Eq. (4a).
(ii) The value ${E_{x}}^{0}$ is obtained from the initial actual value ${H_{y}}^{0}$ and ${\tilde{ε}}_{x}^{} = {\tilde{ε}}_{x}^{0} = ε_{x}^{}$ via Eq. (6), and then the value ${H_{y}}^{0}$ is calculated from the relation ${H_{y}}^{0} = γ^{0} {{\bar{H}}_{y}}^{0}$ through Eq. (5), where $γ^{0} = \sqrt{(2 P {\tilde{ε}}_{x}^{0}) / (Z_{0} n_{e}^{0})}$ .
(iii) Substitute ${\tilde{ε}}_{x}^{1}$ into Eq. (3) to find a new eigenvalue $n_{e}^{1}$ and eigenvector ${{\bar{H}}_{y}}^{1}$ (the superscripts denote the first iterative results of n_e and ${{\bar{H}}_{y}}$ ).
(iv) Repeat steps (i)–(iii) until the criterion $(n_{e}^{s + 1} - n_{e}^{s}) < 10^{- 6}$ is satisfied.

However, the conventional iterative approach (called the n_e-based approach in the following section), which is used properly in tackling linear waveguide modes, fails to find the unstable region with negative slope in the n_e versus P curve. This is because the effective index may be a multi-valued function of the power under certain conditions. Therefore, an alternative approach [18] that uses the power as an eigenvalue has been proposed to obtain both the stable and unstable stationary solutions of nonlinear waveguides. In this study, this alternative approach is also combined into the pseudospectral scheme to determine a complete n_e – P curve. This alternative iterative procedure (called the P-based approach in the following section) is described as follows:

(i) For a given n_e, specify the initial input power value $P^{0} \to 0$ , ${\tilde{ε}}_{x}^{0} = ε_{x}^{}$ , and ${{\bar{H}}_{y}}^{0}$ from the computed mode field of the corresponding linear condition.
(ii) Modify Eq. (3) to $\frac{d^{2} H_{y}}{d y^{2}} + k_{0}^{2} (ε_{x} - n_{e}^{2}) H_{y} = - k_{0}^{2} a f (E_{x}) H_{y},$

and it becomes a nonlinear generalized eigenvalue problem as follows:

\tilde{D} {{\bar{H}}_{y}} = P \tilde{B} {{\bar{H}}_{y}} .

The elements of the matrices $\tilde{D}$ and $\tilde{B}$ are

{\tilde{D}}_{i j} = θ_{j}^{(2)} (y_{i}) + k_{0}^{2} (ε_{x} (y_{i}) - n_{e}^{2}) δ_{i j},

{\tilde{B}}_{i j} = \frac{- k_{0}^{2} a}{P} f (E_{x} (y_{i})) δ_{i j} .

The value ${E_{x}}^{0}$ is obtained from ${H_{y}}^{0}$ and ${\tilde{ε}}_{x}^{0} = ε_{x}^{}$ via Eq. (6), in which the value ${H_{y}}^{0}$ here is calculated from the relation ${H_{y}}^{0} = η^{0} {{\bar{H}}_{y}}^{0}$ through Eq. (5), where $η^{0} = \sqrt{(2 P^{0} {\tilde{ε}}_{x}^{0}) / (Z_{0} n_{e})}$ .

(iii) Solve Eq. (26) to obtain the new iterative values P¹ and ${{\bar{H}}_{y}}^{1}$ , and then calculate the value ${\tilde{ε}}_{x}^{1}$ .
(iv) Repeat steps (i)–(iii) until the desired convergent value of P of the (s + 1)-th iterative step is achieved (satisfy the criterion $(P^{s + 1} - P^{s}) < 10^{- 4}$ ).

The computational procedure for solving the TM mode is similar to that for solving the TE mode except for the following two parts: 1) Eq. (7) relies on the relative permittivities of ${\tilde{ε}}_{y}^{}$ and ${\tilde{ε}}_{z}^{}$ , and they are affected by both the electric field components E_y and E_z because of the biaxial nature as formulated in Eqs. (8)a) and (8b); 2) the right interface condition in Eq. (24) involves ${\tilde{ε}}_{z}^{}$ and needs to be updated together with Eq. (7).

For channel waveguides with 2D cross sections, the optical field ξ(x, y) is expanded by a product of two separable basis sets φ_i(x) and θ_j(y) and the corresponding grid point values ψ_i,j as follows:

ξ (x, y) = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} φ_{i} (x) θ_{j} (y) ξ_{i, j},

where n_x and n_y denote the number of terms of basis functions in the x and y directions, respectively. Note that different sorts of basis sets φ_i(x) and θ_j(y) can be chosen to expand the optical field along distinct directions in each subdomain mainly according to their field features. For subdomains with finite extent, the unknown fields are expanded using Chebyshev polynomials, because they are robust and can represent the non-periodic region. For subdomains with semi-infinite extent, the unknown fields are expanded using LG functions because of the good match between the exponential decay fields of guided modes and mathematical features of LG functions. Substituting Eq. (29) into Eq. (11), we have

A {\bar{ξ}} = {(k_{0} n_{e})}^{2} {\bar{ξ}},

where

{\bar{ξ}}

denotes the normalized grid point vector (eigenvector) that satisfies the relation

\int_{- \infty}^{\infty} {| \bar{ξ} (x, y) |}^{2} d x d y = 1.

where

ξ = ν \bar{ξ}

and

ν = \sqrt{(2 P Z) / (n_{e})}

. The elements of matrix A are

\begin{array}{l} A = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} {[\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + k_{0}^{2} \tilde{ε}]}_{x}_{= x_{i}, y = y_{j}} \\ = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {φ_{p}^{(2)} (x) θ_{q} (y) + φ_{p} (x) θ_{q}^{(2)} (y) + k_{0}^{2} \tilde{ε} (x, y) φ_{p} (x) θ_{q} (y)}] |_{x = x_{i}, y = y_{j}}, \end{array}

where

φ_{p}^{(h)} (x)

and

θ_{q}^{(h)} (y)

denote, respectively, the h-th order derivatives of basis functions φ_p(x) and θ_q(y). Similarly, the contributions from subdomain 1 to t are added to form the global matrix eigenvalue equation as follows:

[\begin{matrix} A_{1} & 0 & 0 & 0 \\ 0 & A_{2} & 0 & 0 \\ 0 & 0 & ⋱ & 0 \\ 0 & 0 & 0 & A_{t} \end{matrix}] [\begin{matrix} {\bar{ξ}}_{1} \\ {\bar{ξ}}_{2} \\ ⋮ \\ {\bar{ξ}}_{t} \end{matrix}] = {(k_{0} n_{e})}^{2} [\begin{matrix} {\bar{ξ}}_{1} \\ {\bar{ξ}}_{2} \\ ⋮ \\ {\bar{ξ}}_{t} \end{matrix}] .

In the scalar wave approximation, the horizontal interface conditions are

\bar{ξ} (y_{r}^{+}) = \bar{ξ} (y_{r}^{-}), \partial \bar{ξ} (y_{r}^{+}) / \partial y = \partial \bar{ξ} (y_{r}^{-}) / \partial y,

where

y_{r}^{+}

and

y_{r}^{-}

denote positions infinitesimally close to interfaces y_r (r = 1, 2, 3, … t − 1) from the top and bottom sides, respectively. In addition, the vertical interface conditions are

\bar{ξ} (x_{r}^{+}) = \bar{ξ} (x_{r}^{-}), \partial \bar{ξ} (x_{r}^{+}) / \partial x = \partial \bar{ξ} (x_{r}^{-}) / \partial x,

where

x_{r}^{+}

and

x_{r}^{-}

denote positions infinitesimally close to interfaces x_r (r = 1, 2, 3, … t − 1) from the right and left sides, respectively. Similarly to the case of slab waveguides, Eq. (33) must be modified by adding the interface conditions in Eqs. (34) and (35) to form the final eigenvalue problem. Moreover, Eq. (11) becomes

\frac{\partial^{2} ξ}{\partial x^{2}} + \frac{\partial^{2} ξ}{\partial y^{2}} + k_{0}^{2} (ε - n_{e}^{2}) ξ = - k_{0}^{2} c f (ξ) ξ,

which becomes a nonlinear generalized eigenvalue problem as follows:

\tilde{A} {\bar{ξ}} = P \tilde{Q} {\bar{ξ}},

where the elements of matrices

\tilde{A}

and

\tilde{Q}

are

\begin{array}{l} \tilde{A} = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} {[\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + k_{0}^{2} ε]}_{x}_{= x_{i}, y = y_{j}} \\ = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {φ_{p}^{(2)} (x) ϕ_{q} (y) + φ_{p} (x) ϕ_{q}^{(2)} (y) + k_{0}^{2} ε (x, y) φ_{p} (x) ϕ_{q} (y)}] |_{x = x_{i}, y = y_{j}} \end{array}

and

\tilde{Q} = \sum_{i = 0}^{n_{x}} \sum_{j = 0}^{n_{y}} [\frac{- k_{0}^{2} c}{P} f (\bar{ξ} (x_{i}, y_{j}))] δ_{i j} .

The iterative procedures are similar to those presented above for slab waveguides.

4. Simulation results and discussion

To verify the numerical accuracy of the proposed scheme, we calculate the dispersion relations for stationary TE mode guided by a slab dielectric nonlinear waveguide. Next the LR- and SR-SPP mode features of a metal film bounded by a nonlinear cladding and a linear substrate are analyzed in detail. Finally, the proposed scheme is further extended to a two-dimensional dielectric nonlinear strip waveguide.

4.1 A slab dielectric waveguide with a Kerr-like nonlinear cladding

A three-layer dielectric nonlinear waveguide is illustrated in Fig. 1 . The Corning 7059 waveguide film composed of a linear medium with permittivity ε_f = 1.57² is bounded by a linear substrate with permittivity ε_s = 1.55² and a Kerr-like nonlinear liquid crystal MBBA cladding with permittivity ${\tilde{ε}}_{c} = ε_{c} + a {| E_{x} |}^{2}$ , where ε_c = 1.55² and $\bar{n} = 10^{- 9}$ m²/W. The operating wavelength of the argon ion laser is λ = 0.515 μm, and the film thickness is d = 2 μm. A benchmark example is often used to examine the accuracy of newly developed TE mode solvers [14–16,21] of nonlinear waveguides because of the availability of an exact solution [8].

Fig. 1 Schematic diagram of a three-layer dielectric waveguide with a Kerr-like nonlinear cladding.

Download Full Size | PDF

In the proposed scheme, the computational domain is divided into three subdomains with homogeneous or continuous refractive index profiles, and the mode field in each subdomain is then expanded by a set of suitable basis functions. The most efficient choice is that the field in the interior region with finite extent is expanded using Chebyshev polynomials and the fields in the exterior regions with semi-infinite extents are expanded using LG functions. Some input powers are given to solve the effective indices for examining the accuracy and efficiency of the proposed scheme. We first analyze the convergences of n_e of the fundamental TE mode versus the iteration time by comparing the results obtained from the proposed scheme using 20 terms of the basis functions for each subdomain with the exact solutions [8], as shown in Fig. 2(a) .

Fig. 2 Convergence of effective index versus (a) iteration time and (b) number of terms of the basis function.

Download Full Size | PDF

As the input power rises, the iteration time to achieve stable convergence increases, because the guided mode profile at lower power is closer to the given initial linear mode characteristic. Note that several iteration times (< 10 times) are enough to achieve stable convergence with an accuracy on the order of 10⁻⁶, even for power as high as P = 100 mW/mm. In addition, the convergence versus the number of terms of the basis function is shown in Fig. 2(b). By increasing the number of terms of the basis function, the accuracy can be further improved to the required tolerance. These results indicate that the proposed scheme for solving a nonlinear waveguide problem is computationally efficient and capable of achieving high accuracy. In the following examples, only 20 terms of the basis functions of each subdomain are used to calculate the dispersion curves. Figure 3 shows the variation of the effective index versus the power for both the n_e-based (lower blue and upper green lines) and P-based (red line) approaches.

Fig. 3 Effective index versus input power for both the n_e- and P-based approaches.

Download Full Size | PDF

For the n_e-based approach, the bistable phenomenon is clearly observed. As the input power increases, the effective index of the guided mode changes slightly (lower branch); however, when the input power exceeds the threshold value of P_th = 111.8 mW/mm, the film-guided mode confined mostly in the film region switches abruptly to a surface-guided mode confined in the nonlinear cladding because of the sufficiently strong self-focusing effect of the cladding nonlinearity. The blue arrow indicates the switching direction. For the upper branch, the self-focusing effect of the cladding region becomes weaker as the mode power decreases. The refractive index in the cladding region is smaller than that in the film region when the power is lower than a threshold value of P = 35.2 mW/mm; thus, the surface-guided mode switches abruptly back to a film-guided mode, as shown in Fig. 3. The green arrow indicates the corresponding switching direction. Consequently, the bistability of nonlinear waveguides can be used to design an all-optical switch or logical gate if we assign different initial input powers. The inset in Fig. 3 magnifies the unstable region (red line with negative slope). However, this part cannot be obtained from the n_e-based approach, because the mode profiles become unstable as the iterations progress; in addition, the effective index is a multivalued function of the guided mode power. Therefore, to obtain the complete power dispersion curve, we need to adopt the P-based approach to determine both the stable (lower and upper branches) and unstable regions. In Fig. 3, the effective indices of guided modes in stable regions obtained from the two approaches are identical, which suggests that the two approaches for solving stable guided mode characteristics are equivalent. To observe the two different modes, the mode profiles for input powers P = 111 mW/mm (labeled as A) and P = 112 mW/mm (labeled as B) are shown in Fig. 4 ; the dashed lines denote the film boundaries. Note that the mode field of P = 111 mW/mm, which is smaller than P_th, is confined most in the film region, similar to its linear counterpart. Beyond the threshold power, the mode field of P = 112 mW/mm shifts to the nonlinear cladding region, and it indeed resembles a surface mode guided by the interface between the nonlinear cladding and film regions.

Fig. 4 Mode profiles for the input powers P = 111 mW/mm < P_th (blue curve labeled A) and P = 112 mW/mm > P_th (green curve labeled B).

Download Full Size | PDF

To demonstrate the accuracy of the proposed scheme, the complete power dispersion curve for the fundamental TE mode obtained from the proposed scheme with the P-based iteration approach and that obtained from the exact solution [8] are shown in Fig. 5 . Note the excellent agreement between the results from the proposed scheme and those from the exact solution.

Fig. 5 Guided power versus effective index for the fundamental TE mode of a three-layered dielectric nonlinear waveguide.

Download Full Size | PDF

4.2 A slab nonlinear plasmonic waveguide

Low-loss plasmonic components are often created in a thin metal film bounded by linear dielectric media, called an insulator–metal–insulator (IMI) structure, because this structure can support a LR-SPP mode characterized by weak field confinement. The surrounding materials may also be composed of nonlinear media, but only a few studies [37,45,46] have discussed the propagation characteristics of a power-dependent configuration. In [45], [46] the power dispersion relation is derived under the uniaxial approximation, because only TM polarizations are supported by a metal-dielectric interface. In [37] commercial finite-element software is used to investigate stationary LR-SPP modes. For similar configurations, called metal–insulator–metal (MIM) structures, approximate relations between the angular frequency and propagation constant for a given power are derived in [31–35]. Here we study the mode characteristics of an IMI nonlinear plasmonic waveguide that is composed of a thin metal (Au) film with permittivity ε_f = –132 bounded by a linear substrate with permittivity ε_s = 1.75² and a Kerr-like nonlinear cladding with permittivity ${\tilde{ε}}_{r} = ε_{c} +$ a(|E_y|² + |E_z|²), where r = y, z, ε_c = 1.75², and $\bar{n} = 10^{- 9}$ m²/W. (For simplicity, we focus mainly on the nonlinear response of the real part of the propagation constant versus the guided power.) The operating wavelength is at a telecommunication wavelength λ = 1.55 μm, and the film thickness is assumed to be d = 50 nm. The power dispersion curves of the LR-SPP and SR-SPP obtained from the proposed scheme using 20 terms of the basis functions each subdomain are shown in Fig. 7 below. Compared with the results in [37] that use the theory developed in [45,46], our results show good agreement. In Fig. 6 , note that the cut-off power of the LR-SPP mode is larger than P = 981.2 mW/mm (corresponding to n_e = 1.77040); however, the curve of the SR-SPP mode shows that the effective index increases monotonically with the guided power. In this case, no bistable phenomenon can be observed in either mode.

Fig. 6 Guided power versus effective index for the long-range and short-range modes of a thin metal film surrounded by a Kerr-like nonlinear cladding and a linear substrate.

Download Full Size | PDF

To observe the evolution of the mode profiles with the guided power, the relative magnetic field profiles of the LR-SPP modes at several powers (P = 0 (n_e = 1.76415), P = 71.8 (n_e = 1.76805), and P = 981.2 (n_e = 1.77040) mW/mm) are shown in Fig. 7 . As the guided power increases, the maximum of the LR-SPP mode profile occurs at the interface between the linear substrate and the metal film, and most of the power is focused in substrate region. The increase of refractive index in the cladding region by the nonlinear response resembles that the guided modes are supported in asymmetric structures [47]. This shifted phenomenon of mode profile is opposite to that observed in the pure dielectric nonlinear waveguide discussed in Section 4.1.

Fig. 7 Relative magnetic profiles of the LR-SPP mode for the guided powers: (a) P = 0, (b) P = 71.8, and (c) P = 981.2 mW/mm.

Download Full Size | PDF

In addition, the relative magnetic field profiles of the SR-SPP modes at several powers (P = 0 (n_e = 1.78014), P = 19.7 (n_e = 1.80025), and P = 30.9 (n_e = 1.81905) mW/mm) are shown in Fig. 8 . Note that the SR-SPP mode profiles exhibit a tendency opposite to that of the LR-SPP mode profiles: the maximum of the SR-SPP mode profile shifts to the interface between the nonlinear cladding and the metal film as the guided power increases.

Fig. 8 Relative magnetic profiles of the SR-SPP mode for the input powers: (a) P = 0, (b) P = 19.7, and (c) P = 30.9 mW/mm.

Download Full Size | PDF

For analyzing the effect of varying the thickness on the power level, the power dispersion relations of the LR-SPP mode at the thicknesses, d = 20nm, 30nm, 50nm, and 80nm, are shown in Fig. 9 . The power scales of the thicknesses, d = 20nm, 30nm, and 50nm, are magnified by 10 times in order to clearly observe the power dispersion curves. We can see that the cut-off power of the LR-SPP mode decreases significantly as the thickness reduces, and the calculated values of d = 30nm, 50nm and 80nm are P = 194.3, P = 981.2 and P = 15092.5 mW/mm, respectively. In particular, the LR-SPP mode for the thickness d = 20nm exhibits the bistability as that occurred in pure dielectric materials while the input power exceeds P = 107.7 mW/mm. However, no bistable phenomena can be observed for the thicker metal films.

Fig. 9 Guided power versus effective index for the LR-SPP modes of a thin metal film with different thicknesses surrounded by a Kerr-like nonlinear cladding and a linear substrate. The power scales of the thicknesses, d = 20nm, 30nm, and 50nm, have been magnified by x10.

Download Full Size | PDF

To observe the evolution of the mode profiles, the relative magnetic field profiles of the LR-SPP modes at power P = 0 mW/mm for different thicknesses are shown in Fig. 10(a) and that at the maximum powers of different thicknesses (P = 110 (d = 20nm), P = 194 (d = 30nm), P = 981 (d = 50nm) and P = 15092 (d = 80nm) mW/mm) are shown in Fig. 10(b). For the condition d = 20nm, the mode profile penetrates into the nonlinear cladding region by self-focusing while the input power exceeds the threshold power P = 107mW/mm as shown in Fig. 10(b), and the peak of mode profile is not close to the interface of the metal film and the nonlinear cladding as that observed in Fig. 4. In contrast, the penetration of the mode profile into cladding region reduces as the thickness of metal film increases, and most power is focused to the interface between the metal and linear substrate (i. e., no bistability phenomenon is seen as the thickness of metal film exceeds a certain value). The results show that the metal thickness is a critical factor for accomplishing the all-optical switching and bistability devices based on the nonlinear IMI structure.

Fig. 10 Relative magnetic profiles |H_x| of the LR-SPP mode at different thicknesses for the input powers: (a) P = 0 mW/mm for all thicknesses and (b) P = 110 (d = 20nm), P = 194.3 (d = 30nm), P = 981.2 (d = 50nm), and P = 15092.5 (d = 80nm) mW/mm.

Download Full Size | PDF

4.3 A nonlinear strip waveguide with a nonlinear substrate

Now we consider an interesting practical application: we extend the proposed scheme to solve a 2D strip waveguide with a nonlinear substrate. Figures 11(a) and (b) show schematic diagrams of the nonlinear strip waveguide and the division of the computational domain, respectively.

Fig. 11 Schematics showing (a) the cross section of a nonlinear strip waveguide with a nonlinear substrate and (b) the division of its computational domain.

Download Full Size | PDF

The geometric parameters of the core region are w = 2 μm and h = 1.2 μm, and the refractive indices of the core and cladding regions are n_co = 1.57 and n_cl = 1.55, respectively, at the operating wavelength λ = 0.515 μm. The nonlinear MBBA liquid crystal substrate with Gaussian saturation is characterized by

\tilde{ε} = n_{s}^{2} + 2 n_{s} Δ n_{s a t} [1 - \exp (- \frac{a {| ψ |}^{2}}{2 n_{s} Δ n_{s a t}})],

where n_s = 1.55,

a = c_{0} ε_{0} n_{s}^{2} \bar{n}

,

\bar{n} = 10^{- 9}

m²/W, and

Δ

n_sat = 0.1. The computational window is divided by the proposed scheme into nine subdomains as illustrated in Fig. 11(b). The basis functions used depend on whether the extent is finite or semi-infinite. For example, for subdomains 1, 3, 7, and 9, the subdomains possess semi-infinite extents in both the x and y directions; therefore, they are all expanded with LGFs in both the x and y directions. For subdomains 2 and 8, Chebyshev polynomials are used to expand the field profiles in the x direction, and LGFs are used to expand them in the y direction. Note that guided mode profiles are well matched to the mathematical characteristics of LGFs; therefore, the proposed scheme does not require extra effort to determine the computational boundary conditions. The power dispersion relations determined from the present scheme using 20 basis functions for each subdomain (only 60 × 60 unknowns are required) and from the FEM [18] are shown in Fig. 12 . The peak power calculated from the proposed scheme is approximately P = 82.7 μW.

Fig. 12 Power dispersion curve as a function of the effective index of the strip waveguide.

Download Full Size | PDF

The guided mode profiles for input powers P = 82.0 μW (< P_th) and P = 86.0 μW (> P_th) are shown in Figs. 13(a) and 13(b), respectively. The bistable phenomenon of the 2D nonlinear strip waveguide is similar to that found in the slab waveguide (Fig. 4). Namely, the guided mode field is mostly confined in the core region when the input power is smaller than the threshold value, and it moves to the nonlinear substrate when the input power is larger than the threshold value. This example suggests that the proposed scheme is an accurate and efficient mode solver for 2D nonlinear waveguides. For analyzing more complicated 2D waveguide structures, more rigorous consideration including the polarization dependence is certainly necessary. The full vector nonlinear mode solver based on the pseudospectral method will be reported in the forthcoming study.

Fig. 13 Mode contours for input powers (a) P = 82 μW (<P_th) and (b) P = 86 μW (>P_th) for the 2D strip waveguide with a nonlinear substrate.

Download Full Size | PDF

5. Conclusion

This study develops a novel pseudospectral mode solver for solving various nonlinear waveguides, including slab and 2D strip waveguides. In this study, two iterative approaches to the proposed scheme are implemented, and identical results are obtained. To determine the complete power dispersion curve including both the stable and unstable solutions, the guided mode power is required as the eigenvalue to be extracted. For a slab dielectric nonlinear waveguide, the accuracy of the effective index is at least of the order of 10⁻⁶ when using only 20 terms of the basis functions for each subdomain and executing ~10 iterations. For a nonlinear plasmonic waveguide composed of a thin metal film surrounded by a Kerr-like nonlinear cladding and a linear substrate, the LR- and SR-SPP modes are analyzed. The mode bistability of the LR-SPP depends on the thickness of the metal film. Therefore, the thickness of the metal film is a critical parameter for achieving the all-optical switching. Finally, the proposed scheme based on the scalar wave equation is extended to solve a 2D strip nonlinear waveguide, and the calculated results are compared with those determined by the FEM to verify the accuracy of the proposed scheme.

Acknowledgments

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. The author would also like to thank Enago (www.Enago.tw) for the English language review.

References and links

1. H. Vach, C. T. Seaton, G. I. Stegeman, and I. C. Khoo, “Observation of intensity-dependent guided waves,” Opt. Lett. 9, 238–240 (1984). [CrossRef] [PubMed]

2. I. Bennion, M. J. Goodwin, and W. J. Stewart, “Experimental nonlinear optical waveguide device,” Electron. Lett. 21(1), 41–42 (1985). [CrossRef]

3. F. Lederer, U. Langbein, and H. E. Ponath, “Nonlinear waves guided by a dielectric slab: I. TE polarization,” Appl. Phys. B 31(2), 69–73 (1983). [CrossRef]

4. F. Lederer, U. Langbein, and H. E. Ponath, “Nonlinear waves guided by a dielectric slab: II. TM–polarization,” Appl. Phys. B 31(3), 187–190 (1983). [CrossRef]

5. D. J. Robbins, “TE modes in a slab waveguide bounded by nonlinear media,” Opt. Commun. 47(5), 309–312 (1983). [CrossRef]

6. G. I. Stegeman, C. T. Seaton, J. Chilwell, and S. D. Smith, “Nonlinear waves guided by thin films,” Appl. Phys. Lett. 44(9), 830–832 (1984). [CrossRef]

7. C. T. Seaton, J. D. Valera, B. Svenson, and G. I. Stegeman, “Comparison of solutions for TM-polarized nonlinear guided waves,” Opt. Lett. 10(3), 149–150 (1985). [CrossRef] [PubMed]

8. C. T. Seaton, J. D. Valera, R. L. Shoemaker, G. I. Stegeman, J. Chilwell, and S. D. Smith, “Calculations of nonlinear TE waves guided by thin dielectric films bounded by nonlinear media,” IEEE J. Quantum Electron. 21(7), 774–783 (1985). [CrossRef]

9. A. D. Boardman, A. A. Maradudin, G. I. Stegeman, T. Twardowski, and E. M. Wright, “Exact theory of nonlinear p-polarized optical waves,” Phys. Rev. A 35(3), 1159–1164 (1987). [CrossRef] [PubMed]

10. R. I. Joseph and D. N. Christodoulides, “Exact field decomposition for TM waves in nonlinear media,” Opt. Lett. 12(10), 826–828 (1987). [CrossRef] [PubMed]

11. G. I. Stegeman, E. M. Wright, C. T. Seaton, J. V. Moloney, T. P. Shen, A. A. Maradudin, and R. F. Wallis, “Nonlinear slab-guided waves in non-Kerr-like media,” IEEE J. Quantum Electron. 22(6), 977–983 (1986). [CrossRef]

12. K. Hayata, M. Nagai, and M. Koshiba, “Finite-element formalism for nonlinear slab-guided waves,” IEEE Trans. Microw. Theory Tech. 36(7), 1207–1215 (1988). [CrossRef]

13. K. Hayata and M. Koshiba, “Full vectorial analysis of nonlinear-optical waveguides,” J. Opt. Soc. Am. B 5(12), 2494–2501 (1988). [CrossRef]

14. K. Ogusu, “TM waves guided by nonlinear planar waveguides,” IEEE Trans. Microw. Theory Tech. 37(6), 941–946 (1989). [CrossRef]

15. M. R. Ramadas, R. K. Varshney, K. Thyagarajan, and A. K. Ghatak, “A matrix approach to study the propagation characteristics of a general nonlinear planar waveguide,” J. Lightwave Technol. 7(12), 1901–1905 (1989). [CrossRef]

16. B. M. A. Rahman, J. R. Souza, and J. B. Davies, “Numerical analysis of nonlinear bistable optical waveguides,” IEEE Photon. Technol. Lett. 2(4), 265–267 (1990). [CrossRef]

17. R. D. Ettinger, F. A. Fernandez, B. M. A. Rahman, and J. D. Davies, “Vector finite element solution of saturable nonlinear strip-loaded optical waveguides,” IEEE Photon. Technol. Lett. 3(2), 147–149 (1991). [CrossRef]

18. Q. Y. Li, R. A. Sammut, and C. Pask, “Variational and finite element analyses of nonlinear strip optical waveguides,” Opt. Commun. 94(1-3), 37–43 (1992). [CrossRef]

19. A. P. Zhao and S. R. Cvetkovic, “Finite-element solution of nonlinear TM waves in multiple-quantum-well waveguides,” IEEE Photon. Technol. Lett. 4(11), 1231–1234 (1992). [CrossRef]

20. X. H. Wang and G. K. Cambrell, “Full vectorial simulation of bistability phenomena in nonlinear-optical channel waveguides,” J. Opt. Soc. Am. B 10(6), 1090–1095 (1993). [CrossRef]

21. T. Rozzi and L. Zappelli, “Modal analysis of nonlinear propagation in dielectric slab waveguide,” J. Lightwave Technol. 14(2), 229–235 (1996). [CrossRef]

22. J. G. Ma and Z. Chen, “Numerically determining the dispersion relations of nonlinear TE slab-guided waves in non-Kerr-like media,” IEEE Trans. Microw. Theory Tech. 45(7), 1113–1117 (1997). [CrossRef]

23. R. K. Varshney, I. C. Goyal, and A. K. Ghatak, “A simple and efficient numerical method to study propagation characteristics of nonlinear optical waveguides,” J. Lightwave Technol. 16(4), 697–702 (1998). [CrossRef]

24. S. S. A. Obayya, B. M. A. Rahman, K. T. V. Grattan, and H. A. El-Mikati, “Full vectorial finite–element solution of nonlinear bistable optical waveguides,’,” IEEE J. Quantum Electron. 38(8), 1120–1125 (2002). [CrossRef]

25. C. W. Kuo, S. Y. Chen, M. H. Chen, C. F. Chang, and Y. D. Wu, “Analyzing multilayer optical waveguide with all nonlinear layers,” Opt. Express 15(5), 2499–2516 (2007). [CrossRef] [PubMed]

26. S. Roy and P. Roy Chaudhuri, “Analysis of nonlinear multilayered waveguides and MQW structures: a field evolution approach using finite–difference formulation,” IEEE J. Quantum Electron. 45(4), 345–350 (2009). [CrossRef]

27. T. H. Wood, “Multiple quantum well (MQW) waveguide modulators,” J. Lightwave Technol. 6(6), 743–757 (1988). [CrossRef]

28. S. Radic, N. George, and G. P. Agrawal, “Optical switching in λ/4-shifted nonlinear periodic structures,” Opt. Lett. 19(21), 1789–1791 (1994). [CrossRef] [PubMed]

29. Y. D. Wu, “New all-optical wavelength auto-router based on spatial solitons,” Opt. Express 12(18), 4172–4177 (2004). [CrossRef] [PubMed]

30. 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]

31. A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar, “Nonlinear plasmonic slot waveguides,” Opt. Express 16(26), 21209–21214 (2008). [CrossRef] [PubMed]

32. A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar, “Quadratic phase matching in nonlinear plasmonic nanoscale waveguides,” Opt. Express 17(22), 20063–20068 (2009). [CrossRef] [PubMed]

33. 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]

34. 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]

35. J. R. Salgueiro and Y. S. Kivshar, “Nonlinear plasmonic directional couplers,” Appl. Phys. Lett. 97(8), 081106 (2010). [CrossRef]

36. 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]

37. A. Degiron and D. R. Smith, “Nonlinear long-range plasmonic waveguides,” Phys. Rev. A 82(3), 033812 (2010). [CrossRef]

38. J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springer–Verlag, 2001).

39. 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]

40. P.-J. Chiang, C.-L. Wu, C.-H. Teng, C.-S. Yang, and H.- Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008). [CrossRef]

41. 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]

42. 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]

43. C. C. Huang, “Numerical investigation of mode characteristics of nanoscale surface plasmon-polaritons using a pseudospectral scheme,” Opt. Express 18(23), 23711–23726 (2010). [CrossRef] [PubMed]

44. 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). [CrossRef] [PubMed]

45. G. I. Stegeman and C. T. Seaton, “Nonlinear surface plasmons guided by thin metal films,” Opt. Lett. 9(6), 235–237 (1984). [CrossRef] [PubMed]

46. 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]

47. P. Berini, “Plasmon-polariton waves guided by thin lossy metal films of finite width: bound modes of asymmetric structures,” Phys. Rev. B 63(12), 125417 (2001). [CrossRef]

Pseudospectral mode solver for analyzing nonlinear optical waveguides

Abstract

1. Introduction

2. Mathematical formulations

3. Computational schemes

4. Simulation results and discussion

4.1 A slab dielectric waveguide with a Kerr-like nonlinear cladding

4.2 A slab nonlinear plasmonic waveguide

4.3 A nonlinear strip waveguide with a nonlinear substrate

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (13)

Equations (43)

Optics Express