The FDTD Method — Simulating Maxwell's Equations on a Grid

1. Why Solve Maxwell's Equations Directly?

Most introductory electromagnetism courses solve Maxwell's equations analytically: plane waves in free space, reflection at a single interface, resonance in a rectangular cavity. These closed-form solutions exist only for a small set of idealised, highly symmetric geometries. Real engineering problems — a phone antenna next to a curved metal chassis, light scattering off a irregular nanoparticle, a radar pulse hitting an aircraft — have no analytical solution.

The FDTD method sidesteps this entirely. Instead of solving for a closed-form field pattern, it discretises space into a grid of cells and time into small steps, then directly evaluates Maxwell's curl equations at every grid point using finite-difference approximations of the spatial and temporal derivatives. The field values are simply numbers stored in arrays, updated every time step according to their neighbours. No matrix inversion, no eigenvalue problem, no Green's function — just a local, explicit update rule applied over and over.

This makes FDTD extraordinarily general: any geometry that can be represented on a grid (arbitrarily complex conductors, dielectrics, anisotropic or dispersive materials) can be simulated with the same core algorithm. The price is computational cost — fine enough grids to resolve wavelength-scale features in 3D can require billions of grid points and millions of time steps for electrically large problems.

The starting point is Maxwell's two curl equations in a source-free, linear, isotropic medium:

These are a coupled pair of first-order partial differential equations: the time derivative of H depends on the spatial curl of E, and the time derivative of E depends on the spatial curl of H (plus a conduction-loss term). FDTD approximates every derivative — in time and in all three spatial directions — with a centred finite difference, then solves for the field at the next time step.

2. The Yee Grid: Staggering E and H in Space

Kane Yee's key insight, published in a 1966 IEEE paper, was to not store E and H at the same grid points. Instead, the electric and magnetic field components are staggered by half a grid cell in every direction — a scheme now called the Yee cell or Yee lattice.

In a 3D Yee cell, each E-field component sits at the centre of a cube edge, and each H-field component sits at the centre of a cube face, pointing normal to that face:

This staggering is not a cosmetic choice — it is exactly what makes centred spatial differences of E naturally line up on the H locations, and vice versa. A centred difference of E_y and E_z around an H_x point yields a second-order-accurate approximation to (∇ × E)_x precisely at the location where H_x lives. The result is that Faraday's law and Ampère's law are each satisfied to second-order accuracy without any interpolation.

The staggering also has a deep physical meaning: it automatically enforces the divergence conditions ∇·B = 0 and ∇·D = ρ at every grid point for all time, because the discrete curl operator on a Yee grid has the same null-space property as the continuous curl operator. This is why FDTD does not need a separate divergence-cleaning step, unlike many other PDE solvers.

In two dimensions (the common case for teaching and for many planar problems), the Yee grid reduces to two independent polarisations. The TM_z mode couples E_z, H_x, H_y; the TE_z mode couples H_z, E_x, E_y. A 2D simulation only needs to update three field components instead of six, which is exactly what a browser-based interactive FDTD demo typically implements.

3. The Update Equations: Leapfrog Time-Stepping

Time is discretised similarly: E-field values are computed at integer time steps (t = 0, Δt, 2Δt, ...) while H-field values are computed at half-integer time steps (t = Δt/2, 3Δt/2, ...). This is the leapfrog scheme: E and H alternate, each one updated using the most recent values of the other, which is already known from the previous half-step.

For the simplest possible case — a 1D simulation of a plane wave travelling along z, with E_x and H_y components only — the update equations reduce to:

Each update is completely explicit: the new field value is computed from a simple weighted sum of already-known values, with no system of equations to solve. This is what makes FDTD fast per time step and trivially parallel — every grid cell's update depends only on its immediate neighbours, making the method ideal for GPU acceleration and domain-decomposition on clusters.

One full simulation cycle is: update all H components from the current E field, then update all E components from the new H field, then repeat. Because H always uses the freshest E and vice versa, the scheme is second-order accurate in both space and time despite being fully explicit — a property that implicit methods usually need much more machinery to achieve.

4. The Courant Stability Limit

An explicit time-marching scheme is only useful if it is numerically stable — if small errors do not grow exponentially with each time step. For FDTD, stability requires that the time step Δt be small enough that information cannot propagate more than one grid cell per step, since the algorithm only exchanges information between adjacent Yee-grid neighbours.

This is the Courant–Friedrichs–Lewy (CFL) condition, and for FDTD it takes the specific form known as the Courant limit:

where c is the speed of light in the medium. In practice, engineers apply a safety margin — typically choosing Δt at 90–99% of the maximum allowed value (the Courant number S = c·Δt/Δx is kept just below 1 in 1D, or below 1/√dimensions in higher dimensions).

Violating the Courant limit is unforgiving: even a single time step above the threshold causes the highest-frequency grid modes to grow without bound, and within a handful of steps the entire simulation diverges into visible numerical noise or overflow. This coupling between spatial resolution and the maximum stable time step means that refining the grid to resolve fine geometric detail automatically forces smaller (more numerous) time steps — a major cost driver in fine-featured 3D FDTD models.

5. Numerical Dispersion and Grid Resolution

Real electromagnetic waves in vacuum are non-dispersive: every frequency travels at exactly c, regardless of wavelength or direction. The discretised FDTD grid, however, is not perfectly isotropic — a wave's numerical phase velocity depends slightly on its frequency and on the direction it travels relative to the grid axes. This artefact is called numerical dispersion.

For a 1D grid, the exact dispersion relation of the discrete update equations is:

which only reduces to the ideal ω = ck in the limit Δx, Δt → 0. For finite grid spacing, waves with shorter wavelength (higher k relative to the grid) travel slightly slower than the true speed of light, and waves travelling diagonally across the grid propagate at a slightly different speed than waves travelling axis-aligned. Over long propagation distances, this accumulates into visible pulse spreading and phase error.

The standard engineering rule of thumb to keep dispersion error acceptably small (typically under a few percent in phase) is:

where λ_min is the shortest wavelength present in the simulation (the highest frequency of interest, or the wavelength inside the highest-index material present). Because memory and runtime scale with the number of grid cells — and in 3D that means the cube of the linear resolution — the choice of cells-per-wavelength is usually the single biggest lever on simulation cost.

6. Injecting Sources: Hard, Soft, and TFSF

A simulation needs a way to introduce energy into the grid. The simplest approach, a hard source, simply overwrites a field value at one grid point every time step — for example forcing E_x[source] = sin(ωt). This is easy to implement but acts as a perfect scatterer: any wave that later returns to that point (a reflection from a boundary or object) is reflected back into the grid rather than passing through, which contaminates the results with spurious re-reflections.

A soft source instead adds the source term to the existing field value rather than overwriting it, which lets outgoing and returning waves superpose naturally at that point without the artificial scattering of a hard source. This is the preferred method for most continuous-wave and pulsed excitation.

For scattering problems — computing how a wave scatters off an object while cleanly separating the scattered field from the huge incident field — the standard technique is the total-field/scattered-field (TFSF) formulation. The grid is divided into an inner "total field" region (incident + scattered wave) and an outer "scattered field" region (scattered wave only), connected by a thin virtual boundary where correction terms inject exactly the incident field. This lets a plane wave illuminate a target cleanly from one direction while the outer region records only the object's scattered response — essential for radar cross-section (RCS) and far-field pattern calculations.

The temporal shape of the source also matters. A common choice is a Gaussian-modulated pulse, which has a smooth, band-limited frequency spectrum and avoids exciting spurious high-frequency grid noise that a discontinuous step-function source would generate:

7. Absorbing Boundaries: The Perfectly Matched Layer

Since a computer grid must be finite, but most electromagnetic problems represent radiation into unbounded free space, the grid's outer edges need to behave as if space continued forever — absorbing outgoing waves with no reflection. A naive truncation (simply stopping the grid) reflects waves back into the domain like a perfect mirror, corrupting the entire simulation.

The modern solution, introduced by Jean-Pierre Bérenger in 1994, is the Perfectly Matched Layer (PML): an artificial absorbing region placed around the simulation domain in which the wave impedance is mathematically matched to free space at every angle of incidence and every frequency, so that in the continuous (unmeshed) limit the reflection coefficient is exactly zero regardless of incidence angle. Inside the PML, the field amplitude decays exponentially via an artificial conductivity profile that is ramped smoothly from zero (at the PML's inner edge) to a maximum value (at the outer, perfectly-conducting termination):

The polynomial grading (rather than a sharp jump to full conductivity) is essential: an abrupt conductivity discontinuity itself causes reflection, so σ is ramped up gradually enough that the discretised grid sees a nearly continuous impedance match. A typical PML is 8–20 cells thick and, well tuned, achieves reflection coefficients below -60 to -80 dB — far below the level needed for accurate far-field and antenna simulations. Because a real PML is implemented on a discrete grid rather than a truly continuous medium, some small residual reflection always remains, but it is negligible for practical engineering purposes.

8. A Minimal 1D FDTD Implementation

The core of an FDTD solver is remarkably short. Here is a complete 1D simulation of a Gaussian pulse propagating in free space, including a simple lossy-layer absorbing boundary, written in JavaScript:

const N = 400;          // number of grid cells
const dx = 1e-3;       // spatial step (m)
const c0 = 3e8;        // speed of light (m/s)
const dt = dx / (2 * c0); // Courant limit, S = 0.5

let Ex = new Float64Array(N);
let Hy = new Float64Array(N);
const eps0 = 8.854e-12, mu0 = 1.2566e-6;

function step(t) {
  // update H from E (leapfrog half-step)
  for (let k = 0; k < N - 1; k++) {
    Hy[k] += (dt / (mu0 * dx)) * (Ex[k + 1] - Ex[k]);
  }
  // update E from H (leapfrog half-step)
  for (let k = 1; k < N; k++) {
    Ex[k] += (dt / (eps0 * dx)) * (Hy[k] - Hy[k - 1]);
  }
  // soft Gaussian source at cell 50
  const t0 = 40, tau = 12;
  Ex[50] += Math.exp(-((t - t0) ** 2) / (2 * tau * tau));
}

Running this loop for a few hundred time steps and plotting Ex against grid position at each frame shows the Gaussian pulse splitting into two counter-propagating wavefronts that travel outward at (very close to) the speed of light — exactly what Maxwell's equations predict, obtained without ever writing down a wave-equation solution by hand. Extending this same skeleton to 2D and 3D, adding materials with different ε and μ per cell, and adding a PML at the edges is precisely how production FDTD solvers such as MEEP, Lumerical FDTD, and openEMS are built, just with far more optimisation and many more field components.

9. Applications: Antennas, Photonics, and Radar

Because FDTD handles arbitrary geometry and arbitrary material properties (including frequency-dependent, anisotropic, and nonlinear media, via extensions of the basic update equations), it has become one of the two dominant methods in computational electromagnetics, alongside the frequency-domain Finite Element Method (FEM).

• Antenna design: predicting radiation patterns, input impedance, and bandwidth of phone antennas, patch antennas, and antenna arrays embedded in realistic device chassis.
• Photonics and nanophotonics: simulating light propagation through waveguides, ring resonators, photonic crystals, and plasmonic nanostructures where feature sizes are comparable to or smaller than the optical wavelength.
• Radar cross-section (RCS): computing how much radar energy an aircraft, ship, or vehicle scatters back toward the radar, using the TFSF technique described above.
• Bioelectromagnetics: modelling how radio-frequency energy from mobile phones or MRI coils is absorbed in tissue, expressed as the Specific Absorption Rate (SAR).
• Microwave and RF circuits: extracting S-parameters of filters, couplers, and connectors directly from time-domain field data via a Fourier transform of the recorded port voltages and currents.

FDTD's chief advantage over frequency-domain methods is that a single broadband time-domain simulation, followed by a Fourier transform of the recorded fields, yields the frequency response across an entire band of interest in one run — whereas a frequency-domain solver must be re-run separately at every frequency point. Its chief disadvantage is memory and runtime for electrically large 3D problems, since grid resolution must resolve the shortest wavelength present everywhere in the domain, even in regions far from any fine geometric detail — a limitation that has motivated non-uniform and adaptive grid variants.