A finite-difference frequency-domain method for full-vectroial mode solutions of anisotropic optical waveguides with arbitrary permittivity tensor

doi:10.1364/OE.17.005965

« Show journal navigation

A finite-difference frequency-domain method for full-vectroial mode solutions of anisotropic optical waveguides with arbitrary permittivity tensor

Ming-yun Chen, Sen-ming Hsu, and Hung-Chun Chang »View Author Affiliations

Optics Express, Vol. 17, Issue 8, pp. 5965-5979 (2009)
http://dx.doi.org/10.1364/OE.17.005965

View Full Text Article

Acrobat PDF (971 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. Formulation
3. Numerical examples
4. 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 (19)
Cited By
Figures (17)
Metrics

Abstract

A new finite-difference frequency-domain (FDFD) method based eigenvalue algorithm is developed for analyzing anisotropic optical waveguides with an arbitrary permittivity tensor. Yee’s mesh is employed in the FD formulation along with perfectly matched layer (PML) absorption boundary conditions. A standard eigenvalue matrix equation is successfully derived through considering simultaneously four transverse field components. The new algorithm is first applied to the mode solution of a proton-exchanged LiNbO₃ optical waveguide and the results agree with those obtained using a full-vectorial finite-element beam propagation method. Then, the algorithm is used to study modes on a liquid-crystal optical waveguide with arbitrary molecular director orientation. This arbitrary orientation may cause the loss of transverse-axis symmetries of the waveguide with symmetric background structure. Asymmetric mode-field profiles under such situations are clearly demonstrated in the numerical examples.

1. Introduction

Solving full-vector modes of optical waveguides has been an essential task in developing and designing state-of-the-art waveguide devices and lightwave circuits. Among different numerical solution methods, the finite-difference (FD) method is attractive due to its advantage of simple formulation and numerical implementation. The FD algorithm can be formulated either through the wave equations [1

C. L. Xu, W. P. Huang, M. S. Stern, and S. K. Chaudhuri, “Full-vectorial mode calculations by finite difference method,” Proc. Inst. Electr. Eng. 141, 281–286 (1994). [CrossRef]

] or the Helmholtz equation [2

P. Lüsse, P. Stuwe, J. Schule, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–494 (1994). [CrossRef]

]. Also high-accuracy algorithms have been pursued [3–6

G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis I: Uniform regions and dielectric interfaces,” J. Lightwave Technol. 20, 1210–1218 (2002). [CrossRef]

]. Besides the traditional approach in which either the electric field components or the magnetic field components are considered at grid points, such as in the above-mentioned algorithms, schemes based on the Yee’s mesh, as in the finite-difference time-domain (FDTD) method [7

K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. AP-14, 302–307 (1966).

], and derived directly from Maxwell’s equations have been formulated [8–10

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, 1627–1634 (2002). [CrossRef]

]. Such schemes have the advantage that the obtained mode fields can easily be cast into the FDTD computation. We name such Yee-mesh-based FD approach the finite-difference frequency-domain (FDFD) method.

Relatively fewer FD works deal with structures including anisotropic materials [11–14

L. Thylen and D. Yevick, “Beam propagation method in anisotropic media,” Appl. Opt. 21, 2751–2754 (1982). [CrossRef] [PubMed]

] and they at most considered anisotropy in the transverse plane, i.e., in the waveguide cross-section. Most recently, Fallahkhair et al. [15

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

] developed a more rigorous vector FD mode solver based on the traditional approach for anisotropic waveguides, but again assuming non-diagonal anisotropy only in the transverse plane. In this paper, we present an FDFD method for solving full-vector modes of optical waveguides with arbitrary permittivity tensor, i.e., with general three-dimensional (3-D) anisotropy. In the finite-element formulations, mode solutions for such general waveguides have been reported by Saitoh and Koshiba [16

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

] but based on the beam propagation approach. Here, we derive a matrix standard eigenvalue problem from the FDFD method, which is easier to solve for the modes. As in [16

], we incorporate perfectly matched layers (PMLs) for anisotropic media [17

F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media”, IEEE Microwave Guid. Wave Lett. 8, 223V–225 (1998). [CrossRef]

] into our formulation as the absorbing boundary condition at the outer boundaries of the computational domain. In so doing, leaky waveguide modes having complex propagation constants can also be analyzed.

After presenting the derivation of the new eigenmode solution algorithm in Section 2, numerical examples are presented in Section 3. The algorithm is first applied to solve a proton exchanged LiNbO₃ optical waveguide that has been analyzed using a full-vectorial finite-element beam propagation method [16

]. Then, modes on a liquid-crystal (LC) optical waveguide with arbitrary molecular director orientation are obtained and discussed. This arbitrary orientation may cause the loss of transverse-axis symmetries of the waveguide with symmetric background structure. Asymmetric mode-field profiles under such situations are clearly demonstrated. The conclusion is drawn in Section 4.

Fig. 1. Cross-section of a general optical waveguide and the computational domain with PML regions.

Fig. 2. Two-dimensional Yee’s mesh for the FDFD method.

2. Formulation

Figure 1 shows the computational domain for the analysis of a dielectric waveguide where the waveguide cross-section in the transverse x-y plane is truncated and surrounded by PML regions of thickness d. The incorporation of PML regions allows the analysis of leaky modes [10

C. P. Yu and H. C. Chang, “Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Opt. Express , 12, 6165V–6177 (2004). http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-25-6165. [CrossRef]

]. For frequency-domain mode solutions, Maxwell’s curl equations are written as

\nabla \times \vec{E} = - jω [μ] \vec{H}

(1)

\nabla \times \vec{H} = jω \vec{D}

(2)

with the assumption that the electromagnetic fields have the time dependence of exp (jωt). Here, ω is the angular frequency and [μ] represents a diagonal permeability tensor with diagonal elements μ_x , μ_y , and μ_z . We are mainly interested in non-magnetic waveguide structures with μ_x = μ_y = μ_z = μ ₀, where μ ₀ is the permeability of free space. But in the PML regions, this more general permeability tensor needs to be adopted. Then, after the assumption of the exp(-jβz) dependence of the fields propagating along the z direction, with β being the modal propagation constant, Eqs. (1) and (2) can be written into the following six component equations in terms of the six electric and magnetic field components:

- jω μ_{x} H_{x} = \frac{\partial E_{z}}{\partial y} + jβ E_{y}

(3)

- jω μ_{y} H_{y} = - jβ E_{x} - \frac{\partial E_{z}}{\partial y}

(4)

- jω μ_{z} H_{z} = \frac{\partial E_{y}}{\partial x} + \frac{\partial E_{x}}{\partial y}

(5)

jω D_{x} = \frac{\partial H_{z}}{\partial y} jβ H_{y}

(6)

jω D_{y} = - jβ H_{x} - \frac{\partial H_{z}}{\partial x}

(7)

jω D_{z} = \frac{\partial H_{y}}{\partial x} - \frac{\partial H_{x}}{\partial y} .

(8)

Applying the Yee’s mesh, as depicted in Fig. 2, and the central difference scheme, Eqs (3)–(8) become

- jω μ_{x} H_{x, (i, j + \frac{1}{2})} = \frac{E_{z, (i, j + 1)} - E_{z, (i, j)}}{Δ y} + jβ E_{y, (i, j + \frac{1}{2})}

(9)

- jω μ_{y} H_{y, (i + \frac{1}{2}, j)} = - jβ E_{x, (i + \frac{1}{2}, j)} - \frac{E_{z, (i + 1, j)} - E_{z, (i, j)}}{Δ x}

(10)

- jω μ_{z} H_{z, (i + \frac{1}{2}, j + \frac{1}{2})} = \frac{E_{y, (i + 1, j + \frac{1}{2})} - E_{y, (i, j - \frac{1}{2})}}{Δ x} + \frac{E_{x, (i + \frac{1}{2}, j + 1)} - E_{x, (i, + \frac{1}{2}, j)}}{Δ y}

(11)

jω D_{x, (i + \frac{1}{2}, j)} = \frac{H_{z, (i + \frac{1}{2}, j + \frac{1}{2})} - H_{z, (i + \frac{1}{2}, j - \frac{1}{2})}}{Δ y} + jβ H_{y, (i + \frac{1}{2}, j)}

(12)

jω D_{y, (i, j + \frac{1}{2})} = - jβ H_{x, (i, j + \frac{1}{2})} - \frac{H_{z, (i + \frac{1}{2}, j + \frac{1}{2})} - H_{z, (i - \frac{1}{2}, j + \frac{1}{2})}}{Δ x}

(13)

jω D_{z, (i, j)} = \frac{H_{y, (i + \frac{1}{2}, j)} - H_{y, (i - \frac{1}{2}, j)}}{Δ x} - \frac{H_{x, (i, j + \frac{1}{2})} - H_{x, (i, j + \frac{1}{2})}}{Δ y}

(14)

which can be written in the matrix form as

- jω μ_{0} [\begin{matrix} μ_{x} H_{x} \\ μ_{y} H_{y} \\ μ_{z} H_{z} \end{matrix}] = [\begin{matrix} 0 & jβ I & U_{y} \\ - jβ I & 0 & - U_{x} \\ - U_{y} & U_{x} & 0 \end{matrix}] [\begin{matrix} E_{x} \\ E_{y} \\ E_{z} \end{matrix}]

(15)

jω [\begin{matrix} D_{x} \\ D_{y} \\ D_{z} \end{matrix}] = [\begin{matrix} 0 & jβ I & V_{y} \\ - jβ I & 0 & - V_{x} \\ - V_{y} & V_{x} & 0 \end{matrix}] [\begin{matrix} H_{x} \\ H_{y} \\ H_{z} \end{matrix}] .

(16)

In Eqs. (15)–(16), I is the identity matrix, U_x , U_y , V_x and V_y are the matrices given as below

U_{x} = \frac{1}{Δ x} [\begin{matrix} - 1 & 1 \\ - 1 & 1 \\ ⋱ & ⋱ \\ ⋱ & ⋱ \\ - 1 & 1 \\ - 1 \end{matrix}]

(17)

U_{y} = \frac{1}{Δy} [\begin{matrix} - 1 & 1 \\ - 1 & ⋱ \\ ⋱ & 1 \\ ⋱ \\ - 1 \\ - 1 \end{matrix}]

(18)

V_{x} = \frac{1}{Δ x} [\begin{matrix} 1 \\ - 1 & 1 \\ - 1 & ⋱ \\ ⋱ & ⋱ \\ - 1 & 1 \\ - 1 & 1 \end{matrix}]

(19)

V_{y} = \frac{1}{Δ y} [\begin{matrix} 1 \\ 1 \\ ⋱ \\ - 1 & ⋱ \\ ⋱ & 1 \\ - 1 & 1 \end{matrix}]

(20)

and each of E_x , E_y , E_z , H_x , H_y , and H_z represents a column vector with its elements being the corresponding Yee-mesh field components over the whole computational domain. Please notice that our notations for Yee-mesh arrangement in Fig. 2 and in Eqs. (9)–(20) are those used in [10

], as adopted from what popularly utilized in the FDTD method, and a little different from those presented in [9

Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10, 853–864 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853. [PubMed]

]. Now, for the general anisotropic dielectric medium, the relation between E⃗ and D⃗ is expressed at each spatial point by

D_{x} = ε_{xx} E_{x} + ε_{xy} E_{y} + ε_{xz} E_{z}

(21)

D_{y} = ε_{yx} E_{x} + ε_{yy} E_{y} + ε_{yz} E_{z}

(22)

D_{z} = ε_{zx} E_{x} + ε_{zy} E_{y} + ε_{zz} E_{z}

(23)

which define the nine ε_ij (i, j = x, y, z) elements for the permittivity tensor at the given point. By including Eqs. (21)–(23) into the formulation of Eqs. (15) and (16), we can obtain

- jω μ_{x} H_{x} = jβ E_{y} + U_{y} E_{z}

(24)

- jω μ_{y} H_{y} = - jβ E_{x} + U_{x} E_{z}

(25)

- jω μ_{z} H_{z} = - U_{y} E_{x} + U_{x} E_{y}

(26)

jω (ε_{xx} E_{x} + ε_{xy} E_{y} + ε_{xz} E_{z}) = jβ H_{y} + V_{y} H_{z}

(27)

jω (ε_{yx} E_{x} + ε_{yy} E_{y} + ε_{yz} E_{z}) = - j β H_{x} - V_{x} H_{z}

(28)

jω (ε_{zx} E_{x} + ε_{zy} E_{y} + ε_{zz} E_{z}) = - V_{y} H_{x} + V_{x} H_{y} .

(29)

Note that in Eqs. (27)–(29), each ε_ijE_j (i, j = x, y, z) should be considered as a column vector with each element representing the value of ε_ijE_j at a corresponding grid point. Please also note that at and near the grid point (i, j), referring to Fig. 2, the values of ε_xj ’s, ε_yj ’s, and ε_zj ’s are assigned according to the material properties at the grid points (i + 1/2, j), (i, j + 1/2), and (i, j), respectively, which are the locations of D_x , D_y , and D_z , respectively. By substituting the H_z and E_z expressions from Eqs. (26) and (29), respectively, into Eqs. (24), (25), (27), and (28), an eigenvalue matrix equation for the four transverse field components can be derived as in the following form:

β [\begin{matrix} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{matrix}] = [\begin{matrix} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \\ A_{41} & A_{42} & A_{43} & A_{44} \end{matrix}] [\begin{matrix} E_{x} \\ E_{y} \\ H_{x} \\ H_{y} \end{matrix}]

(30)

where the sub-matrices A_ij (i, j = 1,⋯ ,4) are defined as below:

A_{11} = - j (\frac{ε_{zx}}{ε_{zz}}) U_{x}

(31)

A_{12} = - j (\frac{ε_{zy}}{ε_{zz}}) U_{x}

(32)

A_{13} = - (\frac{1}{{ωε}_{zz}}) U_{x} V_{y}

(33)

A_{14} = - (\frac{1}{{ωε}_{zz}}) U_{x} V_{y} + ω μ_{y} I

(34)

A_{21} = - j (\frac{ε_{zx}}{ε_{zz}}) U_{y}

(35)

A_{22} = - j (\frac{ε_{zy}}{ε_{zz}}) U_{y}

(36)

A_{23} = - (\frac{1}{{ωε}_{zz}}) U_{y} V_{y} - ω μ_{x} I

(37)

A_{24} = - (\frac{1}{{ωε}_{zz}}) U_{y} V_{x}

(38)

A_{31} = - ω ε_{yx} I + (\frac{1}{ω μ_{z}}) V_{x} U_{y} + \frac{ε_{yz} ω ε_{zx}}{ε_{zz}} I

(39)

A_{32} = - ω ε_{yy} I - (\frac{1}{ω μ_{z}}) V_{x} U_{x} + \frac{ε_{yz} ω ε_{zy}}{ε_{zz}} I

(40)

A_{33} = - j (\frac{ε_{yz}}{ε_{zz}}) V_{y}

(41)

A_{34} = j (\frac{ε_{yz}}{ε_{zz}}) V_{x}

(42)

A_{41} = ω ε_{xx} I + (\frac{1}{ω μ_{z}}) V_{y} U_{y} - \frac{ε_{xz} ω ε_{zx}}{ε_{zz}} I

(43)

A_{42} = ω ε_{xy} I - (\frac{1}{ω μ_{z}}) V_{y} U_{x} - \frac{ε_{xz} ω ε_{zy}}{ε_{zz}} I

(44)

A_{43} = j (\frac{ε_{xz}}{ε_{zz}}) V_{y}

(45)

A_{44} = - j (\frac{ε_{xz}}{ε_{zz}}) V_{x} .

(46)

For the considered problem involving non-diagonal permittivity tensor, we adopt the more general PML for anisotropic media proposed by Teixeira and Chew [17

] as the absorbing boundary condition, as in the formulation in [16

]. Therefore, in the PML regions, the permittivity and permeability tensors are taken to be

[ε_{PML}] = [\begin{matrix} \frac{s_{y} s_{z}}{s_{x}} ε_{xx} & s_{z} ε_{xy} & s_{y} ε_{xz} \\ s_{z} ε_{yx} & \frac{s_{z} s_{x}}{s_{t}} ε_{yy} & s_{x} ε_{yz} \\ s_{y} ε_{zx} & s_{x} ε_{zy} & \frac{s_{x} s_{y}}{s_{z}} ε_{zz} \end{matrix}]

(47)

[μ_{PML}] = μ_{0} [\begin{matrix} \frac{s_{y} s_{z}}{s_{x}} & 0 & 0 \\ 0 & \frac{s_{z} s_{x}}{s_{y}} & 0 \\ 0 & 0 & \frac{s_{x} s_{y}}{s_{z}} \end{matrix}]

(48)

where s_x , s_y , and s_z are the complex PML parameters defined as

s_{j} = 1 - j α_{j}

(49)

with j = x,y,z and α_j ’s are assigned appropriate values to control the field attenuation in PML regions. In our application, we choose

α_{j} = α_{j, max} {(\frac{ρ}{d})}^{2}

(50)

for j = x and y, and α_z = 0, where ρ represents the distance in the j-direction from the beginning of the PML region and α _j,max is determined by the assumed reflectivity value from the PML layer [18

R. Mittra and U. Pekel, “A new look an the perfectly matched layer PML concept for the reflectionless absorption of electromagnetic waves,” IEEE Microwave Guid. Wave Lett. 5, 84–86 (1995). [CrossRef]

3. Numerical examples

To demonstrate the accuracy and applications of the proposed FDFD mode solver algorithm for waveguides with arbitrary permittivity tensor, we present two numerical examples, one with LiNbO₃ (LN) materials and the other involving LCs with arbitrary director orientation.

We first analyze an LN waveguide considered in [16

] with the FE-BPM method. Figure 3 shows its cross-section with surrounded PMLs, where the substrate is that of LN and the core region is made of proton-exchanged LN (PE-LN). The structure and computational-domain parameters are t = w = 3μm, W_x = 9μm, W_y = 7 μm, d = 2μm, and the thickness of the air-cladding region up to the PML region is 1μm. These parameter values are the same as those assumed in [16

]. The ordinary and extraordinary refractive indices of the substrate are n_o = 2.250 and n_e = 2.172, respectively. The refractive index difference between the PE-LN core region and the substrate is assumed to be Δn_e = 0.01 for the extraordinary index and Δn_o = 0 for the ordinary index. The relative permittivity tensor elements for LN are given by

ε_{xx} / ε_{0} = n_{o}^{2}

(51)

ε_{yy} / ε_{0} = n_{e}^{2} \cos^{2} θ_{c} + n_{o}^{2} \sin^{2} θ_{c}

(52)

ε_{zz} / ε_{0} = n_{o}^{2} \cos^{2} θ_{c} + n_{e}^{2} \sin^{2} θ_{c}

(53)

ε_{yz} / ε_{0} = ε_{zy} = (n_{o}^{2} - n_{e}^{2}) \cos θ_{c} \sin θ_{c}

(54)

ε_{xy} / ε_{0} = ε_{yx} = ε_{zx} = ε_{xz} = 0

(55)

where ε ₀ is the permittivity of free space, and θ_c is the angle between the crystalline c-axis and the y-axis in the y-z plane. To compare with [16

], the fundamental transverse-magnetic (TM) like mode is computed at the operating wavelength λ =0.84 μm. As in [16

], by making use of the symmetry of the problem, only the right half of the domain in Fig. 7 is considered in the calculation. The circles in Fig. 4 are the FDFD calculated effective indices at different θ_C ’s, which are seen to agree very well with the solid line obtained using the FE-BPM method. The effective index is defined as β/[ω(μ ₀ ε ₀)^1/2]. Compared with the BPM, the developed standard eigenvalue algorithm has the advantage of directly obtaining solutions from available eigenvalue matrix equation solvers without needing to examine the dependence of numerical convergence on the beam propagation distance [16

]. For the case of θ_C = 50°, we plot the calculated effective index versus the FDFD grid size Δx = Δy in Fig. 5. It is seen that with Δx = Δy = 0.4μm the calculated effective index has reaches three converged digits after the decimal point. The number of field-component unknowns corresponding to Δx = Δy = 0.4μm is 1904, and the computation time based on an Intel Core 2 Duo 1.86-GHz personal computer is 11.24 seconds. When Δx = Δy is reduced to 0.05 μm, the number of unknowns increases to 114400 and the computation time becomes 908.3 seconds. Please note that in this work we use the simple stair-case approximation for modelling the permittivity discontinuity across the dielectric interface. As studied in [10

], for the grid sizes used here, the numerical accuracy in the calculated effective index could be at most on the order of 10^-4 under such approximation. To improve the numerical accuracy, more elaborate treatment of the field continuity conditions across the interface, such as the index-averaging scheme or the proper boundary-condition matching scheme, as detailed in [10

], could be employed. Since the main goal in this work is to provide a solution method for analyzing waveguides with a more general permittivity characteristic, and for practical application purposes, the above-mentioned accuracy would be good enough the issue of further numerical accuracy improvement will not be investigated. The mode-field profiles of the moduli of the three electric field phasor components, ∣E_x ∣, ∣E_y ∣, and ∣E_Z ∣, at θ_C = 30° are shown in Fig. 6(a), (b), and (c), respectively. The presentation is normalized with respect to the peak value of the major field component, ∣E_y ∣. Note that the ∣E_y ∣ profile lines are illustrated down to 10^-4 and the field strength of ∣E_x ∣ and \∣E_Z ∣ is found to be smaller than ~ 10^-4.

Fig. 3. Computational domain and structure parameters for the LN waveguide with PE-LN core.

Fig. 4. Calculated effective index as a function of the crystal angle for the LN waveguide of Fig. 3.

Fig. 5. Numerical convergence versus the grid size for the calculated effective index at θ_C = 50° in Fig. 4.

Fig. 6. Mode-field profiles of (a) ∣E_x ∣, (b) ∣E_y ∣, and (c) ∣E_z ∣ for the TM-like mode of the LN waveguide with θ_c = 30°.

Fig. 7. Computational domain and structure parameters for the LC-core waveguide.

Next, we study an LC optical waveguide. As shown in Fig. 7, we consider a similar kind of structure as the above LN waveguide but with the substrate being glass with the refractive index n_g = 1.45 and the core region being filled with nematic LCs (5CB) [19

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

]. The elements of the relative permittivity tensor of the nematic LCs are given as

Fig. 8. Schematic definition of rotation angles for the LC molecular or director.

Fig. 9. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for modes 1 and 2 of the LC waveguide with θ_c =0°.

ε_{xx} / ε_{0} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \cos^{2} ϕ_{c}

(56)

ε_{xy} / ε_{0} = ε_{yx} = (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \cos ϕ_{c} \cos ϕ_{c}

(57)

ε_{xz} / ε_{0} = ε_{zx} = (n_{e}^{2} - n_{o}^{2}) \sin θ_{c} \cos θ_{c} \cos ϕ_{c}

(58)

ε_{yy} / ε_{0} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \sin^{2} θ_{c} \cos^{2} ϕ_{c}

(59)

ε_{yz} / ε_{0} = ε_{zy} = (n_{e}^{2} - n_{o}^{2}) \sin θ_{c} \cos θ_{c} \sin ϕ_{c}

(60)

ε_{zz} / ε_{0} = n_{o}^{2} + (n_{e}^{2} - n_{o}^{2}) \cos^{2} θ_{c}

(61)

where n_o = 1.5292 and n_e = 1.7072 are, respectively, the ordinary and extraordinary refractive indices of the nematic LCs, θ_C is the angle between the crystal c-axis and the z-axis, and θ_c is the angle between the projection of the crystal c-axis on the x-y plane and the x-axis, as defined in Fig. 8. The waveguide and computational-domain parameters, w, t, W_x , W_y , and d, shown in Fig. 7, are considered to possess the same values as in the above LN waveguide case. The wavelength is taken to be 1.55 μm.

Based on the fact that the LC molecules can be rotated through an applied voltage on suitably designed electrodes, making the change of the permittivity tensor in the LC-core region, we calculate and investigate the guided modes under different LC-rotated situations. We consider θ_c = 0°, 30°, 60°, and 90°, and at each θ_C the angle ϕ_c is varied from 0° to 90°.

For θ_c = 0°, it is obvious that different ϕ_c ’s correspond to the same situation and we obtain two modes with the effective indices of 1.504400 (mode 1) and 1.504382 (mode 2). Note that no spurious-mode problem was encountered in the analysis, including the situations to be discussed in the following. The mode-field profiles of the moduli of the two transverse electric field phasor components, ∣E_x ∣ and ∣E_y ∣, for the two modes are shown in Fig. 9. The presentation is normalized with respect to the peak value of the larger field component. Such normalization in the presentation of mode-field profiles will be similarly adopted in the following discussions of other cases. It is seen that modes 1 and 2 are basically y- and x-polarized, respectively.

Fig. 10. Calculated effective indices versus ϕ_c for the first four modes of the LC waveguide with θ_c = 30°.

Fig. 11. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for the first four modes of the LC waveguide with θ_c = 30° at ϕ_c = 0°, 30°, 60°, and 90°.

Then, we discuss the situation when θ_c = 30°. Now, different ϕ_c corresponds to different core property. In the θ_c = 0° case, due to geometric symmetry and thus the modal symmetry (anti-symmetry), the computational resources can be reduced by considering only the left or right half of the computational domain in Fig. 7. But now, although the waveguide is with symmetric background structure, the θ_c = 30° orientation has caused the loss of transverse-axis symmetries of the problem and we have to include the whole cross-section as shown in Fig. 7 in the computation. Asymmetric mode-field profiles are obtained as will be presented below.

Fig. 12. Calculated effective indices versus ϕ_c for the first seven modes of the LC waveguide with θ_C = 60°.

Fig. 13. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for the first four modes of the LC waveguide with θ_c = 60° at ϕ_c = 0°, 30°, 60°, and 90°.

The calculated effective indices for different ϕ_c ’s for the obtained four modes are shown in Fig. 10. The ∣E_x ∣ and ∣E_y ∣ profiles of these four modes for ϕ_c = 0°, 30°, 60°, and 90° are illustrated in Fig. 11. Except for ϕ_c = 90°, which corresponds to both structure geometry and permittivity tensor symmetries, the field profiles appear to possess asymmetric characteristics of various degrees. For mode 1, it is almost x-polarized at ϕ_c = 0° and 30°, and changes to be almost y-polarized at ϕ_c = 60° and 90°. For modes 2 and 3, the ∣E_x ∣ component is more significant than the ∣E_y ∣ component at ϕ_c = 0° and 30°, and the ∣Ey∣ component becomes more significant than the ∣E_x ∣ component at ϕ_c = 60° and 90°. For mode 4, such polarization significance situation reverses at the corresponding ϕ_c ’s.

The results for θ_c = 60° are given in Figs. 12–14. Seven modes are obtained as seen in Fig.12. Field-profile plots for modes 1–4 are shown in Fig. 13 and those for modes 5–7 are given in Fig. 14. The off-x-y-plane anisotropic nature of the LC-core makes the correspondence between the mode order and the mode pattern feature a little more complicated. For example, at ϕ_c = 0° in Fig. 14, the ∣E_x ∣ pattern of mode 5 appears to have two nodal positions along the y-direction and that of mode 6 has two nodal positions along the x-direction. However, at ϕ_c = 30°, these two respective patterns correspond to reversed mode orders, purely resulting from the different permittivity tensors. Note that we have ordered the modes according to the magnitude of the effective index. Again, at ϕ_c = 90°, the obtained mode profiles possess the required symmetry or anti-symmetry.

Fig. 14. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for modes 5–7 of the LC waveguide with θ_c = 60° at ϕ_c = 0°, 30°, 60°, and 90°.

Fig. 15. Calculated effective indices versus ϕ_c for the first seven modes of the LC waveguide with θ_c = 90°.

Finally, the results for θ_c = 90° are presented in Figs. 15–17. The LC director is now rotated about the z-axis in the x-y plane as ϕ_c is varied. Seven modes are obtained again. As seen in Fig. 15, the relative magnitude differences as well as the ϕ_c dependence of the calculated effective indices of the first four modes are quite similar to those of the first four modes in Fig. 12, except that the effective index is about 0.04 up shifted in Fig. 15 compared to in Fig. 12. And it is observed that one-to-one similarity between Figs. 13 and 16 occurs in most of the panels. But the behaviors of modes 5–7 in Fig. 15 appear to be of more significant difference from those of modes 5–7 in Fig. 12. Such more significant difference also appears in the panel-by-panel comparison between Figs. 14 and 17. It can be seen that at ϕ_c = 30° mode patterns of modes 5 and 6 in Fig. 17 are closer to those of modes 6 and 5, respectively, in Fig. 14, and at ϕ_c = 90° mode patterns of modes 6 and 7 in Fig. 17 are closer to those of modes 7 and 6, respectively, in Fig. 14.

Fig. 16. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for the first four modes of the LC waveguide with θ_c = 90° at ϕ_c = 0°, 30°, 60°, and 90°.

Fig. 17. ∣E_x ∣ and ∣E_y ∣ mode-field profiles for modes 5–7 of the LC waveguide with θ_c = 90° at ϕ_c = 0°, 30°, 60°, and 90°.

4. Conclusion

We have presented a new FDFD method based eigenvalue algorithm for computing guided modes of anisotropic optical waveguides with arbitrary permittivity tensor. Yee’s mesh is employed in the formulation and a standard eigenvalue matrix equation involving four transverse field components simultaneously has been successfully derived. The PML is utilized as the absorption boundary condition for the computational domain. The new algorithm has been first examined by its application to the analysis of a proton-exchanged LiNbO₃ optical waveguide and its agreement with the analysis by a full-vectorial finite-element beam propagation method was demonstrated. Eigenmodes can be more efficiently solved from the standard eigenvalue matrix equation, compared with the beam propagation method that needs an additional consideration about its numerical convergence with the propagation distance. Then, the algorithm has been used to solve guided modes on a liquid-crystal optical waveguide with arbitrary molecular director orientation. The calculated effective indices and mode-field profiles were presented for different θ_c ’s and ϕ_c ’s of the director orientation. This arbitrary orientation may cause the loss of transverse-axis symmetries of the waveguide with symmetric background structure, as was clearly demonstrated by the asymmetric mode-field patterns in the situations with θ_c ≠ 0° and ϕ_c ≠ ±90°. This established mode solver provides an efficient tool for studying and designing waveguides with complicated dielectric materials such as liquid crystals.

Acknowledgments

This work was supported in part by the National Science Council of the Republic of China under grants NSC96-2221-E-002-097 and NSC97-2221-E-002-043-MY2, in part by the Excellent Research Projects of National Taiwan University under grant 97R0062-07, and in part by the Ministry of Education of the Republic of China under “The Aim of Top University Plan” grant.

References and links

1.	C. L. Xu, W. P. Huang, M. S. Stern, and S. K. Chaudhuri, “Full-vectorial mode calculations by finite difference method,” Proc. Inst. Electr. Eng. 141, 281–286 (1994). [CrossRef]
2.	P. Lüsse, P. Stuwe, J. Schule, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–494 (1994). [CrossRef]
3.	G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis I: Uniform regions and dielectric interfaces,” J. Lightwave Technol. 20, 1210–1218 (2002). [CrossRef]
4.	G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners,” J. Lightwave Technol. 20, 1219–1231 (2002). [CrossRef]
5.	P. J. Chiang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44, 56–66, (2008). [CrossRef]
6.	N. Thomas, P. Sewell, and T. M. benson, “A new full-vectorial higher order finite-difference scheme for the modal analysis of rectangular dielectric waveguides,” J. Lightwave Technol. 25, 2563–2570 (2007). [CrossRef]
7.	K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. AP-14, 302–307 (1966).
8.	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, 1627–1634 (2002). [CrossRef]
9.	Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10, 853–864 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853. [PubMed]
10.	C. P. Yu and H. C. Chang, “Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Opt. Express , 12, 6165V–6177 (2004). http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-25-6165. [CrossRef]
11.	L. Thylen and D. Yevick, “Beam propagation method in anisotropic media,” Appl. Opt. 21, 2751–2754 (1982). [CrossRef] [PubMed]
12.	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, 1926–1931 (1994). [CrossRef]
13.	C. L. D. S. Sobrinho and A. J. Giarola, “Analysis of biaxially anisotropic dielectric waveguides with Gaussian-Gaussian index of refraction profiles by the finite-difference method,” IEE Proc.-H 140, 224–230 (1993).
14.	P. Lüsse, K. Ramm, and H.-G. Unger, “Vectorial eigenmode calculation for anisotropic planar optical waveguides,” Electron. Lett. 32, 38–39 (1996). [CrossRef]
15.	A. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector finite difference modesolver for anisotropic dielectric waveguides,” J. Lightwave Technol. 26, 1423–1431 (2008). [CrossRef]
16.	K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19, 405–413 (2001). [CrossRef]
17.	F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media”, IEEE Microwave Guid. Wave Lett. 8, 223V–225 (1998). [CrossRef]
18.	R. Mittra and U. Pekel, “A new look an the perfectly matched layer PML concept for the reflectionless absorption of electromagnetic waves,” IEEE Microwave Guid. Wave Lett. 5, 84–86 (1995). [CrossRef]
19.	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.7370) Optical devices : Waveguides
(260.2110) Physical optics : Electromagnetic optics

ToC Category:
Physical Optics

History
Original Manuscript: January 28, 2009
Revised Manuscript: March 18, 2009
Manuscript Accepted: March 18, 2009
Published: March 30, 2009

Citation
Ming-yun Chen, Sen-ming Hsu, and Hung-Chun Chang, "A finite-difference frequency-domain method for full-vectroial mode solutions of anisotropic optical waveguides with arbitrary permittivity tensor," Opt. Express 17, 5965-5979 (2009)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-8-5965

Sort: Author | Year | Journal | Reset

References

C. L. Xu, W. P. Huang, M. S. Stern, and S. K. Chaudhuri, "Full-vectorial mode calculations by finite difference method," Proc. Inst. Electr. Eng. 141, 281-286 (1994). [CrossRef]
P. L¨usse, P. Stuwe, J. Schule, and H.-G. Unger, "Analysis of vectorial mode fields in optical waveguides by a new finite difference method," J. Lightwave Technol. 12, 487-494 (1994). [CrossRef]
G. R. Hadley, "High-accuracy finite-difference equations for dielectric waveguide analysis I: Uniform regions and dielectric interfaces," J. Lightwave Technol. 20, 1210-1218 (2002). [CrossRef]
G. R. Hadley, "High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners," J. Lightwave Technol. 20, 1219-1231 (2002). [CrossRef]
P. J. Chiang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, "Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations," IEEE J. Quantum Electron. 44, 56-66, (2008). [CrossRef]
N. Thomas, P. Sewell, and T. M. benson, "A new full-vectorial higher order finite-difference scheme for the modal analysis of rectangular dielectric waveguides," J. Lightwave Technol. 25, 2563-2570 (2007). [CrossRef]
K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media," IEEE Trans. Antennas Propag. AP-14, 302-307 (1966).
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, 1627-1634 (2002). [CrossRef]
Z. Zhu and T. G. Brown, "Full-vectorial finite-difference analysis of microstructured optical fibers," Opt. Express 10, 853-864 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-17-853. [PubMed]
C. P. Yu, and H. C. Chang, "Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers," Opt. Express, 12, 6165V-6177 (2004). http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-25-6165. [CrossRef]
L. Thylen and D. Yevick, "Beam propagation method in anisotropic media," Appl. Opt. 21, 2751-2754 (1982). [CrossRef] [PubMed]
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, 1926-1931 (1994). [CrossRef]
C. L. D. S. Sobrinho and A. J. Giarola, "Analysis of biaxially anisotropic dielectric waveguides with Gaussian-Gaussian index of refraction profiles by the finite-difference method," IEE Proc.-H 140, 224-230 (1993).
P. Lüsse, K. Ramm, and H.-G. Unger, "Vectorial eigenmode calculation for anisotropic planar optical waveguides," Electron. Lett. 32, 38-39 (1996). [CrossRef]
A. B. Fallahkhair, K. S. Li, and T. E. Murphy, "Vector finite difference modesolver for anisotropic dielectric waveguides," J. Lightwave Technol. 26, 1423-1431 (2008). [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, 405-413 (2001). [CrossRef]
F. L. Teixeira and W. C. Chew, "General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media," IEEE Microwave Guid. Wave Lett. 8, 223V-225 (1998). [CrossRef]
R. Mittra and U. Pekel, "A new look an the perfectly matched layer PML concept for the reflectionless absorption of electromagnetic waves," IEEE Microwave Guid. Wave Lett. 5, 84-86 (1995). [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

Figures


Fig. 1.	Fig. 2.	Fig. 3.


Fig. 4.	Fig. 5.	Fig. 6.


Fig. 7.	Fig. 8.	Fig. 9.


Fig. 10.	Fig. 11.	Fig. 12.


Fig. 13.	Fig. 14.	Fig. 15.


Fig. 16.	Fig. 17.