Yee FDTD Scheme Overview

• The Yee FDTD scheme is a numerical method that discretizes Maxwell's equations on staggered spatial and temporal grids, ensuring discrete charge-conservation and divergence constraint satisfaction.
• It employs a leap-frog update rule under a strict CFL condition to maintain stability and accurately simulate electromagnetic wave propagation despite numerical dispersion challenges.
• Advanced extensions include interface-corrected updates, metasurface modeling via GSTCs, and SBP-SAT subgridding, which enhance accuracy and efficiency in complex multi-physics applications.

The Yee finite-difference time-domain (FDTD) scheme is a cornerstone computational technique for solving the time-dependent Maxwell equations on regular grids. Introduced in 1966, it forms the basis for numerous modern solvers addressing electromagnetics, photonics, and plasma physics. The Yee scheme is characterized by its use of staggered (collocated) spatial and temporal grids for electric and magnetic field components, ensuring centered approximations of spatial derivatives and inherent preservation of discrete charge-conservation and divergence constraints. Over the past decades, the scheme has undergone significant refinements and extensions, supporting stable subgridding, explicit treatment of thin metasurfaces via generalized sheet transition conditions, interface-corrected updates, and dispersion-optimized formulations. This article reviews the numerical structure of the Yee scheme, its stability properties, convergence and error analysis, advanced interface and metasurface modeling, and extensions for high-order and multi-physics applications.

1. Discretization Structure and Update Rules

The Yee scheme discretizes the time-dependent Maxwell equations by collocating electric field components on the edges of cubic cells and magnetic field components on the cell faces. Time-stepping is leap-frogged: electric fields are known at integer time levels $n\Delta t$, while magnetic fields are advanced at $(n+\tfrac12)\Delta t$.

In three dimensions, the core update equations at grid location $(i,j,k)$ are: $$E_x^{n+1/2}(i,j,k) = E_x^{n-1/2}(i,j,k) + \frac{\Delta t}{\epsilon_0} \left[ \frac{H_z^n(i,j,k) - H_z^n(i,j-1,k)}{\Delta y} - \frac{H_y^n(i,j,k) - H_y^n(i,j,k-1)}{\Delta z} \right]$$ Magnetic field updates are similar, e.g.

$$H_y^{n+1/2}(i,j,k) = H_y^{n-1/2}(i,j,k) - \frac{\Delta t}{\mu_0} \left[ \frac{E_z^n(i,j,k) - E_z^n(i,j,k-1)}{\Delta x} - \frac{E_x^n(i,j,k) - E_x^n(i-1,j,k)}{\Delta z} \right]$$

This highly regular, staggered stencil ensures divergence-free updates provided initial data are compatible, and supports second-order accuracy in both space and time (Smy et al., 2017, Usuki, 2013, Blinne et al., 2017).

2. Stability and Numerical Dispersion

The classical Yee FDTD scheme is conditionally stable under a Courant–Friedrichs–Lewy (CFL) condition: $$c\,\Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}}}$$ Here, $c$ is the speed of light and $\Delta x$, $\Delta y$, $\Delta z$ are the mesh spacings (Smy et al., 2017, Blinne et al., 2017). For explicit time-stepping, one typically uses $\Delta t \approx \Delta / (c\sqrt{D})$, with $D$ being the dimensionality.

A key feature is numerical dispersion: discrete phase and group velocities deviate from the physical value $c$, especially for high wavenumbers (near the grid's Nyquist limit). The discrete dispersion relation is

$$\left( \frac{\sin(\frac{\omega \Delta t}{2})}{c \Delta t} \right)^2 = \sum_{\alpha = x,y,z} \left( \frac{\sin(\frac{k_\alpha \Delta x_\alpha}{2})}{\Delta x_\alpha} \right)^2$$

This leads to $v_\text{ph}, v_g < c$ for typical stencils, affecting wave propagation and leading to artifacts such as numerical Cherenkov radiation in particle-in-cell (PIC) simulations (Blinne et al., 2017). Optimized stencils with extended coefficients can significantly reduce these effects (Blinne et al., 2017).

3. Accuracy at Material Interfaces

The staggered nature of Yee grids requires careful analysis of field continuity at material interfaces. For normally incident harmonic plane waves, FDTD interface models spread the discontinuity over a transition layer of thickness $\Delta_x$, leading to discrete Fresnel reflection and transmission coefficients (Makarov et al., 16 Dec 2025). The FDTD analogs are: $$\widetilde t = \frac{2\eta_2\cos(\frac{\widetilde k_1\Delta_x}{2})}{\eta_2\cos(\frac{\widetilde k_2\Delta_x}{2}) + \eta_1\cos(\frac{\widetilde k_1\Delta_x}{2})}$$

$$\widetilde r = \frac{\eta_2\cos(\frac{\widetilde k_2\Delta_x}{2}) - \eta_1\cos(\frac{\widetilde k_1\Delta_x}{2})}{\eta_2\cos(\frac{\widetilde k_2\Delta_x}{2}) + \eta_1\cos(\frac{\widetilde k_1\Delta_x}{2})}$$

with

$$\sin(\tfrac{\widetilde k_i\Delta_x}{2}) = \frac{n_{r,i}}{S} \sin(\tfrac{\omega\Delta_t}{2})$$

where $n_{r,i}$ is the refractive index and $S = c\Delta t / \Delta_x$ is the Courant number. The leading-order error is $O((k\Delta_x)^2)$, with sign and magnitude determined by the relative impedances of the two media (Makarov et al., 16 Dec 2025). Practical error estimates predict reflection-coefficient errors are always of one sign (over- or underestimation depending on $\eta_1 \gtrless \eta_2$) and that the FDTD “transition layer” effect is fundamental unless advanced interface-corrected schemes are used.

4. Advanced Interface and Metasurface Modeling

Extensions of the Yee-FDTD scheme support physically rigorous modeling of zero-thickness interfaces and complex surfaces. Generalized sheet transition conditions (GSTCs) define discontinuities for tangential fields at metasurfaces in terms of electric and magnetic polarization densities (Smy et al., 2017, Vahabzadeh et al., 2017): $$\hat n\times(\mathbf H_2-\mathbf H_1) = j\omega\,\mathbf P_\parallel + \nabla_\parallel M_n$$

$$\hat n\times(\mathbf E_2-\mathbf E_1) = -j\omega\mu_0\,\mathbf M_\parallel - \nabla_\parallel P_n$$

GSTCs can be embedded in the FDTD update through symmetric, asymmetric, or “tight asymmetric” cell configurations, each with specific discretizations of the field jumps and associated auxiliary differential equations for broadband Lorentzian dispersion. The tight asymmetric cell (TAC) achieves fastest convergence and minimal spurious reflection (Smy et al., 2017). Auxiliary differential equation (ADE) methods are required for dispersive surfaces (Smy et al., 2017).

For stepwise interfaces, high-order accuracy and oscillation-free fields are realized via the correction function method (CFM), which solves local PDEs on patches straddling the interface and injects high-order correction terms to the Yee updates. The method is compatible with both standard and high-order FDTD schemes and achieves up to fourth-order convergence in $L^2$ norm for smooth, complex interfaces (Law et al., 2021).

5. Subgridding, Stability, and High-Order Extensions

For multiscale problems, energy-stable subgridding within the FDTD framework is accomplished using summation-by-parts (SBP) finite difference operators and simultaneous approximation terms (SATs). These techniques enforce boundary and interface conditions weakly and yield a discrete energy estimate, ensuring provable long-term stability even with nonconforming mesh blocks (Cheng et al., 2021, Wang et al., 2022). The interpolation matrices used at block interfaces are constructed to maintain SBP properties and ensure non-dissipative coupling.

Standard and high-order extensions retain the leap-frog update for the bulk grid, while modifications are confined to interface or boundary-adjacent nodes. When SBP-SAT operators are employed, accuracy and stability are preserved for arbitrarily refined subgrids (Cheng et al., 2021, Wang et al., 2022).

6. Frequency-Domain and S-Matrix Formulations

The Yee spatial discretization is also employed in frequency-domain transfer-matrix and S-matrix formulations. Here, Maxwell’s equations are discretized in space and solved (typically slice-by-slice) for modal scattering coefficients. The method guarantees unitarity (exact power conservation) to machine precision and quantifies residuals in eigenmode orthogonality and S-matrix unitarity below $10^{-8}$ (Usuki, 2013). This frequency-domain approach avoids stability issues inherent to direct inversion or explicit time-stepping and is robust to high contrast and complex geometries.

7. Recommended Practices and Limitations

For accurate simulation of wave propagation and interfaces using the Yee scheme:

• Employ at least 30–50 points per wavelength for interface accuracy better than 5% in engineered systems, and 200–400 points per wavelength for optical-frequency metasurfaces (Smy et al., 2017, Makarov et al., 16 Dec 2025).
• Use a causal Lorentzian (or Debye/Drude) dispersion model for frequency-dependent susceptibilities to avoid instability and nonphysical group delay (Smy et al., 2017).
• CFL-limited time-stepping should be based on the smallest grid cell, and the stability constraint remains unchanged for embedded GSTC or corrected-interface schemes (Smy et al., 2017, Law et al., 2021).
• For metasurfaces and sharp interfaces, prefer specialized GSTC or CFM-corrected FDTD rather than naive thin-slab approximation. Avoid frequency-independent surface susceptibilities in time-domain ADE models.

The Yee FDTD family combines regular structure with extensibility for accurate, stable, and efficient solving of Maxwell’s equations across a wide array of electromagnetic, photonic, and multi-physics scenarios. Advances in interface treatment and dispersion optimization continue to refine its applicability to cutting-edge computational research (Smy et al., 2017, Blinne et al., 2017, Makarov et al., 16 Dec 2025, Law et al., 2021, Wang et al., 2022).