## Abstract

This work presents a novel full-vectorial imaginary-distance beam propagation method based on the multidomain pseudospectral scheme, for the first time, to study the modal characteristics of dielectric optical waveguides. The proposed method divides the transverse plane into several subdomains with uniform refractive indices, and expands the optical field in each subdomain in terms of a set of suitable basis functions. Accordingly, the complicated cross-coupling terms, which are required by the finite difference or finite element schemes, can be removed from the full-vectorial formulations. However, the coupling effect can be restored by matching the physical interface conditions. Moreover, to identify the higherorder modes, the residual traces of the preceding lower-order modes are subtracted from the calculated optical fields in each propagation step to suppress the rapid growth of the lower-order modes. Numerical examples of the two-dimensional slab waveguides demonstrate that the present approach yields a highly accurate propagation constant. When applied to three-dimensional rib waveguides the present scheme yields results that remarkably agree with the reliable values obtained from the modal transverse resonance method and finite element scheme.

©2008 Optical Society of America

## 1. Introduction

Developing an efficient and accurate computational scheme for analyzing the modal characteristics of optical waveguides is a crucial subject in the design of diverse photonic integrated circuits. Such an extensively adopted scheme is the beam propagation method (BPM), which solves the paraxial wave equation by applying the fast Fourier transformation [1], finite difference (FD) [2,3] or finite element (FE) [4,5] methods. As well as addressing the propagation of light in longitudinally variant waveguide structures, Yevick and Hermansson [6] proved that when light propagates along an imaginary axis in BPM over a sufficiently long distance, the fundamental mode can be determined from the largest amplification factor. The proposed approach, which is based on the fast Fourier transformation, is called imaginary-distance BPM (ID-BPM).

Based on the idea in another work [6], Xu *et al.* [7] generalized the ID-BPM to a full-vectorial FD-ID-BPM to compute the fundamental vector modes of optical fiber. By assessing the stability of their scheme and calculating the percentage errors in effective refractive index and mode distribution, Xu *et al.* [7] indicated that the ID-BPM yields accurate solutions for arbitrary waveguide structures. To extend the application of ID-BPM [6], Chen and Jüngling [8] proposed the subtraction of the previously obtained lower-order modes from the onset of the starting field and the end of the simulation to determine the higher-order modes. Nowadays, many versions of full-vectorial ID-BPM including ID-BPM, based on the Yee-mesh FD [9,10] and FE [11,12] approaches have been developed to analyze many complex structures, including sloped-side rib waveguides, light-guiding metal lines, electrooptic modulators, and photonic crystal fibers.

In recent years, a novel pseudospectral method (or called spectral collocation method) with high-order accuracy and fast convergence originally developed to solve computational fluid dynamics problems [13], has been introduced to the area of computational electromagnetics [14–16]. For optical waveguide problems, the authors recently developed the semivectorial and full-vectorial pseudospectral mode solvers based on a matrix eigenvalue equation [17,18]. To achieve the same order of accuracy as the FD and FE schemes, the pseudospectral method requires considerably less memory for storage than the FD and FE ones and computational time because the required matrix equations are fewer. Additionally, the pseudospectral method also yields the sparse matrix pattern as those formed by the FD and FE methods. Thereafter, Chiang *et al.* [19] further combined the pseudospectral method with a curvilinear mapping technique to solve waveguides with general curved interfaces. However, the BPM based on the pseudospectral method can be used only to analyze simple slab waveguides [20,21] and no pseudospectral-based ID-BPM has yet been developed.

This work aims to develop a new full-vectorial multidomain pseudospectral-based ID-BPM (MP-ID-BPM) for the eigenmode analysis of complex three-dimensional waveguide structures. The nonuniform mesh in various subdomains makes the present approach more efficient than the one [19] using a single set of basis functions in resolving the strongly varying field, and the better-matched basis functions in each subdomain can be selected flexibly to reduce markedly the number of grid points required. The advantage of ID-BPM over the mode solvers formulated based on a matrix eigenvalue equation [17–19] is that the modal fields calculated by the former can offer a starting field for the propagation beam analysis. Therefore, both the eigenmode and propagation beam analyses ban be performed in a single BPM simulation. The rest of this paper is organized as follows. Section II formulates the paraxial wave equations used in MP-ID-BPM. Section III presents the numerical approach and solution approach. Section IV demonstrates the high accuracy and superior convergence behavior of the proposed method. Finally, Section V draws conclusions.

## 2. Formulations

For light that propagates in the *z* direction, under the paraxial approximation, the full-vectorial wave equations based on the transverse magnetic field vector **H** are given in compact matrix form as follows;

where

*k*
_{0} is the wave number in free space; *n*
_{0} is the reference index, and *H _{x}* and

*H*are the

_{y}*x*and

*y*components of

**H**, respectively. The operators in

**P**are [22],

In the present scheme, since the computational window is partitioned into several subdomains with uniform refractive indices, *P _{xy}*=

*P*=0 and

_{yx}*P*=

_{xx}*P*=

_{yy}*P*in each subdomain, where

*P*is the scalar operator form as follows;

hence, (1) reduces to the following.

In the ID procedure [6], the propagation axis is taken as the imaginary axis (*z*=*jτ*), which yields the matrix equation

The contributions from subdomain 1 to *r*, yield the global matrix equation in diagonal block form:

where

The two components *H _{x}* and

*H*in (7) are decoupled. However, the coupling effect can be restored by considering the physical interface conditions, as shown below.

_{y}First, the normal and tangential components of magnetic fields at each intra-element boundary are continuous. Moreover, the continuities of the longitudinal components *H _{z}* and

*E*across dielectric interfaces provide the coupling effect of

_{z}*H*and

_{x}*H*through the relations ∇×

_{y}**H**=

*jωε*

_{0}

*n*

^{2}(

*x*,

*y*)

**E**and ∇·

**H**=0. For a horizontal interface between the materials, since the derivatives of

*H*and

_{x}*H*with respect to

_{y}*x*on both sides of the interface are equal, the continuity of

*E*yields

_{z}and the continuity of *H _{z}* yields

where *y*
^{+} and *y*
^{-} refer to the locations infinitesimally above and below the horizontal interface, respectively. Likewise, for a vertical interface,

and

where *x*
^{+} and *x*
^{-} refer to the locations infinitesimally to the right and to the left of the vertical interface, respectively. For further dealing with arbitrarily curved interface structure, one can introduce the curvilinear mapping technique [19]. Notably, the formulation (7) must be modified by applying the interface conditions of the equal normal and tangential components of the magnetic fields at all intra-element boundaries and the relations (9)–(12)

## 3. Numerical approaches

#### 3.1 The pseudospectral scheme

In the pseudospectral approach, the unknown fields of a differential equation are represented by a complete set of Lagrange-type interpolation functions and the grid points of fields. This differential equation must then be satisfied exactly at the collocation points [23]. For waveguides with discontinuous refractive index profiles, to preserve the exponential convergence behavior of the pseudospectral scheme, the computational window is partitioned into several subdomains with uniform refractive indices and then the subdomains are patched by matching the physical interface conditions. In each subdomain, the optical fields can be flexibly expanded not only using suitable basis functions but also using a non-uniform mesh to improve the computational efficiency. For subdomain *k*, the transverse magnetic fields (*H _{x}* and

*H*) are represented as a product of the basis function

_{y}*θ*(

*x*) in the

*x*-direction, the basis function

*ψ*(

*y*) in the

*y*-direction, and the grid point values of the

*H*and

_{x}*H*fields at (

_{y}*N*+1)×(

_{x}*N*+1) collocation points, denoted as

_{y}*H*

^{x,k}

_{p,q}and

*H*

^{y,k}

_{p,q}:

and

where *θ ^{k}_{j}*(

*x*)=

_{i}*δ*,

_{ij}*ψ*(

^{k}_{j}*y*)=

_{i}*δ*, and

_{ij}*δ*denotes the Kronecker delta. Equations (13) and (14) are substituted into (7), and (7) must be perfectly satisfied at these (

_{ij}*N*-1)×(

_{x}*N*-1) collocation points, excluding the boundary points of the subdomain

_{y}*k*. Here the operator

*P*is as given below

_{k}$$=\sum _{i=1}^{{N}_{x}-1}\sum _{j=1}^{{N}_{y}-1}[\sum _{p=0}^{{N}_{x}}\sum _{q=0}^{{N}_{y}}\{{\theta}_{p}^{k\left(2\right)}\left(x\right){\psi}_{q}^{k}\left(y\right)+{\theta}_{p}^{k}\left(x\right){\psi}_{q}^{k\left(2\right)}\left(y\right)+\dots $$

$${k}_{0}^{2}\left({n}_{k}^{2}(x,y)-{n}_{0}^{2}\right){\theta}_{p}^{k}\left(x\right){\psi}_{q}^{k}\left(y\right)\left\}\right]{\mid}_{x={x}_{i},y={y}_{j}}.$$

In (15), *θ ^{k}*

^{(h)}

*(*

_{p}*x*) and

*ψ*

^{k}^{(h)}

*(*

_{q}*y*) represent, respectively, the

*h*-th order derivatives of basis functions

*θ*(

_{p}*x*) and

*ψ*(

_{p}*x*) in subdomain

*k*. Notably, the boundary points in the

*x*-directions are replaced by (11) and (12), and those in the

*y*-directions are replaced by (9) and (10). Finally, combining the contributions from all subdomains and considering the relations (9)–(12), yield,

where

, and matrix **S** represents the final matrix to be solved herein and is no longer in diagonal block form, as in (7). For concision, the explicit form of **S** is not given here. A similar implementation to that used to obtain **S** can be found elsewhere [24].

#### 3.2 Basis functions of Lagrange-type

Various basis functions can be used to represent the optical fields in the various subdomains. For the interior subdomains with finite boundaries, Chebyshev polynomials are favored because of their robustness for non-periodic structures [23]. For the exterior subdomains including semi-infinite boundary, Laguerre-Gaussian functions (LGFs) suffice because they exponentially decay to zero at infinity. For Chebyshev polynomials, the explicit form of the Lagrange-type function is as follows [23];

where *T _{n}*(

*x*) is the general Chebyshev polynomial of order

*n*. For LGFs, the explicit form is as follows [23];

where *L _{n}*(

*αx*) is the Laguerre polynomial of order

*n*. Parameter

*α*, called the scaling factor, affects the accuracy for a given term of LGFs. The definite procedure for determining

*α*has been derived in the authors’ earlier work [18].

#### 3.3 Imaginary-distance procedure

The general solution of (16) after propagation through a distance *τ* can be expressed as,

where **Λ**(*x*,*y*,0) is the starting field. Theoretically, the starting field can be expressed as a weighted sum of eigenmodes (including the guided and radiation modes):

where *c _{i}* and

**Ψ**

*(*

_{i}*x*,

*y*) are the amplitude and field distribution of the

*i*-th order eigenmode, respectively. From the eigenvalue equation, the eigenvalue

*λ*and eigenvector

_{i}**Ψ**

*of*

_{i}**S**must satisfy the relation,

where

and *n ^{i}_{e}* is the effective refractive index of the

*i*-th order eigenmode. Substituting (21) and (22) into (20), yield

Applying the generalized Crank-Nicolson scheme to (16), yield,

where *l* is the *l*-th step along the propagation axis τ; Δ*τ* is the step size, and *ρ* is a scheme parameter. Using (24), we write

where **Λ**
* ^{i}_{l}* denotes the portion of the eigenmode

*i*involved in the field

**Λ**

*l*. For the amplification factor

*F*of the eigenmode

^{i}*i*, we have

There are three parameters *ρ*, *n*
_{0}, and Δ*τ* can be optimized to obtain the eigenmode *i* [25,26]. However, the main interesting of this work is to combine the multidomain pseudospectral method and ID-BPM. As a result, the readers can refer to the analytical relationships for determining the three optimal parameters in [25,26]. After a sufficiently long distance, the mode profile with the largest eigenvalue *λ*
_{0} dominates other higher-order modes, and therefore the starting field converges to the fundamental mode. The effective index of the fundamental mode is calculated from the growth of mode profile as follows [27].

Notably, the value *n*
_{0} at step *τ*+ Δ*τ* is updated the value *n _{e}*(

*τ*+Δ

*τ*). Because the value

*n*

_{0}is renewed in each step, only the two parameters

*ρ*and Δ

*τ*are needed to be chosen and will be discussed in the following numerical examples. To find the first-order mode, an arbitrary starting field must be subtracted from the fundamental mode in advance. However, the fundamental mode accumulates continually throughout the propagation even though its component is weak. This work proposes that the subtraction be performed not only at the onset of the input field but also in each step to eliminate perfectly the residual traces of the fundamental mode. Therefore, the new starting field

**Λ**̂(

*x*,

*y*,

*τ*) for finding the (

*N*+1)-th order mode is as follows.

The expansion coefficient *c _{i}* can be easily determined because of the orthogonality of the eigenmodes [8]:

where * means the transpose. The computational procedure of the proposed MP-ID-BPM is as follows.

Step 1: The starting field **Λ**(*x*,*y*,0), initial guess *n*
_{0}(0), step size Δ*τ*, and scheme parameter **α** must be determined in advance to yield the propagation field in the first step **Λ**(*x*, *y*, Δ*τ*) using (25).

Step 2: The effective index in the first step *n*
_{e}(Δ*τ*) is calculated using (28) and the value of *n*
_{e}(Δ*τ*) is then assigned to the reference index *n*
_{0} in the next step, *τ*=Δ*τ*.

Step 3: Substituting the updated value of *n*
_{0}(Δ*τ*) and the obtained field of **Λ**(*x*,*y*,Δ*τ*) into (25) yields the optical field in the second step **Λ**(*x*,*y*,2Δ*τ*).

Step 4: Steps 2 and 3 repeated until the difference between the effective indices in adjacent propagation steps satisfies the required tolerance.

Step 5: To find higher-order modes, Steps 1 to 4 are repeated but the starting field **Λ**(*x*,*y*,0) is replaced by the new starting field **Λ**̂(*x*, *y*,0) based on the equation (29) and (30). Importantly, the subtraction in (29) is performed in each propagation step in the simulation.

## 4. Simulation results and discussion

First, a three-layer waveguide with exact data is used to study in detail the convergence behaviors of the proposed scheme. Next, a directional coupler is analyzed to demonstrate the strong discrimination by method for almost degenerate modes. For the complicated rib waveguide, the fundamental modes calculated using the proposed approach are compared with those obtained using the modal transverse resonance method (MTRM) and FEM with high-order elements. Finally, the higher-order modes of a wider rib waveguide are solved.

#### 4.1 The two-dimensional waveguides with step-index

The refractive indices of the core and cladding are *n _{c}*=1.5 and

*n*=1.3, respectively. The width of the core is 2µm and the optical wavelength is

_{cl}*λ*=1.3µm. The waveguide supports three transverse electric (TE) and transverse magnetic (TM) modes, which were solved by Xu

*et al.*[7] using FD-ID-BPM. To find whole eigenmodes, a Gaussian beam, with beam waist

*w*

_{0}=1µm, centered on the interface between the core and cladding, was selected as the starting field. In practice, the starting field can be chosen arbitrarily. Of course, the propagation distance is strongly affected by the similarity of the starting field and the eigenmode to be found. In particular, the higher-order modes cannot be found if the starting field is centered symmetrically in the core. In the proposed approach, the transverse computational window is divided into three subdomains. The optical field in the interior subdomain is expanded by Chebyshev polynomials, and those in the exterior subdomains are expanded by LGFs with the accurately determined scaling factor, as derived elsewhere [18]. The benefit of applying LGFs to exterior subdomains is that no additional boundary conditions have to be incorporated into the approach because of the close matching of the LGFs and the guided modes. In solving the leaky modes, the absorbing boundary conditions have to be introduced to deal with the traveling waves, and LGFs must be replaced by Chebyshev polynomials to represent the optical field in exterior subdomains [26].

In this example, the reference index *n*
_{0}=1.4 and the scheme parameter *ρ*=0.65 are adopted in the generalized Crank-Nicolson scheme. Fig. 1(a) and (b) present the relative error in the effective refractive index of the fundamental TE and TM modes, respectively, as a function of the propagation distance, with step size Δ*z*=1µm, where N is the number of terms in the basis function for each subdomain. Notably, the reference index can be arbitrarily chosen because of the introduction of the renewed procedure, as in (26).

The relative errors for both TE and TM modes are smaller than 10^{-7} while only eight terms of the basis functions are used in each subdomain (giving a total of 24 unknowns). Increasing the number of terms in the basis functions to N=12 (for a total of 36 unknowns), yields errors of the order of 10^{-12}. Additionally, the effect of the scheme parameter *ρ* is also studied to determine an optimal range for it while the step size Δ*τ*=1 is fixed. Fig. 2(a) and (b) plot the relative error of the fundamental TE and TM modes, respectively, as a function of the propagation distance for various *ρ*. The range *ρ*=0.6 to 0.65 yields the most accurate and efficient results. In the numerical experiments, at *ρ*>0.65, the relative errors increase somewhat.

In the subsequent examples, *ρ*=0.65 is set. After the fundamental mode has solved, this approach is applied to deal with the higher-order modes. Subtracting the obtained fundamental mode in each propagation step to filter out perfectly the fundamental mode causes the starting field to converge to the first-order mode. The same procedure yields the second-order mode by removing the fundamental and first-order modes. Even if the propagation distance is short, the growth of the lower-order modes may dominate the convergent field. Therefore, the subtraction must be performed in each propagation step to prevent efficiently the fast growth of the lower-order modes. Fig. 3(a) and (b) plot the relative errors in the first- and second-order TE modes, respectively as functions of the propagation distance with step size Δ*z*=1 µm. Figure 1 and 3 indicate that the higher-order modes require longer propagation distances and slightly more terms in the basis functions to achieve the same relative error of the effective index.

To validate further the resolution of the proposed approach, a directional coupler with two identical 5µm-wide waveguides, separated by 11µm is considered. The refractive indices of the core and substrate at wavelength *λ*=1.32µm are *n _{c}*=1.60 and

*n*=1.59, respectively. The reference index used is the same as that of the substrate. The directional coupler supports four guided modes that result from the splitting of each individual waveguide that supports two guided modes. The difference between the effective indices of the first two modes is extremely small because of the far separation of the waveguides. The transverse domain is divided into five subdomains. The optical field in each interior subdomain is expanded by Chebyshev polynomials with 20 terms and that in each exterior subdomain is expanded by LGFs with ten terms. The two advantages merits, flexible nonuniform mesh and the use of more matched basis functions in various subdomains- make the proposed approach more efficient in terms of computational time and memory storage. In the propagation direction, to find the symmetric modes (the fundamental and second-order modes) the selected starting field is the superposition of the two identical Gaussian beams (as used in the first example) that are individually centered in each waveguide. Unlike for the symmetric modes, the starting field used for the asymmetric modes (the first- and third-order modes) is the superposition of two inverse Gaussian beams. Table 1 presents the effective indices of the four TE and TM guided modes calculated using the proposed approach. In particular, the difference of the effective refractive index between the fundamental mode TE

_{sub}_{0}(TM

_{0}) and that of the first-order one TE

_{1}(TM

_{1}) are of the order of 10

^{-10}. Fig. 4(a) to (d) display the four TE modal fields obtained herein. The above calculated results demonstrate that the proposed approach can discriminate among the almost degenerate modes, even though a coarse grid with a total of 80 unknowns was used.

#### 4.2 Semiconductor rib waveguides

For the three-dimensional waveguides, the fundamental modes of a benchmark rib waveguide whose cross section is as shown in Fig. 5(a) are analyzed to demonstrate the superior performances of the proposed approach. The geometric parameters are *w*=3µm, *h*+*t*=1µm, and the refractive indices are *n _{c}*=3.44,

*n*=3.4 and

_{s}*n*=1.0 at optical wavelength

_{a}*λ*=1.15µm. The transverse plane of the rib waveguide is divided into 12 subdomains, as displayed in Fig. 5(b).

The reference index is *n*
_{0}=3.4, and the starting field centered in subdomain 5 is a Gaussian beam with beam waist *w*
_{x}=*w*
_{y}=1µm in the *x* and *y* directions. In finding the fundamental quasi-TE mode *H*
^{y}_{11} (for which the major field is *H _{y}* and the minor field is

*H*) and quasi-TM mode

_{x}*H*

^{x}_{11}(for which the major field is

*H*and the minor field is

_{x}*H*), the starting field is assigned to the major field, and the minor field set to zero. The optical field in each subdomain is expanded using basis functions with ten terms in both

_{y}*x*and

*y*directions. Fig. 6 plots the effective indices of

*H*

^{y}_{11}and

*H*

^{x}_{11}calculated by taking

*t*=0.1µm using the proposed method against the propagation distance with step size Δ

*z*=1 µm. According to the calculations, the difference between the effective indices at adjacent steps is smaller than 10

^{-6}at distance of

*z*>150µm.

The advantage of the non-uniform mesh in the proposed approach can be exploited to adjust the distribution of grid points, according to the mode profiles, to reduce efficiently the computational time and required memory for storage. In Fig. 5(b), the optical fields expanded in the *x*-direction in subdomains 1, 3, 4, 6, 7, 9, 10, and 12 and the *y*-direction in subdomains 1, 2, 3, 10, 11, and 12 are fixed using basis functions with ten terms because the mode profile is concentrated in subdomains 5 and 8. Refining the mesh in the transverse plane yields the convergence results of the two modes that are presented in Table 2. The parameter *b* in Table 2 is the normalized propagation constant, defined as *b*=(*n*
^{2}
_{eff}-*n*
^{2}
* _{s}*)/(

*n*

^{2}

*-*

_{c}*n*

^{2}

*), and the N is the number of basis functions used in the*

_{s}*x*- direction in subdomains 2, 5, 8, and 11, and the

*y*-direction in subdomains 4 to 9.

The total number of terms in the basis functions used in the example is of (10+N+10)×(10+N+N+10) in the *x*- and *y*- directions. The achieved convergence values of effective indices are the order of 10^{-6} when the basis functions used have N>14 for both *H ^{y}*

_{11}and

*H*

^{x}_{11}modes. For various values of

*t*, the effective indices in the

*H*

^{y}_{11}mode obtained using the present approach are compared with those presented elsewhere [27] and [28], and presented in Table 3.

In another work [27], Hadley derived highly accurate solutions based on the full vector finite difference method by dealing correctly with the behavior of dielectric corners. The accuracy of the effective index claimed in that work [27] is of the order of 10^{-6} in effective index. In another work [28], the four-digit values of *b* based on the full vector MTRM are believed to be exact. Table 3 indicates that the results obtained in this work agree with those obtained using the highly accurate finite difference method [27] and MTRM [28]. Table 4 present the values for the *H ^{x}*

_{11}mode calculated in this work along with those obtained using the MTRM [28] and FEM [29]. In another work [29], the computational window of the rib waveguide consisted of 2500 triangles, and high-order T

_{15}edge elements with six tangential and two inner normal unknowns were used. Table 4 reveals good agreement between the calculated results for the

*H*

^{x}_{11}mode obtained using the three approaches.

Next, the eigenmode analysis is extended to the higher-order modes of the rib waveguide, as shown in Fig. 5(a) at *t*=0.5µm, but the width is now changed to *w*=6µm. The waveguide supports three quasi-TE and TM modes. Yamauchi *et al.* [30] considered this example using Yee-based finite difference full-vectorial BPM. The computational window adopted in their work [30] was 10µm in the *x*-direction and 6µm in the *y*-direction, and the sample widths in the transverse directions were 0.05µm, yielding their presented convergence results [30]. Therefore, the grid mesh is 200×120 in the *x*- and *y*- directions, respectively. In the scheme herein, however, the starting field and reference index used are the same as were used in the last case but the location of the starting field is moved moderately from the *y*-axis in this case by 4µm to ensure the convergence of the higher-order modes. Table 5 presents the calculated effective indices of the six modes for the mesh distribution with N=20 (N is defined in Table 2). The corresponding grid mesh is 40×60 in the *x*- and *y*- directions, respectively. In this example, the convergent effective indices, which reveal that the differences between the adjacent steps are smaller than 10^{-6}, are obtained at distances of *z*=200µm and *z*=600µm for the fundamental and higher-order modes, respectively.

Figures 7(a) and (b) present the major and minor fields of the second-order TE mode (*H ^{y}*

_{31}), respectively, for comparison with the field profiles obtained by Yamauchi

*et al.*[30]. The fields in the

*H*

^{y}_{31}mode are accurately computed and agree with those of Yamauchi

*et al.*[30]. Notably, the computational efficiency of the proposed approach (with 40×60) markedly exceeds that of Yamauchi

*et al.*[30] (with 200×120). Fig. 8(a) and (b) also present the major and minor fields of the second-order TM mode (

*H*

^{x}_{31}), respectively. Importantly, the starting field must be shifted from the

*y*-axis because the symmetric starting field has no component of a higher-order mode.

Accordingly, if only the fundamental modes are solved, then the position of the starting field can be located arbitrarily. On the contrary, the position of the starting field for finding the higher-order modes can be determined only after some trial computations. In summary, the proposed approach is highly accurate but requires many fewer grid points than the extensively used FE [29] and FD [30] approaches.

## 5. Conclusion

This work combines the full-vectorial multidomain pseudospectral discretization and the imaginary-distance beam propagation method to compute the modal fields and propagation constants of dielectric optical waveguides. The multidomain technique makes the proposed approach highly efficient and accurate because it is inherently of high-order. Moreover, the non-uniform mesh and flexible choice of the matched basis functions improve further the computational efficiency and memory storage. The proposed approach also offers the high resolution to discriminate among the almost degenerate modes with a difference in the effective indices of the order of 10^{-10}. In solving the rib waveguides, the pseudospectral discretization requires much less grid mesh but provides faster convergence than the finite difference and finite element approaches. With respect to the propagation direction, the orthogonality of the eigenmodes is exploited and the preceding calculated lower-order modes are subtracted from the starting field to yield the next higher-order mode. Importantly, the subtraction is performed in each propagation step of the simulation to ensure that the residual traces are completely eliminated. The procedure can efficiently prevent field growth by the accumulation of the lower-order modes even when the propagation distance is short. Future work will develop a leaky mode solver by combining the proposed approach with the efficient perfectly matched layer (PML) absorbing boundary condition.

## Acknowledgments

The authors would like to thank the National Science Council of the Republic of China, Taiwan, for financially supporting this research under Contract No. NSC 97-2221-E-275-001. Ted Knoy is appreciated for his editorial assistance.

## References and links

**1. **M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fiber,” Appl. Opt. **17**, 3990–3998 (1978). [CrossRef] [PubMed]

**2. **W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite difference beam propagation method: analysis and assessment,” J. Lightwave Technol. **10**, 295–304 (1992). [CrossRef]

**3. **W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. **29**, 2639–2649 (1993). [CrossRef]

**4. **Y. Tsuji, M. Koshiba, and T. Shiraishi, “A finite element beam propagation method for strongly guiding and longitudinally varying optical waveguides,” J. Lightwave Technol. **14**, 217–222 (1996). [CrossRef]

**5. **Y. Tsuji, M. Koshiba, and T. Shiraishi, “Finite element beam propagation method for three-dimensional optical waveguide structures,” J. Lightwave Technol. **15**, 1728–1734 (1997). [CrossRef]

**6. **D. Yevick and B. Hermansson, “New formulation of the matrix beam propagation method: Application to rib waveguides,” IEEE J. Quantum Electron. **25**, 221–229 (1989). [CrossRef]

**7. **C. L. Xu, W. P. Huang, and S. K. Chaudhuri, “Efficient and accurate vector mode calculations by beam propagation method,” J. Lightwave Technol. **11**, 1209–1215 (1993). [CrossRef]

**8. **J. C. Chen and S. Jüngling, “Computation of high-order waveguide modes by imaginary-distance beam propagation method,” Opt. Quantum Electron. **26**, S199–S205 (1994). [CrossRef]

**9. **T. Ando, H. Nakayama, S. Numata, J. Yamauchi, and H. Nakano, “Eigenmode analysis of optical waveguides by a Yee-mesh-based imaginary-distance propagation method for an arbitrary dielectric interface,” J. Lightwave Technol. **20**, 1627–1634 (2002). [CrossRef]

**10. **J. Shibayama, T. Yamazaki, J. Yamauchi, and H. Nakano, “Eigenmode analysis of a light-guiding metal line loaded on a dielectric substrate using the imaginary-distance beam-propagation method,” J. Lightwave Technol. **23**, 1533–1539 (2005). [CrossRef]

**11. **S. S. A. Obayya, B. M. A. Rahman, and H. A. El-Mikati, “New full- vectorial numerically efficient propagation algorithm based on the finite element method,” J. Lightwave Technol. **18**, 409–415 (2000). [CrossRef]

**12. **K. Saitoh and M. Koshiba, “Full-vectorial imaginary-distance beam propagation method based on a finite element scheme: application to photonic crystal fibers,” IEEE J. Quantum Electron. **38**, 927–933 (2002). [CrossRef]

**13. **C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, *Spectral methods in fluid dynamics, Springer Series in Computational Physics* (Springer Verlag, New York. 1988).

**14. **B. Yang and J. S. Hesthaven, “A pseudospectral method for time-domain computation of electromagnetic scattering by bodies of revolution,” IEEE Trans. Antennas Propagat. **47**, 132–141 (1999). [CrossRef]

**15. **Q. H. Liu, “A pseudospectral frequency-domain (PSFD) method for computational electromagnetics,” IEEE Antennas Wireless Propagat. Lett. **1**, 131–134 (2002). [CrossRef]

**16. **J. H. Lee and Q. H. Liu, “A 3-D Spectral-Element Time-Domain Method for Electromagnetic Simulation,” IEEE Trans. Microw. Theory Tech. **55**, 983–991 (2007). [CrossRef]

**17. **C. C. Huang and C. C. Huang, “An Efficient and Accurate Semivectorial Spectral Collocation Method for Analyzing Polarized Modes of Rib Waveguides,” J. Lightwave Technol. **23**, 2309–2317 (2005). [CrossRef]

**18. **C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Select. Topics Quantum Electron. **11**, 457–465 (2005). [CrossRef]

**19. **P. J. Chiang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. **44**, 56–66 (2008). [CrossRef]

**20. **I. Deshmukh and Q. H. Liu, “Pseudospectral beam-propagation method for optical waveguides,” IEEE Photon. Technol. Lett. **15**, 60–62 (2003). [CrossRef]

**21. **C. C. Huang and C. C. Huang, “A novel wide-angle beam propagation method based on the spectral collocation scheme for computing tilted waveguides,” IEEE Photon. Technol. Lett. **17**, 1872–1874 (2005). [CrossRef]

**22. **W. P. Huang, C. L. Xu, and S. K. Chaudhuri, “A vector beam propagation method based on H fields,” IEEE Photo. Technol. Lett. **3**, 1117–1120 (1991). [CrossRef]

**23. **J. P. Boyd, *Chebyshev and Fourier Spectral methods. Lecture Notes in Engineering* (New York, Springer, Verlag, 2nd Edition, 2000).

**24. **P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E **75**, 026703 (2007). [CrossRef]

**25. **S. Jungling and J. C. Chen, “A study and optimization of eigenmode calculations using the imaginary-distance beam propgation method,” IEEE J. Quantum Electron. **30**, 2098–2105 (1994). [CrossRef]

**26. **M. M. Spuhler, D. Wiesmann, P. Freuler, and M. Diergardt, “Direst computation of higher-order propagation modes using the imaginary-distance beam propagation method,” Opt. Quantum Electron. **31**, 751–761 (1999). [CrossRef]

**27. **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. **81**, 9–16 (1998).

**28. **C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express **14**, 11631–11652 (2006). [CrossRef] [PubMed]

**29. **G. R. Hadley, “High-accuracy finite-difference equations for dielectric waveguide analysis II: Dielectric corners,” J. Lightwave Technol. **20**, 1219–1232 (2002). [CrossRef]

**30. **C. Vassallo, “1993–1995 optical mode solvers,” Opt. Quantum Electron. **29**, 95–114 (1997). [CrossRef]

**31. **S. Selleri and J. Petracek, “Modal analysis of rib waveguide through finite element and mode matching methods,” Opt. Quantum Electron. **33**, 373–386 (2001). [CrossRef]

**32. **J. Yamauchi, T. Mugita, and H. Nakano, “Implicit Yee-based finite difference full-vectorial beam propagation method,” J. Lightwave Technol. **23**, 1947–1955 (2005). [CrossRef]