Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Full-vectorial meshless finite cloud method for an anisotropic optical waveguide analysis

Xiaoer Wu and Jinbiao Xiao

Open Access Open Access

Abstract

By using the meshless finite cloud method, an efficient full-vectorial mode solver based on the transverse-magnetic-field components is developed to analyze the optical waveguides made of anisotropic materials, in which the waveguide cross-section enclosed by the perfectly matched layers is divided into an appropriate number of homogeneous clouds. The point collocation technique is utilized to create a scattered set of nodes over the cloud, and then the continuity conditions of the longitudinal field components are imposed to appropriately deal with the discrete nodes at the interfaces shared by the adjacent clouds. In comparison with conventional mesh-based numerical techniques, the distributions of solution nodes of the present method can be applied to the area of complexity in a completely free manner. Moreover, an interior nodal distribution adaptively updating along the propagation direction is adopted to accomplish higher computational efficiency while improving numerical accuracy. To validate the proposed method, an anisotropic square waveguide, a magneto-optical raised strip waveguide, and a nematic liquid-crystal channel waveguide are considered as numerical examples, and their modal field distribution and corresponding effective refractive indexes are presented. Results are in good agreement with those published earlier, showing the effectiveness of the present method.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Anisotropic optical waveguides have been widely used in photonic integrated circuits (PICs) or planar lightwave circuits (PLCs) to construct various kinds of optical devices like polarization converters/polarizers [1–3], optical modulators/switches [4,5], optical isolators or polarization splitters [6–8]. As a general rule, for an anisotropic medium, two factors are mandatorily considered, namely, the different polarization behavior and the coupling among the field components. Therefore, developing an efficient full-vectorial mode solver of anisotropic optical waveguides is necessary for clearly analyzing and accurately designing the desired devices. However, very few analytic solutions for computing the guided modes in anisotropic waveguides can be obtained. Therefore, the use of numerical technique is important and perhaps indispensable.

Mesh-based and meshless methods are among the tools normally applied in the numerical computation to investigate modal characteristics of a wide range of anisotropic optical waveguides. The mesh-based methods, which require the meshing of the considered computational domain and the information of the connectivity among the elements, such as the finite difference (FD) method [9–11], the finite element (FE) method [12–15], and the finite-difference frequency-domain (FDFD) method [16,17], have been proposed and applied to analyze the guided modes of the anisotropic optical waveguides. Among them, the FD method is one of the most popular and powerful mesh-based techniques, mostly owing to the explicit, closed-form expression for its differential operator, and visualized, manipulated treatment of solutions obtained by the FD schemes. However, the FD method is usually based on low-order basis functions, normally resulting in large matrix with poor convergence and accuracy. Moreover, the computational time dramatically increases with the increment of the meshed nodes, especially for complex waveguide structures or those waveguides with high refractive index contrast. On the basis of the above fact, meshless methods have emerged as alternatives to the most often employed mesh-based methods, in which the domain of interest is only represented by preselected basis sets. Among them, the most frequently used are the Galerkin [18,19] and collocation methods [20–26]. As a general rule, the collocation method is more superior to the Galerkin method, mainly due to the latter's requirement of considerable evaluations of laborious integral elements. Broadly defined, the commonly used collocation schemes can be divided into two categories. The first category is schemes based on spectral algorithm, namely, the so-called spectral collocation method, in which a preselected orthogonal basis function is used to approximate the unknown functions. The second category is schemes based on free nodal distribution, namely, the point collocation method, where the objective equations are satisfied at each of the nodes covering the computational domain. Currently, several spectral techniques have been proposed, including the multidomain pseudospectral (MPS) method [27], the multidomain collocation (MC) method [28] and the pseudospectral (PS)-based eigenvalue method [29]. In spectral schemes, the redundant integral procedure is avoided because the unknown fields are applied at a set of interpolation nodes to satisfy the wave equation. To further avoid the Gibbs phenomenon, the computational domain in spectral scheme is partitioned into a number of subdomains according to its geometry, and then the mode field profiles in different subdomains are expanded by imposing distinct basis functions. Previous reports showed that the spectral scheme using single basis functions requires more computational time and memory versus the one using different basis functions. In comparison with the FD and FE methods, the spectral approach achieves considerably less memory and better efficiency. However, the scaling factors of the basis functions used in distinct subdomains, which may increase the computational complexity, should be carefully chosen. Moreover, the spectral techniques are often based on mesh or background mesh, which strongly affect the computational efficiency. Consequently, for eigenmode analysis, the point collocation method is more appropriate and efficient.

Recently, a novel technique that was originally proposed to solve the two-dimensional (2-D) elastostatics problem, namely, the finite cloud (FC) method, combining with fast convergence and high-order accuracy, has been introduced to the area of computational electromagnetics [30–33]. The FC method, which uses a point collocation technique to create a scattered set of nodes over the computational domain without requiring geometric knowledge or connectivity information for linking any of the nodes, is a true meshless technique for solving the partial differential equations (PDEs). Thus, instead of regular grids as in FD scheme and background cells as in spectral scheme, the FC method uses an unconstrained nodal distribution. This means that the nodes may not be given in a regular distribution and may vary from a high-density area with a high solution gradient to a low-density area with a low solution gradient. Moreover, the method is also quite suitable for adaptive meshing, mostly owing to the easy addition of extra nodes in any desired areas. As a result, its computational efficiency is superior to the traditional FD and spectral collocation methods. On the other hand, the simplicity of the mathematical bases is also a major consideration for adopting this scheme. In the construction of meshless interpolation functions, the most often used technique is the fixed kernel approximations. Because their major parts, including quadratic basis function, correction function and kernel function, are known priors, major parts of problems can be tackled analytically. However, previous works only considered the isotropic medium [30,32]. Therefore, the aim of the study is to develop a full-vectorial mode solver for anisotropic optical waveguides by using the FC method.

In this paper, a full-vectorial mode solver in terms of transverse magnetic field components for optical waveguides with transverse anisotropy is developed by using the FC method. The robust perfectly matched layers (PMLs) for anisotropic media are incorporated into the full-vectorial wave equations for effectively avoiding the nonphysical reflection from the computational window edges. Moreover, a combination of the point collocation with the modified Gaussian function is used to construct meshless shape functions. Therefore, the computational efficiency can be improved without degrading numerical accuracy. In order to validate the present method, three typical optical dielectric waveguides with anisotropy are considered as numerical examples, and the results are compared with those from the FD, MC and MPS methods. The rest of this paper is organized as follows. The mathematical formulations for the full-vectorial anisotropic wave equations by using the FC method are described in Section II in detail. Section III presents numerical results and discussion, analyzing the modal characteristics for three kinds of anisotropic waveguides. The conclusion is drawn in Section IV.

2. Mathematical formulation

To incorporate the PML boundary conditions into the derivation of the wave equation, the following change of variables is introduced into the Maxwell equations [19,24]:

$$\xi = \int_0^\xi {{s_\xi }({\xi ^{\prime}})} d{\xi ^{\prime}}$$
where ${s_\xi }(\xi )$ represents the complex stretching variables, and $\xi $ stands for $p,\textrm{ }\theta $ or $z$ in the CCS. The typical definition of ${s_\xi }(\xi )$ can be expressed by [19,24]
$${s_\xi } = \left\{ {\begin{array}{lr} {1 - j{\tau_\xi }(\xi )/\omega {\varepsilon_0},}&{\textrm{ in the PML region}}\\ {1,}&{\textrm{ in the non - PML region}} \end{array}} \right.$$
with
$${\tau _\xi } = \frac{{m + 1}}{{150\pi \Delta \xi }}\frac{{|{\xi - {\xi_0}} |}}{{{d^2}}},$$
where d is the thickness of the PML, $\Delta \xi $ is the mesh size along $\xi $ direction, ${\xi _0}$ is the PML interface, and the integral parameter m is taken as 2.

Assuming monochromatic electromagnetic fields that is propagating along the z-direction in an inhomogeneous medium with permittivity tensor $\tilde{\varepsilon }$ and has angular frequency $\omega $, thus, the full-vectorial wave equation in terms of the magnetic field can be written as

$$\tilde{\nabla } \times ({\tilde{\varepsilon }^{ - 1}}\tilde{\nabla } \times \overrightarrow H ) - {\omega ^2}{\mu _0}\overrightarrow H = 0,$$
where ${\mu _0}$ is the magnetic permeability in vacuum and the operator $\tilde{\nabla }$ is constructed as
$$\tilde{\nabla }\textrm{ = }{\alpha _x}\frac{\partial }{{\partial x}}\tilde{x} + {\alpha _y}\frac{\partial }{{\partial y}}\tilde{y} + {\alpha _z}\frac{\partial }{{\partial z}}\tilde{z},$$
where ${\alpha _x}$, ${\alpha _y}$ and ${\alpha _z}$ represents the complex PML parameters. In this study, the transverse anisotropic dielectric medium is considered so that the tensor $\tilde{\varepsilon }$ is taken as follows
$$\tilde{\varepsilon } = \left[ {\begin{array}{ccc} {{\varepsilon_{xx}}}&{{\varepsilon_{xy}}}&0\\ {{\varepsilon_{yx}}}&{{\varepsilon_{yy}}}&0\\ 0&0&{{\varepsilon_{zz}}} \end{array}} \right].$$

In the calculation of the optical waveguides that are assumed to have uniform index profiles along the z-direction, solutions of guided modes can be written as $exp ( - j\beta z)$, where $\beta = {k_0}{n_{\textrm{eff}}},$ ${k_0}$ is the wave number in vacuum and ${n_{\textrm{eff}}}$ is the effective refractive index. Using the divergence relation $\nabla \cdot \overrightarrow H = 0$, the coupled full-vectorial eigenvalue equations based on the transverse field components ${H_x}$ and ${H_y}$ can be derived as

$$\left[ {\begin{array}{cc} {{A_{xx}}}&{{A_{xy}}}\\ {{A_{yx}}}&{{A_{yy}}} \end{array}} \right]\left[ {\begin{array}{c} {{H_x}}\\ {{H_y}} \end{array}} \right] = {\beta ^2}\left[ {\begin{array}{c} {{H_x}}\\ {{H_y}} \end{array}} \right].$$

The ${A_{xx}}$, ${A_{xy}}$, ${A_{yx}}$ and ${A_{yy}}$ are the differential operators, given as follows

$${A_{xx}} = \frac{{{\partial ^2}(I \cdot )}}{{\partial {{\tilde{x}}^2}}} + {\varepsilon _{yy}}\frac{\partial }{{\partial \tilde{y}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}}) + {\varepsilon _{yx}}\frac{\partial }{{\partial \tilde{y}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}}) + {\omega ^2}{\mu _0}{\varepsilon _{yy}},$$
$$\begin{aligned} {A_{xy}} = \frac{{{\partial ^2}(I \cdot )}}{{\partial \tilde{y}\partial \tilde{x}}} - {\varepsilon _{yy}}\frac{\partial }{{\partial \tilde{y}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}}) - {\varepsilon _{yx}}\frac{\partial }{{\partial \tilde{x}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}}) - {\omega ^2}{\mu _0}{\varepsilon _{yx}}, \\ {A_{yx}} = \frac{{{\partial ^2}(I \cdot )}}{{\partial \tilde{x}\partial \tilde{y}}} - {\varepsilon _{xx}}\frac{\partial }{{\partial \tilde{x}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}}) - {\varepsilon _{xy}}\frac{\partial }{{\partial \tilde{y}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}}) - {\omega ^2}{\mu _0}{\varepsilon _{xy}}, \\ {A_{yy}} = \frac{{{\partial ^2}(I \cdot )}}{{\partial {{\tilde{y}}^2}}} + {\varepsilon _{xx}}\frac{\partial }{{\partial \tilde{x}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}}) + {\varepsilon _{xy}}\frac{\partial }{{\partial \tilde{x}}}(\frac{1}{{{\varepsilon _{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}}) + {\omega ^2}{\mu _0}{\varepsilon _{xx}},\end{aligned}$$
where I is an identity matrix and ${\partial / {\partial \tilde{\gamma } \to {\alpha _\gamma }}}{\partial / {\partial \gamma (\gamma = x,y)}}.$ In order to utilize the FC method to solve the derived full-vectorial wave equations for anisotropic optical waveguides, the cross-section of the considered waveguide is divided into several clouds according to its geometry, and then assembled from each cloud while satisfying the physical boundary conditions. This will be followed by a description of the derivation of the full-vector mode solver based on the FC method.

In a fixed kernel approach, for an arbitrary node $({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ of a single cloud $r,$ the approximate solution $H_{xi}^r({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ to the unknown field value $H_{xi}^{ur}({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ is given by

$$H_{xi}^r({\tilde{x}_{ri}},{\tilde{y}_{ri}})\textrm{ = }\int_\Omega {\zeta ({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}},{{\tilde{x}}_{rk}} - s,{{\tilde{y}}_{rk}} - t)} \phi ({\tilde{x}_{rk}},{\tilde{y}_{rk}})H_{xi}^{ur}({\tilde{x}_{ri}},{\tilde{y}_{ri}})dsdt,$$
where $i = 1,2, \cdots ,NP$ and $NP$ is the whole number of nodes in the cloud. $({\tilde{x}_{rk}},{\tilde{y}_{rk}})$ is the kernel node. $\zeta ({\tilde{x}_{ri}},{\tilde{y}_{ri}},{\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t)$ is the correction function and is given by
$$\zeta ({\tilde{x}_{ri}},{\tilde{y}_{ri}},{\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t) = {P^T}({x_{rk}} - s,{\tilde{y}_{rk}} - t)c({\tilde{x}_{ri}},{\tilde{y}_{ri}}),$$
where ${P^T}({x_{rk}} - s,{\tilde{y}_{rk}} - t)$ is a quadratic basis function and is given by
$${P^T}({\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t) = [1,{\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t,{({\tilde{x}_{rk}} - s)^2},({\tilde{x}_{rk}} - s)({\tilde{y}_{rk}} - t),{({\tilde{y}_{rk}} - t)^2}].$$

And $c({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ is a $6 \times 1$ vector of correction function coefficients, which can be determined by satisfying the reproducing or the consistency conditions, i.e.

$$\int_\Omega {{P^T}({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}})} c({\tilde{x}_{ri}},{\tilde{y}_{ri}})\phi ({\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t){p_n}(s,t)dsdt = {p_n}({\tilde{x}_{ri}},{\tilde{y}_{ri}}),n = 1,2, \cdots ,6.$$

In discrete form, the above equation can be written as

$$\sum\limits_{I = 1}^{NP} {{P^T}({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}})c({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}})\phi ({{\tilde{x}}_{rk}} - {{\tilde{x}}_{rI}},{{\tilde{y}}_{rk}} - {{\tilde{y}}_{rI}}){p_n}({{\tilde{x}}_{rI}},{{\tilde{y}}_{rI}})} \Delta {V_{rI}} = {p_n}({\tilde{x}_{ri}},{\tilde{y}_{ri}}),$$
where $\Delta {V_I}$ is the nodal volume associated with node I, and Eq. (13) can be written as
$$\overrightarrow M c({\tilde{x}_{ri}},{\tilde{y}_{ri}})\textrm{ = }p({\tilde{x}_{ri}},{\tilde{y}_{ri}}),$$
where $\overrightarrow M $ is a $6 \times 6$ moment matrix and is a constant. The entries in the moment matrix are given by
$${M_{jl}} = \sum\limits_{I = 1}^{NP} {\{{{p_j}({{\tilde{x}}_{rI}},{{\tilde{y}}_{rI}})\zeta ({{\tilde{x}}_{rk}} - {{\tilde{x}}_{rI}},{{\tilde{y}}_{rk}} - {{\tilde{y}}_{rI}}){p_l}({{\tilde{x}}_{rI}},{{\tilde{y}}_{rI}})\Delta {V_I}} \}} ,j,l = 1,2, \cdots ,6.$$

From Eq. (14), the correction function coefficients are computed as

$$c({\tilde{x}_{ri}},{\tilde{y}_{ri}})\textrm{ = }{\overrightarrow M ^{ - 1}}p({\tilde{x}_{ri}},{\tilde{y}_{ri}}).$$

Substituting the above correction function coefficients into Eq. (10), the correction function is rewritten as

$$\zeta ({\tilde{x}_{ri}},{\tilde{y}_{ri}},{\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t)\textrm{ = }{P^T}({\tilde{x}_{ri}},{\tilde{y}_{ri}}){\overrightarrow M ^{ - 1}}p({\tilde{x}_{ri}},{\tilde{y}_{ri}}).$$

The kernel function $\phi ({\tilde{x}_{rk}},{\tilde{y}_{rk}})$ is centered at $({\tilde{x}_{rk}},{\tilde{y}_{rk}}).$ In the present work, $\phi $ is taken as a modified Gaussian function, i.e.,

$$\phi ({\tilde{x}_{rk}} - {\tilde{x}_{rI}}) = \frac{{\omega ({{\tilde{x}}_{rk}} - {{\tilde{x}}_{rI}})}}{{1 - \omega ({{\tilde{x}}_{rk}} - {{\tilde{x}}_{rI}}) + \varepsilon ^{\prime}}},$$
where $\varepsilon ^{\prime}$ is a small number for avoiding the singularity of $\phi ({\tilde{x}_{rk}} - {\tilde{x}_{rI}}).$ $\omega ({\tilde{x}_{rk}} - {\tilde{x}_{rI}})$ is a normalized Gaussian function given by
$$\omega\left(\tilde{x}_{r k}-\tilde{x}_{r I}\right)= \begin{cases}0, & t>d_{r} \\ \frac{e^{-\left(d_{r}\left(\tilde{x}_{i k}-\tilde{x}_{I}\right) / 2\right)^{2}}-e^{-4}}{1-e^{-4}}, t \leq d_{r}\end{cases},$$
where ${d_r}$ is the support size of the cloud $r.$ In 2-D, the kernel function is constructed as
$$\phi ({\tilde{x}_{rk}} - {\tilde{x}_{rI}},{\tilde{y}_{rk}} - {\tilde{y}_{rI}}) = \phi ({\tilde{x}_{rk}} - {\tilde{x}_{rI}})\phi ({\tilde{y}_{rk}} - {\tilde{y}_{rI}}).$$

Substituting Eqs. (17) and (20) into Eq. (9), the approximate field value $H_{xi}^r({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ of the node $({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ is obtained as

$$H_{xi}^r({\tilde{x}_{ri}},{\tilde{y}_{ri}}) = \sum\limits_{I = 1}^{NP} {{N_{rI}}({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}})} H_{xI}^r,$$
similarly,
$$H_{yi}^r({\tilde{x}_{ri}},{\tilde{y}_{ri}}) = \sum\limits_{I = 1}^{NP} {{N_{rI}}({{\tilde{x}}_{ri}},{{\tilde{y}}_{ri}})} H_{yI}^r,$$
where $H_{xI}^r$ and $H_{yI}^r$ are the field values of node $I.$ ${N_{rI}}({\tilde{x}_{ri}},{\tilde{y}_{ri}})$ is the fixed kernel shape function and represents the single nodal parameter, defined as
$${N_{rI}}({\tilde{x}_{ri}},{\tilde{y}_{ri}}) = {P^T}({\tilde{x}_{ri}},{\tilde{y}_{ri}}){\overrightarrow M ^{ - 1}}p({\tilde{x}_{ri}},{\tilde{y}_{ri}})\zeta ({\tilde{x}_{ri}},{\tilde{y}_{ri}},{\tilde{x}_{rk}} - s,{\tilde{y}_{rk}} - t)\Delta {V_{rI}}.$$

Substituting Eqs. (21) and (22) into Eq. (7), Eq. (7) is demanded to perfectly be satisfied at each node. Equation (7) is thus converted to a system of linear equations to form a matrix eigenvalue equation with an eigenvalue of ${\beta ^2}$ given by

$$\left[ {\begin{array}{cc} {N_{xx}^r}&{N_{xy}^r}\\ {N_{yx}^r}&{N_{yy}^r} \end{array}} \right]\left[ {\begin{array}{c} {H_x^r}\\ {H_y^r} \end{array}} \right] = {\beta ^2}\left[ {\begin{array}{c} {H_x^r}\\ {H_y^r} \end{array}} \right],$$
where the operators $N_{xx}^r$, $N_{xy}^r$, $N_{yx}^r$ and $N_{yy}^r$ are defined as
$$\begin{array}{l} N_{xx}^rH_x^r = {\left. {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial {{\tilde{x}}^2}}} + \frac{{\varepsilon_{yy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} + {\omega^2}{\mu_0}\varepsilon_{yy}^r} \right]H_x^r} \right|_{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}\\ = {\sum\limits_{i = 1}^{NP} {\sum\limits_{j = 1}^{NP} {\left. {\left\{ {\sum\limits_{l = 1}^{NP} {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial {{\tilde{x}}^2}}} + \frac{{\varepsilon_{yy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} + {\omega^2}{\mu_0}\varepsilon_{yy}^r} \right]{N_{rI}}[{H_{xl}^r} ]} } \right\}} \right|} } _{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}, \end{array}$$
$$N_{xy}^rH_y^r = {\left. {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial \tilde{y}\partial \tilde{x}}} - \frac{{\varepsilon_{yy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} - \frac{{\varepsilon_{yx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} - {\omega^2}{\mu_0}\varepsilon_{yx}^r} \right]H_y^r} \right|_{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}$$
$$= {\sum\limits_{i = 1}^{NP} {\sum\limits_{j = 1}^{NP} {\left. {\left\{ {\sum\limits_{l = 1}^{NP} {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial \tilde{y}\partial \tilde{x}}} - \frac{{\varepsilon_{yy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} - \frac{{\varepsilon_{yx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} - {\omega^2}{\mu_0}\varepsilon_{yx}^r} \right]{N_{rI}}[{H_{yl}^r} ]} } \right\}} \right|} } _{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}},$$
$$N_{yx}^rH_x^r = {\left. {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial \tilde{x}\partial \tilde{y}}} - \frac{{\varepsilon_{xx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} - \frac{{\varepsilon_{xy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} - {\omega^2}{\mu_0}\varepsilon_{xy}^r} \right]H_x^r} \right|_{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}$$
$$= {\sum\limits_{i = 1}^{NP} {\sum\limits_{j = 1}^{NP} {\left. {\left\{ {\sum\limits_{l = 1}^{NP} {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial \tilde{x}\partial \tilde{y}}} - \frac{{\varepsilon_{xx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} - \frac{{\varepsilon_{xy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{y}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} - {\omega^2}{\mu_0}\varepsilon_{xy}^r} \right]{N_{rI}}[{H_{xl}^r} ]} } \right\}} \right|} } _{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}},$$
$$N_{xy}^rH_y^r = {\left. {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial {{\tilde{y}}^2}}} + \frac{{\varepsilon_{xx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{1}{{{\varepsilon_{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} + \frac{{\varepsilon_{xy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{1}{{{\varepsilon_{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} + {\omega^2}{\mu_0}\varepsilon_{xx}^r} \right]H_y^r} \right|_{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}$$
$$= {\sum\limits_{i = 1}^{NP} {\sum\limits_{j = 1}^{NP} {\left. {\left\{ {\sum\limits_{l = 1}^{NP} {\left[ {\frac{{{\partial^2}(I \cdot )}}{{\partial {{\tilde{y}}^2}}} + \frac{{\varepsilon_{xx}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{1}{{{\varepsilon_{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{x}}} + \frac{{\varepsilon_{xy}^r}}{{\varepsilon_{zz}^r}}\frac{\partial }{{\partial \tilde{x}}}\frac{1}{{{\varepsilon_{zz}}}}\frac{{\partial (I \cdot )}}{{\partial \tilde{y}}} + {\omega^2}{\mu_0}\varepsilon_{xx}^r} \right]{N_{rI}}[{H_{yl}^r} ]} } \right\}} \right|} } _{\tilde{x} = {{\tilde{x}}_{ri}},\tilde{y} = {{\tilde{y}}_{ri}}}}.$$

Assuming there are m clouds, after assembling the matrix eigenvalue equations in all clouds, the global matrix equation can be obtained, i.e.,

$$\left[ {\begin{array}{ccc} {[{{N_1}} ]}& \cdots &0\\ \vdots & \ddots & \vdots \\ 0& \cdots &{[{{N_m}} ]} \end{array}} \right]\left[ {\begin{array}{c} {[{{H_1}} ]}\\ \vdots \\ {[{{H_m}} ]} \end{array}} \right] = {\beta ^2}\left[ {\begin{array}{c} {[{{H_1}} ]}\\ \vdots \\ {[{{H_m}} ]} \end{array}} \right],$$
where
$$[{{N_r}} ]= \left[ {\begin{array}{cc} {N_{xx}^r}&{N_{xy}^r}\\ {N_{yx}^r}&{N_{yy}^r} \end{array}} \right],[{{H_r}} ]= \left[ {\begin{array}{c} {H_x^r}\\ {H_y^r} \end{array}} \right],r = 1,2, \cdots ,m.$$

It is noted that the nodes shared by the adjacent clouds should be replaced by satisfying the physical interface conditions. Owing to the divergence relation $\nabla \cdot \overrightarrow H = 0,$ the longitudinal field component ${H_Z}$ can be expressed as

$${H_Z} ={-} \frac{j}{\beta }(\frac{{\partial {H_x}}}{{\partial \tilde{x}}} + \frac{{\partial {H_y}}}{{\partial \tilde{y}}}),$$
and ${E_Z}$ can be expressed by
$${E_Z} ={-} \frac{j}{{\omega {\varepsilon _{zz}}}}(\frac{{\partial {H_x}}}{{\partial \tilde{y}}} - \frac{{\partial {H_y}}}{{\partial \tilde{x}}}).$$

Considering a horizontal interface as depicted in Fig. 1(a), the continuity of ${H_Z}$ gives

$${\left. {\frac{{\partial {H_y}}}{{\partial \tilde{y}}}} \right|_{{y_ + }}} = {\left. {\frac{{\partial {H_y}}}{{\partial \tilde{y}}}} \right|_{{y_ - }}}$$
and as well, continuity of ${E_Z}$ yields
$$\varepsilon _{zz}^{{y_ + }}{\left. {\frac{{\partial {H_x}}}{{\partial \tilde{y}}}} \right|_{{y_ - }}} - \varepsilon _{zz}^{{y_ - }}{\left. {\frac{{\partial {H_x}}}{{\partial \tilde{y}}}} \right|_{y + }} = (\varepsilon _{zz}^{{y_ + }} - \varepsilon _{zz}^{{y_ - }})\frac{{\partial {H_y}}}{{\partial \tilde{x}}},$$
where ${y_ + }$ and ${y_ - }$ are respectively referred to the locations at the infinitesimally upper and lower sides of the horizontal interface. Similarly, for a vertical interface as depicted in Fig. 1(b), we have
$${\left. {\frac{{\partial {H_x}}}{{\partial \tilde{x}}}} \right|_{{x_ + }}} = {\left. {\frac{{\partial {H_x}}}{{\partial \tilde{x}}}} \right|_{{x_ - }}}$$
and
$$\varepsilon _{zz}^{{x_ + }}{\left. {\frac{{\partial {H_y}}}{{\partial \tilde{x}}}} \right|_{{x_ - }}} - \varepsilon _{zz}^{{x_ - }}{\left. {\frac{{\partial {H_y}}}{{\partial \tilde{x}}}} \right|_{{x_ + }}} = (\varepsilon _{zz}^{{x_ + }} - \varepsilon _{zz}^{{x_ - }})\frac{{\partial {H_x}}}{{\partial \tilde{y}}},$$
where ${x_ + }$ and ${x_ - }$ are respectively referred to the locations at the infinitesimally right and left of the vertical interface.

 figure: Fig. 1.

Fig. 1. Interfaces of dielectric materials: (a) Horizontal and (b) Vertical interfaces.

Download Full Size | PDF

The boundary conditions can be used to replace the corresponding eigenrow vectors corresponding to the boundary nodes in the sub-matrix. After each node shared by the adjacent clouds is replaced by satisfying the above continuity conditions, a generalized matrix eigenvalue equation is obtained. With the MATLAB subroutines, the obtained matrix is then solved, in which the eigenvectors and eigenvalues are, respectively, related to the modal field distributions and the propagation constants.

3. Results and discussion

A typical anisotropic square waveguide with high index contrast between the core and cladding shown in Fig. 2 is firstly taken as a numerical example to demonstrate the convergence and efficiency of the FC method. The cross-section of the waveguide enclosed by the PMLs is partitioned into several uniform clouds. The computational parameters are: $w = h = 1.0\mathrm{\mu }\textrm{m,}$ $\lambda = 1.3\mathrm{\mu }\textrm{m,}$ and $d = 1.0\mathrm{\mu }\textrm{m}$. The computational window is taken as $4.5\mathrm{\mu }\textrm{m} \times 4.5\mathrm{\mu }\textrm{m}\textrm{.}$ The cladding and core are set to be isotropic and anisotropic, respectively, in which the refractive index of the cladding is taken as ${n_\textrm{c}} = 1.0.$

 figure: Fig. 2.

Fig. 2. Cross section of an anisotropic square waveguide enclosed by the PMLs.

Download Full Size | PDF

In this paper, the angle that the core is homogeneously polarized is fixed at $\theta = {\pi / 4},$ thus, the relative permittivity tensor is given as

$$\tilde{\varepsilon } = \left[ {\begin{array}{ccc} {n_\textrm{o}^2 + (n_\textrm{e}^2 - n_\textrm{o}^2){{\cos }^2}\theta }&{(n_\textrm{e}^2 - n_\textrm{o}^2)\cos \theta \sin \theta }&0\\ {(n_\textrm{e}^2 - n_\textrm{o}^2)\cos \theta \sin \theta }&{n_\textrm{o}^2 + (n_\textrm{e}^2 - n_\textrm{o}^2){{\sin }^2}\theta }&0\\ 0&0&{n_\textrm{o}^2} \end{array}} \right]$$
where ${n_\textrm{e}} = {n_\textrm{f}} + {2 / 3}\Delta n,$ ${n_\textrm{o}} = {n_\textrm{f}} - {1 / 3}\Delta n,$ the refractive index ${n_\textrm{f}}$ is taken as 1.65, and the birefringence of the core is represented as $\Delta n.$ In addition, the small number $\varepsilon ^{\prime}$ is taken as ${10^{ - 6}}$ to effectively avoid the singularity of the kernel function. Considerable computational cost is saved since the mode profiles in the interior clouds and their neighboring clouds are expanded with a modified Gaussian function. In this example and the following two examples, the non-uniform nodal distribution is utilized. For example, if the nodal number used in the x-direction is chosen as 50, the interval among the nodes in the core layer is set to be $0.05\mu m$ and then the nodal number is fixed at 20. Therefore, the corresponding nodal interval in the cladding is chosen as $0.1\mu m$ and the nodal number is fixed at 30.

To effectively evaluate the convergent characteristics of the FC method, the effective indexes of the first- and second-order modes for the anisotropic square waveguide with $\Delta n = 0.015$ versus the transverse nodal number $({N_{{\textrm{T}_\textrm{X}}}})$ are plotted in Fig. 3(a). In particular, it should be mentioned that the vertical nodal number $({N_{{\textrm{T}_\textrm{Y}}}})$ of the computational window is only taken as 37. As is clearly seen, results for the two lowest-order modes approach to the convergent solutions as the transverse nodal numbers of the computational window increase. It also can be observed that the convergent solutions can be obtained when ${N_{{\textrm{T}_\textrm{X}}}}$ is taken as 74, and the results for the two lowest-order modes almost keep constant with further increment of the transverse nodal numbers. Figure 3(b) presents the effective indexes of the first- and second-order modes for the anisotropic square waveguide versus the core birefringence $\Delta n.$ Note that the symbols $n_{\textrm{eff}}^1$ and $n_{\textrm{eff}}^2$ respectively, denote the effective indexes of the first- and second-order modes. The transverse nodal number is set to be ${N_{{\textrm{T}_\textrm{X}}}} = 74.$ In the core layer, the transverse node interval is taken as $0.05\mu m$ and the nodal number is 20, correspondingly, the transverse node interval used in the rest clouds is set to be $0.065\mu m$ and the nodal number is 54. In order to make the comparison more intuitive, results from the FD method described in [9], where a full-vectorial FD discretization is used to deal with more general or non-diagonal anisotropic waveguides, and the MC method described in [28], where the guided mode profiles in the computational subdomains are expanded using the Chebyshev polynomials, are also plotted in the same figure. As can be observed, the effective index of the first- and second-order modes have a noticeable difference since the former linearly increases with the increment of the core birefringence while the latter shows the opposite behavior. It also can be seen that the results from the FC method accord well with those from the MC method. However, the FC method converge apparently faster than the MC method since the former with only 2738 computational nodes can obtain convergent results than the latter with more than 7938 ones. In addition, to obtain the effective index of the field patterns of the two lowest order modes, the computational time of the FC method and the MC method is 9.8 and 12.4 seconds, respectively. It indicates that the FC method is more efficient and practical.

 figure: Fig. 3.

Fig. 3. Effective indexes of the first- and second-order modes for an anisotropic square waveguide versus the transverse nodal number and the core birefringence: (a)${N_{{\textrm{T}_\textrm{X}}}}$ and (b)$\Delta n.$

Download Full Size | PDF

Figure 4 shows the distributions of field intensity of the two lowest-order modes for the anisotropic square waveguide, in which the birefringence of the core and the transverse nodal number are taken as $\Delta n = 0.015$ and ${N_{{\textrm{T}_\textrm{X}}}} = 74,$ respectively. Note that the values of each transverse magnetic field component are normalized to its maximum value. It can be found that, for the first-order mode or the second-order mode, both the transverse magnetic field components ${H_x}$ and ${H_y}$ have the amplitudes of the same order. Therefore, for such waveguide structure, a full-vectorial analysis is indispensable because both the two transverse components are not negligible.

 figure: Fig. 4.

Fig. 4. Distributions of field intensity of the two lowest-order modes for an anisotropic square waveguide: (a)${H_x}$ and (b)${H_y}$ for the first-order mode; (c)${H_x}$ and (d) ${H_y}$ for the second-order mode.

Download Full Size | PDF

Next, the FC method is also applied to compute the eigenmodes of a magneto-optical raised strip waveguide. The cross-section of the waveguide surrounded by the PMLs is divided into nine uniform clouds, as shown in Fig. 5. The waveguide is made of three regions, namely, a core region that is filled with an yttrium iron garnet (YIG), a substrate region that is made of anisotropic gadolinium gallium garnet (GGG) and a cover region that is made of air. The computational parameters are: ${n_\textrm{s}} = 1.95,$ ${n_\textrm{a}} = 1.0,$ $w = 0.8\mathrm{\mu }\textrm{m,}$ $t = 0.7\mathrm{\mu }\textrm{m,}$ $d = 1.0\mathrm{\mu }\textrm{m,}$ and $\lambda = 1.3\mathrm{\mu }\textrm{m}\textrm{.}$ In addition, the computational window is set to be $3.2\mathrm{\mu }\textrm{m} \times 3\mathrm{\mu }\textrm{m}\textrm{.}$ When a static magnetic field is imposed along the propagation direction, the relative permittivity tensor of YIG is defined as

$$\tilde{\varepsilon } = \left[ {\begin{array}{ccc} {n_\textrm{f}^2}&{j\zeta }&0\\ { - j\zeta }&{n_\textrm{f}^2}&0\\ 0&0&{n_\textrm{f}^2} \end{array}} \right]$$
where $\zeta $ stands for the first-order magneto-optic effect and ${n_\textrm{f}} = 2.302.$

 figure: Fig. 5.

Fig. 5. Cross section of a magneto-optical raised strip waveguide surrounded by the PMLs

Download Full Size | PDF

The modal characteristics of this magneto-optical raised strip waveguide are studied in detail. For brevity, the results for the condition $\zeta = 0.005$ are calculated. Figure 6(a) shows the effective indexes of the first- and second-order modes for a magneto-optical raised strip waveguide with respect to the transverse nodal number. Note that the vertical nodal number of the computational window in this example is fixed at the number of ${N_{{\textrm{T}_\textrm{Y}}}} = 24.$ It can be seen that the order of the difference between the effective index calculated by ${N_{{\textrm{T}_\textrm{X}}}} = 20$ and ${N_{{\textrm{T}_\textrm{X}}}} = 35$ is ${10^{ - 4}};$ by increasing the transverse nodal number of the computational window to ${N_{{\textrm{T}_\textrm{X}}}} = 30,$ the order of difference is reduced to ${10^{ - 5}}.$ Therefore, the FC method is more efficient and accurate because it only used $35 \times 24 = 840$ nodes to achieve an accuracy of the order of ${10^{ - 5}}.$ Fig. 6(b) shows the effective indexes of the first- and second-order modes for a magneto-optical raised strip waveguide with respect to the first-order magneto-optic effect. It should be mentioned that the transverse and vertical nodal numbers of the computational window are taken as ${N_{{\textrm{T}_\textrm{X}}}} = 35$ and ${N_{{\textrm{T}_\textrm{Y}}}} = 24,$ respectively. In the core layer, the node interval is taken as $0.065\mu m$ and the nodal number is 16, correspondingly, the node interval used in the rest clouds is set to be $0.1\mu m$ and the nodal number is 19. In order to examine the accuracy of the present FC method, results from the full-vectorial FD method and the full-vectorial MPS method described in [27], where the terms of the basis functions for expanding the exterior and interior subdomains are, respectively, fixed at 10 and 20, are plotted in the same figure. It can be observed that the results of the first-order mode from the FC method have noticeable differences with those from the FD method as $\zeta $ increases, whereas that the differences of the results of the second-order mode from the FC method versus those of the FD method decrease with the increasing of the condition. In addition, the meshed points used in the FD method are $320 \times 350$ over a computational window of $3.2\mathrm{\mu }\textrm{m} \times 2.9\mathrm{\mu }\textrm{m}.$ Moreover, for both two modes, results from the FC method are in good agreement with those from the MPS method. However, to obtain the convergent solutions, the FC method only needs $35 \times 24 = 840$ computational nodes, while the MPS method needs nearly $45 \times 45 = 2025$ ones. Moreover, to obtain the effective index of the field patterns of the first- and second-order modes, the computational time of the FC method and the MPS method is 10.5 and 13.4 seconds, respectively. Comparison of the different result validates that the computational efficiency of the FC method is superior to that of the MPS method. Therefore, the high numerical accuracy and relatively low computational burden of the FC method exhibit its robustness and versatility. For the condition $\zeta = 0.005,$ Fig. 7 shows the field patterns of the two lowest-order modes for a magneto-optical raised strip waveguide. It can be found that the $|{{H_x}} |$ and $|{{H_y}} |$ fields for both modes have the amplitudes of the same order. One can also see that, for the two lowest-order modes, the profiles and magnitudes of the $|{{H_x}} |$ field are almost equal to those of the $|{{H_y}} |$ field.

 figure: Fig. 6.

Fig. 6. Effective indexes of the first- and second-order modes for a magneto-optical raised strip waveguide with respect to the transverse nodal number and the first-order magneto-opic effect: (a)${N_{{\textrm{T}_\textrm{X}}}}$ and (b) $\zeta .$

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Field patterns of the two lowest-order modes for a magneto-optical raised strip waveguide for $\zeta = 0.005:$(a) $|{{H_x}} |$ and (b) $|{{H_y}} |$ for the first-order mode; (c) $|{{H_x}} |$ and (d) $|{{H_y}} |$ for the second-order mode.

Download Full Size | PDF

Currently, in integrated optics, optical devices based on the liquid crystals (LCs) are popular, mainly owing to their birefringence and electro-optical properties. Among them, the nematic LC channel waveguide has been widely used to form various devices for its low manufacturing cost and high electro-optic effect. This will be followed by a further modal analysis of the nematic LC channel waveguide based on the FC method. Figure 8 shows the cross-section of the nematic LC channel waveguide, which is made of a core region filled with nematic LC, embedded in a glass substrate and covered by the air. The computational parameters are: ${n_\textrm{s}} = 1.45,$ ${n_\textrm{a}} = 1.0,$ $w = t = 3\mathrm{\mu }\textrm{m},$ $\lambda = 1.55\mathrm{\mu }\textrm{m}$ and $d = 0.5\mathrm{\mu }\textrm{m}.$ The computational window is set to be $2.5\mathrm{\mu }\textrm{m} \times 2.1\mathrm{\mu }\textrm{m}.$ In addition, Fig. 8 also presents the twist angle $\varphi $ of the LC among the x-axis and the direction ${n_\mathrm{\gamma }}.$ Note that the value of $\varphi $ can be altered by imposing the voltage which is applied by the designed electrodes. In this work, for simplicity, the twist angle $\varphi $ is taken as ${0^ \circ },$ ${30^ \circ },$ ${45^ \circ },$ ${60^ \circ },$ and ${90^ \circ }.$

 figure: Fig. 8.

Fig. 8. Cross section of a nematic liquid crystal (LC) channel waveguide surrounded by the PMLs.

Download Full Size | PDF

The ordinary and extraordinary refractive indexes of the nematic LC are, respectively, set to be ${n_\textrm{e}} = 1.7072$ and ${n_\textrm{o}} = 1.5292.$ The relative permittivity tensor of the core is given by

$$\tilde{\varepsilon } = \left[ {\begin{array}{ccc} {n_\textrm{o}^2 + \Delta \varepsilon {{\cos }^2}\theta }&{\Delta \varepsilon \cos \theta \sin \theta }&0\\ {\Delta \varepsilon \cos \theta \sin \theta }&{n_\textrm{o}^2 + \Delta \varepsilon {{\sin }^2}\theta }&0\\ 0&0&{n_\textrm{o}^2} \end{array}} \right],\Delta \varepsilon = n_\textrm{e}^2 - n_\textrm{o}^2.$$

Firstly, Fig. 9 presents the effective indexes versus the transverse nodal number for the first five modes of the nematic LC channel waveguide with $\varphi = {30^ \circ }$ and $\varphi = {60^ \circ }.$ As is clearly seen that, the effective indexes decrease as the number of mode increases for the two twist angles. It also can be observed that the effective indexes of the first five modes linearly decrease as the transverse nodal numbers increase. In addition, the order of the difference between the effective indexes calculated by ${N_{{\textrm{T}_\textrm{X}}}} = 30$ and ${N_{{\textrm{T}_\textrm{X}}}} = 35$ is ${10^{ - 5}}$ for both five modes. Therefore, 30 transverse nodes and 24 vertical nodes are used in the following analysis. In the core layer, the node interval is taken as $0.1\mu m$ and the nodal number is 10, correspondingly, the node interval used in the rest clouds is set to be $0.1\mu m$ and the nodal number is 20.

 figure: Fig. 9.

Fig. 9. Effective indexes versus the transverse nodal number for the first five modes of the nematic LC channel waveguide with: (a)$\varphi = {30^ \circ }$ and (b) $\varphi = {60^ \circ }$

Download Full Size | PDF

Tables 1 and 2 show the effective indexes at various values of the twist angle of the first five modes for the nematic LC channel waveguide. In order to make the comparison more intuitive, solutions from the MP method are also presented in the same tables. Numerical investigation determined that the FC and MPS methods have the nearly identical convergent behaviors since comparison of the different results shows a high degree of agreement for the two methods. It also can be seen that the FC method converge faster than the MPS method since the former with only $30 \times 24 = 720$ computational nodes can obtain more accurate results than the latter with more than $38 \times 38 = 1440$ grid points. On the other hand, to obtain the effective index of the field patterns of the mode 1, 2 and 5, the computational time of the FC method is 13.6 seconds, respectively. Therefore, the efficiency of the FC method is superior to the MPS method. Furthermore, to verify the variation of the effective index versus the twist angle, Fig. 10 shows the mode patterns of the mode 1, 2 and 5 of a nematic LC channel waveguide for four twist angles. Note that the symbols $|{{H_{nx}}} |$ and $|{{H_{ny}}} |,(n = 1,2,5)$ represent the corresponding mode fields. It can be seen that the variations in the effective indexes of the fundamental guided mode (mode 1) are fairly small for different values of the twist angles, which can be verified by the major field patterns of mode 1 for four twist angles. However, the variations of the effective indexes versus the twist angle show the opposite behavior for high-order modes as shown in Fig. 10. It is clearly observed that the field distributions corresponding to the y-axis are symmetric when the twist angles are taken as ${0^ \circ }$ and ${90^ \circ }.$ However, different from the above two values of the twist angles; more complicated symmetries of the other values may not exist in high-order modes. In contrast, for $\varphi = {45^ \circ },$ the field values of the two transverse magnetic components of the guided modes are almost the same due to the two polarizations experienced the same propagation in the core region. In present work, for brevity, the field patterns for $\varphi = {45^ \circ }$ are not presented.

 figure: Fig. 10.

Fig. 10. $|{{H_x}} |$ and $|{{H_y}} |$ mode patterns of the mode 1, 2 and 5 of a nematic LC channel waveguide for four twist angles.

Download Full Size | PDF

Tables Icon

Table 1. Effective indexes at various values of the twist angle of the first five modes for a nematic LC channel waveguide computed by different methods

View Table | View all tables in this article
Tables Icon

Table 2. Effective indexes at various values of the twist angle of the first five modes for a nematic LC channel waveguide computed by different methods

View Table | View all tables in this article

Finally, the FC method, which bases on irregular nodal distribution in each cloud, can be used to solve the anisotropic optical waveguides with non-square geometries, including the triangle and curvilinear geometries. The major steps for dealing with this work are shown as follows. The first step is the coordinate transformation of nodes and wave equations. Owing to irregular waveguide geometries, extra coordinate transformation is indispensable for dealing with different parts, including the construction of interpolated nodes, the derivation of the full-vector wave equations and the discretization of the equations based on the FC method. The commonly used techniques in the coordinate transformation are the coordinate mapping techniques, including triangle and curvilinear coordinate mapping techniques. The second step is the construction of the nodes in each cloud. In this work, different from the waveguides with square geometries, the irregular nodal distribution in each cloud is used. To do this, a new FC approximation scheme is utilized. As a general rule, this approximation falls into three kinds of steps, which are shown as follows. The first one is the refinement of cloud. Each cloud is divided into several sub-clouds, and each sub-cloud can be chosen as circular, rectangular, triangle and trapezoid. The second one is the construction of the nodes in each sub-cloud. Among these sub-clouds, the nodal distributions are chosen in different regular manners, in which the kernel node is limited only one and further the locations of the kernel nodes in each cloud are fixed at the center. The third one is selecting a new kernel node and removing the sub-clouds. Among these selected kernel nodes, a new kernel node, namely, the so-called star node, should be selected. After the star node has been selected, the sub-clouds are removed but the nodes of each sub-cloud are retained. It should be noted that the above approximation is only suitable for interior nodal distributions. In order to deal with boundaries among clouds, an extra set of nodes should be added. In this work, it is remarkable that, in the guiding layer, the interior nodes increase exponentially from the top to the bottom of the cloud and the extra boundary nodes are transferred to Cartesian coordinate system via the coordinate mapping technique. The third step is the treatment of boundaries. By utilizing the coordinate mapping technique, the continuity conditions of the longitudinal field components used in Eqs. (30-33) can be expanded to non-square waveguide geometries. The final step is the carefully selections of several factors used in the discretization of the objective equations by using the FC method, including the nodal volume, the small number for avoiding the singularity of the kernel function, the support size of the cloud and the dilation parameter. Due to page limitation, more detailed works and results about the transverse anisotropic waveguides with non-square geometries will be shown in the future.

4. Conclusions

A full-vector mode solver for transverse anisotropic optical waveguides by using the FC method in terms of transverse magnetic field components is established, where the robust PML absorbing boundary conditions are adopted for effectively avoiding the nonphysical reflection from the computational window edges. The full-vector wave equations are discretized using a point collocation method and then with the help of the fixed reproducing kernel technique and the continuity conditions of the longitudinal field components, a standard matrix eigenvalue equation is obtained. The traditional numerical methods, including the MC and MPS methods, are usually based on background mesh, which may increase the computational complexity and then strongly affect the computational efficiency. In addition, the computational domain in such two schemes is partitioned into a number of subdomains according to its geometry, and then the mode field profiles in different subdomains are expanded by imposing distinct basis functions, which may increase the computational complexity. The FC method requires a set of scattered nodes without the connectivity information among nodes and additional efforts for choosing the optimum scaling factor are avoided since only one kind of basis functions is utilized. Therefore, the FC method is more efficient since it can achieve the equivalent convergence with much less computational effort. Moreover, the proposed method uses a modified Gaussian function to construct the meshless interpolation function, significantly improving the accuracy of the mode solver. The numerical results for a typical anisotropic square waveguide, a magneto-optical raised strip waveguide and a nematic LC channel waveguide indicate that the FC method shows superior convergent behavior to the FD, MC and MP methods, with high accuracy and efficiency. The proposed FC method can be used in the future to establish a more efficient mode solver for more complex optical dielectric waveguides with full ${3 \times 3}$ anisotropy.

Funding

Natural Science Foundation of Jiangsu Province (BK20211163); National Natural Science Foundation of China (11574046).

Disclosures

The authors declare no conflicts of interest.

Data availability

No data were generated or analyzed in the presented research.

References

1. J. W. Wu, “Birefringent and electro-optic effects in poled polymer films: steady-state and transient properties,” J. Opt. Soc. Am. B 8(1), 142–152 (1991). [CrossRef]  

2. K. Srinivasan and B. J. H. Stadler, “Magneto-optical materials and designs for integrated TE- and TM-mode planar waveguide isolators: a review,” Opt. Mater. Express 8(11), 3307–3318 (2018). [CrossRef]  

3. U. Hempelmann, H. Herrmann, G. Mrozynski, V. Reimann, and W. Sohler, “Integrated optical proton exchanged TM-pass polarizers in LiNbO3: modeling and experimental performance,” J. Lightwave Technol. 13(8), 1750–1759 (1995). [CrossRef]  

4. N. Anwar, C. Themistos, B. M. A. Rahman, and K. T. V. Grattan, “Design consideration for an electrooptic directional coupler modulator,” J. Lightwave Technol. 17(4), 598–605 (1999). [CrossRef]  

5. P. Daggar and S. S. Verma, “Verifying faraday’s magneto-optical effect for some materials,” Eur. J. Biochem. 7(1), 8–14 (2019). [CrossRef]  

6. R. Wolfe, R. A. Lieberman, V. J. Fratello, R. E. Scotti, and N. Kopylov, “Etch-tuned ridged waveguide magneto-optic isolator,” Appl. Phys. Lett. 56(5), 426–428 (1990). [CrossRef]  

7. O. Zhuromoskyy, M. Lohmeyer, N. Bahlmann, H. Dotsch, P. Hertel, and A. F. Popkov, “Analysis of polarization independent Mach-Zehnder-type integrated optical isolator,” J. Lightwave Technol. 17(7), 1200–1205 (1999). [CrossRef]  

8. W. Zaets and K. Ando, “Optical waveguide isolator based on nonreciprocal loss/gain of amplifier converted by ferromagnetic layers,” IEEE Photonics Technol. Lett. 11(8), 1012–1014 (1999). [CrossRef]  

9. P. Lusse, K. Ramm, and H. G. Unger, “Vectorial eigenmodes calculation for anisotropic planar optical waveguides,” Electron. Lett. 32(1), 38–39 (1996). [CrossRef]  

10. C. L. da Silva Souza Sobrinho and A. Giarola, “Analysis of biaxially anisotropic dielectric waveguides with Gaussian-Gaussian index or refraction profiles by the finite-difference method,” Proc. Inst. Elect. Eng. 140(3), 224–230 (1993). [CrossRef]  

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

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

13. 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]  

14. Y. Lu and F. A. Fernandez, “An efficient finite element solution of inhomogeneous anisotropic and lossy dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 41(6), 1215–1223 (1993). [CrossRef]  

15. K. Hayata, K. Miura, and M. Koshiba, “Full vectorial finite element formalism for lossy anisotropic waveguides,” IEEE Trans. Microwave Theory Tech. 37(5), 875–883 (1989). [CrossRef]  

16. M. Y. Chen and H. C. Chang, “Finite-difference frequency-domain analysis of nematic liquid crystal optical waveguides in silicon V-grooves,” 2009 Conference on Lasers & Electro Optics & The Pacific Rim Conference on Lasers and Electro-Optics.

17. 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 an arbitrary permittivity tensor,” Opt. Express 17(8), 5965–5979 (2009). [CrossRef]  

18. X. L. Li, H. Wang, Z. H. Yang, and B. Li, “A new implicit hybridizable discontinuous Galerkin time-domain method for solving the 3-D electromagnetic problems,” Appl. Math. Lett. 93, 124–130 (2019). [CrossRef]  

19. J. Xiao, M. Zhang, and X. Sun, “Full-vectorial analysis of optical waveguides by the mapped Galerkin method based on H-field,” 2005 IEEE International Symposium on Microwave, Antenna, Propagation and EMC Technologies for Wireless Communications.

20. 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]  

21. C. Y. Wang, H. H. Liu, S. Y. Chung, C. H. Teng, C. P. Chen, and H. C. Chang, “High-accuracy waveguide leaky-mode analysis using a multidomain pseudospectral frequency-domain method incorporated with stretched coordinate PML,” J. Lightwave Technol. 31(14), 2347–2360 (2013). [CrossRef]  

22. F. Alharbi, “Full-vectoral meshfree spectral method for optical-waveguide analysis,” IEEE Photonics J. 5(1), 6600315 (2013). [CrossRef]  

23. 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]  

24. J. Xiao, H. B. Ma, N. F. Bai, X. Liu, and X. Sun, “Full-vectorial mode solver for bending waveguides using multidomain pseudospectral method in a cylindrical coordinate system,” IEEE Photonics Technol. Lett. 21(23), 1779–1781 (2009). [CrossRef]  

25. X. Z. Jin, G. Li, and N. R. Aluru, “New approximations and collocation schemes in the finite cloud method,” Comput. Struct. 83(17-18), 1366–1385 (2005). [CrossRef]  

26. P. J. Chiang and H. C. Chang, “A high-accuracy pseudospectral full-vectorial leaky optical waveguide mode solver with carefully implemented UPML absorbing boundary conditions,” Opt. Express 19(2), 1594–1608 (2011). [CrossRef]  

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

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

29. C. C. Huang, “Solving the full anisotropic liquid crystal waveguides by using an iterative pseudospectral-based eigenvalue method,” Opt. Express 19(4), 3363–3378 (2011). [CrossRef]  

30. D. R. Burke and T. J. Smy, “A meshless based solution to vectorial mode fields in optical micro-structured waveguides using leaky boundary conditions,” J. Lightwave Technol. 31(8), 1191–1197 (2013). [CrossRef]  

31. N. R. Aluru and G. Li, “Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation,” Int. J. Numer. Meth. Engng. 50(10), 2373–2410 (2001). [CrossRef]  

32. D. R. Burke and T. J. Smy, “Optical mode solving for complex waveguides using a finite cloud method,” Opt. Express 20(16), 17783–17796 (2012). [CrossRef]  

33. N. R. Aluru, “A point collocation method based on reproducing kernel approximations,” Int. J. Numer. Meth. Engng. 47(6), 1083–1121 (2000). [CrossRef]  

Data availability

No data were generated or analyzed in the presented research.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)
Fig. 1.
Fig. 1. Interfaces of dielectric materials: (a) Horizontal and (b) Vertical interfaces.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Cross section of an anisotropic square waveguide enclosed by the PMLs.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Effective indexes of the first- and second-order modes for an anisotropic square waveguide versus the transverse nodal number and the core birefringence: (a)${N_{{\textrm{T}_\textrm{X}}}}$ and (b)$\Delta n.$

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Distributions of field intensity of the two lowest-order modes for an anisotropic square waveguide: (a)${H_x}$ and (b)${H_y}$ for the first-order mode; (c)${H_x}$ and (d) ${H_y}$ for the second-order mode.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. Cross section of a magneto-optical raised strip waveguide surrounded by the PMLs

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Effective indexes of the first- and second-order modes for a magneto-optical raised strip waveguide with respect to the transverse nodal number and the first-order magneto-opic effect: (a)${N_{{\textrm{T}_\textrm{X}}}}$ and (b) $\zeta .$

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Field patterns of the two lowest-order modes for a magneto-optical raised strip waveguide for $\zeta = 0.005:$(a) $|{{H_x}} |$ and (b) $|{{H_y}} |$ for the first-order mode; (c) $|{{H_x}} |$ and (d) $|{{H_y}} |$ for the second-order mode.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. Cross section of a nematic liquid crystal (LC) channel waveguide surrounded by the PMLs.

View in Article | Download Full Size | PDF

Fig. 9.
Fig. 9. Effective indexes versus the transverse nodal number for the first five modes of the nematic LC channel waveguide with: (a)$\varphi = {30^ \circ }$ and (b) $\varphi = {60^ \circ }$

View in Article | Download Full Size | PDF

Fig. 10.
Fig. 10. $|{{H_x}} |$ and $|{{H_y}} |$ mode patterns of the mode 1, 2 and 5 of a nematic LC channel waveguide for four twist angles.

View in Article | Download Full Size | PDF

Tables (2)

Tables Icon

Table 1. Effective indexes at various values of the twist angle of the first five modes for a nematic LC channel waveguide computed by different methods

View Table | View all tables in this article
Tables Icon

Table 2. Effective indexes at various values of the twist angle of the first five modes for a nematic LC channel waveguide computed by different methods

View Table | View all tables in this article

Equations (43)

Equations on this page are rendered with MathJax. Learn more.

ξ = ∫ 0 ξ s ξ ( ξ ′ ) d ξ ′
s ξ = { 1 − j τ ξ ( ξ ) / ω ε 0 ,  in the PML region 1 ,  in the non - PML region
τ ξ = m + 1 150 π Δ ξ | ξ − ξ 0 | d 2 ,
∇ ~ × ( ε ~ − 1 ∇ ~ × H → ) − ω 2 μ 0 H → = 0 ,
∇ ~  =  α x ∂ ∂ x x ~ + α y ∂ ∂ y y ~ + α z ∂ ∂ z z ~ ,
ε ~ = [ ε x x ε x y 0 ε y x ε y y 0 0 0 ε z z ] .
[ A x x A x y A y x A y y ] [ H x H y ] = β 2 [ H x H y ] .
A x x = ∂ 2 ( I ⋅ ) ∂ x ~ 2 + ε y y ∂ ∂ y ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ y ~ ) + ε y x ∂ ∂ y ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ x ~ ) + ω 2 μ 0 ε y y ,
A x y = ∂ 2 ( I ⋅ ) ∂ y ~ ∂ x ~ − ε y y ∂ ∂ y ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ x ~ ) − ε y x ∂ ∂ x ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ x ~ ) − ω 2 μ 0 ε y x , A y x = ∂ 2 ( I ⋅ ) ∂ x ~ ∂ y ~ − ε x x ∂ ∂ x ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ y ~ ) − ε x y ∂ ∂ y ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ y ~ ) − ω 2 μ 0 ε x y , A y y = ∂ 2 ( I ⋅ ) ∂ y ~ 2 + ε x x ∂ ∂ x ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ x ~ ) + ε x y ∂ ∂ x ~ ( 1 ε z z ∂ ( I ⋅ ) ∂ y ~ ) + ω 2 μ 0 ε x x ,
H x i r ( x ~ r i , y ~ r i )  =  ∫ Ω ζ ( x ~ r i , y ~ r i , x ~ r k − s , y ~ r k − t ) ϕ ( x ~ r k , y ~ r k ) H x i u r ( x ~ r i , y ~ r i ) d s d t ,
ζ ( x ~ r i , y ~ r i , x ~ r k − s , y ~ r k − t ) = P T ( x r k − s , y ~ r k − t ) c ( x ~ r i , y ~ r i ) ,
P T ( x ~ r k − s , y ~ r k − t ) = [ 1 , x ~ r k − s , y ~ r k − t , ( x ~ r k − s ) 2 , ( x ~ r k − s ) ( y ~ r k − t ) , ( y ~ r k − t ) 2 ] .
∫ Ω P T ( x ~ r i , y ~ r i ) c ( x ~ r i , y ~ r i ) ϕ ( x ~ r k − s , y ~ r k − t ) p n ( s , t ) d s d t = p n ( x ~ r i , y ~ r i ) , n = 1 , 2 , ⋯ , 6.
∑ I = 1 N P P T ( x ~ r i , y ~ r i ) c ( x ~ r i , y ~ r i ) ϕ ( x ~ r k − x ~ r I , y ~ r k − y ~ r I ) p n ( x ~ r I , y ~ r I ) Δ V r I = p n ( x ~ r i , y ~ r i ) ,
M → c ( x ~ r i , y ~ r i )  =  p ( x ~ r i , y ~ r i ) ,
M j l = ∑ I = 1 N P { p j ( x ~ r I , y ~ r I ) ζ ( x ~ r k − x ~ r I , y ~ r k − y ~ r I ) p l ( x ~ r I , y ~ r I ) Δ V I } , j , l = 1 , 2 , ⋯ , 6.
c ( x ~ r i , y ~ r i )  =  M → − 1 p ( x ~ r i , y ~ r i ) .
ζ ( x ~ r i , y ~ r i , x ~ r k − s , y ~ r k − t )  =  P T ( x ~ r i , y ~ r i ) M → − 1 p ( x ~ r i , y ~ r i ) .
ϕ ( x ~ r k − x ~ r I ) = ω ( x ~ r k − x ~ r I ) 1 − ω ( x ~ r k − x ~ r I ) + ε ′ ,
ω ( x ~ r k − x ~ r I ) = { 0 , t > d r e − ( d r ( x ~ i k − x ~ I ) / 2 ) 2 − e − 4 1 − e − 4 , t ≤ d r ,
ϕ ( x ~ r k − x ~ r I , y ~ r k − y ~ r I ) = ϕ ( x ~ r k − x ~ r I ) ϕ ( y ~ r k − y ~ r I ) .
H x i r ( x ~ r i , y ~ r i ) = ∑ I = 1 N P N r I ( x ~ r i , y ~ r i ) H x I r ,
H y i r ( x ~ r i , y ~ r i ) = ∑ I = 1 N P N r I ( x ~ r i , y ~ r i ) H y I r ,
N r I ( x ~ r i , y ~ r i ) = P T ( x ~ r i , y ~ r i ) M → − 1 p ( x ~ r i , y ~ r i ) ζ ( x ~ r i , y ~ r i , x ~ r k − s , y ~ r k − t ) Δ V r I .
[ N x x r N x y r N y x r N y y r ] [ H x r H y r ] = β 2 [ H x r H y r ] ,
N x x r H x r = [ ∂ 2 ( I ⋅ ) ∂ x ~ 2 + ε y y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ x ~ + ω 2 μ 0 ε y y r ] H x r | x ~ = x ~ r i , y ~ = y ~ r i = ∑ i = 1 N P ∑ j = 1 N P { ∑ l = 1 N P [ ∂ 2 ( I ⋅ ) ∂ x ~ 2 + ε y y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ x ~ + ω 2 μ 0 ε y y r ] N r I [ H x l r ] } | x ~ = x ~ r i , y ~ = y ~ r i ,
N x y r H y r = [ ∂ 2 ( I ⋅ ) ∂ y ~ ∂ x ~ − ε y y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ x ~ − ε y x r ε z z r ∂ ∂ x ~ ∂ ( I ⋅ ) ∂ x ~ − ω 2 μ 0 ε y x r ] H y r | x ~ = x ~ r i , y ~ = y ~ r i
= ∑ i = 1 N P ∑ j = 1 N P { ∑ l = 1 N P [ ∂ 2 ( I ⋅ ) ∂ y ~ ∂ x ~ − ε y y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ x ~ − ε y x r ε z z r ∂ ∂ x ~ ∂ ( I ⋅ ) ∂ x ~ − ω 2 μ 0 ε y x r ] N r I [ H y l r ] } | x ~ = x ~ r i , y ~ = y ~ r i ,
N y x r H x r = [ ∂ 2 ( I ⋅ ) ∂ x ~ ∂ y ~ − ε x x r ε z z r ∂ ∂ x ~ ∂ ( I ⋅ ) ∂ y ~ − ε x y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ y ~ − ω 2 μ 0 ε x y r ] H x r | x ~ = x ~ r i , y ~ = y ~ r i
= ∑ i = 1 N P ∑ j = 1 N P { ∑ l = 1 N P [ ∂ 2 ( I ⋅ ) ∂ x ~ ∂ y ~ − ε x x r ε z z r ∂ ∂ x ~ ∂ ( I ⋅ ) ∂ y ~ − ε x y r ε z z r ∂ ∂ y ~ ∂ ( I ⋅ ) ∂ y ~ − ω 2 μ 0 ε x y r ] N r I [ H x l r ] } | x ~ = x ~ r i , y ~ = y ~ r i ,
N x y r H y r = [ ∂ 2 ( I ⋅ ) ∂ y ~ 2 + ε x x r ε z z r ∂ ∂ x ~ 1 ε z z ∂ ( I ⋅ ) ∂ x ~ + ε x y r ε z z r ∂ ∂ x ~ 1 ε z z ∂ ( I ⋅ ) ∂ y ~ + ω 2 μ 0 ε x x r ] H y r | x ~ = x ~ r i , y ~ = y ~ r i
= ∑ i = 1 N P ∑ j = 1 N P { ∑ l = 1 N P [ ∂ 2 ( I ⋅ ) ∂ y ~ 2 + ε x x r ε z z r ∂ ∂ x ~ 1 ε z z ∂ ( I ⋅ ) ∂ x ~ + ε x y r ε z z r ∂ ∂ x ~ 1 ε z z ∂ ( I ⋅ ) ∂ y ~ + ω 2 μ 0 ε x x r ] N r I [ H y l r ] } | x ~ = x ~ r i , y ~ = y ~ r i .
[ [ N 1 ] ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ [ N m ] ] [ [ H 1 ] ⋮ [ H m ] ] = β 2 [ [ H 1 ] ⋮ [ H m ] ] ,
[ N r ] = [ N x x r N x y r N y x r N y y r ] , [ H r ] = [ H x r H y r ] , r = 1 , 2 , ⋯ , m .
H Z = − j β ( ∂ H x ∂ x ~ + ∂ H y ∂ y ~ ) ,
E Z = − j ω ε z z ( ∂ H x ∂ y ~ − ∂ H y ∂ x ~ ) .
∂ H y ∂ y ~ | y + = ∂ H y ∂ y ~ | y −
ε z z y + ∂ H x ∂ y ~ | y − − ε z z y − ∂ H x ∂ y ~ | y + = ( ε z z y + − ε z z y − ) ∂ H y ∂ x ~ ,
∂ H x ∂ x ~ | x + = ∂ H x ∂ x ~ | x −
ε z z x + ∂ H y ∂ x ~ | x − − ε z z x − ∂ H y ∂ x ~ | x + = ( ε z z x + − ε z z x − ) ∂ H x ∂ y ~ ,
ε ~ = [ n o 2 + ( n e 2 − n o 2 ) cos 2 θ ( n e 2 − n o 2 ) cos ⁡ θ sin ⁡ θ 0 ( n e 2 − n o 2 ) cos ⁡ θ sin ⁡ θ n o 2 + ( n e 2 − n o 2 ) sin 2 θ 0 0 0 n o 2 ]
ε ~ = [ n f 2 j ζ 0 − j ζ n f 2 0 0 0 n f 2 ]
ε ~ = [ n o 2 + Δ ε cos 2 θ Δ ε cos ⁡ θ sin ⁡ θ 0 Δ ε cos ⁡ θ sin ⁡ θ n o 2 + Δ ε sin 2 θ 0 0 0 n o 2 ] , Δ ε = n e 2 − n o 2 .
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use