Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis

Amgad Abdrabou; A. M. Heikal; S. S. A. Obayya

doi:10.1364/OE.24.010495

1. Introduction

Optical waveguides are the fundamental building blocks around which many photonic devices are developed. Prior to the fabrication of these devices, thorough understanding of the light propagation in optical waveguides is essential. Among the methods which provide analytical or semi-analytical solutions are the mode-matching method [1, 2] and the eigenmode expansion method [3]. However, the most commonly used modal solution techniques are based on the finite difference method [4, 5] and the finite element method [6, 7]. Although successful in handling many problems, they produce matrices of huge size and have a relatively slow convergence rate of order 𝒪(1/N^p), where N is the number grid points and p is a positive constant. On the other hand, multi-domain pseudo-spectral methods (MDPSMs) [8, 9, 10] emerged recently and proved to be very accurate and efficient with a significant reduction in memory usage and computation time in addition to the fast convergence rate. The MDPSM handles semi-infinite extensions of computational domains through different approaches. One technique was proposed in [11] using Laguerre-Gauss basis functions to mimic the semi-infinite extension of the computational domain. Although this approach was successful, it required prior estimation of a very important scaling factor which is related to the decay rate of the field in outer subdomains, and significantly affects the convergence of the method. Another alternative approach was proposed in [12] by using Chebyshev functions in all subdomains while truncating the computational domain. Moreover, perfectly-matched layers (PMLs) are also commonly used to truncate the computational domain. However, there are still serious simulation problems in which PMLs collapse, confront inevitable reflections or even exponential growth [13]. Other major drawbacks are the high sensitivity to the selection of the PML parameters such as conductivity profile, thickness of the layer, permittivity, and the distance to the structure. Moreover, non-physical Berenger PML modes appear in the numerical solution that cannot be wiped out [14].

In this paper, we introduce an efficient numerical technique that can alleviate the problems associated with the existing treatments of semi-infinite computation domains in the context of MDPSMs. Our approach is based on using rational Chebyshev functions in semi-infinite subdomains. The rational basis functions are constructed by composing cardinal Chebyshev functions with a carefully selected conformal map. The proposed approach renders itself to have advantages summarized as:

Unlike Laguerre functions [11], rational Chebyshev functions are bounded i.e. no exponential function with unknown scaling is required to ensure convergence. Hence, convergence is independent of any prior estimates.
Due to the properties of rational Chebyshev functions, we don’t need to employ Perfectly-Matched Layers (PMLs) or PML-like absorbing boundary conditions.
In addition to efficient handling of guided modes, the proposed scheme is very efficient in handling leaky modes. This is accomplished by deriving a nonlinear radiation (Robin-like) boundary condition at the guide-substrate interface, which is linearized using a Taylor series expansion. Through this treatment, exact analytical representation of leaky modes is newly proposed.

2. Rational Chebyshev multi-domain pseudo-spectral method

The wave equation for 2D optical waveguides is

\frac{\partial^{2} u}{\partial y^{2}} - θ \frac{1}{n^{2}} \frac{\partial n^{2}}{\partial y} \frac{\partial u}{\partial y} + k_{0}^{2} (n^{2} (y) - n_{eff}^{2}) u = 0,

where

u (y) = {\begin{array}{l} E_{x} (y) and θ = 0 & for TE modes; \\ H_{x} (y) and θ = 1 & for TM modes . \end{array}

To deal with the discontinuity of n²(y), the whole domain is first divided into a finite number of subdomains at interfaces of discontinuity denoted by d_i, i = 1, 2,...,M for a structure having M interfaces of different materials as shown in Fig. 1. This leads to a computational domain with (M + 1) subdomains, Ω_i, i = 1, 2,...,M + 1. The interface boundary conditions at the points of discontinuity, d_i are given by

E_{x} (d_{i}^{-}) = E_{x} (d_{i}^{+}), \frac{\partial E_{x}}{\partial y} (d_{i}^{-}) = \frac{\partial E_{x}}{\partial y} (d_{i}^{+}),

for TE modes, and

H_{x} (d_{i}^{-}) = H_{x} (d_{i}^{+}), (\frac{1}{n^{2}} \frac{\partial H_{x}}{\partial y}) (d_{i}^{-}) = (\frac{1}{n^{2}} \frac{\partial H_{x}}{\partial y}) (d_{i}^{+}),

for TM modes, where

d_{i}^{-}

and

d_{i}^{+}

refer to locations infinitesimally to the left and to the right of d_i, respectively.

Fig. 1 Schematic diagram of a multi-layer 2D planar waveguide.

Download Full Size | PDF

2.1. Rational Chebyshev functions with algebraic maps

The proposed method is based on the cardinal Chebyshev function expansion [11, 15] defined for ỹ ∈ [−1, 1] as follows

C_{j} (\tilde{y}) = \frac{{(- 1)}^{N - j + 1} (1 - {\tilde{y}}^{2}) T_{N}^{'} (\tilde{y})}{σ_{j} N^{2} (\tilde{y} - {\tilde{y}}_{j})}, \tilde{y} \neq {\tilde{y}}_{j},

where T_N(ỹ) = cos(N cos⁻¹(ỹ)) is the Chebyshev polynomial of degree N, σ₀ = σ_N = 2, and σ_j = 1 for j = 1,...,N − 1. The basis functions in Eq. (5) are defined based on the Gauss-Lobatto collocation points ỹ_k = −cos(kπ/N), k = 0, 1,..,N. For simplicity of reference, we define the subdomain Ω₁ = (−∞, d₁), Ω_i = (d_i−1, d_i), i = 2,...,M, and Ω_M+1 = (d_M, ∞).

Since the cardinal Chebyshev basis functions are defined on the interval [−1, 1], conformal map are needed in order to define suitable basis functions in each subdomain, Ω_i shown in Fig. 1.

Referring to Fig. 1, we have two types of subdomains. First, bounded inner subdomains, Ω₂,...,Ω_M and two semi-infinite subdomains Ω₁ and Ω_M+1.

For the bounded subdomains, the maps are simply linear and take the form $ϕ_{i} (y) = 2 (\frac{y - d_{i - 1}}{d_{i} - d_{i - 1}}) - 1, 2 \leq i \leq M,$ and their mapped Gauss-Lobatto collocation points are given by $y_{j}^{(i)} = ϕ_{i}^{- 1} ({\tilde{y}}_{j}) = d_{i - 1} + (\frac{d_{i} - d_{i - 1}}{2}) (1 - cos (j π / N_{i})) .$
For the semi-infinite subdomains Ω₁ ≡ (−∞, d₁) and Ω_M+1 ≡ (d_M, ∞), we define the algebraic maps ϕ₁(y; N₁) and ϕ_M+1(y; N_M+1) in Ω₁ and Ω_M+1, respectively. In Ω₁ the map ϕ₁(y; N₁) is given by [16] $ϕ_{1} (y; N_{1}) = \frac{τ (N_{1}) + (y - d_{1})}{τ (N_{1}) - (y - d_{1})},$ where τ(N₁) = κ₁ N₁ is an linear function with κ₁ > 0 and N₁ + 1 is the number of basis functions used in Ω₁. The new collocation points, $y_{j}^{(1)}$ in the subdomain Ω₁ are given by $y_{j}^{(1)} = d_{1} - τ (N_{1}) \frac{1 + cos (j π / N_{1})}{1 - cos (j π / N_{1})},$ at which we obtain the physical values of the electric or magnetic field components. Similarly, for the semi-infinite subdomain Ω_M+1, we employ the map ϕ_M+1(y; N_M+1) given by $ϕ_{M + 1} (y; N_{M + 1}) = \frac{(y - d_{m}) - α (N_{M + 1})}{(y - d_{m}) + α (N_{M + 1})},$ where α(N_M+1) = κ₂ N_M+1, κ₂ > 0. The corresponding collocation points, $y_{j}^{(M + 1)}$ in the subdomain Ω_M+1 are then given by $y_{j}^{(M + 1)} = d_{m} + α (N_{M + 1}) \frac{1 - cos (j π N_{M + 1})}{1 + cos (j π N_{M + 1})} .$

The solution of the wave Eq. (1) is obtained by expanding the field in terms of mapped basis functions obtained by composing Chebyshev functions by the suitable conformal map for each subdomain Ω_i, as we have described and imposing Dirichlet zero-boundary conditions at y = ±∞ as shown in Fig. (2).

Fig. 2 Application of zero-boundary conditions for guided modes at y = ±∞.

Download Full Size | PDF

The field expansion in a subdomain Ω_i takes the form

u^{(i)} (y) = \sum_{j = 0}^{N_{i}} u_{j}^{(i)} (C_{j} ○ ϕ_{i}) (y), y \in Ω_{i},

where N_i + 1 denotes the number of basis functions in a subdomain Ω_i and the circle “○” in Eq. (12) denotes a compositions of the functions C_j(ỹ) with the map ϕ_i(y), namely (C_j ○ ϕ_i)(y) = C_j (ϕ_i(y)). The maps ϕ_i(y) all satisfy ϕ_i(Ω_i) = [−1, 1] i.e. the function ϕ_i maps the subdomain Ω_i onto the interval [−1, 1]. The unknown coefficients

u_{j}^{(i)}

in Eq. (12) represent the field values in the subdomain Ω_i at the mapped Gauss-Lobatto collocation points for that subdomain, given by

y_{j}^{(i)} = ϕ_{i}^{- 1} ({\tilde{y}}_{j})

. Cardinality of the basis functions in Eq. (12) is preserved i.e. we have

(C_{j} ○ ϕ_{i}) (y_{k}^{(i)}) = δ_{k j}

, where δ_kj is the Kronecker delta.

It is worth noting that composing the algebraic maps ϕ₁(y) and ϕ_M+1(y) defined by Eqs. (8) and (10), respectively, leads to rational Chebyshev basis functions defined on the semi-infinite subdomains Ω₁ and Ω_M+1. The proposed treatment accounts for semi-infinite extensions of the outer subdomains, hence non-physical modes have no chance to appear in the solution spectrum. It also leads to an efficient replacement of the non-physical perfectly-matched layers (PMLs).

The multi-domain pseudo-spectral method reported [12] applies a domain truncation approach. In that approach, the outer semi-infinite subdomains Ω₁ and Ω_M+1 shown in Fig. 1 are truncated by replacing the original computational domain (−∞, ∞) by a computational window [−W, W] for W > max{|d₁|, |d_M|} and imposing Dirichlet zero-boundary conditions at y = ±W. Truncation of the domain assumes rapid decay rate of the field in outer subdomains. Although this approach is successful, the appropriate value of W should vary according to the refractive index profile (RIP) to accurately describe the field in outer regions. Further, the value of W should increase as the number of basis functions in the outermost subdomains increases to maintain high convergence rate (see [16] and references therein).

Another approach to handle the subdomains Ω₁ and Ω_M+1 was proposed in [11]. It was based on using Laguerre-Gauss basis functions. However, that approach requires a priori estimate for a scaling parameter whose value significantly affects both the accuracy and rate of convergence.

2.2. Leaky modes treatment

For multi-layer planar waveguides where the substrate has a high refractive index, leaky modes arise and accurate calculation of the effective index imaginary part is challenging. Instead of expanding the field in the outer subdomain representing the substrate with high refractive index, we use the analytic expression of the field u(y) in the substrate region. From this expression we derive the suitable boundary conditions at the guide-substrate interface. This approach can be described follows

For a substrate with refractive index n_s in the subdomain Ω_M+1, the field u(y) denoted by u^(M+1) can be expressed as

u^{(M + 1)} (y) = u_{0} e^{- j γ_{s} (y - y_{M})},

where

γ_{s} = \sqrt{k_{0}^{2} n_{s}^{2} - β^{2}}

. Differentiating Eq. (13) with respect to y and evaluating at

y = d_{M}^{-}

, yield

{\frac{d u^{(M + 1)}}{d y} |}_{y = d_{M}^{-}} = - j \sqrt{k_{0}^{2} n_{s}^{2} - β^{2}} {u^{(M + 1)} |}_{y = d_{M}^{-}}

Since the condition

| β^{2} / k_{0}^{2} n_{s}^{2} | < < 1

, the following approximation is valid

γ_{s} = \sqrt{k_{0}^{2} n_{s}^{2} - β^{2}} ≃ k_{0} n_{s} (1 - \frac{1}{2} \frac{β^{2}}{k_{0}^{2} n_{s}^{2}} + \dots)

Substituting Eq. (15) into Eq. (14) yields the following condition imposed at the guide-substrate interface as shown in Fig. (3)

\frac{d u^{(M + 1)}}{d y} + {j k_{0} n_{s} u^{(M + 1)} |}_{y = d_{M}^{-}} = \frac{j}{2} \frac{β^{2}}{k_{0}^{2} n_{s}^{2}} {u^{(M + 1)} |}_{y = d_{M}^{-}}

The treatment of the case of a substrate in the subdomain Ω₁ is similar and yields the following condition

\frac{d u^{(1)}}{d y} - j k_{0} n_{s} {u^{(1)} |}_{y = d_{1} +} = - \frac{j}{2} \frac{β^{2}}{k_{0}^{2} n_{s}^{2}} {u^{(1)} |}_{y = d_{1}^{+}}

Unlike the treatment based on Mur’s absorbing boundary conditions [17], our approach leads to a linear problem while maintaining an accurate estimation of the complex propagation constants. However, if we have

| β^{2} / k_{0}^{2} n_{s}^{2} | < 1

, i.e. the value of

| β^{2} / k_{0}^{2} n_{s}^{2} |

is not small enough, we may need to consider higher-order terms of Eq. (15) yielding a polynomial eigenvalue problem of an arbitrary degree p > 1

A_{0} - λ A_{1} u + λ^{2} A_{2} u + \dots + λ^{p} A_{p} u = 0, λ = β^{2} .

Equation (18) can be easily solved by existing software subroutines (see also [18, 19]).

Fig. 3 Application of leaky modes condition of Eq. (16) at the guide-substrate interface.

Download Full Size | PDF

2.3. Formation of the linear eigenvalue problem

A linear eigenvalue problem for modal solution is deduced by substituting Eq. (12) into the wave equation Eq. (1), and demanding that Eq. (1) is satisfied in each subdomain Ω_i at the corresponding set of collocation points $y_{j}^{(i)}$ except for the boundary points of each subdomain which are used to impose the conditions Eq. (3) and Eq. (4). For guided modes, zeros-boundary conditions are imposed at y = ±∞, while the conditions Eq. (16) or Eq. (17) hold for leaky modes. Define the first order differentiation matrix D ≡ D⁽¹⁾ such that its kj-th entry d_kj = dC_j/dỹ(ỹ_k), defined in [20]. The differentiation matrix of the second order derivative is then given by

D^{(2)} = D D = D^{2} .

For a subdomain Ω_i, we define the diagonal matrix Λ_i(f) for an arbitrary function f by

Λ_{i} (f) = diag (f (y_{0}^{(j)}), \dots, f (y_{N_{i}}^{(i)})) .

Then, differentiation matrices for Ω_i are given by

D_{i} = Λ_{i} (\frac{d ϕ_{i}}{d y}) D,

and

D_{i}^{(2)} = Λ_{i} ({(\frac{d ϕ_{i}}{d y})}^{2}) D^{2} + Λ_{i} (\frac{d^{2} ϕ_{i}}{d y^{2}}) D .

In each subdomain Ω_i, we define the restricted matrices

{\hat{I}}_{i} = {[I_{i}]}_{k j}, 1 \leq k \leq N_{i} - 1, 0 \leq j \leq N_{i},

{\hat{D}}_{i} = {[D_{i}]}_{k j}, 1 \leq k \leq N_{i} - 1, 0 \leq j \leq N_{i},

and

{\hat{D}}_{i}^{(2)} = {[D_{i}^{(2)}]}_{k j}, 1 \leq k \leq N_{i} - 1, 0 \leq j \leq N_{i},

where the matrix I_i in Eq. (23) is an identity matrix of size (N_i + 1) × (N_i + 1). Using the previous notations and definitions, we get for each subdomain Ω_i the eigenvalue problem

[K] u^{(i)} = β^{2} [{\hat{I}}_{i}] u^{(i)}

where the (N_i − 1) × (N_i + 1) matrix K is given by

K = {\hat{D}}_{i}^{(2)} - θ Λ_{i} (\frac{1}{n^{2}} \frac{\partial n^{2}}{\partial y}) {\hat{D}}_{i} + Λ_{i} (k_{0}^{2} n^{2} (y)) {\hat{I}}_{i} .

By assembling the discrete systems of all subdomains and imposing the interface conditions we get a linear eigenvalue problem which is solved for the propagation constants and the field values.

It remains to point out that extension of the proposed approach to handle 3D problems has three main points that should be considered. First, rational Chebyshev basis functions would be incorporated in subdomains with semi-infinite extensions. Second, assuming full-vectorial formulation in terms of the magnetic field components H_x and H_y, the form of approximate expansion of the field components is a double sum of products of basis functions of the form

H_{τ} (x, y) = \sum_{i = 0}^{N} \sum_{j = 0}^{M} h_{i j}^{τ} C_{i} (ϕ (x)) C_{j} (ψ (y)), τ = x or y,

where

h_{i j}^{τ} = H_{τ} (x_{i}, y_{j})

are the field values at the set of collocation points,(x_i, y_j), and (N +1) and (M + 1) denote the number of basis functions in the direction of the x-axis and y-axis, respectively. Finally, the conditions imposed at horizontal and vertical interfaces are the continuity of the field components E_z and H_z (see [10] and references therein).

3. Numerical results

To demonstrate the efficiency of our approach, we used our suggested approach to perform modal solution of some challenging optical waveguides.

3.1. Symmetric step-index slab waveguide

The first example is a three layer symmetric single-mode step-index slab waveguide [21]. It consists of core with refractive index n_c = 1.5 sandwiched between a cladding with refractive index n_cl = 1.48. The guiding layer width is 2 μm and the operating wavelength is 1.5 μm. The analytic value of the effective refractive index, $n_{eff}^{exact} = 1.489305274606916$ is obtained by solving the characteristic equation for TM modes using the analytical method described in [22]. The normalized propagation constant b is defined by $b = (n_{eff}^{2} - n_{cl}^{2}) / (n_{c}^{2} - n_{cl}^{2})$ [23]. Its exact value is b^exact = 0.455715434937231. The relative error in the normalized propagation constant is given by

Δ b = \frac{| b - b^{exact} |}{b^{exact}}

The values of Δb are listed Table 1 as the total number of basis function, N is increased.

Table 1. Comparison of Normalized Propagation Constant and Relative Error for Step-index Slab Waveguide at λ = 1.5μm

View Table | View all tables in this article

The proposed method outperforms the Matrix Method [21], where the minimum obtained relative error is only 1.818E–04 using a truncated series of N = 400 terms, while the proposed method provides a relative error of 4.909E–14 at N = 80 (five times less than) matrix method. The semi-log plot of Δb versus N is shown in Fig. (4b) where a very high convergence rate is achieved. This assures that the interface boundary conditions are considered accurately. The normalized fundamental TM field, |H_x| is shown in Fig. (4a).

Fig. 4 (a) Normalized fundamental TM field, H_x; (b) Convergence of the normalized propagation constant.

Download Full Size | PDF

3.2. Planar waveguide with continuous Gaussian RIP

We consider a planar waveguide with continuous Gaussian RIP [12, 24] shown in Fig. (6a). It is described by

n^{2} (y) = n_{s}^{2} + 2 n_{s} C e^{- y^{2} / t^{2}} + C^{2} e^{- 2 y^{2} / t^{2}}

where n_s = 1.512, C = 0.0833, t = 2.66μm and the operating wavelength λ = 0.6328μm. To apply our technique, we divide the computation domain (−∞, ∞) into two semi-infinite subdomains as shown in Fig. (5).

Fig. 5 Collocation points distribution for proposed method for M = 1, N₁ = N_M+1 = 8.

Download Full Size | PDF

Fig. 6 (a) Continuous Gaussian refractive index profile.; (b) First four TE mode profiles obtained by the proposed method.

Download Full Size | PDF

It is worth noting that the collocation points have high concentration near the center of the waveguide where field is more confined. The treatment in [12] is based on single domain expansion where a domain truncation technique was used to truncate the computational domain from (−∞, ∞) to a finite computational window [−W, W]. Although this approach was successful, the appropriate value of L should vary as different graded-index profiles are considered. This problem is handled here by dividing the computational domain into two semi-infinite subdomains. In each subdomain, the field is expanded using rational Chebyshev functions, taking into account the semi-infinite extension. Convergence of the effective refractive index of the first TE mode is shown in Table 2 with the third column representing an approximate relative error defined by $| n_{eff}^{new} - n_{eff}^{old} | / n_{eff}^{new}$ and the first ten TE modes are listed in Table 3 and the first four TE modes are shown in Fig. (6b).

Table 2. Convergence of the Effective Refractive Index of the First TE Mode

View Table | View all tables in this article

Table 3. Comparison of the Effective Indices for the TE Modes at λ = 0.6328 μm

View Table | View all tables in this article

The results reported in Table 3 show a great agreement of the present method with method in [12]. However, the former method based on domain truncation needs appropriate estimation of the size of the truncated domain i.e. the value of the parameter W, especially for cases when the field decays slowly as y → ±∞.

3.3. Discontinuous planar waveguide with strong discontinuity

The third example studies the case of a discontinuous RIP given by

n^{2} (y) = {\begin{array}{l} n_{s}^{2} + 2 n_{s} C e^{- y^{2} / t^{2}} + C^{2} e^{- 2 y^{2} / t^{2}}, & y \geq 0; \\ 1, & y < 0, \end{array}

where n_s = 1.512, C = 0.0833, t = 2.66μm and λ = 0.6328μm. This refractive index profile is in Fig. (7a). The field profile of the fundamental TE mode obtained by the proposed method is shown in Fig. 7b. This example was studied in [11] using a pseudo-spectral domain decomposition with Laguerre-Gauss functions used in the outer subdomains and using the same approach with domain truncation in [12]. The results shown in Table 4 depict the convergence behaviour of the proposed method as the total number of basis functions increases. The approximate relative error is also given. The obtained convergent value by the proposed method is in good agreement with that obtained in [11] which is 1.579107.

Fig. 7 (a) Discontinuous refractive index profile.; (b) Normalized fundamental TE field for the discontinuous RIP.

Download Full Size | PDF

Table 4. Convergence of the Effective Index for the Discontinuous Planar Waveguide at λ = 0.6328μm

View Table | View all tables in this article

3.4. Slow plasmon resonant nano-structure

We now turn to studying plasmonic structures. First, we study slow surface plasmon-polariton (SPP) modes for a symmetric insulator-metal-insulator (IMI) structure shown in Fig. 8. The propagation constants of long (short) range L(S)R-SPP modes, k_l(s)rsp are given, respectively by the exact formula given in [25]. The dielectric functions of the metal and dielectric are given by ε_d = 1 and ε_m = −23 + 1.69 j, respectively. The effective refractive index for long (short) range SPP modes n_l(s) and the propagation length L_l(s) are then given by

n_{l (s)} = ℜ (k_{l (s) r s p}) / k_{0}, L_{l (s)} = {[2 ℑ (k_{l (s) r s p})]}^{- 1} .

The relative error in the mode effective index δn_l(s) and the relative error in the propagation length δL_l(s) are defined by

δ n_{l (s)} = \frac{| n_{l (s)} - n_{l (s)}^{exact} |}{n_{l (s)}^{exact}}, δ L_{l (s)} = \frac{| L_{l (s)} - L_{l (s)}^{exact} |}{L_{l (s)}^{exact}} .

We consider a case of a film with thickness t = 40 nm. The effective index and propagation length for the long range SPP mode for are reported in Table 5 and those for short range SPP mode are given in Table 6. The field profiles for the normalized ℜ(H_x(y)) for the LR-SPP and SR-SPP modes are shown in Fig. 8. The rapidly decaying relative errors for both the effective index and propagation length reflect the excellent performance and accuracy of the proposed method. In Table 7, films with different thickness values are reported, where N is the number of basis functions in each subdomain.

Table 5. Convergence of n_l and L_l Values for t = 40nm

View Table | View all tables in this article

Table 6. Convergence of n_s and L_s Values for Film Thickness t = 40nm

View Table | View all tables in this article

Fig. 8 Long and short SPP modes in the symmetric IMI nano-structure.

Download Full Size | PDF

Table 7. Relative Error for Different Metal Thickness at N = 100

View Table | View all tables in this article

3.5. The gap SPP mode

We now consider the calculation of the gap SPP (G-SPP) mode in the symmetric metal-insulator metal(MIM) structure shown in Fig. 9 with the field plot of the G-SPP mode. This mode is the only surviving mode for all values of gap thickness, t (nm). The G-SPP propagation constant k_gsp is given by the exact expression [26, 27]

tanh (\frac{α_{d} t}{2}) = - \frac{ε_{d} α_{m}}{ε_{m} α_{d}}

where

α_{m, d} = \sqrt{k_{gsp}^{2} - k_{0}^{2} ε_{m, d}}

. The effective mode index and propagation length are defined as in Eq. (32). The convergence of the effective index, n_gsp and the propagation length L_gsp with the corresponding relative errors given, respectively by δn_gsp and δL_gsp are reported in Table 8 for gap of thicknesses t = 20nm and those for t = 1000nm are reported in Table 9. The exact values were obtained by solving Eq. (34).

Fig. 9 G-SPP mode in the symmetric MIM nano-structure.

Download Full Size | PDF

Table 8. Convergence of n_gsp and L_gsp for Gap Thickness t = 20nm

View Table | View all tables in this article

Table 9. Convergence of n_gsp and L_gsp for Gap Thickness t = 1000nm

View Table | View all tables in this article

In spite of the strong damping of the field inside the metal region in the MIM nano-structure, the convergence of the domain truncation method [12] depends on the size of the computational window [−W, W]. In addition, the numerical results shown in Fig. (10a) and (10b) confirm that convergence of the domain truncation approach cannot be ensured by just using a large value of the computational window parameter W.

Fig. 10 Comparison of (a) the effective refractive index and (b) the propagation length; calculated by the domain truncation method [12] with variable computational window [−W, W] and the proposed method with fixed computational window, W = ∞ for gap of thickness t = 20 nm.

Download Full Size | PDF

3.6. Anti-resonant reflecting optical waveguide (ARROW)

Finally, we consider an anti-resonant reflecting optical waveguide (ARROW) structure [28, 29] as shown in Fig. 11a. The operating wavelength λ = 0.6328 μm, the refractive indices are n_a = 1 for air and n_s = 3.85 for the substrate. The refractive indices for the first and second claddings are n₁ = 2.3 and n₂ = 1.46, respectively.

Fig. 11 (a) Anti-resonant reflecting optical waveguide structure; (b) Normalized TE₁ and TE₁₈ modes at λ = 0.6328μm.

Download Full Size | PDF

Their thicknesses are d₁ = 0.142λ and d₂ = 3.15λ, respectively. The refractive index of the core layer is n_c = 1.46 and its thickness d_c = 6.3λ. The results of the normalized propagation constant, β/k₀ are reported in Table 10. The normalized transverse electric field E_x corresponding to TE₁ and TE₁₈ leaky modes are shown in Fig. 11b. The inset figures demonstrate the effective handling of interface boundary conditions.

Table 10. Normalized Propagation Constant, β/k₀ for TE Modes in ARROW

View Table | View all tables in this article

4. Conclusion

In conclusion, an efficient multi-domain pseudo-spectral method based on rational Chebyshev functions is introduced. Careful selection of the rational maps and employing the rational Chebyshev functions had a significant improvement of the accuracy over traditional domain truncation techniques. The efficiency of our approach was demonstrated through accurate calculation of plasmonic and leaky modes propagation constants. The reported results showed the precision and accuracy of the proposed method with an excellent agreement with other methods in the literature.

References and links

1. Y. C. Shih, “The mode-matching method,” in Numerical Techniques for Microwave and Millimeter-Wave Passive Structures, T. Itoh, ed. (Wiley, 1989), pp. 592–621.

2. A. S. Sudbo, “Film mode matching: a versatile numerical method for vector mode field calculations in dielectric waveguides,” Pure Appl. Opt. 2(3), 211–233 (1993). [CrossRef]

3. P. Bienstman and R. Baets, “Optical modelling of photonic crystals and VCSELs using eigenmode expansion and perfectly matched layers,” Opt. Quantum Electron. 33, 327–341 (2001). [CrossRef]

4. Y. P. Chiou, Y. C. Chiang, C. H. Lai, C. H. Du, and H. C. Chang, “Finite difference modeling of dielectric waveguides with corners and slanted facets,” J. Lightwave Technol. 27, 2077–2086 (2009). [CrossRef]

5. N. Thomas, P. Sewell, and T. M. Benson, “A new full-vectorial higher order finite-difference scheme for the modal analysis of rectangular dielectric waveguides,” J. Lightwave Technol. 25, 2563–2570 (2007). [CrossRef]

6. S. Selleri, L. Vincetti, A. Cucinotta, and M. Zoboli, “Complex FEM modal solver of optical waveguides with PML boundary conditions,” Opt. Quantum Electron. 33, 359–371 (2001). [CrossRef]

7. Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11, 1545–1549 (1993). [CrossRef]

8. C. C. Huang, “Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method,” Opt. Express 14, 3363–3378 (2011). [CrossRef]

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

10. 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, 457–465 (2005). [CrossRef]

11. C. C. Huang, C. C. Huang, and J. Y. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profile using spectral collocation method with domain decomposition,” J. Lightwave Technol. 20, 2284–2296 (2003). [CrossRef]

12. J. Zhu, X. Zhang, and R. Song, “A unified mode solver for optical waveguides based on mapped barycentric rational chebyshev differentiation matrix,” J. Lightwave Technol. 28, 1802–1810 (2010). [CrossRef]

13. A. F. Oskooi, L. Zhang, Y. Avniel, and S. G. Johnson, “The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers,” Opt. Express 16, 11376–11392 (2008). [CrossRef] [PubMed]

14. H. Rogier and D. D. Zutter, “Berenger and leaky modes in optical fibers terminated with a perfectly matched layer,” J. Lightwave Technol. 20, 1141–1148 (2002). [CrossRef]

15. D. B. Haidvogel and T. Zang, “The accurate solution of Poisson equation by expansion in Chebyshev polynomials,” J. Comput. Phys. 30, 167–180 (1979). [CrossRef]

16. J. P. Boyd, “The optimization of convergence for Chebyshev polynomial methods in an unbounded domain,” J. Comput. Phys. 45, 42–79 (1982). [CrossRef]

17. C. C. Huang, “Numerical calculations of ARROW structures by pseudo-spectral approach with Murs absorbing boundary conditions,” Opt. Express 14, 11631–11652 (2006). [CrossRef] [PubMed]

18. F. Tisseur and Karl Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. 43, 235–286 (2001). [CrossRef]

19. J. -P. Dedieu and F. Tisseur, “Perturbation theory for homogeneous polynomial eigenvalue problems,” Linear Algebra Appl. 358, 71–94 (2003). [CrossRef]

20. G.-X. Fan, Q. H. Liu, and J. S. Hesthaven, “Multidomain pseudospectral time-domain simulations of scattering by objects buried in lossy media,” IEEE Trans. Geosci. Remote Sens. 40, 1366–1373 (2002). [CrossRef]

21. L. Wang and C. S. Hsiao, “A matrix method for studying TM modes of optical planar waveguides with arbitrary index profiles,” IEEE J. Quantum Electron. 37, 1654–1660 (2001). [CrossRef]

22. K. Kawano and T. Kitoh, “Analytical Methods,” in Introduction to Optical Waveguide Analysis: Solving Maxwell’s Equations and the Schrödinger Equation (John Wiley & Sons, 2004), pp. 15–18.

23. C. Vassallo, “1993–1995 Optical mode solvers,” Opt. Quantum Electron. 29, 95–114 (1997). [CrossRef]

24. S. Banerjee and A. Sharma, “Propagation characteristics of optical waveguiding structure by direct solution of the helmholtz equation for total fields,” J. Opt. Soc. Amer. A 6, 1884–1894 (1989). [CrossRef]

25. S. I. Bozhevolnyi and T. Sondergaard, “General properties of slow-plasmon resonant nanostructures: nano-antennas and resonators,” Opt. Express 15, 10869–10877 (2007). [CrossRef] [PubMed]

26. E. N. Economou, “Surface plasmons in thin films,” Phys. Rev. 182, 539–554 (1969). [CrossRef]

27. R. Zai, M. D. Selker, P. B. Catrysee, and M. L. Brongersma, “Geometries and materials for subwavelength surface plasmon modes,” J. Opt. Soc. Am. A (21), 2442–2446 (2004).

28. E. Anemogiannis, E. N. Glytsis, and T. K. Gaylord, “Determination of guided and leaky modes in lossless and lossy planar multilayer optical waveguides: reflection pole method and wavevector density method,” J. Lightwave Technol. 17, 929–941 (1999). [CrossRef]

29. E. I. Golant and K. M. Golant, “New method for calculating the spectra and radiation losses of leaky waves in multilayer optical waveguides,” Tech. Phys. 51, 1060–1068 (2006). [CrossRef]

b^exact	N	b	Δb
0.455715434937231	30	0.45570782546737	1.670E–05
	40	0.45571561539621	3.960E-07
	50	0.45571543880946	8.497E-09
	60	0.45571543474805	4.151E-10
	70	0.45571543493613	2.420E-12
	80	0.45571543493721	4.909E-14

N	n_eff	Approx. Relative Error
20	1.58912647924864	–
30	1.58934919934532	1.4013E–04
40	1.58935080298308	1.0090E–06
50	1.58934827515103	1.5905E–06
60	1.58934847814690	1.2772E–07
70	1.58934847308759	3.1833E–09
80	1.58934847224956	5.2728E–10
90	1.58934847242315	1.0922E–10
100	1.58934847240226	1.3144E–11
200	1.58934847240402	1.1075E–12

TE_i	Ref.[24]	Ref.[12]	This work
0	1.58935	1.5893484724038	1.58934847240942
1	1.57775	1.5777480795564	1.57774807955646
2	1.56682	1.5668208288289	1.56682082894205
3	1.55661	1.5566114709170	1.55661147091710
4	1.54717	1.5471718524290	1.54717185318148
5	1.53856	1.5385637047797	1.53856370477617
6	1.53086	1.5308633148130	1.53086331716541
7	1.52417	1.5241700562622	1.52417005624192
8	1.51862	1.5186240734879	1.51862407710320
9	1.51444	1.5144512674959	1.51445126746200

N	n_eff	Approx. Relative Error
30	1.57908597337832	—
40	1.57911013095775	1.52982E–05
50	1.57910781285803	1.46798E–06
60	1.57910748581268	2.07108E–07
70	1.57910746230120	1.48891E–08
80	1.57910746088772	8.95113E–10
90	1.57910746079261	6.02302E–11
100	1.57910746078507	4.77497E–12

N	n_l	δn_l		L_l	δL_l
20	1.00971418429127	4.791E–05		216.716256272217	2.531E–02
30	1.00976098182220	1.568E–06		222.073693701710	1.217E–03
40	1.00976249406630	7.070E–08		222.330033609942	6.459E–05
50	1.00976256206769	3.356E–09		222.343642016321	3.388E–06
60	1.00976256529302	1.615E–10		222.344356651303	1.743E–07
70	1.00976256544818	7.817E–12		222.344393451271	8.822E–09
80	1.00976256545569	3.776E–13		222.344395314153	4.439E–10
90	1.00976256545609	1.297E–14		222.344395405464	3.322E–11
100	1.00976256545607	4.398E–16		222.344395412391	2.066E–12

LR exact values [25]	$n_{l}^{exact}$		$L_{l}^{exact}$

	1.00976256545607		222.34439541285

Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis

Abstract

1. Introduction

2. Rational Chebyshev multi-domain pseudo-spectral method

2.1. Rational Chebyshev functions with algebraic maps

2.2. Leaky modes treatment

2.3. Formation of the linear eigenvalue problem

3. Numerical results

3.1. Symmetric step-index slab waveguide

3.2. Planar waveguide with continuous Gaussian RIP

3.3. Discontinuous planar waveguide with strong discontinuity

3.4. Slow plasmon resonant nano-structure

3.5. The gap SPP mode

3.6. Anti-resonant reflecting optical waveguide (ARROW)

4. Conclusion

References and links

Cited By

Figures (11)

Tables (10)

Equations (34)

Optics Express

N	n_s	δn_s	L_s	δL_s
20	1.04781815201404	2.243E–04	10.9179294176588	6.555E–03
30	1.04804985115533	3.211E–06	10.8448529229272	1.816E–04
40	1.04805329347435	7.325E–08	10.8467650925917	5.355E–06
50	1.04805321766548	9.128E–10	10.8468267199563	3.266E–07
60	1.04805321661594	8.868E–11	10.8468231842769	6.168E–10
70	1.04805321670942	5.201E–13	10.8468231729868	4.240E–10
80	1.04805321670896	8.369E–14	10.8468231776607	6.866E–12
90	1.04805321670892	4.873E–14	10.8468231775793	6.333E–13
100	1.04805321670887	7.203E–15	10.8468231775848	1.225E–13

SR exact values [25]	$n_{s}^{exact}$		$L_{s}^{exact}$

	1.04805321670887		10.8468231775862

t (nm)	Long-range	SPP mode	Short-range	SPP mode

	δn_l	δL_l	δn_s	δL_s
20	3.453E – 14	8.206E – 11	1.325E – 14	2.639E – 13
60	1.312E – 15	7.170E – 13	8.831E – 15	3.370E – 13
100	4.353E – 15	4.157E – 12	8.679E – 15	3.124E – 12
140	8.391E – 14	1.462E – 11	9.255E – 14	1.214E – 11

N	n_gsp	δn_gsp	L_gsp	δL_gsp
20	1.93475631429775	4.591E–05	2.01347998994764	3.834E–03
30	1.93466312749892	2.254E–06	2.02148985865666	1.286E–04
40	1.93466765839929	8.825E–08	2.02122179565734	4.028E–06
50	1.93466748157260	3.147E–09	2.02123017772140	1.195E–07
60	1.93466748786526	1.052E–10	2.02122992948593	3.350E–09
70	1.93466748766451	1.431E–12	2.02122993643811	8.976E–11
80	1.93466748765857	1.639E–12	2.02122993624751	4.541E–12

G-SPP exact values [26]	$n_{gsp}^{exact}$		$L_{gsp}^{exact}$

	1.93466748766174		2.02122993625668

N	n_gsp	δn_gsp		L_gsp	δL_gsp
20	1.03468285611851	9.509E–05		36.0933811311884	1.261E–02
30	1.03458428715564	1.790E–07		36.5580369319735	9.641E–05
40	1.03458449527045	2.216E–08		36.5543697314078	3.908E–06
50	1.03458447169997	6.240E–10		36.5545148706538	6.233E–08
60	1.03458447236412	1.792E–11		36.5545125901131	5.584E–11
70	1.03458447234510	4.672E–13		36.5545125899287	6.089E–11
80	1.03458447234559	1.459E–14		36.5545125922657	3.044E–12

G-SPP exact values [26]	$n_{gsp}^{exact}$		$L_{gsp}^{exact}$

	1.03458447234558		36.5545125921544

Mode	Ref. [28]	Ref. [29]	This work

TE₁	1.45794 – 5.4189E–08 j	1.45794 – 5.4196E–08 j	1.45794 – 5.4030E–08 j
TE₂	1.45192 – 5.2871E–05 j	1.45192 – 5.1991E–05 j	1.45192 – 5.2731E–05 j
TE₃	1.45117 – 1.9203E–04 j	1.45118 – 1.8852E–04 j	1.45118 – 1.9146E–04 j
TE₉	1.37900 – 1.7263E–03 j	1.37901 – 1.6885E–03 j	1.37900 – 1.7220E–03 j
TE₁₀	1.35567 – 1.2860E–04 j	1.35567 – 1.2933E–04 j	1.35567 – 1.2831E–04 j
TE₁₁	1.32510 – 8.9350E–04 j	1.32510 – 8.8172E–04 j	1.32510 – 8.9187E–04 j
TE₁₇	1.13592 – 2.1513E–03 j	1.13593 – 2.1576E–03 j	1.13593 – 2.1492E–03 j
TE₁₈	1.10700 – 7.4694E–03 j	1.10702 – 7.2763E–03 j	1.10700 – 7.4618E–03 j
TE₁₉	1.06241 – 1.6648E–03 j	1.06241 – 1.7338E–03 j	1.06241 – 1.6633E–03 j