2D FDTD Simulator: Yee Grid, TMz Waves and PML Boundary

The original implementation here follows the TMz Yee staggering explained in EMPossible, 2D FDTD formulation (attached s3), and McMaster ECE757 (s2, pp. 1–10). With vacuum outside the block, the fastest speed is $c_0$, and

The spatial grid has 80×80 Ez nodes; outer PEC nodes are 79 cell widths apart. The source is at cell coordinates (39.5,39.5), shared equally between the four nearest electric nodes. Hx sits at (x,y+Δy/2), Hy at (x+Δx/2,y), and H is advanced half a time step ahead of E.

For split-field PML, write $E_z=E_{zx}+E_{zy}$. The differential equations are

In the vacuum PML, $\sigma^*/\mu_0=\sigma/\epsilon_0$. Centering the loss term gives $q=\sigma\Delta t/(2\epsilon_0)$, $a=(1-q)/(1+q)$ and $b=1/(1+q)$. Conductivity is sampled at the actual E and H locations; x and y losses both act in the corners. Inside the dielectric, losses vanish and electric curl increments are divided by εr. The solver stores Ez/Edrive and η₀H/Edrive. This is Berenger split-field PML, not a uniform lossy sponge and not CPML.

For distance $s$ into a layer, $\sigma(s)=\sigma_{\max}(s/d)^m$. At normal incidence, a PEC-terminated continuum layer has round-trip amplitude

Therefore $\sigma_{\max}=(m+1)\ln[1/R(0)]/(2\eta_0d)$, separately for $d_x=p\Delta x$ and $d_y=p\Delta y$. The thickness chart holds the current x-normal σmax fixed and compares m with m+1. Amplitude dB uses 20 log₁₀; energy dB uses 10 log₁₀.

Source correction: Wang et al., arXiv:2509.20801, attached s5 p. 4 Eq. (13), prints ln[R(0)] with a positive prefactor. That would give negative loss for 0<R<1; the attenuation integral requires −ln R = ln(1/R). Its Eqs. (11)–(12) also include /Δ. Here σmax is explicitly in S/m and is applied only once, without an additional division by cell size. Corbino, arXiv:2603.19056v2, s6 p. 8 Eq. (29), supplies a separate conductivity-rule cross-check; the paper’s mimetic damping code is not copied or treated as a split-field benchmark.

At every step the soft increment $g_n=\exp[-((n-4\tau_s)/\tau_s)^2]$ is added equally to the two electric splits at the source. It stops after n=8τs. The setup is symmetric about y = 39.5. The solver advances the lower half-domain with even Ez and Hy and odd Hx. For vacuum it also uses x symmetry, advancing one quarter. The block breaks x symmetry, so that case advances the full width. Mirroring reconstructs the full 80×80 field exactly for the specified source and geometry. Energy sums are scaled to represent the full domain. Colour frames are precomputed every eight steps; energy is sampled every four steps. For the energy diagnostic, average E at n and n−1 and combine it with H at n−½, using $\tfrac12\Delta x\Delta y[\sum\epsilon |E|^2+\sum\mu_0|H|^2]$ per unit z length. Only nodes and edges between the inner PML interfaces are summed; in PEC mode the entire domain is used. This simple quadrature is a diagnostic, not an exact discrete conserved energy. The residual is final sampled energy divided by the largest sampled energy of the same run.

Independent check of absorption

Boito et al. (2020), attached s10 p. 6 Eqs. (43)–(44), provide the retarded 2D scalar-wave Green function. With normalization $(\partial_t^2-c_0^2\nabla^2)K=\delta$, its kernel is

It is applied to the TMz scalar wave equation here; the reference uses no finite boundary or PML equations. Each soft increment $g_k$ is represented as constant Ampere drive $g_k/\Delta t$ during $[(k-1)\Delta t,k\Delta t]$. Its time derivative gives two impulses. After the pulse has passed, their exact convolution is

where $\rho_\ell^2=[(p_x-\ell_x)/C_x]^2+[(p_y-\ell_y)/C_y]^2$ and the four source nodes are (39,39), (40,39), (39,40), (40,40). This evaluates a published analytical solution for the tool’s specified excitation; the numbers are not claimed as a published numerical table.

The frozen benchmark uses vacuum, Δx=Δy=10 mm, safety 0.95, τs=8, 20-cell cubic PML and target R=10⁻⁶. At probe nodes (39,39) and (49,39), the exact normalized fields are respectively −0.000304718749 and −0.000310405656 at step 160, and −0.000176765461 and −0.000178660114 at step 200. Both times allow an undamped return from the outer wall, so the comparison tests absorption. Each case allows 2% relative error for the finite grid and discrete boundary. A read-only mutation removing both PML loss factors makes these late-field comparisons fail, while leaving the early stencil check intact. This verifies absorption at the benchmark’s accuracy; it does not certify the continuum reflection target as a measured discrete reflection.

TMz only, invariant along z. Fields begin at zero. Vacuum has c₀ = 299792458 m/s and η₀ = 376.730313668 Ω; ε₀ is derived consistently as 1/(η₀c₀). The soft source has a fixed peak increment of 1 V/m, so all displayed energies scale with its square.

The optional block uses constant εr ≥ 1 and μr = 1, occupying electric nodes x = 46…52 and y = 33…46, with cell edges shown on the field image. Setting εr = 1 removes it. There are no material losses, dispersion or nonlinearities. The PML is entirely in vacuum and is backed by a PEC edge.

Reference cases distinguish their scope: the Courant case checks a published convention; the conductivity case is algebraic consistency; and the early scalar-wave case checks the interior stencil before PML arrival. Two additional late-time cases compare the actual absorbing run against an independent analytical unbounded-space Green function. The reference numbers are evaluated from the published exact solution, not copied from a paper’s numerical table. Physical invariants supplement these checks.

Use cell aspect ratios from 0.5 to 2, 4–24 PML cells, grading orders 1–5 and Courant safety factors 0.4–0.99. With a dielectric block, leave at least four vacuum cells before the PML. Simulation length must be a multiple of eight. Playback shows the preceding stored frame (eight-step spacing); colour limits stay fixed at ±0.04 of the source increment and larger values are colour-clipped.

The Gaussian spectrum reaches −20 dB in amplitude at f₂₀ = √ln(10)/(πτ). The reported wavelength resolution uses the highest permittivity and the larger cell size. Ten cells per wavelength is only a screening rule: all results still have directional numerical dispersion, staircased material interfaces and finite-grid PML error. Refine and test convergence in a full solver for quantitative design.

The continuum target is an amplitude ratio, with reflected power R². It predicts attenuation of an ideal normally incident return from the outer wall; it does not include discrete interface mismatch, oblique or grazing waves, the point source’s low-frequency tail, or resonances. More aggressive σmax is not guaranteed to reduce the simulated residual. Energy sampled every four steps can miss a peak and is shown with a −120 dB chart floor. The last diagnostic is at (Nrun−½)Δt. A residual is not a measured reflection coefficient.

The centred source and fixed block are mirror-symmetric in y; arbitrary off-centre geometry is outside this tool. The fixed 80×80 electric-node mesh is intentional. The article’s 200×200 grid is too heavy for this browser’s interactive computation and storage budget. Use a dedicated FDTD package for larger grids, 3D, conductive/dispersive media, CPML, calibrated scattering parameters or radiation quantities.