## Abstract

I extend the full vector pseudospectral-based eigenvalue scheme, based on the transverse magnetic field components, to analyze the mode behaviors of dielectric optical waveguides with transverse, nondiagonal anisotropy. One of the principal axes of the anisotropic materials is thus constrained to point in the longitudinal direction of the waveguide. I expand the guided mode fields in the interior subdomains with finite extent by using Chebyshev polynomials and those in the exterior subdomains with semi-infinite extent by using Laguerre–Gaussian functions with an accurately determined scaling factor. This study analyzes two examples: (1) the circularly-polarized modes of a magneto-optical raised strip waveguide and (2) the guided mode patterns of a nematic liquid-crystal channel waveguide under different orientations of the liquid-crystal molecule. The comparison of the numerical results with those from the vector finite difference approach demonstrates that my numerical approach has a higher computational efficiency and requires less computer memory.

©2010 Optical Society of America

## 1. Introduction

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

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

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

## 2. Mathematical formulations

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

*ω*is the angular frequency and

*μ*

_{0}is the permeability in vacuum. This study assumes that one of the material principal axes is aligned in the waveguide direction; thus, the permittivity tensor [

**ε**] yields the transverse, nondiagonal anisotropic form as

*ε*

_{0}is the permittivity in vacuum and [

*ε*] is the relative permittivity tensor. For a mode problem, all electromagnetic fields are assumed to have a

_{r}*z*dependence of exp(−

*jßz*), where

*β*=

*k*

_{0}

*n*, the wave number in vacuum is

_{eff}*k*

_{0}=

*ω*

^{2}

*μ*

_{0}

*ε*

_{0}, and

*n*is the effective refractive index. Using the divergence relation$\nabla \cdot H=0$, the magnetic field component

_{eff}*H*in Eq. (1) can be eliminated to yield the full vector eigenvalue equation presented by the transverse magnetic field components

_{z}*H*and

_{x}*H*as follows:

_{y}*P*,

_{xx}*P*,

_{xy}*P*, and

_{yx}*P*are defined as

_{yy}*H*and

_{x}*H*components for an anisotropic material with nondiagonal terms of [

_{y}**ε**], the coupling is determined from both the wave equations (

*P*and

_{xy}*P*terms) and the physical interface conditions. In the proposed approach, the computational domain is first divided at the interfaces between different materials into several subdomains with uniform or continuous refractive index profiles, and the whole computational domain is then assembled from the subdomains while satisfying the physical interface conditions. These interface conditions include the continuous normal and tangential components of magnetic fields at each intra-element boundary, and the continuities of the longitudinal components

_{yx}*H*and

_{z}*E*from the relations$\nabla \times H=$

_{z}*jωε*

_{0}

*n*

^{2}(

*x*,

*y*)

**E**and$\nabla \cdot H=0$, respectively. For a mode problem, all electromagnetic fields are assumed to have a

*z*dependence of exp(−

*jßz*), thus the continuity of the longitudinal component

*H*yields

_{z}*E*yields

_{z}*et al*. [20] introduced the curvilinear mapping technique to the pseudospectral scheme.

## 3. Numerical approach

In the pseudospectral approach [16,24], the computational domain is divided into several subdomains with uniform or continuous refractive index profiles. For an arbitrary interior rectangle subdomain *r*, it involves (*n _{x}* + 1) × (

*n*+ 1) collocation points in the corresponding (

_{y}*n*+ 1) × (

_{x}*n*+ 1) field unknowns as shown in Fig. 1(a) (where

_{y}*n*=

_{x}*n*= 10 is illustrated here).

_{y}The transverse magnetic field components to be solved in this subdomain are thus expanded as a product of suitable basis functions comprising the Lagrange-type interpolation *θ*(*x*) in the *x*-direction, *ψ*(*y*) in the *y*-direction, and the *H _{x}* and

*H*field values denoted as

_{y}*H*and

_{x,pq}*H*at the (

_{y,pq}*n*+ 1) × (

_{x}*n*+ 1) collocation points as follows:

_{y}*δ*denotes the Kronecker delta. Substituting the Eqs. (8a) and (8b) into Eq. (3), and then the pseudospectral scheme requires that the Eq. (3) has to be satisfied exactly at the specific (

_{mp}*n*+ 1) × (

_{x}*n*+ 1) collocation points that depend on the selected basis functions. The wave equations are thus converted to a system of linear equations. Note that the boundary collocation points (

_{y}*i*.

*e*., (

*x*

_{0},

*y*), (

_{j}*x*

_{10},

*y*), (

_{j}*x*,

_{i}*y*

_{0}), and (

*x*,

_{i}*y*

_{10}) as shown in Fig. 1(a), where

*i*= 0 to 10 and

*j*= 0 to 10) of the rectangle subdomain

*r*are shared by the adjacent subdomains, and the linear equations for the boundary collocation points have to be replaced by the physical interface conditions. These physical interface conditions are the continuous normal and tangential components of magnetic fields, and Eqs. (6) and (7). Removing the linear equations at these boundary collocation points, the matrix size is reduced to 2×(

*n*−1)×(

_{x}*n*− 1), and a matrix eigenvalue equation for the subdomain

_{y}*r*is given by

*P*,

_{xx}*P*,

_{xy}*P*, and

_{yx}*P*are

_{yy}*h-*th order derivatives of

*θ*(

_{p}*x*) to

*x*and

*ψ*(

_{q}*y*) to

*y*, respectively. After obtaining the matrix eigenvalue equation for the subdomain

*r*, a global matrix equation can be obtained by assembling all the matrix equations. Assuming there are total

*m*subdomains in whole computational window, the pattern of the matrix elements forms the following block diagonal matrix equation:

*P*includes the field unknowns

_{xx}H_{x}*H*

_{x}_{,}

*at interface boundary points (where the subscripts of summations are from*

_{pq}*p*= 0 to

*n*and

_{x}*q*= 0 to

*n*), although the linear equations satisfied at interface boundary points have been removed (where the subscripts of summations are only from

_{y}*i*= 1 to

*n*-1 and

_{x}*j*= 1 to

*n*-1). Therefore, these interface field unknowns ({

_{y}*H*

_{x}_{,0,0},

*H*

_{x}_{,0,1,…}

*H*

_{x}_{,0,}

*},{*

_{ny}*H*

_{x}_{,}

_{nx}_{,0},

*H*

_{x}_{,}

_{nx}_{,1,…}

*H*

_{x}_{,}

_{nx}_{,}

*},{*

_{ny}*H*

_{x}_{,0,0},

*H*

_{x}_{,1,0,…}

*H*

_{x}_{,}

_{nx}_{,0}},{

*H*

_{x}_{,0,}

*,*

_{ny}*H*

_{x}_{,1,}

_{ny}_{,…}

*H*

_{x}_{,}

_{nx}_{,}

*}) at interface boundary points involved in Eq. (11) have to be replaced by the field unknowns at all collocation points other than the interface boundary points by using the physical boundary conditions. The similar results are also found in other elements*

_{ny}*P*,

_{xy}H_{y}*P*, and

_{yx}H_{x}*P*.

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

*H*in Eq. (6) gives

_{z}*E*in Eq. (7) gives

_{z}*i*= 0, 1, 2,…

*n*. In Eqs. (13) and (15), the ${H}_{x}^{(1)}$ and ${H}_{y}^{(1)}$denote the magnetic field components of

_{y}**H**

_{x}and

**H**

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

We now determine the appropriate basis functions by expanding the *H _{x}* and

*H*components. The reports [18,19] show that the Chebyshev polynomials are better for expanding the optical fields in interior subdomains because of their mathematical robustness to non-periodic structures. The other choice to expand the field profiles in the interior regions is Legendre polynomials. The comparison of the convergence between Chebyshev and Legendre basis functions for analyzing a circular waveguide was made by Chiang

_{y}*et al*. [20]. They calculated the relative errors in effective index relative to the exact values of HE

_{11}mode. The results showed that Chebyshev polynomials achieve the better accuracy and faster convergence than that based on Legendre polynomials. This result is resulted from the different mesh division profiles (

*i.e*., collocation point profiles) between Chebyshev and Legendre polynomials [16,24]. In contrast to the exterior subdomains, the LGFs with exponential decay characteristics match the guided field profiles with exponential decay. In general, there are three ways to expand the unknown fields in a semi-infinite interval [0,∞): (1) expand in LGFs, (2) map the semi-infinite extent into a finite one and then expand in Chebyshev polynomials, and (3) truncate the semi-infinite domain to finite one [0,

*x*

_{max}] and use a Chebyshev expansion [16]. For the third choice, only one way by increasing

*x*

_{max}, as the number of terms of basis functions can be used to achieve the exponential convergence. Therefore, the computational efforts in CPU time and memory storage are much heavy. The second approach has several types of map including the algebraic, logarithmic, and exponential maps to be considered. In addition, these maps also need to efficiently determine their optimal scale factors. In contrast, LGFs with exponential decay characteristics (in mathematical) efficiently match the exponential field profiles of guided modes (in physical) on waveguide problems, and the scaling factor in LGFs has been definitely determined by the report for the isotropic media [19]. Besides, Boyd [24] has noticed that uniform spectral accuracy of LGFs requires at least exponential decay of the function as

*x→*∞, the feature is similar to the characteristic of the general guided field profiles. If Chebyshev polynomials are used to represent the guided mode profiles, more terms of basis functions than that of LGFs are required. In mathematics, various basis functions can be converted to their Lagrange-type form. For Chebyshev polynomials, the explicit form of the Lagrange-type interpolation function is as follows [24]:

*T*(

_{v}*x*) is the general Chebyshev polynomial of order

*v*and

*x*denotes the collocation points for the Chebyshev polynomials. For the LGFs, the explicit form is as follows [24]:

_{p}*L*(

_{v}*αx*) is the Laguerre polynomial of order

*v*and

*x*denotes the collocation points for the LGFs. In Eq. (18), the scaling factor

_{p}*α*significantly affects the numerical accuracy of results obtained with the LGFs. The definite value of

*α*for isotropic materials has been derived in previous studies [19]; hence, results obtained with LGFs are sufficiently accurate.

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

*M*is the spreading of the guided fields and

*x*is the collocation point at position

_{i}*i*. The choice of

*α*is thus converted to determine the value of

*M*. Using EIM, the effective index and the decay rate of the desired guided mode has to be obtained in advance, and then

*M*can be estimated subsequently. For finding out the effective index and the decay rate, a rib waveguide shown in Fig. 2(a) is used to describe the determination of

*M*, where, [

**ε**

_{c}], [

**ε**

_{s}], and [

**ε**

_{a}] denote the permittivity tensors of core, substrate, and cover, respectively.

The first process is to divide the rib waveguide by the vertical dotted lines into three subregions labeled as 1, 2, and 3, and then the dispersion relation [27] of three-layer waveguide along the *y*-direction in each subregion computes the effective indices of subregions 1, 2 and 3 as *n _{y}*

_{1},

*n*

_{y}_{2}and

*n*

_{y}_{3}, respectively, as shown in Fig. 2(b). Importantly, different from the condition of isotropic media [19], the refractive index for each material to be used in the dispersion relation [27] is (

*ε*)

_{xx}^{1/2}in this study. The dispersion relation is utilized again in the

*x*-direction to obtain the equivalent effective index

*n*to the original waveguide structure. The decay rates

_{eq}*γ*along the

*y*-direction extending to the substrate for the three subregions are given by

*y*-direction extending to air areand the decay rates along the

*x*-direction arewhere

*n*is the effective refractive index for the original waveguide. Accordingly, the spreadings of the magnetic fields along the

_{eq}*y*-direction extending to the substrate

*M*are given

_{ysi}*y*-direction extending to the air

*M*are

_{yai}*x*-direction

*M*are

_{xi}*x*and

_{t}*y*are the turning points in

_{t}*x*- and

*y*-directions, respectively. In a step-index waveguide, the turning points are located at the interfaces between different materials. As shown in Eqs. (21a)–(21c), the parameter

*κ*denotes the ratio of the field amplitude at the positions of turning points plus the spreadings

*M*relative to that at turning points. For simplicity, we let

*κ*be equal in both the

*x*- and the

*y*-directions. In the report [19], the ratio of

*κ*= 10

^{−2}is sufficient to involve the significant contribution of the guided field profiles.

## 4. Simulation results and discussion

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

#### 4.1 A magneto-optical raised strip waveguide

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

*n*= 1.

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

*n*= 2.302 and

_{f}*ζ*is the first-order magneto-optic effect and proportional to static magnetization. The width

*w*of the core is 800 nm, and its thickness

*t*is set to be 607.6 nm, which was carefully adjusted to achieve the nearly degenerate mode between the first and the second order modes when

*ζ*= 0 [8,28], for the operating wavelength

*λ*= 1300 nm in free space.

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

For examining the convergence of the proposed scheme, I calculated the results for the condition *ζ* = 0.005. Note that the terms of the basis functions for expanding the exterior subdomains were fixed at the number *N _{ex}* = 10. Table 1
shows the convergence of effective refractive indices of the first (labeled as${n}_{eff}^{(1)}$) and the second (labeled as${n}_{eff}^{(2)}$) order modes and the computational time versus the term of the basis functions

*N*expanding the interior subdomains. Throughout this study, equal terms of the basis function are used in both the

_{int}*x*- and

*y*-directions.The calculations were performed on a personal computer with an Intel core i7 CPU 920 at 2.67 GHz. The number of grid points in this example was (

*N*+ 2

_{int}*N*− 2) × (

_{ex}*N*+ 2

_{int}*N*− 2), and the computational window was approximately 3.2 μm × 3 μm. The difference of the calculated effective index between

_{ex}*N*= 10 and

_{int}*N*= 25 was in the order of 10

_{int}^{−4}; by increasing the terms of basis functions to

*N*= 20, this difference was reduced to the order of 10

_{int}^{−5}. Thus the proposed scheme obtains a sufficiently accurate result (it achieved accuracy of the order of 10

^{−4}using only 28 × 28 grid points). The growth of the computational time listed in Table 1 was fairly slow due to the use of LGFs for exterior subdomains. If Chebyshev polynomials are used to expand the field profiles for all subdomains, the computational times were 9.9 s, 44.8 s, and 281.9s while using 15, 20, and 25 terms, respectively. Therefore, the proposed scheme is more efficient than that using the unified set of basis functions all over the whole domain. For different values of

*ζ*, Table 2 shows the results obtained by the proposed scheme using

*N*= 25 and

_{int}*N*= 10 (the grid points were now 45 × 45) and those obtained with the vector FD method [8]. The results obtained from the two schemes are in good agreement. In the vector FD method [8], the grid points were 320 × 348 over a computational window of 3.2 μm × 2.9 μm. This comparison demonstrates that the computational efficiency of the proposed method is superior to that of the vector FD method. For the condition

_{ex}*ζ*= 0.005, Fig. 4 shows the |

*H*| and

_{x}*|H*| field profiles of the first order mode, and Fig. 5 shows the |

_{y}*H*| and

_{x}*|H*| field profiles of the second order mode. The magnitudes and profiles of the |

_{y}*H*| and

_{x}*|H*| fields are nearly equal for the two modes.

_{y}To clearly represent the polarization states of the two modes, Figs. 6
and 7
plot the field profiles of |*H _{±}* |

*= |H*/2of the first and the second order modes, respectively. For the first order mode, the power of the |

_{x}± jH_{y}|*H*−

_{x}*jH*| /2field is suppressed by approximately 17 dB compared with the peak value of the |

_{y}*H*+

_{x}*jH*| /2field, a suppression amount equal to that calculated by [8]. Figure 6 indicates that the phase difference between the complex

_{y}*H*and

_{x}*H*components was approximately π/2, and the first order mode is in the left-hand circularly-polarized state. In contrast, Fig. 7 indicates that the second order mode is in the right-hand circularly-polarized state. In addition, the calculated Faraday rotation rate for the condition

_{y}*ζ*= 0.005 is ${\Theta}_{F}=(\pi /\lambda )({n}_{eff}^{(1)}-{n}_{eff}^{(2)})={2397}^{\text{o}}/\text{cm}$.

#### 4.2 A nematic liquid-crystal channel waveguide

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

The permittivity tensor of the core region is given by

*n*= 1.5292 and

_{o}*n*= 1.7072 are the ordinary and extraordinary refractive indices of the nematic LC, respectively, at the wavelength

_{e}*λ*= 1.55 μm. In Fig. 8(b),

*φ*denotes the twist angle of the LC between the

*x*-axis and director$\widehat{n}$. The voltage applied by the designed electrodes is applied to alter the value of

*φ*. By the proposed scheme, the division of the computational window is the same as the first example shown in Fig. 3(b), and all the calculations used

*N*= 10 and

_{ex}*N*= 20 in each subdomain. Table 3 shows the calculated effective indices at various values of

_{int}*φ*for the first seven modes.For the mode 1, the variation of effective index versus

*φ*is insignificant. The result emerges in the major field patterns as shown in Fig. 9 . The small differences of effective index only result in the different minor field patterns.

However, the variations of effective index versus *φ* are more significant for the higher order modes as shown in Fig. 10
. The symbols |*H _{m,x}*| and |

*H*| denote the modului of the

_{m,y}*H*and

_{x}*H*components, respectively, of the

_{y}*m*-th order guided mode. In Fig. 10, the symmetric field patterns corresponding to the

*y*-axis are observed for the conditions

*φ*= 0° and

*φ*= 90° because the elements of the permittivity tensor only appear in the principal axes. For other values of

*φ*, more complicated symmetries exist only in the lower order modes. For the conditions

*φ*= 0° and

*φ*= 90°, both of the minor powers corresponding to the peak values are suppressed by more than 40 dB, and those for

*φ*= 30° and

*φ*= 60° are suppressed by approximately 5 dB.

In contrast, the field modului of the *H _{x}* and

*H*components are nearly identical for

_{y}*φ*= 45° because the two polarizations experience identical propagation in the core region, except for the small asymmetries in the

*x*- and

*y*- directions. Hence, the calculated field patterns of |

*H*| and |

_{x}*H*| are not perfectly identical. For brevity, the patterns for

_{y}*φ*= 45° are not shown. In particular, the magnitudes of the major |

*H*| and minor |

_{x}*H*| fields are reversed for the seventh mode (not shown here) at these values of

_{y}*φ*. This reverse phenomenon appears only in some of the higher order modes.

## 5. Conclusion

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

## Acknowledgements

I thank the National Science Council of the Republic of China, Taiwan, for financially supporting this research under Contract No. NSC 99-2112-M-005-005-MY3.

## References and links

**1. **C. L. Xu, W. P. Huang, J. Chrostowski, and S. K. Chaudhuri, “A full-vectorial beam propagation method for anisotropic waveguides,” J. Lightwave Technol. **12**(11), 1926–1931 (1994). [CrossRef]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

**19. **C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron. **11**(2), 457–465 (2005). [CrossRef]

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

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

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

**23. **J. B. Xiao and X. H. Sun, “Full-vectorial mode solver for anisotropic optical waveguides using multidomain spectral collocation method,” Opt. Commun. **283**(14), 2835–2840 (2010). [CrossRef]

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

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

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

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

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

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

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

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

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

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