Open Access
Add to BibSonomyAbstract
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 LiNbO3 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.
©2009 Optical Society of America
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] or the Helmholtz equation [2]. Also high-accuracy algorithms have been pursued [3–6]. 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], and derived directly from Maxwell’s equations have been formulated [8–10]. 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] and they at most considered anisotropy in the transverse plane, i.e., in the waveguide cross-section. Most recently, Fallahkhair et al. [15] 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] 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] 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 LiNbO3 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.
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]. For frequency-domain mode solutions, Maxwell’s curl equations are written as
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 = μ 0, where μ 0 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:
Applying the Yee’s mesh, as depicted in Fig. 2, and the central difference scheme, Eqs (3)–(8) become
which can be written in the matrix form as
In Eqs. (15)–(16), I is the identity matrix, Ux, Uy, Vx and Vy are the matrices given as below
and each of Ex, Ey, Ez, Hx, Hy, and Hz 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]. Now, for the general anisotropic dielectric medium, the relation between E⃗ and D⃗ is expressed at each spatial point by
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
Note that in Eqs. (27)–(29), each εijEj (i, j = x, y, z) should be considered as a column vector with each element representing the value of εijEj 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 Dx, Dy, and Dz, respectively. By substituting the Hz and Ez 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:
where the sub-matrices Aij (i, j = 1,⋯ ,4) are defined as below:
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
where sx, sy, and sz are the complex PML parameters defined as
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
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].
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 LiNbO3 (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, Wx = 9μm, Wy = 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 no = 2.250 and ne = 2.172, respectively. The refractive index difference between the PE-LN core region and the substrate is assumed to be Δne = 0.01 for the extraordinary index and Δno = 0 for the ordinary index. The relative permittivity tensor elements for LN are given by
where ε 0 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 β/[ω(μ 0 ε 0)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, ∣Ex∣, ∣Ey∣, and ∣EZ∣, 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, ∣Ey∣. Note that the ∣Ey∣ profile lines are illustrated down to 10-4 and the field strength of ∣Ex∣ and \∣EZ∣ is found to be smaller than ~ 10-4.
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) ∣Ex∣, (b) ∣Ey∣, and (c) ∣Ez∣ for the TM-like mode of the LN waveguide with θc = 30°.
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 ng= 1.45 and the core region being filled with nematic LCs (5CB) [19]. The elements of the relative permittivity tensor of the nematic LCs are given as
where no = 1.5292 and ne = 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, Wx, Wy, 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, ∣Ex∣ and ∣Ey∣, 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. ∣Ex∣ and ∣Ey∣ 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. ∣Ex∣ and ∣Ey∣ 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 ∣Ex∣ and ∣Ey∣ 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 ∣Ex∣ component is more significant than the ∣Ey∣ component at ϕc = 0° and 30°, and the ∣Ey∣ component becomes more significant than the ∣Ex∣ 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 ∣Ex∣ 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. ∣Ex∣ and ∣Ey∣ 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. ∣Ex∣ and ∣Ey∣ mode-field profiles for the first four modes of the LC waveguide with θc = 90° at ϕc = 0°, 30°, 60°, and 90°.
Fig. 17. ∣Ex∣ and ∣Ey∣ 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 LiNbO3 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).
