This paper discusses a numerical method for computing the electromagnetic modes supported by multilayer planar optical waveguides constructed from lossy or active media, having in general a diagonal permittivity tensor. The method solves the dispersion equations in the complex plane via the Cauchy integration method. It is applicable to lossless, lossy and active waveguides, and to AntiResonant Reflecting Optical Waveguides (ARROW’s). Analytical derivatives for the dispersion equations are derived and presented for what is believed to be the first time, and a new algorithm that significantly reduces the time required to compute the derivatives is given. This has a double impact: improved accuracy and reduced computation time compared to the standard approach. A different integration contour, which is suitable for leaky modes is also presented. Comparisons are made with results found in the literature; excellent agreement is noted for all comparisons made.
© Optical Society of America
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) , , 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 , Newton’s method and the one-dimensional scan method  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 ,  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 , . 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 (µ=µ0), is shown in Figure 1. The refractive index tensor, , of the i-th layer can be in general complex, i.e., ñjji =njji-jkjji where njji and kjji 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.
The TMM is a well-known approach for generating the dispersion equations governing the modes supported by a multilayer planar optical waveguide . 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  for source-free, time-harmonic fields in anisotropic media are:
where ε0 is the free space permittivity, ω is the angular frequency, and ε̿r is a tensor of relative permittivity having the form:
For a TE mode (Ex, Ez, Hy =0) propagating in the +ẑ direction in the i-th layer, (xi ≤x≤x i+1), the non-zero electric and magnetic field components are:
where x̂, ŷ, ẑ are the unit vectors in the x, y, z directions, respectively, =K 0(β-jα) is the complex propagation constant with β and α the normalized phase and attenuation constants, respectively, and k0 =ω/c=2π/λ0 , c is the speed of light in free space and λ0 is the free space wavelength. From equations (1) and (3), we can derive:
where . 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:
where xi 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=xi ), can be expressed as a function of the field and its derivative within that layer:
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 (Eys, dEys/dx) and at the boundary at the bottom of the cover layer (Eyc, dEyc/dx) are related via the matrix product:
The Mi are the transfer matrices for all of the r layers of thickness di . For guided modes in an open structure, the tangential fields in the substrate and cover must be exponentially decaying:
where , 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:
The zeros of Equation (12) correspond to the complex propagation constants .
2.3 Transfer Matrix and the Dispersion Equation for the TM Modes
For the TM modes (Hx, Hz, Ey =0) the following equations hold:
Following a derivation similar to the TE case the transfer matrices of each layer are determined:
where and the dispersion equation is:
3. Cauchy Integration Method
3.1 Summary of the Method
The Cauchy integration method  is based on the argument principle  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:
where Nz is the number of zeros and Np is number of poles inside the region enclosed by the contour C. If there are no poles in the region enclosed by C, then S0 is the number of zeros, and from the residue theorem we have:
where zi, i=1, 2, …, S0 are the roots of f(z) inside C and Sm is the sum of with m=1, 2, …, S0 . Equation (18) leads to a system of equations that can be used to evaluate the coefficients of a polynomial p(z) of degree S 0, which has the same roots z 1, … as the function f(z) inside C. The approximation polynomial p(z) can be written as:
with . The coefficients Ck are given via Newton’s recursive formula:
The polynomial p(z) can be solved by standard techniques such as Laguerre’s method . 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 C1 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 C1 . The integral contour C2 is selected for leaky modes. These contours do not enclose any of the poles.
(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 , which states that the first derivative of a function f(z) at z=z0 is given by:
where D is any closed path which encloses the point z0 and f(z) is analytic inside and on D.
Reference  gives an algorithm for solving Equation (21) numerically. The algorithm is based on the parametric variable transformation z=z0 +Rej2πx where R is the radius of a circle centered at z0 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:
The computation of the numerical derivative according to Equation (22) is the most intensive computational task in the standard CIM . When z0 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  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 .
To solve Equations (17) and (18), integration is performed numerically by applying an adaptive integration method  in a manner similar to reference . This approach is summarized here for completeness. The quadrature integration of a function F(z) is represented by the summation:
where Wi are the weights and a≤zi ≤b are the integration nodes. A 7-point Gaussian rule G7f and 15-point Kronrod  rule K15f 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 K15f since the Kronrod rule is more accurate than the Gaussian rule. If the local error estimate |K15f-G7f|/|K15f| 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 f,(zi )/f(zi ) for m=1, …, S 0 can be incorporated such that f(zi) and f’(zi ) are evaluated only once for the specific node zi for all the summations S 1, …, .
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  with the initial guess being the roots of p(z).
As discussed in reference , the region enclosed by the contour of interest C ought to be subdivided into smaller regions such that S0 ≤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 =k0ũ, , and . Then from Equation (9) we derive for each layer i:
The derivative of the transfer matrix describing the multilayer structure for i=1, …, r is then:
The derivatives of s and c with respect to ũ are: and , thus from Equation (12) we obtain the derivative of the dispersion equation as:
4.2 TM Modes
Similarly, for the TM modes where . Then and , and from Equation (15) we derive for each layer i:
From Equation (16), the derivative of the dispersion equation is:
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 1 and are computed. For the second layer (i=2), M 2 and are computed, then and M=M 1 M 2 are accumulated. For the third layer (i=3), M 3 and are computed, then and M=M 1 M 2 M 3 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 .
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  was also implemented and its computation times are included for comparison. The parameters of the structure analyzed were taken from reference :
where ns =2.177, nc =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 TE0 mode and β=2.17930 for the TE1 mode . 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.
Figure 3 shows the execution time required to compute the effective index of the TE0 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.
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.
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 s and c must be chosen. If we assume s= + and c= + , for β<ns (β<nc ), we should choose <0 and >0 ( <0 and >0). The integral contour is selected as C2 , shown in Fig. 2(a) for ns > nc . 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 , as shown.
5.4 ARROW Waveguides
ARROW (AntiResonant Reflecting Optical Waveguide) waveguides ,  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 , where the effective index of the first 6 TE and TM ARROW modes found using our CIM are also given.
Figure 4 shows the spatial distribution of the Hy 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 .
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 . The structure of interest is defined in Table V and the first 4 ARROW modes for the isotropic (AR=1) and anisotropic cases (AR=nxxi/nzzi=1.03) are given. Our results for the TE modes are in complete agreement with those reported in  for the case AR=1.
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.
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]