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

Optics Express, Vol. 7, Issue 8, pp. 260-272 (2000)

http://dx.doi.org/10.1364/OE.7.000260

Acrobat PDF (164 KB)

### Abstract

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

## 1. Introduction

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]

## 2. Transfer Matrix Method

_{0}), is shown in Figure 1. The refractive index tensor,

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

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.1 Maxwell’s Equations and TE Field Solutions

*ε*

_{0}is the free space permittivity,

*ω*is the angular frequency, and

*ε̿*

_{r}is a tensor of relative permittivity having the form:

*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:

*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

*k*

_{0}=

*ω/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:

*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

*i*-th layer, (

*x*=

*x*

_{i}), can be expressed as a function of the field and its derivative within that layer:

*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:

*ñ*

_{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:

*.*γ ˜

### 2.3 Transfer Matrix and the Dispersion Equation for the TM Modes

*H*

_{x}

*, H*

_{z}

*, E*

_{y}=0) the following equations hold:

## 3. Cauchy Integration Method

### 3.1 Summary of the Method

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]

*f(z)*is analytic and does not go to zero over a closed integral contour, then the argument principle is of the form:

*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*

_{0}is the number of zeros, and from the residue theorem we have:

*z*

_{i}

*, i*=1, 2, …,

*S*

_{0}are the roots of

*f(z)*inside

*C*and

*S*

_{m}is the sum of

*m*=1, 2, …,

*S*

_{0}. 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}, …

*f(z)*inside

*C*. The approximation polynomial

*p(z)*can be written as:

*C*

_{k}are given via Newton’s recursive formula:

*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

*β*=

*ñ*

_{c}and

*β*=

*ñ*

_{s}as shown in Figure 2, parts (a) and (b). The rectangular integral contour

*C*

_{1}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*

_{1}. The integral contour

*C*

_{2}is selected for leaky modes. These contours do not enclose any of the poles.

*ñ*

_{c}

*<ñ*

_{s}, (b)

*ñ*

_{c}=

*ñ*

_{s}.

*f(z)*at

*z*=

*z*

_{0}is given by:

*D*is any closed path which encloses the point

*z*

_{0}and

*f(z)*is analytic inside and on

*D*.

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

*z*=

*z*

_{0}+

*Re*

^{j2πx}where

*R*is the radius of a circle centered at

*z*

_{0}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:

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

*z*

_{0}is very close to the zero of the function, the numerical derivative converges very slowly and inaccurately.

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]

**10**, 1344–1351, (1992) [CrossRef]

**10**, 1344–1351, (1992) [CrossRef]

*F(z)*is represented by the summation:

*W*

_{i}are the weights and

*a*≤

*z*

_{i}≤

*b*are the integration nodes. A 7-point Gaussian rule

*G*

_{7}

*f*and 15-point Kronrod [13] rule

*K*

_{15}

*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

_{15}f since the Kronrod rule is more accurate than the Gaussian rule. If the local error estimate |

*K*

_{15}

*f*-

*G*

_{7}

*f*|/|

*K*

_{15}

*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

*f*,(

*z*

_{i})/

*f*(

*z*

_{i}) for

*m*=1, …,

*S*

_{0}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*

_{1}, …,

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

**10**, 1344–1351, (1992) [CrossRef]

*C*ought to be subdivided into smaller regions such that

*S*

_{0}≤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

*=*γ ˜

*k*

_{0}

*ũ*,

*i*=1, …,

*r*is then:

*and*γ ˜

_{s}*with respect to*γ ˜

_{c}*ũ*are:

### 4.2 TM Modes

*i*:

### 4.3 Algorithm for Computing the Derivative of the Transfer Matrices

*i*=1),

*M*

_{1}and

*i*=2),

*M*

_{2}and

*M*=

*M*

_{1}

*M*

_{2}are accumulated. For the third layer (

*i*=3),

*M*

_{3}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 [14].

## 5. Numerical Results and Discussion

### 5.1 Guided Modes in a Lossless Waveguide

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]

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]

*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

_{0}mode and β=2.17930 for the TE

_{1}mode [14

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]

_{0}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

### 5.3 Leaky Modes in Lossless Waveguides

*and*γ ˜

_{s}*must be chosen. If we assume*γ ˜

_{c}*=*γ ˜

_{s}*=*γ ˜

_{c}*β*<

*n*

_{s}(

*β*<

*n*

_{c}), we should choose

*C*

_{2}, 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

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]

### 5.4 ARROW Waveguides

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]

### 5.5 Anisotropic ARROW Waveguides

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]

*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

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]

*AR*=1.

## 6. Conclusion

## 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 |

2. | L. M. Walpita, “Solutions for planar optical waveguide equations by selecting zero elements in a characteristic matrix,” J. Opt. Soc. Am. A |

3. | K. H. Schlereth and M. Tacke, “The complex propagation constant of multilayer waveguides: An algorithm for a personal computer,” IEEE J. Quantum Electron. , |

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 |

5. | L. M. Delves and J. N. Lyness, “A numerical method for locating the zeros of an analytic function,” Math. Comp. , |

6. | L. C. Botten and M. S. Craig, “Complex zeros of analytic functions”, Comput. Phys. Commun. , |

7. | E. Anemogiannis and E. N. Glytsis, “Multilayer waveguides: efficient numerical analysis of general structures”, J. Lightwave Tech. , |

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

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

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

16. | W. Huang, R. M. Shubiar, A. Nathan, and Y. L. Chow, “The modal characteristics of ARROW structures”, J. Lightwave Tech. , |

17. | J. Deng and Y. Huang,“A novel hybrid coupler based on antiresonant reflecting optical waveguides,” J. Lightwave Tech. , |

18. | B. Ray and G W. Hanson, “Some effects of anisotropy on planar antiresonant reflecting optical waveguides”, J. Lightwave Tech. , |

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

**OCIS Codes**

(130.2790) Integrated optics : Guided waves

(230.7390) Optical devices : Waveguides, planar

**ToC Category:**

Research Papers

**History**

Original Manuscript: September 5, 2000

Published: October 9, 2000

**Citation**

Chengkun Chen, Pierre Berini, Dazeng Feng, Stoyan Tanev, and Velko Tzolov, "Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media," Opt. Express **7**, 260-272 (2000)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-7-8-260

Sort: Journal | Reset

### References

- 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]
- 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]
- 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]
- 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]
- 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]
- L. C. Botten and M. S. Craig, "Complex zeros of analytic functions," Comput. Phys. Commun. 29, 245- 259, (1983) [CrossRef]
- E. Anemogiannis, and E. N. Glytsis, "Multilayer waveguides: efficient numerical analysis of general structures," J. Lightwave Tech. 10, 1344-1351, (1992) [CrossRef]
- 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]
- Hermann A. Haus, Waves and Fields in Optoelectronics, (New Jersey, Prentice-Hall Inc., 1984). Ch. 11
- J. W. Brown and R. V. Churchill, Complex Variables and Applications, (Sixth Edition, New York: McGraw-Hill, 1996)
- W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in C, (Second Edition, Cambridge, 1994).
- J. R. Rice, Numerical Methods, Software, and Analysis, (IMSL Reference Edition. New York: McGraw-Hill, 1983)
- A. S. Kronrod, Nodes and Weights of Quadrature Formulas, (NewYork: Consultants Bureau, 1965)
- 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]
- 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]
- 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]
- J. Deng and Y. Huang,"A novel hybrid coupler based on antiresonant reflecting optical waveguides," J. Lightwave Tech. 16, 1062-1069, (1998) [CrossRef]
- B. Ray and G W. Hanson, "Some effects of anisotropy on planar antiresonant reflecting optical waveguides," J. Lightwave Tech. 14, 202-208, (1996) [CrossRef]
- 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]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

OSA is a member of CrossRef.