Smooth approximation of a varying refractive-index profile and its application in the computation of light waves

Jianxin Zhu; Jianxin Zhu; Peipei Yu; Zhaochen Zhu; Zhaochen Zhu

doi:10.1364/OSAC.418185

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

n = \left\{ {\begin{array}{lc} {{n_1}},& {x < 0},\\ {{n_0}(x)},&{0 \le x \le d,}\\ {{n_2}},& {x > d,} \end{array}} \right.

where ${n_0}(x)$ is a continuous function with ${n_0}(x) > \max ({n_1},{n_2})$, and $x = 0$ and $x = d$ are two discontinuous points. The Sturm-Liouville eigenvalue equation is

(1)$$\left\{ \begin{array}{l} \mathrm{\rho} \frac{\rm d }{{{\rm d} x}}\left( {\frac{1}{\rho }\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}}} \right) + \mathrm{\kappa}_0^2n_0^2(x)\mathrm{\phi} = {\mathrm{\beta}^2}\mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - \infty < x < \infty ,\\ \mathop {\lim }\limits_{|x |\to + \infty } \mathrm{\phi} = 0. \end{array} \right.$$

with $\mathrm{\rho} = 1$ for TE case and $\mathrm{\rho} = n_0^2(x)$ for TM case, where ${\mathrm{\kappa} _0}$ is the free space wavenumber, ${\mathrm{\beta} ^2}$ is an eigenvalue, and $\mathrm{\phi} (x)$ is an eigenfunction corresponding to ${\mathrm{\beta} ^2}$. For simplification, we only study the TE case

Fig. 1. Subdivision diagram of a three-layers’ waveguide (Ref. [23], Fig. 1), where $d = k \cdot h, k = 3$, and $h$ is a step width.

Download Full Size | PDF

In the TE case, Eq. (1) in $[0,d]$ is given by

(2)$$\frac{{{\rm{d} ^2}\mathrm{\phi} }}{{{\rm d} {x^2}}} + \mathrm{\kappa} _0^2n_0^2(x)\mathrm{\phi} = {\mathrm{\beta} ^2}\mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 < x < d,$$

(3)$$\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}} ={-} i \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}} \mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x = {0^ + },$$

(4)$$\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}} = i \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}} \mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{{cc}} {}&{} \end{array}{\kern 1pt} x = {\rm{d}^ - }{\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt}$$

where $i = \sqrt { - 1}$.

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:

(5)$$\frac{{{\textrm{\rm{d}}^2}{y_j}}}{{\textrm{\rm{d}}{x^2}}} + ({a_j}{x^2} + {b_j}x + {c_j}){y_j} = 0,{\kern 1pt}$$

where $(j = 1,2, \cdot{\cdot} \cdot ,k)$.

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

(6)$$V_{{x_i}}^{{x_{i + 1}}}({n_0^2(x)} )= \frac{1}{k}V_0^d({n_0^2(x)} ),$$

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

(7)$$\left\{ \begin{array}{l} {a_j} = \frac{{2{{\kappa}_0}^2}}{{{h^2}}}[{n_0}^2({t_0}) - 2{n_0}^2({t_1}) + {n_0}^2({t_2})],\\ {b_j} = \frac{{{{\kappa}_0}^2}}{{{h^2}}}[(1 - 4j){n_0}^2({t_0}) + (8j - 4){n_0}^2({t_1}) + (3 - 4j){n_0}^2({t_2})],\\ {c_j} = {{\kappa}_0}^2[{(2{j^2} - j){n_0}^2({t_0}) + 4(j - {j^2}){n_0}^2({t_1}) + (2{j^2} - 3j + 1){n_0}^2({t_2})} ]- {{\beta}^2}. \end{array} \right.$$

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

(8)$$\left\{ \begin{array}{l} {a_j} = \frac{{{M_{j + 1}} - {M_j}}}{{2({x_{j + 1}} - {x_j})}},\\ {b_j} = \frac{{{M_j}{x_{j + 1}} - {M_{j + 1}}{x_j}}}{{{x_{j + 1}} - {x_j}}},\\ {c_j} = \frac{{{M_{j + 1}}{x_j}^2 - {M_j}{x_{j + 1}}^2}}{{2({x_{j + 1}} - {x_j})}} + \frac{1}{2}{M_j}({x_{j + 1}} - {x_j}) + {\kappa}_0^2n_0^2({x_j}) - {{\beta}^2}, \end{array} \right.$$

where ${M_j} = S^{\prime}({x_j})$, $(j = 1,2, \cdots ,k)$, and $S(x)$ is the quadratic spline interpolation function.

Furthermore, ${M_j}$ and ${M_{j + 1}}$ satisfies

(9)$${M_j} + {M_{j + 1}} = \frac{{2{{\kappa} _0}^2({n_0}^2({x_{j + 1}}) - {n_0}^2({x_j}))}}{{{x_{j + 1}} - {x_j}}}.$$

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:

(10)$${y_j}(x) = {A_j}{\alpha _j}(x) + B{}_j{{\beta} _j}(x)\;\;\;\;,\;\;\;\;\;\;(j = 1,2, \cdot \cdot \cdot ,k).$$

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

(11)$${y_j}({x_j}) = {y_{j + 1}}({x_j})\;\;\;\;,\;\;\;\;\;\;{y^{\prime}_j}({x_j}) = {y^{\prime}_{j + 1}}({x_j}).$$

And the outgoing boundary conditions at $x = 0$ and $x = d$ become

(12)$${y^{\prime}_1}(0) ={-} i \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}} {y_1}(0), $$

(13)$${y^{\prime}_k}(d) = i \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}} {y_k}(d). $$

According to these conditions, a linear system of ${A_j}$ and ${B_j}(j = 1,2, \cdot{\cdot} \cdot ,k)$ can be obtained.

(14)$$\left\{ \begin{array}{l} {{\mathrm{\alpha}^{\prime}}_1}(0){A_1} + {{{\mathrm{\beta}}^{\prime}}_1}(0){B_1} ={-} i \sqrt {\mathrm{\kappa}_0^2n_1^2 - {{\mathrm{\beta}}^2}} [{{\mathrm{\alpha}_1}(0){A_1} + {{\mathrm{\beta}}_1}(0){B_1}} ],\\ {A_j}{\mathrm{\alpha}_j}({x_{j + 1}}) + {B_j}{{\mathrm{\beta}}_j}({x_{j + 1}}) = {A_{j + 1}}{\mathrm{\alpha}_{j + 1}}({x_{j + 1}}) + {B_{j + 1}}{{\mathrm{\beta}}_{j + 1}}({x_{j + 1}}),\\ {A_j}{{\mathrm{\alpha}^{\prime}}_j}({x_{j + 1}}) + {B_j}{{{\mathrm{\beta}}^{\prime}}_j}({x_{j + 1}}) = {A_{j + 1}}{{\mathrm{\alpha}^{\prime}}_{j + 1}}({x_{j + 1}}) + {B_{j + 1}}{{{\mathrm{\beta}}^{\prime}}_{j + 1}}({x_{j + 1}}),\\ {{\mathrm{\alpha}^{\prime}}_k}(d){A_k} + {{{\mathrm{\beta}}^{\prime}}_k}(d){B_k} = i \sqrt {\mathrm{\kappa}_0^2n_2^2 - {{\mathrm{\beta}}^2}} [{{\mathrm{\alpha}_k}(d){A_k} + {{\mathrm{\beta}}_k}(d){B_k}} ]. \end{array} \right.$$

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

(15)$$\left\{ \begin{array}{l} {T_1} ={-} \frac{{{{\alpha^{\prime}}_1}(0) + i {\alpha_1}(0){\gamma_1}}}{{{{{\beta}^{\prime}}_1}(0) + i {{\beta}_1}(0){\gamma_1}}},\\ {T_{j + 1}} ={-} \frac{{{\alpha_{j + 1}}(jh)[{{\alpha^{\prime}}_j}(jh) + {{{\beta}^{\prime}}_j}(jh){T_j}] + {{\alpha^{\prime}}_{j + 1}}(jh)[{\alpha_j}(jh) + {{\beta}_j}(jh){T_j}]}}{{{{\beta}_{j + 1}}(jh)[{{\alpha^{\prime}}_j}(jh) + {{{\beta}^{\prime}}_j}(jh){T_j}] + {{{\beta}^{\prime}}_{j + 1}}(jh)[{\alpha_j}(jh) + {{\beta}_j}(jh){T_j}]}},\\ \begin{array}{{ccc}} {}&{(j = 1,2, \cdots ,k - 1),}&{} \end{array}\\ {T_k} ={-} \frac{{{{\alpha^{\prime}}_k}(d) + i {\alpha_k}(d){\gamma_2}}}{{{{{\beta}^{\prime}}_k}(d) + i {{\beta}_k}(d){\gamma_2}}}. \end{array} \right.$$

where ${\gamma _1} = \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}}$ and ${\gamma _2} = \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}}$, and the evaluation on ${T_k}$ depends on ${T_1}$ recursively. Furthermore, ${T_k}$ is determined by ${\gamma _1}$. In addition, ${T_k}$ is given by the last equation, which depends on ${\gamma _2}$. Hence, the dispersion relation with respect to ${{\beta} ^2}$ is as follow:

(16)$$F({\beta} ) = {T_k} + \frac{{{{\alpha ^{\prime}}_k}(d) + i {\alpha _k}(d){\gamma _2}}}{{{{{\beta} ^{\prime}}_k}(d) + i {{\beta} _k}(d){\gamma _2}}}.$$

2.3 Asymptotic solutions for the TE case

By the WKB method, we can obtain the following formulas [18].

(17)$$\sqrt {{\kappa} _0^2n_0^2({0^ + }) - {{\beta} ^2}} \approx \frac{W}{{{a_0}}} - \frac{{a{a_2}}}{{{W^2}}} - \frac{{{a^2}{a_3}}}{{{W^3}}} - \frac{{{a^3}{a_4}}}{{{W^4}}}\textrm{ - } \cdot{\cdot} \cdot ,$$

where

{a_0} ={-} \frac{{i d}}{2}\;\;;{a_2} = - \frac{{{\delta _0}\textrm{ + }{\delta _1}\textrm{ + }{\delta _2}}}{8} ;{a_3} = \frac{i }{{4d}}({{\delta_0}\textrm{ + }{\delta_1}\textrm{ + }{\delta_2}} );{a_4} = - \frac{{ - i d{b_4} - 2{a_3}}}{{2{a_0}}} - \frac{{{a_2}^2}}{2};

{b_4} = - \frac{1}{{128}}({5{\delta_1}^2 + 5{\delta_2}^2 + {\delta_0}^2 + 2{\delta_1}{\delta_2} + 2{\delta_0}{\delta_2} + 2{\delta_0}{\delta_1}} );{\delta _0} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_0^2(d^-)

{\delta _1} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_1^2\;\;\;;{\delta _2} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_2^2;{\delta _3} = {\delta _2} - {\delta _0},

W = {W_m} = \textrm{LambertW} (p,\frac{{ - {i ^{ - q}}}}{2} \cdot {({\delta _1}{\delta _3})^{\frac{1}{4}}}); m ={-} 4(p + 1) + q\;\;\;;(q = 0,1,2,3;p = - 1, - 2, \cdot \cdot \cdot )

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.

Fig. 2. The fitting function of the spline interpolation with $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).

Download Full Size | PDF

Fig. 3. The fitting function of the spline interpolation with non-equidistant partition as $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).

Download Full Size | PDF

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.

Fig. 4. The propagation constants of the piecewise interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of ${\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PDF

Fig. 5. The propagation constants of the spline interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PDF

Fig. 6. The propagation constants of the spline interpolation with non-equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PDF

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.

Fig. 7. The relative error of the piecewise interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PDF

Fig. 8. The relative error of the spline interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PDF

Fig. 9. the relative error of spline interpolation with non-equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PDF

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.

Fig. 10. The real part of eigenfunction with $k = 4$. Where the solid line (marked in blue) stands for the spline interpolation, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.

Download Full Size | PDF

Fig. 11. The imaginary part of eigenfunction with $k = 4$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.

Download Full Size | PDF

Fig. 12. The real part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.

Download Full Size | PDF

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.

Fig. 13. The imaginary part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.

Download Full Size | PDF

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:

(18)$$({\hat{\lambda }\textrm{ - }\lambda } )\left( {\int_{{t_1}}^{{t_2}} {{\mathrm{\phi}^2}d x} + \sum\limits_{i = 1}^2 {\frac{{{\mathrm{\phi}^2}({t_i})}}{{2{\gamma_i}(\lambda )}}} } \right) \approx{-} \sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {R(f(x)) \cdot {\mathrm{\phi} ^2}d x} } ,$$

where ${{\beta} ^2} = \lambda$, ${\hat{{\beta} }^2} = \hat{\lambda }$, ${\gamma _i}(\lambda ) = \sqrt {\lambda - {{\kappa} _0}^2{n_i}^2}$, $R(f(x))$ is the remainder of interpolation $t{}_0 = 0 \,\, \textrm{and} \,\, t{}_1 = d$

Further,

(19)$$|{\hat{\lambda }\textrm{ - }\lambda } |\le \frac{{2\sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {|{R(f(x))} |\cdot {{|\mathrm{\phi} |}^2}d x} } }}{{\left|{2\int_{{t_1}}^{{t_2}} {{\mathrm{\phi}^2}d x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Case 1: For the piecewise interpolation, we can obtain the remainder formula of the piecewise quadratic interpolation:

(20)$$R(f(x))\textrm{ = }{R_1}(f(x)) = \frac{{f^{\prime\prime}{^{\prime}}(\xi )}}{6}{\omega _3}(x),$$

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:

(21)$$|{\hat{\lambda }\textrm{ - }\lambda } |\le \frac{{\sqrt 3 {h^3}}}{{108}} \cdot \frac{{\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |\cdot \int_0^d {{{|\mathrm{\phi} |}^2}d x} }}{{\left|{2\int_0^d {{\mathrm{\phi}^2}\rm{d} x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Case 2: For the quadratic spline interpolation, we can obtain the remainder formula as follow [26]:

(22)$$R(f(x)) = {R_2}(f(x)) = \int_0^d {K(x,t)f^{\prime\prime}{^{\prime}}(t)} d t,$$

where

(23)$$K(x,t) = {R_2}\left\{ {\frac{{{{(x - t)}^2}_ + }}{{2!}}} \right\},\begin{array}{{cc}} {}&{\begin{array}{{cc}} {}&{} \end{array}} \end{array}$$

here ${R_2}\, (\cdot)$ is represented as the remainder of the quadratic spline interpolation.

The integral result of the function $K(x,t)$ on $x \in [{x_i},{x_{i + 1}}]$ becomes

(24)$$\int_{{x_i}}^{{x_{i + 1}}} {|{K(x,t)} |d t} = \frac{1}{6}{(x - {x_i})^2}({x_{i + 1}} - x).$$

By the method of enlarging and reducing, we have

(25)$$\int_0^d {|{K(x,t)} |d t} = \frac{1}{6}(x - {x_i})({x_{i + 1}} - x)(x - {x_i} + \sum\limits_{j = 1}^{i - 1} {{h_j}/{h_i}} ) < \frac{1}{{24}}{h_i}\sum\limits_{j = 1}^i {{h_j}^2} .$$

Thus, the form of remainder formula is

(26)$$R(f(x)) = {R_2}(f(x)) = \frac{1}{{24}}\int_0^d {{h_i}\sum\limits_{j = 1}^i {{h_j}^2} |{f^{\prime\prime}{^{\prime}}(t)} |} d t.\begin{array}{{cc}} {}&{} \end{array}$$

Let $h = \max \{ {h_1},{h_2}, \cdot{\cdot} \cdot ,{h_k}\}$, we have

(27)$$|{\hat{\lambda }\textrm{ - }\lambda } |< \frac{1}{{12}} \cdot \frac{{\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |\cdot \sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {{h_i}\sum\limits_{j = 1}^i {{h_j}^2} } } |{{\mathrm{\phi}^2}} |d x}}{{\left|{2\int_0^d {{\mathrm{\phi}^2}d x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Particularly, if the third derivative of the original function is a bounded variation function, and the interval partition satisfies

(28)$$\frac{1}{{{h_i}}}\sum\limits_{j = 1}^{i - 1} {|{{h_j} - {h_{j - 1}}} |} \le \frac{h}{\rm{d}},(i = 1,2, \cdot{\cdot} \cdot ,k),$$

a better residual estimate of the quadratic spline interpolation can be deduced.

(29)$$|{{R_2}(f(x))} |\le \frac{1}{8}{h^3}\left( {\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |+ V_0^d({f^{\prime\prime}{^{\prime}}(x)} )} \right).$$

Appendix 2

When $a \ne 0$, let’s begin to consider the second order differential equation as follow

(30)$$y^{\prime\prime}(x) + (a{x^2} + bx + c)y(x) = 0,$$

by the change of variables $t = x + ({b / {2a}})$, we have

(31)$$y^{\prime\prime}(t) + (a{t^2} + C)y(t) = 0,$$

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

(32)$$z\frac{{{\textrm{\rm{d}}^2}w}}{{\textrm{\rm{d}}{z^2}}} + \left( {\frac{1}{2} - z} \right)\frac{{\textrm{\rm{d}}w}}{{\textrm{\rm{d}}z}} + (\frac{1}{4} - \frac{C}{{4\sqrt { - a} }})w(z) = 0,$$

which is in the form of the confluent hypergeometric equation

(33)$$z\frac{{{\textrm{\rm{d}}^2}w}}{{\textrm{\rm{d}}{z^2}}} + ({B - z} )\frac{{\textrm{\rm{d}}w}}{{\textrm{\rm{d}}z}} + Aw(z) = 0,$$

with $A = \frac{1}{4}\textrm{ - }\frac{C}{{4\sqrt { - a} }}$ and $B = {1 / 2}$. Two linearly independent solutions are given by the Kummer functions $M(A,B,z)$ and $U(A,B,z)$. However, $U(A,B,z)$ shows a branch point $z = 0$. Hence, two linearly independent solutions of Eq. (30) are given by

(34)$$\left\{ \begin{array}{l} \alpha (x) = \textrm{exp} (\textrm{ - }{{z(x)} / 2}) \cdot M\left( {A,\frac{1}{2},z(x)} \right),\\ {\beta} (x) = \textrm{exp} (\textrm{ - }{{z(x)} / 2}) \cdot {[z(x)]^{\frac{1}{2}}}M\left( {A + \frac{1}{2},\frac{3}{2},z(x)} \right) \end{array} \right.,{\kern 1pt}$$

where $z(x) = \sqrt { - a} {[{x + {b / {(2a)}}} ]^2}$, $A = \frac{C}{{4\sqrt { - a} }} = \frac{{{b^2} - 4ac}}{{16a\sqrt { - a} }}.$

Recall the following differential formulas for the Kummer function:

(35)$$\left\{ \begin{array}{l} \alpha^{\prime}(x) ={-} \frac{1}{2}z^{\prime}(x)\alpha (x) + 2Az^{\prime}(x)\textrm{exp} (\textrm{ - }{{z(x)} / 2})M\left( {A + 1,\frac{3}{2},z(x)} \right),\\ {\beta}^{\prime}(x) ={-} \frac{{1 - z(x)}}{{2z(x)}}z^{\prime}(x){\beta} (x)\\ \begin{array}{{cc}} {\begin{array}{{cc}} {}&{} \end{array}}&{} \end{array} + \frac{{2A + 1}}{3}\sqrt {z(x)} z^{\prime}(x)\textrm{exp} [\textrm{ - }{{z(x)} / 2}]M\left( {A + \frac{3}{2},\frac{5}{2},z(x)} \right). \end{array} \right.$$

Special Case: when $a = 0$,and $b \ne 0$, Eq. (30) becomes

(36)$$y^{\prime\prime}(x) + (bx + c)y(x) = 0.$$

and two linearly independent solutions are given by the Airy functions

(37)$$\alpha (x) = {\rm Ai} (w(x)),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} (x) = {\rm Bi} (w(x)),$$

where $w(x) ={-} \frac{{bx + c}}{{{b^{2/3}}}}$.

Therefore,

(38)$$\alpha ^{\prime}(x) = \textrm{ - }{b^{1/3}}{\rm{A i}}^{\prime}(w(x)),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} ^{\prime}(x) = \textrm{ - }{b^{1/3}}{\rm{B i}}^{\prime}(w(x)).$$

The Airy functions and their derivatives can be evaluated in Matlab as

\textrm{airy} (z )= {\rm{Ai}} (z )\;\;\;,\;\;\;\textrm{airy}({1,z} )= {\rm{A i}}^{\prime}(z ),\\

\textrm{airy} ({2,z} )= {\rm{B i}} (z )\;\;,\;\;\;\textrm{airy} ({3,z} )= B i^{\prime}(z )

When $a = b = 0$, Eq. (30) becomes

(39)$$y^{\prime\prime}(x) + cy(x) = 0,$$

which has two linearly independent solutions

(40)$$\alpha (x) = {e^{ - i\sqrt c x}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} (x) = {e^{i\sqrt c x}}.$$

Funding

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).

Smooth approximation of a varying refractive-index profile and its application in the computation of light waves

Abstract

Corrections

1. Introduction

2. Mathematical modelling

2.1 Interval partitions

2.1.1 Equidistant partition

2.1.2 Adaptive partition

2.2 Interpolation methods

2.2.1 By the piecewise quadratic interpolation

2.2.2 By the quadratic spline interpolation

2.3 Asymptotic solutions for the TE case

3. Numerical example and comparison

3.1 Comparison of the fitting effect by two kinds of spline interpolation methods

3.2 Comparison of eigenvalues obtained by different spline interpolation methods

3.3 Numerical verification of the convergence rate

3.4 Comparison of eigenfunctions obtained by two kinds of interpolation methods

4. Conclusions

Appendix 1

Appendix 2

Funding

Disclosures

References

Cited By

Figures (13)

Equations (47)

OSA Continuum