Finite-difference time-domain (FDTD) method

The finite-difference time-domain method discretizes space into a grid of cells and time into small steps, and updates the electric and magnetic fields in every cell from their neighbors at each step using the curl equations of Maxwell. A source injects a field, usually a short pulse, and monitors record the fields as they evolve; Fourier transforming the recorded fields gives transmission, reflection and field profiles over the whole pulse bandwidth from one run. Because it makes no assumption about the geometry beyond the mesh, FDTD handles scattering, strong reflection, radiation out of plane and resonances, which makes it the standard tool for grating couplers, photonic crystals, metasurfaces and plasmonic structures. The article on choosing a photonics simulation method compares it with the mode solver, eigenmode expansion and the beam propagation method, and gives a full three-dimensional cost estimate.

The Yee grid

Yee's 1966 scheme places the six field components at different points of each cell: electric field components on the cell edges and magnetic field components on the face centers, so that every component of $\mathbf{E}$ is surrounded by four components of $\mathbf{H}$ that form its curl, and vice versa. $\mathbf{E}$ and $\mathbf{H}$ are also offset by half a time step, and the algorithm leapfrogs: update $\mathbf{H}$ from the curl of $\mathbf{E}$, then $\mathbf{E}$ from the curl of $\mathbf{H}$. The central differences that result are second-order accurate in both space and time, and the staggering enforces the divergence conditions automatically. Materials enter through the permittivity assigned to each $\mathbf{E}$ point; dispersive materials such as metals or silicon over a wide band are modeled with Drude or Lorentz terms solved alongside the field equations.

Courant limit

An explicit time step is stable only if a wave cannot cross more than about one cell per step. On a uniform grid of spacing $\Delta x$ in $D$ dimensions the Courant condition is

$$\Delta t \le \frac{\Delta x}{c\sqrt{D}},$$

with $c$ the vacuum speed of light (the fastest wave anywhere in the domain). For $\Delta x$ = 10 nm, the limit is 33.4 as in 1D, 23.6 as in 2D and 19.3 as in 3D. Codes usually run slightly below the limit. Refining the mesh therefore increases the cost twice over: more cells, and more steps to cover the same physical time.

Mesh and numerical dispersion

The finite-difference grid makes the simulated phase velocity depend slightly on wavelength and direction (numerical dispersion). In one dimension at half the Courant limit, 10 cells per wavelength gives a phase velocity error of about 1.3% and 20 cells about 0.3%, falling with the square of the mesh. The wavelength that matters is the one inside the highest-index material: at 1550 nm in silicon ($n$ ≈ 3.48) it is 445 nm, so 20 cells per wavelength means a 22 nm mesh. Thin layers and gaps need a finer local mesh, and staircasing of curved or slanted boundaries on a rectangular grid adds error unless the code averages permittivity at interfaces.

Boundaries and run length

The domain is surrounded by absorbing boundaries, almost always the perfectly matched layer (PML) introduced by Berenger in 1994: a lossy layer, several to tens of cells thick, whose impedance matches the interior at every angle and frequency so that outgoing waves enter it without reflection and decay inside. Periodic and Bloch boundaries model gratings and photonic crystal unit cells; symmetry boundaries halve or quarter the domain when the structure and source allow it.

The run must last until the fields have decayed. Spectral resolution is set by the time window, $\Delta f \approx 1/T$: resolving 1 nm at 1550 nm (124.8 GHz) needs about 8.0 ps of simulated time. As a worked example, a 2D domain of 20 × 10 µm at a 10 nm mesh has two million cells, and 8.0 ps at the 2D Courant limit takes about 340,000 steps. A resonator with quality factor 10⁴ at 1550 nm has a photon lifetime of 8.2 ps, so its ring-down needs several times that; for high-Q devices the spectrum is usually obtained by fitting the decaying field or by simulating only the coupling region.

Pitfalls

• Unconverged mesh. A result is reliable only when it stops changing as the mesh is refined.
• PML too close. Evanescent fields that reach the PML are absorbed artificially, adding false loss, especially for high-index-contrast modes and leaky structures.
• Stopping too early. Truncating a slowly decaying field produces ripples and artificial broadening in the computed spectrum.

Common questions

When is FDTD the right choice?

When the structure changes along the propagation direction on the scale of the wavelength, scatters or radiates strongly, or needs a broadband response from a single run. For a waveguide cross-section uniform along its length a mode solver is far cheaper, and long, slowly varying structures are better handled by eigenmode expansion.

What is 2.5D FDTD?

A planar device is collapsed in the vertical direction with the effective index method, leaving a 2D FDTD problem that runs orders of magnitude faster at some loss of accuracy, especially for out-of-plane radiation.

References: K. S. Yee, IEEE Trans. Antennas Propag. 14, 302 (1966); J.-P. Berenger, J. Comput. Phys. 114, 185 (1994); A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005).