Open Access
Add to BibSonomy
Abstract
A three-point finite difference method with an arbitrary order of accuracy is proposed for the modal analysis of chiral planar waveguides. The elaborate application of a nonuniform grid, compact finite difference technique, and boundary conditions results in an efficient, easily implemented, and versatile tool for the modal analysis of chiral planar waveguides with an arbitrarily discontinuous profile of permittivity, permeability, and chirality. In particular, this method efficiently resolves the fine structures in plasmon and photonic crystal waveguides. For the test model of a chiral-metallic plasmon waveguide, stable convergence up to a sixteenth order of accuracy can be obtained, which produces a relative error on the effective index that approaches the machine precision with only eighty grid points.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Recently, renewed interest has arisen in chiral planar waveguides that consist of optically active media owing to their potential applications in chiral detection [1–3] and separation [4,5]. These applications take advantage of the enhanced light-matter interaction of surface plasmon polaritons (SPPs) and rely on the unique properties of guided modes in waveguide structures. Hence, modal analysis is a basic procedure in device design and detection processing. At present, the modal analysis of chiral plasmon waveguides is based on analytic characteristic equations and limited to simple double-layered [1] and three-layered [2–4] cases. A lack of efficient numerical tools hinders the analysis of more interesting chiral multilayered waveguides.
In the past few decades, the treatment of discontinuities in refractive index profiles and the search for high truncation order schemes have been two main objectives in the development of the finite difference (FD) method, which is popular for its easy implementation. Following Stern’s pioneering work [6–16], researchers have developed arbitrary truncation order FD methods for the modal analysis of both one-dimensional [13,14] and two-dimensional [15] waveguides. Sujecki [16] built a remarkable three-point scheme with an arbitrary truncation order for planar waveguides using the compact FD technique [17].
However, the high-order numerical methods mentioned above are restricted to filling materials without chirality. As for chiral waveguides with an arbitrarily shaped cross section, finite element methods that enforce boundary conditions without a higher-order extension have been developed [18,19]. For chiral circular optical fibers, a mode solver with a fourth-order of accuracy has been proposed [20], utilizing the immersed interface method (IIM) and a five-point stencil. Reference [21] presents a second-order FD method for the modal analysis of chiral planar waveguides that is also based on the IIM. Although this method may be used to simulate chiral plasmon waveguides, the efficiency is limited by its low order of accuracy. According to Ref. [20], in principle, the IIM is capable of generating arbitrary higher-order schemes; however, the complexity of formulation appears to rapidly increase as the order of accuracy increases.
In this paper, an arbitrary-order three-point FD method is reported for the modal analysis of chiral planar waveguides with an arbitrarily discontinuous profile of permittivity, permeability, and chirality. In this method, several tricks are used to improve simplicity, efficiency, and versatility. First, the FD scheme is directly built on a nonuniform grid, which is a strategy rarely employed by previous researchers, especially in the construction of high-order schemes. However, for planar waveguides, a nonuniform grid has a crucial advantage over a uniform grid. A nonuniform grid not only leads to a simpler FD scheme for irregular grid points in the presence of interfaces, but also provides a flexible allocation of grid points to adapt the fine structures encountered in field patterns or waveguide structures, which usually accelerates the convergence. Second, the compact FD technique is applied to obtain a three-point scheme with an arbitrary order of accuracy, which guarantees rapid convergence and simplified formulation and implementation. Third, the procedure of enforcing the boundary conditions and applying the compact FD technique is elaborately designed such that only the jump conditions of the field derivatives up to second order are required to construct arbitrary-order schemes. This final trick significantly simplifies the formulation and has not been previously reported according to the author’s knowledge.
These tricks are coherent, resulting in an efficient, easily implemented, and versatile tool for the modal analysis of chiral planar waveguides with an arbitrary profile of electromagnetic parameters. Because the chiral medium (characterized by $\mathbf {D}=\epsilon \mathbf {E}+i\xi \mathbf {H}$ and $\mathbf {B}=\mu \mathbf {H}-i\xi \mathbf {E}$) includes the isotropic medium (achiral, characterized by $\mathbf {D}=\epsilon \mathbf {E}$ and $\mathbf {B}=\mu \mathbf {H}$) as its special case of null chirality ($\xi =0$), one formulation set that retains all of the above advantages is immediately acquired for common achiral planar waveguides with an arbitrary refractive index profile if the chirality is allowed to vanish.
This paper is organized as follows: In Section 2, the formulation is described in detail, where the derivation of an arbitrary-order FD approximation of the first and second-order derivatives of the fields on a nonuniform grid makes up the core content. In Section 3, a numerical example of a multi-quantum-well waveguide without chirality is first analyzed to demonstrate the more flexible and efficient features of this method compared to the method based on a uniform grid. Then, two numerical examples of a chiral-metallic plasmon waveguide and chiral photonic crystal (PC) waveguide supporting surface waves are presented to validate this method for chiral waveguides. A summary is given in Section 4.
2. Formulation
2.1 Modal eigenvalue problem and jump conditions
For a chiral medium with the constitutive relations $\mathbf {D}=\epsilon \mathbf {E}+i\xi \mathbf {H}$, and $\mathbf {B}=\mu \mathbf {H}-i\xi \mathbf {E}$, the normalized Maxwell’s curl equations read as
where $k_0=\omega \sqrt {\epsilon _0\mu _0}$ is the vacuum wave number, $\epsilon _r=\epsilon /\epsilon _0$, $\mu _r=\mu /\mu _0$, and $\xi _r=\xi /\sqrt {\epsilon _0\mu _0}$ denote the relative permittivity, permeability, and chirality, respectively, and the magnetic field $\mathbf {H}$ is normalized as $\mathbf {F}=i\sqrt {\mu _0/\epsilon _0}\mathbf {H}$. A general chiral planar multilayered waveguide is considered, each layer of which can be filled with a chiral medium. The $x$-axis is laid normal to the interfaces, and modal propagation along the $z$-direction, that is, the propagation term $\exp \left [i\left (\beta z-\omega t\right )\right ]$ for the fields, is assumed. Through this setup of the coordinate axis system, the fields are suggested as independent of the $y$-coordinate. The matrix Helmholtz equation expressed as an eigenvalue problem, can be derived for vector $\mathbf {g}=[E_y;F_y]$, where $\bar {x}=k_0x$ is the normalized $x$-coordinate, matrix $\mathbf {N}=[\xi _r,\mu _r;\epsilon _r,\xi _r]$ is a generalization of the refractive index of an achiral medium, and $n_\mathrm {eff}=\beta /k_0$ is the effective index of the mode. When chirality vanishes, namely $\xi _r=0$, Eq. (3) decouples into two scalar Helmholtz equations for the ordinary transverse electric (TE) and transverse magnetic (TM) modes, whereas in the chiral case, the TE and TM fields are coupled to form a hybrid mode.Unlike previous research, in which jump conditions of higher-order derivatives of field $\phi$ are used to generate higher-order FD approximations of $\phi ^{(2)}$, in this study, only the jump conditions of $\mathbf {g}$, $\mathbf {g}^{(1)}$, and $\mathbf {g}^{(2)}$ are required to approximate $\mathbf {g}^{(2)}$ to an arbitrary order of accuracy; this is achieved by taking full advantage of the compact FD technique. The jump conditions of $\mathbf {g}$ and $\mathbf {g}^{(1)}$ directly originate from the physical boundary conditions expressed as
where $\left [f\right ]=0$ indicates that $f$ is continuous across the interface. Combining Eqs. (3) and (4), the jump condition of $\mathbf {g}^{(2)}$ is derived as2.2 Accurate expressions for first and second-order derivatives
Making use of the flexible allocation of grid points owing to the nonuniform grid, a grid point can be placed at the interface, as shown in Fig. 1. Physical quantities on the left and right sides of the interface are labelled with subscripts ‘$-$’ and ‘$+$’, respectively. Grid point $i$ is observed on the left side of interface, namely,
Because Eqs. (4) and (8) imply $\mathbf {g}_i=\mathbf {g}_-=\mathbf {g}_+$, grid point $i+1$ is effectively regular in the three-point stencil; in other words, no special treatment with respect to the interface is required for grid point $i+1$. Only grid point $i$ is irregular, at which the FD approximation of $\mathbf {g}^{(2)}$ must be modified taking the jump conditions into account. If a uniform grid is used, the interface is generally located between grid points $i$ and $i+1$, which are both irregular. As a result, the corresponding FD schemes are highly complex [13,16], because the distances from the interface to the two irregular grid points enter into the formulation as additional parameters.Taking a Taylor series expansion of $\mathbf {g}_{i+1}$ and $\mathbf {g}_{i-1}$ about grid point $i$, we obtain
where it should be noted that $\mathbf {g}_+^{(n)}\neq \mathbf {g}_-^{(n)}$ for a general $n$ owing to the presence of the interface. Multiplying Eq. (9) by $\mathbf {N}_+^{-1}$ and Eq. (10) by $\mathbf {N}_-^{-1}$ separately, and then eliminating $\mathbf {g}_+^{(1)}$ and $\mathbf {g}_-^{(1)}$ via the explicit expression of Eq. (5), namely $\mathbf {N}_-^{-1}\mathbf {g}^{(1)}_-=\mathbf {N}_+^{-1}\mathbf {g}^{(1)}_+$, we obtain an equation without any terms of the first-order derivative. In this equation, by further eliminating $\mathbf {g}_+^{(2)}$ through the explicit expression of Eq. (6), namely $\mathbf {g}_-^{(2)}+\mathbf {N}_-^2\mathbf {g}_-=\mathbf {g}_+^{(2)}+\mathbf {N}_+^2\mathbf {g}_+$, and finally noting that $\mathbf {g}_-=\mathbf {g}_+=\mathbf {g}_i$ and $\mathbf {g}_-^{(2)}=\mathbf {g}_i^{(2)}$, we arrive at an accurate expression for $\mathbf {g}^{(2)}_i$In the accurate expression of $\mathbf {g}_i^{(2)}$, namely Eq. (11), the summation on the right-hand side contains both odd and even higher-order terms. As will be shown below, to develop a higher-order FD approximation of $\mathbf {g}_i^{(2)}$, these odd and even higher-order terms must be replaced by the FD approximations of $\mathbf {g}_i^{(1)}$ and $\mathbf {g}_i^{(2)}$, respectively, with certain orders. Thus, a parallel set of accurate expressions and higher-order FD approximations of $\mathbf {g}_i^{(1)}$ must also be derived.
Using a similar procedure and eliminating all terms of the second-order derivative from Eqs. (9) and (10), an accurate expression for $\mathbf {g}_i^{(1)}$ can be derived.
2.3 Arbitrary-order FD approximations of first and second-order derivatives
Truncating the infinite summation in Eq. (11) after $s=\gamma$, where $\gamma$ is an integer $\geq 1$, we obtain an approximate expression for $\mathbf {g}_i^{(2)}$ with a $(2\gamma +1)^\text {th}$-order accuracy. Note that the order of accuracy of this expression is not altered after approximating $\mathbf {g}_{\pm }^{(2s+1)}$ using $\delta ^{2s+1}_{2(\gamma -s)+2}\mathbf {g}_{\pm }$ and $\mathbf {g}_{\pm }^{(2s+2)}$ using $\delta ^{2s+2}_{2(\gamma -s)+1}\mathbf {g}_{\pm }$, where $\delta ^n_m\mathbf {g}_\pm$ denotes the $m^\text {th}$-order approximation of $\mathbf {g}_{\pm }^{(n)}$. We then have a $(2\gamma +1)^\text {th}$-order FD approximation of $\mathbf {g}_i^{(2)}$, expressed as
To make Eq. (21) numerically computable, $\delta ^{2s+1}_{2(\gamma -s)+2}\mathbf {g}_{\pm }$ must be related to $\delta ^{1}_{2(\gamma -s)+2}\mathbf {g}_{\pm }$, and $\delta ^{2s+2}_{2(\gamma -s)+1}\mathbf {g}_{\pm }$ to $\delta ^{2}_{2(\gamma -s)+1}\mathbf {g}_{\pm }$. Using the compact FD technique [17], this is accomplished through the application of Eq. (3). Successively taking a second-order differentiation of Eq. (3) $s$ times produces
Using a parallel procedure, we can derive the numerically computable $(2\gamma +2)^\text {th}$-order FD approximation of $\mathbf {g}_i^{(1)}$, that is, $\delta ^1_{2\gamma +2}\mathbf {g}_i$,
The structure of Eqs. (26) and (27) is analyzed as follows: It is always the case that a higher-order approximation of $\mathbf {g}_i^{(2)}$ or $\mathbf {g}_i^{(1)}$ is composed of a linear combination of all the lower-order approximations of $\mathbf {g}_i^{(2)}$ and $\mathbf {g}_i^{(1)}$. For example, $\delta ^2_1\mathbf {g}_i$ and $\delta ^1_2\mathbf {g}_i$ form $\delta ^2_3\mathbf {g}_i$ or $\delta ^1_4\mathbf {g}_i$; and $\delta ^2_1\mathbf {g}_i$, $\delta ^1_2\mathbf {g}_i$, $\delta ^2_3\mathbf {g}_i$, and $\delta ^1_4\mathbf {g}_i$ form $\delta ^2_5\mathbf {g}_i$ or $\delta ^1_6\mathbf {g}_i$. Thus, a higher-order approximation of $\mathbf {g}_i^{(2)}$ or $\mathbf {g}_i^{(1)}$ eventually resolves into a linear combination of $\delta ^2_1\mathbf {g}_i$ and $\delta ^1_2\mathbf {g}_i$: the two basic terms given in Eqs. (13) and (18).
For the sake of program implementation, we rewrite the right-hand sides of Eqs. (26) and (27) as polynomials of $n_{\mathrm {eff}}^2$, namely
2.4 Matrix eigenvalue problem
Definitions (13) and (18) of the two basic FD approximations $\delta ^2_1\mathbf {g}_i=\mathbf {p}_{10}$ and $\delta ^1_2\mathbf {g}_i=\mathbf {N}_-\mathbf {p}_{20}$ can be turned into
Grouping all the FD equations from Eqs. (35) for $i=1,\ldots,N$, where $N$ is the total number of grid points, we obtain a polynomial eigenvalue problem for sparse matrices, which can be reduced to a generalized eigenvalue problem and efficiently solved using a sparse matrix technique [22]. The compact nature of the FD scheme guarantees that the matrices of the final eigenvalue problem are always heptadiagonal for any convergence order. This feature greatly eases the algorithm implementation and reduces the demand for data memory and computation time.
If the grid is uniform, $\delta ^2_{2\gamma +1}\mathbf {g}_i$ upgrades to $\delta ^2_{2\gamma +2}\mathbf {g}_i$, and the uniform grid scheme appears to be more efficient in convergence than the corresponding nonuniform grid scheme. However, the global convergence behavior of the nonuniform grid scheme is not solely determined by the order of the local truncation error. It is also influenced by the allocation of grid points. A layer-wise uniform grid could be used, where only the irregular grid points at the interfaces suffer first-order lower accuracy. In such a case, the nonuniform grid scheme will gain the same global order accuracy as the uniform grid scheme because the codimension-one of the interfaces in the computational domain will exactly complement the local first-order lower accuracy. Additionally, it is possible to further improve the convergence property by utilizing a more adaptive allocation of grid points.
3. Numerical examples
3.1 Achiral multi-quantum-well waveguide
First, we compare the performance between the present formulation and that developed by Sujecki [16] based on a uniform grid for achiral waveguides. The present formulation is relatively general; the program implementation can directly analyze the achiral waveguides as a special case of $\xi =0$ for all layers. However, the TE and TM modes are computed simultaneously. Alternatively, we can derive individual formulations for the TE and TM modes from the general formulation by setting $\xi =0$.
The test model is a multi-quantum-well waveguide consisting of $56$ barriers with a width of $0.012\,\mathrm{\mu}\mathrm {m}$ and a refractive index of $3.2874$, $55$ wells with a width of $0.007\,\mathrm{\mu}\mathrm {m}$ and a refractive index of $3.3704$, and two claddings with a refractive index of $3.2224$, as shown in Fig. 2(a). The operating wavelength is $1.55\,\mathrm{\mu}\mathrm {m}$. For a better comparison, a modal analysis is performed based on the individual formulation of the TM modes; as in [23], the computational window edges are $16\,\mathrm{\mu}\mathrm {m}$ away from the corresponding outermost interfaces. $\mathbf {g}^{(2)}_i$ is approximated using $\delta ^2_1\mathbf {g}_i$, which leads to a standard eigenvalue problem and global second order of accuracy for the layer-wise uniform grid.
Fig. 2. Test model of a multi-quantum-well waveguide without chirality. (a) Illustration of the waveguide structure. (b) Relative error on $n_{\mathrm {eff}}$ versus the number of grid points for the first TM mode.
Figure 2(b) shows the convergence curves of the effective index of the first TM mode, where the relative errors on $n_{\mathrm {eff}}$ are plotted against the number of grid points used to sample the computational domain; for a nonuniform grid scheme, the number of grid points rather than the many-valued step size determines the computational cost. The uniform grid scheme cannot work until the number of grid points is sufficiently large that the step size is comparable to the minimum layer thickness [23]. For the nonuniform grid scheme, the range of applications and the accuracy can be improved through various approaches: (1) A layer-wise uniform (LWU) grid with a fixed allocation ratio of $128:1$ between the cladding layer and fine barrier/well layer; (2) the same allocation ratio as in (1), with the grid points in the cladding layers generated from quadratic stretching (QS) of the coordinates, starting from the outermost interfaces to the corresponding computational window edges; (3) two grid points in each fine barrier/well layer, and the same quadratic stretching (QS*) of the coordinates in the cladding layers as in (2).
The LWU approach has the expected second-order accuracy and nearly the same accuracy level as the uniform grid scheme. The QS approach also has a second-order accuracy, with a steady improvement in accuracy over an order of magnitude compared with the LWU approach. This accuracy enhancement may be due to the better adaptive allocation of the grid points to the evanescent field in the cladding layer. Using the QS* approach, the range of applications extends to a remarkable 230 grid points for this multilayered waveguide with 113 layers; two grid points in the barrier/well layer, and four grid points in the cladding layer. Additionally, the convergence is no longer of the second order, but becomes spectral-like.
From this numerical example, it is clear that compared to the three-point method introduced by Sujecki [16,23], our method has the advantage of extending the usability of the algorithm to the conditions of a lower sampling density for waveguides with fine structures and improving the accuracy for a fixed computational cost. Additionally, the convergent performance of our method could be further improved by learning from data to allocate the grid points.
3.2 Chiral-metallic plasmon waveguide
Subsequently, the validity of the proposed method is verified for the efficient modal analysis of chiral plasmon waveguides. Consider the chiral-metallic single interface plasmon waveguide shown in Fig. 3(a) as the test model, whose modal effective indices could be easily computed with a sufficiently high accuracy by solving the analytic characteristic equation [1]. The relative permittivity and chirality of the chiral layer are $\epsilon _{r,c}=2$ and $\xi _r=0.02$, respectively, the relative permittivity of the metallic layer is $\epsilon _{r,m}=-80+6i$, and both layers are non-magnetic. The ratio of skin depths in the chiral and metallic layers is approximately $40$. The computational window edges are placed where the evanescent fields have decayed to $10^{-16}$ of their initial magnitudes at the interface. An equal number of grid points are placed in each layer according to the QS of the coordinates.
Fig. 3. Test model of a chiral-metallic single interface plasmon waveguide. (a) Illustration of the waveguide structure. (b) Relative error on $n_{\mathrm {eff}}$ versus the number of grid points.
Figure 3(b) shows the convergence curves of the effective index of the surface plasmon mode, approximating $\mathbf {g}_i^{(2)}$ using $\delta ^2_{2\gamma +1}\mathbf {g}_i$ for $\gamma =0,1,\ldots,7$. The orders of accuracy are steadily and nearly precise at $2,4,\ldots,16$, namely the local truncation error of order $2\gamma +1$ corresponds to the global error of order $2\gamma +2$. This demonstrates the arbitrary superconvergent ability of the proposed compact FD method.
Tab. 1 presents the actual orders of accuracy fitted from the convergence curves and data on the relative error and elapsed time for $80$ grid points. The computation was performed on an Intel Core i7 3.40 GHz personal computer. Higher-order schemes reduce the error by orders of magnitude, whereas the elapsed time is approximately proportional to the order of accuracy. When $\gamma =7$, only $80$ grid points can provide a relative error on the effective index below $10^{-15}$, which approaches the limit of double precision. Inevitably, higher-order schemes maintain faster convergence at the expense of a greater accumulation of round-off error. As shown in Fig. 3(b), when $\gamma =7$, the accuracy decreases below that for lower $\gamma$ values with $10$ grid points. A further boost in the convergence rate relies on the improvement of machine precision.
Table 1. Convergence data for Fig. 3(b)
3.3 Chiral PC waveguide
Finally, we consider a test model of a chiral multilayered waveguide. Current research on chirality sensing by means of optical surface waves focuses on metallic plasmon waveguides [1–3]. Nevertheless, the surface waves guided along the boundary of PCs are also attractive for sensing applications owing to their low attenuation at optical frequencies [24,25]. The structure of a waveguide is achiral/chiral-$\text {H}'(\text {LH})^5$-air; for the chiral layer, $\epsilon _{r,c}=1.5$ and $\xi _r=0$ (achiral contrast) or $\xi _r=10^{-4}$, the one-dimensional PC consists of five pairs of low (L) and high (H) index layers, with $\epsilon _{r,L}=2$, $\epsilon _{r,H}=13$ and normalized widths $d_L=1.571$, $d_H=0.453$, and an extra high-index layer with $d_{H'}=0.1d_H$ is inserted between the chiral layer and PC. The waveguide structure is illustrated in Fig. 4(a).
Fig. 4. Test model of an achiral/chiral-PC-air waveguide. (a) Illustration of the waveguide structure. (b) Field distributions of a surface mode in the achiral-PC-air waveguide: TE polarization; the inset demonstrates the continuity of $E_y$ at the first and second interfaces. (c) Field distributions of a surface mode in the chiral-PC-air waveguide: hybrid polarization; the inset demonstrates the continuity of $F_y$ at the third and fourth interfaces. The inner vertical lines in (b) and (c) indicate the interfaces.
The modal fields of the surface modes are shown in Fig. 4(b) and 4(c). The number of grid points in the inner layers are set as the same to effectively resolve the thin H and $\text {H}'$ layers. When $\xi _r=0$, the surface mode is purely TE polarized (see Fig. 4 (b)), and the TM component increases significantly for $\xi _r=10^{-4}$ (see Fig. 4 (c)). The insets in Fig. 4(b) and 4(c) demonstrate that the continuity of $E_y$ and $F_y$ is well enforced.
4. Conclusion
In summary, a three-point FD method with an arbitrary truncation order is proposed for chiral planar waveguides. The combination of a nonuniform grid and compact FD technique not only simplifies the formulation, but also provides a flexible and efficient solution for the modal analysis of chiral waveguides with fine structures for SPP or PC-related applications. This proposal may also be suitable for other one-dimensional wave systems, such as circular fibers and matter waves in quantum-well structures.
Funding
Natural Science Foundation of Ningxia Province (2020AAC03071).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. G. Mi and V. Van, “Characteristics of surface plasmon polaritons at a chiral-metal interface,” Opt. Lett. 39(7), 2028–2031 (2014). [CrossRef]
2. Q. Zhang and J. Li, “Characteristics of surface plasmon polaritons in a dielectrically chiral-metal-chiral waveguiding structure,” Opt. Lett. 41(14), 3241–3244 (2016). [CrossRef]
3. Q. Zhang, J. Li, X. Liu, and D. J. Gelmecha, “Dispersion, propagation, and transverse spin of surface plasmon polaritons in a metal-chiral-metal waveguide,” Appl. Phys. Lett. 110(16), 161114 (2017). [CrossRef]
4. Q. Zhang, J. Li, X. Liu, D. J. Gelmecha, and W. Zhang, “Optical screwdriving induced by the quantum spin hall effect of surface plasmons near an interface between strongly chiral material and air,” Phys. Rev. A 97(1), 013822 (2018). [CrossRef]
5. Q. Zhang, J. Li, and X. Liu, “Optical lateral forces and torques induced by chiral surface-plasmon-polaritons and their potential applications in recognition and separation of chiral enantiomers,” Phys. Chem. Chem. Phys. 21(3), 1308–1314 (2019). [CrossRef]
6. M. Stern, “Semivectorial polarised finite difference method for optical waveguides with arbitrary index profiles,” IEE Proc. J Optoelectron. UK 135(1), 56–63 (1988). [CrossRef]
7. C. Vassallo, “Improvement of finite difference methods for step-index optical waveguides,” IEE Proc. J Optoelectron. UK 139(2), 137–142 (1992). [CrossRef]
8. J. Yamauchi, M. Sekiguchi, O. Uchiyama, J. Shibayama, and H. Nakano, “Modified finite-difference formula for the analysis of semivectorial modes in step-index optical waveguides,” IEEE Photonics Technol. Lett. 9(7), 961–963 (1997). [CrossRef]
9. J. Yamauchi, G. Takahashi, and H. Nakano, “Modified finite-difference formula for semivectorial h-field solutions of optical waveguides,” IEEE Photonics Technol. Lett. 10(8), 1127–1129 (1998). [CrossRef]
10. Chiou Yih-Peng, Chiang Yen-Chung, and Chang Hung-Chun, “Improved three-point formulas considering the interface conditions in the finite-difference analysis of step-index optical devices,” J. Lightwave Technol. 18(2), 243–251 (2000). [CrossRef]
11. Y.-C. Chiang, Y.-P. Chiou, and H.-C. Chang, “Improved full-vectorial finite-difference mode solver for optical waveguides with step-index profiles,” J. Lightwave Technol. 20(8), 1609–1618 (2002). [CrossRef]
12. Y. Lu, L. Yang, W. Huang, and S. Jian, “Improved full-vector finite-difference complex mode solver for optical waveguides of circular symmetry,” J. Lightwave Technol. 26(13), 1868–1876 (2008). [CrossRef]
13. Y.-P. Chiou and C.-H. Du, “Arbitrary-order interface conditions for slab structures and their applications in waveguide analysis,” Opt. Express 18(5), 4088–4102 (2010). [CrossRef]
14. C.-H. Du and Y.-P. Chiou, “Higher-order full-vectorial finite-difference analysis of waveguiding structures with circular symmetry,” IEEE Photonics Technol. Lett. 24(11), 894–896 (2012). [CrossRef]
15. Y.-P. Chiou and C.-H. Du, “Arbitrary-order full-vectorial interface conditions and higher order finite-difference analysis of optical waveguides,” J. Lightwave Technol. 29(22), 3445–3452 (2011). [CrossRef]
16. S. Sujecki, “Arbitrary truncation order three-point finite difference method for optical waveguides with stepwise refractive index discontinuities,” Opt. Lett. 35(24), 4115–4117 (2010). [CrossRef]
17. W. F. Spotz and G. F. Carey, “High-order compact finite difference methods,” in ICOSAHOM.95: Proceedings of the Third International Conference on Spectral and High Order Methods, A. V. Ilin and L. R. Scott, eds. (Houston Journal of Mathematics, 1996).
18. J. A. M. Svedin, “Propagation analysis of chirowaveguides using the finite-element method,” IEEE Trans. Microwave Theory Technol. 38(10), 1488–1496 (1990). [CrossRef]
19. L. Valor and J. Zapata, “A simplified formulation to analyze inhomogeneous waveguides with lossy chiral media using the finite-element method,” IEEE Trans. Microwave Theory Technol. 46(2), 185–187 (1998). [CrossRef]
20. Y. Cao, “Efficient full-vectorial modal analysis based on immersed interface method for dielectric chiral optical fibers,” J. Opt. Soc. Am. A 36(12), 1957–1967 (2019). [CrossRef]
21. Y. Cao, “Efficient mode solver based on immersed interface method for dielectric chiral planar waveguide,” in International Photonics and OptoElectronics Meeting 2019 (OFDA, OEDI, ISST, PE, LST, TSA), (Optical Society of America, 2019), p. JW4A.79.
22. Y. Saad, Numerical Methods for Large Eigenvalue Problems (SIAM, 2011).
23. S. Sujecki, “Accuracy of three-point finite difference approximations for optical waveguides with step-wise refractive index discontinuities,” Opto-Electron. Rev. 19(2), 145–150 (2011). [CrossRef]
24. W. Robertson, “Experimental measurement of the effect of termination on surface electromagnetic waves in one-dimensional photonic bandgap arrays,” J. Lightwave Technol. 17(11), 2013–2017 (1999). [CrossRef]
25. F. Villa, L. E. Regalado, F. Ramos-Mendieta, J. Gaspar-Armenta, and T. Lopez-Ríos, “Photonic crystal sensor based on surface waves for thin-film characterization,” Opt. Lett. 27(8), 646–648 (2002). [CrossRef]
