1. Introduction

Mode calculation is a key task for simulation and analysis of optical waveguides and photonic integrated circuits. Several powerful numerical methods such as the finite difference (FD) [1–5] and the finite element (FE) [6,7] methods have been developed and widely adopted for mode calculations of practical optical waveguides. In the conventional mode calculation method, the open optical waveguides are discretized over the transverse section and the boundary-value problem is reduced to a matrix eigen-value problem defined on finite computation domain facilitated by artificial numerical boundary conditions. As such, the accuracy, computational efficiency as well as memory requirements for the conventional mode solvers depend critically on the size of the computation widow and the number of meshes. Thanks to the powerful matrix algorithm, the solutions of the eigen-value problem associated with the mode computation are normally highly efficient with an approximately linear scaling factor with the total number of meshes. One of the shortcomings of the boundary-value eigen solution is that the memory requirement scales drastically with the size of the matrix and hence become a bottleneck especially for dealing with optical waveguide problem that need large number of mesh point to achieve sufficient accuracy.

An alterative approach to mode calculation is to formulate the problem as an initial-value problem. By launching an initial field at the input of the waveguide, one may track its propagation along the waveguide axis by a field propagation algorithm such as the beam propagation method (BPM) [8–10] in frequency domain or the finite-difference time-domain (FDTD) [11,12] method in time-domain. Once the propagation reaches its steady state, a signal processing method such as the Fourier transform (FT) method [8,13] may be used to extract the modal parameters.

In this paper, the optical mode solver is based on wave equation formulation in time-domain. The compact 2D-FDTD method [12] combined with matrix pencil method (MPM) [14] is adopted for calculation of the scalar, semi-vector and full-vector modes. Uniaxial perfectly matched layer (UPML) [15–17] absorbing boundary conditions (ABCs) are employed for the FDTD scheme in order to reduce the computation region. By launching an initial optical field into the waveguide, the light is propagating along the direction (z) through the FDTD method for a certain time (t). When the steady-state is reached, a series of instantaneous response field values at the end of the propagation are adopted by MPM algorithm for mode parameter extractions. Since the MPM algorithm is much more efficient and accurate for the mode parameter extraction compared with the traditional FT method, only limited propagation distance is required in FDTD simulation. Since the FDTD method is an inherently data-parallel algorithm, we implement it in Graphics Processing Units (GPUs) by using Nvidia’s Compute Unified Device Architecture (CUDA) [18]. With the help of the high computing power brought by GPUs [19,20], the simulation time of the FDTD method is further reduced to less than 3% of that used for the implementation on a standard 3.0GHz CPU.

This paper is organized as follows. The wave Eq. (-)based compact 2D-FDTD method and the matrix pencil method for optical waveguide modal calculation are described in section 2 and 3 followed by its validation in section 4.1. The convergence test is presented in section 4.2. The performance of the enhanced method implemented on GPUs using CUDA is discussed in section 4.3. Comparisons are made among the standard FDTD, the accelerated FDTD and the conventional eigen mode solvers.

2. Wave Eq. (-)based compact 2D-FDTD method with UPML ABCs

For waveguide structures with uniform refractive index, derivatives of the mode field with respect to longitudinal direction can be replaced with $\mp j\beta$and consequently reduce one dimension complexity in space which can save significant computing time and memory consumption. This method is called compact 2D-FDTD method which was introduced by Xiao et al. [12] and has been widely used in microwaves. To avoid the complex-variables in Maxwell’s Eq. (-)based compact 2D-FDTD method, a wave Eq. (-)based compact 2D-FDTD method is introduced in [21].

In this paper, the wave Eq. (-)based compact 2D-FDTD method with UPML absorbing boundary conditions is employed. In Appendix I, the detailed formulations of the full-vector, semi-vector, and scalar schemes are summarized for the sake of completeness. Compared with Benenger’s PML [22], UPML is more preferable for parallel computing due to the reason that it avoids the nonphysical field splitting by composing a uniaxial anisotropic medium.

3. Extraction of mode parameters by matrix pencil method

When the wave equations are solved in frequency domain, the operation frequency is treated as an input parameter whereas the propagation constant β is an output in the mode calculation. In our algorithm, the wave equations are solved in time domain and hence β is selected as an input parameter so as to the mode frequencies ${\omega }_m$, the mode losses ${\xi }_m$, and the optical fields (${\Phi }_m$) for different modes are extracted in the post processing after the FDTD simulation. Since the optical waveguide has a uniform index along the propagation direction, an arbitrary transverse field (${\Psi }_t\left(x,y,t\right)$) in the waveguide can be represented by the summation of the eigenmodes ${\Phi }_m\left(x,y\right)$

In order to extract the mode parameter${\omega }_m$ for different modes, the Fourier transform can be employed to obtain the frequency spectrum from the steady-state time domain optical field captured in a window with proper width. The mode parameters ${\omega }_m$ are the positions of resonant peaks in the frequency spectrum. Since the accuracy of the mode frequency ${\omega }_m$depends on the resolution of the frequency spectrum, a large time window in FDTD simulation is required to achieve a high accuracy, which is extremely computation intensive. It is also very difficult to obtain the mode loss ${\xi }_m$ accurately by fitting Lorentz shapes centered at mode frequencies in the frequency spectrum. The matrix pencil method (MPM) is proposed to overcome the drawback of the Fourier transform method so that the mode frequency ${\omega }_m$ and mode losses ${\xi }_m$can be extracted accurately with reasonable computation resources.

In the MPM, the time-dependent signal is represented by a summation of a series of complex exponentials [14] as follow.

$R_m$is a real number and $s_m$is a complex number. By mapping Eq. (1) to Eq. (2), the mode frequency ${\omega }_m$and mode loss factor ${\xi }_m$can both be obtained by

The loss factor can be easily converted to attenuation coefficient ${\alpha }_m$by the relation ${\alpha }_m={\xi }_m/\left(d{\omega }_m/d{\beta }_m\right)$. It will be illustrated in the following sections that only limited number of field samples in time domain is required by MPM to obtain mode parameters accurately, which is very robust and more efficient compared with Fourier transform method.

4. Simulation results

4.1. Assessment and validation

In order to validate the mode solver presented in this paper, we calculated the real guided mode of a rectangular channel waveguide and the complex leaky mode of a ridge waveguide, respectively.

4.1.1. Guided mode analysis

The square channel waveguide structure is shown in Fig. 1 . The width ($w_1$) of the square core region is chosen as $2\mu m$ and the indices of the core region and the cladding region are 1.52 and 1.46, respectively. The number of PML layers is chosen to be 20 in our simulation.

 

Fig. 1 Schematic diagram of the square channel waveguide structure, $n_1=1.52$, $n_2=1.46$, $w_1=2\mu m$, PML Layers = 20.

Download Full Size | PDF

An initial optical field with the preset propagation constant β and the Gaussian distribution is launched at the position $\left(x_0,y_0\right)$ inside the core region in FDTD simulation. The excitation source is a Gaussian modulated sinusoidal waveform in time domain. The waveform in the frequency domain also follows a Gaussian distribution with a central frequency of${\omega }_0$, formulated by

When the steady-state is reached after the optical field propagates in the waveguide for a period of time, proper numbers of temporal samples of the optical field are taken at one position inside the core region in MPM analysis from which the resonant frequency for each mode can be extracted. The accuracy of extracted values by the full-vectorial solver is validated through comparisons made with benchmark results obtained by Goell’s model [23] in Fig. 3 . Note that for the given FDTD time steps and the MPM sampling numbers, the relatively larger error is observed near the cut-off. This is due to the reason that the effective width of each mode is increasing with the decrease of the normalized frequency V (horizontal label in Fig. 3). By using a larger computation window (in space), this error can be reduced.

 

Fig. 3 Guided modal analysis by Goell’s model and solver in this paper.

Download Full Size | PDF

Once the mode frequency ${\omega }_m$is extracted, the field distribution of each mode (${\Phi }_m$) can be readily obtained from the instantaneous field (${\Psi }_t$) by the formula

Figure 2 shows the calculated field pattern of the fundamental mode for the square channel waveguide depicted in Fig. 1.

 

Fig. 2 Fundamental mode profile of the square channel waveguide.

Download Full Size | PDF

4.1.2. Leaky mode analysis

In this section, the proposed method is applied to a deep-etched GaAs-AlGaAs optical waveguide structure [24] shown in Fig. 4 . The indices are 3.4519, 3.4434, 3.3955, 3.3675 for 5%, 6.5%, 15% and 20% AlGaAs material, respectively, and 3.4804 for GaAs cap and substrate. Note that the refractive index of the 5% AlGaAs core region is higher than those of the upper and lower cladding layers. The fundamental mode can be well confined in the rib. However high order modes are much lossier and can strongly leak through the lower 15% AlGaAs cladding layer and radiate into the substrate [24].

 

Fig. 4 Schematic diagram of a deep-etched GaAs-alGaAs optical waveguide structure.

Download Full Size | PDF

The TE20 mode frequency ${\omega }_m$ and attenuation coefficient ${\alpha }_m$ obtained by the proposed method with full-vectorial scheme are compared with those obtained from FD mode solver [4] in Table 1 . It also shows that results are well agreed with each other. The accuracy of the extracted results can be controlled by several factors such as the FDTD time steps and the number of sampling points for the MPM. For critical applications, the accuracy can be further enhanced with the cost of extra computation efforts. The TE20 mode profile of the ridge waveguide under investigation is calculated and shown in Fig. 5 . It is observed that the field is radiative towards to the substrate.

Table 1. Leaky modal analysis by FD mode solver and solver in this paper

View Table

 

Fig. 5 TE20 mode profile of the ridge waveguide.

Download Full Size | PDF

4.2. Error analysis and convergence test

The accuracy of the extracted parameters for the optical waveguide modal analysis depends on several factors, namely, the accuracy of the propagating field simulated by the FDTD and that of the MPM post processing. The former is related to the transverse mesh size whereas the latter on the FDTD running time (i.e., related to the sampling window in time) and MPM sample numbers. Assume that the FDTD mesh size is chosen in accordance with the numerical dispersion condition which is $\delta \le \lambda /12$. We may focus exclusively on the impact of the accuracy associated with the post data processor. Under this assumption, the convergence of the mode solver as functions of the FDTD running time and the MPM sample size is simulated and analyzed for the square channel waveguide structure described previously in section 4.1.1.

4.2.1. Convergence with FDTD time steps

The extracted wavelength of the fundamental mode is calculated with 480x480 mesh points and 7.810043rad/um propagation constant. The percentage error is defined with respect to the benchmark results obtained from the FD mode solver and shown in Fig. 6 as a function of FDTD time steps for a fixed number (1000) of MPM samples. It is observed that the error of the extracted parameter converges rapidly toward zero with the increase of FDTD time steps.

 

Fig. 6 Percentage error of extracted wavelength of fundamental mode versus different FDTD time steps.

Download Full Size | PDF

4.2.2. Convergence with MPM sample number

To study the convergence of this method with MPM sample numbers, the FDTD time steps is chosen to be 8000. The mode calculation is then carried out for the square channel waveguide with different numbers of MPM samples. The percentage error of extracted wavelength of the fundamental mode is shown in Fig. 7 . It is noted that the extracted parameter converges well with the increase of MPM sample numbers.

 

Fig. 7 Percentage error of extracted wavelength of fundamental mode versus MPM sample numbers.

Download Full Size | PDF

4.3. Performance of the mode solver on GPUs

In the mode solver based on FDTD, the most computational intensive work is for the FDTD simulation. Since the index of the waveguide is uniform along the propagation direction (z), the compact FDTD method is employed so that significant computation time can be saved. However, in the modal analysis of an optical waveguide with relatively large dimensions, a large number of meshes are required to ensure adequate accuracy. This leads to drastic increase in computation time and render this approach not applicable as a practical tool for mode calculation on CPUs of today’s desktop computers. In order to improve the computation efficiency, this method is implemented on Nvidia GTX 295 GPUs with parallel calculation algorithms. Comparisons of the memory and time consumptions are made in implementing this method on GPUs, CPU as well as the FD mode solver, respectively.

4.3.1. Memory consumption

When the FDTD mode solver is implemented on CPU and GPU, the memory requirements are virtually the same. The memory consumptions for the scalar mode calculation by the standard and the accelerated FDTD and the conventional FD mode solvers are estimated and shown in Fig. 8 . It is clearly shown that the memory consumption for the FDTD methods scale much slowly than that of the FD eigen solver. A saving of 94.6% memory can be realized for a problem of total mesh points of 480x480. In addition, we also examined the memory requirements of the scalar, semi-vector and full-vector mode calculations by the accelerated FDTD solvers showing in Fig. 9 . It is noted that there are moderate increases of memory consumption for the full and semi-vector modes (e.g., 20% and 35% increase, respectively) in comparison with the scalar mode. We conclude that the FDTD method is far more scalable than the FD scheme in terms of memory requirement and hence is capable of solving waveguide problems with large number of mesh points.

 

Fig. 8 Memory consumptions of the scalar mode calculation by the accelerated FDTD and the conventional FD mode solvers.

Download Full Size | PDF

 

Fig. 9 Memory consumptions of the scalar, semi-vector and full-vector mode calculations by the accelerated FDTD solvers.

Download Full Size | PDF

4.3.2. Time consumption

The computational efficiency is measured by the computation time required to achieve a desirable accuracy. In Fig. 10 , the computation time for scalar mode by the standard FDTD, accelerated FDTD and conventional FD solvers are plotted. Significant improvement in computation efficiency (e.g., reduction of computation time as much as 97% is realized for the scalar mode for a problem of total mesh points of 640x640) is achieved by the accelerated FDTD mode solver in comparison with the standard FDTD mode solver running on CPU of a standard desktop computer. It also shows that the computational efficiency of the accelerated FDTD method is in the same order of magnitude of the standard FD eigen mode solver. In Fig. 11 , the computation time for the scalar, semi-vector and full-vector modes by the accelerated FDTD mode solvers are plotted. It is noted that there are moderate increases of time consumption for the full and semi-vector modes in comparison with the scalar mode.

 

Fig. 10 Time consumption of the scalar mode calculation by the standard FDTD, accelerated FDTD and the conventional FD mode solvers.

Download Full Size | PDF

 

Fig. 11 Time consumption of the scalar, semi-vector and full-vector mode calculations by accelerated FDTD mode solvers on GPUs.

Download Full Size | PDF

5. Summary

It is proposed and demonstrated that the numerical mode solver based on the compact FDTD method in combination of the matrix pencil method (MPM) can be accelerated significantly by commercially available, low-cost and easy-to-implement high-performance computing techniques such as the CUDA algorithm on GPUs. In the meanwhile, the relative low memory requirements and simple numerical scheme inherent in the time explicit FDTD method is preserved. This idea is implemented and validated for calculations of scalar, semi-vector and full-vector modes and hence has general applicability to a wide range of optical waveguide problems. The advantages of this new approach become even more pronounced when the total number of meshes is large, in which the standard FDTD method is prohibitively time-consuming and the conventional FD solver is severely limited by its memory consumption.