Improved Three-Point Formulas Considering the
Interface Conditions in the Finite-Difference Analysis
of Step-Index Optical Devices
Yih-Peng Chiou, Yen-Chung Chiang, and Hung-Chun Chang, Member, IEEE, Member, OSA
Abstract—A general relation, considering the interface condi-tions, between a sampled point and its nearby points is derived. Making use of the derived relation and the generalized Douglas scheme, the three-point formulas in the finite-difference modeling of step-index optical devices are extended to fourth order accuracy irrespective of the existence of the step-index interfaces. With nu-merical analysis and nunu-merical assessment, several frequently used formulas are investigated.
Index Terms—Finite-difference method, generalized Douglas scheme, step-index optical waveguides.
I. INTRODUCTION
W
ITH the rapid progress of computers, in both software and hardware, simulation programs or computer-aided design (CAD) tools for the design of optoelectronic devices have become more and more convenient and important. Among these CAD tools, the finite-difference method (FDM) is one of the most well-known numerical methods, which is widely used in the mode solvers and beam propagation methods (BPM’s). Comprehensive reviews can be found in [1] and [2].Since the differential equations are directly approximated with their corresponding difference equations in the FDM, the efficiency and accuracy of the FDM are greatly affected by its finite-difference (FD) formulas. Various formulations have been proposed to elevate the accuracy and efficiency in the modeling. The simplest FD formula is based on the scalar approximation. Stern [3] derived vectorial formulas based on graded index approximation, in which the dielectric interface conditions between different refractive indexes were matched by means of averaging the permittivity over meshes. As will be shown later, these formulas have truncation errors where is the grid spacing and denotes that the order is th power of It should be noted that the accuracy is not elevated
Manuscript received August 9, 1999. This work was supported by the National Science Council of the Republic of China under Grant NSC88-2215-E-002-015.
Y.-P. Chiou was with the Department of Electrical Engineering, National Taiwan University. He is now with Taiwan Semiconductor Manufacturing Co., Ltd., Hsin-Chu 300, Taiwan, R.O.C.
Y.-C. Chiang is with the Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan 106-17, R.O.C.
H.-C. Chang is with the Department of Electrical Engineering, the Graduate Institute of Electro-Optical Engineering, and the Graduate Institute of Com-munication Engineering, National Taiwan University, Taipei, Taiwan 106-17, R.O.C.
Publisher Item Identifier S 0733-8724(00)01321-9.
when finer grid spacings are used. Besides, the interface is required to be in the middle between the sampled points, i.e., in Fig. 1. If the interface is not in the middle between the sampled points, the truncation errors would be large than Vassallo [4] provided improvement of the FDM for step-index optical waveguides without averaging the permittivity over meshes. The resulting formulas were derived from the Taylor series expansion and from matching of interface conditions, and gave more accurate results. The truncation errors are usually irrespective of the location of the interface with respect to the sampled points. If the interface is in the middle between the sampled points, then the truncation errors in Vassallo’s formulas are Usually, the truncation error of the commonly used three-point FD schemes can only be at best. Formulations with higher order truncation errors can be obtained when higher order terms are retained in the derivation. The generalized Douglas (GD) scheme [5]–[7] was used in the BPM’s to increase the accuracy of the FD formulas. The accuracy was elevated to when the medium is homogeneous, i.e., the refractive
indexes in Fig. 1. However, interface
conditions were not treated. The accuracy was still at best and was reduced to when the interface is not in the middle between the sampled points or the discretization is
nonuniform [7], i.e., or in Fig. 1. Yamauchi
et al. [8] derived formulas under nonuniform discretization for the GD scheme and improved Vassallo’s formulas to accuracy irrespective of the location of the interface by means of evaluating higher order terms through the BPM. The above derivations all placed the interface between the sampled points. Lüsse et al. [9] derived another formulation by placing the interfaces exactly at the sampled points, i.e., in Fig. 1. Hadley [10] derived quasifourth order equations by taking the interface conditions and nonuniformity of grid spacings into consideration. The higher order terms were evaluated through the BPM and a complicated averaging operator.
In this paper, formulas similar to Vassallo’s [4] are derived simply by the Taylor series expansion and matching the inter-face conditions. The interinter-faces need not be located exactly at the sampled points, that is, they can be placed at the grids or elsewhere. Higher order terms are retained and the GD scheme is adopted. It is also found that higher order terms are not nec-essarily evaluated through the BPM and the complicated aver-aging operator as in [8] and [10], but they can be included in a more general and simple algebraic form. The resulting eigen-value problems can be solved directly and efficiently, instead
Fig. 1. Sketch of interfaces between sampled points.
of indirectly and inefficiently using the BPM. The derived for-mulas can also be applied to the BPM with little extra compu-tation efforts, compared with the compucompu-tation in the eigenvalue problem, in finding the coefficients.
Several formulations are derived in Section II, followed by a description of numerical implementation in Section III. In Sec-tion IV the formulaSec-tion is assessed and applicaSec-tions to modal so-lutions for slab waveguides and multiple-quantum-well waveg-uides are numerically demonstrated. Section V gives the con-clusion.
II. FORMULATION
A general relation between the sampled field and the fields nearby as shown in Fig. 1 is derived in this section. It is found that Stern’s [3] and Vassallo’s [4] formulations are lower order cases of our derived formulation. Retaining the higher order terms and making use of the GD scheme, the truncation error of our formulation can be extended to irrespective of the existence of the interfaces.
A. General Formulation Considering the Interface Conditions Consider the magnetic field at a sampled point and the
fields nearby, and as shown in Fig. 1, where and
represent fields at just to the left and and just to the right sides of the interface, respectively. Using the Taylor series expansion,
is expressed as
(1)
When (1) is differentiated times successively and multiplied with we can express the derivatives of in terms of the derivatives of as
(2)
When the first terms are retained in (1) and higher order terms (H.O.T.) are ignored, we have
(3)
Equation (3) can be rewritten in a matrix form as
.. . .. . ... ... ... . .. ... .. . (4) or denoted as (5) where and denote the first and second derivatives, respec-tively, and denotes the th derivative. Similarly, and its derivatives can be expressed in terms of and its deriva-tives as
(6)
or denoted as
(7) The interface conditions require that
(8)
(9) where
for transverse electric (TE) case for transverse magnetic (TM) case.
(10) Making use of the Helmholtz equation, the relation between the higher order derivatives of and can be obtained. The Helmholtz equation for is
where is the propagation constant, is the wave number in free space, and is the refractive index. From (8) and (11), we have
(12) or
(13)
where and and are the refractive
indexes of two adjacent regions, as shown in Fig. 1. Repeating the similar process successively, higher order derivatives of can be expressed in terms of and its derivatives as
(14) (15) (16) .. . or denoted as (17) From (5), (7), and (17), and its derivatives can be expressed in terms of and its derivatives as
(18) where
(19) Note that (17) is exact without any truncation error while (5) and (7) have truncation errors, which reveals that the trunca-tion errors will not increase with the interface conditrunca-tions being included and that the truncation errors are from the neglect of higher order terms of the Taylor series expansions. Also, when (17) is multiplied with 1/θ, we have the -field formulation, which reveals that the -field and the -field formulations are equivalent, and the accuracy and efficiency are the same for these two formulations.
The relations between the above-mentioned fields can be il-lustrated as
(20)
where TSE denotes Taylor Series Expansion and MBC denotes Matching the Boundary Condition. Similarly, and its derivatives can be expressed in terms of and its derivatives as (21)
Considering the first rows of we have
(22)
where and are the first rows of and ,
respec-tively.
B. Scalar Approximation and Graded-Index Approximation The simplest three point formula is based on the scalar ap-proximation, and the vectorial nature is ignored. The second derivative is the same as that in a homogeneous medium and is expressed as
(23) and reduced to
(24)
when the grid spacing is uniform Equations
(23) and (24) are good approximations when the index contrast is low but fail when the index contrast is high. Besides, the vec-torial nature is not included.
Stern [3] derived a vectorial formulation and similar formula-tion can be derived from the graded-index approximaformula-tion, which was frequently used later in the optical waveguide simulation [11]–[13]. In their derivation,
and the second derivatives are assumed to be continuous. Thus, the formulation is the same as (24) for TE cases and
(25)
for TM cases. Since the second derivatives are assumed to be continuous, the truncation error in (22) is for
and thus the truncation errors of (24) and (25) are
which means that the truncation errors would not de-crease with the grid size with the existence of the interfaces. Equations (24) and (25) may be a good approximation in the modeling of structures with graded index or low index contrast but they are not suitable in the modeling of step-index structures. C. Improved Formulation with Truncation Errors
Vassallo [4] derived an improved formulation by making use of the Helmholtz equation. Equivalently, in the derivation. From (22), we have
(26)
where the coefficients ’s and ’s are and are given in the Appendix. is obtained by ignoring and eliminating
(28)
Since the coefficients of and are and
respectively, the approximation in (28) is Interestingly, when
(29) and
or (30)
the coefficients would have such relation
and as shown in the Appendix.
Eliminating would result in eliminating simultaneously and the approximation in (28) is thus
D. Improved Formulation of Higher Order Accuracy
In this subsection the improved formulation is derived based on the GD scheme. Again, can be expressed in terms of and its derivatives as
(31)
(32) where the coefficients ’s and ’s are and given in the Appendix. If the higher order terms containing and in (31) and (32) are ignored, then and can be solved as
(33)
(34)
Following the procedure in Section II-C, it can be found that,
usually, and Also,
when (29) and (30) are satisfied. Elimi-nating in (31) and (32) by the linear combination
we have (35) shown at the bottom of the page or denoted as
(36) which is approximated with
(37) or
(38)
Since is and
the approximation in (37) and (38) are, at worst, Similarly, when (29) and (30) are satisfied, the coefficients would have such relation
and as shown in the Appendix.
Eliminating would result in eliminating simultaneously.
Besides, and
and are Therefore, formulas with accuracy is
obtained when (29) and (30) are satisfied. III. IMPLEMENTATION
Substituting (38) into the Helmholtz equation
(39) leads to
(40) or
(41) Combining all sampled points together from (41), we have an algebraic equation
(42) or denoted as
(43)
where diag
results from the operator results from the
operator and is an identity matrix. The
matrices , and are all tridiagonal. Note that when the GD scheme is not adopted, and the formulation reduces to Vassallo’s [4]. The eigenvalue problem in (43) can be solved efficiently with the shifted-inverse power method [14].
Substituting with in (41) and assuming the
enve-lope approximation where is the
refer-ence refractive index, we have
(44) This formulation can also be applied to the beam propagation method directly as described in [7] and [8].
IV. NUMERICALRESULTS
A. Assessment of the Formulation
The investigated structure is a symmetric slab waveguide.
The core and cladding indexes are and
for GaAs and Ga Al As, respectively.
The waveguide width is m and the wavelength is
m. The effective index for the TE mode of this
wave-guide is and the transverse wave
numbers in the core and cladding are and
respectively. In the assessment of the
formulation the magnitudes of the fields, and are
assigned with exact values and then or are evaluated using (24), (25), (28), and (38). For example
(45) for our formulation, where in the core region and
in the cladding region. The relative error in the transverse wave number is defined as
(46) in the core region, or
(47) in the cladding region.
Fig. 2 shows for the TE case with respect to the grid
size for uniform grids, i.e., and
in Fig. 1, where 3pt. denotes that three-point formulation is adopted and 5pt. formulation is adopted. All the sampled points are in the core region. It can be seen that the
slopes are two and four for and formulations,
re-spectively. Note that calculated using the five-point scheme is slightly larger than that calculated using the three-point scheme. in the cladding region is of the same order as (not shown). The truncation error of our formulation is exactly
Fig. 2. Relative error in the transverse wave number with respect to the grid size: uniform grid size.
Fig. 3. Relative error in the transverse wave number with respect to the grid size: nonuniform grid size.
The distorted behavior of for is due
to finite digits in the computer.
Fig. 3 shows the relative error with respect to the average grid
size for nonuniform grids, i.e., and
in Fig. 1. The three sampled points again are all in the core region. The average grid size is defined as
(48) where the ratio in this case. is larger than that in the case of uniform grids, but still behaves as It can be seen that calculated using the five-point scheme is slightly smaller than that calculated using the three-point scheme.
in the cladding region is again of the same order as Fig. 4 shows the relative error with respect to the grid size for an interface lying between sampled points, that is,
and in
Fig. 1. It can be seen that the truncation error of our formulation is also even if the interfaces exist. The accuracy of our formulation is close to that of the uniform grid size as shown in Fig. 2. The transverse wave vector cannot be well approximated using the graded-index approximation whose truncation error is
Fig. 4. Relative error in the transverse wave number with respect to the grid size: interfaces lying between sampled points.
only and cannot be reduced with smaller grid size due to the existence of the interface.
As shown in Figs. 2–4, the accuracy is greatly enhanced when higher-order formulation is adopted. Or on the other side, the efficiency is elevated because coarser grids can be used. For example, the grid size of the formulation can be about ten times larger that of the formulation at
All the above results are for TE cases. The TM cases have also been assessed, the results being close to those for the TE cases for the same parameters.
From the above assessment of the numerical error, the deriva-tion in Secderiva-tion II is confirmed and the truncaderiva-tion error of our derived formulation is indeed We have applied the de-rived formulation to various parameters, which also validates our derivation.
B. Application to Mode Solvers: Slab Waveguides
First, a weakly guiding waveguide is investigated. The struc-ture is the same as that considered in [10]. The wavelength is 1.55µm and the waveguide width is m. The core and
cladding relative permittivities are and
respectively. The exact effective index is
The transparent boundary condition [15] is adopted.
Fig. 5(a) and (b) shows the field profiles calculated using the formulation for large grid sizes
and 1 m, respectively. The solid curves are the analyt-ical field profiles and the circles are the results calculated using our formulation. The numerical effective indexes are
and
for and 1 m, respectively. In this case the propagation constant and the modal profile can be calculated accurate to four digits even though only one sampled point is placed in the core. It seems that our formulation may be more “economic” than that of Hadley’s [10] which requires a sampled point at the each interface. Our results are accurate to the same order as those of [10] with slightly smaller error when the grid sizes are the same.
Fig. 5. Field profiles for grid size1x = (a) 2 m and (b) 1 m.
Fig. 6 shows the relative error in the reduced propagation stant with respect to the grid size. The reduced propagation con-stant is defined as [10]
(49) The dashed curve is calculated using an scheme across the interface and an scheme elsewhere. It can be seen that such rough approximation yields worse results for large grid sizes, while the effect of inaccurate sampled points near the in-terface is less weighted due to much more sampled points for small grid sizes. If the interface conditions are considered, the results become better for large grid sizes. It can be seen that the truncation error of our formulation is indeed The grid size with truncation errors can be ten times larger or more than that with truncation errors, resulting in great reduction of the computation effort.
C. Application to Mode Solvers: Multiple-Quantum-Well Waveguides
The modal characteristics of multiple-quantum-well (MQW) optical waveguides have attracted considerable attention be-cause of their distinctive features and potential applications in
Fig. 6. Relative error in the reduced propagation constant with respect to the grid size for different three-point formulations.
photonic devices. Unfortunately, those complicated structures are difficult to analyze using conventional FD techniques since the refractive index varies significantly and the widths of wells and barriers are very small and usually not uniform. The equivalent index methods [16] or the finite element methods [17] are often employed. The situation changes when the formulation of high accuracy is applied. The grid size is much more flexible and the accuracy can be much higher than those of conventional FD schemes.
The modal characteristics of a complicated MQW waveguide is investigated to show the excellent performance of our formu-lation. Fig. 7 shows the index distribution. The wavelength is
m, and the widths of wells and barriers are
m and m, respectively. The refractive
indexes of the wells, barriers, and claddings are
and for GaAs, Ga Al As,
and Ga Al As, respectively. There are 55 wells and 56 barriers in the MQW structure. The exact effective propagation
constants are and
for TE and TM modes, respectively. Only one sampled point is placed at the center of each well or barrier. The grid size is 0.0095 m in the MQW region and 0.01 m in the cladding region. Fig. 8 shows the calculated field distribution which is indistinguishable from the analytical one. The calculated effec-tive propagation constants are
and for TE and TM modes, respectively, which
are both correct to ten digits even nonuniform grid sizes are used under such strongly guiding structure.
V. CONCLUSION
We have derived improved three-point formula for the finite-difference analysis of step-index optical waveguides. A gen-eral relation is derived from the Taylor series expansion and matching the interface conditions. From the relation, the trun-cation errors are not from matching the boundary conditions but from the truncation of the higher order terms in the Taylor se-ries expansion. Also, the -field and the -field formulations are actually equivalent. Making use of the the general relation and the Generalized Douglas scheme the improved formulation
Fig. 7. Refractive index profile of the investigated MQW optical waveguide.
Fig. 8. Field profile in the MQW optical waveguide.
is derived. It is found that the frequently adopted formulations of Stern’s [3] and Vassallo’s [4] are lower-order cases of our de-rived formulation.
We have assessed the frequently used formulations and our derived formulation for various parameters. Our formulation is indeed in the uniform discretization cases irrespective of the existence of the interfaces. The graded-index approximation
is not suitable for cases with high refractive index contrast ratio since the accuracy is only When the discretization is nonuniform, the accuracy is reduced, especially in cases with high refractive contrast ratio.
The improved formulation has been applied to the mode solver of slab waveguides and multiquantum waveguides, which shows its preference in accuracy and efficiency. The derived formulation can be applied not only to the analysis of optical waveguides, but also to other electromagnetic simula-tions for structures with step index.
APPENDIX
The parameters ’s in (27) are as follows:
(50)
Replacing and with and respectively,
’s can be obtained. Specifically, when and ’s
become
(51)
Or when and ’s become
(52) The parameters ’s in (27) are as follows:
(53)
Replacing and with and respectively,
’s can be obtained. Specifically, when and ’s
become
(54)
Or when and ’s become
REFERENCES
[1] C. Vassallo, “1993–1995 optical mode solvers,” Opt. Quantum
Elec-tron., vol. 29, pp. 95–114, 1997.
[2] D. Yevick, “Some recent advances in field propagation techniques,”
Proc. SPIE, pp. 502–511, 1996.
[3] M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” Inst. Elect. Eng. Proc.-J., vol. 135, pp. 56–63, 1988.
[4] C. Vassallo, “Improvement of finite difference methods for step-index optical waveguides,” Inst. Elect. Eng. Proc.-J., vol. 139, pp. 137–142, 1992.
[5] L. Sun and G. L. Yip, “Modified finite-difference beam-propaga-tion method based on the Douglas scheme,” Opt. Lett., vol. 18, pp. 1229–1231, 1993.
[6] 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, pp. 961–963, 1997.
[7] C. Vassallo, “Interest of improved three-point formulas for finite-dif-ference modeling of optical devices,” J. Opt. Soc. Amer., vol. 14, pp. 3273–3284, 1997.
[8] J. Yamauchi, J. Shibayama, O. Saiti, O. Uchiyama, and H. Nakano, “Im-proved finite-difference beam propagation method based on the gener-alized Douglas scheme and its application to semivectorial analysis,” J.
Lightwave Technol., vol. 14, pp. 2401–2406, 1996.
[9] P. Lüsse, K. Ramm, and H.-G. Unger, “Comparison of a vectorial and new semiconductor finite-difference approach for optical waveguides,”
Opt. Quantum Electron., vol. 29, pp. 115–120, 1997.
[10] G. R. Hadley, “Low-truncation-error finite difference equations for pho-tonic simulation I: Beam propagation,” J. Lightwave Technol., vol. 16, pp. 134–141, 1998.
[11] Y. Chung and N. Dagli, “Analysis ofz-invariant and z-variant semi-conductor rib waveguides by explicit finite difference beam propagation method with nonuniform mesh configuration,” IEEE J. Quantum
Elec-tron., vol. 27, pp. 2296–2305, 1991.
[12] P.-L. Liu, S. L. Yang, and D. M. Yuan, “The semivectorial beam prop-agation method,” IEEE J. Quantum Electron., vol. 29, pp. 1205–1211, 1993.
[13] 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, pp. 2639–2649, 1993.
[14] A. Jennings, Matrix Computation for Engineers and Scientists. New York: Wiley, 1977.
[15] C. Vassallo and J. M. van der Keur, “Comparison of a few transparent boundary conditions for finite-difference optical mode-solvers,” J.
Lightwave Technol., vol. 15, pp. 397–402, 1997.
[16] M. Saini and E. K. Sharma, “Equivalent refractive index of MQW waveguides,” IEEE J. Quantum Electron., vol. 32, pp. 1383–1390, 1996.
[17] B. M. A. Rahman, Y. Liu, and K. T. V. Grattan, “Finite-element modeling of one- and two-dimensional MQW semiconductor optical waveguides,”
IEEE Photon. Technol. Lett., vol. 8, pp. 928–931, 1996.
Yih-Peng Chiou was born in Taoyuan, Taiwan,
R.O.C., on October 10, 1969. He received the B.S. degree in electrical engineering from the National Taiwan University, Taipei, in 1992 and the Ph.D. degree from the College of Electrical Engineering, National Taiwan University in 1998.
He is currently with the Taiwan Semiconductor Manufacturing Co., Ltd., Hsin-Chu, Taiwan. His research interests include new finite difference methods for solving the propagation characteristics of optical wave in dielectric waveguides.
Yen-Chung Chiang was born in Hualien, Taiwan,
R.O.C., on March 10, 1970. He received the B.S. and M.S. degrees from the Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan, in 1992 and 1994, respectively. He is currently working towards the Ph.D. degree at the same university.
Hung-Chun Chang (S’78–M’83) was born in
Taipei, Taiwan, R.O.C., on February 8, 1954. He received the B.S. degree from National Taiwan University, Taipei, Taiwan, in 1976 and the M.S. and Ph.D. degrees from Stanford University, Stanford, CA, in 1980 and 1983, respectively, all in electrical engineering.
From 1978 to 1984, he was with the Space, Telecommunications, and Radioscience Laboratory of Stanford University. In August 1984, he joined the Faculty of the Electrical Engineering Department of National Taiwan University, where he is currently a Professor. He served as Vice-Chairman of the Electrical Engineering Department from 1989 to 1991, and Chairman of the newly established Graduate Institute of Electro-Optical Engineering at the same university from 1992 to 1998. His current research interests include the theory, design, and application of guided-wave structures and devices for fiber optics, integrated optics, optoelectronics, and microwave and millimeter-wave circuits.
Dr. Chang is a member of Sigma Xi, the Phi Tau Phi Scholastic Honor Society, the Chinese Institute of Engineers, the Photonics Society of Chinese-Ameri-cans, the Optical Society of America (OSA), and China/SRS (Taipei) National Committee (a Standing Committee member during 1988–1993) and Commis-sion H of U.S. National Committee of the International Union of Radio Science (URSI). In 1987, he was among the recipients of the Young Scientists Award at the URSI XXIInd General Assembly. In 1993, he was one of the recipients of the Distinguished Teaching Award sponsored by the Ministry of Eduction of the Republic of China.
