Abstract
In this paper, the smooth approximation of light waves is studied for an open optical waveguide with a distinct refractive-index profile, which involves high-precision computation of the eigenmodes and corresponding eigenfunctions. During analysis, the refractive-index function is first approximated by a quadratic spline interpolation function. Since the quadratic spline function is a polynomial of degree two in every sub-interval (sub-layer), it is equivalent to a piecewise polynomial of degree two, based on which, the corresponding Sturm-Liouville eigenvalue problem of the Helmholtz operator in sub-layer can be solved analytically by the Kummer functions. Finally, the approximate dispersion equation is established to the TE case. Obviously, the approximate dispersion equations converge fast to the exact ones, as the maximum value of the sub-interval sizes tends to zero. Furthermore, eigenmodes may be obtained by Müller’s method with suitable initial values. To refine the accuracy, the equidistant partition and the non-equidistant partition are applied to divide the interval. Numerical simulations show that the eigenfunctions of the spline interpolation are much smoother than the ones with piecewise interpolation. In addition, the non-equidistant partition can help improve the accuracy and the order of convergence of general solutions reaches the third.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Corrections
11 October 2021: A typographical correction was made to the author affiliations.
1. Introduction
The propagation characteristics (called as propagation constants) of optical waveguides play an important role in the optimal designs of optical devices, including fiber [1,2] and the computation of wave propagation with the eigenmode expansion method. Furthermore, the wave field can be decomposed into a linear combination of a few discrete guided modes and an infinite integral of a continuum of radiation modes [3]. The infinite integral in the decomposition of a field can be approximately expressed as an infinite sum of leaky modes through discretization [4,5], which requires to develop some efficient methods to obtain high-precision propagation constants [6,7]. In general, numerical methods, such as the finite difference method [8,9], the finite element method [10], the multi-domain pseudo-spectral method [11,12], and the B-spline modal method [13], can help simplify the problem into a matrix eigenvalue problem. However, for the reason that these methods all arise complex matrices, it is too complicated to calculate the related eigenvalues, let alone satisfying high precision requirements. Now, two problems are encountered. For one thing, the propagation constants of the guided modes are easier to get by common methods compared with the unreal or unbounded ones, since they are real and bounded [14]. For another, the propagation constants of the leaky modes are complex numbers [15], implying the computation difficulty, especially those with large modules. To solve above tasks, an alternative approach is proposed, that is to transform the eigen problem into a nonlinear dispersion equation, whose roots correspond to the eigenvalues or the propagation constants. Some methods can be adopted to deal with the piecewise constant slab waveguides [14–16] and rib waveguides [17]. For the varying refractive-index waveguides, part of approximate dispersion equations are constructed by the Wentzel-Kramers-Brillouin (WKB) method [18] and the differential transfer matrix method (DTMM) [19–22]. However, these approximations cannot guarantee the modes’ accuracy. Recently, some efficient methods are provided to compute the propagation constants, but their corresponding eigenfunctions are not smooth [23], but continuous, which changes the physical meaning of the line.
In this article, we consider the eigenmode problems of varying refractive-index waveguides in TE case. First, we approximate the refractive-index function by a quadratic spline interpolation function. Second, since the quadratic spline function is a polynomial of degree two in every sub-interval (sub-layer), the corresponding Sturm-Liouville eigenvalue problem of the Helmholtz operator in sub-layer can be solved analytically by the Kummer functions. Finally, the approximate dispersion equation is established to the TE case. The obtained dispersion equations can be solved by the secant method or Muller’s method. Numerical simulations show that our method reaches the third order accuracy. The error estimate originally proved for the real eigenvalues in [24], is still valid for the complex eigenvalues. The rest of this paper is organized as follows. In Section 2, the mathematical model and treatment are introduced, and the dispersion equation is established, where asymptotic solutions of the propagation constants (corresponding to leaky modes) are provided for TE case. In Section 3, numerical results and corresponding analysis are presented. Finally, the conclusions are given in Section 4.
2. Mathematical modelling
In this section, we consider a general two-dimensional waveguide with a refractive-index profile ${n_0}(x)$ shown in Fig. 1. For simplicity, we assume that the medium is homogeneous for $x < 0$ and $x > d$ with d being a constant. That is
In the TE case, Eq. (1) in $[0,d]$ is given by
First, as $x \in [{x_j},{x_{j + 1}}] \subset [0,d]$, the function $\mathrm{\kappa} _0^2n_0^2(x)$ is interpolated by a polynomial of degree two, that is $\mathrm{\kappa} _0^2n_0^2(x)\textrm{ - }{{\beta} ^2} \approx {a_j}{x^2} + {b_j}x + {c_j}$. Denote ${y_j}(x)$ as the approximation of the $\mathrm{\phi} (x)$. Thus, the form of approximated equation is as following form:
2.1 Interval partitions
We adopt following two partitions as follows.
2.1.1 Equidistant partition
The interval $[0,d]$ is divided into k subintervals with width $h = d / k$, then the subintervals are denoted by ${I_j} = [{x_{j - 1}},{x_j}]$ related to ${x_j} = jh,(j = 1,2, \cdot{\cdot} \cdot ,k)$.
2.1.2 Adaptive partition
The interval $[0,d]$ is arbitrarily divided into k subintervals $[{x_1},{x_2}],[{x_2},{x_3}], \cdot{\cdot} \cdot ,[{x_k},{x_{k + 1}}]$. Assume $n_0^2(x)$ is a bounded variation function on the interval $[0,d]$, we have
where $V_0^d({n_0^2(x)} )$ and $V_{{x_i}}^{{x_{i + 1}}}({n_0^2(x)} )$ are represented as the total variations of $n_0^2(x)$ on the intervals $[0,d]$ and $[{x_i},{x_{i + 1}}]$, respectively, $(i = 1,2, \cdot{\cdot} \cdot ,k).$2.2 Interpolation methods
We have two forms as follows.
2.2.1 By the piecewise quadratic interpolation
The coefficients of Eq. (5) are given by
Here ${t_0} = (j - 1)h$, ${t_1} = (j - 1/2)h$,and ${t_2} = jh,(j = 1,2, \cdot{\cdot} \cdot ,k)$.
2.2.2 By the quadratic spline interpolation
The coefficients of Eq. (5) are given by
Furthermore, ${M_j}$ and ${M_{j + 1}}$ satisfies
To fix ${M_j},(j = 1,2, \cdot{\cdot} \cdot ,k)$, we choose the below boundary condition: $S^{\prime\prime}(0) = {M_2} - {M_1} = 0$.
Moreover, Eq. (5) can be changed to the confluent hypergeometric equation, and a pair of linearly independent solutions $\{ {\alpha _j}(x),{{\beta} _j}(x)\}$ as $x \in [{x_j},{x_{j + 1}}]$ are given by the Kummer functions. Then the approximated solution of $\mathrm{\phi}$ as $x \in [{x_j},{x_{j + 1}}]$ can be expressed in the following form:
To obtain the constants of ${A_j}$ and ${B_j}$, $(j = 1,2, \cdots ,k)$, we require the solutions ${y_j}(x),(j = 1,2, \cdots ,k)$ of field $\mathrm{\phi} (x)$ and their first order derivatives are continuous. Thus, the interface conditions at ${x_j},(j = 1,2, \cdots ,k)$ are
And the outgoing boundary conditions at $x = 0$ and $x = d$ become
According to these conditions, a linear system of ${A_j}$ and ${B_j}(j = 1,2, \cdot{\cdot} \cdot ,k)$ can be obtained.
In order to simplify these equations, let ${T_j} = {B_j}/{A_j},(j = 1,2, \cdot{\cdot} \cdot ,k)$, the dispersion relations of ${\beta}$ are
2.3 Asymptotic solutions for the TE case
By the WKB method, we can obtain the following formulas [18].
For the sake of unification, let m be the rearrangement index (see [14]). In order to find the roots of the nonlinear equation $F({\beta} ) = 0$, we adopt Müller’s method [25] for implementation. Müller’s method requires three initial values, which are taken by the asymptotic solutions. The error analysis of the approximation is listed in the APPENDIX 1.
3. Numerical example and comparison
At last, the theory of this paper has been implemented and tested on a number of examples. Let ${n_0}(x) = 3.5 \cdot [1 - 0.08{(x - 1)^2}]$, ${n_1} = 2.1$, ${n_2} = 3.17$, $d = 2 \mu \textrm{m}$, $\lambda = 1.55 \mu \textrm{m}$, ${{\kappa} _0} = \frac{{2\pi }}{\lambda }$. The iteration stopping criteria is chosen as $|F |< {10^{ - 10}}$.
3.1 Comparison of the fitting effect by two kinds of spline interpolation methods
The fitting diagrams of the spline interpolation functions are shown in Figs. 2 and 3, respectively.
It is demonstrated that the approximation effect of the non-equidistant and adaptive spline interpolation (see Fig. 3) is better than that of the equidistant spline interpolation (see Fig. 2).
We have found that two fitting curves obtained by the quadratic spline interpolation (equidistant partition) with $k = 16$ and the piecewise quadratic interpolation with $k = 2$ are almost consistent with the original curve $y = {\kappa} _0^2n_0^2(x)$, which matches the reality. The result is not a coincidence, for the original curve is a quadratic polynomial. Since the additional boundary condition at $x = 0$ of the quadratic spline interpolation is artificially assigned, rather than the boundary condition of the original curve, its spline interpolation function is not equal to the original one. However, the spline interpolation function can perfectly approximate the original function as long as the value of k is large enough.
3.2 Comparison of eigenvalues obtained by different spline interpolation methods
The diagrams of approximate eigenvalues as $k = 8$ and $k = 16$ are shown in Figs. 4 and 5, respectively. It is clear that the solutions ${\beta}$ of the spline interpolation are generally similar to the ones of piecewise interpolation except the accuracy as $k = 8$. To address this problem, we adopt non-equidistant partition to improve the accuracy, which can help achieve good results with $k = 8$ and $k = 16$ shown in Fig. 6.
3.3 Numerical verification of the convergence rate
For the error analysis, we choose the absolute value of the relative error, which is denoted by ${E_r}$, and plot ${\log _2}({E_r})$ against ${\log _2}(h)$ in Figs. 7, 8 and 9, respectively. In these figures (Figs. 7–9), as k increases, namely ${\log _2}(2/k)$ decreases, the slopes between the points are not less than 3. That is, our method meet the accuracy of order 3, and the accuracy can be improved by non-equidistant partition. That is to say, the order of convergence of general solutions (eigenmodes) reach the third.
3.4 Comparison of eigenfunctions obtained by two kinds of interpolation methods
The diagram of eigenfunctions are shown in Figs. 10, 11, 12 and 13, respectively. The simulation results show that, the eigenfunction of the spline interpolation is smoother than the one of the piecewise interpolation.
Obviously, the smoothness of the eigenfunction obtained by the spline interpolation can be well preserved regardless of the value of k. However, the piecewise interpolation cannot work like the spline interpolation, since if we apply the piecewise interpolation, there may be some jumping points. For example, in Figs. 10–13, the jumping points appears near the beginning or the end of the eigenfunction. The reason is that piecewise interpolation only keeps the function continuous, but it cannot guarantee the continuity of the first derivative. Therefore, the deviation of numerical calculation appears, which in reality, is not consistent with the wave propagation.
4. Conclusions
For a variety of refractive-index waveguide, the function ${\kappa} _0^2n_0^2(x)$ is approximated by a quadratic spline interpolation function. Then, it can establish to the approximate dispersion equation, only for TE case. Since the approximate equation is equivalent to the confluent hypergeometric equation, its’ solutions can be expressed analytically by the Kummer functions in each layer (sub-interval). In the numerical example, we find the roots of the dispersion equation (i.e. the propagation constants corresponding to the leaky modes) by Müller’s method, where three different asymptotic solutions play the roles of initial values. Obviously, the approximate dispersion equation converges fast to the exact one for the continuous refractive-index function, as the maximum value of the subinterval sizes tends to zero. In addition, to refine the accuracy, we applied the equidistant partition and the non-equidistant partition to divide the interval. We apply the piecewise quadratic interpolation and the quadratic spline interpolation to solve the approximate function, respectively. Numerical simulations demonstrate that the eigenfunction of the spline interpolation is smoother, compared with the one with the piecewise interpolation. Moreover, the non-equidistant partition can help promote the accuracy, and the order of convergence of general solutions meet the third.
Appendix 1
Error analysis: the error of the coefficient in Eq. (5) will inevitably lead to the error of the exact solution of the eigenvalue, and the error of eigenvalue can be obtained as follow:
Further,
Case 1: For the piecewise interpolation, we can obtain the remainder formula of the piecewise quadratic interpolation:
where $|{\omega _3}(x)|= |(x - {x_j})(x - \frac{{{x_j} + {x_{j + 1}}}}{2})(x - {x_{j + 1}})|$, which can reach the maximum value when $x = {{({x_j} + {x_{j + 1}})} / 2} + {3^{ - 1/2}}{{({x_{j + 1}} - {x_j})} / 2}$.Let $h = \max \{ {h_1},{h_2}, \cdot{\cdot} \cdot ,{h_k}\}$, then the error estimation of the piecewise quadratic interpolation is shown as follow:
Case 2: For the quadratic spline interpolation, we can obtain the remainder formula as follow [26]:
whereThe integral result of the function $K(x,t)$ on $x \in [{x_i},{x_{i + 1}}]$ becomes
By the method of enlarging and reducing, we have
Thus, the form of remainder formula is
Let $h = \max \{ {h_1},{h_2}, \cdot{\cdot} \cdot ,{h_k}\}$, we have
Particularly, if the third derivative of the original function is a bounded variation function, and the interval partition satisfies
Appendix 2
When $a \ne 0$, let’s begin to consider the second order differential equation as follow
by the change of variables $t = x + ({b / {2a}})$, we have where $C = c - [{{{b^2}} / {(4a)}}]$.Further, the change of variables $z = \sqrt { - a} {t^2}$ and $\textrm{exp} (\textrm{ - }{z / 2})w(z) = y(t)$, leads to
Recall the following differential formulas for the Kummer function:
Therefore,
The Airy functions and their derivatives can be evaluated in Matlab as
When $a = b = 0$, Eq. (30) becomes
which has two linearly independent solutionsFunding
Ministry of Science and Technology of the People's Republic of China (2018YFB1107600); Natural Science Foundation of Zhejiang Province (Q20A010007); First Class Discipline of Zhejiang - A (Zhejiang Gongshang University - Statistics).
Disclosures
The authors declare no conflicts of interest.
References
1. J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985). [CrossRef]
2. R. März, Integrated Optics: Design and Modeling (Artech House, 1995).
3. J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992).
4. A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985). [CrossRef]
5. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photon. 1(1), 58–106 (2009). [CrossRef]
6. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994). [CrossRef]
7. A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014). [CrossRef]
8. A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).
9. 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]
10. Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993). [CrossRef]
11. P. J. Ciang, 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]
12. S. F. Chiang, B. Y. Lin, C. H. Teng, C. Y. Wang, and S. Y. Chung, “A multidomain pseudospectral mode solver for optical waveguide analysis,” J. Lightwave Technol. 30(13), 2077–2087 (2012). [CrossRef]
13. M. Walz, T. Zebrowski, J. Kuchenmeister, and K. Büsch, “B-spline modal method: A polynomial approach compared to the Fourier modal method,” Opt. Express 21(12), 14683–14697 (2013). [CrossRef]
14. J. Zhu and Y. Lu, “Leaky modes of slab waveguides-asymptotic solutions,” J. Lightwave Technol. 24(3), 1619–1623 (2006). [CrossRef]
15. L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005). [CrossRef]
16. C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14(24), 11631–11652 (2006). [CrossRef]
17. D. Song and Y. Lu, “Pesudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32(8), 1624–1630 (2014). [CrossRef]
18. J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008). [CrossRef]
19. S. Khorasani and K. Mehrany, “Differential transfer-matrix method for solution of one-dimensional linearnonhomogeneous optical structures,” J. Opt. Soc. Am. B 20(1), 91–96 (2003). [CrossRef]
20. M. Eghlidi, K. Mehrany, and B. Rashidian, “Modified differential transfer-matrix method for solution ofonedimensional linear inhomogeneous optical structures,” J. Opt. Soc. Am. B 22(7), 1521–1528 (2005). [CrossRef]
21. J. Zhu and Z. Shen, “Dispersion relation of leaky modes in nonhomogeneous waveguides and its applications,” J. Lightwave Technol. 29(21), 3230–3236 (2011). [CrossRef]
22. J. Zhu and J. Zheng, “Nonlinear equation of the modes in circular slab waveguides and its application,” Appl. Opt. 52(33), 8013–8023 (2013). [CrossRef]
23. Y. Li and J. Zhu, “Efficient approximations of dispersion relations in optical waveguides with varying refractive-index profiles,” Opt. Express 23(9), 11952–11964 (2015). [CrossRef]
24. S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973). [CrossRef]
25. D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956). [CrossRef]
26. J. Peetre, Lecture Notes: A Theory of Interpolation of Normed Spaces (Brazilia, 1963).