# Causal Materials in FDTD and EM Simulation *Dispersive materials in time-domain electromagnetic solvers must be implemented as K–K-consistent pole expansions; otherwise naive frequency-table lookups produce acausal wavefronts, late-time instability, and non-passive energy growth.* ## Why naive dispersion breaks FDTD The finite-difference time-domain (FDTD) method solves Maxwell's equations by stepping E and H fields on a staggered spatial grid at discrete times. In vacuum the update equation is ``` E^{n+1} = E^n + (Δt/ε₀) · ∇ × H^{n+1/2} H^{n+1/2} = H^{n−1/2} − (Δt/μ₀) · ∇ × E^n ``` and time-steps at roughly Δt ≈ Δx/(c√D) for a D-dimensional grid. The method is intrinsically a time-domain update; it has no frequency variable to plug a table of ε(ω) into. The **wrong** way to add a dispersive material is to multiply the frequency-domain FFT of E by ε(ω) each step, or to evaluate ε at the instantaneous carrier frequency of a narrowband pulse. Such schemes - violate causality: fields at t = 0⁺ depend on frequency components evaluated with future information; - are in general non-passive: small rounding errors can feed energy back into the fields, producing late-time instability that only shows up after millions of time-steps; - cause acausal wavefronts: the leading edge of a pulse can appear *before* t = z/c in the grid, a dead giveaway. The right fix is to represent ε(ω) as a sum of **simple causal poles** whose inverse Fourier transforms give simple time-domain convolutions, and then turn those convolutions into **auxiliary differential equations** or **recursive convolutions** that march in lockstep with the Yee update. ## Pole-fit representations Every passive dispersive material has ε(ω) that is analytic in the UHP (physics convention `e^{−iωt}`) with poles strictly in the lower half plane. The natural basis of K–K-consistent pole terms is a mix of Debye, Drude, and Lorentz contributions: ``` ε(ω) = ε_∞ + Σ_k Δε_k / (1 − iω τ_k) (Debye) + Σ_k −ω_p,k² / (ω² + iω γ_k) (Drude) + Σ_k A_k / (ω_k² − ω² − iγ_k ω) (Lorentz) ``` The combined form most commonly implemented in simulators is the Lorentz sum ``` ε(ω) = ε_∞ + Σ_k A_k / (ω_k² − ω² − iγ_k ω) ``` with each term contributing two complex-conjugate poles at ω = ±√(ω_k² − γ_k²/4) − iγ_k/2, strictly in the LHP as required by passivity. A Debye term is the degenerate case ω_k → 0 with one overdamped pole at ω = −i/τ_k; a Drude term is ω_k = 0 with a pole at ω = 0 and one at ω = −iγ_k. Because each summand is an analytic rational function with LHP poles only, the total ε(ω) **automatically satisfies K–K**. There is no consistency check to add — the pole expansion *is* the causality constraint, built into the ansatz. ## Fitting measured tables A measurement table of ε(ω) (e.g. from ellipsometry or from the K–K reflectance workflow of Chapter 23) is fitted to the pole sum above by minimising a weighted least-squares error over ω, subject to the **passivity constraints** ``` A_k, ω_k², γ_k ≥ 0, Im ε(ω) ≥ 0 for all real ω ≥ 0 ``` These constraints are not cosmetic. An unconstrained fit can in principle reproduce the table with negative γ_k and pass it to the solver; the resulting simulation will blow up exponentially as energy flows out of material poles into the fields. Vector-fitting algorithms (Gustavsen & Semlyen, 1999) enforce the LHP-pole constraint by construction and are the numerical workhorse for extracting pole models from measured data; commercial solvers HFSS and CST include proprietary equivalents. ## Auxiliary differential equation (ADE) update For each Lorentz term, introduce a polarisation vector P_k(t) obeying ``` d²P_k/dt² + γ_k dP_k/dt + ω_k² P_k = ε₀ A_k E(t) ``` which is a damped driven oscillator whose frequency-domain solution reproduces the k-th pole term. Discretise with central differences, yielding an update equation for P_k^{n+1} in terms of P_k^n, P_k^{n−1}, E^n. The electric displacement is then ``` D(t) = ε₀ ε_∞ E(t) + Σ_k P_k(t) ``` and the Yee update for E becomes ``` E^{n+1} = (ε_∞)^{−1} [D^{n+1}/ε₀ − Σ_k P_k^{n+1}/ε₀] ``` with D advanced by ∇×H as in vacuum. Storage cost is two extra vector fields per Lorentz pole per grid cell; cost per time-step is a few multiply-adds per pole. The resulting scheme is provably stable under a Courant condition that reduces to the vacuum one as the pole contributions vanish. ## Piecewise-linear recursive convolution (PLRC) An alternative formulation, due to Luebbers et al. (1990), writes the dispersive displacement as a running convolution ``` D(t) = ε₀ ε_∞ E(t) + ε₀ ∫₀^t χ(t − τ) E(τ) dτ ``` and approximates E(τ) as piecewise linear between samples. Each pole contributes a complex recursion of the form ``` Ψ_k^{n+1} = e^{−γ_k Δt} Ψ_k^n + α_k E^n + β_k E^{n+1} ``` where Ψ_k stores the running-convolution state. The update for E then folds Ψ_k^{n+1} into the usual Yee step. PLRC uses one complex state variable per pole (vs two real variables for ADE), roughly the same cost, and is popular in commercial codes. Both ADE and PLRC reproduce the *exact* K–K-consistent pole expansion in the limit Δt → 0, and both inherit the passivity of the underlying pole model. ## Passivity ↔ stability The stability criterion for the combined Yee + ADE/PLRC update is subtler than the vacuum Courant bound. A sufficient condition is - The pole expansion is passive (LHP poles, positive residues in an appropriate sense; formally, Im ε(ω) ≥ 0 for ω ≥ 0). - The Yee time-step satisfies the vacuum Courant bound. Pass either of these and the simulation is unconditionally stable out to arbitrarily late times. Violate passivity, even by floating-point rounding of a marginal fit, and the simulation grows exponentially — often after 10^5 to 10^7 steps, long after the user has stopped paying attention. This makes K–K-consistent pole fitting not just a niceness but an engineering necessity. ## Beyond FDTD The same machinery applies to any time-domain EM solver: FETD, DGTD, transmission-line matrix (TLM), and circuit-level time-domain SPICE macromodels for dispersive transmission lines. In every case the challenge is to approximate a measured ε(ω) (or impedance, or S-parameter block) by a **sum of causal poles** whose time-domain realisation is cheap and stable. The Kramers–Kronig relation is the built-in feature of the ansatz, not a cross-check applied afterwards. ## See also - [01_introduction.md](01_introduction.md) - [21_dielectric_spectroscopy.md](21_dielectric_spectroscopy.md) - [22_optical_constants.md](22_optical_constants.md) - [23_reflectance_spectroscopy.md](23_reflectance_spectroscopy.md) ## References - A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed., Artech House, 2005 — Chapter 9 (dispersive, nonlinear, and anisotropic media). - D. M. Sullivan, *Electromagnetic Simulation Using the FDTD Method*, 2nd ed., IEEE Press / Wiley, 2013. - R. J. Luebbers et al., "A frequency-dependent finite-difference time-domain formulation for dispersive materials", *IEEE Trans. Electromagn. Compat.* 32, 222 (1990). - B. Gustavsen and A. Semlyen, "Rational approximation of frequency domain responses by vector fitting", *IEEE Trans. Power Delivery* 14, 1052 (1999). - [Wikipedia: Finite-difference time-domain method](https://en.wikipedia.org/wiki/Finite-difference_time-domain_method)