## Abstract

The dynamic-thermal electron-quantum medium finite-difference time-domain (DTEQM-FDTD) method is used for efficient analysis of mode profile in elliptical microcavity. The resonance peak of the elliptical microcavity is studied by varying the length ratio. It is observed that at some length ratios, cavity mode is excited instead of whispering gallery mode. This depicts that mode profiles are length ratio dependent. Through the implementation of the DTEQM-FDTD on graphic processing unit (GPU), the simulation time is reduced by 300 times as compared to the CPU. This leads to an efficient optimization approach to design microcavity lasers for wide range of applications in photonic integrated circuits.

© 2013 OSA

## 1. Introduction

In the last twenty years, photonic devices had emerged strongly and played an important part in every aspect of our daily life. Some of these include laser sensing security identification and verification, communication using fibre optics and entertainment using light emitting diode (LED) and laser [1–4]. This huge development is owning to the enhanced interest in research and design of photonic devices. The first step in the research and design path is numerical simulation. Numerical simulation leads to preliminary prediction of functionality and reliability of a photonic device before manufacturing and marketing.

For numerical simulations of photonic devices and to predict interaction of lightwaves, Maxwell’s equations [5] are backbone. Maxwell equations represent fundamental combination of electric and magnetic fields, and predict the electromagnetic wave phenomena. Usually, Maxwell equations are simulated by either using the frequency and time domain method. Frequency domain methods are efficient for narrow bandwidth applications. Some of the frequency domain techniques are finite element method (FEM), finite different frequency domain (FDFD) and method of moment (MOM) [6–8]. Time domain methods are efficient for wide band applications. Finite-difference time-domain (FDTD) method is the most common numerical technique in time domain, with some unconditionally stable techniques such as locally one dimensional FDTD (LOD-FDTD) and alternating direction implicit FDTD (ADI-FDTD) [9–11]. FDTD has attracted much attention from computational scientist because it offers many advantages such as high accuracy, robustness and systematic computing. The FDTD method, also known as the Yee’s algorithm performs calculation by using Maxwell’s equations [12].

However, it is essential to divide photonic structures into smaller cells to achieve optimum accuracy in simulation results. The dimensions of each cell are often more than one-tenth of the functional wavelength so as to achieve optimum accuracy in the simulation results. This require smaller time step so as to achieve numerical stability and consistency. This strict requirement causes a major problem in simulation and calculation. The large number of cells and time step would require a larger computational memory to run the simulation especially for real time three dimensional structures. This heavy load results in long computational time. In addition, many current photonic devices consist of many complex shapes and media. Complex shapes and thin media require much smaller meshing size and this adds extra burden to the existing long computational time. Hence, the simulation time can increase tremendously. Additionally, large and extensive computational resources are required and may cause inconvenient to users. To solve these problems, the FDTD algorithm is implemented on GPU to accelerate the simulation speed. These implementations have been reported in refs [13–16]. In the reported work, the FDTD method is implemented on GPU for simple media definition. However, this work is different because the DTEQM-FDTD algorithm is implemented on GPU, which is used to simulate media by taking into account the complex electrons dynamics.

The objective of this paper is to study the different design cases of the microcavity. The results of the case studies are obtained by implementing complex media FDTD algorithm on graphic processor unit (GPU). The spectra and field distribution results are obtained by using a CPU or GPU and results are compared to check the accuracy. Then the spectra properties and field distribution of different microcavity designs are investigated. Subsequently, the extra-ordinary results of the different designed phenomena are discussed. Finally, the improvement in simulation speed because of the implementation of the complex media FDTD algorithm on CPU and GPU is demonstrated

## 2. Theoretical model: DTEQM-FDTD

The DTEQM-FDTD algorithm [17,18] is a computational model that is used for the simulation of optical devices with complex material properties. The model includes the electrodynamics and thermal process that is incorporated into Maxwell equations, and is modeled with FDTD. For electrodynamics process, the quantum structure of material is used to simulate the electron dynamics between different bands and energy levels. Several physical principles such as the Pauli Exclusion Principle and Fermi-Dirac thermalization [19,20] governed electrodynamics equations are incorporated in such a way that the computational model can be used to treat different types of media. These principles considered in the DTEQM-FDTD algorithm allow the modeling of active media with reasonable physical realism and accuracy. It is also sophisticated enough to encompass the essential atomic principles in the model. The presence of atomic parameters allows versatile applications of the model for a wide range of active photonic devices. In additions, due to these interesting features this model is also implemented for some interesting nanophotonics applications such as active plasmonics for waveguides, antennae, plasmonics source and light extraction [21–23].

In the DTEQM-FDTD method, electron dynamics and thermal processes contribute to the dielectric polarization of the electric field. The Maxwell equations with polarization term *P* are written as

*P*(where

_{i,r}*r*represent the

*x*,

*y*and

*z*components) that is different from the conventional electric field equations for the FDTD method. The parameter,

*i*represent the number of electron energy levels and

*M*is the maximum number of energy levels considered in the simulation. The energy level separation is represented by Δλ and is defined according to the number of energy level defined in the DTEQM-FDTD algorithm. For the magnetic field, the equations are identical to the equations used in conventional FDTD because the material is non-magnetic and is not affected by any external factors. However, if the magnetic properties of the material are affected by external magnetic field, then the magnetization density would need to be considered in the FDTD algorithm to model the effect of magnetic torque [24].

The polarization densities along the *x*, *y*, and *z* directions are defined in [17,18] and are written as

*μ*is the atomic dipole moment of electrons involved in the interband transition between the

_{ir}*i*th level in the conduction and valence bands for

*x*,

*y*and

*z*directions. The expression for the atomic dipole moment

*μ*, is given by ${\left|{\mu}_{ir}\right|}^{2}=\frac{3\pi \hslash {\epsilon}_{0}{c}^{3}}{{\left({\omega}_{ai}\right)}^{3}{\tau}_{i}}$ .

_{ir}The parameter, τ* _{i}* is the interband transition time between the

*i*th energy level. The interband frequency is denoted by ω

_{αi}. The parameter

*A*is the vector potential, where r =

_{r}*x*,

*y*or

*z*direction and

*γ*is atomic dipole dephrasing time at

_{i}*i*th energy level after excitation. The polarization density,

*P*is dependent on the atomic properties of the medium. It relates the quantized electromagnetic field and atomic dynamics.

_{ir}The electron and hole density at the *i*th energy level in the conduction band and valence band respectively is defined by *N _{Ci}* and

*N*. The rate equation [17] for the carrier densities which describes the intraband transition is given as

_{Vi}*represent the intraband transition time within the energy levels in conduction and valence band respectively. For the interband, the electron density that is being transferred from the valence to conduction band is given by*

_{C/V}_{pump}is the electron pumping rate to the active semiconductor media. The carrier density at the initial energy level is given by ${N}_{C/Vi}^{0}\left(r\right)$.

Although there are many transitions between and within the energy levels and band structure, the overall total number of carrier must be maintained at constant during the different simulation time to ensure charge conservation. The overall carrier density must satisfy the following conditions

The introduction of the 3D DTEQM-FDTD model to simulate active solid state semiconductor devices includes the complex electron dynamics in the band structure with the FDTD method. The inclusion of the electron dynamics Eqs. (2) – (6) would increase computation time. For efficient simulation, in the next section, the implementation of the method on GPU is given and discussed.

## 3. GPU implementation for DTEQM-FDTD

The graphics processing unit (GPU) has recently been the platform of choice for scientific simulations [24]. This is especially good for scientific simulation algorithm that requires plenty of parallelism. It is not difficult to understand the rationale behind this phenomenon. The current generation of GPU has hundreds of processing cores. For example, Nvidia Tesla C2050 has 14 multiprocessors of 32 processing cores each. Therefore, these 448 processing cores in the GPU provide the computation power that the scientific simulation can utilize. In order for the GPU to be utilized efficiently, thousands of processing threads have to be executed.

Apart from the numerous processing cores, memory access is another important factor in determining the computation speed. In order for the GPU to perform computation without much CPU intervention, each GPU is packaged on a board with its own memory, known as device memory. For example, Nvidia Tesla C2050 has 3 gigabyte of device memory. In order for the GPU to maximize the memory access throughput, it is important to make full use of coalescing. Coalescing is a mechanism to combine and reduce the memory access transaction if a set of consecutive neighboring threads is accessing a set of consecutive neighboring memory locations.

In the device memory, there is a section of memory that is classified as constant memory. For the Nvidia Tesla C2050, there is 64 kilobyte of constant memory. As the name implies, it is used for storing constant with the constraint of read-only access. To further reduce the memory access latency, there is also a 8 kilobyte cache on Nvidia Tesla C2050 per multiprocessor.

Within each multiprocessor, there is 48 kilobyte of on-chip memory on the Nvidia Tesla C2050, known as the shared memory. As mentioned earlier, there are a lot of threads running on each multiprocessor. However, the threads do not share their memory space. Therefore, the shared memory is a mechanism for these threads to access a common pool of memory. Furthermore, the shared memory is on the same chip, it is as fast as the register. It could also be programmed as a form of cache to speed up the memory access.

In the implementation of DTEQM-FDTD algorithm, all the characteristics of the GPU are taken into consideration and the implementation is designed to make full use of these processing cores and its corresponding memory access mechanism. Figure 1 show the flow chart of the implementation, which consists of three parts: pre-processing, GPU kernel execution and post-processing. Pre-processing, post-processing and the conditional check are executed on the CPU, while the kernels are executed on the GPU.

## 4. Design, simulation and results of DTEQM-FDTD implemented on GPU

In this section, an example of semiconductor elliptical microcavity is demonstrated. The accuracy and efficiency of the GPU implemented DTEQM-FDTD is also discussed. The results obtained with the DTEQM-FDTD method by using CPU and GPU are compared to confirm the correctness and accuracy. The optical spectrum and field distribution of the elliptical microcavity are compared to confirm the accuracy. Subsequently, more comprehensive study of the optical properties of the elliptical microcavity is presented by varying its dimensions.

#### 4.1 Elliptical microcavity parameters and setup

The simulation of three dimensional (3D) elliptical shaped microcavity is studied in this paper. The elliptical shaped microcavity has a semi-major and semi-minor length of 400 and 200 nm respectively. The thickness of the elliptical microcavity is 200 nm. The schematic layout of the elliptical microcavity is shown in Fig. 2
. The microcavity is made up of III-V material and is used as a semiconductor laser. Since the elliptical microcavity is made up of III-V material, the band structure is assumed to have direct bandgap and have parabolic energy distribution. This makes the modeling and simulation easier with DTEQM-FDTD method in this paper. The refractive index of the elliptical microcavity is 3.5. The elliptical microcavity is surrounded by air which has a refractive index of n_{air} = 1. An optical pulse is injected into the elliptical microcavity to generate resonance for lasing [17,18,20]. After 100000 simulation time step, the field distribution for the elliptical microcavity is taken.

The number of the electronic energy level *i*-states in the conduction and valence band is then defined. There are infinite numbers of energy level *i*-states in the conduction and valence band due to the large number of atoms in the semiconductor structure. For more practical modeling, the infinite energy levels are grouped together due to short dipole dephrasing time of the electron as compared to the decay lifetime of the photon. In this paper, 20 energy levels are defined in the conduction and valence band for sufficient modeling accuracy. The energy level separation is set to Δλ = 20 nm. The energy discrete levels (in electron volt units) are given by the Eq., Ei = 1.65 + 0.129**i* (eV), where *i* is the *i*th energy level and *i* = 1, 2, 3…..20, for 20 energy levels in the conduction and valence band. The value of 1.65 eV for E_{i} comes from the designated resonance wavelength at 752 nm. The interband transition wavelength between the conduction and valence bands is 850 nm. The effective mass of the electrons and holes in the valence and conduction band are given as 0.047 m_{e} and 0.36 m_{e}, where m_{e} is the free electron mass. This effective mass is obtained by considering the internal atomics binding and repulsion forces.

The transition time between conduction and valence band is defined as the interband transition time and is taken as 1 nanosecond (ns). The downward intraband transition time between the sub-energy levels (i.e. E_{c(i)} to E_{c(i-1)}) in the conduction and valence band is given as 1 picosecond (ps) and 100 femtosecond (fs) respectively. The upward intraband transition time can be obtained from the Fermi-Dirac thermalization steady state equations and is given as [25,26]

^{9}/s. For electron dynamic model in this paper, the spontanenous noise equations are not considered. Hence, a short optical pulse is manually added into the computational algorithm to excite the whispering gallery mode in the elliptical microcavity to instigate lasing. The excitation pulse is a 18 fs plane wave.

#### 4.2 Field distribution and resonance spectrum for elliptical microcavity

The normalized electric field distribution with respect to the incident electric field in the *z*-direction, *E _{z}* is shown in Fig. 3(a)
. The

*E*field distribution is observed at the middle of the three dimensional elliptical microcavity by slicing along the

_{z}*x*-

*y*plane at thickness of 100 nm from the top. Figures 3(b) and 3(c) shows the

*E*field distribution by slicing along the

_{z}*x-z*and

*y-z*direction at the major and minor length path. The whispering gallery mode [27] of the elliptical microcavity is the dominate mode as shown in Figs. 3(a)-3(c). Figure 4 shows the wavelength spectrum of the elliptical microcavity. The wavelength spectrum, represented by the bold line, is obtained by measuring the power at the edge of the elliptical microcavity. It shows that the resonance wavelength is at 752 nm. Hence, the field distribution pattern shown in Fig. 3(a) has wavelength of 752 nm.

Figures 3(a)-3(c) and the bold line spectrum in Fig. 4 are obtained by using the DTEQM-FDTD method on CPU platform. These results are used as a reference to check correctness and accuracy of the method when implemented on GPU. Figures 3(d)-3(f) and the scattered circle plot in Fig. 4 show the respective electric field distributions and spectrum that are obtained by using the GPU platform. It is observed that the *E _{z}* field distributions and the resonance spectrum match very well with the results that are obtained from CPU simulation. This supports the facts that GPU implementation is correct and accurate.

The *E _{z}* field distribution in Figs. 3(a)-3(c) and the spectrum in Fig. 4 are obtained at fixed semi-major and semi-minor lengths of 400 and 200 nm respectively. In the next section, the semi-major and semi-minor lengths are varied separately to investigate the effect on the resonance wavelength and electric field distribution in the elliptical microcavity.

#### 4.3 Case study 1: Varying the semi-minor length and its effects on wavelength shift and the electric field distribution

The length ratio, *L _{r}*, is defined as the ratio of the semi-major length to semi-minor length and is used to plot the results obtained from the simulation. The length ratio,

*L*is varied between 1.6 and 2.8 (i.e. the minor length is varied between 142 and 250). In Fig. 5(a) , the spectra for selected length ratio, L

_{r}_{r}is plotted on CPU and GPU. The plot in Fig. 5(b) shows the resonance peak wavelength for both CPU and GPU simulation. It is observed that the respective spectra and the resonance peak wavelength obtained from CPU and GPU agreed very well. In Fig. 5(a), the spectra for different L

_{r}1.739, 1.818 and 1.86 is shown for both CPU and GPU. The line and scatter symbol represent CPU and GPU implementation respectively and both results match very well.

The plot in Fig. 5(b) shows that the resonance wavelength is dependent on the length ratio, L_{r} of the elliptical microcavity. The decrease in length ratio results in the increase of size of the elliptical microcavity, which increases optical path length of the whispering gallery mode. Hence, the resonance wavelength increases because longer wavelength can be supported in the larger microcavity. It also satisfies the condition of the appearance of whispering gallery mode. In general, the resonance peak is red shifted with decrease in L_{r} (or increase in semi-minor length). However, there are two distinguish points in Figs. 5(a) and 5(b), at which the wavelength is blue shifted significantly. At the length ratio of 1.818 or minor length of 220 nm, the resonance wavelength drops to 568 nm. Similarly, at length ratio of 2.5, the resonance wavelength drops to 591 nm. In order to investigate this abnormal drop in resonance wavelength, we plot the only resonance peak wavelength at selected length ratios as shown in Fig. 5(b). It is shown clearly that the resonance peak wavelength steadily increases with decrease in the length ratio except for L_{r} values at 1.818 and 2.5.

To explain the blue shift phenomenon at these two length ratios, we plot the electric field distribution at each of the highest resonance peak for L_{r} = 1.86, 1.818 and 1.739. The electric field (*E _{z}*) distribution in the x-y plane is plotted for the length ratio of 1.818 as shown in Fig. 6(b)
, whereas for the neighboring length ratios 1.86 and 1.739 is also plotted in Figs. 6(a) and 6(c) respectively to demonstrate the difference in resonance wavelengths. Figures 6(a) and 6(c) show the whispering gallery mode in the elliptical microcavity, where field is mainly found on the outer part of the micorcavity and there is no field presence at the centre. However, the E

_{z}field distribution in Fig. 6(b) shows presence of more than one mode in the elliptical microcavity. It is observed that the field distribution consists of mixture of fundamental whispering gallery and cavity modes. The fields are present on the circumference as well as on the internal area of the elliptical microcavity. The whispering gallery mode in the elliptical microcavity is much weaker whereas the cavity mode is stronger. The cavity mode is the peak resonance mode as shown in Fig. 5(a). It has lower resonance wavelength compared to the whispering gallery mode. This is the reason why lower resonance wavelength is observed. If there is only whispering gallery mode, then the field distribution would be similar to that in Figs. 6(a) or 6(c). Since the cavity mode dominates in the elliptical microcavityat 568 nm wavelength, therefore, a change in the resonance wavelength is detected. This explains the blue shift at length ratio, L

_{r}= 1.818. The second blue-shift in the wavelength peak occurs when the length ratio is 2.5. The

*E*field distribution at length ratio of 2.5 and 2.353 is shown in Figs. 7(a) and 7(b) respectively. It is observed that in Fig. 7(a), the

_{z}*E*field distribution consists of mixture of cavity and whispering gallery modes. Again the cavity mode is a dominate mode and it causes a blue shift in the resonance peak.

_{z}By varying the semi-minor length of the elliptical microcavity, it is observed that the whispering gallery mode can be weakened. Hence, it is possible to control the mode profile in the elliptical microcavity. When the semi-minor length is at particular values, the condition for normal cavity mode is satisfied. Lightwaves reflect to and fro in the elliptical microcavity. This results in the presence of strong cavity mode. The cavity mode has higher amplitude compared to the whispering gallery mode in the elliptical microcavity. Hence, the cavity mode becomes the dominate mode. The resonance wavelength of the cavity mode is shorter. Therefore, a large blue shift occurs due to the dominant cavity mode in the elliptical microcavity. The weaker whispering gallery mode contributes much less to the resultant mode distribution in the elliptical microcavity. In the next section, the semi-major length is varied to investigate its effect on the field distribution pattern and resonance wavelength peak.

#### 4.4 Case study 2: Varying the semi-major length and its effects on wavelength shift and the electric field distribution

In the previous two sections, it is demonstrated that the implementation of the DTEQM-FDTD method on GPU provides results that are in very good agreement with CPU results. Therefore, all the results shown in this section are obtained by using GPU. The semi-major length is varied at length ratio, L_{r} from 1.75 to 2.5 (between 350 to 500 nm), while the semi-minor length is fixed at 200 nm. Figure 8(a)
shows the wavelength spectra for different length ratios by varying the semi-major length. Figure 8(b) shows the different resonance peaks for different length ratios. It is observed that the peak resonance wavelength varies with the length ratio in a “stair-step” manner. For example, the peak resonance wavelength increases from length ratio of 1.75 to 1.85. Then it drops to wavelength of 632 nm at L_{r} = 1.9. This trend continues until the length ratio of 2.2. After L_{r} = 2.2, the resonance wavelength trend increases up to length ratio of 2.5.

To explain the trend of the resonance peak wavelength, L_{r} = 2.05, 2.1 and 2.15 are selected to investigate the drop in resonance peak wavelength. The *E _{z}* field distribution pattern at selected L

_{r}is shown in Figs. 9(a) -9(c) respectively. The

*E*field distribution in Fig. 9(a) is for L

_{z}_{r}= 2.05 and the whispering gallery mode is present in the elliptical microcavity. However, Fig. 9(b) shows that the cavity mode is excited along the

*x*axis and oscillating in the elliptical microcavity. This cavity mode is the dominate mode, hence this explain the blue shift in resonance wavelength. As the length ratio increases to 2.15, the whispering gallery mode is excited again as shown in Fig. 9(c).

These results demonstrate that by varying the semi-major length of the elliptical microcavity, different peak resonance wavelengths can be obtained as compared to by varying the semi-minor length of the elliptical microcavity. Therefore, it is crucial that the correct semi-major and semi-minor lengths should be selected so that desired resonance modes and wavelengths can be excited in the elliptical microcavity.

## 5. Simulation speed improvement with GPU

The DTEQM-FDTD method is implemented on GPU (hardware accelerator) to improve the simulation speed. The DTEQM-FDTD consists of sophisticated electrons transition mechanism, and many parameters to encompass the essential physics. The executions of all these sophisticated physical equations require long simulation speed and large computational resources.

In Fig. 10(a) , simulation time comparison of the DTEQM-FDTD method is made with GPU and CPU. The GPU card used is Nvidia Tesla C2050, and the computer hardware specifcations used in the numerical experiments is Intel Core 2 Quad 3.2 GHz workstation with 4GB RAM. The simulation time using CPU is obviously much longer, whereas by using GPU the simulation time is reduced drastically. Figure 10(b) shows the relative improvement in the simulation time with respect to semi-minor length. The speed ratio is defined as the ratio of CPU time with respect to GPU time. It is observed that the speed up ratio increases at longer semi-minor length. The simulation time by using GPU is very short due to built in multiprocessing features. Hence, increasing the size of the elliptical microcavity only increases a very small amount of the GPU simulation time. However, this affects the simulation of the CPU by a considerable amount as shown in Fig. 10(a). This explains why the speed up ratio increases at longer semi-minor length. GPU implemented DTEQM-FDTD benefits longer semi-minor length more than shorter semi-minor length. It is observed that the simulation time by using GPU is reduced by approximately 300 times as compared to the by using CPU. This drastic reduction in simulation time makes the optimization and variation of different elliptical microcavity dimensions more convenient and meaningful. In addition, the shorter simulation time would allow the inclusion of other theories into the DTEQM-FDTD method for higher accuracy.

## 6. Conclusions

In conclusion, various cases are studied and discussed for different dimensions of the elliptical microcavity by implementing the DTEQM-FDTD method on GPU. The length ratios, L_{r} of elliptical microcavity were varied and some changes in the reasonance wavelength were observed. These include shifting of resonance peak wavelength at some particular length ratios, increasing the complexity of the mode profiles with different modes existing together, and mixing of the different modal field distribution in the elliptical microcavity. At few selected L_{r} values, the cavity mode was excited and becomes dominate mode. As a result, sharp blue shift in the resonance wavelength was observed. At other values of L_{r}, only whispering gallery mode was observed. GPU based DTEQM-FDTD approach drastically reduced the simulation time by approximately 300 times as compared to CPU. This achievement also results in more efficient and comprehensive study and understanding of the mode properties in the elliptical microcavity at different length ratios. This approach is believed to be useful for the simulation of large complex optical system with different complex materials. It also encourages the inclusion of more complex physical theories to the DTEQM-FDTD method for higher accuracy with minimal demand of computational resource. It is believed that the achievement in this paper is useful for simulation of nanophotonics devices such as surface plasmonic waveguides, optical antennae, nanolaser and exciton-plasmon interaction systems by changing the material conditions and formulations.

## Acknowledgment

The authors would like to express their sincerest gratitude to Prof Ho Seng Tiong of Northwestern University for his useful discussion. This work was supported by A*STAR Joint Council Organization and A*STAR-Mindef with Grant Nos. 12302FG012 and 1021740172.

## References and links

**1. **J. Rezac and A. Rosenberger, “Locking a microsphere whispering-gallery mode to a laser,” Opt. Express **8**(11), 605–610 (2001). [CrossRef]

**2. **M. Ohtsu, *Principles of Nanophotonics* (CRC Press/Taylor & Francis NW, 2008).

**3. **S. L. McCall, A. F. J. Levi, R. E. Slusher, S. J. Pearton, and R. A. Logan, “Whispering-gallery mode microdisk laser,” Appl. Phys. Lett. **60**(3), 289–291 (1992). [CrossRef]

**4. **R. E. Slusher, A. F. J. Levi, U. Mohideen, S. L. McCall, S. J. Pearton, and R. A. Logan, “Threshold characteristics of semiconductor microdisk laser,” Appl. Phys. Lett. **63**(10), 1310–1312 (1993). [CrossRef]

**5. **J. D. Jackson, *Classical Electrodynamic* (Wiley Press, 1999).

**6. **S. M. Hsu and H. C. Chang, “Full-vectorial finite element method based eigenvalue algorithm for the analysis of 2D photonic crystals with arbitrary 3D anisotropy,” Opt. Express **15**(24), 15797–15811 (2007). [CrossRef]

**7. **Z. H. Liu, E. K. Chua, and K. Y. See, “Accurate and efficient evaluation of method of moments matrix based on a generalized analytical approach,” PIERS **94**, 367–382 (2009). [CrossRef]

**8. **G. Strang and G. Fix, *An Analysis of The Finite Element Method* (Prentice Hall Press, 1973).

**9. **I. Ahmed, E. H. Khoo, and E. P. Li, “Development of the CPML for three-dimensional unconditionally stable LOD-FDTD method,” IEEE Trans. Antenn. Propag. **58**(3), 832–837 (2010). [CrossRef]

**10. **F. Zheng, Z. Chen, and J. Zhang, “A finite-difference time-domain method without the Courant stability conditions,” IEEE Microw. Guided Wave Lett. **9**(11), 441–443 (1999). [CrossRef]

**11. **I. Ahmed, E. K. Chua, E. P. Li, and Z. Chen, “Development of the three dimensional unconditionally stable LOD-FDTD method,” IEEE Trans. Antenn. Propag. **56**(11), 3596–3600 (2008). [CrossRef]

**12. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antenn. Propag. **14**(3), 302–307 (1966). [CrossRef]

**13. **S. E. Krakiwsky, L. E. Turner, and M. M. Okoniewski, “Acceleration of finite different time domain (FDTD) using graphics processor units (GPU),” IEEE Int. Microw. Sym. Digest **2**, 1033–1036 (2004).

**14. **S. Adam, J. Payne, and R. Boppana, “Finite different time domain (FDTD) simulations using graphics processor,” Proceedings of the Department of Defense High Performance Computing Modernization Program Users Group Conference, 334–338 (2007).

**15. **R. Sypek, A. Dziekonski, and M. Mrozowski, “How to render FDTD computations more effective using a graphics accelerator,” IEEE Trans. Magn. **45**(3), 1324–1327 (2009). [CrossRef]

**16. **K. H. Lee, I. Ahmed, R. S. M. Goh, E. H. Khoo, E. P. Li, and T. G. G. Hung, “Implementation of the FDTD method based on Lorentz-Drude model on GPU for plasmonics applications,” PIERS **116**, 441–456 (2011).

**17. **Y. Huang and S. T. Ho, “Computational model of solid-state, molecular, or atomic media for FDTD simulation based on a multi-level multi-electron system governed by Pauli exclusion and Fermi-Dirac thermalization with application to semiconductor photonics,” Opt. Express **14**(8), 3569–3587 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-14-8-3569. [CrossRef]

**18. **E. H. Khoo, S. T. Ho, I. Ahmed, E. P. Li, and Y. Huang, “Light energy extraction from the minor surface arc of an electrically pumped elliptical microcavity laser,” IEEE J. Quantum Electron. **46**(1), 128–136 (2010). [CrossRef]

**19. **R. K. Chang and A. J. Campillo, *Optical processes in Microcavities, Advanced series in Applied Physics* (World Scientific, Singapore 1996).

**20. **E. H. Khoo, I. Ahmed, and E. P. Li, “Enhancement of light energy extraction from elliptical microcavity using external magnetic field for switching applications,” Appl. Phys. Lett. **95**(12), 121104 (2009). [CrossRef]

**21. **I. Ahmed, E. H. Khoo, O. Kurniawan, and E. P. Li, “Modeling and simulation of plasmonic with FDTD method by using solid state and Lorentz -Drude dispersion model,” J. Opt. Soc. Am. B **28**(3), 352–359 (2011). [CrossRef]

**22. **E. H. Khoo, I. Ahmed, and E. P. Li, “Investigation of light energy extraction efficiency using surface plasmonics in electrically pumped semiconductor microcavity,” Proc. SPIE **7764**, 7764B (2010).

**23. **O. Kurniawan, I. Ahmed, and E. P. Li, “Generation of surface plasmon polariton using plasmonic resonant cavity based on microdisk laser,” IEEE Photon. J. **3**, 344–352 (2011).

**24. **R. Shams and P. Sadeghi, “On optimization of finite-difference time-domain (FDTD) computation on heterogeneous and GPU clusters,” J. Parallel Distrib. Comput. **71**(4), 584–593 (2011). [CrossRef]

**25. **S. Noda, M. Fujita, and T. Asano, “Spontaneous-emission control by photonic crystals and nanocavities,” Nat. Photonics **1**(8), 449–458 (2007). [CrossRef]

**26. **W. Fang, J. Y. Xu, A. Yamilov, H. Cao, Y. Ma, S. T. Ho, and G. S. Solomon, “Large enhancement of spontaneous emission rates of InAs quantum dots in GaAs microdisks,” Opt. Lett. **27**(11), 948–950 (2002). [CrossRef]

**27. **Y. Kogami, Y. Tomabechi, and K. Matsumura, “Resonance characteristic of whispering-gallery mode in an elliptic disk resonator,” IEEE Trans. Microw. Theory Tech. **44**(3), 473–475 (1996). [CrossRef]