Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method

doi:10.1364/OE.19.003363

« Show journal navigation

Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method

Chia-Chien Huang »View Author Affiliations

Optics Express, Vol. 19, Issue 4, pp. 3363-3378 (2011)
http://dx.doi.org/10.1364/OE.19.003363

View Full Text Article

Acrobat PDF (1515 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Computational schemes
4. Simulation results and discussion
5. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (34)
Cited By
Figures (11)
Metrics

Abstract

This study develops an efficient mode solver based on pseudospectral eigenvalue algorithm to analyze liquid crystal waveguides with full 3 × 3 anisotropic permittivity tensors. Present formulation yields a cubic eigenvalue matrix equation with an eigenvalue of the propagation constant, and they are solved using an iterative approach following the transformation of the matrix equation to a standard linear eigenvalue equation. The proposed scheme significantly reduces the memory storage and computational time by using only transverse magnetic field components. Although the proposed scheme requires an iterative procedure, the convergent eigenvalues are achieved after performing only four iterations. Therefore, for this scheme, computational efforts remain greatly lower than those for other reported schemes that used at least three field components. For solving the modes of nematic liquid crystal waveguides, the numerical results obtained by the proposed scheme are in good agreement with those calculated by using the finite-element and the finite-difference frequency-domain schemes, thus verifying the applicability of the proposed approach. Furthermore, the mode patterns of liquid crystal waveguides under arbitrary molecular orientations are also characterized in detail.

1. Introduction

Dielectric waveguides are essential building blocks of diverse photonic integrated circuits such as polarizers, switches, filters, rotators, and modulators. To accomplish the desired characteristics, these functional devices are mainly designed by altering the material characteristics by using electro-optic, magneto-optic, and thermo-optic modulations. Recently, liquid crystals (LCs) have been reconsidered as versatile optoelectronic materials for integrated optics [1

J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect,” Opt. Quantum Electron. 37(1-3), 95–106 (2005). [CrossRef]

–11

D. Donisi, B. Bellini, R. Beccherelli, R. Asquini, G. Gilardi, M. Trotta, and A. Dálessandro, “A switchable liquid-crystal optical channel waveguide on silicon,” IEEE J. Quantum Electron. 46(5), 762–768 (2010). [CrossRef]

] because of their low scattering loss at near-infrared wavelengths [12

M. Kawachi, N. Shibata, and T. Edahiro, “Possibility of use of liquid crystals as optical waveguide material for 1.3μm and 1.55μm bands,” Jpn. J. Appl. Phys. 21(Part 2, No. 3), L162–L164 (1982). [CrossRef]

] and low power requirement resulting from their appreciable birefringence. The birefringence can be controlled using an applied voltage from the designed electrode. For LCs with arbitrary molecular orientations, the nine elements of the permittivity tensor are all nonzero. In consequence, for accurately analyzing the LC device performance, it is indispensible to develop a numerical scheme that considers an arbitrary permittivity tensor.

For solving anisotropic LC waveguides with full 3 × 3 permittivity tensors by using full-wave techniques, a number of numerical schemes have been proposed. These approaches are the finite-element frequency-domain (FEFD)-based eigenvalue mode solvers [6

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]

,13

M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]

–15

S. M. Hsu and H. C. Chang, “Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy,” Opt. Express 15(24), 15797–15811 (2007). [CrossRef] [PubMed]

], the finite-element beam propagation methods (FE-BPMs) [3

G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finite-element beam propagation method,” Opt. Quantum Electron. 40(10), 733–748 (2008). [CrossRef]

P. J. M. 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]

,16

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]

,17

K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19(3), 405–413 (2001). [CrossRef]

], the finite-difference beam propagation method (FD-BPM) [4

E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707–5714 (2000). [CrossRef]

], and the finite-difference frequency-domain (FDFD)-based eigenvalue mode solvers [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]

]. In the FE-BPMs [3

,16

,17

], three electric (or magnetic) field components with 3N unknowns, are used to result in a matrix size of 9N ², where N is the number of grid points. To avoid introducing Pade approximation, the FE-BPM algorithm developed by Vanbrabant et al. [9

] used the transverse electric field components and their derivatives with respect to the propagation direction z to determine the second-order derivative term; however, this leads to an increase of up to 4N in the number of unknowns, resulting in a matrix size of 16N ². In the FD-BPM [4

E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707–5714 (2000). [CrossRef]

], a system of coupled differential equations, which involves an electric and a magnetic field component with 2N unknowns (the matrix size is 4N ²), is proposed for solving the twisted nematic pixels found in microdisplays. The required computational effort is relatively low. These BPM-based schemes [3

E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707–5714 (2000). [CrossRef]

,16

,17

] can study field profiles along the longitudinally variant waveguides but they need to consider the extra convergence conditions involving the step size of propagation, the two-point recurrence scheme, and the reference refractive index while solving the mode problems. Moreover, the BPM-based schemes cannot determine higher-order modes unless they transform the real propagation axis to the imaginary axis as in the imaginary-distance BPMs [18

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

,19

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]

]. Thus far, to the best of my knowledge, imaginary-distance BPM schemes have not been developed to solve full anisotropy waveguide problems. In contrast, mode solver algorithms based on eigenvalues, such as the FEFD [6

,13

M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]

–15

] and FDFD [10

] schemes, can deal with problems including all high-order modes. In [13

M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]

–15

], the FEFD schemes are formulated in term of all three magnetic or electric components with 3N unknowns to study the gyrotropic rectangular waveguides [13

M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]

], tilted channel waveguides [14

G. Tartarini and H. Renner, “Efficient finite-element analysis of tilted open anisotropic optical channel waveguides,”, ” IEEE Microw. Guid. Wave Lett. 9(10), 389–391 (1999). [CrossRef]

], and 2D photonic crystals [15

]. To maintain the symmetry of matrices in the FEFD method, Beeckman et al. [6

] used 6N unknowns (the matrix size is 36N ²) including three electric field components (3N) and the product (3N) of three electric field components and the wavevector component along z. To formulate a standard eigenvalue matrix equation in the FDFD [10

], the authors adopted 4N unknowns including two transverse electric field and two transverse magnetic field components to study LC waveguides with arbitrary molecular director orientations.

In recent years, the multidomain pseudospectral method [20

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

] has demonstrated fast convergence and high-order accuracy computational performance for solving propagation characteristics of various dielectric waveguides and optical fibers [21

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]

–27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

]. To achieve the same order of accuracy as the FDFD and FE algorithms based on eigenvalues, the pseudospectral scheme requires considerably less memory storage because it requires fewer unknown variables. However, in previous work, the authors considered only the isotropic media [21

–25

C. C. Huang, “Improved pseudospectral mode solver by prolate spheroidal wave functions for optical waveguides with step-index,” J. Lightwave Technol. 27(5), 597–605 (2009). [CrossRef]

] or at most the transverse anisotropic media [26

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]

,27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

]. The aim of this study is to develop an algorithm based on pseudospectral eigenvalues for determining the mode characteristics of anisotropic LCs with an arbitrary orientation of the LC director (i.e., having a full 3 × 3 anisotropic permittivity tensor). The proposed scheme is based on transverse magnetic field components. Thus, it significantly reduces the unknowns to only 2N, resulting in a matrix size of only 4N ², making the proposed scheme more efficient than approaches mentioned above. The extended Jones matrix method has also confirmed that only two of the electric or magnetic field components are independent [28

T. Scharf, Polarized Light in Liquid Crystals and Polymers (John Wiley and Sons. Inc., NJ, 2007, Chap 4).

]. However, the associated penalty is that one needs to deal with a cubic nonlinear eigenvalue equation corresponding to an eigenvalue of the propagation constant. Two methods exist for converting this cubic eigenvalue equation to a linear one. The first method involves extending the size of matrices by three times but the memory requirement is too heavy to be solved using a personal computer. The second method is an iterative approach that follows the arrangement of the nonlinear eigenvalue equation to a linear one with an eigenvalue of the square of propagation constant. In particular, the iterative scheme significantly reduces the computational resources, and thus, is able to be undertaken using a personal computer. The remainder of the paper is organized as follows. Section 2 presents the computational scheme including mathematical formulations and a numerical approach. Section 3 investigates mode characteristics of nematic LC waveguides with various orientations of the LC director; the numerical result demonstrates the superior convergence and computational efficiency of the proposed method. Finally, Section 4 draws conclusions.

2. Computational schemes

A. Mathematical formulations

Assuming a monochromatic electromagnetic wave that has a time dependence of exp(jωt) and is propagating along the z direction in a medium with the permittivity tensor [ε]. The vector wave equation based on the magnetic field vector H is given by

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

(1)

where ω is the angular frequency and μ ₀ is the permeability in vacuum. In this work, the permittivity tensor [ε] with the general form is considered as follows:

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

(2)

where ε ₀ is the permittivity in free space and [ε_r ] is the relative permittivity tensor. Considering a z-invariant waveguide structure, all electromagnetic field components are assumed to have a z dependence of exp(−jßz), where the propagation constant is β = k ₀ n_eff , the wave number in free space is k ₀ = ω ² μ ₀ ε ₀, and n_eff is the effective refractive index. In the study, the governing wave equations based on transverse magnetic field components H_x and H_y are developed, and the magnetic field component in the z-direction H_z is thus replaced by transverse magnetic field components H_x and H_y through applying the divergence-free of the magnetic field vector. Under the consideration of full anisotropy of [ε], a cubic eigenvalue matrix equation is formulated. By transforming the cubic eigenvalue matrix equation to a linear eigenvalue equation with an eigenvalue of β ², the full vector eigenvalue equation represented using the transverse magnetic field components H_x and H_y is obtained 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}],

(3)

where the differential operators P_xx (β), P_xy (β), P_yx (β), and P_yy (β) are defined as a function of β by the following.

\begin{array}{l} 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} + j β (η_{y z} + η_{z y}) \frac{\partial H_{x}}{\partial y} ... \\ + \frac{j}{β} (η_{z x} \frac{\partial}{\partial y} - η_{z y} \frac{\partial}{\partial x}) \frac{\partial^{2} H_{x}}{\partial y \partial x}, \end{array}

(4a)

\begin{array}{l} 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}} - j β (η_{y z} \frac{\partial H_{y}}{\partial x} + η_{z x} \frac{\partial H_{y}}{\partial y}) ... \\ + \frac{j}{β} (η_{z x} \frac{\partial}{\partial y} - η_{z y} \frac{\partial}{\partial x}) \frac{\partial^{2} H_{y}}{\partial y^{2}}, \end{array}

(4b)

\begin{array}{l} 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}} - j β (η_{x z} \frac{\partial H_{x}}{\partial y} + η_{z y} \frac{\partial H_{x}}{\partial x}) ... \\ + \frac{j}{β} (η_{z y} \frac{\partial}{\partial x} - η_{z x} \frac{\partial}{\partial y}) \frac{\partial^{2} H_{x}}{\partial x^{2}}, \end{array}

(4c)

\begin{array}{l} 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} + j β (η_{x z} + η_{z x}) \frac{\partial H_{y}}{\partial x} ... \\ + \frac{j}{β} (η_{z y} \frac{\partial}{\partial x} - η_{z x} \frac{\partial}{\partial y}) \frac{\partial^{2} H_{y}}{\partial x \partial y}, \end{array}

(4d)

where

η_{x x} = (ε_{y y} ε_{z z} - ε_{y z} ε_{z y}) / Δ,

η_{x y} = (ε_{x z} ε_{z y} - ε_{x y} ε_{z z}) / Δ,

η_{x z} = (ε_{x y} ε_{y z} - ε_{x z} ε_{y y}) / Δ,

η_{y x} = (ε_{y z} ε_{z x} - ε_{y x} ε_{z z}) / Δ,

η_{y y} = (ε_{x x} ε_{z z} - ε_{x z} ε_{z x}) / Δ,

η_{y z} = (ε_{x z} ε_{y x} - ε_{x x} ε_{y z}) / Δ,

η_{z x} = (ε_{y x} ε_{z y} - ε_{y y} ε_{z x}) / Δ,

η_{z y} = (ε_{x y} ε_{z x} - ε_{x x} ε_{z y}) / Δ,

η_{z z} = (ε_{x x} ε_{y y} - ε_{x y} ε_{y x}) / Δ,

and

Δ = ε_{x x} ε_{y y} ε_{z z} + ε_{x y} ε_{y z} ε_{z x} + ε_{x z} ε_{y x} ε_{z y} - ε_{x z} ε_{y y} ε_{z x} - ε_{x y} ε_{y x} ε_{z z} - ε_{x x} ε_{y z} ε_{z y} .

(5)

Compared with that in the transverse anisotropic media [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

], the operators in Eqs. (4)a) to (4d) including the eigenvalue β lead to a nonlinear eigenvalue matrix equation rather than a simple linear one as shown in [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

]. Therefore, an approach for solving a nonlinear eigenvalue problem is required. In the proposed scheme, the computational window is first divided at interfaces between different materials into several subdomains having uniform refractive index profiles and is then assembled using the physical interfacial conditions between these subdomains. Interfacial conditions include continuous normal and tangential components of the magnetic fields H_x and H_y at each intra-element boundary, and the other two interfacial conditions can be derived from the divergence condition of the magnetic field vector and from the Ampere’s law. In the divergence condition of the magnetic field

\nabla \cdot H = 0

, the longitudinal component H_z is expressed by

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

(6)

and that the E_z is expressed by

E_{z} = η_{z x} [(\frac{- 1}{β}) \frac{\partial}{\partial y} (\frac{\partial H_{x}}{\partial x} + \frac{\partial H_{y}}{\partial y}) + β H_{y}] + η_{z y} [(\frac{1}{β}) \frac{\partial}{\partial x} (\frac{\partial H_{x}}{\partial x} + \frac{\partial H_{y}}{\partial y}) - β H_{x}] + j η_{z z} (\frac{\partial H_{x}}{\partial y} - \frac{\partial H_{y}}{\partial x}),

(7)

through the Ampere’s law

\nabla \times H =

jω[ε]E. Accordingly, Eqs. (6) and (7) are used as the two interfacial conditions because of the continuities of H_z and E_z at interfaces between different materials. Note that the first three continuous interface conditions, the normal and tangential components of the magnetic fields (H_x and H_y ) and H_z , are the same as those applied in the transverse anisotropic media [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

]. The main difference of interface conditions between the full 3 × 3 anisotropic here and the transverse anisotropic media [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

] is appeared in Eq. (7). Equation (7) involving the eigenvalue β makes the proposed work significantly complex while patching the subdomains than that in [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

] which has only the last term of Eq. (7). In the proposed approach, only horizontal and vertical interfaces are considered because the aim of this study is to develop an efficient waveguide mode solver for anisotropic waveguides with arbitrary permittivity tensors. If the waveguides, such as optical fibers, tapered waveguides, or photonic crystals, with curved boundaries are considered, the transfinite interpolation schemes with bending functions [29

W. J. Gordon and C. A. Hall, “Transfinite element methods: blending function interpolation over arbitrary curved element domains,” Numer. Math. 21(2), 109–129 (1973). [CrossRef]

] are the most often used approaches in various numerical schemes such as the FD, FE, and pseudospectral methods [30

D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations : Algorithms for Scientists and Engineers (Springer–Verlag, 2009).

]. The schemes use linear interpolations to map a reference orthogonal square in computational space to an arbitrarily curvilinear quadrilateral region that represent the physical space. Under the mapping, the Eqs. (3), (6), and (7) are transformed by the result of chain rule. The implementations of the transfinite interpolation schemes in the pseudospectral approach for analyzing the optical fibers and photonic crystals with curved interfaces, readers are referred to references [22

P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75(2 Pt 2), 026703 (2007). [CrossRef] [PubMed]

,23

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]

B. Numerical approach

In the pseudospectral scheme [20

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

], the computational window is divided at interfaces between different materials into several subdomains having uniform 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 The mesh division of an arbitrary interior subdomain r

In the subdomain r, the magnetic field components H_x and H_y are expanded using a set of suitable orthogonal basis functions 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},

(8a)

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},

(8b)

where θ(x) and ψ(y) are Lagrange-type basis functions, H_x,pq and H_y,pq are grid point values at (n_x + 1) (n_y + 1) collocation points, and

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

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

, δ_mp denotes the Kronecker delta. Substituting Eqs. (8)a) and (8b) into Eq. (3), Eq. (3) is demanded to perfectly be satisfied at these (n_x − 1) × (n _y − 1) interior collocation points (namely, (x ₁, y_q ), (x ₂, y_q ),… (x ₉, y_q ), as shown in Fig. 1, where q = 1 to 9). Equation (3) is thus converted to a system of linear equations to form a matrix eigenvalue equation with an eigenvalue of β ² 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}]

(9)

where the operators

P_{x x}^{r} (β)

P_{x y}^{r} (β)

P_{y x}^{r} (β)

, and

P_{y y}^{r} (β)

are

P_{u v}^{r} (β) = [\begin{matrix} {({\tilde{D}}_{00}^{n_{x} - 1, 1})}_{_{u v}}^{r} & {({\tilde{D}}_{01}^{n_{x} - 1, 1})}_{_{u v}}^{r} & \dots & {({\tilde{D}}_{0 n_{y}}^{n_{x} - 1, 1})}_{_{u v}}^{r} \\ {({\tilde{D}}_{00}^{n_{x} - 1, 2})}_{_{u v}}^{r} & {({\tilde{D}}_{01}^{n_{x} - 1, 2})}_{_{u v}}^{r} & \dots & {({\tilde{D}}_{0 n_{y}}^{n_{x} - 1, 2})}_{_{u v}}^{r} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ {({\tilde{D}}_{00}^{n_{x} - 1, n_{y} - 1})}_{_{u v}}^{r} & {({\tilde{D}}_{01}^{n_{x} - 1, 2})}_{_{u v}}^{r} & \dots & {({\tilde{D}}_{0 n_{y}}^{n_{x} - 1, 2})}_{_{u v}}^{r} \end{matrix}], u v = x x, x y, y x, and y y,

(10a)

{({\tilde{D}}_{0 s}^{n_{x} - 1, t})}_{_{u v}}^{r} = [\begin{matrix} {(D_{0 s}^{1, t})}_{_{u v}}^{r} & {(D_{1 s}^{1, t})}_{_{u v}}^{r} & \dots {(D_{n_{x} s}^{1, t})}_{_{u v}}^{r} \\ {(D_{0 s}^{2, t})}_{_{u v}}^{r} & {(D_{1 s}^{2, t})}_{_{u v}}^{r} & \dots {(D_{n_{x} s}^{2, 1})}_{_{u v}}^{r} \\ {(D_{0 s}^{n_{x} - 1, t})}_{_{u v}}^{r} & {(D_{1 s}^{n_{x} - 1, t})}_{_{u v}}^{r} & \dots {(D_{n_{x} s}^{n_{x} - 1, t})}_{_{u v}}^{r} \end{matrix}], s = 0, 1, ..., n_{y}, t = 0, 1, ..., n_{x},

(10b)

\begin{array}{l} {(D_{p q}^{i, j})}_{x x}^{r} = η_{y y}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) + η_{z z}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r}) - η_{y x}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) + k_{0}^{2} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) ... \\ + j β (η_{y z}^{r} + η_{z y}^{r}) θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) + \frac{j}{β} (η_{z x}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r}) - η_{z y}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r})), \end{array}

(10c)

\begin{array}{l} {(D_{p q}^{i, j})}_{x y}^{r} = (η_{y y}^{r} - η_{z z}^{r}) θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) - η_{y x}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r}) - j β (η_{y z}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) ... \\ + η_{z x}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r})) + \frac{j}{β} (η_{z x}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (3)} (y_{_{j}}^{r}) - η_{z y}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r})), \end{array}

(10d)

\begin{array}{l} {(D_{p q}^{i, j})}_{y x}^{r} = (η_{x x}^{r} - η_{z z}^{r}) θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) - η_{x y}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) - j β (η_{x z}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) ... \\ + η_{z y}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r})) + \frac{j}{β} (η_{z y}^{r} θ_{p}^{r (3)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) - η_{z x}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r})), \end{array}

(10e)

\begin{array}{l} {(D_{p q}^{i, j})}_{y y}^{r} = η_{z z}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) + η_{x x}^{r} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r}) + η_{x y}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) + k_{0}^{2} θ_{p}^{r} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) ... \\ + j β (η_{x z}^{r} + η_{z x}^{r}) θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r} (y_{_{j}}^{r}) + \frac{j}{β} (η_{z y}^{r} θ_{p}^{r (2)} (x_{_{i}}^{r}) ψ_{q}^{r (1)} (y_{_{j}}^{r}) - η_{z x}^{r} θ_{p}^{r (1)} (x_{_{i}}^{r}) ψ_{q}^{r (2)} (y_{_{j}}^{r})), \end{array}

(10f)

where

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

and

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

represent the h-th order derivatives of θ_p (x) to x and ψ_q (y) to y, respectively. Once the matrix eigenvalue equation for each subdomain is obtained, a global matrix equation can be formed by assembling the matrix equations in all subdomains. Assuming there are m subdomains, the pattern of the matrix elements is given by

[\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}]

(11)

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) .

(12)

Note that the interfacial conditions between subdomains have not yet been enforced to Eq. (11) to form the final matrix equation. The exponential convergence behavior can be preserved as the same as that in the single domain version by explicitly enforcing the physical interface conditions. First, the continuities of the normal and tangential components of the magnetic field vector H are implicitly imposed. In addition, the other two interface conditions are the continuities of H_z in Eq. (6) and E_z in Eq. (7). For a specific vertical interface as shown in Fig. 1, the interface collocation points (x ₀, y_j ) and (x_nx , y_j ) belong to the right (denoted by + ) and left (denoted by −) subdomains, respectively, where j = 1, 2, 3… n_y . The continuity of H_z gives

[\begin{matrix} A_{^{x}}^{+} & A_{^{y}}^{+} \end{matrix}] [\begin{matrix} H_{x}^{+} \\ H_{y}^{+} \end{matrix}] = [\begin{matrix} A_{^{x}}^{-} & A_{^{y}}^{-} \end{matrix}] [\begin{matrix} H_{x}^{-} \\ H_{y}^{-} \end{matrix}],

(13)

where

A_{k}^{\pm} = [\begin{matrix} A_{k}^{\pm} (ψ_{0} (y_{0}^{\pm})) & A_{k}^{\pm} (ψ_{1} (y_{0}^{\pm})) & \dots & A_{k}^{\pm} (ψ_{n_{y}} (y_{0}^{\pm})) \\ A_{k}^{\pm} (ψ_{0} (y_{1}^{\pm})) & A_{k}^{\pm} (ψ_{1} (y_{1}^{\pm})) & \dots & A_{k}^{\pm} (ψ_{n_{y}} (y_{1}^{\pm})) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ A_{k}^{\pm} (ψ_{0} (y_{n_{y}}^{\pm})) & A_{k}^{\pm} (ψ_{1} (y_{n_{y}}^{\pm})) & \dots & A_{k}^{\pm} (ψ_{n_{y}} (y_{n_{y}}^{\pm})) \end{matrix}], k = x, y,

(14a)

A_{x}^{+} (ψ_{i} (y_{j}^{+})) = [\begin{matrix} θ_{0}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) & θ_{1}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) & \dots θ_{n_{x}}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) \end{matrix}],

(14b)

A_{y}^{+} (ψ_{i} (y_{j}^{+})) = [\begin{matrix} θ_{0}^{} (x_{_{0}}^{+}) ψ_{i}^{(1)} (y_{j}^{+}) & θ_{1}^{} (x_{_{0}}^{+}) ψ_{i}^{(1)} (y_{j}^{+}) & \dots θ_{n_{x}} (x_{_{0}}^{+}) ψ_{i}^{(1)} (y_{j}^{+}) \end{matrix}],

(14c)

A_{x}^{-} (ψ_{i} (y_{j}^{-})) = [\begin{matrix} θ_{0}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) & θ_{1}^{(1)} (x_{n_{x}}^{-}) ψ_{i} (y_{j}^{-}) & \dots θ_{n_{x}}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) \end{matrix}],

(14d)

A_{y}^{-} (ψ_{i} (y_{j}^{-})) = [\begin{matrix} θ_{0}^{} (x_{_{n_{x}}}^{-}) ψ_{i}^{(1)} (y_{j}^{-}) & θ_{1}^{} (x_{n_{x}}^{-}) ψ_{i}^{(1)} (y_{j}^{-}) & \dots θ_{n_{x}} (x_{n_{x}}^{-}) ψ_{i}^{(1)} (y_{j}^{-}) \end{matrix}],

(14e)

and the continuity of E_z gives

[\begin{matrix} B_{^{x}}^{+} & B_{^{y}}^{+} \end{matrix}] [\begin{matrix} H_{x}^{+} \\ H_{y}^{+} \end{matrix}] = [\begin{matrix} B_{^{x}}^{-} & B_{^{y}}^{-} \end{matrix}] [\begin{matrix} H_{x}^{-} \\ H_{y}^{-} \end{matrix}],

(15)

where

B_{k}^{\pm} = [\begin{matrix} B_{k}^{\pm} (ψ_{0} (y_{0}^{\pm})) & B_{k}^{\pm} (ψ_{1} (y_{0}^{\pm})) & \dots B_{k}^{\pm} (ψ_{n_{y}} (y_{0}^{\pm})) \\ B_{k}^{\pm} (ψ_{n_{y}} (y_{0}^{\pm})) & B_{k}^{\pm} (ψ_{1} (y_{1}^{\pm})) & \dots B_{k}^{\pm} (ψ_{n_{y}} (y_{1}^{\pm})) \\ ⋮ & ⋮ & ⋮ \\ B_{k}^{\pm} (ψ_{0} (y_{n_{y}}^{\pm})) & B_{k}^{\pm} (ψ_{1} (y_{n_{y}}^{\pm})) & \dots B_{k}^{\pm} (ψ_{n_{y}} (y_{n_{y}}^{\pm})) \end{matrix}], k = x, y,

(16a)

B_{x}^{+} (ψ_{i} (y_{j}^{+})) = {[\begin{matrix} η_{z x} (\frac{- 1}{β}) θ_{0}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) + η_{z y} (\frac{1}{β}) θ_{0}^{(2)} (x_{0}^{+}) ψ_{i} (y_{j}^{+}) - β η_{z y} δ_{i j} + j η_{z z} θ_{0} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) \\ η_{z x} (\frac{- 1}{β}) θ_{1}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) + η_{z y} (\frac{1}{β}) θ_{1}^{(2)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) - β η_{z y} δ_{i j} + j η_{z z} θ_{1} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) \\ ⋮ \\ η_{z x} (\frac{- 1}{β}) θ_{n_{x}}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) + η_{z y} (\frac{1}{β}) θ_{n_{x}}^{(2)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) - β η_{z y} δ_{i j} + j η_{z z} θ_{n_{x}} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) \end{matrix}]}^{T},

(16b)

B_{y}^{+} (ψ_{i} (y_{j}^{+})) = {[\begin{matrix} η_{z x} (\frac{- 1}{β}) θ_{0} (x_{_{0}}^{+}) ψ_{_{i}}^{(2)} (y_{j}^{+}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{0}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) - j η_{z z} θ_{_{0}}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) \\ η_{z x} (\frac{- 1}{β}) θ_{1} (x_{_{0}}^{+}) ψ_{_{i}}^{(2)} (y_{j}^{+}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{1}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) - j η_{z z} θ_{_{1}}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) \\ ⋮ \\ η_{z x} (\frac{- 1}{β}) θ_{n_{x}} (x_{_{0}}^{+}) ψ_{_{i}}^{(2)} (y_{j}^{+}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{n_{x}}^{(1)} (x_{_{0}}^{+}) ψ_{_{i}}^{(1)} (y_{j}^{+}) - j η_{z z} θ_{_{n_{x}}}^{(1)} (x_{_{0}}^{+}) ψ_{i} (y_{j}^{+}) \end{matrix}]}^{T},

(16c)

B_{x}^{-} (ψ_{i} (y_{j}^{-})) = {[\begin{matrix} η_{z x} (\frac{- 1}{β}) θ_{0}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) + η_{z y} (\frac{1}{β}) θ_{0}^{(2)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) - β η_{z y} δ_{i j} + j η_{z z} θ_{0} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) \\ η_{z x} (\frac{- 1}{β}) θ_{1}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) + η_{z y} (\frac{1}{β}) θ_{1}^{(2)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) - β η_{z y} δ_{i j} + j η_{z z} θ_{1} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) \\ ⋮ \\ η_{z x} (\frac{- 1}{β}) θ_{n_{x}}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) + η_{z y} (\frac{1}{β}) θ_{n_{x}}^{(2)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) - β η_{z y} δ_{i j} + j η_{z z} θ_{n_{x}} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) \end{matrix}]}^{T},

(16d)

B_{y}^{-} (ψ_{i} (y_{j}^{-})) = {[\begin{matrix} η_{z x} (\frac{- 1}{β}) θ_{0} (x_{n_{x}}^{-}) ψ_{_{i}}^{(2)} (y_{j}^{-}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{0}^{(1)} (x_{n_{x}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) - j η_{z z} θ_{_{0}}^{(1)} (x_{n_{x}}^{-}) ψ_{i} (y_{j}^{-}) \\ η_{z x} (\frac{- 1}{β}) θ_{1} (x_{n_{x}}^{-}) ψ_{_{i}}^{(2)} (y_{j}^{-}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{1}^{(1)} (x_{n_{x}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) - j η_{z z} θ_{_{1}}^{(1)} (x_{n_{x}}^{-}) ψ_{i} (y_{j}^{-}) \\ ⋮ \\ η_{z x} (\frac{- 1}{β}) θ_{n_{x}} (x_{n_{x}}^{-}) ψ_{_{i}}^{(2)} (y_{j}^{-}) + β η_{z x} δ_{i j} + η_{z y} (\frac{1}{β}) θ_{n_{x}}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{_{i}}^{(1)} (y_{j}^{-}) - j η_{z z} θ_{_{n_{x}}}^{(1)} (x_{_{n_{x}}}^{-}) ψ_{i} (y_{j}^{-}) \end{matrix}]}^{T},

(16e)

where i = 0,1 2, …n_x and T denotes the transpose of a matrix. The derived formulations are similar for a horizontal interface, if the vertical interface collocation points (x ₀, y_j ) and (x_nx , y_j ) are altered by the horizontal ones (x_i , y ₀) and (x_i , y_ny ). For brevity, these formulations are not shown here. Following the global matrix equation of Eq. (11), in addition to the continuous normal and tangential components of the magnetic fields at all interface points, the final matrix equation is formed by also imposing the interfacial conditions in Eqs. (13) and (15) to Eq. (11), which no longer has a block diagonal form but turns into a matrix with an approximate sparsity of 47%. A detailed description for clearly recognizing the coupling of fields between subdomains can be found in the appendix of the work [22

], which solved the band structures of photonics crystals by scalar pseudospectral method. Notably, solving the full 3 × 3 anisotropic waveguides is approximately 8.5 times the memory required and 1.9 times the computational time of those for solving the transverse anisotropic ones [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

] having the approximate sparsity of 38%. The dependent variables in the final matrix include only the interior collocation points of each subdomain; however, the dependent variables in interface points can be obtained once the solutions of the final matrix are solved. In addition, the computational boundary conditions are automatically satisfied because of the use of LGFs to represent the exterior subdomains.

We now determine the appropriate basis functions by expanding the H_x and H_y components. Chebyshev polynomials are used for expanding the optical fields in interior subdomains having finite intervals because of their mathematical robustness to non-periodic structures. In contrast, LGFs are used for expanding the exterior subdomains having semi-infinite intervals because the exponential decay characteristic of LGFs matches the guided field profiles. For Chebyshev polynomials, the explicit form of the Lagrange-type interpolation function is represented as follows [20

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

θ_{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}

(17)

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

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

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

(18)

where L_v (αx) is the Laguerre polynomial of order v, and x_p denotes the collocation points for the LGFs. The scaling factor α in Eq. (18) significantly affects the numerical accuracy of the results obtained using LGFs. For a given number of terms of basis functions, α is determined by the derivation of Tang [31

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

] and the idea of the effective index method (EIM) [32

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

]. From the derivation [31

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

], the optimum α is determined by the following formula:

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

(19)

where M is the spreading of the guided fields to be solved, and x_i is the collocation point at position i. Now, the determination of M is become a main concern for various physical problems while using LGFs. For a light beam propagating in waveguides, the M of guided modes can be computed by a simple EIM [32

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

]. The detailed derivation of determining M has been shown in the previous work [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

] for solving optical waveguides with transverse anisotropic media. The processes for finding α is the same for the full 3 × 3 anisotropic permittivity tensors. The resultant global matrix is a cubic eigenvalue equation with respect to β, which results from the wave equations of Eq. (11) and the interfacial condition of Eq. (15). By transforming the final matrix to a linear eigenvalue equation with an eigenvalue of β ², a simple iterative scheme such as the following process is used to discover the mode eigenvalues. (1) Assign an initially guessed value of n_eff by that of the substrate (in fact, this method should be followed only if a reasonable range of values of n_eff is required). (2) Calculate the operators P _xx (β), P _xy (β), P _yx (β), and P _yy (β). (3) A standard linear eigenvalue problem with an eigenvalue of β ² is obtained and solved using the Arnoldi iteration method [33

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]

], which efficiently solves generalized eigenvalue problems. (4) After obtaining a new n_eff in step (3), repeatedly keep executing the processes (2) and (3) to achieve the desired convergence (the difference of n_eff between adjacent iterations is set at 10⁻⁸ in this study). Particularly, four iterations are needed to achieve the accepted convergence solutions due to the nonlinearity of the eigenvalue equation for the full anisotropy, and thus the whole computational time is 8 times of solving a transverse anisotropic problem [27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

4. Simulation results and discussion

In this section, a square LC waveguide with five nonzero elements of the permittivity tensor is first calculated to examine the convergence of the proposed mode solver compared with that by using the FEFD-based eigenvalue approach with higher-order elements [6

]. Furthermore, to show the full capability of the proposed scheme, the second example considers a nematic LC waveguide with nine nonzero elements of the permittivity tensor. The calculated effective indices are then compared with those obtained using the FDFD-based eigenvalue approach [10

]. The mode patterns under various orientations of the LC director are also characterized in detail.

Figure 2(a) shows a cross-sectional view of the square LC waveguide surrounded by an isotropic cladding with refractive index n_s = 1.5. The structure parameters of its core region are: width w = 1 μm and thickness t = 1 μm.

Fig. 2 (a) Cross-sectional view of the square LC waveguide with the permittivity tensor [ε] in the core region and refractive index n_s of the surrounding substrate. (b) Schematic diagram of the orientation

\hat{n}

of the LC director. (c) Division of the computational domain for the LC waveguide.

Assuming that the orientation of the square LC director

\hat{n}

tilts at an angle of θ_c = 45° with respect to the z-axis (the propagation direction) as shown in Fig. 2(b), and thus, the permittivity tensor [ε] of the LC is given by

[ε] = ε_{0} [\begin{matrix} n_{o}^{2} & 0 & 0 \\ 0 & n_{o}^{2} + Δ ε \sin^{2} θ_{c} & Δ ε \sin θ_{c} \cos θ_{c} \\ 0 & Δ ε \sin θ_{c} \cos θ_{c} & n_{o}^{2} + Δ ε \cos^{2} θ_{c} \end{matrix}]

(20)

where n_o = 1.55, n_e = 1.8, and

Δ ε = n_{e}^{2} - n_{o}^{2}

. In this example, the operating wavelength λ = 1μm in free space. To preserve the fast convergence of the pseudospectral method, the proposed scheme divides the computational window at material interfaces into 9 subdomains (Fig. 2(c)). The guided mode profiles in the interior subdomains having finite intervals are expanded using Chebyshev polynomials, and those in the exterior subdomains having semi-infinite intervals are expanded using LGFs. For example, the mode profiles in both the x- and y-directions of subdomains 1, 3, 7, and 9 are expanded using LGFs. In contrast, the Chebyshev polynomials expand the mode profiles in both the x- and y-directions of the subdomain 5. For subdomains 2 and 8, the mode profiles in the x-direction are expanded using Chebyshev polynomials, and those in the y-direction are expanded using LGFs.

For examining the convergence of the proposed scheme, Fig. 3 shows the effective indices of the i-th iteration of the fundamental mode versus the iteration time under different number of terms of basis functions n. We note that equal number of terms of the basis function in both the x- and y- directions for all subdomains are used throughout the paper, and the chosen initial guess of the effective index is set as n _guess = n_o = 1.55. In particular, the semi-infinite intervals are fixed at 10 terms of the basis functions due to the well-matched behaviors of the exponential decay fields of guided modes and the mathematical behavior of LGFs. In Fig. 3, the ns denote the number of terms required to expand the finite intervals. Thus, the total number of unknowns used are 2 × [(10 − 1) + (n − 2) + (10 − 1)]².

Fig. 3 Convergent effective index versus iteration time for different terms of the basis functions n.

The differences of the effective indices between the fourth and fifth iterations are achieved in the order of 10⁻⁸, and the calculated effective indices are 1.604108 and 1.604152 for n = 20 (the number of unknowns is 2592) and n = 25 (the number of unknowns is 3362) at the fourth iteration, respectively. Compared with the effective index n_eff = 1.6055 obtained using the FEFD-based eigenvalue approach with higher-order elements [6

] (the degrees of freedom used are ~10000), the result n_eff = 1.6041 calculated using the proposed scheme shows good agreement. However, the proposed scheme has the advantage that it requires fewer degrees of freedom (less memory storage). In addition, Fig. 4 shows the mode profiles of |H_x | and |H_y |.

Fig. 4 Mode profiles of (a) |H_x | and (b) |H_y | of the fundamental mode of the square nematic LC waveguide.

The following example considers an LC dielectric waveguide comprising a core region filled with nematic LCs (5CB) [34

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

], surrounded by a glass substrate having refractive index n_g = 1.45 and air with refractive index n_a , (Fig. 5(a) ). In addition, this example investigates arbitrary orientations of the LC director shown in Fig. 5(b).

Fig. 5 (a) Cross-sectional view of the nematic LC-core channel waveguide with the core permittivity tensor [ε], substrate refractive index n_g , and air refractive index n_a . (b) Schematic diagram of the twist angle φ_c and the tilt angle θ_c of the LC director

\hat{n}

The permittivity tensor [ε] of the LC-core region at the operating wavelength λ = 1.55 μm is given by

[ε] = ε_{0} [\begin{matrix} n_{o}^{2} + Δ ε \sin^{2} θ_{c} \cos^{2} φ_{c} & Δ ε \sin^{2} θ_{c} \sin φ_{c} \cos φ_{c} & Δ ε \sin θ_{c} \cos θ_{c} \cos φ_{c} \\ Δ ε \sin^{2} θ_{c} \sin φ_{c} \cos φ_{c} & n_{o}^{2} + Δ ε \sin^{2} θ_{c} \sin^{2} φ_{c} & Δ ε \sin θ_{c} \cos θ_{c} \sin φ_{c} \\ Δ ε \sin θ_{c} \cos θ_{c} \cos φ_{c} & Δ ε \sin θ_{c} \cos θ_{c} \sin φ_{c} & n_{o}^{2} + Δ ε \cos^{2} θ_{c} \end{matrix}],

(21)

where

Δ ε = n_{e}^{2} - n_{o}^{2}

, the ordinary refractive index is n_o = 1.5292, the extraordinary refractive index is n_e = 1.7072, θ_c is the tilt angle between the LC director

\hat{n}

and the z-axis, and φ_c is the twist angle between the projection of

\hat{n}

on the xy plane and the x-axis. The width and thickness of the core are w = 3 μm and t = 3 μm, respectively. The orientation

\hat{n}

of the LC director can be controlled using an applied voltage on the designed electrode. For the condition θ_c = 90°, the elements of [ε] under different angles of φ_c are reduced to those of a waveguide with transverse anisotropy [26

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]

,27

C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]

]. For the condition θ_c = 0°, nonzero elements of [ε] are further reduced to be only in diagonal positions. This study accordingly considers the conditions θ_c = 30° and θ_c = 60° in which nine elements of [ε] are all nonzero. The computational window is divided into nine subdomains, which is similar to that in Fig. 2(c). Furthermore, all calculations use n = 20, and results are obtained at the fourth iteration. The calculated effective indices for θ_c = 30° and θ_c = 60° under different φ_c ’s along with the results obtained using the FDFD [10

] are shown in Figs. 6(a) and 6(b), respectively. On comparison with the calculated effective indices by the FDFD [10

] (indicated by circles in Fig. 6), the results obtained using the proposed scheme (indicated by crosses in Fig. 6) are found to be in good agreement.

Fig. 6 Calculated effective indices versus different values of φ_c for tilt angles (a) θ_c = 30° and (b) θ_c = 60°.

For each individual value of the tilt angles θ_c , the differences in the effective indices of the fundamental guided mode (mode 1) are fairly small (of the order of 10⁻⁴) for different values of twist angles φ_c . This can be understood by the opposite variances of the diagonal elements ε_xx and ε_yy as well as the non-diagonal elements ε_xz and ε_yz as the twist angles φ_c alter. The major field patterns of mode 1 for four values of φ_c shown in Fig. 7 (for θ_c = 30°) and Fig. 8 (for θ_c = 60°) also verify these results. The symbols |H_1,x | and |H_1,y | on the left sides of Figs. 7 and 8 denote the moduli of the H_x and H_y components of mode 1, respectively. The effects of altering φ_c emerge in the minor field patterns and the relative amplitude between |H_x | and |H_y |. In addition, Figs. 7 and 8 also show the effect of transforming the major field H_y to H_x beyond the twist angle φ_c = 45°.

Fig. 7 |H_x | and |H_y | mode patterns of mode 1 for the tilt angle θ_c = 30°

Fig. 8 |H_x | and |H_y | mode patterns of mode 1 for the tilt angle θ_c = 60°.

For higher-order modes, the differences in the values of the calculated effective indices for different values of φ_c are greater than 0.001; thus, the major field patterns show modest variations. Figure 9 shows the |H_x | and |H_y | components of higher-order modes 2–4 for four values of φ_c at θ_c = 30°; Fig. 10 shows the same for modes 2−4 at θ_c = 60°.

Fig. 9 |H_x | and |H_y | mode patterns of the higher-order modes for the tilt angle θ_c = 30°.

Fig. 10 |H_x | and |H_y | mode patterns of the higher-order modes for the tilt angle θ_c = 60°.

As shown in Figs. 9 and 10, mode symmetries relative to the y-axis are observed only at φ_c = 0° and φ_c = 90°, and symmetric features also appear at θ_c = 45°. Particularly, the major field does not always have the same polarization for all modes for an individual value of φ_c . For instance, at the orientation of the LC director θ_c = 30° and φ_c = 90°, the major field is H_x polarization for the first three modes but turns into H_y polarization for mode 4. Contrary to the condition of φ_c = 90°, only mode 4 has major H_x field while the orientation is at φ_c = 0°. For the case of θ_c = 60° and φ_c = 90°, the major field H_y polarizations for the first six modes are converted to the H_x polarization for mode 7 (not shown here). Further observations were made of the higher-order modes; however, the reverse phenomenon appears merely in some modes. In addition, for examining the influence of θ_c at certain angles of φ_c in more detail, the calculated effective indices by the proposed approach versus different values of θ_c at φ_c = 30° and φ_c = 60° are shown in Fig. 11(a) and (b) , respectively. The effective index increases monotonously as the angle θ_c increases. For a specific angle of φ_c , it can be realized that the elements ε_xx , ε_yy , and ε_xy increase as the angle of θ_c increases.

Fig. 11 Calculated effective indices versus different values of θ_c for twisted angles (a) φ_c = 30° and (b) φ_c = 60°.

5. Conclusion

This study has developed an efficient pseudospectral mode solver for solving dielectric waveguides with full 3 × 3 anisotropy. The proposed scheme was formulated as an eigenvalue matrix equation by applying only the transverse magnetic field components with 2N unknowns (matrix size of 4N ²), where N is the number of unknowns. Thus, it requires fewer grid points than both the finite-element frequency-domain (FEFD)-based eigenvalue scheme with 3N unknowns (matrix size of 9N ²) and the finite-difference frequency-domain (FDFD)-based eigenvalue scheme with 4N unknowns (matrix size of 16N ²). In consequence, the computational time and memory required for the proposed scheme are significantly reduced. Having to treat a cubic eigenvalue equation creates a computational penalty because of the complicated wave equations and the interfacial conditions; however, convergence of the effective indices can be achieved by performing only four iterations. As a result, the computational efforts for the proposed scheme are still much lesser than those required for reported FEFD and FDFD approaches. The first example considered a square LC waveguide with the orientation of the director at an tilt angle of θ_c = 45° with respect to the propagation direction (having five nonzero elements in the permittivity tensor), and the calculated results of the proposed scheme are compared with those obtained using the FEFD scheme. A good agreement of the effective index of the fundamental mode is obtained to validate the numerical accuracy of the proposed scheme. To further demonstrate the full capability of the proposed scheme, in this study, I analyzed a nematic liquid crystal with arbitrary orientations of the director embedded in a glass substrate. Good agreement was also shown for the calculated results when compared with those from the FDFD scheme. In addition, mode profiles at arbitrary orientations of the LC director were presented. On examining these mode profiles, the proposed scheme was shown to offer not only an accurate numerical scheme but also a solution method that is more efficient than the previously reported approaches. Finally, the proposed mode solver can be used in the future to investigate and design more complex photonic devices while using arbitrary anisotropic materials with nine nonzero elements of the permittivity tensor.

Acknowledgements

The author would like to thank the National Science Council of China, Taiwan, for financially supporting this research under Contract No. NSC 99-2112-M-005-005-MY3.

References and links

1.	J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect,” Opt. Quantum Electron. 37(1-3), 95–106 (2005). [CrossRef]
2.	A. D’Alessandro, B. D. Donisi, R. Beccherelli, and R. Asquini, “Nematic liquid crystal optical channel waveguides on silicon,” IEEE J. Quantum Electron. 42(10), 1084–1090 (2006). [CrossRef]
3.	G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finite-element beam propagation method,” Opt. Quantum Electron. 40(10), 733–748 (2008). [CrossRef]
4.	E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707–5714 (2000). [CrossRef]
5.	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]
6.	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]
7.	P. J. M. 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]
8.	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]
9.	P. J. M. 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]
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.	D. Donisi, B. Bellini, R. Beccherelli, R. Asquini, G. Gilardi, M. Trotta, and A. Dálessandro, “A switchable liquid-crystal optical channel waveguide on silicon,” IEEE J. Quantum Electron. 46(5), 762–768 (2010). [CrossRef]
12.	M. Kawachi, N. Shibata, and T. Edahiro, “Possibility of use of liquid crystals as optical waveguide material for 1.3μm and 1.55μm bands,” Jpn. J. Appl. Phys. 21(Part 2, No. 3), L162–L164 (1982). [CrossRef]
13.	M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]
14.	G. Tartarini and H. Renner, “Efficient finite-element analysis of tilted open anisotropic optical channel waveguides,”, ” IEEE Microw. Guid. Wave Lett. 9(10), 389–391 (1999). [CrossRef]
15.	S. M. Hsu and H. C. Chang, “Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy,” Opt. Express 15(24), 15797–15811 (2007). [CrossRef] [PubMed]
16.	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]
17.	K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19(3), 405–413 (2001). [CrossRef]
18.	J. C. Chen and S. Jüngling, “Computation of high-order waveguide modes by imaginary-distance beam propagation method,” Opt. Quantum Electron. 26, 199–205 (1994). [CrossRef]
19.	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]
20.	J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springer–Verlag, 2nd edition, 2001).
21.	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]
22.	P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75(2 Pt 2), 026703 (2007). [CrossRef] [PubMed]
23.	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]
24.	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]
25.	C. C. Huang, “Improved pseudospectral mode solver by prolate spheroidal wave functions for optical waveguides with step-index,” J. Lightwave Technol. 27(5), 597–605 (2009). [CrossRef]
26.	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]
27.	C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]
28.	T. Scharf, Polarized Light in Liquid Crystals and Polymers (John Wiley and Sons. Inc., NJ, 2007, Chap 4).
29.	W. J. Gordon and C. A. Hall, “Transfinite element methods: blending function interpolation over arbitrary curved element domains,” Numer. Math. 21(2), 109–129 (1973). [CrossRef]
30.	D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations : Algorithms for Scientists and Engineers (Springer–Verlag, 2009).
31.	T. Tang, “The Hermite spectral method for Gauss-type functions,” SIAM J. Sci. Comput. 14(3), 594–605 (1993). [CrossRef]
32.	T. Tamir, Guides-Wave Optoelectronics (Springer-Verlag, 1988).
33.	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]
34.	P. Yeh and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons. Inc., New York, 1999).

OCIS Codes
(130.2790) Integrated optics : Guided waves
(230.3720) Optical devices : Liquid-crystal devices
(230.7380) Optical devices : Waveguides, channeled

ToC Category:
Integrated Optics

History
Original Manuscript: December 7, 2010
Revised Manuscript: January 25, 2011
Manuscript Accepted: January 28, 2011
Published: February 4, 2011

Citation
Chia-Chien Huang, "Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method," Opt. Express 19, 3363-3378 (2011)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-4-3363

Sort: Author | Year | Journal | Reset

References

J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D lateral light propagation in nematic-liquid-crystal cells with tilted molecules and nonlinear reorientational effect,” Opt. Quantum Electron. 37(1-3), 95–106 (2005). [CrossRef]
A. D’Alessandro, B. D. Donisi, R. Beccherelli, and R. Asquini, “Nematic liquid crystal optical channel waveguides on silicon,” IEEE J. Quantum Electron. 42(10), 1084–1090 (2006). [CrossRef]
G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finite-element beam propagation method,” Opt. Quantum Electron. 40(10), 733–748 (2008). [CrossRef]
E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707–5714 (2000). [CrossRef]
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]
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]
P. J. M. 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]
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]
P. J. M. 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]
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]
D. Donisi, B. Bellini, R. Beccherelli, R. Asquini, G. Gilardi, M. Trotta, and A. Dálessandro, “A switchable liquid-crystal optical channel waveguide on silicon,” IEEE J. Quantum Electron. 46(5), 762–768 (2010). [CrossRef]
M. Kawachi, N. Shibata, and T. Edahiro, “Possibility of use of liquid crystals as optical waveguide material for 1.3μm and 1.55μm bands,” Jpn. J. Appl. Phys. 21(Part 2, No. 3), L162–L164 (1982). [CrossRef]
M. Koshiba, K. Hayata, and M. Suzuki, “Finite element solution of anisotropic waveguides with arbitrary tensor permittivity,” J. Lightwave Technol. 4(2), 121–126 (1986). [CrossRef]
G. Tartarini and H. Renner, “Efficient finite-element analysis of tilted open anisotropic optical channel waveguides,”, ” IEEE Microw. Guid. Wave Lett. 9(10), 389–391 (1999). [CrossRef]
S. M. Hsu and H. C. Chang, “Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy,” Opt. Express 15(24), 15797–15811 (2007). [CrossRef] [PubMed]
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]
K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19(3), 405–413 (2001). [CrossRef]
J. C. Chen and S. Jüngling, “Computation of high-order waveguide modes by imaginary-distance beam propagation method,” Opt. Quantum Electron. 26, 199–205 (1994). [CrossRef]
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]
J. P. Boyd, Chebyshev and Fourier Spectral Methods (Springer–Verlag, 2nd edition, 2001).
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]
P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75(2 Pt 2), 026703 (2007). [CrossRef] [PubMed]
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]
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]
C. C. Huang, “Improved pseudospectral mode solver by prolate spheroidal wave functions for optical waveguides with step-index,” J. Lightwave Technol. 27(5), 597–605 (2009). [CrossRef]
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]
C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). [CrossRef] [PubMed]
T. Scharf, Polarized Light in Liquid Crystals and Polymers (John Wiley and Sons. Inc., NJ, 2007, Chap 4).
W. J. Gordon and C. A. Hall, “Transfinite element methods: blending function interpolation over arbitrary curved element domains,” Numer. Math. 21(2), 109–129 (1973). [CrossRef]
D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers (Springer–Verlag, 2009).
T. Tang, “The Hermite spectral method for Gauss-type functions,” SIAM J. Sci. Comput. 14(3), 594–605 (1993). [CrossRef]
T. Tamir, Guides-Wave Optoelectronics (Springer-Verlag, 1988).
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]
P. Yeh and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons. Inc., New York, 1999).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper