Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach

Chia-Chien Huang

doi:10.1364/OE.18.026583

1. Introduction

Simulations of dielectric waveguides provide crucial information for designing photonic integrated circuits. In particular, most functional photonic devices, such as isolators, switches, modulators, couplers, and polarizers, use the birefringence of the material to accomplish the desired characteristics. The directionally dependent permittivity by electro-optic, magneto-optic, and thermo-optic modulations can produce such birefringence in various crystals. The off-diagonal terms of the permittivity tensor are nonzero when the three mutually orthogonal principal axes of the crystals are not aligned with the waveguide coordinates (namely the x, y, and z axes). Therefore, developing a numerical scheme according to the anisotropic wave equations provides a vital tool for modeling diverse and practical photonic components. The commonly used numerical schemes can be divided into two categories. The first category is schemes based on the beam propagation method (BPM), such as the finite difference (FD) [1–3] and finite element (FE) [4–7]. The second category is schemes based on the eigenvalue algorithm, such as the FD [8–10], FE [11,12], and Fourier decomposition [13] methods. The mode field profiles offer a starting field for the BPM-based method; however, this method needs to consider extra convergence conditions such as the step size of propagation, the two-point recurrence scheme, and the reference refractive index. In addition, the real-axis BPM cannot determine the higher order modes, unless the real propagation axis is transformed to the imaginary axis as in the imaginary-distance BPM [14,15]. Consequently, for the eigenmode analysis, the eigenvalue-based algorithm is more appropriate and efficient and can solve all higher order modes simultaneously.

Recently, a pseudospectral method with high-order accuracy and fast convergence that was originally developed to solve computational fluid dynamics problems [16], has been introduced to the area of computational electromagnetics [17]. In order to analyze the guided modes in waveguides, I developed semivectorial [18] and full-vectorial [19] pseudospectral mode solvers based on the matrix eigenvalue equation. In previous studies [18,19], I used different basis functions to expand the mode field profiles in distinct subdomains (I used Chebyshev polynomials to expand interior subdomains and Laguerre–Gaussian functions (LGFs) to expand exterior subdomains). However, Chiang et al. [20] used a single kind of Chebyshev polynomials to expand the mode field profiles in distinct subdomains. My numerical analyses of guided modes showed that the pseudospectral scheme using different basis functions [19] is better than the one using single basis functions [20] because it requires less computational time and memory. This efficiency is achieved mainly because LGFs are well-matched to the exponential decay of guided modes. In addition, the pseudospectral scheme is applied to solving the leaky modes of an antiresonant reflecting optical waveguide [21] and a one-ring triangular holey fiber structure [22]. To achieve the same order of accuracy as the FD and FE schemes, the pseudospectral approach requires considerably less memory because it has fewer unknown variables [19,20].

However, these reports [16–20] only consider isotropic media. Although a recent report proposed an anisotropic waveguide mode solver based on the multidomain spectral collocation method [23], it used single Chebyshev polynomials to expand the mode fields in all subdomains; thus, the more efficient LGFs were not used to expand the mode fields in the exterior subdomains. This study extends the previous pseudospectral scheme [19] to investigate dielectric optical waveguides with transverse, nondiagonal anisotropy using hybrid basis functions. Accordingly, one of the principal axes of the anisotropic materials is thus constrained in the waveguide direction. Numerical examples include the circularly-polarized modes of a magneto-optical raised strip waveguide and the mode patterns of a nematic liquid-crystal (LC) waveguide under different orientations (twist angles) of the liquid crystal molecule (director). The remaining sections in this paper are organized as follows. Section 2 formulates the wave equations for the transverse anisotropic materials; Section 3 presents the numerical approach; Section 4 presents two examples to demonstrate the computational performances in convergence and accuracy of the proposed method; Section 5 draws the conclusions.

2. Mathematical formulations

Consider a monochromatic electromagnetic wave, which has a time dependence of exp(jωt), propagating along the z direction in an inhomogeneous medium with permittivity tensor [ε], and the vector wave equation based on the magnetic field vector H in the frequency domain is given by

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

where ω is the angular frequency and μ ₀ is the permeability in vacuum. This study assumes that one of the material principal axes is aligned in the waveguide direction; thus, the permittivity tensor [ε] yields the transverse, nondiagonal anisotropic form as

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

where ε ₀ is the permittivity in vacuum and [ε_r] is the relative permittivity tensor. For a mode problem, all electromagnetic fields are assumed to have a z dependence of exp(−jßz), where β = k ₀ n_eff, the wave number in vacuum is k ₀ = ω ² μ ₀ ε ₀, and n_eff is the effective refractive index. Using the divergence relation

\nabla \cdot H = 0

, the magnetic field component H_z in Eq. (1) can be eliminated to yield the full vector eigenvalue equation presented by the transverse magnetic field components H_x and H_y as follows:

[\begin{matrix} P_{x x} & P_{x y} \\ P_{y x} & P_{y y} \end{matrix}] [\begin{array}{l} H_{x} \\ H_{y} \end{array}] = β^{2} [\begin{matrix} η_{y y} & - η_{y x} \\ - η_{x y} & η_{x x} \end{matrix}] [\begin{array}{l} H_{x} \\ H_{y} \end{array}],

where the differential operators P_xx, P_xy, P_yx, and P_yy are defined as

P_{x x} H_{x} = η_{y y} \frac{\partial^{2} H_{x}}{\partial x^{2}} + η_{z z} \frac{\partial^{2} H_{x}}{\partial y^{2}} - η_{y x} \frac{\partial^{2} H_{x}}{\partial y \partial x} + k_{0}^{2} H_{x},

P_{x y} H_{y} = (η_{y y} - η_{z z}) \frac{\partial^{2} H_{y}}{\partial x \partial y} - η_{y x} \frac{\partial^{2} H_{y}}{\partial y^{2}} .

P_{y x} H_{x} = (η_{x x} - η_{z z}) \frac{\partial^{2} H_{x}}{\partial y \partial x} - η_{x y} \frac{\partial^{2} H_{x}}{\partial x^{2}}, and

P_{y y} H_{y} = η_{z z} \frac{\partial^{2} H_{y}}{\partial x^{2}} + η_{x x} \frac{\partial^{2} H_{y}}{\partial y^{2}} + η_{x y} \frac{\partial^{2} H_{y}}{\partial x \partial y} + k_{0}^{2} H_{y},

where

η_{x x} = \frac{ε_{y y} ε_{z z}}{Δ}, η_{x y} = η_{y x} = \frac{- ε_{x y} ε_{z z}}{Δ}, η_{y y} = \frac{ε_{x x} ε_{z z}}{Δ}, η_{z z} = \frac{ε_{x x} ε_{y y} - ε_{x y} ε_{y x}}{Δ}, Δ = ε_{z z} (ε_{x x} ε_{y y} - ε_{x y} ε_{y x}) .

The H_x and H_y components for an anisotropic material with nondiagonal terms of [ε], the coupling is determined from both the wave equations (P_xy and P_yx terms) and the physical interface conditions. In the proposed approach, the computational domain is first divided at the interfaces between different materials into several subdomains with uniform or continuous refractive index profiles, and the whole computational domain is then assembled from the subdomains while satisfying the physical interface conditions. These interface conditions include the continuous normal and tangential components of magnetic fields at each intra-element boundary, and the continuities of the longitudinal components H_z and E_z from the relations

\nabla \times H =

jωε ₀ n ²(x,y)E and

\nabla \cdot H = 0

, respectively. For a mode problem, all electromagnetic fields are assumed to have a z dependence of exp(−jßz), thus the continuity of the longitudinal component H_z yields

H_{z} = \frac{1}{j β} (\frac{\partial H_{x}}{\partial x} + \frac{\partial H_{y}}{\partial y}) .

and the continuity of E_z yields

E_{z} = \frac{j}{ω ε_{0} ε_{z z}} (\frac{\partial H_{x}}{\partial y} - \frac{\partial H_{y}}{\partial x}) .

In this case, only the horizontal and vertical interfaces are considered because the main aim of this study is to develop an anisotropic mode solver. For dealing with arbitrarily curved interfaces, Chiang et al. [20] introduced the curvilinear mapping technique to the pseudospectral scheme.

3. Numerical approach

In the pseudospectral approach [16,24], the computational domain is divided into several subdomains with uniform or continuous refractive index profiles. For an arbitrary interior rectangle subdomain r, it involves (n_x + 1) × (n_y + 1) collocation points in the corresponding (n_x + 1) × (n_y + 1) field unknowns as shown in Fig. 1(a) (where n_x = n_y = 10 is illustrated here).

Fig. 1 (a) The mesh division of an arbitrary interior subdomain r and (b) The mesh division of a problem with two subdomains labeled as 1 and 2.

Download Full Size | PDF

The transverse magnetic field components to be solved in this subdomain are thus expanded as a product of suitable basis functions comprising the Lagrange-type interpolation θ(x) in the x-direction, ψ(y) in the y-direction, and the H_x and H_y field values denoted as H_x,pq and H_y,pq at the (n_x + 1) × (n_y + 1) collocation points as follows:

H_{x}^{r} (x, y) = \sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} θ_{p}^{r} (x) ψ_{q}^{r} (y) H_{x, p q}^{r},

H_{y}^{r} (x, y) = \sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} θ_{p}^{r} (x) ψ_{q}^{r} (y) H_{y, p q}^{r},

where

θ_{p}^{r} (x_{m}) = δ_{m p}

,

ψ_{p}^{r} (y_{m}) = δ_{m p}

, and δ_mp denotes the Kronecker delta. Substituting the Eqs. (8a) and (8b) into Eq. (3), and then the pseudospectral scheme requires that the Eq. (3) has to be satisfied exactly at the specific (n_x + 1) × (n_y + 1) collocation points that depend on the selected basis functions. The wave equations are thus converted to a system of linear equations. Note that the boundary collocation points (i.e., (x ₀, y_j), (x ₁₀, y_j), (x_i, y ₀), and (x_i, y ₁₀) as shown in Fig. 1(a), where i = 0 to 10 and j = 0 to 10) of the rectangle subdomain r are shared by the adjacent subdomains, and the linear equations for the boundary collocation points have to be replaced by the physical interface conditions. These physical interface conditions are the continuous normal and tangential components of magnetic fields, and Eqs. (6) and (7). Removing the linear equations at these boundary collocation points, the matrix size is reduced to 2×(n_x −1)×(n_y − 1), and a matrix eigenvalue equation for the subdomain r is given by

[\begin{matrix} P_{x x}^{r} & P_{x y}^{r} \\ P_{y x}^{r} & P_{y y}^{r} \end{matrix}] [\begin{array}{l} H_{x}^{r} \\ H_{y}^{r} \end{array}] = β^{2} [\begin{matrix} η_{y y}^{r} & - η_{y x}^{r} \\ - η_{x y}^{r} & η_{x x}^{r} \end{matrix}] [\begin{array}{l} H_{x}^{r} \\ H_{y}^{r} \end{array}],

where the operators P_xx, P_xy, P_yx, and P_yy are

\begin{matrix} P_{x x}^{r} H_{x}^{r} = [η_{y y}^{r} \frac{\partial^{2} H_{x}^{r}}{\partial x^{2}} + η_{z z}^{r} \frac{\partial^{2} H_{x}^{r}}{\partial y^{2}} - η_{y x}^{r} \frac{\partial^{2} H_{x}^{r}}{\partial y \partial x} + k_{0}^{2} H_{x}^{r}] |_{x = x_{i}, y = y_{j}}, \\ = \sum_{i = 1}^{n_{x} - 1} \sum_{j = 1}^{n_{y} - 1} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {η_{y y}^{r} θ_{p}^{r (2)} (x) ψ_{q}^{r} (y) + η_{z z}^{r} θ_{p}^{r} (x) ψ_{q}^{r (2)} (y) - η_{y x}^{r} θ_{p}^{r (1)} (x) ψ_{q}^{r (1)} (y) ... \\ + k_{0}^{2} θ_{p}^{r} (x) ψ_{q}^{r} (y)}] |_{x = x_{i}, y = y_{j}} \cdot [H_{x, p q}^{r}], \end{matrix}

\begin{matrix} P_{x y}^{r} H_{y}^{r} = [(η_{y y}^{r} - η_{z z}^{r}) \frac{\partial^{2} H_{y}^{r}}{\partial x \partial y} - η_{y x}^{r} \frac{\partial^{2} H_{y}^{r}}{\partial y^{2}}] |_{x = x_{i}, y = y_{j}}, \\ = \sum_{i = 1}^{n_{x} - 1} \sum_{j = 1}^{n_{y} - 1} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {(η_{y y}^{r} - η_{z z}^{r}) θ_{p}^{r (1)} (x) ψ_{q}^{r (1)} (y) - η_{y x}^{r} θ_{p}^{r} (x) ψ_{q}^{r (2)} (y)}] |_{x = x_{i}, y = y_{j}} \cdot [H_{y, p q}^{r}], \end{matrix}

\begin{matrix} P_{y x}^{r} H_{x}^{r} = [(η_{x x}^{r} - η_{z z}^{r}) \frac{\partial^{2} H_{x}^{r}}{\partial x \partial y} - η_{x y}^{r} \frac{\partial^{2} H_{x}^{r}}{\partial x^{2}}] |_{x = x_{i}, y = y_{j}}, \\ = \sum_{i = 1}^{n_{x} - 1} \sum_{j = 1}^{n_{y} - 1} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {(η_{x x}^{r} - η_{z z}^{r}) θ_{p}^{r (1)} (x) ψ_{q}^{r (1)} (y) - η_{x y}^{r} θ_{p}^{r (2)} (x) ψ_{q}^{r} (y)}] |_{x = x_{i}, y = y_{j}} \cdot [H_{x, p q}^{r}], \end{matrix}

\begin{matrix} P_{y y}^{r} H_{y}^{r} = [η_{z z}^{r} \frac{\partial^{2} H_{y}^{r}}{\partial x^{2}} + η_{x x}^{r} \frac{\partial^{2} H_{y}^{r}}{\partial y^{2}} + η_{x y}^{r} \frac{\partial^{2} H_{y}^{r}}{\partial y \partial x} + k_{0}^{2} H_{y}^{r}] |_{x = x_{i}, y = y_{j}}, \\ = \sum_{i = 1}^{n_{x} - 1} \sum_{j = 1}^{n_{y} - 1} [\sum_{p = 0}^{n_{x}} \sum_{q = 0}^{n_{y}} {η_{z z}^{r} θ_{p}^{r (2)} (x) ψ_{q}^{r} (y) + η_{x x}^{r} θ_{p}^{r} (x) ψ_{q}^{r (2)} (y) + η_{x y}^{r} θ_{p}^{r (1)} (x) ψ_{q}^{r (1)} (y) ... \\ + k_{0}^{2} θ_{p}^{r} (x) ψ_{q}^{r} (y)}] |_{x = x_{i}, y = y_{j}} \cdot [H_{y, p q}^{r}], \end{matrix}

and

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

and

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

indicate the h-th order derivatives of θ_p(x) to x and ψ_q(y) to y, respectively. After obtaining the matrix eigenvalue equation for the subdomain r, a global matrix equation can be obtained by assembling all the matrix equations. Assuming there are total m subdomains in whole computational window, the pattern of the matrix elements forms the following block diagonal matrix equation:

[\begin{matrix} Q_{1} & 0 & 0 & 0 \\ 0 & Q_{2} & 0 & 0 \\ 0 & 0 & ⋱ & 0 \\ 0 & 0 & 0 & Q_{m} \end{matrix}] [\begin{matrix} H_{1} \\ H_{2} \\ ⋮ \\ H_{m} \end{matrix}] = β^{2} [\begin{matrix} η_{1} & 0 & 0 & 0 \\ 0 & η_{2} & 0 & 0 \\ 0 & 0 & ⋱ & 0 \\ 0 & 0 & 0 & η_{m} \end{matrix}] [\begin{matrix} H_{1} \\ H_{2} \\ ⋮ \\ H_{m} \end{matrix}],

where

Q_{r} = [\begin{matrix} P_{x x}^{r} & P_{x y}^{r} \\ P_{y x}^{r} & P_{y y}^{r} \end{matrix}], Η_{r} = [\begin{matrix} H_{x}^{r} \\ H_{y}^{r} \end{matrix}], η_{r} = [\begin{matrix} η_{y y}^{r} & - η_{y x}^{r} \\ - η_{x y}^{r} & η_{x x}^{r} \end{matrix}], (r = 1, 2, 3... m) .

However, the final matrix to be solved must include other constrains at interface boundary points which satisfy the physical interface conditions. In Eq. (10a), one observes that the element P_xxH_x includes the field unknowns H_x _, _pq at interface boundary points (where the subscripts of summations are from p = 0 to n_x and q = 0 to n_y), although the linear equations satisfied at interface boundary points have been removed (where the subscripts of summations are only from i = 1 to n_x-1 and j = 1 to n_y-1). Therefore, these interface field unknowns ({H_x _,0,0, H_x _,0,1,… H_x _,0, _ny},{H_x _, _nx _,0, H_x _, _nx _,1,… H_x _, _nx _, _ny},{H_x _,0,0, H_x _,1,0,… H_x _, _nx _,0},{H_x _,0, _ny, H_x _,1, _ny _,… H_x _, _nx _, _ny}) at interface boundary points involved in Eq. (11) have to be replaced by the field unknowns at all collocation points other than the interface boundary points by using the physical boundary conditions. The similar results are also found in other elements P_xyH_y, P_yxH_x, and P_yyH_y.

For clearly, a problem with only two subdomains shown in Fig. 1(b), labeled as 1 and 2, is adopted to illustrate the coupling between subdomains. The interface boundary points ( $(x_{n_{x}}^{(1)}, y_{i}) = (x_{0}^{(2)}, y_{i}),$ where i = 0, 1, 2,…n_y, labeled as the red circles) located in the vertical line are shared by the two subdomains. By the physical interface conditions, the continuities of the normal and tangential components of the magnetic field vector have to be implicitly satisfied, and the continuity of H_z in Eq. (6) gives

[\begin{matrix} R_{^{x}}^{(1)} & R_{^{y}}^{(1)} \end{matrix}] [\begin{matrix} H_{x}^{(1)} \\ H_{y}^{(1)} \end{matrix}] = [\begin{matrix} R_{^{x}}^{(2)} & R_{^{y}}^{(2)} \end{matrix}] [\begin{matrix} H_{x}^{(2)} \\ H_{y}^{(2)} \end{matrix}],

where

R_{s}^{(1)} = [\begin{matrix} R_{s}^{(1)} (ψ_{0} (y_{0}^{(1)})) & R_{s}^{(1)} (ψ_{1} (y_{0}^{(1)})) & \begin{matrix} \begin{matrix} . & . & . \end{matrix} \end{matrix} & R_{s}^{(1)} (ψ_{n_{y}} (y_{0}^{(1)})) \\ R_{s}^{(1)} (ψ_{0} (y_{1}^{(1)})) & R_{s}^{(1)} (ψ_{1} (y_{1}^{(1)})) & \begin{matrix} . & . & . \end{matrix} & R_{s}^{(1)} (ψ_{n_{y}} (y_{1}^{(1)})) \\ \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . & . & . \\ . & . & . \\ . & . & . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} \\ R_{s}^{(1)} (ψ_{0} (y_{n_{y}}^{(1)})) & R_{s}^{(1)} (ψ_{1} (y_{n_{y}}^{(1)})) & \begin{matrix} . & . & . \end{matrix} & R_{s}^{(1)} (ψ_{n_{y}} (y_{n_{y}}^{(1)})) \end{matrix}], s = x, y,

\begin{array}{l} R_{x}^{(1)} (ψ_{i} (y_{0}^{(1)})) = [\begin{matrix} θ_{0}^{(1)} (x_{_{n_{x}}}^{(1)}) ψ_{i} (y_{0}^{(1)}) & θ_{1}^{(1)} (x_{_{n_{x}}}^{(1)}) ψ_{i} (y_{0}^{(1)}) & . & . & . & θ_{n_{x}}^{(1)} (x_{_{n_{x}}}^{(1)}) ψ_{i} (y_{0}^{(1)}) \end{matrix}], \\ R_{y}^{(1)} (ψ_{i} (y_{0}^{(1)})) = [\begin{matrix} θ_{0}^{} (x_{n_{x}}^{(1)}) ψ_{i}^{(1)} (y_{0}^{(1)}) & θ_{1}^{} (x_{n_{x}}^{(1)}) ψ_{i}^{(1)} (y_{0}^{(1)}) & . & . & . & θ_{n_{x}} (x_{_{n_{x}}}^{(1)}) ψ_{i}^{(1)} (y_{0}^{(1)}) \end{matrix}], \\ R_{s}^{(2)} = [\begin{matrix} R_{s}^{(2)} (ψ_{0} (y_{0}^{(2)})) & R_{s}^{(2)} (ψ_{1} (y_{0}^{(2)})) & \begin{matrix} \begin{matrix} . & . & . \end{matrix} \end{matrix} & R_{s}^{(2)} (ψ_{n_{y}} (y_{0}^{(2)})) \\ R_{s}^{(2)} (ψ_{0} (y_{1}^{(2)})) & R_{s}^{(2)} (ψ_{1} (y_{1}^{(2)})) & \begin{matrix} . & . & . \end{matrix} & R_{s}^{(2)} (ψ_{n_{y}} (y_{1}^{(2)})) \\ \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . & . & . \\ . & . & . \\ . & . & . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} \\ R_{s}^{(2)} (ψ_{0} (y_{n_{y}}^{(2)})) & R_{s}^{(2)} (ψ_{1} (y_{n_{y}}^{(2)})) & \begin{matrix} . & . & . \end{matrix} & R_{s}^{(2)} (ψ_{n_{y}} (y_{n_{y}}^{(2)})) \end{matrix}], s= x, y, \\ R_{x}^{(2)} (ψ_{i} (y_{0}^{(2)})) = [\begin{matrix} θ_{0}^{(1)} (x_{_{0}}^{(2)}) ψ_{i} (y_{0}^{(2)}) & θ_{1}^{(1)} (x_{_{0}}^{(2)}) ψ_{i} (y_{0}^{(2)}) & . & . & . & θ_{n_{x}}^{(1)} (x_{_{0}}^{(2)}) ψ_{i} (y_{0}^{(2)}) \end{matrix}], \\ R_{y}^{(2)} (ψ_{i} (y_{0}^{(2)})) = [\begin{matrix} θ_{0}^{} (x_{0}^{(2)}) ψ_{i}^{(1)} (y_{0}^{(2)}) & θ_{1}^{} (x_{0}^{(2)}) ψ_{i}^{(1)} (y_{0}^{(2)}) & . & . & . & θ_{n_{x}}^{} (x_{_{0}}^{(2)}) ψ_{i}^{(1)} (y_{0}^{(2)}) \end{matrix}], \end{array}

and the continuity of E_z in Eq. (7) gives

[\begin{matrix} U_{^{x}}^{(1)} & U_{^{y}}^{(1)} \end{matrix}] [\begin{matrix} H_{x}^{(1)} \\ H_{y}^{(1)} \end{matrix}] = [\begin{matrix} U_{^{x}}^{(2)} & U_{^{y}}^{(2)} \end{matrix}] [\begin{matrix} H_{x}^{(2)} \\ H_{y}^{(2)} \end{matrix}],

where

U_{s}^{(1)} = [\begin{matrix} U_{s}^{(1)} (ψ_{0} (y_{0}^{(1)})) & U_{s}^{(1)} (ψ_{1} (y_{0}^{(1)})) & \begin{matrix} \begin{matrix} . & . & . \end{matrix} \end{matrix} & U_{s}^{(1)} (ψ_{n_{y}} (y_{0}^{(1)})) \\ U_{s}^{(1)} (ψ_{0} (y_{1}^{(1)})) & U_{s}^{(1)} (ψ_{1} (y_{1}^{(1)})) & \begin{matrix} . & . & . \end{matrix} & U_{s}^{(1)} (ψ_{n_{y}} (y_{1}^{(1)})) \\ \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . & . & . \\ . & . & . \\ . & . & . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} \\ U_{s}^{(1)} (ψ_{0} (y_{n_{y}}^{(1)})) & U_{s}^{(1)} (ψ_{1} (y_{n_{y}}^{(1)})) & \begin{matrix} . & . & . \end{matrix} & U_{s}^{(1)} (ψ_{n_{y}} (y_{n_{y}}^{(1)})) \end{matrix}], s = x, y,

\begin{array}{l} U_{x}^{(1)} (ψ_{i} (y_{0}^{(1)})) = \frac{1}{ε_{_{z z}}^{(1)}} [\begin{matrix} θ_{0} (x_{_{n_{x}}}^{(1)}) ψ_{_{i}}^{(1)} (y_{0}^{(1)}) & θ_{1} (x_{_{n_{x}}}^{(1)}) ψ_{_{i}}^{(1)} (y_{0}^{(1)}) & . & . & . & θ_{n_{x}} (x_{_{n_{x}}}^{(1)}) ψ_{_{i}}^{(1)} (y_{0}^{(1)}) \end{matrix}], \\ U_{y}^{(1)} (ψ_{i} (y_{0}^{(1)})) = \frac{- 1}{ε_{_{z z}}^{(1)}} [\begin{matrix} θ_{0}^{(1)} (x_{n_{x}}^{(1)}) ψ_{i} (y_{0}^{(1)}) & θ_{1}^{(1)} (x_{n_{x}}^{(1)}) ψ_{i} (y_{0}^{(1)}) & . & . & . & θ_{n_{x}}^{(1)} (x_{_{n_{x}}}^{(1)}) ψ_{i} (y_{0}^{(1)}) \end{matrix}], \\ U_{s}^{(2)} = [\begin{matrix} U_{s}^{(2)} (ψ_{0} (y_{0}^{(2)})) & U_{s}^{(2)} (ψ_{1} (y_{0}^{(2)})) & \begin{matrix} \begin{matrix} . & . & . \end{matrix} \end{matrix} & U_{s}^{(2)} (ψ_{n_{y}} (y_{0}^{(2)})) \\ U_{s}^{(2)} (ψ_{0} (y_{1}^{(2)})) & U_{s}^{(2)} (ψ_{1} (y_{1}^{(2)})) & \begin{matrix} . & . & . \end{matrix} & U_{s}^{(2)} (ψ_{n_{y}} (y_{1}^{(2)})) \\ \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} & \begin{matrix} . & . & . \\ . & . & . \\ . & . & . \end{matrix} & \begin{matrix} . \\ . \\ . \end{matrix} \\ U_{s}^{(2)} (ψ_{0} (y_{n_{y}}^{(2)})) & U_{s}^{(2)} (ψ_{1} (y_{n_{y}}^{(2)})) & \begin{matrix} . & . & . \end{matrix} & U_{s}^{(2)} (ψ_{n_{y}} (y_{n_{y}}^{(2)})) \end{matrix}], s = x, y, \\ U_{x}^{(2)} (ψ_{i} (y_{0}^{(2)})) = \frac{1}{ε_{_{z z}}^{(2)}} [\begin{matrix} θ_{0} (x_{_{0}}^{(2)}) ψ_{_{i}}^{(1)} (y_{0}^{(2)}) & θ_{1} (x_{_{0}}^{(2)}) ψ_{_{i}}^{(1)} (y_{0}^{(2)}) & . & . & . & θ_{n_{x}} (x_{_{0}}^{(2)}) ψ_{_{i}}^{(1)} (y_{0}^{(2)}) \end{matrix}], \\ U_{y}^{(2)} (ψ_{i} (y_{0}^{(2)})) = \frac{- 1}{ε_{_{z z}}^{(2)}} [\begin{matrix} θ_{0}^{(1)} (x_{0}^{(2)}) ψ_{i} (y_{0}^{(2)}) & θ_{1}^{(1)} (x_{0}^{(2)}) ψ_{i} (y_{0}^{(2)}) & . & . & . & θ_{n_{x}}^{(1)} (x_{_{0}}^{(2)}) ψ_{i} (y_{0}^{(2)}) \end{matrix}], \end{array}

and i = 0, 1, 2,…n_y. In Eqs. (13) and (15), the

H_{x}^{(1)}

and

H_{y}^{(1)}

denote the magnetic field components of H _x and H _y, respectively, at the collocation point including the boundary points in subdomain 1. Using Eqs. (13) and (15), the field unknowns at interface boundary points are represented by the interior collocation points inside the subdomains 1 and 2. Extending to all subdomains in computational window, the final matrix includes only the collocation points except the interface boundary points. The final matrix pattern to be solved is no longer a block diagonal form but turns to a matrix with approximately 38% sparsity. In this study, the final matrix equation is solved using the Arnoldi iteration method [25], which efficiently solves generalized eigenvalue problems. Once the field values are obtained by solving the final matrix, the field values at interface boundary points are then found through Eqs. (13) and (15).

We now determine the appropriate basis functions by expanding the H_x and H_y components. The reports [18,19] show that the Chebyshev polynomials are better for expanding the optical fields in interior subdomains because of their mathematical robustness to non-periodic structures. The other choice to expand the field profiles in the interior regions is Legendre polynomials. The comparison of the convergence between Chebyshev and Legendre basis functions for analyzing a circular waveguide was made by Chiang et al. [20]. They calculated the relative errors in effective index relative to the exact values of HE₁₁ mode. The results showed that Chebyshev polynomials achieve the better accuracy and faster convergence than that based on Legendre polynomials. This result is resulted from the different mesh division profiles (i.e., collocation point profiles) between Chebyshev and Legendre polynomials [16,24]. In contrast to the exterior subdomains, the LGFs with exponential decay characteristics match the guided field profiles with exponential decay. In general, there are three ways to expand the unknown fields in a semi-infinite interval [0,∞): (1) expand in LGFs, (2) map the semi-infinite extent into a finite one and then expand in Chebyshev polynomials, and (3) truncate the semi-infinite domain to finite one [0, x _max] and use a Chebyshev expansion [16]. For the third choice, only one way by increasing x _max, as the number of terms of basis functions can be used to achieve the exponential convergence. Therefore, the computational efforts in CPU time and memory storage are much heavy. The second approach has several types of map including the algebraic, logarithmic, and exponential maps to be considered. In addition, these maps also need to efficiently determine their optimal scale factors. In contrast, LGFs with exponential decay characteristics (in mathematical) efficiently match the exponential field profiles of guided modes (in physical) on waveguide problems, and the scaling factor in LGFs has been definitely determined by the report for the isotropic media [19]. Besides, Boyd [24] has noticed that uniform spectral accuracy of LGFs requires at least exponential decay of the function as x→∞, the feature is similar to the characteristic of the general guided field profiles. If Chebyshev polynomials are used to represent the guided mode profiles, more terms of basis functions than that of LGFs are required. In mathematics, various basis functions can be converted to their Lagrange-type form. For Chebyshev polynomials, the explicit form of the Lagrange-type interpolation function is as follows [24]:

θ_{p} (x) = \frac{{(- 1)}^{p + 1} (1 - x^{2}) T_{v}^{'} (x)}{c_{p} n^{2} (x - x_{p})}, c_{p} = {\begin{cases} 2, if p =0, N \\ 1,  if 1 \leq p \leq N -1 \end{cases}

where T_v(x) is the general Chebyshev polynomial of order v and x_p denotes the collocation points for the Chebyshev polynomials. For the LGFs, the explicit form is as follows [24]:

θ_{p} (α x) = \frac{x L_{v} (α x)}{α (x - x_{p}) [x L_{v}^{'} (α x)] |_{x = x_{p}}} e^{- α (x - x_{p}) / 2},

where L_v(αx) is the Laguerre polynomial of order v and x_p denotes the collocation points for the LGFs. In Eq. (18), the scaling factor α significantly affects the numerical accuracy of results obtained with the LGFs. The definite value of α for isotropic materials has been derived in previous studies [19]; hence, results obtained with LGFs are sufficiently accurate.

For a given number of terms of basis functions, α is determined by combining two techniques. One is from the derivation of Tang [26] and the other is from the idea of effective index method (EIM) [27]. Firsyt, the optimum α is obtained through the following formula [26]

α = \max_{0 \leq i \leq n} \frac{{x_{i}}}{M}

where M is the spreading of the guided fields and x_i is the collocation point at position i. The choice of α is thus converted to determine the value of M. Using EIM, the effective index and the decay rate of the desired guided mode has to be obtained in advance, and then M can be estimated subsequently. For finding out the effective index and the decay rate, a rib waveguide shown in Fig. 2(a) is used to describe the determination of M, where, [ε _c], [ε _s], and [ε _a] denote the permittivity tensors of core, substrate, and cover, respectively.

Fig. 2 (a) Cross section of a rib waveguide. (b) Equivalent slab waveguide to original rib waveguide after using EIM along the y direction.

Download Full Size | PDF

The first process is to divide the rib waveguide by the vertical dotted lines into three subregions labeled as 1, 2, and 3, and then the dispersion relation [27] of three-layer waveguide along the y-direction in each subregion computes the effective indices of subregions 1, 2 and 3 as n_y ₁, n_y ₂ and n_y ₃, respectively, as shown in Fig. 2(b). Importantly, different from the condition of isotropic media [19], the refractive index for each material to be used in the dispersion relation [27] is (ε_xx)^1/2 in this study. The dispersion relation is utilized again in the x-direction to obtain the equivalent effective index n_eq to the original waveguide structure. The decay rates γ along the y-direction extending to the substrate for the three subregions are given by

γ_{y s i} = k_{0} \sqrt{n_{y i}^{2} - ε_{x x, s}}, i =1, 2, 3,

and those along the y-direction extending to air are

γ_{y a i} = k_{0} \sqrt{n_{y i}^{2} - ε_{x x, a}}, i =1, 2, 3,

and the decay rates along the x-direction are

γ_{x i} = k_{0} \sqrt{n_{e q}^{2} - n_{y i}^{2}}, i =1, 3,

where n_eq is the effective refractive index for the original waveguide. Accordingly, the spreadings of the magnetic fields along the y-direction extending to the substrate M_ysi are given

\frac{H (y_{t} + M_{y s i})}{H (y_{t})} = \exp [- \int_{y_{t}}^{y_{t} + M_{y s i}} γ_{y s i} d y] = κ, i =1, 2, 3,

those along the y-direction extending to the air M_yai are

\frac{H (y_{t} + M_{y a i})}{H (y_{t})} = \exp [- \int_{y_{t}}^{y_{t} + M_{y a i}} γ_{y a i} d y] = κ, i =1, 2, 3,

and the spreadings along the x-direction M_xi are

\frac{H (x_{t} + M_{x i})}{H (x_{t})} = \exp [- \int_{x_{t}}^{{\begin{matrix} x \end{matrix}}_{t} + M_{x i}} γ_{x i} d x] = κ, i =1, 3,

where x_t and y_t are the turning points in x- and y-directions, respectively. In a step-index waveguide, the turning points are located at the interfaces between different materials. As shown in Eqs. (21a)–(21c), the parameter κ denotes the ratio of the field amplitude at the positions of turning points plus the spreadings M relative to that at turning points. For simplicity, we let κ be equal in both the x- and the y-directions. In the report [19], the ratio of κ = 10⁻² is sufficient to involve the significant contribution of the guided field profiles.

4. Simulation results and discussion

This section demonstrates the capabilities of the proposed scheme for simulating anisotropic waveguides through a magneto-optical raised strip waveguide and a nematic LC optical waveguide. The calculated results of the magneto-optical raised strip waveguide are compared with those obtained from the rigorous vector FD method [8], and the convergences of the two schemes are also discussed. For the second example, the mode characteristics of the nematic LC optical waveguide versus the twist angle of the director are presented in detail.

4.1 A magneto-optical raised strip waveguide

A polarization Faraday rotator based on the phase matching of the first and the second order modes can be achieved by applying longitudinally directed magnetization in a magneto-optic channel waveguide [28]. Figure 3(a) shows a schematic diagram of the lossless magneto-optical raised strip waveguide, which consists of a core region filled with a yttrium iron garnet (YIG) over an isotropic gadolinium gallium garnet (GGG) substrate with a refractive index n_s = 1.95, and the cover is air with n_a = 1.

Fig. 3 (a) Schematic showing cross section of magneto-optical raised strip waveguide with permittivity tensor [ε] in the core region and the refractive indices of GGG substrate n_s and air n_a, and (b) division of computational domain for magneto-optical raised strip waveguide.

Download Full Size | PDF

Ignoring the differences of the diagonal components of the x, y, and z positions, which result from the growth- and stress-induced anisotropies of growing a YIG-film on a (111)-oriented GGG-substrate, the permittivity tensor of YIG applied by a static magnetic field along the propagation direction (z-axis) is given as [28]:

[ε] = ε_{0} [\begin{matrix} n_{f}^{2} & j ς & 0 \\ - j ς & n_{f}^{2} & 0 \\ 0 & 0 & n_{f}^{2} \end{matrix}]

where n_f = 2.302 and ζ is the first-order magneto-optic effect and proportional to static magnetization. The width w of the core is 800 nm, and its thickness t is set to be 607.6 nm, which was carefully adjusted to achieve the nearly degenerate mode between the first and the second order modes when ζ = 0 [8,28], for the operating wavelength λ = 1300 nm in free space.

The proposed scheme divides the computational window into nine subdomains as illustrated in Fig. 3(b). The mode profiles in the interior subdomains were expanded with Chebyshev polynomials, and those in the exterior subdomains were expanded with LGFs with the determined scaling factor α. For example, the field profiles of the subdomains 1, 3, 7, and 9 with semi-infinite extents were expanded with LGFs in both the x- and y-directions. For subdomains 2 and 8, Chebyshev polynomials expanded the field profiles in the x-direction and LGFs expanded them in the y-direction. For subdomain 5 with finite extents both in the x- and the y-directions, a single kind of Chebyshev polynomial expanded them. Because of the well-matched profiles of the guided-mode profile and the mathematical characteristic of LGFs, the proposed scheme requires no extra effort for the computational boundary conditions. In [23], the field profiles in the exterior subdomains were expanded with Chebyshev polynomials incorporating the perfectly matched layer (PML) condition; however, there are two disadvantages to this approach. The first is that the zero boundary conditions must be enforced; this thus significantly enlarges the computational window and thus consumes extra computational efforts. The other one is that although adding a PML to absorb the field profiles can reduce the computational window, such efforts are not required for solving the guided modes unless the leaky modes are considered. However, the report [21] has demonstrated that the desired accuracy can be achieved by using the simple second-order Mur’s absorbing boundary condition.

For examining the convergence of the proposed scheme, I calculated the results for the condition ζ = 0.005. Note that the terms of the basis functions for expanding the exterior subdomains were fixed at the number N_ex = 10. Table 1 shows the convergence of effective refractive indices of the first (labeled as $n_{e f f}^{(1)}$ ) and the second (labeled as $n_{e f f}^{(2)}$ ) order modes and the computational time versus the term of the basis functions N_int expanding the interior subdomains. Throughout this study, equal terms of the basis function are used in both the x- and y-directions.The calculations were performed on a personal computer with an Intel core i7 CPU 920 at 2.67 GHz. The number of grid points in this example was (N_int + 2N_ex − 2) × (N_int + 2N_ex− 2), and the computational window was approximately 3.2 μm × 3 μm. The difference of the calculated effective index between N_int = 10 and N_int = 25 was in the order of 10⁻⁴; by increasing the terms of basis functions to N_int = 20, this difference was reduced to the order of 10⁻⁵. Thus the proposed scheme obtains a sufficiently accurate result (it achieved accuracy of the order of 10⁻⁴ using only 28 × 28 grid points). The growth of the computational time listed in Table 1 was fairly slow due to the use of LGFs for exterior subdomains. If Chebyshev polynomials are used to expand the field profiles for all subdomains, the computational times were 9.9 s, 44.8 s, and 281.9s while using 15, 20, and 25 terms, respectively. Therefore, the proposed scheme is more efficient than that using the unified set of basis functions all over the whole domain. For different values of ζ, Table 2 shows the results obtained by the proposed scheme using N_int = 25 and N_ex = 10 (the grid points were now 45 × 45) and those obtained with the vector FD method [8]. The results obtained from the two schemes are in good agreement. In the vector FD method [8], the grid points were 320 × 348 over a computational window of 3.2 μm × 2.9 μm. This comparison demonstrates that the computational efficiency of the proposed method is superior to that of the vector FD method. For the condition ζ = 0.005, Fig. 4 shows the |H_x| and |H_y| field profiles of the first order mode, and Fig. 5 shows the |H_x| and |H_y| field profiles of the second order mode. The magnitudes and profiles of the |H_x| and |H_y| fields are nearly equal for the two modes.

Table 1. Effective refractive indices of the first and the second order modes and the computational time versus the term of the basis functions for expanding the interior subdomains N_int, while ζ = 0.005

View Table | View all tables in this article

Table 2. The calculated effective indices obtained by this study and the vector FD method [8] for different values of ζ

View Table | View all tables in this article

Fig. 4 Mode profiles of (a) |H_x| and (b) |H_y| of the first order mode for ζ = 0.005.

Download Full Size | PDF

Fig. 5 Mode profiles of (a) |H_x| and (b) |H_y| of the second order mode for ζ = 0.005.

Download Full Size | PDF

To clearly represent the polarization states of the two modes, Figs. 6 and 7 plot the field profiles of |H _± | = |H_x ± jH_y| /2of the first and the second order modes, respectively. For the first order mode, the power of the |H_x − jH_y| /2field is suppressed by approximately 17 dB compared with the peak value of the |H_x + jH_y| /2field, a suppression amount equal to that calculated by [8]. Figure 6 indicates that the phase difference between the complex H_x and H_y components was approximately π/2, and the first order mode is in the left-hand circularly-polarized state. In contrast, Fig. 7 indicates that the second order mode is in the right-hand circularly-polarized state. In addition, the calculated Faraday rotation rate for the condition ζ = 0.005 is $Θ_{F} = (π / λ) (n_{e f f}^{(1)} - n_{e f f}^{(2)}) = 2397^{o} / cm$ .

Fig. 6 Mode profiles of (a) |H_x − jH_y |/2and (b) |H_x + jH_y |/2of the first order mode for ζ = 0.005.

Download Full Size | PDF

Fig. 7 Mode profiles of (a) |H_x − jH_y|/2and (b) |H_x + jH_y|/2of the second order mode for ζ = 0.005.

Download Full Size | PDF

4.2 A nematic liquid-crystal channel waveguide

The following points should be considered when using LCs to design integrated optics components: the low scattering loss near the infrared wavelength, the high electro-optic effect, the low power requirement, and the low manufacturing cost [29]. D’Alessandro et al. [30] experimentally demonstrated the feasibility of an integrated optics polarizer, which comprises LCs, with extinction ratios of 25 dB. Therefore, many photonic devices comprising LCs have been manufactured. Recently, some numerical schemes, such as FD-BPM [3], FDFD [9,10], FE [31], and FE-BPM [32], were developed to analyze the propagation characteristics of LC waveguides. Figure 8(a) shows a schematic diagram of a nematic LC channel optical waveguide considered in this study. A core of width w = 3 μm and height t = 3 μm filled with LCs (5CB) [33] is embedded in a glass substrate with a refractive index n_s = 1.45, and the cover is air.

Fig. 8 (a) Schematic showing cross section of the nematic LC channel optical waveguide with permittivity tensor of the core region [ε] and the refractive indices of glass substrate n_s and air n_a, (b) Schematic showing the twist angle φ of LC molecule and orientation of director $\hat{n}$ .

Download Full Size | PDF

The permittivity tensor of the core region is given by

[ε] = ε_{0} [\begin{matrix} ε_{x x} & ε_{x y} & 0 \\ ε_{y x} & ε_{y y} & 0 \\ 0 & 0 & ε_{z z} \end{matrix}] = ε_{0} [\begin{matrix} n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \cos^{2} φ & (n_{e}^{2} - n_{o}^{2}) \cos φ \sin φ & 0 \\ (n_{e}^{2} - n_{o}^{2}) \cos φ \sin φ & n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \sin^{2} φ & 0 \\ 0 & 0 & n_{o}^{2} \end{matrix}],

where n_o = 1.5292 and n_e = 1.7072 are the ordinary and extraordinary refractive indices of the nematic LC, respectively, at the wavelength λ = 1.55 μm. In Fig. 8(b), φ denotes the twist angle of the LC between the x-axis and director

\hat{n}

. The voltage applied by the designed electrodes is applied to alter the value of φ. By the proposed scheme, the division of the computational window is the same as the first example shown in Fig. 3(b), and all the calculations used N_ex = 10 and N_int = 20 in each subdomain. Table 3 shows the calculated effective indices at various values of φ for the first seven modes.For the mode 1, the variation of effective index versus φ is insignificant. The result emerges in the major field patterns as shown in Fig. 9 . The small differences of effective index only result in the different minor field patterns.

Table 3. Calculated effective indices at φ = 0°, 30°, 45°, 60°, and 90° for the first seven modes using N_ex = 10 and N_int = 20 in each subdomain

View Table | View all tables in this article

Fig. 9 |H_x| and |H_y| mode patterns of the mode 1 for four twist angles φ.

Download Full Size | PDF

However, the variations of effective index versus φ are more significant for the higher order modes as shown in Fig. 10 . The symbols |H_m,x| and |H_m,y| denote the modului of the H_x and H_y components, respectively, of the m-th order guided mode. In Fig. 10, the symmetric field patterns corresponding to the y-axis are observed for the conditions φ = 0° and φ = 90° because the elements of the permittivity tensor only appear in the principal axes. For other values of φ, more complicated symmetries exist only in the lower order modes. For the conditions φ = 0° and φ = 90°, both of the minor powers corresponding to the peak values are suppressed by more than 40 dB, and those for φ = 30° and φ = 60° are suppressed by approximately 5 dB.

Fig. 10 |H_x| and |H_y| mode patterns of the mode 2 to mode 5 for four twist angles φ.

Download Full Size | PDF

In contrast, the field modului of the H_x and H_y components are nearly identical for φ = 45° because the two polarizations experience identical propagation in the core region, except for the small asymmetries in the x- and y- directions. Hence, the calculated field patterns of |H_x| and |H_y| are not perfectly identical. For brevity, the patterns for φ = 45° are not shown. In particular, the magnitudes of the major |H_x| and minor |H_y| fields are reversed for the seventh mode (not shown here) at these values of φ. This reverse phenomenon appears only in some of the higher order modes.

5. Conclusion

This study has presented a vector anisotropic mode solver based on a pseudospectral scheme with hybrid basis functions. The field profiles in the interior subdomains were expanded with Chebyshev polynomials and those in the exterior subdomains were expanded with Laguerre-Gaussian functions with an accurately determined scaling factor. The proposed scheme has considered the dielectric waveguides with transverse anisotropy, and thus one of the principal axes of materials is limited to point in the propagation direction of the waveguides. The simulations of the circularly-polarized modes in a magneto-optical raised strip waveguide with the method proposed in this paper have demonstrated the higher efficiency of the proposed method compared with the vector finite difference approach. In addition, this study has presented the detailed mode patterns in a nematic liquid-crystal channel waveguide at various twist angles of the liquid-crystal molecule. The proposed method offers an efficient modeling tool for designing functional photonic devices of considerable interest. Further analysis of altering the tilt angle of the liquid-crystal molecule, which has nine nonzero elements in the permittivity tensor, will be reported in oncoming day.

Acknowledgements

I 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. C. L. Xu, W. P. Huang, J. Chrostowski, and S. K. Chaudhuri, “A full-vectorial beam propagation method for anisotropic waveguides,” J. Lightwave Technol. 12(11), 1926–1931 (1994). [CrossRef]

2. P. Lüsse, K. Ramm, and H. G. Unger, “Vectorial eigenmode calculation for anisotropic planar optical waveguides,” Electron. Lett. 32(1), 38–39 (1996). [CrossRef]

3. Q. Wang, G. Farrell, and Y. Semenova, “Modeling liquid-crystal devices with the three-dimensional full-vector beam propagation method,” J. Opt. Soc. Am. A 23(8), 2014–2019 (2006). [CrossRef]

4. Y. Tsuji, M. Koshiba, and N. Takimoto, “Finite element beam propagation method for anisotropic optical waveguides,” J. Lightwave Technol. 17(4), 723–728 (1999). [CrossRef]

5. S. Selleri, L. Vincetti, and M. Zoboli, “Full-vector finite-element beam propagation method for anisotropic optical device analysis,” IEEE J. Quantum Electron. 36(12), 1392–1401 (2000). [CrossRef]

6. K. Saitoh and M. Koshiba, “Approximate scalar finite element beam-propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19(5), 786–792 (2001). [CrossRef]

7. J. P. da Silva, H. E. Hernandez-Figueroa, and A. M. F. Frasson, “Improved vectorial finite-element BPM analysis for transverse anisotropic media,” J. Lightwave Technol. 21(2), 567–576 (2003). [CrossRef]

8. A. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector finite difference modesolver for anisotropic dielectric waveguides,” J. Lightwave Technol. 26(11), 1423–1431 (2008). [CrossRef]

9. M. F. O. Hameed, S. S. A. Obayya, K. Al-Begain, M. I. Abo el Maaty, and A. M. Nasr, “Modal properties of an index guiding nematic liquid crystal based photonic crystal fiber,” J. Lightwave Technol. 27(21), 4754–4762 (2009). [CrossRef]

10. M. Y. Chen, S. M. Hsu, and H. C. Chang, “A finite-difference frequency-domain method for full-vectorial mode solutions of anisotropic optical waveguides with arbitrary permittivity tensor,” Opt. Express 17(8), 5965–5979 (2009). [CrossRef] [PubMed]

11. V. Schulz, “Adjoint high-order vectorial finite elements for nonsymmetric transversally anisotropic waveguides,” IEEE Trans. Microw. Theory Tech. 51(4), 1086–1095 (2003). [CrossRef]

12. B. G. Ward, “Finite element analysis of photonic crystal rods with inhomogeneous anisotropic refractive index tensor,” IEEE J. Quantum Electron. 44(2), 150–156 (2008). [CrossRef]

13. R. Pashaie, “Fourier decomposition analysis of anisotropic inhomogeneous dielectric waveguide structures,” IEEE Trans. Microw. Theory Tech. 55(8), 1689–1696 (2007). [CrossRef]

14. J. C. Chen and S. Jüngling, “Computation of high-order waveguide modes by imaginary-distance beam propagation method,” Opt. Quantum Electron. 26(3), S199–S205 (1994). [CrossRef]

15. T. Ando, H. Nakayama, S. Numata, J. Yamauchi, and H. Nakano, “Eigenmode analysis of optical waveguides by a Yee-mesh-based imaginary-distance propagation method for an arbitrary dielectric interface,” J. Lightwave Technol. 20(8), 1627–1634 (2002). [CrossRef]

16. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics (Springer-Verlag, New York, 1988).

17. Q. H. Liu, “A pseudospectral frequency-domain (PSFD) method for computational electromagnetics,” IEEE Antennas Wirel. Propag. Lett. 1(6), 131–134 (2002). [CrossRef]

18. C. C. Huang and C. C. Huang, “An efficient and accurate semivectorial spectral collocation method for analyzing polarized modes of rib waveguides,” J. Lightwave Technol. 23(7), 2309–2317 (2005). [CrossRef]

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

20. 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(1), 56–66 (2008). [CrossRef]

21. C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14(24), 11631–11652 (2006). [CrossRef] [PubMed]

22. P. J. Chiang and Y. C. Chiang, “Pseudospectral frequency-domain formulae based on modified perfectly matched layers for calculating both guided and leaky modes,” IEEE Photon. Technol. Lett. 22(12), 908–910 (2010). [CrossRef]

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

24. J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springer-Verlag, 2nd edition, 2001).

25. R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly re-started Arnoldi iteration,” SIAM J. Matrix Anal. Appl. 17(4), 789–821 (1996). [CrossRef]

26. T. Tang, “The Hermite spectral method for Gauss-type functions,” SIAM J. Sci. Comput. 14(3), 594–605 (1993). [CrossRef]

27. T. Tamir, Guides-Wave Optoelectronics (Springer-Verlag, 1988).

28. M. Loymeyer, N. Bahlmann, O. Zhuromskyy, H. Dotsch, and P. Hertel, “Phase-matched rectangular magnetooptic waveguides for applications in integrated optics isolators: numerical assessment,” Opt. Commun. 158(1-6), 189–200 (1998). [CrossRef]

29. B. Bellini and R. Beccherelli, “Modeling, design and analysis of liquid crystal waveguides in preferentially etched silicon grooves,” J. Phys. D Appl. Phys. 42(4), 045111 (2009). [CrossRef]

30. A. D’Álessandro, B. Bellini, D. Donisi, R. Beccherelli, and R. Asquini, “Nematic liquid crystal optical channel waveguides on silicon,” IEEE J. Quantum Electron. 42, 1084–1090 (2006). [CrossRef]

31. J. Beeckman, R. James, F. A. Í. Fernandez, W. De Cort, P. J. M. Vanbrabant, and K. Neyts, “Calculation of fully anisotropic liquid crystal waveguide modes,” J. Lightwave Technol. 27(17), 3812–3819 (2009). [CrossRef]

32. P. J. Vanbrabant, J. Beeckman, K. Neyts, R. James, and F. A. Fernandez, “A finite element beam propagation method for simulation of liquid crystal devices,” Opt. Express 17(13), 10895–10909 (2009). [CrossRef] [PubMed]

33. P. Yeh, and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons. Inc., New York, 1999).

N_int	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$	Time (sec)
5	2.049644	2.047895	0.8
10	2.048921	2.047191	1.7
15	2.048805	2.047074	3.6
20	2.048759	2.047028	6.8
25	2.048734	2.047003	13.4

	This study		Vector FD method [8]
ζ	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$
0.0001	2.047887	2.047850	2.047909	2.047873
0.0002	2.047904	2.047833	2.047927	2.047855
0.0005	2.047955	2.047782	2.047982	2.047801
0.0010	2.048042	2.047695	2.048072	2.047710
0.0020	2.048215	2.047522	2.048253	2.047530
0.0050	2.048734	2.047003	2.048795	2.046988

Mode number	The twist angle φ
Mode number	0°	30°	45°	60°	90°
1	1.674320	1.674065	1.673798	1.673512	1.673203
2	1.629201	1.626711	1.626975	1.628761	1.630778
3	1.620445	1.621682	1.620035	1.616712	1.612987
4	1.574271	1.574378	1.573683	1.572310	1.569640
5	1.553208	1.543317	1.546410	1.552893	1.560291
6	1.533944	1.540108	1.533665	1.523996	1.514011
7	1.501704	1.501848	1.501991	1.502157	1.502471

N_int	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$	Time (sec)
5	2.049644	2.047895	0.8
10	2.048921	2.047191	1.7
15	2.048805	2.047074	3.6
20	2.048759	2.047028	6.8
25	2.048734	2.047003	13.4

	This study		Vector FD method [8]
ζ	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$	$n_{e f f}^{(1)}$	$n_{e f f}^{(2)}$
0.0001	2.047887	2.047850	2.047909	2.047873
0.0002	2.047904	2.047833	2.047927	2.047855
0.0005	2.047955	2.047782	2.047982	2.047801
0.0010	2.048042	2.047695	2.048072	2.047710
0.0020	2.048215	2.047522	2.048253	2.047530
0.0050	2.048734	2.047003	2.048795	2.046988

Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach

Abstract

1. Introduction

2. Mathematical formulations

3. Numerical approach

4. Simulation results and discussion

4.1 A magneto-optical raised strip waveguide

4.2 A nematic liquid-crystal channel waveguide

5. Conclusion

Acknowledgements

References and links

Cited By

Figures (10)

Tables (3)

Equations (36)

Optics Express