## Abstract

A novel bidirectional operator marching method based on the Dirichlet-to-Neumann (DtN) mapping for three-dimensional optical waveguide structures is developed and implemented using iterative methods. The backward propagation wave is integrated into the classical operator marching approach which represents the forward propagating wave. The bidirectional range marching formulas are exact for each range-independent piece and a large range step is possible in both directions. The validity and effectiveness of our proposed method are verified by analyzing uniform waveguides and longitudinal waveguides with varying refractive indices.

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

## 1. Introduction

An efficient computational method for simulation of optical waveguide structures is of fundamental importance in the design of new optical devices. Since a certain number of grid points are needed for each wavelength, a direct numerical approach by the finite difference (FD) method or the finite element (FE) method that discretize the whole wave-guiding structure is prohibitively expensive [1–3]. For the slowly varying waveguides or piecewise invariant structures, some efficient numerical methods, such as the beam propagation met hod (BPM) [4–6] and the operator marching method (OMM) [7–9], have been extremely successful in computation of the forward propagating optical fields.

However, many practical optical devices, such as laser facets, grating structures, multiple dielectric layers or junctions of different waveguides, contain refractive index discontinuities along the direction of propagation of the optical field. Therefore, a bidirectional technique which considers coupling of the forward and backward waves is needed for accurate modeling of reflection or transmission [10–13]. The OMM, based on the Dirichlet-to-Neumann (DtN) mapping with the exact reformulations for pairs of operators, is used to solve the propagation problem in three dimensional optical waveguides that is formulated in a domain with just one direction (i.e. $z$) having a particularly large length, due to its excellent simulation speed with sufficient accuracy. In [14], a three-dimensional OMM of the scalar Helmholtz equation in the three-dimensional rectangular box with the radiation boundary was developed. In [15], a three-dimensional adjoint-based OMM (AOMM) of the complex Helmholtz-PML equation for long-range inhomogeneous optical waveguide devices was proposed. Hence, there is an understandable desire to incorporate the bidirectional technique and OMM into developing the so-called bidirectional OMM (Bi-OMM) in three dimensional optical waveguides since it can simultaneously consider the forward and backward waves.

In this paper, we propose a new Bi-OMM based on the bidirectional DtN mapping for three-dimensional optical waveguide structures. The bidirectional reflection operator is firstly determined, subsequently the exact one-way re-formulation is constructed through the bidirectional DtN map [7,16,17], and finally the local base transformation is implemented by the series-expansion technique [18,19]. To validate the present approach, a typical buried waveguide and the surface plasmon stripe waveguide with complex varying refractive indices are considered as the numerical examples [5,20]. Results show that the proposed Bi-OMM can deal with evanescent waves properly with accuracy comparable to that of the well-established bidirectional beam propagation method (Bi-BPM) [4,10]. And in the inhomogeneous waveguides, our proposed Bi-OMM also shows more effective with comparison of the AOMM.

The rest of this paper is organized as follows. Section 2 presents the basic mathematical model in unbounded waveguides. In Section 3 and 4 explain the formulation of the bidirectional operator marching scheme in detail, which is the main goal of this paper. To show the effectiveness and usefulness of the proposed approach, some numerical results are shown in Section 5. The conclusion is drawn in Section 6.

## 2. Basic problem

Considering the optical waves propagating in a three-dimensional waveguide structure, the propagation direction is supposed to be $z$-direction and the transverse cross-section is set to be $xy$-plane. In the Bi-OMM, the longitudinally varying integrated optical waveguide is first discretized and replaced by a suitable number of longitudinally invariant inhomogeneous segments. Then, by utilizing the governing wave equations in frequency-domain for a monochromatic wave, the electric field satisfies the following three-dimensional complex Helmholtz equation

The transverse operator on the $xy$ plane is defined as

At $z =0$, the incident wave is given as

And an outgoing condition can be written down at $z=L$,

Therefore, the propagation problem is formulated in a domain with just the main propagation direction (i.e., $z$ direction), having a particularly large length, and a relatively small transverse region.

## 3. Bidirectional operator marching

The total wave field can be split into right and left propagating components, denoted by $E_R$ and $E_L$, each of which satisfies a one-way re-formulation based on the DtN map. We write

and the right and left propagating fields in each $z$-invariant segment are, respectively, identified with the two formal solutions to Eqs. (6) and (7). The electric field at any point along the propagation direction $z$ can be represented by the right ("$E_R$") and left ("$E_L$") waves. And it follows that the wave field splitting in a transversely inhomogeneous environment, is accomplished by the DtN map, where $Q_1$ and $Q_2$ are the two operator solutions of the quadratic operator equationSubstituting Eq. (8) into the Eq. (1), we have

Since the above formulas Eqs. (10) and (11) are true for any solution of Eq. (1) (satisfying the boundary conditions at (3), and the outgoing condition at $z = L$), we obtain

andThe above Eqs. (12) and (13) should be solved for decreasing $z$ with an exact radiation condition at $+\infty$. Then, the fundamental solution operators $G_1$ and $G_2$ have been introduced for stability reasons,

andThe initial condition at $z=L$ is $G_1(L) = I$ and $G_2(L) = I$, where $I$ is the identity operator. When $G_1(0)$ and $G_2(0)$ are calculated, the solution at $z =L$ can be generated by the multiplications with the starting field $E_R(x, y, 0)$ and $E_L(x, y, 0)$, namely

The formulation above allows us to evaluate $E(x, y, L)$, given $E(x, y, 0)$. However, typically we are given $E_R(x, y, 0)$, and should determine $E_L(x, y, 0)$ such that the boundary condition that no backward field will be present at the output, i.e., $E_L(x, y, L) = 0$, is satisfied.

Now, we consider the multiple-layered structure shown in Fig. 1, which consists of $N$ interfaces, separating $N-1$ transversely inhomogeneous layers. Let the discrete points $z_i$ ($1\leq i < N$) satisfy

Each interval $(z_i, z_{i+1})$ corresponds to a range-independent segment and $h = L/N$.

The approximated Helmholtz equation by piecewise range-independent segments using the medium values at $(z_i, z_{i+1})$ is

where $z_{i+1/2} = z_i + h/2 = (z_i + z_{i+1})/2$. Exact solutions in each segment can be written down and continuity conditions at $z_{i+1/2}$ are used to match the solutions. We take $i=0$, for example, to describe the bidirectional marching form $z_0$ to $z_1$.On the interval $(z_0, z_1)$, the Eqs. (12) and (13) are approximated by

and where $\mathscr {B}_{1/2}=\mathscr {B}({z_{1/2}})$. Given the initial conditions $Q_1^{(1)}\approx Q_1(z_1)$ and $Q_2^{(1)}\approx Q_2(z_1)$, the exact solutions of $Q_1^{(0)}$ and $Q_2^{(0)}$ are used to obtain $Q_1^{(0)}\approx Q_1(z_0)$ and $Q_2^{(0)}\approx Q_2(z_0)$.Then the total wave field can be decomposed into the right and left going waves, satisfying

andThe solution to this general problem can be written symbolically, in terms of the introduced reflection operators $\mathscr {R}$ for right and left propagating incident fields. At the point $z_i$, we define the reflection operator $\mathscr {R}_i$ by

And it is easy to obtain from $E_z(x, y, z_1) = Q_1^{(1)} E_R + Q_2^{(1)} E_L$ and the above relationship that

Similarly, at the point $z_0$, we can obtain

Compare the definition of $\mathscr {R}_1$ and $\mathscr {R}_0$ and notice that

The formulas (24), (25) and (27) reveal the exact relationship between $Q_1^{(0)}$, $Q_2^{(0)}$ and $Q_1^{(1)}$, $Q_2^{(1)}$ for Eqs. (19) and (20). In the DtN formulation above, the continuity condition is implicit in the requirement that $Q_1^{(1)}$ and $Q_2^{(1)}$ obtained from the previous calculation in $(z_1, z_2)$ be the same $Q_1^{(1)}$ and $Q_2^{(1)}$ used to calculate $Q_1^{(0)}$ and $Q_2^{(0)}$ in $(z_0, z_1)$.

Therefore, on the interval $(z_0, z_1)$, we have

At this point, we know how to construct the bidirectional DtN operators $Q_1$ and $Q_2$, and the fundamental solution operators $G_1$ and $G_2$. Meanwhile, solving such operators from $z=L$ to $z =0$ with the formulas (24), (25) and (27) are used in the bidirectional marching step from $z_{i+1}$ to $z_i$, where $1\leq i\leq N$.

## 4. Local base transformation and iteration schemes

In this section, a new treatment for local base transformation is introduced [21], then the above bidirectional operator marching formulas are considered into numerical implementation through the matrices.

The eigenvalue problem of the transverse operator $\mathscr {B}$ is important, and it is defined as

corresponding to the boundary condition (3). The square root operator,We will also take the segment $(z_0, z_1)$ for example. For integers $n_x$ and $n_y$, the transverse operator $\mathscr {B}$ are expanded by bivariate Chebyshev coefficients such that

The derivatives of two transverse directions $x$ and $y$ can be expressed in terms of the grid values $\mathscr {B}(x(\xi _i),y(\xi _j),z_{1/2})$. That is, the $p$th derivatives

In consideration of the boundary conditions (3), we can make use of the $(n_x-1) \times (n_x-1)$ matrix ${\mathscr {D}}^2_{n_x-1}$ obtained by stripping $\tilde {\mathscr {D}}_{n_x}^2$ of its first and last rows and columns. In MATLAB notation:

where $\tilde {\mathscr {D}}$ is the first derivative Chebyshev differentiation matrix of size $(n_x +1) \times (n_x +1)$. Similarly Thus we can get the matrix form of the transverse operator $\mathscr {B}(z_{1/2})$ is where $n_x = n_y$, $\otimes$ is the kronecker product, $K_{1/2}=\kappa _0^2n^2(x(\xi _k),y(\xi _t),z_{1/2})I \otimes I$ and $I$ is the identity matrix.Suppose $\nu (x, y)_m$ and $\hat {\nu }(x,y)_m$ are the normalized eigenvectors of $\mathscr {M}$ and $\mathscr {M}^{\textrm{H}}$, respectively, and assume their corresponding eigenvalues are $\lambda _m$ and $\bar {\lambda }_m$, where $\mathscr {M}^{\textrm{H}}$ is the conjugate transverse of $\mathscr {M}$. Then

where $\left \langle {*,*}\right \rangle$ is the Hilbert space $L^2(\Omega )$ inner product $\left \langle {\phi ,\psi }\right \rangle = \int _\Omega \phi \bar {\psi }d\Omega$, and the overline indicates complex conjugation.Otherwise, the eigenvector decomposition of $\mathscr {M}_{1/2}$ and $\mathscr {M}^{\textrm{H}}_{1/2}$ are

Thus the bidirectional operator marching algorithm (Bi-OMM) can be outlined as follows:

(1) $E_R^{(0)}$ is given, $E_L^{(0)}$ is guessed (e.g., set to 0) and choose a step size $h >0$;

(2) for $i = 0, 1, 2, \ldots$ do:

a) $\mathscr {O}_1 = \mathscr {N}\hat {\mathscr {O}}_1\mathscr {N}^{-1}$,

$\mathscr {O}_2 = \mathscr {N}\hat {\mathscr {O}}_2\mathscr {N}^{-1},$

$\mathscr {P}_1 := ( \textrm{i} \sqrt { \Lambda _{1/2, k} } + \mathscr {O}_2 )^{-1}(\textrm{i}\sqrt {\Lambda _{1/2, k}}-\mathscr {O}_1),$

$\mathscr {P}_0 = e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}}\mathscr {P}_1e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}},$

$\mathscr {S}_0 := \mathscr {N}\mathscr {S}_1\mathscr {N}^{-1}(I + \mathscr {P}_1)e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}}(I + \mathscr {P}_0)^{-1},$

$\left ( {\begin {array}{c} {{E_R}^{(i)}}\\ {{E_L}^{(i)}} \end {array}} \right ) = \hat {\nu }_{1/2}\mathscr {S}_0\hat {\nu }^{\textrm{H}}_{1/2}\left ( {\begin {array}{c} {{E_R}^{(0)}}\\ {{E_L}^{(i)}} \end {array}} \right )$;

b) if $\left | {{E_L}^{(k)}} \right | < 1e - 6$, stop; otherwise,

c) $\hat {\mathscr {O}}_1 = \mathscr {N}^{-1}\mathscr {O}_1\mathscr {N}$,

$\hat {\mathscr {O}}_2 = \mathscr {N}^{-1}\mathscr {O}_2\mathscr {N}$

$\mathscr {P}_0 := ( \textrm{i} \sqrt { \Lambda _{1/2, k} } + \hat {\mathscr {O}}_2 )^{-1}(\textrm{i}\sqrt {\Lambda _{1/2, k}}-\hat {\mathscr {O}}_1),$

$\mathscr {P}_1 = e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}}\mathscr {P}_0e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}},$

$\mathscr {S}_1 := \mathscr {N}^{-1}\mathscr {S}_0\mathscr {N}(I + \mathscr {P}_0)e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}}(I + \mathscr {P}_1)^{-1},$

$\left ( {\begin {array}{c} {*}\\ {{E_L}^{(i+1)}} \end {array}} \right ) = \hat {\nu }^{\textrm{H}}_{3/2}\mathscr {S}_1\hat {\nu }_{3/2}\left ( {\begin {array}{c} {{E_R}^{(i)}}\\ {\alpha \cdot {E_L}^{(i)}} \end {array}} \right )$, where $0 < \alpha <1$ and * indicates an unused component.

In brief, at each iteration, the guess for $E_L^{(i)}$ is propagated along with $E_R^{(0)}$ forward through the system. The fields $E_R^{(i)}$ and $E_L^{(i)}$ are propagated backward through the structure to obtain a new guess for $E_R^{(i+1)}$. Then, $E_L^{(i)}$ should be zero when we reach a solution. The damping factor, $\alpha$, is selected to optimize convergence. Compared with the Bi-BPM, our method is developed from the exact formulation of pairs of operators, and no paraxial or slowly varying envelope approximation has been used as in the BPM.

## 5. Numerical results

Two numerical examples will be presented here to illustrate the validity of the proposed Bi-OMM. The relative error (RE) $\varepsilon$ using the discrete $L^2$ norm is defined:

#### 5.1 Buried waveguide

To validate the proposed Bi-OMM, we first apply it to analyze a typical Buried waveguide, as shown in Fig. 2 [4,15]. The core of the buried waveguide is assumed to have the reflective index $n_2 = 3.54$. The reflective index of the cladding material is $n_1 = 3.17$. The core thickness is fixed to $d = 0.4$ $\mu$m, while the core width is set to be $\omega = 0.4$ $\mu$m. The operating wavelength is assumed to be $\lambda = 1.30$ $\mu$m. The propagation distance is $L = 10$ $\mu$m. The width of the transverse computational window is $D_1 = D_2 = 1$ $\mu$m. Meanwhile, the transverse operator is approximated with a grid of $20\times 20$ and $k = 8$ in the local subtraction of the eigenmodes.

The starting field at $z=0$ is $E(x, y, 0) =\sin (2\pi x)\sin (2\pi y)$, which is the third propagating mode, and this case has an analytical solution

The exact fields in Eq. (44) are used for reference.In Fig. 3, the numerical fields $EL_R\approx E_R(x, y, L)$ with the step length $h = 1$ $\mu$m are plotted. Meanwhile, in Fig. 4, the RE is shown for the numerical fields $EL_R$ with $h = 1$ $\mu$m compared to the exact fields for a buried waveguide in the $XY$ plane. In contrast, the RE for the numerical fields and the exact solution is small enough, and the accuracy of the proposed Bi-OMM is pretty good in the buried waveguide. Thus good solutions are already obtained with $h =1$ $\mu$m. In Table 1, we list the relative errors ( compared with the "exact" solution) of the numerical solutions of $E(x, y, L)$ with different choices of the range step size.

For this case, we also compare the efficiency of the Bi-OMM, FD-Bi-BPM, and AOMM for the buried waveguide by finding the number of nodes, CPU time and the required computing memory for two predetermined $\varepsilon$ error goals of $0.01$ and $0.001$. The results in Table 2 show the obvious superiority of the proposed Bi-OMM with the much lower computational cost over the FD-Bi-BPM and the results are especially impressive for the higher accuracy goal.For the buried waveguide, our proposed Bi-OMM can achieve similar accuracy as the adjoint-based OMM (AOMM), meanwhile, avoiding the construction of the adjoint operator for the transverse operator that usually is complicated. When the transverse operator is approximated with a grid of $40 \times 40$, and similarly $k=8$ for the eigenmode local subtraction, the effectiveness of the proposed Bi-OMM in the buried waveguide has no significant change. It is mainly because the differentiation matrices Eq. (37) are dense, and the series-expansion method is a powerful mode solver, with its low memory storage requirements for a small matrix, could obtain highly accurate results [18,22].

#### 5.2 Surface plasmon stripe waveguide

To further confirm the applicability of the proposed Bi-OMM to plasmonic devices, we consider a surface plasmon stripe waveguide structure (as shown in Fig. 5) and discuss the numerical accuracy of Bi-OMM in lossy and leaky inhomogeneous waveguide analysis. The refractive indices of the cladding (air), core and substrate, are $n_0 = 1.0$, $n_1 = (-26.1437 - \textrm{i}1.8479) [1 + 0.05e^{-(z/L-0.5)^2}\sin ^2(\pi x)\sin ^2(\pi y)]$ and $n_2 = 1.46^2$. The core height is fixed to $d = 55$ nm, while the core width is set to be $\omega = 3.3$ $\mu$m. The propagation distance is set to be $L = 2$ mm, and $D_1 = D_2 = 10$ $\mu$m. The operating wavelength is set to be $0.8$ $\mu$m. Meanwhile, the transverse operator is discretized by a grid of $24\times 24$ and $k = 10$ in the local subtraction of the eigenmodes. The starting field at $z = 0$ is set to be $E(x, y, 0) = \sin (2\pi x)\sin (2\pi y)e^{-x-y}.$

The numerical fields calculated by our Bi-OMM with $h = 60$ $\mu$m are shown in Fig. 6. The RE on the numerical fields $EL_R$ with $h = 60$ $\mu$m compared to a much more accurate solution obtained with $h/32$ for a surface plasmon stripe waveguide in the $XY$ plane are thoroughly presented in Fig. 7. We can see from the RE of Fig. 7 that the results are still satisfactory and a reasonably good solution is already obtained with $h$. As above, we list the relative errors of the numerical solutions of $EL_R$ with different choices of the range step size in Table 3. And the numerical fields with a step length $h/128$ are used for reference, because analytic solutions are not available here.

A comparison of the computational cost of the proposed Bi-OMM, FD-Bi-BPM and AOMM for the surface plasmon stripe waveguide with two predetermined $\varepsilon$ error goals of $0.01$ and $0.001$ is also shown in Table 4. The proposed Bi-OMM has much lower computational cost compared to the FD-Bi-BPM, the required computer memory and CPU time needed are also much smaller. For the varying refractive index profiles, it is even harder than the constant refractive index profile case of the Bi-BPM to do rational approximations in each segment. For the surface plasmon stripe waveguide with the varying refractive index profiles, our proposed Bi-OMM is more efficient than the AOMM, and what is more, the backward waves can also be considered by the Bi-OMM. When the transverse operator is approximated with a grid of $48 \times 48$, and similarly $k = 10$ for the eigenmode local subtraction, the effectiveness of the proposed Bi-OMM in a surface plasmon stripe waveguide is still satisfactory as above.

## 6. Conclusion

In this paper, we have developed an efficient Bi-OMM based on the bidirectional DtN mapping and the large range step method for three-dimensional optical waveguide structures. Using a traditional buried waveguide example, the technique proves to be very accurate compared to the analytical result. Numerical results also show that the present Bi-OMM can deal with inhomogeneous waveguide structures properly with accuracy comparable to the well-established FD-Bi-BPM and less demanding in terms of memory and CPU-time resources. Currently, work is being done to extend the proposed method to the full-vectorial formulation.

## Funding

Zhejiang Provincial Natural Science Foundation of China (LQ19A010012).

## 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. **A. W. Snyder and J. D. Love, * Optical Waveguide Theory* (Chapman and Hall, 1983).

**2. **R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. **6**(1), 150–162 (2000). [CrossRef]

**3. **M. Levy, * Parabolic Equation Methods for Electromagnetic Wave Propagation* (IET, 2009).

**4. **S. Obayya, * Computational Photonics* (Wiley, 2011).

**5. **S. Kawai, A. Iguchi, and Y. Tsuji, “Study on high precision and stable finite element beam propagation method based on incomplete third order hybrid edge/nodal element,” J. Lightwave Technol. **36**(11), 2278–2285 (2018). [CrossRef]

**6. **X. Ma, S. Zhu, L. Li, J. Li, X. Shao, and A. Galvanauskas, “Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence,” Opt. Express **26**(4), 3995–4003 (2018). [CrossRef]

**7. **Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. **152**(1), 231–250 (1999). [CrossRef]

**8. **Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. **1**(6), 1056–1075 (2006).

**9. **S. Zhao and J. Zhu, “Efficient marching treatment for optical propagation in the waveguide with a discontinuous interface,” J. Opt. Soc. Am. B **35**(9), 2326–2333 (2018). [CrossRef]

**10. **H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. **11**(7), 830–832 (1999). [CrossRef]

**11. **Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. **14**(11), 1533–1535 (2002). [CrossRef]

**12. **S. Wu and J. Xiao, “An efficient semivectorial bidirectional beam propagation method for 3-D optical waveguide structures,” J. Lightwave Technol. **34**(4), 1313–1321 (2016). [CrossRef]

**13. **A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A **100**(5), 053818 (2019). [CrossRef]

**14. **J. Zhu and H. Yang, “Efficient marching solver for a class of three-dimensional optical waveguides,” Opt. Express **26**(12), 15110–15123 (2018). [CrossRef]

**15. **H. Yang and Z. Tang, “3D adjoint-based marching scheme for optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. B **37**(5), 1298–1306 (2020). [CrossRef]

**16. **L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. **28**(5), 865–876 (1993). [CrossRef]

**17. **L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. **18**(18), 1967–1969 (2006). [CrossRef]

**18. **J. Zhu and H. Yang, “Adaptive mode solver for varying refractive index profile’s waveguides with modified spectral element method,” J. Opt. Soc. Am. B **34**(3), 538–545 (2017). [CrossRef]

**19. **C. C. Huang, C. C. Huang, and J. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profiles using spectral collocation method with domain decomposition,” J. Lightwave Technol. **21**(10), 2284–2296 (2003). [CrossRef]

**20. **H. H. Liu and H. C. Chang, “High-resolution analysis of leaky modes in surface plasmon stripe waveguide,” J. Lightwave Technol. **34**(11), 2752–2757 (2016). [CrossRef]

**21. **J. Zhu and G. Wang, “High-precision computation of optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. A **32**(9), 1653–1660 (2015). [CrossRef]

**22. **L. N. Trefethen, * Approximation Theory and Approximation Practice* (SIAM, 2013).

**23. **A. Abdrabou, A. M. Heikal, and S. S. A. Obayya, “Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis,” Opt. Express **24**(10), 10495–10511 (2016). [CrossRef]