A higher-order finite-difference time-domain (HO-FDTD) numerical method is proposed for the time-domain analysis of planar optical waveguide devices. The anisotropic perfectly matched layer (APML) absorbing boundary condition for the HO-FDTD scheme is implemented and the numerical dispersion of this scheme is studied. The numerical simulations for the parallel-slab directional coupler are presented and the computing results using this scheme are in highly accordance with analytical solutions. Compared with conventional FDTD method, this scheme can save considerable computational resource without sacrificing solution accuracy and especially could be applied in the accurate analysis of optical devices.
©2006 Optical Society of America
The finite-difference time-domain (FDTD) method  has been widely used to model and simulate optical devices [2, 3], since it can readily address complex geometric features and inhomogeneous materials. However, the conventional FDTD formulas exhibit only 2nd-order accuracy both in time and space domain. As a result, its numerical dispersion becomes a dominant factor in limiting applications of the FDTD to electrically small problems. For electrically large problems, the numerical dispersion inherent in the classical FDTD algorithm will introduce significant errors. Various techniques have been developed to improve the FDTD computation efficiency by developing the schemes that have extremely low numerical dispersion errors . One of them is the pseudospectral time-domain (PSTD) method  and another is the multiresolution time-domain (MRTD) method . By using the fast Fourier transform to present spatial derivatives, the PSTD method can reduce the spatial discretization to two steps per wavelength. However, its complex computation inherited from FFT doubles the compute memory requirement and the FFT iteration results in additional accumulating computation errors. By using orthonormal wavelet expansions, MRTD scheme can also reduce the spatial discretization to two steps per wavelength. However, the method is not very straightforward since the computations in MRTD are performed for wavelet expansion coefficients rather than the actual electromagnetic field quantities. In fact, the simplest method to reduce the dispersion error is to retain more terms in the Taylor series approximation of the spatial derivatives than it is done in 2nd-order central difference. Because of the additional large storage requirements when applying a similar approach to the temporal derivatives, the standard 2nd-order central differences are maintained in the approximation of the temporal derivatives [7, 8]. Fang was the first to present this approach in conjunction with solving Maxwell’s equations. He investigated the use of a 2nd-order accurate in time, 4th-order accurate in space FDTD algorithm, which we denote as the FDTD (2, 4) algorithm. The dispersion errors of higher order finite-difference time-domain (HO-FDTD) algorithms are better than those of MRTD algorithms that have equivalent spatial stencil sizes .
In this paper, the FDTD method with higher-order approximation has been introduced to the time-domain analysis of planar optical waveguide devices. The perfectly matched layers (PML) absorbing boundary condition for HO-FDTD schemes has been implemented and the numerical dispersion of this scheme is studied. The numerical simulations for the parallel-slab directional coupler are presented and the computing results using this scheme are in highly accordance with analytical solutions. Compared with conventional FDTD method, this scheme can save considerable computational resource without sacrificing solution accuracy and especially can be applied in computer-aided design and accurate analysis of optical devices.
2. Higher-order FDTD schemes
Maxwell’s equations for the two-dimensional TE wave can be written as
By applying the Taylor series expansion to a function f(x) and doing some manipulations, we can obtain (2M) th-order approximation for the first-order derivation of f(x)
where a(l)’s are the coefficient weights for (2M) th-order finite difference schemes and can be calculated by 
The coefficient weights a(l) for l≥0 can be calculated by Eq. (3), and those for l<0 can be obtained by a(-l)=-a(l-1). Now replacing the first-order time derivative in Eq. (1) with the 2nd-order central difference approximation and the first-order spatial derivatives in Eq. (1) with (2M) th-order central difference approximation, we can obtain the general (2, 2M) th-order FDTD formula of Eq. (1):
The coefficient weights a(l) for some of HO-FDTD schemes are listed in Table.1.
where v max is the maximum of propagation speed in the computational region.
Since the computational domain is limited in space by storage limitations, absorbing boundary condition should be implemented to effectively simulate open space and satisfy the radiation condition. According to the theory of APML [12,13], Maxwell’s equations for TE wave case can be written as
where sx and sy are the stretched coordinate variables in the x and y directions, they can be chosen as
where ση is the conductivity profile in the PML region along the η direction.
By using the frequency-time transform, we can derive their corresponding time-domain equations. Similar to the FDTD approach in Ref. , a two-step approach to derive the HOFDTD update equations in PML region is used
From the update equations above, it is noted that when updating the fields at the points which are close to the boundary the field components outside the domain are required, we can use the image theory to obtain these components . According to the image theory, the tangential components of the image E-field (parallel to the PMC mirror wall) and the normal components of the image H-field (normal to the PMC mirror wall) are always even symmetric about the original fields, while the normal components of the image E-field and the tangential components of the image H-field are always odd symmetric about the original fields. Using the above guidelines, we can express all fields in terms of those in the original region. Figure 1 shows the image theory for two-dimensional TE wave.
3. Numerical dispersion
The numerical dispersion for some of HO-FDTD schemes is investigated and compared with that of the standard FDTD.
By substituting a time-harmonic trial solution into the update equations, the dispersion relation for FDTD(2,2M) can be obtained
Where kx and ky are the x and y components of the numeric wave vector k, respectively, ω is the angular frequency and c is the speed of light.
Assuming Δx=Δy=Δs and ϕ define the direction of propagation, i.e. kx=k cos ϕ, ky=k sin ϕ, the dispersion relation can be expressed as
The numerical dispersion relation is a transcendental equation and the numeric phase velocity vp can be calculated by Newton’s method. The dispersion error ve can be evaluated by the following expression
In order to illustrate the relationship between the numerical dispersion error and the cell size, the numerical dispersion error is calculated with the propagation angle in the range of 0~90 degrees. Figure 2 shows the maximum of numerical dispersion error for different propagation angle versus cells per wavelength for Δt=0.1Δs/c and Δt=0.5Δs/c, respectively. From Fig. 2, we can see that when Δt=0.5Δs/c the maximum of dispersion error is about -6% for the standard FDTD, i.e. FDTD (2,2) and 1% for HO-FDTD. This means that those HO-FDTD schemes can exhibit the excellent dispersion characteristics when compared with the standard FDTD scheme. Figure 3 shows the numerical error versus the direction of propagation for Δs=λ/10 when Δt=0.1Δs/c and Δt=0.5Δs/c, respectively. It is found that the dispersion characteristic of FDTD (2, 4) is better than FDTD(2, 6) and FDTD(2, 8) when Δt=0.5Δs/c, but HO-FDTD schemes less depend on the direction of propagation.
4. Numerical results
Some HO-FDTD schemes are validated by the simulation of the wave propagation in the parallel-slab directional coupler made of two identical slabs (see Fig. 4). The refractive indexes of the coupler are n 1=2.42 and n 2=2. The dimension of the coupler is A=32µm and B=6µm. The thickness of each slab is D=0.4µm, and the gap between the slabs is S=0.3µm. The spatial distribution of excitation source is chosen to be similar to that of the TE0 mode and the wavelength of the excitation source in free space is 1.5 λ=µm. The whole analysis region is surrounded by ten-layer PMLs and the time-step is Δt=0.5Δs/c.
Figure 5 shows the field distribution of the parallel-slab directional coupler calculated by FDTD (2, 4) scheme with cell size of Δs=0.05µm. The power exchange between the two waveguides is clearly illustrated. By calculating the distance between the maximum and the minimum of the field amplitudes in the slab along the propagation direction, the coupling length is obtained as 10.775 Lc=µm. It is in good agreement with the analytical solution of 10.806µm . The relative error is about -0.287%.
Table 2 shows the results of coupling length for the coupler using different HO-FDTD schemes with three discretization levels: coarse-cell (Δs=0.1µm), normal-cell (Δs=0.05µm) and fine-cell (Δs=0.025µm). The calculated relative error reduces as the cell size becomes smaller, and is less than 1% when fine-cell is used. Compared with the standard FDTD scheme, these HO-FDTD schemes demonstrate obvious advantage in calculating efficiency and computing accuracy. As the dispersion characteristics of the HOFDTD scheme is better than the standard FDTD scheme, the HO-FDTD scheme need less grids and can save a lot of computer resource, while the same accuracy is attained. For example, the relative error calculated by FDTD (2, 4) scheme with coarse-cells is -1.906%, while the error calculated by FDTD (2, 2) scheme with normal cells is -2.6%, but both the computational time ratio and the number of cells ratio are about 1:4. Moreover, the increase in the number of order doesn’t provide obvious improvement in the accuracy of the solution, and an additional computational time is used. Thus the FDTD (2, 4) scheme is the best one when both calculating efficiency and computing accuracy are considered.
In this paper, the HO-FDTD schemes have been applied to the time-domain analysis of planar optical waveguide devices. The PML absorbing boundary condition for HO-FDTD schemes has been implemented and the numerical dispersion of this scheme is studied. The numerical simulations for the parallel-slab directional coupler are presented and the computing results using this scheme are in highly accordance with analytical solutions. In addition, the FDTD (2, 4) scheme is the best one when both calculating efficiency and computing accuracy are considered. Compared with conventional FDTD method, this higher-order scheme can save considerable computational resource without sacrificing solution accuracy and especially can be applied to the accurate analysis of arbitrary optical waveguide structures.
References and links
1. K. S. Yee, “Numerical solution of initial boundary value problem involving Maxwell’s equation in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1996).
2. S. Chu and S. K Chaudhuri, “A finite-difference time-domain method for the design and analysis of guidedwave optical structures,” IEEE, J. Lightwave Technol. 7, 2033–2038 (1989). [CrossRef]
3. W. Huang, C. L. Xu, and A. Goss, “A scalar finite-difference time-domain approach to guided-wave optics,” IEEE Photon. Technol. Lett. 3, 524–526 (1991). [CrossRef]
4. K. L. Shlager and J. B. Schneider, “Comparison of the dispersion properties of several low dispersion finite-difference time-domain algorithms,” IEEE Trans. Antennas and Propag. 51, 642–653 (2003). [CrossRef]
5. Q. H. Liu, “The pseudospectral time-domain (PSTD) method: a new algorithm for solution of Maxwell’s equations,” IEEE Antennas Propag. Soc. Int. Symp. 1, 122–125 (1997).
6. M. Krumpholz and L. P. B. Katehi. “MRTD New time-domain schemes based on multiresolution analysis.” IEEE Trans. Microwave Theory Tech. 44, 555–571 (1996). [CrossRef]
7. J. Fang, “Time Domain Finite Difference Computation for Maxwell’s Equations,” PhD thesis, (Univ. California, Berkeley, 1989).
8. C. W. Manry, S. L. Broschat, and J. B. Schneider, “Higher-order FDTD methods for large problems,” Appl. Comput. Electromagn. Soc. J. 10, 17–29 (1995).
9. K. L. Shlager and J. B. Schneider, “Comparison of the dispersion properties of higher order FDTD schemes and equivalent-sized MRTD schemes,” IEEE Trans. Antennas and Propag. 52, 1095–1104 (2004). [CrossRef]
10. M. L. Ghrist, “High order finite difference methods for wave equations,” M. S. dissertation, (University of Colorado, 1997).
11. E. M. Tentzeris, R. L. Robertson, and J. F. Harvey, “Stability and dispersion analysis of Battle -Lemarie-Based MRTD schemes,” IEEE Trans. Microwave Theory Tech. 47, 1004–1012 (1999). [CrossRef]
12. J. P. Berenger, “A perfectly match layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]
13. S. D. Gedney. “An anisotropic perfectly matched layer-absorbing medium for truncation of FDTD lattices,” IEEE Trans. Antennas and Propag. 44, 1630–1639 (1996). [CrossRef]
14. A. Yefet and P. G. Petropoulos, “A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell’s equations,” J. Comput. Phys. 168, 286–315 (2001). [CrossRef]
15. K. Okamoto, Fundamentals of optical waveguides, (Academic Press, 2000).