We show how accuracy in finite-difference time-domain (FDTD) simulations [1] can be significantly improved for interfaces between anisotropic materials by careful selection of a subpixel smoothing scheme based on recent developments in perturbation theory [2]. This extends our previous work for interfaces between isotropic media [3], combined with a recently proposed scheme for stable FDTD in anisotropic media [4], and replaces our previous heuristic proposal for anisotropic-material interfaces [5]. Ordinarily, the presence of discontinuous material interfaces degrades the accuracy of FDTD to first-order $\left[O\left(\Delta x\right)\right]$ from the usual second-order $\left[O\left(\Delta x^2\right)\right]$ accuracy [6], 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 [4]. 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 [7], most prominently to design cloaking metamaterials [8]. 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 [9]).

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 [3]. Advances in perturbation theory have enabled us to extend this scheme to interfaces between anisotropic materials, initially for a plane-wave method [2]. Here, we adapt the technique to FDTD, combined with a recent FDTD scheme with improved stability for anisotropic media [4]. 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 $\mathbf{E}∕\mathbf{H}$.

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 ${\epsilon }^a$ and ${\epsilon }^b$, we showed that the proper smoothed permittivity (in this coordinate frame) at each point is [3]

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 $[i,j,k]=(i\Delta x,j\Delta y,k\Delta z)$, the $E_x$ and $D_x$ components are discretized at $[i+0.5,j,k]$, while $E_y∕D_y$ are at $[i,j+0.5,k]$ and $E_z∕D_z$ are at $[i,j,k+0.5]$ [1]. At each time step, $\mathbf{E}={\epsilon }^{-1}\mathbf{D}$ must be computed, but any off-diagonal parts of ε couple components stored at different locations. For example, a nonzero ${\left({\epsilon }_{xy}\right)}^{-1}$ means that the computation of $E_x$ requires $D_y$, but the value of $D_y$ is not available at the same grid point as $E_x$, as depicted in Fig. 1. One approach is to average the four adjacent $D_y$ values and use them in updating $E_x$, along with ${\left({\epsilon }_{xy}\right)}^{-1}$ at the $E_x$ point [3, 4]. This approach, however, is theoretically unstable and leads to divergences for a long simulation [4]. Instead, a modified technique was recently shown to satisfy a necessary condition for stability with Hermitian ε [4]: as depicted in Fig. 1, one first averages $D_y$ at $[i,j\pm 0.5,k]$ and multiplies by ${\left({\epsilon }_{xy}\right)}^{-1}$ at $[i,j,k]$, and then averages the two results at $[i,j,k]$ and $[i+1,j,k]$ to update $E_x$ at $[i+0.5,j,k]$. (Although [4] 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 [4].) We use this scheme here, and find that it greatly improves stability compared with the simpler scheme from our previous work [3]. The subpixel averaging is performed as follows. At the $E_x$ point $[i+0.5,j,k]$ [light gray dot (orange online) in Fig. 1], the smoothed $\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}$ is computed by Eq. (2), averaging over the pixel centered at that point. Then $\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}$ is inverted to obtain ${\left({\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}}^{-1}\right)}_{xx}$, which is stored at the $E_x$ point. The subpixel averaging $\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}$ is also performed for a pixel centered at the $[i,j,k]$ point [black dot (blue online)] halfway between two $D_y$ points [dark gray dot (red online)], and ${\left({\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}}^{-1}\right)}_{xy}$ is computed and stored at that point. This is similar for other components. (Note that the $\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}$ 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 ${\left({\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}}^{-1}\right)}_{xx}$ at $[i+0.5,j,k]$, ${\left({\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}}^{-1}\right)}_{yy}$ at $[i,j+0.5,k]$, ${\left({\stackrel{\mathrm{̃}}{\mathbf{\epsilon }}}^{-1}\right)}_{zz}$ at $[i,j,k+0.5]$, and all the off-diagonal components are at $[i,j,k]$—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 ${\epsilon }^{-1}$ necessary for that point. Each component of ${\epsilon }^{-1}$ 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 ${\epsilon }^a$ surrounded by ${\epsilon }^b$, a photonic crystal [13]. We choose ${\epsilon }^{a,b}$ to be random positive-definite symmetric matrices with random eigenvalues in the interval [1, 5] (1.45, 2.81, and 4.98) for ${\epsilon }^a$ and in [9, 13] (8.49, 8.78, and 11.52) for ${\epsilon }^b$. We compute the lowest ω for an arbitrary Bloch wave vector $\mathbf{k}=(0.4,0.2,0.3)2\pi ∕a$, 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 [14] to obtain the eigenfrequency ω, obtaining the relative error $|\omega -{\omega }_0|∕{\omega }_0$ by comparison with the “exact” ${\omega }_0$ from a plane-wave calculation [5] 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 [4].

We compare the new smoothing technique of Eq. (2) to the nonsmoothed case as well as to two simple smoothing techniques: using the mean $〈\epsilon 〉$ [15] and also the harmonic mean ${⟨{\epsilon }^{-1}⟩}^{-1}$. We do not compare with a previous heuristic that we had proposed without the benefit of perturbation theory [5], 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 [2], and first-order FDTD accuracy for that heuristic was also shown in [4] [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 [3], 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 $〈\epsilon 〉$ does better than the harmonic mean ${⟨{\epsilon }^{-1}⟩}^{-1}$, probably because it treats roughly two of the three field components correctly [3].) 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 [5] can sometimes be fortuitous [2], 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 [3]. 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 [3]. 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.

 

Fig. 1 Schematic 2D Yee FDTD discretization near a dielectric interface, showing the method [4] used to compute the part of ${\mathit{E}}_x$ that comes from ${\mathit{D}}_y$ and the locations where various ${\epsilon }^{-1}$ components are required.

Download Full Size | PDF

 

Fig. 2 Relative error Δω/ω for an eigenmode calculation with a square lattice (period a) of 2D anisotropic ellipsoids (right inset) versus spatial resolution (units of pixels per vacuum wavelength λ) for a variety of subpixel smoothing techniques. Straight lines for perfect linear (dashed) and perfect quadratic (solid) convergence are shown for reference. Most curves are for the first eigenvalue band (left inset shows $E_z$ in unit cell), with vacuum wavelength $\lambda =4.85a$. Hollow squares show new method for band 15 (middle inset), with $\lambda =1.7a$. Maximum resolution for all curves is $100\text{\hspace{0.17em}pixels}∕a$.

Download Full Size | PDF

 

Fig. 3 Relative error Δω/ω for an eigenmode calculation with a cubic lattice (period a) of 3D anisotropic ellipsoids (right inset) versus spatial resolution (units of pixels per vacuum wavelength λ), for a variety of subpixel smoothing techniques. Straight lines for perfect linear (dashed) and perfect quadratic (solid) convergence are shown for reference. Most curves are for the first eigenvalue band (left inset shows $E_x$ in $xy$ cross-section of unit cell), with vacuum wavelength $\lambda =5.15a$. Hollow squares show new method for band 13 (middle inset), with $\lambda =2.52a$. New method for bands 1 and 13 is shown for resolution up to $100\text{\hspace{0.17em}pixels}∕a$.

Download Full Size | PDF

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]  

3. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. Joannopoulos, S. Johnson, and G. Burr, Opt. Lett. 31, 2972 (2006). [CrossRef]   [PubMed]  

4. G. Werner and J. Cary, J. Comput. Phys. [CrossRef]  226, 1085 (2007).

5. S. G. Johnson and J. D. Joannopoulos, Opt. Express 8, 1732001. [CrossRef]   [PubMed]  

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]  

8. J. B. Pendry, D. Schurig, and D. R. Smith, Science 312, 1780 (2006). [CrossRef]   [PubMed]  

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]