An Efficient Method for Computing Optical Waveguides With Discontinuous Refractive Index Profiles Using Spectral Collocation Method With Domain Decomposition

An Efficient Method for Computing Optical

Waveguides With Discontinuous Refractive Index

Profiles Using Spectral Collocation Method With

Domain Decomposition

Chia-Chien Huang, Chia-Chih Huang, and Jaw-Yen Yang

Abstract—An accurate and efficient solution method using

spec-tral collocation method with domain decomposition is proposed for computing optical waveguides with discontinuous refractive index profiles. The use of domain decomposition divides the usual single domain into a few subdomains at the interfaces of discontinuous refractive index profiles. Each subdomain can be expanded by a suitable set of orthogonal basis functions and patched at these in-terfaces by matching the physical boundary conditions. In addi-tion, a new technique incorporating the effective index method and the Wentzel–Kramers–Brillouin method for the a-priori determi-nation of the scaling factor in Hermite–Gauss or Laguerre–Gauss basis functions is introduced to considerably save computational time without experimenting with it. This method shares the same desirable property of the spectral collocation method of providing a fast and accurate solution but avoids the oscillatory solutions and improves the poor convergence problem of the simple spec-tral collocation method with single domain where regions of dis-continuous refractive index profiles exist. Applications to several two- and three-dimensional waveguide structures having exact or accurate approximate solutions are given to test the accuracy and efficiency; all the results are found to be in excellent agreement.

Index Terms—Discontinuous refractive index profiles, domain

decomposition, optical waveguides, orthogonal basis function, scaling factor, spectral collocation method.

I. INTRODUCTION

R

ECENTLY, there has been growing interest in developing integrated optics devices in optical communication sys-tems. The fundamental operating principle and optimal design of optical devices such as modulators, switches, filters, fibers, and semiconductor lasers are often based on the comprehension of optical waveguide theory [1]. Hence, solving the optical waveguide problems is an essential way to clearly realize these devices. However, only a limited number of waveguides, which have simple geometry and specific refractive index profiles (RIPs), have analytic solutions [1], [2]. Consequently, employing

Manuscript received June 9, 2003. This work was supported by the National Science Council, Taiwan, R.O.C., under Contracts NSC 89-2212-1-002-151 and NSC 90-2212-1-002-203.

C.-C. Huang and J.-Y. Yang are with the Computational Electromagnetics and Plasma Laboratory, Institute of Applied Mechanics, National Taiwan University, 10617 Taipei, Taiwan, R.O.C. (e-mail: [email protected]).

C.-C. Huang is with the Department of Electrical Engineering and Graduate Institute of Electro-Optical Engineering, National Taiwan University, 10617 Taipei, Taiwan, R.O.C.

Digital Object Identifier 10.1109/JLT.2003.816895

approximate or numerical methods to obtain propagation char-acteristics of general waveguides is definitely necessary and of practical importance.

There exist many effective and accurate numerical methods for solving optical waveguide problems. These include the beam propagation method (BPM), [3], [4], finite difference (FD), [5], [6], and finite element (FE) methods [7], which all produce good results. In the past decade, the series expan-sion method has became very popular and has been widely applied in studying the optical waveguide problems due to its superior accuracy and lower computational storage because of the smaller size of matrix equations required. The spectral collocation method (SCM) [8] (i.e., the collocation method in [9] proposed by Sharma and Banerjee) and the Galerkin method [10] are the two most prevalent methods in series expansion. More recently, many researchers [11]–[15] have successfully implemented the Galerkin method with different orthogonal basis functions to investigate optical waveguides and optical fibers but few pursued SCM. Sharma and Banerjee [9], [16] first applied SCM with Hermite–Gauss (HG) basis functions in propagation characteristics of optical waveguiding structures for two kinds of problems as follows. First, for the waveguide modal problems, SCM attained the same degree of convergence compared with the Galerkin method for given terms of orthogonal basis functions but needed much less computational time, resulting from the avoidance of evaluating laborious integral elements as undertaken by the Galerkin method [16]. Second, the propagation problems using SCM also provided better convergence than BPM [16].

These advantages of SCM produced good results only ap-plicable under the condition of the continuous RIP in interest domain. As for discontinuous RIP, the SCM solutions show the oscillatory phenomenon (the so-called Gibbs phenomenon) and induce difficulty to converge, as illustrated in [16]. Causa

et al. [17], [18] also used SCM with HG in some cases, such as

longitudinally nonuniform, Kerr-nonlinear media and nonlinear carrier diffusion in semiconductor optical devices. Because of the existence of discontinuous RIP, even though the difference of refractive index between different mediums is small, too many terms of HG still were needed (the number of terms of HG are ) to achieve desirable results. As mentioned above, although SCM is powerful and simple, there still exists a major drawback when dealing with the discontinuous RIP 0733-8724/03$17.00 © 2003 IEEE

in interest domain. For the solution of a differential equation solved through SCM in terms of a complete set of orthog-onal basis functions, the expansion coefficients are determined by requiring the expansion to satisfy the differential equation exactly at certain discrete points, which are the so-called collo-cation points. Meanwhile, RIP is indicated as an interpolation of the refractive index at these collocation points; therefore, the discontinuous RIP cannot always be described accurately at the discontinuous interfaces. With this notable flaw, the merit of the high accuracy of SCM has not been completely exploited, and it can result in poor convergence due to the existence of discontinuous functions in computational domain [8]. The larger the difference of refractive index at discontin-uous interfaces is, the larger the error in describing the actual RIP in SCM may become.

In this paper, our purpose is first to propose SCM with domain decomposition (DD) to deal with the discontinuous RIP in the interest domain, while SCM with only single domain usually produces oscillatory results and poor convergence [16]; second to decide the suitable set of orthogonal basis functions that naturally describes the electromagnetic characteristics at each subdomain; and finally to derive a technique for choosing the scaling factor when using HG or Laguerre–Gauss (LG) basis functions (LGs are usually used in optical fiber of circular structure [19]) avoiding iterative experimenting with it and enhancing the accuracy. By combining these techniques, the accuracy of solution can be significantly improved and the computational time of finding a scaling factor is substantially reduced.

This paper is organized as follows. In Section II, we describe the mathematical formulations of the wave equations. The de-scription of the SCM with DD is stated in Section III. Section IV analyzes in detail the choice of the scaling factor. In Section V, we first compare two cases calculated by Banerjee and Sharma [16] with the results obtained by the present method; then a few examples that include planar and diffused channel waveguides solved by analytical or highly accurate approximate methods are also tested. After that, we discuss the effect of scaling factor. Section VI concludes this paper.

II. MATHEMATICALFORMULATIONS OFWAVEEQUATIONS Assuming monochromatic electromagnetic fields with a time dependence of , propagating along the -direction in an inhomogeneous medium with refractive index , the vector wave equation can be derived from the Maxwell’s equa-tions for a source free region in frequency domain as follows:

(1) Similarly, the vector wave equation based on the magnetic field vector is given by

(2) where and is the wavelength in free space. Con-sidering the longitudinally (i.e., -direction) invariant structure, namely, , the refractive index depends only on the

two transverse directions, that is, . Here, we are only concerned with the guided mode solutions, with a periodic dependence of the form , where is the propaga-tion constant. When the two polarizapropaga-tions of and (or and ) are weak for three-dimensional waveguides, the cou-pling terms may be neglected [4]. Hence we can individually deal with the uncoupled equations termed as the semivector [7], [8] or quasi-vector [15] analysis, as follows:

(3) for the quasi-TE modes and

(4) for the quasi-TM modes.

Here is the effective refractive index. The second terms of (3) and (4) are polarization dependence. Furthermore, if the weakly guiding approximation is assumed in the whole interest domain, then we can neglect polarization dependence to obtain the scalar wave equation

(5) that considerably reduces the computational time. For three-di-mensional (3-D) structures, we only compute a diffused channel waveguide in this paper, which has been solved by Sharma and Bindal [20] assuming scalar approximation. In two-dimensional (2-D) structures such as planar waveguides, the derivatives with respect to in (3) and (4) are zero if the infinite extension is along the -direction. Hence, we obtain the formulas

(6) for TE modes and

(7) for TM modes.

III. SPECTRALCOLLOCATIONMETHODWITHDOMAIN DECOMPOSITION

In this section, we describe how to solve the optical wave-guide problems mentioned previously (5)–(7) by SCM with DD. Apparently, to deal with the inaccurate description of discontin-uous RIP at interfaces encountered by SCM with single domain [16], we need to divide the whole domain into a few subdomains at these interfaces. According to the feature of each subregion, we may employ different orthogonal basis functions to expand the optical fields in individual subregion. For the internal re-gions, whose boundaries are finite, Chebyshev polynomials are the favored choice because of their robustness for nonperiodic structure [8]. For the outer regions, which include semi-infinite boundaries, LGs are adequate due to the mathematical charac-teristic of the exponentially decay to zero at infinity.

For illustration purposes, we confine our description to planar waveguide structures with the transverse spatial direction . For

any dependent variable such as electric field in (6) or magnetic field in (7) in each internal region, we expand

it by Chebyshev polynomials with

(8)

where are the expansion coefficients and in each outer sub-domain. We expand it by

(9)

where denotes the LGs, i.e.,

, , and is

the Laguerre polynomial of degree . In (9), is a scaling factor to be determined and can influence the accuracy for a given number of terms of basis functions (discussed in the next section). The independent variable is discretized into

1 collocation points ( ). An

alterna-tive expression is that the dependent field variable is directly expanded by the cardinal basis functions [8], [21] and the grid point values of fields

at collocation points as follows:

(10)

where

(11)

Here, is a basis function of order 1 and the

prime denotes the first derivative with respect to . are the collocation points with the conditions , where is the Kronecker delta. Thus the computations are operated in physical space instead of in spectral space . The field profiles then can be directly obtained from the unknowns in (10). In this paper, we use the latter form (10) because of the direct obtainment of the field profiles . The explicit forms of for distinct basis functions and collocation points are represented below. For instance, in (11), if we take Chebyshev polynomials as basis functions and Chebyshev Gauss–Lobatto points as collocation points (i.e., the extrema of Chebyshev polynomials together with the endpoints at ) for internal finite subdomains, then we have

(12)

where denotes Chebyshev polynomial of order ,

, and ( ). Correspondingly,

for outer infinite (or semi-infinite) subdomains, we take LG,

namely, as the basis

functions, and the roots of Laguerre polynomial of

order plus as collocation points to facilitate the incor-poration of boundary conditions. We have

(13)

Finally, using HG (i.e.,

) as basis functions and

the roots of Hermite polynomials of order 1 as collocation points, we have

(14)

The basis functions are adequate and efficient only for the con-tinuous RIP along transverse spatial directions using SCM with usual single domain as illustrated in [16]. Using SCM, we substitute (10), based on either (12) for internal regions or (13) for outer regions, into (6) and (7), respectively. Requiring that (6) and (7) be satisfied at these 1 collocation points, these lead to

(15)

and

(16) Subsequently, (6) and (7) are transformed into 1 algebraic eigenvalues (15) and (16) for each subdomain. We obtain in ma-trix form

(17)

where is an 1 1 square matrix.

are grid point vectors (eigenvectors), where denotes the transpose. The elements of matrix are

(18) for TE modes and

(19) for TM modes. In the above formulas, denotes the th derivative of cardinal basis function of degree with respect to . The numerical examples of planar waveguides in this paper are mainly based on (18) and (19).

After expanding different basis functions in internal and outer subdomains, we need to employ the boundary conditions of electromagnetic field at these interfaces , ({ }, where is the number of discontinuity), between different di-electric constants to patch these subdomains as follows:

for TE modes and

(21) for TM modes. Here, and are referred to the locations at the infinitesimally right and left of the interface , respectively. As a result of special care taken to treat the discontinuity of RIP using DD, the discontinuous RIP in interest domain can be correctly represented and the oscillatory solution in usual SCM with single domain can be efficiently eliminated.

Considering the three-dimensional waveguide structures, the dependent variable in any subdomain is expanded by a product of two separable basis sets and as follows:

(22)

where and denote two cardinal basis functions of order and along - and -directions, respectively, and now are grid point values in two transverse directions. Using the similar procedures as mentioned above,

after substituting (22) into (5), algebraic

eigenvalue equations are obtained.

IV. CHOICE OFSCALINGFACTOR

An important point to note is that any orthogonal basis func-tion for an infinite interval implicitly contains a scaling factor to set the extension of the mapping, such as HG or LG [8]. Thus reasonably choosing can substantially reduce the requiring number of basis functions and evade many additional trial com-putations with it. In SCM, the scaling factor is usually picked up through trial and error with additional expense of a consider-able amount of computational time for a given number of terms of the basis functions [16]. Hence a definite way to determine the optimum value of the scaling factor is thus very desirable and necessary. For given terms of basis functions , decides the spread and the extent of collocation points in the transverse direction. In general, there exists an adequate for a given in optical waveguide problems.

In this paper, we adopt techniques from Tang [22] and Ramanujam et al. [23] and combine them to offer a new novel approach for determining for HG and LG. In [22], Tang proposed an interesting and definite procedure to choose for Gaussian-type functions that is expanded by the first 1 HG of order using the expression of series coefficients as given by (9) in spectral space. Here, another approach, which is based on the expression given by (10) in physical space, is used. A more convenient and direct expression than those given in [22] may be derived to obtain an appropriate . For any dependent variable in outer regions, we use as basis functions. For Chebyshev polynomials used for internal regions, because of the finite interval of , no scaling factor is provided to adjust. In contrast to HG or LG, can be independently accommodated in each outer region. We write

(23)

where , , are the grid point values as in

(10). At each collocation point , we obtain

(24)

which is equivalent to

(25)

For the HG basis functions, are the zeros of

with the order , while for LG, are

the zeros of plus . When solving the differential equations, the th derivatives of dependent variable with respect to are also needed and are expanded by

(26)

where ( ) is the so-called th

differentiation matrix (DM) [24] of the given basis functions. Suppose that the function has a finite support ,

i.e., for . Thus has sufficient

contribution to in (25) or in (26) under the

condition

for all

Accordingly, the following relation determines the scaling factor associated with given terms

(27) The above choice of the scaling factor is due to Tang [22], which he originally applied to known Gaussian-type functions. However, the value for the finite support remains to be determined. The finite support in Tang’s work was assumed to be given by a simple truncation of infinite domain into a finite support. In [23], a novel technique has been proposed in Galerkin method with trigonometric basis functions (the Fourier method) for the a-priori determination of the finite support (or numerical window in [23]) of the Fourier method by using the Wentzel–Kramers–Brillouin (WKB) method and effective index method (EIM). Here we adopt this technique to determine for in (27) in our SCM with DD.

According to the character of guided modes in optical wave-guide theory [1], the mode shapes of all wave-guided modes oscil-late within the core and exponentially decay beyond the turning point. For step-index waveguides, the turning point of all guided modes locates at the same discontinuous interface that is be-tween the outermost and its adjacent layers. As for graded index waveguides, the turning point of each guided mode locates on the different position, where the refractive index is equal to the effective index. If we can obtain the effective index of the highest guided mode, the turning points and the decay rates then can be estimated. Indeed, the characteristic of the optical field is to extend to infinity; thus, we have to set a criterion to decide

the finite support beyond which the field strength is suffi-ciently small and can be neglected. For diffused waveguides, Hocker and Burns [25]–[27] have used WKB method to treat the diffused planar waveguides and combined EIM and WKB to extend to the 3-D diffused channel waveguides.

Suppose we have the following 2-D RIP:

(28) where is the substrate index, is the maximum index at ( ), and is the cover index. The functions

and 2 are the normalized profiles of the graded index. Here and are the parameters of diffusion depth in - and -direction, respectively. First, we use the EIM for the profile in the -direction at each position and define the normalized guide depth and effective index parameters [27] as follows:

(29) and

(30) where is the effective index, which depends on . Then for a specified value of employing the WKB theory, we can find by the dispersion relation [26, (10)] as follows:

(31)

where , , and is the turning point, which

satisfies at a given ; and ,

denotes the order of mode along the -direction. Once we

have , we can obtain at each through (30).

Here, the phase shift at the turning point ( ) is 2, and that at the interface of air and material ( ) is

2 . As shown in [26],

denotes the asymmetry index for a given , and its value is

variant for TE ( ) and TM

( ) modes,

. A set of data ( ) then can be obtained point by point through the repeated computations from (29)–(31), and the equivalent planar wave-guide to the original channel wavewave-guide can be represented by the symmetric index profile along the -direction. The normalized guide depth and effective index of the equivalent planar waveguide can be defined by

(32) and

(33)

where is the effective index of the original channel wave-guide and denotes the value of the symmetric effective index profile at . Apply the WKB method again to the -dependent RIP . The dispersion relation equations [27, (18)] determining for the equivalent symmetric planner waveguides are as follows:

(34)

where , , is the turning point in

-di-rection, which satisfies , and ,

rep-resents the mode number along the -direction. Once we have obtained , can be determined by (33) and is calculated through . Using the same definition as [23], we denote the decay rates for the mode as

(35a) (35b) (35c) Accordingly, we can estimate the decay rates of optical fields , along the -direction into cover and substrate sides, respectively, and is along the -direction into substrate. For determining the finite support or numerical window , we need to specify a criterion to assume that the field strength at is much less than the field strength at turning points. Moreover, the ratio of the field strength at relative to the field strength at turning points has to be adequately and reasonably decided. If the ratio is too small, then the spreading we calculate is too wide to reduce computational time and the number of the basis function. Contrarily, the large ratio represents the field strength at may not be negligible. In this paper, following Ramanujam

et al. [23], we take the ratio as 0.01. Hence, to find the finite

support in each direction, we first let

(36a)

(36b)

(36c)

as indicated in Fig. 1, and determine , , and from (36a)–(36c). Then the finite support along -direction is

. For the -direction, and

are the finite supports along cover and substrate sides, respectively. We want to emphasize here that the a-priori choice of the numerical window only needs the roughly approximate values, because the main function of scaling factor in HG or LG is to redistribute the collocation points corresponding to (no scaling) to more efficiently describe the field profiles. Thus the tolerance of may not be so strict.

After that, we can adjust the density of collocation points by choosing suitable and concerning more collocation points near the regions where the field values change explosively. It

(a)

(b)

Fig. 1. Illustration of modal profiles for the cross-sectional view of a diffused channel waveguide having symmetric and asymmetric graded index profiles alongx- and y-directions in (a) and (b), respectively. In (a), the dotted lines represent the location of the turning points (at6x ), is decay rate, and x + 1x is position of spreading for a given mode. For (b), the turning points are at zero andy , where , denote decay rates and y 0 1y , 1y are spreadings for substrate and air sides, respectively.

is noted that in [23], the numerical window is used in Galerkin method with trigonometric functions. Due to the characteristics of trigonometric basis set, the boundary conditions at infinity of optical field cannot be satisfied; hence the requirement to decide numerical window is severe and the accuracy strongly depends on the parameter . In the present method, the HG or LG can automatically satisfy the boundary conditions at infinity and no enclosed window is needed. Indeed, the a-priori choice of scaling factor can substantially reduce the computational time and efficiently improve the convergence.

V. NUMERICALRESULTS ANDDISCUSSION

We first discuss two numerical examples in detail that allow the existence of discontinuous RIP in interest domain, as studied in [16], to demonstrate the significant improvement of conver-gent behavior of effective refractive index using DD and the a-priori choice of scaling factor in SCM. Second, to examine

Fig. 2. Refractive index profile for a three-layer step-index waveguide divided into three subdomains by the vertical dotted lines, where the refractive indexes of core and cladding aren = 1:45 and n = 1:447 636, respectively, and W = 3 m at  = 1:3 m.

(a)

(b)

Fig. 3. Schematic diagrams of refractive index profiles for asymmetric graded index waveguides divided into two regions aty = 0. (a) Gaussian profile in region II,n = 1:512, n = 1, and the diffused depth d = 2:66 m at wavelength = 0:6328 m. (b) Exponential profile in region II, n = 2:47, n = 1, d = 2:5 m, and the wavelength  = 0:633 m.

the accuracy of this method, several waveguides are investi-gated, including asymmetric graded index waveguide, which possesses high-order modes, metal-clad waveguide, planar di-rectional coupler, and diffused channel waveguide. Finally, we discuss the effect of scaling factor.

TABLE I

CONVERGENCE OF THEEFFECTIVEINDEXn OF ATHREE-LAYERSTEP-INDEXSTRUCTURE AT THEWAVELENGTH = 1:3 m

TABLE II

CONVERGENCE OF THEEFFECTIVEINDEXn OF AASYMMETRICGRADEDINDEXSTRUCTURE AT THEWAVELENGTH = 0:6328 m

A. The Fundamental Mode of Slab Waveguide for Symmetric Step Index and Asymmetric Graded Index Profiles in [16]

We study two examples that have been solved by Banerjee and Sharma in [16]. The first example is a three-layer sym-metric step-index slab waveguide (see Fig. 2), which supports single mode. The refractive indexes of core and claddings are

and , respectively. The half-width

of the core is m, and the optical wavelength is m. We divide the whole domain into three subdomains at the two interfaces of discontinuous RIP as indicated in Fig. 2. Region I is expanded by Chebyshev polynomials, because of the finite boundary; however, region II and region III are both ex-panded by LG due to the semi-infinite extension in each region. Using the transcendental equation of TE mode for three-layer step-index waveguide [2] to find the effective index, and ac-cording to (35) and (36), the decay rates and spreadings are then

obtained. Here, we have the spreadings m and

m for the two outer regions II and III, respec-tively. In Table I, , , and represent the number of terms of orthogonal basis functions in each region. We first fix and simultaneously adjust and because of the sym-metry. When , the convergent mode value only achieves 1.448 978 no matter how many terms of and are used. After increasing to , the result nearly converges to the exact solution 1.448 90 associated with . After

that, (namely, scaling factors for region

II and region III, respectively) can be evaluated by (27). Com-paring to the results of [16, Table 4], which were obtained using

terms of HG in single domain, our solutions obvi-ously eliminate the oscillatory phenomenon and can achieve the exact solution using merely 15 terms of orthogonal basis func-tions altogether.

The second example is an asymmetric graded index profile slab waveguide as depicted in Fig. 3(a). Here, the RIP is a Gaussian profile described as follows:

(37)

where , , , and m

at m. The interest domain is divided into two subdomains, and we obtain the spreadings m for region I and m for region II, which are calcu-lated by the WKB method (through (29)–(31) at ). In [16], only the fundamental mode of the example is solved. Even though the number of the terms of HG is up to in usual single domain, the effective index seems to have difficulty in convergence. Similar to the first example, we set and increase the number of to 8. The exact result is attained with a total of 12 terms of LG, as given in Table II. At the

same time, and are determined. From

these two examples, we can see the superior convergence of the present method and its ability to completely exclude oscillatory solutions. Meanwhile, the a-priori choice of scaling factor has

TABLE III

COMPARSION OF THEEFFECTIVEINDEXESn FORTEANDTM MODES AT = 0:633 m

been implemented. Otherwise, more computational effort must be spent, if determined through trial and error.

B. The High-Order Modes of Asymmetric Planar Waveguide With Exponential RIP

To investigate optical waveguides, including high-order modes is also important in precisely designing many devices such as multimode fiber and multimode interference used in dense wavelength-division multiplexing (DWDM). For examining the accuracy for high-order mode problems, we solved an asymmetric slab waveguide with exponential RIP, as shown in Fig. 3(b), which supports eight TE and TM modes. The exact solutions have been obtained in [28]. The RIP is described as follows:

(38)

where is the substrate index, is the cover

index, is the maximum index change, m is

the depth of diffusion, and m is the wavelength, re-spectively. The interest domain is divided into two subdomains. In our calculation for TE modes, the spreadings are

m and m. The number of LG in region

I is and in region II is , for which and are obtained. Likewise, for TM modes,

m and m are calculated, ,

are used, and we have and . The exact

[28] and approximate [29] solutions of effective indexes of TE and TM modes together with ours are listed in Table III. In ad-dition, the results of TM modes solved by expanding both the index profile and modal field into a Fourier series [30] are also added in. There were a few modifications for the exact effec-tive indexes [28] calculated by Noro and Nakayama displayed in [29] using the characteristic equations [28, (8) and (13)] for TE and TM modes, respectively. Most of the data are correct, but the third order (m ) and the sixth order (m ) of TM modes show slight deviation. We recalculated them using secant method, and found that the effective index of the third order is 2.486 660 2 and sixth order is 2.473 250 7. After rounding off all the data of effective indexes, the values are shown in Table III.

Fig. 4. TE mode profiles obtained by the present method (filled circles) versus exact solutions for modes 0 (solid line), 4 (dashed line), and 7 (dotted line) at the wavelength = 0:633 m.

Fig. 5. Metal-clad guide with exponential index profile in region II, wherea is the depth of diffusion and the index of cover layern = 010:3 0 i1:0 in region I is complex.

TABLE IV

COMPARSION OF THECOMPLEXPROPAGATIONCONSTANTS OFTEANDTM MODES FORASYMMETRICMETAL-CLADWAVEGUIDEWITHEXPONENTIAL

INDEXPROFILE FORa = 5 m

The present method achieves better accuracy than [29] and [30], which used terms of Fourier series, and shows the in-distinguishable result with exact solution within 30 terms of LG. Fig. 4 shows the mode profiles of modes 0 (solid line), 4 (dashed line), 7 (dotted line) of TE by exact method [28] and our solu-tions (filled circles); they are all in excellent agreement. From the results of this case, our method can accurately provide the ef-fective indexes and mode profiles for optical waveguide modes simultaneously, even when they support higher order modes.

C. The Metal-Clad Optical Waveguides

We have illustrated the accuracy obtained by our method for lossless optical waveguides in above cases. Here, the SCM with DD is extended to analyze the lossy optical or the so-called metal-clad optical waveguides. The metal-clad waveguides provide many applications in integrated optics devices, for example, polarizers and mode filters. Many researchers have investigated the metal-clad waveguides theoretically and ex-perimentally for step index [31], [32] and graded index profiles [33]–[35]. In this case, we solve the exponential RIP like case B but the refractive index is replaced by the complex refractive index , as illustrated in Fig. 5. This problem has been studied by an analytical, highly accurate numerical method in [33] and by a Galerkin method with trigonometric basis functions in [35]. The profile is given as follows:

(39)

where the waveguide index parameters are ,

, the gold cladding, and is the depth of dif-fusion at the wavelength m. Regions I and II in Fig. 5 are both expanded by LG. Because of the difficulty of finding the spreading of the imaginary part of the optical field in advance, we only consider the real part of to approxi-mately calculate the spreading of the real field of the highest mode by WKB [(35) and (36)]. Although this idea is not rig-orous, the small error is acceptable for the approximate predic-tion of spreading due to the fact that the spreading of real fields is close to the imaginary fields, as illustrated by [32, Fig. 7]. In addition, because the difference of refractive indexes at the in-terface ( ) of metal and surface of material is large (i.e.,

Re ), the phase shift at may be

assumed as [25], [26]; thus we obtain the same values of

Fig. 6. Schematic diagram of refractive index profile for the two parallel slab waveguide structures (located in regions II and III) divided into five subdomains by the four vertical dotted lines, whereW = 2 m, S = 1:95 m, n = 2:19, andn is variable.

spreadings of TE and TM modes. Accordingly, the spreadings

of the highest mode (m ) are m in region I

and m in region II. We used , ,

, and for m.

In Table IV ( m), our solutions still show highly accu-rate results as compared with the accuaccu-rate numerical solutions in [33] no matter whether the real or image part of propagation constants is encountered. Here, the total number of terms of LG in the present method is within 20. In [35] (only TE modes were solved), even using Fourier series, the solutions were still inferior. This case exhibits that the optimal scaling factors

( and ) are extremely hard to achieve

through trial and error. Furthermore, through the DD, the scaling factor can be independently determined in each subdomain; un-like SCM applied in single domain, is the same in the whole computational domain and is chosen inflexibly. Hence, the con-vergence for our method is achieved more quickly than SCM using only single domain.

D. The Planar Directional Coupler

A directional coupler is an interesting optical device that often is used as a beam splitter, filter, optical switcher, etc. Due to the wide applications in optical communications, it is necessary to understand how a directional coupler works. In general, the

simplest directional coupler consists of two parallel slab waveg-uides, as shown in Fig. 6. When a directional coupler is used as an optical switcher, the modulation fields in the two waveguides (located in regions II and III in Fig. 6) are usually applied with an equal amount of electric field but with opposite sign to change refractive indexes due to the electrooptic effect of crystals, resembling lithium niobate (LiNbO ). The refractive index of the right waveguide (region III) is then increased by 2 and decreased by 2 in the left waveguide (region II). Generally, is often a tensor quantity for anisotropic medium, but here, only an isotropic change of refractive index is considered. The widths of the two waveguides are m. For the zero modulation fields at wavelength m, the refractive indexes of the two waveguides are . The separation between the guides is m. Except the guides, the re-fractive indexes are everywhere. Here, we divide the whole domain into five subdomains at these discontinuous in-terfaces. The spreadings calculated in regions IV and V are ob-tained by regarding the guides as isolated (i.e., weakly coupled approximation) in this structure. However, the spreadings we

obtained are m for all the different due

to the weakly coupled approximation, which are calculated by the same dispersion relation for a three-layer step index wave-guide as adopted in the first example in case A. If the coupling is strong, we need additional analysis to evaluate the spreading. The optical fields under three various refractive indexes induced

by different modulation fields , ,

are illustrated in Fig. 7(a)–(c), respectively. In Fig. 7(a)–(c), the real lines and dashed lines represent the even and odd modes, respectively. Under variant applied fields, the induced change of refractive indexes is indicated in Table V, and the propagation constants of the even ( ) and odd ( ) modes calculated by the present method using 31 basis func-tions altogether are exact compared with analytical solufunc-tions and the solutions obtained by finite element method [36] using 4000 mesh divisions. Here, we used seven Chebyshev

polyno-mials in regions I, II, and III ( ) and five

LG in regions IV and V ( ). The present method can be straightforwardly extended to some devices containing anisotropic and layered waveguide structures.

E. The Diffused Channel Waveguide

In the last case, this method is extended to investigate dif-fused channel waveguide, which is formed by diffusing titanium ( ) into LiNbO in - and -directions and hence produces graded RIP. The RIP is given by (28), where ,

, and

(40)

Here, erf denotes error function and the parameters are

m, m, , , wavelength

m, and is the normalized

frequency. Sharma and Bindal [20] listed a few results solved by scalar variational analysis using different trial fields, and the solutions obtained by [37] were acclaimed as the most accurate variational analysis by them. We tested our method using scalar

(a)

(b)

(c)

Fig. 7. Optical field profiles of the even (solid line) and odd (dotted line) modes of TE supermodes for different change of refractive index: (a)1n = 0:0, (b)1n = 0:0008, (c) 1n = 0:0020.

TABLE V

THEEVEN ANDODDTE SUPERMODES OF THECOUPLEDWAVEGUIDEARECALCULATED BYDIFFERENTMETHODS. INTHISPAPER,THESPREADING OFOUTER

REGIONSAREABOUTM = M = 5 mFORDIFFERENT1n. N = N = N = 7, N = N = 5 AREUSED TOACHIEVE THECONVERGENCEVALUES

Fig. 8. The cross-sectional view of general diffused channel waveguide divided into two regions at the interface ofy = 0, where n(x; y) is diffused refractive index profile alongx- and y-directions.

wave (5) to this case and divided the computational domain into two subdomains, region I ( )and region II ( ), as illustrated in Fig. 8.

In (22), the optical field along the -direction is expanded by HG due to the continuous RIP, and the -directions in region I and region II are individually expanded by LG due to the dis-continuous RIP the same as previous cases. The spreadings ob-tained are m for different values of , which are

nearly the same m and m

for , respectively. For the case ,

the spreading m, and m, the

con-vergence is obtained by using ( ) terms

of LG for region I along -direction into cover,

( ) terms of LG for region II along -direction into

substrate, and ( ) HG along

-direc-tion. The cross-section of modal profile at for the case calculated by this method is plotted in Fig. 9, and the characteristics of rapid decay in the interface (as illustrated by the dashed line in Fig. 9) between air and the surface of material can be seen. Our solutions of normalized propagation constants versus compared with results in [20] are shown in Fig. 10. We can see that our results are closer to the results in [37] than those obtained by evanescent secant hyperbolic (ESH) or secant hyperbolic (SH) trial fields. The propagation constants can be efficiently and accurately ob-tained because the spreadings solved by WKB and EIM yield good approximation with actual optical fields, as depicted in Fig. 9.

Fig. 9. Modal profiles alongy-direction (at x = 0) for the case V = 3:179.

Fig. 10. Normalized propagation constantb versus the normalized frequency V of diffused channel waveguide obtained by this method and a few variational analysis using different trial fields SH, ESH, and the more accurate variational analysis in [37].

Finally, we note that for all the cases studied in this section, we have used the a-priori choice of scaling factor in HG or LG

for optical waveguide problems and obtained efficient conver-gence to effective indexes. In principle, decides the extension of collocation points. Hence, for a given number of terms of basis functions, if is chosen too small, the too-sparse density of collocation points is used to accurately describe the modal profiles, which change explosively. In contrast, if is too large, the narrower distribution of collocation points cannot properly take into account the contribution of long extension of higher order modes. In some previous works, using Galerkin method with HG [13], [15] or with LG [19] for step index or trun-cated parabolic optical fibers, a fixed value of the scaling factor was chosen, where is the refrac-tive index of the core, is the index of cladding, and rep-resents the radius of the core for step index circular fiber or is half the width of channel waveguides. However, for asymmetric diffused waveguides of general profiles, this scaling factor cannot be certainly assigned for the region I in Fig. 3(a) and (b) or Fig. 5 due to the indefinite parameter . Accordingly, the a-priori determination of by the present method not only can be used in step index or symmetric graded index waveguides but also can be used in asymmetric graded index waveguides.

VI. CONCLUSION

We have proposed an efficient solution method that com-bines spectral collocation method and domain decomposition for computing optical waveguides with discontinuous refrac-tive index profile. Through domain decomposition, one can flexibly expand different basis functions in distinct subdomains depending on the extension of each region. In addition, we have presented a new approach to choose the scaling factor to optimally determine the spreading of the highest guided mode of optical waveguides by combining Tang’s work [22] and the technique proposed by Ramanujam et al. [23]. The major factor to improve the convergence is to divide the whole do-main into a few subdodo-mains at the interfaces of discontinuous refractive index profile to enable that interpolation of the refrac-tive index at these collocation points can faithfully represent the discontinuous refractive index profile. Consequently, the oscillatory solutions obtained by spectral collocation method with usual single domain [16] can be completely eliminated and the degree of convergence can be improved. Moreover, the a-priori choice of scaling factor in Hermite–Gauss or La-guerre–Gauss functions materially reduces the computational time without trial and error. Merging these two techniques, namely, the domain decomposition and the optimal choice of scaling factor, into conventional spectral collocation method renders the present solution method both accurate and effi-cient. Our results have demonstrated significant improvement of convergence over the solutions obtained using usual single domain for spectral collocation method. A wide range of wave-guide problems has been tested, such as planar, metal-clad, directional coupler, and diffused channel waveguides, and all the results obtained are found in excellent agreement with exact or other best available accurate approximate solutions. Also, significant saving of computational time has been il-lustrated. For the three-dimensional channel waveguide, only

scalar wave equation is solved here. However, this method can be further extended to semivector or full-vector formulas [4] for investigating more complicated devices, which require the polarization dependence and coupling terms to be considered. Extension of this method of full-vectorial modal analysis of dielectric waveguides will be reported elsewhere.

REFERENCES

[1] A. W. Snyder and J. D. Love, Optical Waveguides Theory: Chapman and Hall, 1983.

[2] T. Tamir, Guides-Wave Optoelectronics. Berlin, Germany: Springer-Verlag, 1988.

[3] W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite difference vector beam propagation method-Analysis and assessment,”

J. Lightwave Technol., vol. 10, pp. 295–305, Mar. 1992.

[4] 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, Oct. 1993.

[5] M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” Proc. Inst. Elect. Eng. J., vol. 135, pp. 56–63, Feb. 1988.

[6] , “Semivectorial polarized H field solutions for dielectric waveg-uides with arbitrary index profiles,” Proc. Inst. Elect. Eng. J., vol. 135, pp. 333–338, Oct. 1988.

[7] B. M. A. Rahman and J. B. Davies, “Finite-element solution of integrated optical waveguides,” J. Lightwave Technol., vol. LT-2, pp. 682–688, Oct. 1984.

[8] J. P. Boyd, Chebyshev and Fourier Spectral Methods. New York: Springer-Verlag, 1989, Lecture Notes in Engineering.

[9] A. Sharma and S. Banerjee, “Method for propagation of total fields or beam through optical waveguides,” Opt. Lett., vol. 14, pp. 94–96, Jan. 1989.

[10] S. G. Mikhlin and K. L. Smolitskiy, Approximate Methods for Solution

of Differential and Integral Equations. New York: Elsevier, 1967, pp. 250–252.

[11] J. P. Meunier, T. Pigeon, and J. N. Massot, “A numerical technique for the determination of propagation characteristics of inhomogeneous planar optical waveguides,” Opt. Quantum Electron., vol. 15, pp. 77–85, 1983.

[12] C. H. Henry and B. H. Verbeek, “Solution of the scalar wave equation for arbitrarily shaped dielectric waveguides by two-dimensional Fourier analysis,” J. Lightwave Technol., vol. 7, pp. 308–313, Feb. 1989. [13] R. L. Gallawa, I. C. Goyal, Y. Tu, and A. K. Ghatak, “Optical

wave-guide modes: An approximate solution using Galerkin’s method with Hermite-Gauss basis functions,” IEEE J. Quantum Electron., vol. 27, pp. 518–522, Mar. 1991.

[14] D. Marcuse, “Solution of vector wave equation for general dielectric waveguides by the Galerkin method,” IEEE J. Quantum Electron., vol. 28, pp. 459–465, Feb. 1992.

[15] A. Weisshaar, J. Li, R. L. Gallawa, and I. C. Goyal, “Vector and quasivector solutions for optical waveguide modes using efficient Galerkin’s method with Hermite-Gauss basis functions,” J. Lightwave

Technol., vol. 13, pp. 1795–1800, Aug. 1995.

[16] S. Banerjee and A. Sharma, “Propagation characteristics of optical waveguiding structure by direct solution of the Helmholtz equation for total fields,” J .Opt. Soc. Amer. A, vol. 6, pp. 1884–1894, Dec. 1989. [17] F. Causa, J. Sarma, and R. Balasubramanyam, “A new method for

com-puting nonlinear carrier diffusion in semiconductor optical devices,”

IEEE Trans. Electron Devices, vol. 46, pp. 1135–1139, June 1999.

[18] F. Causa and J. Sarma, “A versatile method for analyzing paraxial optical propagation in dielectric structures,” J. Lightwave Technol., vol. 18, pp. 1445–1452, Oct. 2000.

[19] J. P. Meunier, J. pigeon, and J. N. Massot, “A general approach to the numerical determination of modal propagation constants and field dis-tribution of optical fibers,” Opt. Quantum Electron., vol. 13, pp. 71–83, 1981.

[20] A. Sharma and P. Bindal, “Variational analysis of diffused planar and channel waveguides and directional couplers,” J. Opt. Soc. Amer. A, vol. 11, pp. 2244–2248, Aug. 1994.

[21] D. Funaro, Polynomials Approximation of Differential Equa-tions. Berlin, Germany: Springer-Verlag, 1992.

[22] T. Tang, “The Hermite spectral method for Gauss-type functions,” SIAM

[23] N. Ramanujam, L. Li, J. J. Burke, and M. A. Gribbons, “Determination of the truncation order and numerical window for modeling general di-electric waveguides by the Fourier method,” J. Lightwave Technol., vol. 14, pp. 500–508, Mar. 1996.

[24] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral

Methods in Fluid Dynamics. New York: Springer Verlag, 1988, Springer Series in Computational Physics.

[25] G. B. Hocker and W. K. Burns, “Modes in diffused optical waveguides of arbitrary index profile,” IEEE J. Quantum Electron., vol. QE-11, pp. 270–276, June 1975.

[26] G. B. Hocker, “Strip-loaded diffused optical waveguides,” IEEE J.

Quantum Electron., vol. QE-12, pp. 232–236, Apr. 1976.

[27] G. B. Hocker and W. K. Burns, “Mode dispersion in diffused channel waveguides by the effective index method,” Appl. Opt., vol. 16, pp. 113–118, Jan. 1977.

[28] E. M. Conwell, “Modes in optical waveguides formed by diffusion,”

Appl. Phys. Lett., vol. 23, pp. 328–329, Sept. 1973.

[29] H. Noro and T. Nakayama, “A new approach to scalar and semivector mode analysis of optical waveguides,” J. Lightwave Technol., vol. 14, pp. 1546–1556, June 1996.

[30] L. Wang and C. S. Hsiao, “A matrix method for studying TM modes of optical planar waveguides with arbitrary index profiles,” IEEE J.

Quantum Electron., vol. 37, pp. 1654–1660, Dec. 2001.

[31] A. Reisinger, “Characteristics of optical guided modes in lossy waveg-uides,” Appl. Opt., vol. 12, pp. 1015–1025, May 1973.

[32] I. P. Kaminow, W. L. Mammel, and H. P. Weber, “Metal-clad waveg-uides: Analytical and experimental study,” Appl. Opt., vol. 13, pp. 396–405, Feb. 1974.

[33] S. J. Al-Bader and H. A. Jamid, “Guided wave characteristics of metal-clad graded-index planar optical waveguides: Analytical ap-proach,” IEEE J. Quantum Electron., vol. QE-23, pp. 539–544, May 1987.

[34] M. Masuda, A. Tanji, Y. Ando, and J. Koyama, “Propagation losses of guided modes in anoptical graded-index slab waveguides with metal-cladding,” IEEE Trans. Microwave Theory Tech., vol. MTT-25, pp. .773–776, Sept. 1977.

[35] W. Y. Lee and S. Y. Wang, “Guided-wave characteristics of optical graded-index planar waveguides with metal-cladding: A simple analysis method,” J. Lightwave Technol., vol. 13, pp. 416–421, Mar. 1995. [36] T. Wongcharoen, B. M. A. Rahman, and K. T. V. Grattan, “Electro-optic

directional coupler switch characterization,” J. Lightwave Technol., vol. 15, pp. 377–382, Feb. 1997.

[37] A. Sharma and P. Bindal, “An accurate variational analysis of single-mode diffused channel waveguides,” Opt. Quantum Electron., vol. 24, pp. 1359–1371, 1992.

Chia-Chien Huang was born in I-Lan, Taiwan, R.O.C., on August 21, 1972. He received the B.S. degree from the Department of Chemical Engi-neering, National Tsing-Hua University, Hsin-chu, Taiwan, in 1994, and the M.S. and Ph.D. degrees from the Institute of Applied Mechanics, National Taiwan University, Taipei, Taiwan, in 1998 and 2003, respectively.

In August 2003, he joined the faculty of the Department of Information Management, Ling Tung College, Taichung, Taiwan, where he is currently an Assistant Professor. His current research interests are numerical methods for optical waveguide and optical fiber devices.

Chia-Chih Huang was born in I-Lan, Taiwan, R.O.C., on July 2, 1971. He received the B.S. degree in electrical engineering from National Taiwan University, Taipei, Taiwan, in 1993, the M.S. degree in electrical engineering from National Tsing-Hua University, Hsin-chu, Taiwan, in 1995, and the Ph.D. degree in the graduate institute of electrooptical engineering from National Taiwan University.

Since August 2003, he has been with the Depart-ment of Electronic Engineering, Tung Nan Institute of Technology, Taipei, Taiwan, where he is an Assis-tant Professor. His current research interest includes modeling and design of passive integrated optical devices.

Jaw-Yen Yang received the Ph.D. degree from Stan-ford University, StanStan-ford, CA.

He was a National Research Council Associate at NASA Ames Research Center. He is currently a Pro-fessor in the Institute of Applied Mechanics at the National Taiwan University. His research focuses on computational modeling of electromagnetic wave in-teractions with complex structures, plasma stability, quantum hydrodynamics and aerodynamics.