method proposed by McKeeâMitchell is applied to full-vectorial paraxial wave equations. ... developed an ADI process for parabolic equations wit...

0 downloads 0 Views 197KB Size

JOURNAL OF LIGHTWAVE TECHNOLOGY, VOL. 16, NO. 12, DECEMBER 1998

Full-Vectorial Beam-Propagation Method Based on the McKee–Mitchell Scheme with Improved Finite-Difference Formulas Junji Yamauchi, Member, IEEE, Member, OSA, Gen Takahashi, and Hisamatsu Nakano, Fellow, IEEE

Abstract—The alternating-direction implicit method proposed by McKee–Mitchell is applied to full-vectorial paraxial wave equations. The high computational efficiency of the present method is demonstrated in comparison with an iterative solver. Novel finite-difference formulas that take into account discontinuities of the fields are proposed and employed to ensure second-order accuracy. Calculations regarding the effective index of rib waveguides show that the present results remarkably agree with values obtained from the modal transverse resonance method. Index Terms—Finite-difference methods, optical beam propagation, optical waveguides.

I. INTRODUCTION

C

ONSIDERABLE effort has been directed to improving a full-vectorial beam-propagation method (VBPM) from various viewpoints. The VBPM’s developed so far have often caused problems in terms of stability and/or accuracy. This is due to the existence of mixed derivatives in coupled wave equations. Since vectorial treatment is absolutely necessary to simulate polarization-dependent and coupling devices, accurate and stable numerical methods should be developed. Pioneer work on a VBPM has been made by Huang and Xu [1] using finite-difference (FD) techniques with an iterative solver. Since the iterative solver is somewhat time-consuming, Mansour et al. [2] introduced the alternating-direction implicit (ADI) technique into the VBPM. Their formulation was based on the Peaceman–Rachford scheme. Although the Peaceman–Rachford splitting is effective in scalar [3] and semivectorial calculations [4], [5], the application to the fullvectorial case often causes instability in numerical results. Yevick et al. [6] introduced the ADI technique using lines at 45 to the principal axis. This technique, the so-called Douglas–Gunn splitting, yields stable results but is extremely complicated. On the other hand, McKee–Mitchell [7] have developed an ADI process for parabolic equations with a mixed derivative, which is known to yield reasonable results with simpler splitting. However, no attempt has been made to apply the McKee–Mitchell scheme to a VBPM. One should also note that inaccuracy in previous VBPM’s stems from improper treatment of the mixed derivatives.

Since the field and its derivatives are often discontinuous, the conventional FD formulas severely deteriorate the numerical results. In the case of a semivectorial analysis, several improved FD formulas have been proposed, which attain considerable accuracy [8], [9]. This fact encourages us to develop improved FD formulas even for the mixed derivatives in the VBPM. In this paper, we propose a novel VBPM based on the McKee–Mitchell (MM) scheme [7], and demonstrate its better stability than the VBPM based on the well-known Peaceman–Rachford (PR) scheme, while maintaining high computational efficiency. Special attention is paid to evaluation of mixed derivatives as well as second derivatives. We propose new FD formulas taking into account discontinuities of refractive indexes. The new FD formulas ensure secondorder accuracy even in full-vectorial calculations. The present VBPM can also be used in imaginary-axis propagation. We assess the accuracy of the present VBPM in the eigenmode analysis in a rib waveguide as a classical benchmark test, since reliable data on propagation constants are available [10]. Calculation shows how appropriate evaluation of the mixed derivatives improves convergence behavior. Negligence of discontinuities of the first derivative in previous works results in poor convergence. Remarkable agreement not obtainable in previous works is observed with the results derived by the modal transverse resonance method [10], [11]. II. MCKEE–MITCHELL SCHEME Assuming that the refractive index varies slowly along , we solve the following paraxial vector wave equations [1]: (1) (2) where

Manuscript received February 25, 1998; revised June 30, 1998. The authors are with the College of Engineering, Hosei University, Koganei, Tokyo 184 Japan (e-mail: [email protected]). Publisher Item Identifier S 0733-8724(98)09120-8. 0733–8724/98$10.00 1998 IEEE

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

(3) (4) (5)

YAMAUCHI et al.: FULL-VECTORIAL BEAM PROPAGATION

2459

with , and , in which is the free-space wavenumber, is the index profile, is the reference refractive index. and Based on the MM scheme, which is known to yield reasonable results with simpler splitting than the Douglas–Gunn scheme [7], we derive FD equations from the coupled paraxial wave equations. Since (1) and (2) have a similar form, we only treat (1) in the following. After integrating (1) over an , we get interval of

We first express the fields is given by expansions.

and

using Taylor-series

(10) , the continuity relations at the interface With respect to and refer to the fields at should be satisfied. Let the infinitesimally right and left of the interface, respectively. is used, can be written as When

(6) (11) where the differential operators placed by difference operators Using the assumption of

, and are re, and , respectively. , we obtain

On the other hand, using

can be expressed as

(12) Since the relation between

and

[8], [9]

(7) The ADI process, derived by McKee–Mitchell, is in unsplit form

can be rewritten in the form (8) Now we introduce the Douglas–Rachford type splitting, so that we obtain (13) and . where The FD formula for the second derivative can be obtained by eliminating the first derivative from (10) and (13) [9] (9a) (14) where

(9b) Equation (9) involves the solution of only two tridiagonal sets of equations at each propagation step. This is in contrast to the Douglas–Gunn splitting, which requires the solution of four tridiagonal sets of equations at each propagation step. III. IMPROVED FD FORMULAS We treat three consecutive mesh points with a discontinuity and derive FD formulas when the between points and interface between different media is located midway between two mesh points.

, and , in which . The coefficient can be approximated by without significant effects. Equation (14), regarded as an extension of Stern’s formula [12], is a special case of (6) in our previous paper [9]. Note that recently Vassallo [13] also developed similar improved FD formulas for the second derivative. Consideration is now given to the treatment of the mixed derivative. No attempt has been made to evaluate the mixed derivatives accurately. Direct discretization of the mixed derivative results in poor convergence behavior, as will be shown in Fig. 5. To evaluate the mixed derivative correctly, we consider the first derivative, since the mixed derivative is composed of the combination of the first derivatives. As in the case of derivation of (14), we can derive the following

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

2460

JOURNAL OF LIGHTWAVE TECHNOLOGY, VOL. 16, NO. 12, DECEMBER 1998

Fig. 2. Comparison in computational effort between the McKee–Mitchell scheme and the iterative solver. Fig. 1. Overlap integral between the numerical and exact fields for a minor component of HE11 mode of a step-index fiber as a function of propagation distance z .

three-point FD formula for the first derivative by eliminating the second derivative from (10) and (13): (15) and . where Equations (14) and (15) are utilized to evaluate the first and second derivatives in (9). For example, (14) is applied to or (3), and to or (4) with . Equation (15) is and in [see (5)] with . used for may be evaluated by another FD formula, in which the interface conditions of electric flux density are satisfied. Following the similar procedure mentioned above, we obtain (16) Corresponding formulas similar to (14)–(16) are also derived and (see when the discontinuity lies between points the Appendix). IV. RESULTS We first compare the MM scheme with the PR scheme in propagation behavior of HE11 mode of a step-index fiber. Fig. 1 shows the overlap integral between the numerical and exact fields for the minor component as a function of propagation distance. The exact fields are chosen as input fields. The transverse mesh size is taken to be m. The total number of mesh points is . The configuration parameters are as m, and the refractive follows. The core radius is and indexes of the core and cladding are . A wavelength of m is used, so that the . This fiber is not realistic, normalized frequency is but is considered to investigate sensitivity of the schemes in a strongly guiding structure. It is found that the MM scheme than that in the allows a larger propagation step length PR scheme. Further calculation shows similar tendency in the major component. The PR scheme is efficient only for scalar [3] and semivectorial cases [4], [5].

Fig. 3. Rib waveguide geometry.

High computational efficiency of the MM scheme is demonstrated in Fig. 2, in which comparison with the VBPM using an iterative solver (Bi-CGSTAB) [14] is made. The computing time in the Bi-CGSTAB depends on a tolerance factor for convergence. The tolerance factor is chosen in such a way that the results obtained by the Bi-CGSTAB are comparable to those obtained by the MM scheme. In the MM scheme, the standard Thomas algorithm can be used, leading to less CPU time and memory. For a step-index fiber, we cannot correctly assess the accuracy of FD formulas due to discretization error of a circular core. Hence, we next consider a rib waveguide (shown in Fig. 3) used as a classical benchmark [10]. The configuration , , m, rib parameters are m, central rib height m, and lateral width varying from 0.1 to 0.9 m. The computational height m above the top of the rib, domain parameters are m below the guiding layer, and m (or m) aside from the rib lateral side. 5.5 m for The fundamental mode can be calculated by the imaginary distance procedure in which the coordinate in the propaga[15]. The effective index is tion direction is changed to evaluated by the growth in the field amplitude [16], as shown in (17) at the bottom of the next page. For evaluation of the , the reference index effective index defined by has to be a value close to the exact one. It should be noted, however, that the exact value is unknown in practice. In this paper, we adopt a new technique of iteratively renewing [17]. The effectiveness of this technique is shown in

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

YAMAUCHI et al.: FULL-VECTORIAL BEAM PROPAGATION

2461

Fig. 4. Convergence behavior of the normalized propagation constant B as a function of propagation .

Fig. 4, in which the normalized propagation constant is presented as a function of . The data m. The transverse are for the quasi-TE mode with and propagation step length are mesh size m and m, respectively. taken to be As the input field, we choose a step function as for otherwise

m and

Fig. 5. Convergence behavior of the normalized propagation constant B as a function of transverse mesh size .

1

TABLE I NORMALIZED PROPAGATION CONSTANT B OBTAINED WITH MTRM AND DEVIATIONS IN SEMIVECTORIAL AND FULL-VECTORIAL CASES

m

for the major component and zero for the minor component. to be or , converges to When we fix iteratively leads a different value. In contrast, renewing to monotone convergence to the same value regardless of an (it has been confirmed that converges initial value of to the exact value in a two-dimensional waveguide [17]). The convergence value is close to a value obtained with the modal transverse resonance method (MTRM) as will be discussed in detail. for the quasi-TE mode as The convergence behavior of is shown in Fig. 5, a function of transverse mesh size is where is chosen to be 0.5 m. In this calculation, (although a larger fixed to be 0.05 m regardless of is available when is large). For comparison, the results obtained using the combination of other FD formulas [1], [12] are also shown. It is revealed that the use of (14)–(16) achieves faster convergence. Although not illustrated, a similar tendency is also observed in the quasi-TM mode. The results regarding the semivectorial (sv) and full-vectorial (fv) cases with the improved FD formulas are found to have almost the same convergence behavior. values obtained from the present Table I tabulates the technique in both the quasi-TE and quasi-TM modes. The is taken to be 0.025 m. For reference, the data mesh size values for the semivectorial case are also presented. The 10 from are expressed with the deviations

the MTRM. The four-digit values of ’s are believed to be exact [10], [11]. Table I also presents extrapolated and indicates that the difference between values at the semivectorial and full-vectorial results is slight. It should be noted, however, that the data on the full-vectorial case ’s. The deviation for fv remarkably agree with is only within 1. It can be said that accurate evaluation of the mixed derivatives greatly contributes to improvement of accuracy in terms of agreement with the MTRM. Typical field distributions in the quasi-TE and quasi-TM modes are illustrated in Figs. 6 and 7, respectively. Singularities and discontinuities of the fields are clearly displayed. Comparison in deviations with other techniques is shown in Fig. 8. All the data are derived using an extrapolation

(17)

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

2462

JOURNAL OF LIGHTWAVE TECHNOLOGY, VOL. 16, NO. 12, DECEMBER 1998

(a) Fig. 6. Field distributions in the quasi-TE mode: (a)

(b)

Ex

and (b)

Ey .

(a) Fig. 7. Field distributions in the quasi-TM mode: (a)

(b)

Ey

and (b)

Ex .

technique and expressed as a function of . One of the data was obtained by Pregla using the method of lines (MOL) [10], [18]. Another data was obtained using the Yee’s mesh VBPM [19], [20]. In general, all the results relatively agree with the data of the MTRM. Strictly speaking, however, the MOL achieves better results in the quasi-TE mode than the quasi-TM mode, while the Yee’s mesh VBPM indicates opposite tendency. It is noteworthy that the present method is successful for both the quasi-TE and the quasi-TM mode. Fig. 9 shows another interesting comparison among three different discretization schemes: Stern type (present), Bierwirth type [21], and Yee type [20]. The rib waveguide configuration to be considered here is the one used in [21]. Note that in contrast to the Stern scheme, the discontinuity lines are just on mesh points in the Bierwirth scheme. Also note that in the Yee scheme, the components of and are interlaced within a unit cell. It is worth mentioning that the present scheme again shows the fastest convergence, although each scheme is second-order accurate. The reason why a Stern-type scheme shows faster convergence than a Bierwirth-type scheme is

probably due to the fact that the Stern-type scheme avoids evaluation of the fields at the corners of the rib waveguide. One should note that quite recently, Hadley succeeded in deriving a high-accuracy formula in the Bierwirth type [22]. It seems, however, that the high-accuracy formula is effective except for corner points. V. CONCLUSIONS We have presented an improved vectorial finitedifference BPM based on the McKee–Mitchell scheme. The present scheme is more stable than the conventional Peaceman–Rachford scheme and achieves high computational efficiency in comparison with an iterative solver. The mixed derivatives in the coupled wave equations are evaluated taking into account discontinuities of refractive indexes. The present BPM ensures second-order accuracy, provided that the discontinuity lies midway between two mesh points. The propagation constants of a rib waveguide for a classical benchmark are evaluated using the imaginary distance procedure in both quasi-TE and quasi-TM modes. The

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

YAMAUCHI et al.: FULL-VECTORIAL BEAM PROPAGATION

2463

(a)

(b)

Fig. 8. Comparison in deviations from MTRM as a function of t: (a) quasi-TE mode and (b) quasi-TM mode.

ACKNOWLEDGMENT The authors would like to thank Prof. W. P. Huang for sending literature regarding an iterative solver. Thanks are also due to Dr. W. W. Lui of NTT Opto-Electronics Laboratories for invaluable discussions on field singularities at the corners of a rib waveguide. The authors also acknowledge the help of N. Morohashi, who provided the data concerning the Yee’s mesh VBPM. REFERENCES

Fig. 9. Comparison in convergence behavior of the normalized propagation constant B among three discretization schemes.

obtained results remarkably agree with the values derived with the modal transverse resonance method. APPENDIX When a discontinuity lies between points and , the FD formulas corresponding to (14)–(16), respectively, are derived as follows: (18) , , in which

where , and

,

, and (19) where

and (20)

can be approximated by The coefficient significant effects.

without

[1] W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron., vol. 29, no. 10, pp. 2639–2649, 1993. [2] I. Mansour, A.-D. Capobianco, and C. Rosa, “Noniterative vector beam propagation method with a smoothing digital filter,” J. Lightwave Technol., vol. 14, no. 5, pp. 908–913, 1996. [3] J. Yamauchi, T. Ando, and H. Nakano, “Beam propagation analysis of optical fibers by alternating direction implicit method,” Electron. Lett., vol. 27, no. 18, pp. 1663–1664, 1991. [4] P.-L. Liu and B. J. Li, “Semivectorial beam-propagation method for analyzing polarized modes of rib waveguides,” J. Quantum Electronics, vol. 28, no. 4, pp. 778–782, 1992. [5] J. Yamauchi, J. Shibayama, O. Saito, O. Uchiyama, and H. Nakano, “Improved finite-difference beam-propagation method based on the generalized Douglas scheme and its application to semivectorial analysis,” J. Lightwave Technol., vol. 14, no. 10, pp. 2401–2406, 1996. [6] D. Yevick, J. Yu, W. Bardyszewski, and M. Glasner, “Stability issues in vector electric field propagation,” Photon. Tech. Lett., vol. 7, no. 6, pp. 658–660, 1995. [7] S. McKee and A. R. Mitchell, “Alternating direction methods for parabolic equations in two space dimensions with a mixed derivative,” Comput. J., vol. 13, no. 1, pp. 81–86, 1970. [8] C. Vassallo, “Improvement of finite difference methods for step-index optical waveguides,” Proc. Inst. Elect. Eng.,, vol. 139, no. 2, pt. J, pp. 137–142, 1992. [9] J. Yamauchi, M. Sekiguchi, O. Uchiyama, J. Shibayama, and H. Nakano, “Modified finite-difference formula for the analysis of semivectorial modes in step-index optical waveguides,” IEEE Photon. Technol. Lett., vol. 9, no. 7, pp. 961–963, 1997. [10] C. Vassallo, “1993–1995 optical mode solvers,” Optical Quantum Electron., vol. 29, no. 2, pp. 95–114, 1997. [11] A. S. Sudbø, “Film mode matching: A versatile numerical method for vector mode field calculations in dielectric waveguides,” J. Eur. Opt. Soc. A, Pure Appl. Opt., vol. 2, pp. 211–233, 1993. [12] M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” Proc. Inst. Elect. Eng., vol. 135, no. 1, pt. J, pp. 56–63, 1988. [13] C. Vassallo, “Interest of improved three-point formulas for finitedifference modeling of optical devices,” J. Opt. Soc. Amer. A, vol. 14, no. 12, pp. 3273–3284, 1997. [14] H. A. Van Der Vorst, “Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear system,” SIAM J. Sci. Statist. Comput., vol. 13, pp. 631–644, 1992.

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

2464

[15] D. Yevick and B. Hermansson, “New formulations of the matrix beam propagation method: Application to rib waveguides,” IEEE J. Quantum Electron., vol. 25, no. 2, pp. 221–229, 1989. [16] C. L. Xu, W. P. Huang, and S. K. Chaudhuri, “Efficient and accurate vector mode calculations by beam propagation method,” J. Lightwave Technol., vol. 11, no. 7, pp. 1209–1215, 1993. [17] J. Shibayama, M. Sekiguchi, J. Yamauchi, and H. Nakano, “Eigenmode analysis of optical waveguides by an improved finite-difference imaginary-distance beam propagation method,” IEICE Trans., vol. J81C-I, no. 1, pp. 9–16, 1998. [18] R. Pregla, “MOL-BPM method of lines based beam propagation method,” in Progress in Electromagnetics Research (PIERS 11), W. P. Huang, Ed. Cambridge, MA: EMW, 1995, ch. 2. [19] S. M. Lee, “Finite-difference vectorial-beam-propagation method using Yee’s discretization scheme for modal fields,” J. Opt. Soc. Amer. A, vol. 13, no. 7, pp. 1369–1377, 1996. [20] J. Yamauchi, N. Morohashi, and H. Nakano, “Techniques of improving the imaginary-distance full-vectorial BPM based on Yee’s mesh,” in Integrated Photonics Research. Victoria, B.C., Canada, 1998, pp. 56–58. [21] P. L¨usse, P. Stuwe, J. Sch¨ule, and H.-G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol., vol. 12, no. 3, pp. 487–494, 1994. [22] G. R. Hadley, “Low-truncation-error finite difference equations for photonics simulation—I: Beam propagation,” J. Lightwave Technol., vol. 16, no. 1, pp. 134–141, 1998.

Junji Yamauchi (M’85) was born in Nagoya, Japan, on August 23, 1953. He received the B.E., M.E., and Dr.E. degrees from Hosei University, Tokyo, Japan, in 1976, 1978, and 1982, respectively. From 1984 to 1988, he was a Lecturer in the Electrical Engineering Department, Tokyo Metropolitan Technical College, Japan. Since 1988, he has been a member of the Faculty of Hosei University, where he has been a Professor of Electronic Informatics. His research interests include optical waveguides and circularly polarized antennas. Dr. Yamauchi is a member of the Institute of Electronics, Information and Communication Engineers (IEICE) of Japan and the Optical Society of America (OSA).

JOURNAL OF LIGHTWAVE TECHNOLOGY, VOL. 16, NO. 12, DECEMBER 1998

Gen Takahashi was born in Kagoshima, Japan, on May 9, 1974. He received the B.E. degree from Hosei University, Tokyo, Japan, in 1997 and is currently pursuing the M.E. degree. Mr. Takahashi is a member of the Institute of Electronics, Information and Communication Engineers (IEICE) of Japan.

Hisamatsu Nakano (M’75–SM’87–F’92) was born in Ibaraki, Japan, on April 13, 1945. He received the B.E., M.E., and Dr.E. degrees in electrical engineering from Hosei University, Tokyo, Japan, in 1968, 1970, and 1974, respectively. Since 1973, he has been a member of the Faculty of Hosei University, where he is now a Professor of electronic informatics. His research topics include numerical methods for antennas, electromagnetic wave scattering problems, and lightwave problems. He has published more than 130 refereed journal papers and 90 international symposium papers on antenna and relevant problems. He is the author of Helical and Spiral Antennas, (New York: Research Studies Press, Wiley, 1987). He published “Antenna Analysis Using Integral Equations,” in Analysis Methods of Electromagnetic Wave Problems, vol. 2 (Norwood, MA: Artech House, 1996). He was a Visiting Associate Professor at Syracuse University, Syracuse, NY, during May–September 1981, a Visiting Professor at the University of Manitoba, Canada, during March–September 1986, and a Visiting Professor at the University of California, Los Angeles, during September 1986–March 1987. Dr. Nakano received an International Scientific Exchange Award from the Natural Sciences and Engineering Research Council of Canada. In 1987, he received the Best Paper Award from the IEEE 5th International Conference on Antennas and Propagation. In 1994, he received the IEEE Antennas and Propagation Society Best Application Paper Award (H. A. Wheeler Award). He is an Associate Editor of IEEE ANTENNAS AND PROPAGATION MAGAZINE.

Authorized licensed use limited to: HOSEI UNIVERSITY KOGANEI LIBRARY. Downloaded on August 31, 2009 at 23:53 from IEEE Xplore. Restrictions apply.

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & close