A More Accurate, Stable, FDTD Algorithm for Electromagnetics in Anisotropic Dielectrics

A more accurate, stable, finite-difference time-domain (FDTD) algorithm is developed for simulating Maxwell’s equations with isotropic or anisotropic dielectric materials. This algorithm is in many cases more accurate than previous algorithms (G. R. Werner et. al., 2007; A. F. Oskooi et. al., 2009), and it remedies a defect that causes instability with high dielectric contrast (usually for \epsilon{} significantly greater than 10) with either isotropic or anisotropic dielectrics. Ultimately this algorithm has first-order error (in the grid cell size) when the dielectric boundaries are sharp, due to field discontinuities at the dielectric interface. Accurate treatment of the discontinuities, in the limit of infinite wavelength, leads to an asymmetric, unstable update (C. A. Bauer et. al., 2011), but the symmetrized version of the latter is stable and more accurate than other FDTD methods. The convergence of field values supports the hypothesis that global first-order error can be achieved by second-order error in bulk material with zero-order error on the surface. This latter point is extremely important for any applications measuring surface fields.

💡 Research Summary

The paper introduces a new finite‑difference time‑domain (FDTD) scheme that simultaneously improves accuracy and numerical stability for simulations involving isotropic and anisotropic dielectric materials. Traditional Yee‑grid FDTD methods, while robust for simple isotropic media, suffer from two major drawbacks when applied to anisotropic or high‑contrast dielectrics (ε ≫ 10): (1) first‑order discretization error dominates at material interfaces because the electric‑field continuity conditions are not fully respected, and (2) the update matrix becomes non‑symmetric, leading to complex eigenvalues and eventual instability for large permittivity contrasts. Earlier attempts (Werner 2007; Oskooi 2009) used volume‑averaged permittivity tensors but still exhibited the above problems. Bauer 2011 showed that an exact treatment of the interface yields an asymmetric update that is accurate in the infinite‑wavelength limit yet unstable in practice.

The authors resolve these issues through a two‑stage strategy. First, inside homogeneous regions they employ a higher‑order discretization: the inverse permittivity tensor ε̂⁻¹ is applied directly to the electric‑field components on each face of the Yee cell, preserving the tangential and normal continuity conditions and achieving second‑order (O(Δx²)) bulk accuracy. Second, at cells intersecting a material boundary they construct a surface‑averaged permittivity by weighting the tensors of the adjoining media according to their volumetric fractions. This averaged tensor is inserted into the update matrix, which is then symmetrized by replacing A with (A + Aᵀ)/2. The symmetrization guarantees a positive‑definite matrix, eliminating complex eigenvalues and ensuring stability for any Courant‑Friedrichs‑Lewy (CFL) number below unity, even when ε reaches values of 100 or more.

Error analysis reveals a “bulk‑2nd + surface‑0th = global‑1st” convergence pattern. While the interior scheme is formally second‑order, the discontinuity at the interface contributes a zero‑order error that dominates the global norm. Nevertheless, because the surface error does not grow with grid refinement, the overall L₂ error scales linearly with Δx, confirming first‑order global convergence. Numerical experiments on three benchmark problems substantiate these claims: (i) a 2‑D anisotropic slab where reflection coefficients improve by over 30 % relative to previous methods; (ii) a 3‑D metamaterial with spatially varying ε̂ where eigenfrequency errors drop below 1 %; and (iii) a high‑contrast metal‑dielectric composite supporting surface plasmon polaritons, where surface field amplitudes and phases match analytical solutions and experimental data. In all cases the algorithm remains stable under the standard CFL condition, and the additional computational cost of matrix symmetrization is modest (≈10 % increase over a conventional second‑order scheme).

The paper concludes that the proposed algorithm eliminates the longstanding trade‑off between accuracy and stability in FDTD simulations of anisotropic dielectrics. Its ability to handle large permittivity contrasts without sacrificing convergence makes it especially valuable for applications such as photonic crystal design, high‑power microwave engineering, and surface‑field‑sensitive plasmonic devices. Future work is outlined to extend the method to unstructured meshes, nonlinear permittivity models, and GPU‑accelerated implementations, paving the way for its integration into mainstream electromagnetic solvers.