Efficient approximations of dispersion relations in optical waveguides with varying refractive-index profiles

Yutian Li; Jianxin Zhu

doi:10.1364/OE.23.011952

1. Introduction

Many practical applications of dielectrical waveguides, such as the optimal designs of optical devices [1], depend on the propagation characteristics of these waveguides. When the eigen-mode expansion method is applied to the computation of wave propagation, the modal analysis and computation are crucial. This raised a strong interest in developing efficient methods to fast determine their characteristics with high-accuracy [2], due to the fact that exact solutions for the fields in dielectric waveguides exist only for few structures [3, 4]. For an open optical waveguide, the total electromagnetic field can be expressed as the sum of two components. One component describes the spatial steady state, and the other is the radiation field, which describes the spatial transient [5]. Mathematically, the field can be decomposed into a linear combination of a few discrete guided modes and an infinite integral of a continuum of radiation modes [6]. Common numerical solvers, such as the finite difference method [7], the finite element method [8–11], the multi-domain pseudo-spectral method [12–14], and the B-spline modal method [15], reduce the problem to a matrix eigenvalue problem. However, it is difficult to obtain eigenvalues with satisfactory accuracy since these methods all produce large and complex matrices and the corresponding eigenvalue problems are too complicated to solve.

The infinite integral in the decomposition of a field can be approximately expressed as an infinite sum of leaky modes through discretization [16,17]. The character constants of the guided modes are relatively easy to obtain by common solvers because they are real and bounded [18]. On the other hand, the character constants of the leaky modes are complex numbers [18], which leads to difficulties in computing the leaky modes, especially those with large modules, with high accuracy by the solvers of discretizing the Helmholtz operator [19]. For this reason, an alternative approach is proposed, that is to transform the characteristic problem to a nonlinear dispersion equation, whose roots correspond to the eigenvalues or the propagation constants. There are several efficient methods to deal with the piecewise constant slab waveguides [18,20,21] and rib waveguides [22]. While for the varying refractive-index waveguides, several approximate dispersion equations are constructed by the Wentzel-Kramers-Brillouin (WKB) method [23] and the differential transfer matrix method (DTMM) [24–27]. However, all these methods can not guarantee the accuracy of the modes. Recently, an exact dispersion equation of transverse electric (TE) modes for nonhomogeneous optical waveguides is established [19]. Unfortunately, this exact equation involves highly oscillatory integrals and derivatives of the refractive-index function, which bring difficulties to the numerical computation in practice.

In the present paper, we consider the eigenmode problems for varying refractive-index waveguides in both TE and TM cases. First, we approximate the refractive-index function by a piecewise polynomial of degree two. Then, the approximate eigenfunctions of the Helmholtz operator are expressed analytically by the Kummer functions. Finally, the approximate and analytical dispersion equations are established. The obtained dispersion equations can be solved by the secant method or Müller’s method. Numerical simulations show that our method has a third order accuracy, that is O(h³) with h standing for the step size. This suggests the error estimate originally proved for real eigenvalue problems in [28] is still valid for complex eigenvalue problems. The rest of this paper is organized as follows. In Section 2, the mathematical model and treatment are introduced; and the dispersion equations are established and the formulas of asymptotic solutions are provided for TE and TM cases, respectively. In Section 3, numerical results and corresponding analysis for the TE case are presented. Finally, the conclusions are given in Section 4.

2. Mathematical model and treatment

In this section, we consider a general two-dimensional waveguide shown in Fig. 1. Let the refractive-index function n(x) be

n (x) = {\begin{matrix} n_{1}, & x < 0; \\ n_{0} (x), & 0 < x < d; \\ n_{2}, & x > d . \end{matrix}

Here, n₀(x) is a given positive-valued function where n₀(x) > max(n₁, n₂), and there are two discontinuous points at x = 0 and x = d. The Sturm-Liouville eigenvalue equation is

ρ \frac{d}{d x} (\frac{1}{ρ} \frac{d ϕ}{d x}) + κ_{0}^{2} n^{2} (x) ϕ = β^{2} ϕ, - \infty < x < \infty,

with ρ = 1 for TE case and ρ = n²(x) for TM case. Here, κ₀ is the free space wavenumber, β² is an eigenvalue, and ϕ(x) is an eigenfunction corresponding to β². Since the solutions in media 1 and 2 are exponential functions that decrease away from the media discontinuity, the boundary conditions are

lim_{| x | \to \infty} ϕ = 0

.

Fig. 1 Subdivision diagram of a waveguide

Download Full Size | PDF

Our strategy is to approximate the function $n_{0}^{2} (x)$ piece-wisely to obtain a simpler equation that can be solved analytically. Given the function $n_{0}^{2} (x)$ , we could interpolate this function by a piecewise polynomial H(x) of degree two. Assume the interval [0, d] is divided into k subintervals [0, h], [h, 2h], ···, [(k − 1)h, kh] with h = d/k. On each subinterval, three nodes are chosen to be the two endpoints and the midpoint. In particular, H(x_j) = n²(x_j), (j = 0, 1, 2, ··· ,k).

2.1. TE case

In the TE case, the Sturm-Liouville problem is given by

\frac{d^{2} ϕ}{d x^{2}} + κ_{0}^{2} n_{0}^{2} (x) ϕ = β^{2} ϕ, 0 < x < d,

\frac{d ϕ}{d x} = - i γ_{1} ϕ, x = 0^{+},

\frac{d ϕ}{d x} = i γ_{2} ϕ, x = d^{-},

where

γ_{1} = {(k_{0}^{2} n_{1}^{2} - β^{2})}^{1 / 2}

,

γ_{2} = {(k_{0}^{2} n_{2}^{2} - β)}^{1 / 2}

, and

i = \sqrt{- 1}

. The boundary conditions in Eqs. (4)–(5) are often referred to as out-going conditions.

On each subinterval, the function $n_{0}^{2} (x)$ is interpolated by a polynomial of degree two, and hence the approximated equation of Eq. (3) is given by

\frac{d^{2} y_{j}}{d x^{2}} + (a_{j} x^{2} + b_{j} x + c_{j}) y_{j} (x) = 0, (j - 1) h < x < j h,

where

{\begin{cases} a_{j} = 2 κ_{0}^{2} / h^{2} [n_{0}^{2} (x_{0}) - 2 n_{0}^{2} (x_{1}) + n_{0}^{2} (x_{2})], \\ b_{j} = κ_{0}^{2} / h [(1 - 4 j) n_{0}^{2} (x_{0}) + (8 j - 4) n_{0}^{2} (x_{1}) + (3 - 4 j) n_{0}^{2} (x_{2})], \\ c_{j} = κ_{0}^{2} [(2 j^{2} - j) n_{0}^{2} (x_{0}) + 4 (j - j^{2}) n_{0}^{2} (x_{1}) + (2 j^{2} - 3 j + 1) n_{0}^{2} (x_{2})] - β^{2} . \end{cases}

Here (j = 1, 2,··· ,k), and x₀ = (j − 1)h, x₁ = (j − 1/2)h, x₂ = jh.

The approximate Eq. (6) is equivalent to the confluent hypergeometric equation and its solution can be given by the Kummer functions; see Appendixes A and B. Let {α_j(x), β_j(x)} be a pair of linearly independent solutions on j-th interval. Then,

y_{j} (x) = A_{j} α_{j} (x) + B_{j} β_{j} (x); j = 1, 2, \dots, k,

where A_j and B_j are constants. To fix these constants, we need to impose two conditions at each interface x = jh, where we require the solution y(x) and its first order derivative are continuous. That is,

y_{j} (j h) = y_{j + 1} (j h), {\frac{d y_{j}}{d x} |}_{x = j h} = {\frac{d y_{j + 1}}{d x} |}_{x = j h}, (j = 1, 2, \dots k - 1) .

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

{\frac{d y_{1}}{d x} |}_{x = 0} = - i {(k_{0}^{2} n_{1}^{2} - β^{2})}^{1 / 2} y_{1} (0)

{\frac{d y_{k}}{d x} |}_{x = d} = i {(k_{0}^{2} n_{2}^{2} - β^{2})}^{1 / 2} y_{k} (d) .

Thus, we obtain a linear system of A_j and B_j:

{\begin{cases} α_{1}^{'} (0) A_{1} + β_{1}^{'} (0) B_{1} = - i {(k_{0}^{2} n_{1}^{2} - β^{2})}^{1 / 2} [α_{1} (0) A_{1} + β_{1} (0) B_{1}], \\ α_{j} (j h) A_{j} + β_{j} (j h) B_{j} - α_{j + 1} (j h) A_{j + 1} - β_{j + 1} (j h) B_{j + 1} = 0, \\ α_{j}^{'} (j h) A_{j} + β_{j}^{'} (j h) B_{j} - α_{j + 1}^{'} (j h) A_{j + 1} - β_{j + 1}^{'} (j h) B_{j + 1} = 0, \\ α_{k}^{'} (d) A_{k} + β_{k}^{'} (d) B_{k} = i {(k_{0}^{2} n_{2}^{2} - β^{2})}^{1 / 2} [α_{k} (d) A_{k} + β_{k} (d) B_{k}], \end{cases}

where j = 1, 2,··· ,k − 1. Denote T_j = B_j/A_j, (j = 1, 2, ··· ,k). The above linear system leads to the dispersion relation of β as follows:

{\begin{cases} T_{1} = - \frac{α_{1}^{'} (0) + i α_{1} (0) γ_{1}}{β_{1}^{'} (0) + i β_{1} (0) γ_{1}}, \\ T_{j + 1} = - \frac{α_{j + 1} (j h) [α_{j}^{'} (j h) + β_{j}^{'} (j h) T_{j}] - α_{j + 1}^{'} (j h) [α_{j} (j h) + β_{j} (j h) T_{j}]}{β_{j + 1} (j h) [α_{j}^{'} (j h) + β_{j}^{'} (j h) T_{j}] - β_{j + 1}^{'} (j h) [α_{j} (j h) + β_{j} (j h) T_{j}]}, \\ T_{k} = - \frac{α_{k}^{'} (d) - i α_{k} (d) γ_{2}}{β_{k}^{'} (d) - i β_{k} (d) γ_{2}}, \end{cases}

where (j = 1,··· ,k − 1). Note that T_k is determined from T₁ recursively, and hence depends on

γ_{1} = {(k_{0}^{2} n_{1}^{2} - β^{2})}^{1 / 2}

. On the other hand, T_k is given by the last equation, and hence depends on

γ_{2} = {(k_{0}^{2} n_{2}^{2} - β^{2})}^{1 / 2}

. Therefore, we obtain an equation (a dispersion relation) for propagation constant β.

Special case (slab waveguide) for TE case: If n₀(x) = n₀ (constant) and k = 1, then α₁(x) = exp(−iγ₀x) and β₁(x) = exp(iγ₀x), with $γ_{0} = {(κ_{0}^{2} n_{0}^{2} - β^{2})}^{1 / 2}$ . By

- \frac{α_{1}^{'} (d) - i α_{1} (d) γ_{2}}{β_{1}^{'} (d) - i β_{1} (d) γ_{2}} = - \frac{α_{1}^{'} (0) + i α_{1} (0) γ_{1}}{β_{1}^{'} (0) + i β_{1} (0) γ_{1}},

we have

exp (- 2 i γ_{0} d) \cdot \frac{γ_{0} + γ_{2}}{γ_{2} - γ_{0}} = \frac{γ_{1} - γ_{0}}{γ_{1} + γ_{0}},

that is

\frac{(γ_{0} + γ_{1}) (γ_{0} + γ_{2})}{(γ_{0} - γ_{1}) (γ_{0} - γ_{2})} = exp (2 i γ_{0} d) .

This dispersion relation is the same as the Eq. (19) in [18].

The dispersion relation given in Eq. (10) for the propagation constant β is nonlinear. In order to evaluate β numerically, we shall adopt Müller’s method for implementation [29]. Müller’s method requires three starting values, for which we use the WKB approximations obtained in [23].

Asymptotic solutions for TE case: By the WKB method, we can obtain the following formulas [23].

{[κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}]}^{1 / 2} \approx \frac{W}{a_{0}} - \frac{a_{0} a_{2}}{W^{2}} - \frac{a_{0}^{2} a_{3}}{W^{3}} - \frac{a_{0}^{3} a_{4}}{W^{4}} + \dots,

where

\begin{array}{l} a_{0} = - \frac{i d}{2}, a_{2} = - \frac{δ_{0} + δ_{1} + δ_{2}}{8}, a_{3} = \frac{i}{4 d} (δ_{0} + δ_{1} + δ_{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 δ_{1}^{2} + 5 δ_{2}^{2} + δ_{0}^{2} + 2 δ_{2} δ_{1} + 2 δ_{2} δ_{0} + 2 δ_{1} δ_{0}), \\ δ_{0} = κ_{0}^{2} n_{0}^{2} (0^{+}) - κ_{0}^{2} n_{0}^{2} (d^{-}), δ_{1} = κ_{0}^{2} n_{0}^{2} (0^{+}) - κ_{0}^{2} n_{1}^{2}, \\ δ_{2} = κ_{0}^{2} n_{0}^{2} (0^{+}) - κ_{0}^{2} n_{2}^{2}, δ_{3} = δ_{2} - δ_{0}; \end{array}

and

W = LambertW (p, - \frac{i^{- q} d}{4} (δ_{1} δ_{3}))

, (q = 0, 1, 2, and 3), (p = −1, −2,···).

In the approximation formula given in Eq. (11), $Im [κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}] < 0$ is required. Notice that the Lambert W function is the inverse function of x = We^W and it is multi-valued. On the complex plane, the equation We^W = x has a countably infinite number of solutions. They are represented by LambertW(p, x) with p ranging over the integers and standing for different branches. Since the leaky modes β are in the first quadrant, the value of p should be a negative integer, and the modulus of Lambert(p, x) increases as p decreases. Therefore, smaller negative integer p gives better approximation to the leaky modes. The integer −q is the power of i, and we simply take q = 0, 1, 2, 3, which give all the four values of i to an integer power.

For Müller’s method, we shall use the first two, first three and first four terms in the WKB approximations as the initial values.

2.2. TM case

For the TM case, the Sturm-Liouville eigenvalue problem becomes

n_{0}^{2} (x) \frac{d}{d x} (\frac{1}{n_{0}^{2} (x)} \frac{d ϕ}{d x}) + κ_{0}^{2} n^{2} (x) ϕ = β^{2} ϕ, 0 < x < d,

\frac{d ϕ}{d x} = - i γ_{1} \frac{n_{0}^{2} (0^{+})}{n_{1}^{2}} ϕ, at x = 0,

\frac{d ϕ}{d x} = i γ_{2} \frac{n_{0}^{2} (d^{-})}{n_{2}^{2}} ϕ, at x = d .

In order to apply the same method as in the TE case, we need to simplify Eq. (12). Let p(x) = ϕ(x)/n₀(x) (see [23]), then Eqs. (12)–(14) become

p_{x x} + (κ_{0}^{2} s^{2} (x) - β^{2}) p = 0, 0 < x < d,

p_{x} = [- i γ_{1} \frac{n_{0}^{2} (0^{+})}{n_{1}^{2}} - \frac{n_{0 x} (0^{+})}{n_{0} (0^{+})}] p ≜ \hat{A} (β) p, at x = 0,

p_{x} = [i γ_{2} \frac{n_{0}^{2} (d^{-})}{n_{2}^{2}} - \frac{n_{0 x} (d^{-})}{n_{0} (d^{-})}] p ≜ \hat{B} (β) p, at x = d,

where

s^{2} (x) = n_{0}^{2} + \frac{n_{0 x x}}{κ_{0}^{2} n_{0}} - \frac{2 n_{0 x}^{2}}{κ_{0}^{2} n_{0}^{2}},

and n_0x and n_0xx denote the first and second derivatives of n₀(x), respectively; see Eqs. (19)–(21) of [23].

To solve the problem in Eqs. (15)–(17), we use the same method as we have done for the TE case in Section 2.1. The interval [0, d] is divided into k equal subintervals, and on each subinterval the function $k_{0}^{2} s^{2} (x) - β^{2}$ is approximated by a polynomial of degree two. The rest is almost straightforward: rewriting the equations in Section 2.1. Only need to notice the differences in the boundary conditions between Eqs. (16)–(17) and Eqs. (4)–(5).

On each subinterval, we also have an equation similar to Eq. (6)

\frac{d^{2} p_{j}}{d x^{2}} + ({\hat{a}}_{j} x^{2} + {\hat{b}}_{j} x + {\hat{c}}_{j}) p_{j} (x) = 0, (j - 1) h < x < j h,

where

{\begin{cases} {\hat{a}}_{j} = \frac{2 κ_{0}^{2}}{h^{2}} [s^{2} (x_{0}) - 2 s^{2} (x_{1}) + s^{2} (x_{2})], \\ {\hat{b}}_{j} = \frac{κ_{0}^{2}}{h} [(1 - 4 j) s^{2} (x_{0}) + (8 j - 4) s^{2} (x_{1}) + (3 - 4 j) s^{2} (x_{2})], \\ {\hat{c}}_{j} = κ_{0}^{2} [(2 j^{2} - j) s^{2} (x_{0}) + 4 (j - j^{2}) s^{2} (x_{1}) + (2 j^{2} - 3 j + 1) s^{2} (x_{2})] - β^{2} . \end{cases}

Here (j = 1, 2,··· ,k), and x₀ = (j − 1)h, x₁ = (j − 1/2)h, x₂ = jh. Solutions of (19) are given by the Kummer functions

p_{j} (x) = {\hat{A}}_{j} α_{j} (x) + {\hat{B}}_{j} β_{j} (x); j = 1, 2, \dots, k .

We again impose the following interface conditions at x = jh:

p_{j} (j h) = p_{j + 1} (j h), {\frac{d p_{j}}{d x} |}_{x = j h} = {\frac{d p_{j + 1}}{d x} |}_{x = j h}, (j = 1, 2, \dots, k - 1) .

The outgoing boundary conditions at x = 0 and x = d now become

{\frac{d p_{1}}{d x} |}_{x = 0} = \hat{A} (β) p_{1} (0),

{\frac{d p_{k}}{d x} |}_{x = d} = \hat{B} (β) p_{k} (d) .

Thus, we obtain a linear system of Â_j and B̂_j:

{\begin{cases} α_{1}^{'} (0) {\hat{A}}_{1} + β_{1}^{'} (0) {\hat{B}}_{1} = \hat{A} (β) [α_{1} (0) {\hat{A}}_{1} + β_{1} (0) {\hat{B}}_{1}], \\ α_{j} (j h) {\hat{A}}_{j} + β_{j} (j h) {\hat{B}}_{j} - α_{j + 1} (j h) {\hat{A}}_{j + 1} - β_{j + 1} (j h) {\hat{B}}_{j + 1} = 0, \\ α_{j}^{'} (j h) {\hat{A}}_{j} + β_{j}^{'} (j h) {\hat{B}}_{j} - α_{j + 1}^{'} (j h) {\hat{A}}_{j + 1} - β_{j + 1}^{'} (j h) {\hat{B}}_{j + 1} = 0, \\ α_{k}^{'} (d) {\hat{A}}_{k} + β_{k}^{'} (d) {\hat{B}}_{k} = \hat{B} (β) [α_{k} (d) {\hat{A}}_{k} + β_{k} (d) {\hat{B}}_{k}], \end{cases}

where (j = 1, 2,··· ,k − 1). Denote T̂_j = B̂_j/Â_j, (j = 1, 2, ··· ,k). We then obtain the dispersion relation of β as follows:

{\begin{cases} {\hat{T}}_{1} = - \frac{α_{1}^{'} (0) - \hat{A} (β) α_{1} (0)}{β_{1}^{'} (0) - \hat{A} (β) β_{1} (0)}, \\ {\hat{T}}_{j + 1} = - \frac{α_{j + 1} (j h) [α_{j}^{'} (j h) + β_{j}^{'} (j h) {\hat{T}}_{j}] - α_{j + 1}^{'} (j h) [α_{j} (j h) + β_{j} (j h) {\hat{T}}_{j}]}{β_{j + 1} (j h) [α_{j}^{'} (j h) + β_{j}^{'} (j h) {\hat{T}}_{j}] - β_{j + 1}^{'} (j h) [α_{j} (j h) + β_{j} (j h) {\hat{T}}_{j}]}, \\ {\hat{T}}_{k} = - \frac{α_{k}^{'} (d) - \hat{B} (β) α_{k} (d)}{β_{k}^{'} (d) - \hat{B} (β) β_{k} (d)}, \end{cases}

where (j = 1,··· ,k − 1).

Special case (slab waveguide) for TM case: If n₀(x) = n₀ (constant) and k = 1, then α₁(x) = exp(−iγ₀x), β₁(x) = exp(iγ₀x), $\hat{A} (β) = - i γ_{1} n_{0}^{2} / n_{1}^{2}$ , and $\hat{B} (β) = i γ_{2} n_{0}^{2} / n_{2}^{2}$ . By

- \frac{α_{1}^{'} (d) - \hat{B} (β) α (d)}{β_{1}^{'} (d) - \hat{B} (β) β_{1} (d)} = - \frac{α_{1}^{'} (0) - \hat{A} (β) α_{1} (0)}{β_{1}^{'} (0) - \hat{A} (β) β_{0}},

we have

exp (- 2 i γ_{0} d) \cdot \frac{\frac{n_{0}^{2}}{n_{2}^{2}} γ_{2} + γ_{0}}{\frac{n_{0}^{2}}{n_{2}^{2}} γ_{2} - γ_{0}} = \frac{\frac{n_{0}^{2}}{n_{1}^{2}} γ_{1} - γ_{0}}{\frac{n_{0}^{2}}{n_{1}^{2}} γ_{1} + γ_{0}},

that is,

\frac{(μ_{0} + μ_{1}) (μ_{0} + μ_{2})}{(μ_{0} - μ_{1}) (μ_{0} - γ_{2})} = exp (2 i γ_{0} d),

where

μ_{j} = γ_{j} / n_{j}^{2}

, (j = 1, 2). This dispersion relation is the same as Eq. (22) in [18].

We shall adopt Müller’s method to solve the dispersion Eq. (22), where the derivatives of n₀(x) can be approxomated by the higher-order finite difference, see Appendix C. For the initial values, we shall make use of the WKB approximations obtained in [23].

Asymptotic solutions for TM case: By the WKB method, we can obtain the following formulas [23].

(1) Leading-term approximation: ${[κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}]}^{1 / 2} \approx \frac{π m}{d} - \frac{i ln \tilde{A_{0}}}{2 d} - \tilde{b_{2}} ≜ K_{m};$
(2) Two-term Approximation: ${[κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}]}^{1 / 2} \approx K_{m} + \tilde{b_{2}} - \frac{\tilde{A_{i}}}{2 d \tilde{A_{0}} K_{m}} ≃ K_{m 1};$
(3) Three-term approximation: ${[κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}]}^{1 / 2} \approx K_{m 1} - \frac{i (2 \tilde{A_{2}} \tilde{A_{0}} - {\tilde{A_{1}}}^{2})}{2 d {\tilde{A_{0}}}^{2} K_{m 1}^{2}};$ where m = 1, 2, 3,···, $\tilde{A_{0}} = \tilde{c_{0}} e_{0}$ , $\tilde{A_{1}} = \tilde{c_{1}} e_{0} + \tilde{c_{0}} e_{1}$ , $\tilde{A_{2}} = \tilde{c_{0}} e_{2} + \tilde{c_{1}} e_{1} + \tilde{c_{2}} e_{0}$ , and $\begin{array}{l} \tilde{a_{2}} = δ_{0} - γ_{x x} (d^{-}) + γ_{x}^{2} (d^{-}), \tilde{b_{2}} = γ_{x}^{2} (0^{+}) - γ_{x x} (0^{+}), γ (x) = ln n_{0} (x), \\ \tilde{c_{2}} = \frac{n_{2}^{2} + n_{0}^{2} (d^{-})}{n_{2}^{2} - n_{0}^{2} (d^{-})}, \tilde{c_{1}} = - \frac{2 i γ_{x} (d^{-}) n_{2}^{4}}{{(n_{2}^{2} - n_{0}^{2} (d^{-}))}^{2}}, \tilde{c_{2}} = - \frac{2 n_{2}^{6} γ_{x}^{2} (d^{-})}{{(n_{2}^{2} - n_{0}^{2} (d^{-}))}^{3}} + \frac{(\tilde{a_{2}} - δ_{2}) n_{0}^{2} (d^{-}) n_{2}^{2}}{{(n_{2}^{2} - n_{0}^{2} (d^{-}))}^{2}}, \\ e_{0} = \frac{n_{1}^{2} + n_{0}^{2} (0^{+})}{n_{1}^{2} - n_{0}^{2} (0^{+})}, e_{1} = \frac{2 i γ_{x} (0^{+}) n_{1}^{4}}{{(n_{1}^{2} - n_{0}^{2} (0^{+}))}^{2}}, e_{2} = - \frac{2 n_{1}^{6} γ_{x}^{2} (0^{+})}{{(n_{1}^{2} - n_{0}^{2} (0^{+}))}^{3}} + \frac{(\tilde{b_{2}} - δ_{1}) n_{0}^{2} (0^{+}) n_{1}^{2}}{{(n_{1}^{2} - n_{0}^{2} (0^{+}))}^{2}} . \end{array}$ Here, $Im [κ_{0}^{2} n_{0}^{2} (0^{+}) - β^{2}] < 0$ is also required.

3. Numerical example

The eigenmodes of an open waveguide consist of two sets, namely the guided modes and the leaky modes. The propagation constants β of the guided modes are all real and satisfy $κ_{0}^{2} max {n_{1}^{2}, n_{2}^{2}} \leq β^{2} < κ_{0}^{2} max_{0 \leq x \leq d} n_{0}^{2} (x)$ , see [18]. Thus, the propagation constants are easy to obtain through by finding the roots of Eq. (10) with Mesh Generation method [30] or Müller’s method with some suitable initial values of β² chosen in this interval. So, we omit finding the propagation constants of the guided modes here.

To verify the effectiveness of our method, for simplicity, we use the initial values of β by asymptotic solutions which correspond to leaky modes.

Example. Let n₀(x) = 3.5[1 − 0.08(x − 1.0)²], n₁ = 2.10, n₂ = 3.17, d = 2.0cm, λ = 1.55μm and κ₀ = 2π/λ.

For simplicity, we only give the numerical simulation for the TE case since the TM case is similar. By Eq. (10), we set

g = T_{k} + \frac{α_{k}^{'} (d) - i α_{k} (d) γ_{2}}{β_{k}^{'} (d) - i β_{k} (d) γ_{2}},

where T_k is obtained from T₁ by a recursive manner. To find the roots of the nonlinear equation g = 0, we use Müller’s method, where three initial values are taken as the two-term, three-term and four-term asymptotic approximations by the WKB method given in Eq. (11). The iteration stopping criteria is chosen as |g| < 10⁻⁸. For the two indices p and q appearing in the Lambert W function, we choose p = −2, −3,···, and q = 0, 1, 2, 3, since p = −1 does not correspond to any actual leaky mode, ant it should be removed; see the paragraph under Eq. (11). For the convenience of labeling, we set m = −4(p + 2) + (q + 1), with p = −2, −3,···, and q = 0, 1, 2, 3. The number of subintervals k is chosen to be 2, 4, 8, 16 for the simulation and part of the numerical results are listed in Table 1.

Table 1. Propagation constants of Example for TE case, k = 2, 4, 8, 16.

View Table

The propagation constants are plotted in Fig. 2. We use the computed values from k = 8 and k = 16. As we can see from Table 1, most of the values are nearly the same for k = 8 and k = 16 when m is large, and so they coincide in the figure. In order to examine the convergence order, we choose three different β values, and compute the dispersion equation with k = 2, 4, 8, 16, 32, 64, 128 respectively. The reference solution is obtained with k = 256 and a higher accuracy |g| < 10⁻¹⁰. These reference values of propagation constants are β̂₁ = 10.4184075001 + 0.8691830321i, β̂₂ = 2.88113263565 + 12.9270707357i and β̂₃ = 2.81221819020 + 26.4431630191i. For the error, we use the absolute value of the relative error, which is denoted by E_r. We plot log₂(E_r) against log₂(h) in Fig. 2, with h = 2/k being the step size. It is clear that our method converges and the numerical results suggest that our method has an O(h³) accuracy.

Fig. 2 Results from Example. Left: the propagation constants: ‘*’ stands for β obtained from the four-term WKB approximations; ‘o’ stands for results of β with k = 8; ‘+’ stands for results of β with k = 16. Right: relative errors for different step sizes: E_r stands for the absolute value of the relative error; h = 2/k is the step size; ‘o’ represents log₂(E_r) for β̂₁; ‘+’ represents log₂(E_r) for β̂₂; ‘*’ represents log₂(E_r) for β̂₃.

Download Full Size | PDF

The obtained propagation constants are all complex and the imaginary part of β is positive, which means that the optical wave decays in the transverse direction. Even our example is a lossless waveguide, there are leaky modes with complex propagation constants, which are not real modes of a lossless medium in the strict sense. Mathematically, these leaky modes appear in the approximation of the radiation modes, and they are useful in the computation of the lightwave propagation when the eigenmode expansion method is applied. The physical meaning of these leaky modes is discussed in detail in [17].

4. Conclusions

For a varying refractive-index waveguide, the refractive-index function n₀(x) is approximated by a piecewise polynomial of degree two. Since the eigenfunctions of the resulting Sturm-Liouville problem of the Helmholtz operator can be expressed analytically by the Kummer functions in each layer, analytical approximate dispersion equations are established for both TE and TM cases. Obviously, the approximate dispersion equations will approach to the exact ones as the number k of subintervals tends to infinity, or equivalently the step size h tends to zero, with a third order convergence. In the numerical example, we find the roots of the dispersion equation (i.e. propagation constants) by Müller’s method, where three different asymptotic solutions play the roles of initial values. Numerical simulations demonstrate that the iteration converges fast and high-precision values for propagation constants may be obtained only if a suitable root-finding method is adopted and good initial values are given. Further research results will appear in the future study.

5. Appendix A. The confluent hypergeometric functions

Let’s consider the second order equation with the quadratic potential

y^{″} (x) + (a x^{2} + b x + c) y (x) = 0,

where a ≠ 0. By the change of variables

t = x + \frac{b}{2 a}

, we have

\frac{d y^{2}}{d t^{2}} + (a t^{2} + C) y = 0,

where

C = c - \frac{b^{2}}{4 a}

. Further change of variables

z = {(- a)}^{1 / 2} t^{2} and w (z) exp (- z / 2) = y (t)

leads to

z \frac{d^{2} w}{d z^{2}} + (\frac{1}{2} - z) \frac{d w}{d z} - (\frac{1}{4} - \frac{C}{4 {(- a)}^{1 / 2}}) w (z) = 0,

which is in the form of the confluent hypergeometric equation

z \frac{d^{2} w}{d z^{2}} + (B - z) \frac{d w}{d z} - A w (z) = 0

with

A = \frac{1}{4} - \frac{C}{4 {(- a)}^{1 / 2}}

and

B = \frac{1}{2}

. Two linearly independent solutions are given by the Kummer functions M(A, B, z) and U(A, B, z), see [31]. However, U(A, B, z) has a branch point z = 0, and a fundamental pair of solutions that is numerically satisfactory near the origin is given by M(A, B, z) and z^1−BM(A + B − 1, 2 − B, z); see §13.2 of [31]. Hence, two linearly independent solutions of Eq. (23) are given by

{\begin{cases} α (x) = exp [- \frac{1}{2} z (x)] M (A, \frac{1}{2}, z (x)) \\ β (x) = {[z (x)]}^{\frac{1}{2}} exp [- \frac{1}{2} z (x)] M (A + \frac{1}{2}, \frac{3}{2}, z (x)), \end{cases}

where

z (x) = {(- a)}^{1 / 2} {(x + \frac{b}{2 a})}^{2}, A = \frac{1}{4} - \frac{C}{4 {(- a)}^{1 / 2}} = \frac{1}{4} + \frac{b^{2} - 4 a c}{16 a {(- a)}^{1 / 2}};

see Eq. (13.2.26) in [31]. Recall the differential formula for the Kummer functions:

\frac{d}{d z} M (a, b, z) = \frac{a}{b} M (a + 1, b + 1, z);

see Eq. (13.3.15) in [31]. The derivatives of α(x) and β(x) can be expressed as

{\begin{cases} α^{'} (x) = - \frac{1}{2} z^{'} (x) α (x) + 2 A z^{'} (x) exp [- \frac{z (x)}{2}] M (A + 1, \frac{3}{2}, z (x)) \\ β^{'} (x) = \frac{1 - z (x)}{2 z (x)} z^{'} (x) β (x) + \frac{2 A - 1}{3} {[z (x)]}^{\frac{1}{2}} z^{'} (x) exp [- \frac{z (x)}{2}] M (A + \frac{3}{2}, \frac{5}{2}, z (x)), \end{cases}

where z′(x) denotes the derivative of z(x).

6. Appendix B. Special case: a = 0

When a = 0, and b ≠ 0, Eq. (23) becomes

y^{″} (x) + (b x + c) y (x) = 0,

and two linearly independent solutions are given by the Airy functions

α (x) = Ai (w (x)), β (x) = Bi (w (x))

with

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

Therefore,

α^{'} (x) = - b^{1 / 3} {Ai}^{'} (w (x)), β^{'} (x) = - b^{1 / 3} {Bi}^{'} (w (x))

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

airy (z) = Ai (z), airy (1, z) = {Ai}^{'} (z), airy (2, z) = Bi (z), airy (3, z) = {Bi}^{'} (z) .

When a = b = 0, Eq. (23) becomes

y^{″} (x) + c y (x) = 0,

which has two linearly indecent solutions

α (x) = sin (c^{1 / 2} x), β (x) = cos (c^{1 / 2} x) .

7. Appendix C. Numerical evaluation of n_x and n_xx

In the TM case we have used a transformation to reduce Eq. (12) to a simpler form given in Eq. (15), and the function s²(x) plays the same role as n²(x) in the TE case. To this end, one has to calculate the values of n_x and n_xx at the nodes when interpolating s²(x); c.f. Eq. (18). Note that the piece-wise interpolation to s²(x) has the accuracy O(h³). In order to evaluate n_x and n_xx numerically with the same accuracy O(h³), we shall use high-order finite difference formulas. Five point formulas for the first and second order derivatives will be used.

Let x₀ denote one of the nodes. Then, the central five point formulas are

n^{'} (x_{0}) = \frac{1}{6 h} [n (x_{0} - h) - 8 n (x_{0} - 0.5 h) + 8 n (x_{0} + 0.5 h) - n (x_{0} + h)] + O (h^{4}),

n^{″} (x_{0}) = \frac{1}{3 h^{2}} [n (x_{0} - h) - 16 n (x_{0} - 0.5 h) - 30 f (x_{0}) + 16 n (x_{0} + 0.5 h) - n (x_{0} + h)] + O (h^{4}) .

In case that x₀ is one of the first two or last two nodes, the central difference can not apply. We shall make use of the forward or backward difference formulas. For example, when x₀ = 0 or x₀ = h/2, we shall use

n^{'} (x_{0}) = \frac{1}{6 h} [- 25 n (x_{0}) + 48 n (x_{0} + 0.5 h) - 36 n (x_{0} + h) + 16 n (x_{0} + 1.5 h) - 3 n (x + 2 h)] + O (h^{4}),

n^{″} (x_{0}) = \frac{1}{3 h^{2}} [35 n (x_{0}) - 104 n (x_{0} + 0.5 h) + 114 n (x_{0} + h) - 56 n (x_{0} + 1.5 h) + 11 n (x + 2 h)] + O (h^{3}) .

When x₀ = d or x₀ = d − h/2, one just needs to replace h by −h in the last two formulas:

n^{'} (x_{0}) = \frac{1}{6 h} [- 25 n (x_{0}) + 48 n (x_{0} - 0.5 h) - 36 n (x_{0} - h) + 16 n (x_{0} - 1.5 h) - 3 n (x - 2 h)] + O (h^{4}),

n^{″} (x_{0}) = \frac{1}{3 h^{2}} [35 n (x_{0}) - 104 n (x_{0} - 0.5 h) + 114 n (x_{0} - h) - 56 n (x_{0} - 1.5 h) + 11 n (x - 2 h)] + O (h^{3}) .

Acknowledgments

Both of the authors would like to thank the reviewers for their valuable comments that help improve the manuscript. The work of Y. Li was partially supported by a General Research Fund from Hong Kong Research Grants Council (No. 201513) and the HKBU Strategic Development Fund. The work of J. Zhu was partially supported by a grant from the Natural Science Foundation of China (NSFC) (No. 11371319), a grant from the Key Project of the Major Research Plan of NSFC (No. 91130004), and a grant from the Zhejiang Provincial Natural Science Foundation of China (No. LY13A010002).

References and links

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, 906–915 (1985). [CrossRef]

2. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26, S113–S134 (1994). [CrossRef]

3. D. Marcuse, Theory of Dielectric Optical Waveguides (Academic Press, 1974).

4. E. F. Keuster and R. C. Pate, “Fundamental mode propagation on dielectric fibres of arbitrary cross-section,” Proc. IEE 127(1), 41–51 (1980).

5. A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).

6. J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992). [CrossRef]

7. R. Ye and D. Yevick, “Noniterative calculation of complex propagation constants in planar waveguides,” J. Opt. Soc. Am. A 18, 2819–2822 (2001). [CrossRef]

8. D. Stowell and J. Tausch, “Variational formulation for guided and leaky modes in multilayer dielectric waveguides,” Commun. Comput. Phys. 7, 564–579 (2010).

9. M. Koshiba, K. Hayata, and M. Suzuki, “Vectorial finite-element method without spurious solutions for dielectric waveguide problems,” Electron. Lett. 20, 409–410 (1984). [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, 1545–1549 (1993). [CrossRef]

11. S. Selleri, L. Vincetti, A. Cucinotta, and M. Zoboli, “Complex FEM modal solver of optical waveguides with PML boundary conditions,” Opt. Quantum Electron. 33, 359–371 (2001). [CrossRef]

12. 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, 457–465 (2005). [CrossRef]

13. 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, 56–66 (2008). [CrossRef]

14. 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, 2077–2087 (2012). [CrossRef]

15. 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, 14683–14697 (2013). [CrossRef] [PubMed]

16. A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17, 311–321 (1985). [CrossRef]

17. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photon. 1, 58–106 (2009). [CrossRef]

18. J. Zhu and Y. Lu, “Leaky modes of slab waveguides-asymptotic solutions,” J. Lightwave Technol. 24, 1619–1623 (2006). [CrossRef]

19. J. Zhu and J. Zheng, “Exact dispersion equation of transverse electric leaky modes for nonhomogeneous optical waveguides,” J. Opt. Soc. Am. B 32, 92–100 (2015). [CrossRef]

20. 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, 230–238 (2005). [CrossRef]

21. C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14, 11631–11652 (2006). [CrossRef] [PubMed]

22. D. Song and Y. Lu, “Pesudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32, 1624–1630 (2014). [CrossRef]

23. J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for mirochip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48, 555–562 (2008). [CrossRef]

24. S. Khorasani and K. Mehrany, “Differential transfer-matrix method for solution of one-dimensional linear non-homogeneous optical structures,” J. Opt. Soc. Am. B 20, 91–96 (2003). [CrossRef]

25. M. Eghlidi, K. Mehrany, and B. Rashidian, “Modified differential transfer-matrix method for solution of one-dimensional linear inhomogeneous optical structures,” J. Opt. Soc. Am. B 22, 1521–1528 (2005). [CrossRef]

26. J. Zhu and Z. Shen, “Dispersion relation of leaky modes in nonhomogeneous waveguides and its applications,” J. Lightwave Technol. 29, 3230–3236 (2011). [CrossRef]

27. J. Zhu and J. Zheng, “Nonlinear equation of the modes in circular slab waveguides and its application,” Appl. Opt. 52, 8013–8023 (2013). [CrossRef]

28. S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10, 55–68 (1973). [CrossRef]

29. R. L. Burden and J. D. Faires, Numerical Analysis (Brooks/Cole, 2010).

30. P. J. Frey and P. George, Mesh Generation (Wiley, 2010).

31. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

m	β for k = 2	β for k = 4	β for k = 8	β for k = 16
1	10.417366 + 0.867180i	10.417863 + 0.868322i	10.418553 + 0.869157i	10.418414 + 0.869183i
2	6.636388 + 2.699835i	8.738858 + 1.501234i	8.739917 + 1.501640i	8.739609 + 1.501915i
3	6.636388 + 2.699835i	6.638866 + 2.700423i	8.739917 + 1.501640i	8.739609 + 1.501915i
4	6.636388 + 2.699835i	6.638866 + 2.700423i	4.518018 + 5.021335i	6.639969 + 2.701231i
5	3.453215 + 7.948715i	4.516944 + 5.021703i	4.518018 + 5.021335i	4.518790 + 5.021659i
6	3.453215 + 7.948715i	3.054062 + 10.577451i	4.518018 + 5.021335i	3.456118 + 7.946856i
7	2.879598 + 12.928845i	3.054062 + 10.577451i	3.054503 + 10.577007i	2.881129 + 12.927022i
8	2.879598 + 12.928845i	2.799926 + 15.096018i	2.880952 + 12.927243i	2.881129 + 12.927022i
9	2.799230 + 15.097164i	2.799926 + 15.096018i	2.800366 + 15.095566i	2.800478 + 15.095377i
10	2.764141 + 17.143554i	2.799926 + 15.096018i	2.765045 + 17.142103i	2.765119 + 17.141943i
11	2.753670 + 19.104243i	2.754224 + 19.103265i	2.754421 + 19.102815i	2.754470 + 19.102683i
12	2.757737 + 21.001673i	2.758165 + 21.000766i	2.758320 + 21.000393i	2.758362 + 21.000289i
13	2.770714 + 22.851069i	2.771071 + 22.850221i	2.771221 + 22.849879i	2.771265 + 22.849785i
14	2.789390 + 24.662743i	2.789722 + 24.661989i	2.789848 + 24.661669i	2.789887 + 24.661581i
15	2.811707 + 26.444330i	2.812027 + 26.443573i	2.812140 + 26.443289i	2.812174 + 26.443209i
16	2.836397 + 28.201307i	2.836659 + 28.200590i	2.836790 + 28.200331i	2.836818 + 28.200259i
17	2.862568 + 29.937992i	2.862794 + 29.937313i	2.862932 + 29.937035i	2.862959 + 29.936972i
18	2.889654 + 31.657597i	2.889870 + 31.656975i	2.889978 + 31.656687i	2.890008 + 31.656635i
19	2.917229 + 33.362775i	2.917441 + 33.362155i	2.917517 + 33.361884i	2.917557 + 33.361839i
20	2.945021 + 35.055551i	2.945200 + 35.054956i	2.945267 + 35.054713i	2.945320 + 35.054661i

Efficient approximations of dispersion relations in optical waveguides with varying refractive-index profiles

Abstract

1. Introduction

2. Mathematical model and treatment

2.1. TE case

2.2. TM case

3. Numerical example

4. Conclusions

5. Appendix A. The confluent hypergeometric functions

6. Appendix B. Special case: a = 0

7. Appendix C. Numerical evaluation of n_x and n_xx

Acknowledgments

References and links

Cited By

Figures (2)

Tables (1)

Equations (63)

Optics Express

Abstract

1. Introduction

2. Mathematical model and treatment

2.1. TE case

2.2. TM case

3. Numerical example

4. Conclusions

5. Appendix A. The confluent hypergeometric functions

6. Appendix B. Special case: a = 0

7. Appendix C. Numerical evaluation of nx and nxx

Acknowledgments

References and links

Cited By

Figures (2)

Tables (1)

Equations (63)

Optics Express

7. Appendix C. Numerical evaluation of n_x and n_xx