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

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 frequencydomain 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. ©2011 Optical Society of America OCIS codes: (130.2790) Guided waves; (230.7380) Waveguides, channeled; (230.3720) Liquid crystal devices. 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). 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). 3. G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finiteelement beam propagation method,” Opt. Quantum Electron. 40(10), 733–748 (2008). 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). 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). 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). 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). 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). 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). 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). 11. D. Donisi, B. Bellini, R. Beccherelli, R. Asquini, G. Gilardi, M. Trotta, and A. Dálessandro, “A switchable liquidcrystal optical channel waveguide on silicon,” IEEE J. Quantum Electron. 46(5), 762–768 (2010). #139244 $15.00 USD Received 7 Dec 2010; revised 25 Jan 2011; accepted 28 Jan 2011; published 4 Feb 2011 (C) 2011 OSA 14 February 2011 / Vol. 19, No. 4 / OPTICS EXPRESS 3363 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 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). 27. C. C. Huang, “Modeling mode characteristics of transverse anisotropic waveguides using a vector pseudospectral approach,” Opt. Express 18(25), 26583–26599 (2010). 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). 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). 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). 34. P. Yeh and C. Gu, Optics of Liquid Crystal Displays (John Wiley and Sons. Inc., New York, 1999).


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][2][3][4][5][6][7][8][9][10][11] because of their low scattering loss at near-infrared wavelengths [12] 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 fullwave techniques, a number of numerical schemes have been proposed.These approaches are the finite-element frequency-domain (FEFD)-based eigenvalue mode solvers [6,[13][14][15], the finite-element beam propagation methods (FE-BPMs) [3,9,16,17], the finite-difference beam propagation method (FD-BPM) [4], and the finite-difference frequency-domain (FDFD)-based eigenvalue mode solvers [10].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 2 , 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 2 .In the FD-BPM [4], a system of coupled differential equations, which involves an electric and a magnetic field component with 2N unknowns (the matrix size is 4N 2 ), is proposed for solving the twisted nematic pixels found in microdisplays.The required computational effort is relatively low.These BPM-based schemes [3,4,9,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,19].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][14][15] and FDFD [10] schemes, can deal with problems including all high-order modes.In [13][14][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], tilted channel waveguides [14], 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 2 ) 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] has demonstrated fast convergence and high-order accuracy computational performance for solving propagation characteristics of various dielectric waveguides and optical fibers [21][22][23][24][25][26][27].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][22][23][24][25] or at most the transverse anisotropic media [26,27].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 2 , 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].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.

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 12 0 where ω is the angular frequency and μ 0 is the permeability in vacuum.In this work, the permittivity tensor [ε] with the general form is considered as follows: where ε 0 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 0 n eff , the wave number in free space is k 0 = ω 2 μ 0 ε 0 , 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 0   H .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 β 2 , the full vector eigenvalue equation represented using the transverse magnetic field components H x and H y is obtained as follows: where the differential operators P xx (β), P xy (β), P yx (β), and P yy (β) are defined as a function of β by the following.
Compared with that in the transverse anisotropic media [27], 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].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 0 and that the E z is expressed by through the Ampere's law   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].The main difference of interface conditions between the full 3 × 3 anisotropic here and the transverse anisotropic media [27] is appeared in Eq. (7).Equation ( 7) involving the eigenvalue β makes the proposed work significantly complex while patching the subdomains than that in [27] 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] are the most often used approaches in various numerical schemes such as the FD, FE, and pseudospectral methods [30].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,23].

B. Numerical approach
In the pseudospectral scheme [20], 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).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: , 00 ( , ) ( ) ( ) ,  3) is demanded to perfectly be satisfied at these (n x  1)  (n y  1) interior collocation points (namely, (x 1 , y q ), (x 2 , y q ),… (x 9 , 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 β 2 given by where the operators ()  (1) (1) where () () and () () h q 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 where 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 0 , 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 , (1) 01 0 0 0 (1) 01 0 0 0 (1) 01 (1) 01 and the continuity of E z gives , where , (1) 0 0 0 , 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 0 , y j ) and (x nx , y j ) are altered by the horizontal ones (x i , y 0 ) 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] 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 semiinfinite 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]: 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]: 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] and the idea of the effective index method (EIM) [32].From the derivation [31], the optimum α is determined by the following formula: where M is the spreading of the guided fields to be solved, and 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].The detailed derivation of determining M has been shown in the previous work [27] 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 β 2 , 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 β 2 is obtained and solved using the Arnoldi iteration method [33], 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 8 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].

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.Assuming that the orientation of the square LC director n tilts at an angle of θ c = 45° with respect to the z-axis (the propagation direction) as shown in Fig. 2 where n o = 1.55, n e = 1.8, and 22 eo nn     .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 semiinfinite 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)] 2 .The differences of the effective indices between the fourth and fifth iterations are achieved in the order of 10 8 , 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 |.The following example considers an LC dielectric waveguide comprising a core region filled with nematic LCs (5CB) [34], 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).
where 22 eo nn     , 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 n and the z-axis, and φ c is the twist angle between the projection of 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 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,27].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.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

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 2 ), 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 2 ) and the finite-difference frequency-domain (FDFD)-based eigenvalue scheme with 4N unknowns (matrix size of 16N 2 ).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.

Fig. 1 .
Fig. 1.The mesh division of an arbitrary interior subdomain r

Fig. 2 .
Fig. 2. (a) Cross-sectional view of the square LC waveguide with the permittivity tensor [ε] in the core region and refractive index ns of the surrounding substrate.(b) Schematic diagram of the orientation n of the LC director.(c) Division of the computational domain for the LC waveguide.
(b), and thus, the permittivity tensor [ε] of the LC is given by

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

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

Fig. 5 .
Fig. 5. (a) Cross-sectional view of the nematic LC-core channel waveguide with the core permittivity tensor [ε], substrate refractive index ng, and air refractive index na.(b) Schematic diagram of the twist angle φc and the tilt angle θc of the LC director n .The permittivity tensor [ε] of the LC-core region at the operating wavelength λ = 1.55 μm is given by 2 2 2 2 2 2 2 2 0 22 sin cos sin sin cos sin cos cos [ ] sin sin cos sin sin sin cos sin , sin cos cos sin cos sin cos

Fig. 6 .
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 4 ) 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°.
(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.