Finite-difference time-domain methods suffer from reduced accuracy when discretizing discontinuous materials. We previously showed that accuracy can be significantly improved by using subpixel smoothing of the isotropic dielectric function, but only if the smoothing scheme is properly designed. Using recent developments in perturbation theory that were applied to spectral methods, we extend this idea to anisotropic media and demonstrate that the generalized smoothing consistently reduces the errors and even attains second-order convergence with resolution.
© 2009 Optical Society of America
We show how accuracy in finite-difference time-domain (FDTD) simulations  can be significantly improved for interfaces between anisotropic materials by careful selection of a subpixel smoothing scheme based on recent developments in perturbation theory . This extends our previous work for interfaces between isotropic media , combined with a recently proposed scheme for stable FDTD in anisotropic media , and replaces our previous heuristic proposal for anisotropic-material interfaces . Ordinarily, the presence of discontinuous material interfaces degrades the accuracy of FDTD to first-order from the usual second-order accuracy , but our work demonstrates how an appropriate choice of subpixel smoothing can both restore second-order asymptotic accuracy and give the lowest errors compared with competing schemes even at modest resolutions. Subpixel smoothing has an additional benefit: it allows the simulation to respond continuously to changes in the geometry, such as during optimization or parameter studies, rather than changing in discontinuous jumps as interfaces cross pixel boundaries. This technique additionally yields much smoother convergence of the error with resolution, which makes it easier to evaluate the accuracy and enables the possibility of extrapolation to gain another order of accuracy . The ability to handle anisotropic materials is becoming increasingly important via the use of anisotropic materials to represent arbitrary coordinate transformations in Maxwell’s equations , most prominently to design cloaking metamaterials . Our smoothing scheme requires preprocessing of the materials and does not otherwise modify the FDTD algorithm. It is therefore particularly simple to implement (free software is available ).
Our basic approach, as described previously [2, 3], is to smooth the structure to eliminate the discontinuity before discretizing, but because the smoothing itself changes the geometry we use first-order perturbation theory to select a smoothing with zero first-order effect. For isotropic materials, this approach made rigorous a smoothing scheme that had previously been proposed heuristically [10, 11, 12] and explained its second-order accuracy . Advances in perturbation theory have enabled us to extend this scheme to interfaces between anisotropic materials, initially for a plane-wave method . Here, we adapt the technique to FDTD, combined with a recent FDTD scheme with improved stability for anisotropic media . Although this Letter focuses on the case of anisotropic electric permittivity ε, exactly the same smoothing and discretization schemes apply to magnetic permeabilities μ owing to the equivalence in Maxwell’s equations under interchange of ε/μ and .
We define an interface-relative coordinate frame as in Fig. 1 , so that the first component “1” is the direction normal to the interface. Previously, for an interface between two isotropic materials and , we showed that the proper smoothed permittivity (in this coordinate frame) at each point is 2]:2] and we will not repeat it here, but we point out that Eq. (1) is now obtained as the special case for isotropic ε.
An additional difficulty for anistropic media occurs in FDTD: to accurately discretize the spatial derivatives, each field component is discretized on a different grid. In the standard Yee discretization for grid coordinates , the and components are discretized at , while are at and are at . At each time step, must be computed, but any off-diagonal parts of ε couple components stored at different locations. For example, a nonzero means that the computation of requires , but the value of is not available at the same grid point as , as depicted in Fig. 1. One approach is to average the four adjacent values and use them in updating , along with at the point [3, 4]. This approach, however, is theoretically unstable and leads to divergences for a long simulation . Instead, a modified technique was recently shown to satisfy a necessary condition for stability with Hermitian ε : as depicted in Fig. 1, one first averages at and multiplies by at , and then averages the two results at and to update at . (Although  derives no sufficient condition for stability with inhomogeneous media once the Yee time discretization is included, this method has been stable in all numerical experiments to date .) We use this scheme here, and find that it greatly improves stability compared with the simpler scheme from our previous work . The subpixel averaging is performed as follows. At the point [light gray dot (orange online) in Fig. 1], the smoothed is computed by Eq. (2), averaging over the pixel centered at that point. Then is inverted to obtain , which is stored at the point. The subpixel averaging is also performed for a pixel centered at the point [black dot (blue online)] halfway between two points [dark gray dot (red online)], and is computed and stored at that point. This is similar for other components. (Note that the tensor from (2) must be rotated from the interface-normal to Cartesian coordinates at each point.) Thus, for each Yee cell in three dimensions, the subpixel averaging is performed four times, obtaining at , at , at , and all the off-diagonal components are at —in other words, we apply the same averaging procedure in Eq. (2) to pixels centered around different points/corners in the Yee cell, and then for each point we store only the components of necessary for that point. Each component of need only be stored at most once per Yee cell, so no additional storage is required compared to other anisotropic FDTD schemes. After this smoothing, the anisotropic FDTD scheme proceeds without modification.
To illustrate the discretization error, we compute an eigenfrequency ω of a periodic (square in 2D or cubic in 3D, period a) lattice of dielectric ellipsoids made of surrounded by , a photonic crystal . We choose to be random positive-definite symmetric matrices with random eigenvalues in the interval [1, 5] (1.45, 2.81, and 4.98) for and in [9, 13] (8.49, 8.78, and 11.52) for . We compute the lowest ω for an arbitrary Bloch wave vector , giving wavelengths comparable with the feature sizes. In an FDTD simulation with Bloch-periodic boundaries and a Gaussian pulse source, we analyze the response with a filter-diagonalization method  to obtain the eigenfrequency ω, obtaining the relative error by comparison with the “exact” from a plane-wave calculation  at a high resolution. We looked at eigenvalue bands 1 and 15 (in 2D) or 1 and 13 (in 3D), where the higher band is clearly non-plane-wave-like (see inset fields), to counter suggestions that subpixel averaging may perform poorly for higher bands .
We compare the new smoothing technique of Eq. (2) to the nonsmoothed case as well as to two simple smoothing techniques: using the mean  and also the harmonic mean . We do not compare with a previous heuristic that we had proposed without the benefit of perturbation theory , since our previous paper already demonstrated that this heuristic (which does not yield zero first-order perturbation) is much less accurate than the new method , and first-order FDTD accuracy for that heuristic was also shown in  [who did not examine the isotropic case where Eq. (1) remains correct].
Results from 2D and 3D simulations are shown in Figs. 2 and 3 , respectively. In both cases, similar to our previous results for isotropic materials , the new smoothing algorithm has the lowest error, often by 1 order of magnitude or more, and it is the only technique that appears to give second-order accuracy in the limit of high resolution. (The simple mean does better than the harmonic mean , probably because it treats roughly two of the three field components correctly .) Similar accuracy is obtained for both lower and higher (non-plane-wave-like) bands at comparable resolutions per wavelength (although higher bands require greater absolute resolution per a, of course, because their wavelengths are smaller). As we have noted, apparent quadratic convergence obtained in a single structure  can sometimes be fortuitous , but we have confidence in these results (obtained now in multiple settings) because they are backed by a clear theory rather than an ad hoc heuristic.
Because the new smoothing scheme greatly improves the accuracy of FDTD simulation for anisotropic materials, without increasing the computational/storage cost (other than a one-time preprocessing step), it should be an attractive technique and subsumes our previous scheme . A remaining challenge is to accurately handle objects with sharp corners, where the resulting field singularities are known to degrade the accuracy to between first- and second-order once the smoothing eliminates the first-order error . We are hopeful that an accurate smoothing can be developed for corners once the corresponding perturbation theory is derived.
This work was supported in part by the MRSEC program of the National Science Foundation (NSF) under award DMR-0819762 and the Army Research Office through the Institute for Soldier Nanotechnologies under contract DAAD-19-02-D0002.
1. A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method, 3rd ed. (Artech, 2005),
2. C. Kottke, A. Farjadpour, and S. Johnson, Phys. Rev. E 77, 036611 (2008). [CrossRef]
4. G. Werner and J. Cary, J. Comput. Phys. [CrossRef] 226, 1085 (2007).
6. A. Ditkowski, K. Dridi, and J. S. Hesthaven, J. Comp. Physiol. 170, 39 (2001).
7. A. J. Ward and J. B. Pendry, J. Mod. Opt. 43, 773 (1996). [CrossRef]
9. Meep FDTD, http://ab-initio.mit.edu/meep.
10. R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, Phys. Rev. B 48, 8434 (1993). [CrossRef]
11. S. G. Johnson, R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, Phys. Rev. B 55, 15942 (1997).
12. J.-Y. Lee and N.-H. Myung, Microwave Opt. Technol. Lett. 23, 245 (1999). [CrossRef]
13. J. D. Joannopoulos, S. G. Johnson, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton U. Press, 2008).
14. V. A. Mandelshtam and H. S. Taylor, J. Chem. Phys. 107, 6756 (1997). [CrossRef]
15. S. Dey and R. Mittra, IEEE Trans. Microwave Theory Tech. 47, 1737 (1999). [CrossRef]