Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media

Chengkun Chen; Pierre Berini; Dazeng Feng; Stoyan Tanev; Velko P. Tzolov

doi:10.1364/OE.7.000260

1. Introduction

The notions of guided and leaky modes are fundamental concepts in optical waveguiding theory. A knowledge of mode propagation characteristics is essential to the design of numerous guided-wave optoelectronic devices, passive and active components such as semiconductor lasers, electro-absorption and electro-optic modulators, switches, photodetectors, filters and couplers, to name but a few. Numerical methods that can efficiently and accurately model planar optical waveguides are thus of obvious importance since they are used as a basic tool in the design process.

The Transfer Matrix Method (TMM) [1], [2], as one of the primary tools for multilayer planar optical waveguide analysis, can generate the dispersion equation of the TE and TM modes supported by such structures in a straightforward manner. The waveguides can consist of any combination of lossless, lossy (dielectric, semiconductor, metallic) and active (including uniaxially anisotropic quantum wells) layers. Solving the dispersion equation yields the propagation constant and the electromagnetic field distribution for a mode. For the modes of lossy or active waveguides, and for leaky modes, the mode propagation constants, which are the roots of the dispersion equation, are complex numbers.

Traditional numerical zero-search algorithms such as the downhill method [3], Newton’s method and the one-dimensional scan method [4] cannot give the number of propagating modes supported by the waveguides, and need an initial guess value close to the actual root. Therefore these methods are not efficient and reliable, especially for a general-purpose mode solver. There is a rigorous mathematical technique [5], [6] which is capable of finding the zeros or poles of any analytic function in the complex plane. This technique is based on integration in the complex plane and can be used to solve the dispersion equation of a multilayer planar optical waveguide [7], [8]. The method is often called the Cauchy Integration Method (CIM) or the Argument Principle Method (APM).

The novel features of the TMM and CIM that we present in this paper include: the derivation and use of analytical derivatives of the dispersion equations in anisotropic media, an algorithm for the rapid calculation of the derivatives, and a different integration contour for locating leaky modes. Using analytical derivatives in conjunction with the algorithm suggested improves the accuracy of the method and greatly reduces the CPU time required to find modes. The new integration contour improves the numerical efficiency of the integration calculation. The method is applicable to lossless, lossy and active waveguides, and to AntiResonant Reflecting Optical Waveguides (ARROW’s). The method can find leaky and guided modes. Anisotropic dielectrics described by a diagonal permittivity tensor are handled by the formulation.

2. Transfer Matrix Method

The transfer matrix method provides a straightforward formulation of the multilayer planar optical waveguide problem. A multilayer nonmagnetic anisotropic slab waveguide structure (µ=µ₀), is shown in Figure 1. The refractive index tensor, ${\bar{\bar{n}}}_{i}$ , of the i-th layer can be in general complex, i.e., ñ_jji =n_jji-jk_jji where n_jji and k_jji are the refractive index and extinction coefficient along the jj direction of the i-th layer, and i=1,…, r is the layer number between the substrate and cover.

Figure 1. Structure of the multilayer planar optical waveguide

Download Full Size | PDF

The TMM is a well-known approach for generating the dispersion equations governing the modes supported by a multilayer planar optical waveguide [1]. The technique is only summarized here and the equations needed to obtain the derivatives analytically are given.

2.1 Maxwell’s Equations and TE Field Solutions

Maxwell’s curl equations [9] for source-free, time-harmonic fields in anisotropic media are:

\nabla \times \overset{⇀}{E} = - j ω μ_{0} \overset{⇀}{H}

\nabla \times \overset{⇀}{H} = - j ω ε_{0} {\bar{\bar{ε}}}_{r} \overset{⇀}{E}

where ε₀ is the free space permittivity, ω is the angular frequency, and ε̿_r is a tensor of relative permittivity having the form:

{\bar{\bar{ε}}}_{r} = [\begin{matrix} {\tilde{ε}}_{r, xx} \\ {\tilde{ε}}_{r, yy} \\ {\tilde{ε}}_{r, zz} \end{matrix}] = [\begin{matrix} {\tilde{n}}_{xx}^{2} \\ {\tilde{n}}_{yy}^{2} \\ {\tilde{n}}_{zz}^{2} \end{matrix}]

For a TE mode (E_x, E_z, H_y =0) propagating in the +ẑ direction in the i-th layer, (x_i ≤x≤x _i+1), the non-zero electric and magnetic field components are:

{\overset{⇀}{E}}_{i} = \hat{y} E_{yi} (x) \exp (j ω t - j \tilde{γ} z)

{\overset{⇀}{H}}_{i} = [\hat{x} H_{xi} (x) + \hat{z} H_{zi} (x)] \exp (j ω t - j \tilde{γ} z)

where x̂, ŷ, ẑ are the unit vectors in the x, y, z directions, respectively, $\tilde{γ}$ =K ₀(β-jα) is the complex propagation constant with β and α the normalized phase and attenuation constants, respectively, and k₀ =ω/c=2π/λ₀ , c is the speed of light in free space and λ₀ is the free space wavelength. From equations (1) and (3), we can derive:

\frac{d^{2} E_{yi} (x)}{dx} + {\tilde{κ}}_{i}^{2} E_{yi} (x) = 0

and

H_{zi} (x) = \frac{j}{ω μ_{0}} \frac{{dE}_{yi} (x)}{dx}

H_{xi} (x) = - \frac{\tilde{γ}}{ω μ_{0}} E_{yi} (x)

where ${\tilde{κ}}_{i} = \pm \sqrt{k_{0}^{2} {\tilde{n}}_{yyi}^{2} - {\tilde{γ}}^{2}}$ . The tangential electric field and its derivative must satisfy equation (4) and are related to the magnetic fields through equation (5) within the i-th layer. Solutions to equation (4) are written as:

E_{yi} (x) = \cos [{\tilde{κ}}_{i} (x - x_{i})] E_{yi} (x_{i}) + \frac{1}{{\tilde{κ}}_{i}} \sin [{\tilde{κ}}_{i} (x - x_{i})] \frac{{dE}_{yi} (x_{i})}{dx}

\frac{{dE}_{yi} (x_{i})}{dx} = {- \tilde{κ}}_{i} \sin [{\tilde{κ}}_{i} (x - x_{i})] E_{yi} (x_{i}) + \cos [{\tilde{κ}}_{i} (x - x_{i})] \frac{{dE}_{yi} (x_{i})}{dx}

where x_i defines the boundary between the i-th and (i+1)-th layer. Equation (6) implies that any sign for the square root of k̃_i is acceptable.

2.2 Transfer Matrix and the Dispersion Equation for the TE Modes

Using Equation (6), the tangential electric field and its derivative at the bottom of the i-th layer, (x=x_i ), can be expressed as a function of the field and its derivative within that layer:

(\begin{matrix} E_{yi} (x_{i}) \\ \frac{{dE}_{yi} (x_{i})}{dx} \end{matrix}) = (\begin{matrix} \cos [{\tilde{κ}}_{i} (x - x_{i})] & - \frac{1}{{\tilde{κ}}_{i}} \sin [{\tilde{κ}}_{i} (x - x_{i})] \\ {\tilde{κ}}_{i} \sin [{\tilde{κ}}_{i} (x - x_{i})] & \cos [{\tilde{κ}}_{i} (x - x_{i})] \end{matrix}) (\begin{matrix} E_{yi} (x) \\ \frac{{dE}_{yi} (x)}{dx} \end{matrix})

Imposing the continuity of tangential fields at any layer interface in the multilayer structure, the fields tangential to the boundary at the top of the substrate layer (E_ys, dE_ys/dx) and at the boundary at the bottom of the cover layer (E_yc, dE_yc/dx) are related via the matrix product:

(\begin{matrix} E_{ys} \\ \frac{{dE}_{ys}}{dx} \end{matrix}) = \prod_{i = 1}^{r} M_{i} (\begin{matrix} E_{yc} \\ \frac{{dE}_{yc}}{dx} \end{matrix})

where

M_{i} = (\begin{matrix} \cos ({\tilde{κ}}_{i} d_{i}) & - \frac{1}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) \\ {\tilde{κ}}_{i} \sin ({\tilde{κ}}_{i} d_{i}) & \cos ({\tilde{κ}}_{i} d_{i}) \end{matrix}) for i = 1, 2, \dots, r

The M_i are the transfer matrices for all of the r layers of thickness d_i . For guided modes in an open structure, the tangential fields in the substrate and cover must be exponentially decaying:

\begin{matrix} E_{ys} (x) = A_{s} \exp ({\tilde{γ}}_{s} x) \\ \frac{{dE}_{ys} (x)}{dx} = {\tilde{γ}}_{s} A_{s} \exp ({\tilde{γ}}_{s} x) \end{matrix}} for x \leq 0

\begin{matrix} E_{yc} (x) = B_{c} \exp [- {\tilde{γ}}_{c} (x - x_{r + 1})] \\ \frac{{dE}_{yc} (x)}{dx} = - {\tilde{γ}}_{c} B_{c} \exp [- {\tilde{γ}}_{c} (x - x_{r + 1})] \end{matrix}} for x \geq x_{r + 1}

where ${\tilde{γ}}_{s} = \pm \sqrt{{\tilde{γ}}^{2} - k_{0}^{2} {\tilde{n}}_{yys}^{2}}$ , ${\tilde{γ}}_{c} = \pm \sqrt{{\tilde{γ}}^{2} - k_{0}^{2} {\tilde{n}}_{yyc}^{2}}$ and ñ_yys and ñ_yyc are the substrate and cover complex refractive indices along the y direction, respectively. Equations (10) and (11) in conjunction with (8) yield the dispersion equation for the TE modes:

F (\tilde{γ}) = {\tilde{γ}}_{s} m_{11} + {\tilde{γ}}_{c} m_{22} - m_{21} - {\tilde{γ}}_{s} {\tilde{γ}}_{c} m_{12} = 0

The zeros of Equation (12) correspond to the complex propagation constants $\tilde{γ}$ .

2.3 Transfer Matrix and the Dispersion Equation for the TM Modes

For the TM modes (H_x, H_z, E_y =0) the following equations hold:

\frac{d^{2} H_{yi} (x)}{dx} + {\tilde{κ}}_{i}^{2} H_{yi} (x) = 0

E_{zi} (x) = - \frac{j}{ω ε_{0} {\tilde{n}}_{zzi}^{2}} \frac{{dH}_{yi} (x)}{dx}

E_{xi} (x) = - \frac{\tilde{γ}}{ω ε_{0} {\tilde{n}}_{xxi}^{2}} H_{yi} (x)

Following a derivation similar to the TE case the transfer matrices of each layer are determined:

M_{i} = (\begin{matrix} \cos ({\tilde{κ}}_{i} d_{i}) & - \frac{{\tilde{n}}_{zzi}^{2}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) \\ \frac{{\tilde{κ}}_{i}}{{\tilde{n}}_{zzi}^{2}} \sin ({\tilde{κ}}_{i} d_{i}) & \cos ({\tilde{κ}}_{i} d_{i}) \end{matrix}) for i = 1, 2, \dots, r

where ${\tilde{κ}}_{i} = \pm \frac{{\tilde{n}}_{zzi}}{{\tilde{n}}_{xxi}} \sqrt{k_{0}^{2} {\tilde{n}}_{xxi}^{2} - {\tilde{γ}}^{2}}$ and the dispersion equation is:

F (\tilde{γ}) = \frac{{\tilde{γ}}_{s}}{{\tilde{n}}_{zzs}^{2}} m_{11} + \frac{{\tilde{γ}}_{c}}{{\tilde{n}}_{zzc}^{2}} m_{22} - m_{21} - \frac{{\tilde{γ}}_{s} {\tilde{γ}}_{c}}{{\tilde{n}}_{zzs}^{2} {\tilde{n}}_{zzc}^{2}} m_{12} = 0

3. Cauchy Integration Method

3.1 Summary of the Method

The Cauchy integration method [5] is based on the argument principle [10] and the residue theorem of complex analysis. If a function f(z) is analytic and does not go to zero over a closed integral contour, then the argument principle is of the form:

S_{0} = \frac{1}{j 2 π} \oint_{c} \frac{f^{’} (z)}{f (z)} dz = N_{z} - N_{p}

where N_z is the number of zeros and N_p is number of poles inside the region enclosed by the contour C. If there are no poles in the region enclosed by C, then S₀ is the number of zeros, and from the residue theorem we have:

S_{m} = \frac{1}{j 2 π} \oint_{c} z^{m} \frac{f^{’} (z)}{f (z)} dz = \sum_{i = 1}^{s_{0}} z_{i}^{m} for m = 1, 2, \dots, S_{0}

where z_i, i=1, 2, …, S₀ are the roots of f(z) inside C and S_m is the sum of $z_{i}^{m}$ with m=1, 2, …, S₀ . Equation (18) leads to a system of equations that can be used to evaluate the coefficients of a polynomial p(z) of degree S ₀, which has the same roots z ₁, … $z_{S_{0}}$ as the function f(z) inside C. The approximation polynomial p(z) can be written as:

p (z) = \prod_{i = 1}^{s_{0}} (z - z_{i}) = \sum_{k = 0}^{s_{0}} C_{k} z^{k}

with $C_{s_{0}} = 1$ . The coefficients C_k are given via Newton’s recursive formula:

C_{k} = \frac{1}{(k - S_{0})} \sum_{j = 1}^{s_{0} - k} S_{j} C_{k + j} for k = S_{0} - 1, \dots, 0

The polynomial p(z) can be solved by standard techniques such as Laguerre’s method [11]. The problem of finding the zeros of an arbitrary function f(z) is thus transformed to the simpler problem of finding the zeros of the polynomial p(z), for which a variety of reliable and efficient numerical methods exist.

3.2 Numerical Implementation

The CIM can be used to solve the dispersion equation of the multilayer planar optical waveguide in a straightforward manner. The poles are first identified at β=ñ_c and β=ñ_s as shown in Figure 2, parts (a) and (b). The rectangular integral contour C₁ is selected for the guided modes. The area between the minimum and the maximum real parts of all the refractive indices is enclosed by the contour C₁ . The integral contour C₂ is selected for leaky modes. These contours do not enclose any of the poles.

Figure 2. The integral contours in the complex plane. The oe-7-8-260-i001 correspond to the poles, ñ _max is the complex refractive index with the maximum real part and is enclosed by C1, for TE modes in anisotropic media: ñ _max=ñ _yymax, ñ_s =ñ_yys and ñ_c =ñ_yyc . For TM modes in anisotropic media: ñ _max=ñ _xxmax, ñ_s =ñ_xxs and ñ_c =ñ_xxc .

Download Full Size | PDF

(a) ñ_c <ñ_s , (b) ñ_c =ñ_s .

To solve Equations (17) and (18) the derivative of the dispersion equation is needed. If the derivative cannot be obtained analytically then it can be estimated numerically via Cauchy’s theorem [10], which states that the first derivative of a function f(z) at z=z₀ is given by:

f^{’} (z_{0}) = \frac{1}{j 2 π} \oint_{D} \frac{f (z)}{{(z - z_{0})}^{2}} dz

where D is any closed path which encloses the point z₀ and f(z) is analytic inside and on D.

Reference [7] gives an algorithm for solving Equation (21) numerically. The algorithm is based on the parametric variable transformation z=z₀ +Re^j2πx where R is the radius of a circle centered at z₀ defining the contour over which the function f(z) is evaluated. The m-point trapezoidal integration rule is then applied to Equation (21), with z transformed as described, yielding the following numerical approximation to the derivative:

R f^{’} (z_{0}) \approx R {[f^{’} (z_{0})]}^{m} = \frac{1}{m} \sum_{l = 1}^{m} \frac{f (z_{0} + {Re}^{j \frac{2 π}{m} l})}{e^{j \frac{2 π}{m} l}}

The computation of the numerical derivative according to Equation (22) is the most intensive computational task in the standard CIM [7]. When z₀ is very close to the zero of the function, the numerical derivative converges very slowly and inaccurately.

In multilayer planar optical waveguides, the derivative of the dispersion equations can be obtained analytically. Reference [8] gives an analytical derivative but in exponential form and for isotropic media. We have obtained the exact analytical expressions for the derivatives of the transfer matrices and of the dispersion equations for the multilayer waveguide problem of interest in this paper. The expressions are given in the next section and an algorithm to efficiently construct them is described. The derivatives computed in this manner are then used to help solve Equations (17) and (18), thus providing a gain in accuracy and reduction in execution time compared to the standard CIM approach [7].

To solve Equations (17) and (18), integration is performed numerically by applying an adaptive integration method [12] in a manner similar to reference [7]. This approach is summarized here for completeness. The quadrature integration of a function F(z) is represented by the summation:

\int_{a}^{b} F (z) = \sum_{i = 1}^{K} W_{i} F (z_{i})

where W_i are the weights and a≤z_i ≤b are the integration nodes. A 7-point Gaussian rule G₇f and 15-point Kronrod [13] rule K₁₅f are used to estimate integration along the straight line connecting the end-points of the interval [a, b] in the complex plane. The local estimate is taken as K₁₅f since the Kronrod rule is more accurate than the Gaussian rule. If the local error estimate |K₁₅f-G₇f|/|K₁₅f| is less than a specified tolerance, the local estimate is accepted as the integration result for the specific interval. Otherwise the interval is bisected and the same procedure is applied to the two new subintervals. In the same integration routine the calculation of the products $z_{i}^{m}$ f,(z_i )/f(z_i ) for m=1, …, S ₀ can be incorporated such that f(z_i) and f’(z_i ) are evaluated only once for the specific node z_i for all the summations S ₁, …, $S_{S_{0}}$ .

The polynomial p(z) can be solved by Laguerre’s method. The roots of polynomial p(z) and f(z) do not coincide exactly due to the integration errors introduced in Equation (23). After the roots of the polynomial p(z) are obtained, a further refinement must be performed to find the roots of f(z) by applying Muller’s Method [11] with the initial guess being the roots of p(z).

As discussed in reference [7], the region enclosed by the contour of interest C ought to be subdivided into smaller regions such that S₀ ≤4. This ensures that the degree and the coefficients of the polynomial p(z) remain small, thus helping to reduce numerical errors when locating the roots. Considering all sub-regions enclosed by the original contour leads to all roots of interest.

4. Analytical Derivatives of the Transfer Matrices and Dispersion Equations

4.1. TE Modes

For the TE modes, let $\tilde{γ}$ =k₀ũ, ${\tilde{κ}}_{i} = \pm k_{0} \sqrt{{\tilde{n}}_{yyi}^{2} - {\tilde{u}}^{2}}$ , $\frac{{d \tilde{κ}}_{i}}{d \tilde{u}} = - \frac{k_{0}^{2}}{{\tilde{κ}}_{i}} \tilde{u}$ and $\frac{d}{d \tilde{u}} (\frac{1}{{\tilde{κ}}_{i}}) = \frac{k_{0}^{2}}{{\tilde{κ}}_{i}^{3}} \tilde{u}$ . Then from Equation (9) we derive for each layer i:

\frac{{dM}_{i}}{d \tilde{u}} = \tilde{u} (\begin{matrix} \frac{k_{0}^{2} d_{i}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) & \frac{k_{0}^{2} d_{i}}{{\tilde{κ}}_{i}^{2}} \cos ({\tilde{κ}}_{i} d_{i}) - \frac{k_{0}}{{\tilde{κ}}_{i}^{3}} \sin ({\tilde{κ}}_{i} d_{i}) \\ - \frac{k_{0}^{2}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) - k_{0}^{2} d_{i} \cos ({\tilde{κ}}_{i} d_{i}) & \frac{k_{0}^{2} d_{i}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) \end{matrix}) for i = 1, \dots, r

The derivative of the transfer matrix describing the multilayer structure for i=1, …, r is then:

\frac{dM}{d \tilde{u}} = (\begin{matrix} \frac{{dm}_{11}}{d \tilde{u}} & \frac{{dm}_{12}}{d \tilde{u}} \\ \frac{{dm}_{21}}{d \tilde{u}} & \frac{{dm}_{22}}{d \tilde{u}} \end{matrix}) = \sum_{j = 1}^{r} (\frac{{dM}_{j}}{d \tilde{u}} \cdot \prod_{\underset{i \neq j}{i = 1}}^{r} M_{i})

The derivatives of $\tilde{γ}$ _s and $\tilde{γ}$ _c with respect to ũ are: $\frac{{d \tilde{γ}}_{s}}{d \tilde{u}} = \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{s}}$ and $\frac{{d \tilde{γ}}_{c}}{d \tilde{u}} = \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{c}}$ , thus from Equation (12) we obtain the derivative of the dispersion equation as:

\frac{dF (\tilde{u})}{d \tilde{u}} = \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{s}} m_{11} + \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{c}} m_{22} + {\tilde{γ}}_{s} \frac{{dm}_{11}}{d \tilde{u}} + {\tilde{γ}}_{c} \frac{{dm}_{22}}{d \tilde{u}} - \frac{{dm}_{21}}{d \tilde{u}} - \frac{k_{0}^{2} \tilde{u} {\tilde{γ}}_{c}}{{\tilde{γ}}_{s}} m_{12} - \frac{k_{0}^{2} \tilde{u} {\tilde{γ}}_{s}}{{\tilde{γ}}_{c}} m_{12} - {\tilde{γ}}_{s} {\tilde{γ}}_{c} \frac{{dm}_{12}}{d \tilde{u}}

4.2 TM Modes

Similarly, for the TM modes ${\tilde{κ}}_{i} = \pm \frac{{\tilde{n}}_{zzi}}{{\tilde{n}}_{xxi}} k_{0} \sqrt{{\tilde{n}}_{xxi}^{2} - {\tilde{u}}^{2}} = \pm k \sqrt{{\tilde{n}}_{xxi}^{2} - {\tilde{u}}^{2}}$ where $k = \frac{{\tilde{n}}_{zzi}}{{\tilde{n}}_{xxi}} k_{0}$ . Then $\frac{{d \tilde{κ}}_{i}}{d \tilde{u}} = - \frac{k^{2}}{{\tilde{κ}}_{i}} \tilde{u}$ and $\frac{d}{d \tilde{u}} (\frac{1}{{\tilde{κ}}_{i}}) = \frac{k^{2}}{{\tilde{κ}}_{i}^{3}} \tilde{u}$ , and from Equation (15) we derive for each layer i:

\frac{{dM}_{i}}{d \tilde{u}} = \tilde{u} (\begin{matrix} \frac{k^{2} d_{i}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) & \frac{k^{2} {\tilde{n}}_{zzi}^{2} d_{i}}{{\tilde{κ}}_{i}^{2}} \cos ({\tilde{κ}}_{i} d_{i}) - \frac{k^{2} {\tilde{n}}_{zzi}^{2}}{{\tilde{κ}}_{i}^{3}} \sin ({\tilde{κ}}_{i} d_{i}) \\ - \frac{k^{2}}{{\tilde{κ}}_{i} {\tilde{n}}_{zzi}^{2}} \sin ({\tilde{κ}}_{i} d_{i}) - \frac{k^{2} d_{i}}{{\tilde{n}}_{zzi}^{2}} \cos ({\tilde{κ}}_{i} d_{i}) & \frac{k^{2} d_{i}}{{\tilde{κ}}_{i}} \sin ({\tilde{κ}}_{i} d_{i}) \end{matrix}) for i = 1, \dots, r

From Equation (16), the derivative of the dispersion equation is:

\frac{dF (\tilde{u})}{d \tilde{u}} = \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{s} {\tilde{n}}_{xxs}^{2}} m_{11} + \frac{k_{0}^{2} \tilde{u}}{{\tilde{γ}}_{c} {\tilde{n}}_{xxc}^{2}} m_{22} + \frac{{\tilde{γ}}_{s}}{{\tilde{n}}_{zzs}^{2}} \frac{{dm}_{11}}{d \tilde{u}} + \frac{{\tilde{γ}}_{c}}{{\tilde{n}}_{zzc}^{2}} \frac{{dm}_{22}}{d \tilde{u}} - \frac{{dm}_{21}}{d \tilde{u}} - \frac{{\tilde{γ}}_{s} {\tilde{γ}}_{c}}{{\tilde{n}}_{zzs}^{2} {\tilde{n}}_{zzc}^{2}} \frac{{dm}_{12}}{d \tilde{u}}

4.3 Algorithm for Computing the Derivative of the Transfer Matrices

The transfer matrices, and their derivatives as expressed by Equation (25), can be computed in an efficient manner. The idea is to accumulate the derivative of the transfer matrix as the latter is being constructed. Consider accumulation over the first 3 layers of a waveguide structure. For the first layer (i=1), M ₁ and $\frac{{dM}_{1}}{d \tilde{u}}$ are computed. For the second layer (i=2), M ₂ and $\frac{{dM}_{2}}{d \tilde{u}}$ are computed, then $\frac{dM}{d \tilde{u}} = \frac{d (M_{1} M_{2})}{d \tilde{u}} = M_{1} \frac{{dM}_{2}}{d \tilde{u}} + \frac{{dM}_{1}}{d \tilde{u}} M_{2}$ and M=M ₁ M ₂ are accumulated. For the third layer (i=3), M ₃ and $\frac{{dM}_{3}}{d \tilde{u}}$ are computed, then $\frac{dM}{d \tilde{u}} = \frac{d (M_{1} M_{2} M_{3})}{d \tilde{u}} = \frac{d (M_{1} M_{2})}{d \tilde{u}} M_{3} + (M_{1} M_{2}) \frac{{dM}_{3}}{d \tilde{u}}$ and M=M ₁ M ₂ M ₃ are accumulated. The accumulation process continues until all layers in the structure have been handled. Clearly, one matrix product is required per layer to accumulate the transfer matrix and two products are required to accumulate its derivative. This algorithm thus accumulates the transfer matrix for r layers in r matrix products and the derivative of the transfer matrix in 2r matrix products. The CIM based on our analytical derivatives and algorithm is thus an order r method, which is a substantial improvement over a similar approach of order r! discussed in [14].

5. Numerical Results and Discussion

5.1 Guided Modes in a Lossless Waveguide

An inhomogeneous waveguide having an index of refraction that follows an exponential distribution was analyzed in order to compare the computational efficiency of the CIM using numerical and analytical derivatives. The ADR method [14] was also implemented and its computation times are included for comparison. The parameters of the structure analyzed were taken from reference [14]:

n (x) = \sqrt{n_{s}^{2} + 2 Δ n_{s} e^{(x ⁄ α)}} for x < 0

n (x) = n_{c} for x > 0

where n_s =2.177, n_c =1.0, Δ=0.043, α=0.931 µm and the thickness of the inhomogeneous layer is set to d=4 µm. The normalized phase constant (or effective index) of the two TE modes supported by this structure can be obtained analytically. They are β=2.19075 for the TE₀ mode and β=2.17930 for the TE₁ mode [14]. The effective indices of both modes, computed using our approach, are shown in Table I. The computed values are seen to converge to the analytical values as the number of layers used to approximate the continuous profile is increased. Agreement with the analytical values is excellent for a large number of layers.

TABLE I. Effective Index of Guided Modes Supported by an Inhomogeneous Waveguide.

View Table | View all tables in this article

Figure 3 shows the execution time required to compute the effective index of the TE₀ mode supported by the structure as a function of the number of layers used to approximate Equation (29a). The numerical derivative was computed to a relative accuracy of 10^-3 and the tolerance in assessing the local error estimate when evaluating Equation (23) was set to 10^-5 for all three methods. From Figure 3, it is apparent that the CIM based on the analytical derivative is order r in matrix operations, like the CIM based on the numerical derivative and the ADR method. Using an analytical derivative is more accurate and clearly much more computationally efficient than using a numerical one. The computational cost of the CIM with the analytical derivative is only slightly higher than the ADR method, which is less robust and more complex.

Figure 3. Comparison of the normalized execution times required for the ADR and the CIM based on a numerical derivative (ND) and an analytical derivative (AD).

Download Full Size | PDF

5.2 Guided Modes in Quantum Well Active Waveguides

Quantum well (QW) layers used in active photonic devices are in fact anisotropic. Their optical response is different for polarizations parallel and perpendicular to the direction of growth. The anisotropy is uniaxial and appears even though quantum wells are made up from cubic semiconductor materials. Modeling uniaxial media is essential if accurate waveguide propagation constants are desired. We applied our method to a novel InP-based QW laser structure as defined in Table II in order to illustrate its ability to handle anisotropy and gain. The second layer is the active anisotropic quantum well layer. The results are shown in Table II and gain is observed for the first TE and TM modes.

TABLE II. Effective Index of Guided Modes Supported by a 6-layer Quantum Well Active Waveguide.

View Table | View all tables in this article

5.3 Leaky Modes in Lossless Waveguides

The CIM can be used to find leaky modes. The leaky mode propagation constants satisfy the same dispersion equation as the guided modes, but the appropriate sign of $\tilde{γ}$ _s and $\tilde{γ}$ _c must be chosen. If we assume $\tilde{γ}$ _s= $γ_{s}^{r}$ + ${j γ}_{s}^{i}$ and $\tilde{γ}$ _c= $γ_{c}^{r}$ + ${j γ}_{c}^{i}$ , for β<n_s (β<n_c ), we should choose $γ_{s}^{r}$ <0 and $γ_{s}^{i}$ >0 ( $γ_{c}^{r}$ <0 and $γ_{c}^{i}$ >0). The integral contour is selected as C₂ , shown in Fig. 2(a) for n_s > n_c . The reason for selecting the top part of the integral contour along the β axis is that the imaginary part of the leaky mode propagation constants is negative. The integral contour above the β axis will lead to large integration errors and generate roots for the polynomial p(z) far away from the roots of the dispersion equation. For the TE modes supported by the waveguide structure described in Table III, our results are complete agreement with those reported in [1], as shown.

TABLE III. Effective Index of Leaky Modes Supported by a 4-layer Lossless Waveguide.

View Table | View all tables in this article

5.4 ARROW Waveguides

ARROW (AntiResonant Reflecting Optical Waveguide) waveguides [15], [16] are based on Fabry-Perot reflection instead of total internal reflection. The modes supported by such structures can be particularly difficult to find due to the proximity of the corresponding zeros in the complex plane. An example structure is described in Table IV [17], where the effective index of the first 6 TE and TM ARROW modes found using our CIM are also given.

TABLE IV. Effective Index of ARROW Modes Supported by a 9-layer ARROW Waveguide

View Table | View all tables in this article

Figure 4 shows the spatial distribution of the H_y field component related to the first two TM ARROW modes, computed using our method. It is indeed clear from this figure that the structure can be used as a directional coupler [17].

Figure 4. (a) Field distribution for the symmetric TM mode: TM₁. (b) Field distribution for the first anti-symmetric TM mode: TM₂.

Download Full Size | PDF

5.5 Anisotropic ARROW Waveguides

Many anisotropic materials have desirable properties at optical wavelengths, including low losses, and large electro-optic or photo-elastic effects. In addition, anisotropy may be introduced during material growth or device processing. It is thus important for a mutilayer waveguide modeling tool to be able to handle effectively such materials.

As the last example, we have applied our method to obtain the propagation constant of modes supported by an anisotropic ARROW waveguide [18]. The structure of interest is defined in Table V and the first 4 ARROW modes for the isotropic (AR=1) and anisotropic cases (AR=n_xxi/n_zzi=1.03) are given. Our results for the TE modes are in complete agreement with those reported in [19] for the case AR=1.

TABLE V. Effective Index of ARROW Modes Supported by a 3-layer Anisotropic ARROW Waveguide

View Table | View all tables in this article

6. Conclusion

A numerical method based on integration in the complex plane has been developed to characterize the modes supported by multilayer planar optical waveguides constructed from lossy or active anisotropic media. The analytical derivative of the transfer matrices and dispersion equations were obtained and presented for the first time. An algorithm for efficiently computing the derivative was given and a different integral contour to be used for finding leaky modes was presented. The method reported has several important advantages over other methods for finding modes of planar waveguides. An advantage is the ability to find the total number of modes in a region of interest in the complex plane. Another advantage is that the method can be used to characterize lossless, lossy, active and ARROW waveguides in anisotropic media. The method can generate guided modes as well as leaky modes. The combination of high accuracy, low computation time and broad range of applicability makes the method very attractive for integration into commercial computer-aided design and modeling tools for integrated optics.

7. Acknowledgements

The authors gratefully acknowledge Valery I. Tolstikhin, for providing the second example structure (an active anisotropic waveguide) investigated in this work.

References and links

1. J. Chilwell and I. Hodgkinson, “Thin-films field-transfer matrix theory of planar multilayer waveguides and reflection from prism-loaded waveguide,” J. Opt. Soc. Am. A 1, 742–753, (1984) [CrossRef]

2. L. M. Walpita, “Solutions for planar optical waveguide equations by selecting zero elements in a characteristic matrix,” J. Opt. Soc. Am. A 2, 595–602, (1985) [CrossRef]

3. K. H. Schlereth and M. Tacke, “The complex propagation constant of multilayer waveguides: An algorithm for a personal computer,” IEEE J. Quantum Electron. , 26, 627–630, (1990) [CrossRef]

4. L. Sun and E. Marhic, “Numerical study of attenuation in multilayer infrared waveguides by the circle-chain convergence method”, J. Opt. Soc. Am. B 8, 478–483, (1991) [CrossRef]

5. L. M. Delves and J. N. Lyness, “A numerical method for locating the zeros of an analytic function,” Math. Comp. , 21, 543–560, (1967) [CrossRef]

6. L. C. Botten and M. S. Craig, “Complex zeros of analytic functions”, Comput. Phys. Commun. , 29, 245–259, (1983) [CrossRef]

7. E. Anemogiannis and E. N. Glytsis, “Multilayer waveguides: efficient numerical analysis of general structures”, J. Lightwave Tech. , 10, 1344–1351, (1992) [CrossRef]

8. R. E. Smith, S. N. Houde-Walter, and G. W. Forbes, “Mode determination for planar waveguides using the four-sheeted dispersion relation,” IEEE J. Quantum Electron. , 28, 1520–1526, (1992) [CrossRef]

9. Hermann A. Haus, “Waves and Fields in Optoelectronics,” (New Jersey, Prentice-Hall Inc., 1984). Ch. 11

10. J. W. Brown and R. V. Churchill, “Complex Variables and Applications,” (Sixth Edition, New York: McGraw-Hill, 1996)

11. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Numerical Recipes in C,” (Second Edition, Cambridge, 1994)

12. J. R. Rice, “Numerical Methods, Software, and Analysis,” (IMSL Reference Edition. New York: McGraw-Hill, 1983)

13. A. S. Kronrod, “Nodes and Weights of Quadrature Formulas,” (NewYork: Consultants Bureau, 1965)

14. E. Anemogiannis, E. N. Glytsis, and T. K. Gaylord, “Efficient solution of eigenvalue equations of optical waveguiding structures”, J. Lightwave Tech. , 12, 2080–2084, (1994) [CrossRef]

15. T. Baba and Y. Kokubun, “Dispersion and radiation loss characteristics of antiresonsnt reflecting optical waveguides-Numerical Results and Analytical Expressions”, IEEE J. Quantum Electron. , 28, 1689–1700, (1992) [CrossRef]

16. W. Huang, R. M. Shubiar, A. Nathan, and Y. L. Chow, “The modal characteristics of ARROW structures”, J. Lightwave Tech. , 10, 1015–1022. (1992) [CrossRef]

17. J. Deng and Y. Huang,“A novel hybrid coupler based on antiresonant reflecting optical waveguides,” J. Lightwave Tech. , 16, 1062–1069, (1998) [CrossRef]

18. B. Ray and G W. Hanson, “Some effects of anisotropy on planar antiresonant reflecting optical waveguides”, J. Lightwave Tech. , 14, 202–208, (1996) [CrossRef]

19. E. Anemogiannis, E. N. Glytsis, and T. K. Gaylord, “Determination of guided and leaky modes in lossless and lossy planar multiplayer optical waveguides: reflection pole method and wavevector density method”, J. Lightwave Tech. , 17, 929–941, (1999) [CrossRef]

Number of Layers	β - TE₀	β - TE₁
5	2.191567664	2.179485876
10	2.190948646	2.179376481
20	2.190805357	2.179325841
40	2.190775872	2.179315352
60	2.190770683	2.179313499
80	2.190768887	2.179312857
100	2.190768059	2.179312561
200	2.190766958	2.179312167

6-layer Anisotropic Active Waveguide: ñ_s =3.13575-j6.2264×10^-5, ñ ₁=3.16404-j1.0005×10^-5, ñ _xx2=ñ _yy2=3.393856+j0.0069093, ñ _zz2=3.395891+j0.0032438, ñ ₃=3.16987-j1.3991×10^-4, ñ ₄=3.38838-j2.7498×10^-4, ñ ₅=3.17034-j2.829×10^-4, ñ ₆=3.46930-j0.083828, ñ_c =0.59-j12.63, d₁=1.5µm, d₂=0.2µm, d₃=0.2µm, d₄=0.01µm, d₅=1.3µm, d₆=0.2µm, λ₀=1.55µm
Mode	β	α×10^-3
TE₀	3.211912271	-2.295751586
TE₁	3.146335751	1.833386157
TE₂	3.137997320	2.519718051
TM₀	3.204193956	-1.750073292
TM₁	3.144266723	0.1557238150

4-layer Lossless Waveguide: n_s=1.5, n₁=1.66, n₂=1.60, n₃=1.53, n₄=1.66, n_c=1.0, d₁=d₂=d₃=d₄=0.5µm, λ₀=0.6328µm
Present Method			Reference [1]
Mode	β	α	β	α
TE₄	1.461856641	0.007155871	1.46186	0.00716
TE₅	1.382489223	0.018165877	1.38250	0.01817
TE₆	1.281364436	0.035877392	1.28136	0.03588
TE₇	1.142314462	0.052876075	1.14231	0.05288
TE₈	1.003037019	0.070770941	1.00304	0.07077
TM₄	1.451534978	0.011923599
TM₅	1.370664375	0.030142063
TM₆	1.273737061	0.056791773
TM₇	1.157312853	0.087578491
TM₈	1.036950265	0.103078083

9-layer ARROW Waveguide: n_s=3.5, n₁=1.46, n₂=1.50, n₃=1.46, n₄=1.50, n₅= 1.46, n₆=1.50, n₇=1.46, n₈=1.50, n₉=1.46, n_c=1.0, d₁=2.0µm, d₂=0.448µm, d₃= 4.0µm, d₄=0.448µm, d₅=2.0µm, d₆=0.448µm, d₇=4.0µm, d₈=0.448µm, d₉=2.0µm, λ₀=0.6328µm
Mode	β	α×10^-4
TE₁	1.457920191	0.007106242
TE₂	1.457791244	0.009053396
TE₃	1.453780369	0.114698816
TE₄	1.453045406	0.420121480
TE₅	1.451864807	0.693651857
TE₆	1.450269491	0.732515869
TM₁	1.457925423	0.045880488
TM₂	1.457782773	0.057163274
TM₃	1.453795449	0.645756672
TM₄	1.452928430	2.555862981
TM₅	1.451781628	4.567101184
TM₆	1.450247659	4.357488809

3-layer Anisotropic ARROW Waveguide: n_s=3.85, n_zz1=1.46, n_zz2=2.3, n_zz3=1.46, n_c =1, d₁=3.15λ₀, d₂=0.142λ₀, d₃=6.3λ₀, λ₀=0.6328µm, n_xxi=n_yyi, AR=n_xxi/n_zzi, i=1,2,3
AR=1					AR=1.03
Present Method			Reference [19]		Present Method
Mode	β	α×10^-3	β	α×10^-3	β	α×10^-3
TE₁	1.457941265	0.000054189	1.45794	0.000054189	1.501798936	0.000050179
TE₂	1.451919174	0.052870681	1.45192	0.052871	1.495945499	0.053815143
TE₃	1.451174055	0.192035341	1.45117	0.19203	1.495255344	0.184243873
TE₄	1.441371363	0.004374469	1.44137	0.0043745	1.485698165	0.004051178
TM₁	1.457890856	0.002450742			1.501625054	0.002544521
TM₂	1.451754691	0.553891897			1.495287895	0.576101022
TM₃	1.451304282	1.151033285			1.494855078	1.189339701
TM₄	1.440916633	0.319061714			1.484121307	0.197863211

Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media

Abstract

1. Introduction

2. Transfer Matrix Method

2.1 Maxwell’s Equations and TE Field Solutions

2.2 Transfer Matrix and the Dispersion Equation for the TE Modes

2.3 Transfer Matrix and the Dispersion Equation for the TM Modes

3. Cauchy Integration Method

3.1 Summary of the Method

3.2 Numerical Implementation

4. Analytical Derivatives of the Transfer Matrices and Dispersion Equations

4.1. TE Modes

4.2 TM Modes

4.3 Algorithm for Computing the Derivative of the Transfer Matrices

5. Numerical Results and Discussion

5.1 Guided Modes in a Lossless Waveguide

5.2 Guided Modes in Quantum Well Active Waveguides

5.3 Leaky Modes in Lossless Waveguides

5.4 ARROW Waveguides

5.5 Anisotropic ARROW Waveguides

6. Conclusion

References and links

Cited By

Figures (4)

Tables (5)

Equations (37)

Optics Express