8 Computational Electromagnetics
Electromagnetic waves govern a remarkable range of phenomena: radio communication, radar systems, optical fibers, photonic devices, medical imaging, and even the light by which we see. This chapter shows how the finite difference techniques developed throughout this book apply to Maxwell’s equations, the fundamental laws governing electromagnetic fields.
We focus on the Finite-Difference Time-Domain (FDTD) method, introduced by Kane Yee in 1966 (Yee 1966). The FDTD method is distinguished by its use of a staggered grid (the Yee cell) where electric and magnetic field components are offset by half a grid cell in both space and time. This elegant arrangement naturally satisfies the divergence conditions and produces a stable, non-dissipative scheme.
The chapter is organized as follows:
Maxwell’s Equations (Section 8.2): We introduce the governing equations and derive the wave equation for electromagnetic fields, connecting to the scalar wave equation from Chapter 2.
The Yee Scheme (Section 8.3): We present the staggered grid discretization and leapfrog time stepping that form the core of FDTD.
1D Implementation (Section 8.4): We implement a 1D FDTD solver in Devito, demonstrating plane wave propagation and material interfaces.
2D Implementation (Section 8.5): We extend to 2D and introduce Perfectly Matched Layer (PML) absorbing boundary conditions.
Stability and Dispersion (Section 8.6): We analyze the CFL condition and numerical dispersion, showing connections to the wave equation analysis from Section 2.40.
Verification (Section 8.7): We verify our implementations using manufactured solutions and convergence testing.
Applications: We present two practical applications:
- Dielectric Waveguides (Section 8.8): Guided wave propagation in optical structures
- Ground Penetrating Radar (Section 8.9): Subsurface imaging with electromagnetic pulses
8.1 Prerequisites
This chapter assumes familiarity with:
- Finite difference methods for the wave equation (Chapter 2)
- Stability analysis and CFL conditions (Section 2.42)
- Numerical dispersion (Section 2.43)
- Devito programming patterns (Chapter 1)
8.2 Maxwell’s Equations
Maxwell’s equations describe how electric and magnetic fields are generated and altered by each other and by charges and currents. In their differential form, they are:
8.2.1 Faraday’s Law
A time-varying magnetic field induces an electric field: \[ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \tag{8.1}\]
8.2.2 Amp`ere’s Law (with Maxwell’s correction)
A current or time-varying electric field produces a magnetic field: \[ \nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t} \tag{8.2}\]
8.2.3 Gauss’s Laws
\[ \nabla \cdot \mathbf{D} = \rho, \quad \nabla \cdot \mathbf{B} = 0 \tag{8.3}\]
where \(\mathbf{E}\) is the electric field [V/m], \(\mathbf{H}\) is the magnetic field [A/m], \(\mathbf{D}\) is the electric displacement [C/m\(^2\)], \(\mathbf{B}\) is the magnetic flux density [T], \(\mathbf{J}\) is the current density [A/m\(^2\)], and \(\rho\) is the charge density [C/m\(^3\)].
8.2.4 Constitutive Relations
The fields are related through the material properties: \[ \mathbf{D} = \varepsilon \mathbf{E}, \quad \mathbf{B} = \mu \mathbf{H} \tag{8.4}\]
where \(\varepsilon\) is the permittivity [F/m] and \(\mu\) is the permeability [H/m]. In free space: \[ \varepsilon_0 = 8.854 \times 10^{-12} \text{ F/m}, \quad \mu_0 = 4\pi \times 10^{-7} \text{ H/m} \]
The speed of light follows from these constants: \[ c = \frac{1}{\sqrt{\varepsilon_0 \mu_0}} = 299{,}792{,}458 \text{ m/s} \tag{8.5}\]
In a material with relative permittivity \(\varepsilon_r\) and relative permeability \(\mu_r\), the wave speed is: \[ v = \frac{c}{\sqrt{\varepsilon_r \mu_r}} = \frac{c}{n} \tag{8.6}\]
where \(n = \sqrt{\varepsilon_r \mu_r}\) is the refractive index.
8.2.5 Derivation of the Wave Equation
To connect with the wave equation from Chapter 2, we can eliminate either \(\mathbf{E}\) or \(\mathbf{H}\) from Maxwell’s equations. Taking the curl of Faraday’s law: \[ \nabla \times (\nabla \times \mathbf{E}) = -\frac{\partial}{\partial t}(\nabla \times \mathbf{B}) = -\mu \frac{\partial}{\partial t}(\nabla \times \mathbf{H}) \]
Using the vector identity \(\nabla \times (\nabla \times \mathbf{E}) = \nabla(\nabla \cdot \mathbf{E}) - \nabla^2 \mathbf{E}\) and assuming no free charges (\(\nabla \cdot \mathbf{E} = 0\)): \[ -\nabla^2 \mathbf{E} = -\mu \frac{\partial}{\partial t}\left(\mathbf{J} + \varepsilon\frac{\partial \mathbf{E}}{\partial t}\right) \]
For source-free regions (\(\mathbf{J} = 0\)): \[ \nabla^2 \mathbf{E} = \mu \varepsilon \frac{\partial^2 \mathbf{E}}{\partial t^2} \tag{8.7}\]
This is the vector wave equation with wave speed \(c = 1/\sqrt{\mu\varepsilon}\). Similarly for \(\mathbf{H}\): \[ \nabla^2 \mathbf{H} = \mu \varepsilon \frac{\partial^2 \mathbf{H}}{\partial t^2} \tag{8.8}\]
While the wave equations 8.7 and 8.8 are mathematically equivalent to Maxwell’s equations (for linear, homogeneous media), the FDTD method works directly with the first-order system 8.1 and 8.2. This preserves both \(\mathbf{E}\) and \(\mathbf{H}\), which is essential for:
- Computing power flow (Poynting vector \(\mathbf{S} = \mathbf{E} \times \mathbf{H}\))
- Handling material interfaces correctly
- Applying boundary conditions on either field
8.2.6 Dimensional Reduction: 1D TM Mode
For our first implementation, we consider the simplest case: fields varying only in the \(x\)-direction with electric field polarized in the \(z\)-direction. This is the transverse magnetic (TM) mode with:
- \(E_z(x, t)\): electric field (z-component only)
- \(H_y(x, t)\): magnetic field (y-component only)
Maxwell’s equations reduce to: \[ \frac{\partial E_z}{\partial t} = \frac{1}{\varepsilon}\frac{\partial H_y}{\partial x} \tag{8.9}\] \[ \frac{\partial H_y}{\partial t} = \frac{1}{\mu}\frac{\partial E_z}{\partial x} \tag{8.10}\]
These are two coupled first-order PDEs. Compare with the scalar wave equation: \[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]
which is a single second-order PDE. The first-order system 8.9 and 8.10 can be converted to this form by differentiating and substituting, confirming that \(c = 1/\sqrt{\varepsilon\mu}\).
8.2.7 Lossy Media
In conductive media, Ohm’s law adds a current density term: \[ \mathbf{J} = \sigma \mathbf{E} \]
where \(\sigma\) is the electrical conductivity [S/m]. The 1D equations become: \[ \frac{\partial E_z}{\partial t} = \frac{1}{\varepsilon}\frac{\partial H_y}{\partial x} - \frac{\sigma}{\varepsilon} E_z \tag{8.11}\]
The additional term \(-(\sigma/\varepsilon) E_z\) causes exponential decay of the electric field amplitude, representing absorption of electromagnetic energy as heat in the medium.
8.2.8 Wave Impedance
The ratio of electric to magnetic field amplitudes in a plane wave defines the wave impedance: \[ \eta = \frac{|E|}{|H|} = \sqrt{\frac{\mu}{\varepsilon}} \tag{8.12}\]
In free space, \(\eta_0 = \sqrt{\mu_0/\varepsilon_0} \approx 377\ \Omega\). The impedance is crucial for:
- Computing reflection and transmission at interfaces
- Verifying simulation results (E and H must be related by \(\eta\))
- Designing matched absorbing boundaries
8.2.9 Boundary Conditions
Perfect Electric Conductor (PEC): The tangential electric field vanishes: \[ \mathbf{n} \times \mathbf{E} = 0 \quad \text{(PEC)} \tag{8.13}\]
For our 1D case with \(E_z\), this means \(E_z = 0\) at the boundary.
Perfect Magnetic Conductor (PMC): The tangential magnetic field vanishes: \[ \mathbf{n} \times \mathbf{H} = 0 \quad \text{(PMC)} \tag{8.14}\]
Absorbing Boundary Conditions: For open-domain problems, we need boundaries that absorb outgoing waves without reflection. The simplest approach is Mur’s first-order ABC (Mur 1981), based on the one-way wave equation: \[ \frac{\partial E_z}{\partial t} + c \frac{\partial E_z}{\partial x} = 0 \]
for a wave traveling in the \(+x\) direction. More sophisticated approaches include the Perfectly Matched Layer (PML) (Berenger 1994), covered in Section 8.5.
8.3 The Yee Scheme
The Finite-Difference Time-Domain (FDTD) method was introduced by Kane Yee in 1966 (Yee 1966). Its key innovation is the staggered grid, where electric and magnetic field components are offset by half a grid cell in both space and time. This arrangement, known as the Yee cell, naturally represents the curl operations in Maxwell’s equations.
8.3.1 The Staggered Grid Concept
Unlike the collocated grid used for the scalar wave equation in Chapter 2, where all quantities are defined at the same grid points, the Yee scheme places \(E\) and \(H\) at different locations:
1D Staggered Grid:
- \(E_z\) is defined at integer grid points: \(x_i = i \cdot \Delta x\)
- \(H_y\) is defined at half-integer points: \(x_{i+1/2} = (i + 1/2) \cdot \Delta x\)
E_z[0] E_z[1] E_z[2] E_z[3] E_z[4]
| | | | |
*---------*---------*---------*---------*
H_y[0] H_y[1] H_y[2] H_y[3]Similarly for time:
- \(E_z\) is computed at integer time steps: \(t^n = n \cdot \Delta t\)
- \(H_y\) is computed at half-integer times: \(t^{n+1/2} = (n + 1/2) \cdot \Delta t\)
This staggering means that every derivative in Maxwell’s equations becomes a centered difference, achieving second-order accuracy automatically.
8.3.2 Leapfrog Time Integration
The time stepping alternates between updating \(E\) and \(H\) fields:
- Use \(H^{n-1/2}\) to compute \(E^n\)
- Use \(E^n\) to compute \(H^{n+1/2}\)
- Use \(H^{n+1/2}\) to compute \(E^{n+1}\)
- Repeat…
This is a leapfrog scheme where each field “leaps over” the other in time. The method is:
- Explicit: No linear system to solve
- Non-dissipative: Energy is conserved (in lossless media)
- Second-order accurate: In both space and time
Compare this with the leapfrog scheme for the wave equation in Section 2.6, which also achieves second-order accuracy through centered differences.
8.3.3 1D FDTD Discretization
Starting from the 1D Maxwell equations 8.9 and 8.10:
\(H_y\) update (from \(t^{n-1/2}\) to \(t^{n+1/2}\)):
Using Faraday’s law \(\partial H_y/\partial t = (1/\mu) \partial E_z/\partial x\): \[ \frac{H_y|_{i+1/2}^{n+1/2} - H_y|_{i+1/2}^{n-1/2}}{\Delta t} = \frac{1}{\mu} \frac{E_z|_{i+1}^{n} - E_z|_{i}^{n}}{\Delta x} \tag{8.15}\]
Solving for the new value: \[ H_y|_{i+1/2}^{n+1/2} = H_y|_{i+1/2}^{n-1/2} + \frac{\Delta t}{\mu \Delta x}\left(E_z|_{i+1}^{n} - E_z|_{i}^{n}\right) \tag{8.16}\]
\(E_z\) update (from \(t^{n}\) to \(t^{n+1}\)):
Using Amp`ere’s law \(\partial E_z/\partial t = (1/\varepsilon) \partial H_y/\partial x\): \[ \frac{E_z|_{i}^{n+1} - E_z|_{i}^{n}}{\Delta t} = \frac{1}{\varepsilon} \frac{H_y|_{i+1/2}^{n+1/2} - H_y|_{i-1/2}^{n+1/2}}{\Delta x} \tag{8.17}\]
Solving for the new value: \[ E_z|_{i}^{n+1} = E_z|_{i}^{n} + \frac{\Delta t}{\varepsilon \Delta x}\left(H_y|_{i+1/2}^{n+1/2} - H_y|_{i-1/2}^{n+1/2}\right) \tag{8.18}\]
8.3.4 Update Coefficients
For uniform materials, we can precompute the update coefficients: \[ C_h = \frac{\Delta t}{\mu \Delta x}, \quad C_e = \frac{\Delta t}{\varepsilon \Delta x} \tag{8.19}\]
The update equations become: \[ H_y|_{i+1/2}^{n+1/2} = H_y|_{i+1/2}^{n-1/2} + C_h \left(E_z|_{i+1}^{n} - E_z|_{i}^{n}\right) \] \[ E_z|_{i}^{n+1} = E_z|_{i}^{n} + C_e \left(H_y|_{i+1/2}^{n+1/2} - H_y|_{i-1/2}^{n+1/2}\right) \]
For spatially varying materials, \(C_h\) and \(C_e\) become arrays indexed by position.
8.3.5 Lossy Media Update
For media with conductivity \(\sigma\), the \(E_z\) update becomes: \[ E_z|_{i}^{n+1} = C_a E_z|_{i}^{n} + C_b \left(H_y|_{i+1/2}^{n+1/2} - H_y|_{i-1/2}^{n+1/2}\right) \tag{8.20}\]
where: \[ C_a = \frac{1 - \sigma\Delta t/(2\varepsilon)}{1 + \sigma\Delta t/(2\varepsilon)}, \quad C_b = \frac{\Delta t/(\varepsilon\Delta x)}{1 + \sigma\Delta t/(2\varepsilon)} \tag{8.21}\]
This semi-implicit treatment of the loss term maintains second-order accuracy.
8.3.6 The 2D Yee Cell
In two dimensions, a common pedagogical choice is the TM polarization (with \(E_z\), \(H_x\), \(H_y\)), where the electric field is out of the simulation plane. This reduces the field update to a single electric component and two magnetic components, which is convenient for introductory FDTD and matches the implementation in src.em.maxwell2D_devito.
H_y[i,j] H_y[i,j+1]
| |
| E_z |
--------|------ * ------|--------
| [i,j] |
H_x[i,j]| |H_x[i+1,j]
| |
--------|---------------|--------
| |
H_y[i+1,j] H_y[i+1,j+1]The field components are located at:
- \(E_z\): cell centers \((i, j)\)
- \(H_x\): cell edges parallel to \(x\)-axis \((i+1/2, j)\)
- \(H_y\): cell edges parallel to \(y\)-axis \((i, j+1/2)\)
This arrangement ensures that the discrete curl terms are built from centered differences between staggered components.
8.3.7 2D Update Equations (TM Mode)
\(H_x\) update: \[ H_x|_{i,j}^{n+1/2} = H_x|_{i,j}^{n-1/2} - \frac{\Delta t}{\mu}\frac{E_z|_{i,j+1}^{n} - E_z|_{i,j}^{n}}{\Delta y} \tag{8.22}\]
\(H_y\) update: \[ H_y|_{i,j}^{n+1/2} = H_y|_{i,j}^{n-1/2} + \frac{\Delta t}{\mu}\frac{E_z|_{i+1,j}^{n} - E_z|_{i,j}^{n}}{\Delta x} \tag{8.23}\]
\(E_z\) update: \[ E_z|_{i,j}^{n+1} = E_z|_{i,j}^{n} + \frac{\Delta t}{\varepsilon}\left( \frac{H_y|_{i,j}^{n+1/2} - H_y|_{i-1,j}^{n+1/2}}{\Delta x} - \frac{H_x|_{i,j}^{n+1/2} - H_x|_{i,j-1}^{n+1/2}}{\Delta y} \right) \tag{8.24}\]
8.3.8 CFL Stability Condition
The Yee scheme is conditionally stable. Von Neumann analysis (similar to Section 2.42) yields the CFL condition:
1D: \[ c \Delta t \leq \Delta x \tag{8.25}\]
2D: \[ c \Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}} \tag{8.26}\]
For a uniform grid with \(\Delta x = \Delta y\): \[ c \Delta t \leq \frac{\Delta x}{\sqrt{2}} \tag{8.27}\]
3D: \[ c \Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}}} \tag{8.28}\]
Yee’s original 1966 paper (Yee 1966) contained an error in the stability condition for 3D. The correct condition was established by Taflove and Brodwin in 1975 (Taflove and Brodwin 1975).
8.3.9 Boundary Conditions
Perfect Electric Conductor (PEC):
For \(E_z = 0\) at a PEC boundary, we simply set \(E_z = 0\) at the boundary grid points. No special stencil is needed.
Perfect Magnetic Conductor (PMC):
For \(H_y = 0\) at a PMC boundary (normal \(H\) component vanishes), we set the \(H_y\) values at the boundary to zero.
Absorbing Boundaries:
For open-domain problems, simple boundary conditions cause spurious reflections. Options include:
- First-order ABC (Mur) (Mur 1981): \(E^{n+1}_0 = E^n_1 + \frac{c\Delta t - \Delta x}{c\Delta t + \Delta x}(E^{n+1}_1 - E^n_0)\)
- Perfectly Matched Layer (PML) (Berenger 1994): A specially designed absorbing region (see Section 8.5)
8.3.10 Connection to Scalar Wave Equation
The Yee scheme for Maxwell’s equations can be viewed as a first-order splitting of the second-order wave equation. If we define: \[ u = E_z, \quad v = \eta H_y \] where \(\eta = \sqrt{\mu/\varepsilon}\) is the wave impedance, the 1D Maxwell equations become: \[ \frac{\partial u}{\partial t} = c \frac{\partial v}{\partial x}, \quad \frac{\partial v}{\partial t} = c \frac{\partial u}{\partial x} \]
This is a symmetric hyperbolic system. Eliminating \(v\): \[ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} \]
The centered differences in the Yee scheme correspond to the standard second-order scheme for the wave equation, explaining why both achieve the same accuracy and stability properties.
8.4 1D FDTD Implementation
We now implement the 1D Yee scheme in Python. Unlike the scalar wave equation where we could use Devito’s automatic differentiation directly, the staggered grid requires more careful handling of the half-integer offsets.
8.4.1 The Staggered Grid Challenge
Devito’s TimeFunction places all values on a collocated grid. For the Yee scheme, we need:
- \(E_z\) at positions \(x_i = i \cdot \Delta x\) for \(i = 0, 1, \ldots, N_x\)
- \(H_y\) at positions \(x_{i+1/2} = (i + 0.5) \cdot \Delta x\) for \(i = 0, 1, \ldots, N_x - 1\)
We handle this by interpreting the \(H_y\) array with shifted indices: H_y[i] represents \(H_y\) at location \(x_{i+1/2}\).
8.4.2 The Solver Module
The src.em.maxwell1D_devito module provides a complete 1D FDTD solver:
from src.em import solve_maxwell_1d, MaxwellResult1D
import numpy as np
# Gaussian pulse initial condition
def gaussian_pulse(x, x0=0.5, sigma=0.05):
return np.exp(-((x - x0)**2) / (2 * sigma**2))
# Solve 1D Maxwell equations
result = solve_maxwell_1d(
L=1.0, # Domain length [m]
Nx=200, # Number of grid cells
T=5e-9, # Simulation time [s]
CFL=0.9, # Courant number
E_init=gaussian_pulse,
save_history=True,
)
print(f"Wave speed: {result.c:.2e} m/s")
print(f"Time step: {result.dt:.2e} s")
print(f"Grid spacing: {result.dx:.2e} m")8.4.3 Understanding the Implementation
The core of the solver implements the update equations 8.16 and 8.18. Here is a simplified version:
# Update coefficients
Ch = dt / (mu * dx) # H update coefficient
Ce = dt / (eps * dx) # E update coefficient
# Time stepping loop
for n in range(Nt):
# Update H_y (uses current E_z)
# H_y[i] represents H_y at x_{i+1/2}
H_y[:] = H_y[:] + Ch * (E_z[1:] - E_z[:-1])
# Update E_z (uses new H_y)
E_z[1:-1] = E_z[1:-1] + Ce * (H_y[1:] - H_y[:-1])
# Apply boundary conditions
E_z[0] = 0.0 # PEC at left
E_z[-1] = 0.0 # PEC at rightNote how the array slicing naturally handles the staggered indexing:
E_z[1:] - E_z[:-1]computes \(E_z|_{i+1} - E_z|_{i}\) for all \(i\)H_y[1:] - H_y[:-1]computes \(H_y|_{i+1/2} - H_y|_{i-1/2}\) for interior points
8.4.4 Visualizing Wave Propagation
Let’s simulate a Gaussian pulse propagating and reflecting from PEC boundaries:
import matplotlib.pyplot as plt
from src.em import solve_maxwell_1d
# Initial Gaussian pulse
def gaussian(x):
return np.exp(-((x - 0.3)**2) / (0.02**2))
result = solve_maxwell_1d(
L=1.0, Nx=400, T=8e-9, CFL=0.9,
E_init=gaussian,
save_history=True,
)
# Plot snapshots
fig, axes = plt.subplots(2, 3, figsize=(12, 6))
times_idx = [0, 20, 40, 60, 80, 100]
for ax, idx in zip(axes.flat, times_idx):
if idx < len(result.E_history):
ax.plot(result.x_E, result.E_history[idx])
ax.set_title(f't = {result.t_history[idx]*1e9:.2f} ns')
ax.set_xlabel('x [m]')
ax.set_ylabel('E_z [V/m]')
ax.set_ylim(-1.2, 1.2)
plt.tight_layout()8.4.5 Material Interfaces
One of FDTD’s strengths is handling material interfaces. The update coefficients simply use the local material properties:
# Two-layer medium: vacuum (eps_r=1) and glass (eps_r=4)
eps_r = np.ones(Nx + 1)
eps_r[Nx//2:] = 4.0 # Glass in second half
result = solve_maxwell_1d(
L=1.0, Nx=400, T=8e-9,
eps_r=eps_r,
E_init=gaussian,
save_history=True,
)At the interface, the wave partially reflects and partially transmits. The reflection coefficient for normal incidence is: \[ R = \frac{\eta_2 - \eta_1}{\eta_2 + \eta_1} = \frac{\sqrt{\varepsilon_1} - \sqrt{\varepsilon_2}}{\sqrt{\varepsilon_1} + \sqrt{\varepsilon_2}} \]
For our vacuum/glass interface: \[ R = \frac{1 - 2}{1 + 2} = -\frac{1}{3} \]
The negative sign indicates a phase reversal upon reflection.
8.4.6 Source Injection
For many applications, we inject a source rather than using initial conditions. The soft source adds to the existing field:
from src.em import ricker_wavelet
# Ricker wavelet source at 500 MHz
def source(t):
return ricker_wavelet(np.array([t]), f0=500e6)[0]
result = solve_maxwell_1d(
L=1.0, Nx=400, T=10e-9,
source_func=source,
source_position=0.1, # 10 cm from left boundary
bc_left="abc", # Absorbing BC to prevent reflection
bc_right="abc",
save_history=True,
)The Ricker wavelet (Mexican hat) is commonly used in GPR and seismic modeling: \[ r(t) = \left(1 - 2\pi^2 f_0^2 (t - t_0)^2\right) e^{-\pi^2 f_0^2 (t - t_0)^2} \tag{8.29}\]
8.4.7 Absorbing Boundary Conditions
Simple PEC boundaries cause total reflection. For open-domain problems, we need absorbing boundaries. The simplest approach uses a first-order ABC:
# At right boundary: wave traveling in +x direction
# One-way wave equation: dE/dt + c*dE/dx = 0
# Discretized: E[N]^{n+1} = E[N-1]^n + (C-1)/(C+1) * (E[N-1]^{n+1} - E[N]^n)This works well for normal incidence but becomes less effective at oblique angles. The 2D implementation in Section 8.5 introduces the more robust Perfectly Matched Layer.
8.4.8 Field Relationship Verification
A key verification is that \(E_z\) and \(H_y\) maintain the correct amplitude ratio. For a plane wave: \[ \frac{|E_z|}{|H_y|} = \eta = \sqrt{\frac{\mu}{\varepsilon}} \]
In free space, \(\eta_0 \approx 377\ \Omega\). We can verify this:
from src.em import verify_units, EMConstants
const = EMConstants()
E_max = np.max(np.abs(result.E_z))
H_max = np.max(np.abs(result.H_y))
ratio = E_max / H_max
expected = const.eta0
print(f"E/H ratio: {ratio:.1f} Ohm")
print(f"Expected: {expected:.1f} Ohm")
print(f"Error: {100*abs(ratio-expected)/expected:.2f}%")8.4.9 Complete Example: Plane Wave Verification
Here is a complete example that verifies our implementation against the exact plane wave solution:
from src.em import solve_maxwell_1d, exact_plane_wave_1d
import numpy as np
# Physical parameters
L = 1.0 # Domain length [m]
Nx = 200
wavelength = 0.1 # [m]
k = 2 * np.pi / wavelength
# Initial plane wave
def E_init(x):
return np.cos(k * x)
def H_init(x):
from src.em import EMConstants
return np.cos(k * x) / EMConstants().eta0
# Run short simulation (avoid boundary effects)
T = 1e-9 # 1 ns
result = solve_maxwell_1d(
L=L, Nx=Nx, T=T, CFL=0.9,
E_init=E_init, H_init=H_init,
bc_left="abc", bc_right="abc",
)
# Compare with exact solution
E_exact, H_exact = exact_plane_wave_1d(
result.x_E, result.t,
amplitude=1.0, k=k
)
error = np.sqrt(np.mean((result.E_z - E_exact)**2))
print(f"L2 error: {error:.2e}")The error should be small for short simulation times, limited mainly by the ABC’s inability to perfectly absorb non-normal components of the numerical dispersion.
8.5 2D FDTD Implementation
Extending to two dimensions introduces new challenges: the 2D Yee cell has three field components arranged on cell edges and centers, and the stability condition becomes more restrictive. Most importantly, we need absorbing boundary conditions that work at oblique incidence angles.
8.5.1 TE and TM Modes in 2D
In 2D simulations, we have two independent polarizations:
TE Mode (Transverse Electric): \(E_x\), \(E_y\), \(H_z\)
- Electric field lies in the \(xy\)-plane
- Magnetic field is normal to the plane
TM Mode (Transverse Magnetic): \(H_x\), \(H_y\), \(E_z\)
- Magnetic field lies in the \(xy\)-plane
- Electric field is normal to the plane
We implement the TM mode (a single out-of-plane electric component). The governing equations are: \[ \frac{\partial H_x}{\partial t} = -\frac{1}{\mu}\frac{\partial E_z}{\partial y} \tag{8.30}\] \[ \frac{\partial H_y}{\partial t} = \frac{1}{\mu}\frac{\partial E_z}{\partial x} \tag{8.31}\] \[ \frac{\partial E_z}{\partial t} = \frac{1}{\varepsilon}\left(\frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}\right) \tag{8.32}\]
8.5.2 The 2D Solver
The src.em.maxwell2D_devito module provides a 2D TM-mode solver using Devito (Luporini et al. 2020):
from src.em import solve_maxwell_2d, gaussian_source_2d
import numpy as np
# Gaussian initial condition at domain center
def E_init(X, Y):
return gaussian_source_2d(X, Y, x0=0.5, y0=0.5, sigma=0.05)
result = solve_maxwell_2d(
Lx=1.0, Ly=1.0, # Domain size [m]
Nx=100, Ny=100, # Grid points
T=3e-9, # Simulation time [s]
CFL=0.5, # Courant number (<= 1/sqrt(2) for stability)
E_init=E_init,
pml_width=15, # PML absorbing boundary (grid cells)
save_history=True,
)8.5.3 2D CFL Condition
The stability limit in 2D is more restrictive than 1D: \[ c \Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}} \tag{8.33}\]
For a uniform grid (\(\Delta x = \Delta y = h\)): \[ c \Delta t \leq \frac{h}{\sqrt{2}} \approx 0.707 h \]
This means the maximum Courant number is \(1/\sqrt{2} \approx 0.707\), compared to 1.0 in 1D. The solver enforces this:
>>> from src.em import solve_maxwell_2d
>>> solve_maxwell_2d(Lx=1.0, Ly=1.0, Nx=100, Ny=100, T=1e-9, CFL=0.9)
ValueError: CFL=0.9 > 1/sqrt(2) ~ 0.707 violates 2D stability condition8.5.4 Perfectly Matched Layer (PML)
Simple ABCs like Mur’s condition work poorly in 2D because waves approach boundaries at various angles. The Perfectly Matched Layer (PML), introduced by Berenger in 1994 (Berenger 1994), provides a much more effective solution.
The key insight is to create an absorbing region with:
- Matched impedance: No reflection at the interface (for any angle)
- Exponential decay: Waves attenuate as they propagate into the PML
Mathematically, this is achieved through complex coordinate stretching: \[ \tilde{x} = x + \frac{j}{\omega}\int_0^x \sigma_x(x') dx' \]
This makes the wave “see” a lossy medium while maintaining impedance matching. For the acoustic wave equation, an efficient PML formulation is given by Grote and Sim (2010). The Convolutional PML (CPML) variant (Roden and Gedney 2000) offers improved stability for anisotropic and dispersive media.
8.5.5 PML Implementation
Our solver supports two PML implementations, selectable via the pml_type parameter:
Graded-conductivity absorbing layer (pml_type='conductivity'). The simplest approach adds a spatially varying conductivity to the update equations. The conductivity increases polynomially from zero at the PML interface to \(\sigma_{\max}\) at the outer boundary:
def create_pml_profile(N, pml_width, sigma_max, order=3):
"""Create polynomial PML conductivity profile."""
sigma = np.zeros(N)
for i in range(pml_width):
d = (pml_width - i) / pml_width
sigma[i] = sigma_max * (d ** order) # Left PML
for i in range(N - pml_width, N):
d = (i - (N - pml_width - 1)) / pml_width
sigma[i] = sigma_max * (d ** order) # Right PML
return sigmaConvolutional PML (pml_type='cpml', default). The CPML (Roden and Gedney 2000) uses the complex frequency-shifted (CFS) stretching function and implements the PML via recursive convolution. The key update for the auxiliary \(\Psi\) fields is: \[
\Psi_x^{n+1} = b_x \Psi_x^n + a_x \frac{\partial E_z}{\partial x}\bigg|^n,
\tag{8.34}\] where the coefficients \(b_x\) and \(a_x\) are derived from the CFS-PML parameters: \[
b_x = e^{-(\sigma_x/\kappa_x + \alpha_x)\Delta t}, \quad
a_x = \frac{\sigma_x}{\kappa_x(\sigma_x + \kappa_x \alpha_x)}(b_x - 1).
\]
The CPML modifies the spatial derivatives in both the H and E updates: \[ \left.\frac{\partial E_z}{\partial x}\right|_{\text{PML}} = \frac{1}{\kappa_x}\frac{\partial E_z}{\partial x} + \Psi_x. \]
This approach is superior to the simple conductivity method because:
- It provides true impedance matching at the PML interface
- It handles evanescent waves (with \(\alpha > 0\))
- It is stable for long-time simulations and dispersive media
Typical parameters:
pml_width: 10–20 grid cellsorder: 3–4 (polynomial order for \(\sigma\) grading)sigma_max: Computed from optimal formula
The optimal \(\sigma_{\max}\) minimizes total reflection (numerical + PML): \[ \sigma_{\max} = \frac{(m+1)}{150\pi \Delta x} \tag{8.35}\]
where \(m\) is the polynomial order.
8.5.6 Visualizing 2D Propagation
import matplotlib.pyplot as plt
from src.em import solve_maxwell_2d, gaussian_source_2d
import numpy as np
# Point source excitation
def source(t):
from src.em import ricker_wavelet
return ricker_wavelet(np.array([t]), f0=1e9)[0]
result = solve_maxwell_2d(
Lx=0.5, Ly=0.5, Nx=200, Ny=200, T=2e-9, CFL=0.5,
source_func=source,
source_position=(0.25, 0.25),
pml_width=20,
save_history=True,
save_every=10,
)
# Plot snapshots
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
for ax, i in zip(axes.flat, range(0, len(result.E_history), len(result.E_history)//6)):
im = ax.imshow(result.E_history[i].T, origin='lower',
extent=[0, 0.5, 0, 0.5], cmap='RdBu',
vmin=-1, vmax=1)
ax.set_title(f't = {result.t_history[i]*1e9:.2f} ns')
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
plt.tight_layout()The simulation shows:
- Circular wavefront: Expanding from the point source
- No visible reflections: PML absorbs outgoing waves
- Correct propagation speed: Wavefront radius = \(c \times t\)
8.5.7 Scattering from Objects
A classic FDTD application is computing scattering from objects. We can embed a dielectric or conducting scatterer:
from src.em.materials import create_cylinder_model_2d, GLASS
# Create cylinder scatterer
eps_r, sigma = create_cylinder_model_2d(
Nx=200, Ny=200, Lx=0.5, Ly=0.5,
center=(0.25, 0.25),
radius=0.03,
cylinder_material=GLASS,
)
result = solve_maxwell_2d(
Lx=0.5, Ly=0.5, Nx=200, Ny=200, T=2e-9,
eps_r=eps_r,
source_func=source,
source_position=(0.1, 0.25), # Source to left of cylinder
pml_width=20,
save_history=True,
)The scattered field shows:
- Reflection from the front surface
- Transmission through the cylinder
- Diffraction around the edges
- Internal resonances for certain sizes
8.5.8 Connection to Wave Chapter ABCs
The PML can be viewed as a sophisticated extension of the absorbing boundary conditions discussed for the scalar wave equation (Section 2.50). Compare:
| Method | Principle | Angle Dependence | Complexity |
|---|---|---|---|
| First-order ABC | One-way wave equation | Normal incidence only | Simple |
| Higher-order ABC | Multiple angles | Improved | Moderate |
| PML | Impedance matching | All angles | Higher |
The PML achieves angle-independent absorption through the impedance matching condition, which ensures zero reflection at the PML interface regardless of incidence angle.
8.5.9 Practical Considerations
Grid Resolution:
The rule of thumb is 10-20 grid points per wavelength for acceptable accuracy. For a 1 GHz signal in vacuum: \[ \lambda = \frac{c}{f} = \frac{3 \times 10^8}{10^9} = 0.3 \text{ m} \]
So we need \(\Delta x \leq 0.03\) m (30 mm) for 10 points per wavelength.
Computational Cost:
2D FDTD scales as \(O(N_x \times N_y \times N_t)\). For a \(200 \times 200\) grid with 1000 time steps, this is 40 million field updates. Each requires only a few floating-point operations, making FDTD very efficient.
Memory:
We need to store 3 field arrays (\(E_x\), \(E_y\), \(H_z\)) plus material property arrays. For a \(200 \times 200\) grid in double precision: \[ \text{Memory} \approx 6 \times 200 \times 200 \times 8 \text{ bytes} \approx 2 \text{ MB} \]
This is modest by modern standards, allowing much larger simulations.
8.6 Stability and Dispersion Analysis
The Yee scheme shares fundamental properties with the centered-difference scheme for the scalar wave equation. This section analyzes stability and numerical dispersion, drawing parallels to the wave equation analysis in Section 2.40.
8.6.1 Von Neumann Stability Analysis
Following the approach in Section 2.42, we insert a plane wave ansatz into the discrete equations: \[ E_z|_i^n = \hat{E} e^{j(kx_i - \tilde\omega t_n)}, \quad H_y|_{i+1/2}^{n+1/2} = \hat{H} e^{j(kx_{i+1/2} - \tilde\omega t_{n+1/2})} \]
where \(\tilde\omega\) is the numerical angular frequency (which may differ from the physical \(\omega = ck\)).
Substituting into the update equations 8.16 and 8.18 and simplifying yields the numerical dispersion relation: \[ \left(\frac{\sin(\tilde\omega \Delta t/2)}{\Delta t/2}\right)^2 = c^2 \left(\frac{\sin(k\Delta x/2)}{\Delta x/2}\right)^2 \tag{8.36}\]
This can be rewritten as: \[ \sin(\tilde\omega \Delta t/2) = C \sin(k\Delta x/2) \tag{8.37}\]
where \(C = c\Delta t/\Delta x\) is the Courant number.
Stability Condition: For \(\tilde\omega\) to be real (non-growing solutions), we need: \[ |C \sin(k\Delta x/2)| \leq 1 \quad \text{for all } k \]
The maximum occurs at \(k\Delta x = \pi\), requiring \(|C| \leq 1\), which gives the CFL condition 8.25.
8.6.2 Numerical Phase Velocity
The numerical phase velocity is: \[ \tilde{v}_p = \frac{\tilde\omega}{k} \]
From 8.37: \[ \tilde\omega = \frac{2}{\Delta t}\arcsin\left(C \sin\frac{k\Delta x}{2}\right) \]
The phase velocity error is: \[ \frac{\tilde{v}_p - c}{c} = \frac{\tilde\omega}{kc} - 1 \tag{8.38}\]
This error depends on:
- Courant number \(C\)
- Points per wavelength \(N_\lambda = \lambda/\Delta x = 2\pi/(k\Delta x)\)
8.6.3 The “Magic” Time Step
A remarkable property of the 1D Yee scheme: at \(C = 1\), the dispersion vanishes exactly. From 8.37 with \(C = 1\): \[ \sin(\tilde\omega \Delta t/2) = \sin(k\Delta x/2) \]
With \(c\Delta t = \Delta x\) (i.e., \(C = 1\)): \[ \tilde\omega = \frac{2}{\Delta t}\arcsin\left(\sin\frac{k\Delta x}{2}\right) = \frac{k\Delta x}{\Delta t} = kc \]
This is the exact dispersion relation! Waves of all frequencies travel at exactly the correct speed.
The magic time step \(C = 1\) makes the numerical scheme exact (in 1D for uniform media). This is the same phenomenon observed for the scalar wave equation in Section 2.43. It arises because the centered difference stencil, when \(C = 1\), samples the wave at exactly the points where information propagates.
8.6.4 Dispersion in 2D/3D
In higher dimensions, the dispersion becomes anisotropic—it depends on the propagation direction. For 2D with \(\Delta x = \Delta y = h\): \[ \sin^2\left(\frac{\tilde\omega \Delta t}{2}\right) = S_x^2 \sin^2\left(\frac{k_x h}{2}\right) + S_y^2 \sin^2\left(\frac{k_y h}{2}\right) \tag{8.39}\]
where \(S_x = c\Delta t/\Delta x\) and \(S_y = c\Delta t/\Delta y\).
The phase velocity error varies with angle \(\theta\) (where \(k_x = k\cos\theta\), \(k_y = k\sin\theta\)). For a given resolution, dispersion is:
- Minimum along grid axes (\(\theta = 0, 90^\circ\))
- Maximum along diagonals (\(\theta = 45^\circ\))
This grid anisotropy is inherent to the Cartesian Yee grid. It can be reduced by using finer grids or higher-order schemes.
8.6.5 Quantifying Dispersion Error
The src.em.analysis.dispersion_maxwell module provides tools for computing dispersion errors:
from src.em.analysis import compute_dispersion_error, phase_velocity_error_1d
import numpy as np
# Error vs. points per wavelength for various Courant numbers
N_lambda = np.array([5, 10, 20, 40, 80])
for C in [0.5, 0.9, 1.0]:
error = compute_dispersion_error(N_lambda, courant_number=C, dim=1)
print(f"C = {C}:")
for n, e in zip(N_lambda, error):
print(f" N_lambda = {n:2d}: error = {100*e:+.3f}%")Typical output:
C = 0.5:
N_lambda = 5: error = -2.467%
N_lambda = 10: error = -0.617%
N_lambda = 20: error = -0.154%
N_lambda = 40: error = -0.039%
N_lambda = 80: error = -0.010%
C = 0.9:
N_lambda = 5: error = -0.449%
N_lambda = 10: error = -0.112%
N_lambda = 20: error = -0.028%
N_lambda = 40: error = -0.007%
N_lambda = 80: error = -0.002%
C = 1.0:
N_lambda = 5: error = +0.000%
N_lambda = 10: error = +0.000%
N_lambda = 20: error = +0.000%
N_lambda = 40: error = +0.000%
N_lambda = 80: error = +0.000%8.6.6 Practical Implications
Resolution Requirements:
The rule of thumb is 10-20 points per minimum wavelength. With \(N_\lambda = 10\):
- At \(C = 0.5\): ~0.6% phase error per wavelength traveled
- At \(C = 0.9\): ~0.1% phase error per wavelength traveled
For a wave traveling 100 wavelengths, these accumulate to 60% and 10% total phase error, respectively. Long-distance propagation requires either higher resolution or \(C\) closer to 1.
Broadband Signals:
Signals with multiple frequency components (like the Ricker wavelet in GPR) disperse because different frequencies travel at different speeds. The pulse broadens and develops oscillatory tails. Higher resolution mitigates this.
Grid Design:
When possible, design the grid so that the primary propagation direction aligns with grid axes (lower dispersion). For problems with waves traveling in all directions (e.g., scattering), use sufficient resolution to keep anisotropy acceptable.
8.6.7 Comparison with Wave Equation
The dispersion relation 8.36 is identical to that for the scalar wave equation discretized with central differences (compare with Section 2.43). This is not coincidental—the Yee scheme can be viewed as a first-order system reformulation of the second-order wave equation, and both share the same dispersion properties.
| Property | Scalar Wave | Maxwell (Yee) |
|---|---|---|
| Order of accuracy | 2nd in space and time | 2nd in space and time |
| CFL limit (1D) | \(C \leq 1\) | \(C \leq 1\) |
| CFL limit (2D) | \(C \leq 1/\sqrt{2}\) | \(C \leq 1/\sqrt{2}\) |
| Magic time step | \(C = 1\) (exact) | \(C = 1\) (exact) |
| Grid anisotropy | Yes (2D/3D) | Yes (2D/3D) |
The key difference is that Maxwell’s equations preserve both \(\mathbf{E}\) and \(\mathbf{H}\) fields, which is essential for computing quantities like the Poynting vector and handling material interfaces correctly.
8.6.8 Group Velocity
The group velocity determines how wave packets (pulses) propagate: \[ v_g = \frac{\partial\omega}{\partial k} \]
From the numerical dispersion relation: \[ \tilde{v}_g = c \frac{\cos(k\Delta x/2)}{\sqrt{1 - C^2\sin^2(k\Delta x/2)}} \tag{8.40}\]
For \(C = 1\), \(\tilde{v}_g = c\cos(k\Delta x/2)\), which is less than \(c\) for finite wavelengths. This means even at the magic time step, pulses experience some dispersion due to their finite bandwidth.
8.7 Verification
Verification ensures our implementation correctly solves the intended equations. We apply the techniques from Section 1.5: exact solutions, convergence testing, and the Method of Manufactured Solutions (MMS).
8.7.1 Verification Checklist
A complete verification suite for the Maxwell solver should confirm:
- Wave speed: A pulse travels at \(c = 1/\sqrt{\varepsilon\mu}\)
- PEC reflection: Electric field vanishes at PEC boundaries
- Impedance relation: \(|E|/|H| = \eta = \sqrt{\mu/\varepsilon}\)
- Convergence rate: Second-order in space and time
- Energy conservation: Total EM energy is constant (lossless case)
8.7.2 Exact Solution: Plane Wave
The simplest exact solution is a plane wave: \[ E_z(x, t) = E_0 \cos(kx - \omega t), \quad H_y(x, t) = \frac{E_0}{\eta} \cos(kx - \omega t) \tag{8.41}\]
where \(\omega = kc\) and \(\eta = \sqrt{\mu/\varepsilon}\).
from src.em import solve_maxwell_1d, exact_plane_wave_1d, EMConstants
import numpy as np
const = EMConstants()
L = 1.0
Nx = 200
wavelength = 0.1
k = 2 * np.pi / wavelength
# Initial condition from exact solution
def E_init(x):
return np.cos(k * x)
def H_init(x):
return np.cos(k * x) / const.eta0
# Run simulation
T = 0.5e-9 # Short time to avoid boundary effects
result = solve_maxwell_1d(
L=L, Nx=Nx, T=T, CFL=0.9,
E_init=E_init, H_init=H_init,
bc_left="abc", bc_right="abc",
)
# Compare with exact
E_exact, _ = exact_plane_wave_1d(result.x_E, result.t, k=k)
error = np.sqrt(np.mean((result.E_z - E_exact)**2))
print(f"L2 error: {error:.2e}")8.7.3 Wave Speed Verification
To verify the wave travels at the correct speed, we track a pulse peak:
from src.em.verification import verify_wave_speed, EMConstants
import numpy as np
const = EMConstants()
# Gaussian pulse simulation
result = solve_maxwell_1d(
L=1.0, Nx=400, T=3e-9, CFL=0.9,
E_init=lambda x: np.exp(-((x - 0.2)**2) / 0.01**2),
bc_left="abc", bc_right="abc",
save_history=True,
)
passed, measured_c = verify_wave_speed(
result.E_history,
result.x_E,
result.t_history,
expected_c=const.c0,
tolerance=0.05,
)
print(f"Expected speed: {const.c0:.3e} m/s")
print(f"Measured speed: {measured_c:.3e} m/s")
print(f"Verification: {'PASSED' if passed else 'FAILED'}")8.7.4 Convergence Testing
The Yee scheme should exhibit second-order convergence. We verify this by running a grid refinement study:
from src.em import convergence_test_maxwell_1d
import numpy as np
grid_sizes, errors, order = convergence_test_maxwell_1d(
grid_sizes=[50, 100, 200, 400],
T=0.5e-9,
CFL=0.5,
wavelength=0.1,
)
print("Grid refinement study:")
print("-" * 40)
for N, err in zip(grid_sizes, errors):
dx = 1.0 / N
print(f"Nx = {N:3d}, dx = {dx:.4f}, error = {err:.2e}")
print("-" * 40)
print(f"Observed order: {order:.2f}")
print(f"Expected order: 2.00")
print(f"Verification: {'PASSED' if 1.9 <= order <= 2.1 else 'FAILED'}")Expected output:
Grid refinement study:
----------------------------------------
Nx = 50, dx = 0.0200, error = 2.34e-03
Nx = 100, dx = 0.0100, error = 5.85e-04
Nx = 200, dx = 0.0050, error = 1.46e-04
Nx = 400, dx = 0.0025, error = 3.66e-05
----------------------------------------
Observed order: 2.00
Expected order: 2.00
Verification: PASSED8.7.5 Method of Manufactured Solutions
MMS provides a systematic way to verify that the code solves the equations correctly, independent of the existence of analytical solutions.
We choose a smooth manufactured solution: \[ E_z^{mms}(x, t) = \sin(\pi x/L) \cos(\omega t) e^{-\alpha t} \tag{8.42}\]
This satisfies a modified PDE with a source term \(f(x, t)\) that we compute symbolically. By running the solver with this source and comparing to the manufactured solution, we verify correctness.
from src.em.verification import manufactured_solution_1d
import numpy as np
# Generate manufactured solution at t = 0.5e-9
x = np.linspace(0, 1, 201)
t = 0.5e-9
E_mms, H_mms, source = manufactured_solution_1d(
x, t,
omega=2*np.pi*1e9,
alpha=1e8,
)
# The source term is what we'd need to add to the RHS
# to make E_mms the exact solution
print(f"Max manufactured E: {np.max(np.abs(E_mms)):.3e}")
print(f"Max source term: {np.max(np.abs(source)):.3e}")8.7.6 Energy Conservation
In lossless media, the total electromagnetic energy should be constant: \[ U = \int \left(\frac{1}{2}\varepsilon E^2 + \frac{1}{2}\mu H^2\right) dx \tag{8.43}\]
from src.em.verification import verify_energy_conservation
from src.em import solve_maxwell_1d, EMConstants
import numpy as np
const = EMConstants()
# Run with PEC boundaries (closed system)
result = solve_maxwell_1d(
L=1.0, Nx=200, T=5e-9, CFL=0.9,
E_init=lambda x: np.exp(-((x - 0.5)**2) / 0.02**2),
bc_left="pec", bc_right="pec",
save_history=True,
)
passed, max_change, energy = verify_energy_conservation(
result.E_history,
result.H_history,
result.dx,
const.eps0,
const.mu0,
tolerance=0.01,
)
print(f"Initial energy: {energy[0]:.6e} J")
print(f"Final energy: {energy[-1]:.6e} J")
print(f"Max relative change: {100*max_change:.4f}%")
print(f"Verification: {'PASSED' if passed else 'FAILED'}")8.7.7 Verification Against Published Results
For additional confidence, we compare against published benchmarks.
Monk-S"uli Convergence Rates (Monk and Süli 1994):
The seminal paper by Monk and S"uli proves that the Yee scheme achieves second-order convergence in the \(L^2\) norm. Our tests should reproduce their Table 2 results.
from src.em.verification import verify_monk_suli_convergence
def test_function(N):
"""Run solver and return (error, dx)."""
result = solve_maxwell_1d(
L=1.0, Nx=N, T=0.5e-9, CFL=0.5,
E_init=lambda x: np.sin(np.pi * x),
bc_left="pec", bc_right="pec",
)
# Compare with standing wave solution
E_exact = np.sin(np.pi * result.x_E) * np.cos(np.pi * const.c0 * result.t)
error = np.sqrt(np.mean((result.E_z - E_exact)**2))
return error, result.dx
passed, order, errors = verify_monk_suli_convergence(
test_function,
grid_sizes=[25, 50, 100, 200],
expected_order=2.0,
tolerance=0.2,
)
print(f"Observed order: {order:.2f}")
print(f"Monk-Suli expected: 2.0")
print(f"Verification: {'PASSED' if passed else 'FAILED'}")Taflove Dispersion Formula (Taflove 1980):
We verify that our numerical dispersion matches the theoretical formula:
from src.em.verification import taflove_dispersion_formula
from src.em.analysis import numerical_dispersion_relation_1d
import numpy as np
# Test parameters
c = EMConstants().c0
dx = 0.01 # 1 cm grid
dt = 0.9 * dx / c # CFL = 0.9
# Test at various wavenumbers
for k_dx in [0.1, 0.5, 1.0, 2.0]:
k = k_dx / dx
omega = k * c # Physical frequency
# Our implementation
omega_num = numerical_dispersion_relation_1d(k, c, dx, dt)
# Taflove's formula
ratio_taflove = taflove_dispersion_formula(omega, c, dx, dt)
print(f"k*dx = {k_dx:.1f}: omega_num/omega = {omega_num/omega:.6f}, "
f"Taflove = {ratio_taflove:.6f}")8.7.8 Summary of Verification Tests
| Test | What it Verifies | Acceptance Criterion |
|---|---|---|
| Plane wave | Correct wave speed | Error < 1% |
| PEC reflection | Boundary conditions | \(E_z = 0\) at boundary |
| Impedance ratio | \(E\)-\(H\) relationship | \(\|E\|/\|H\| = \eta \pm 1\%\) |
| Convergence | Order of accuracy | Rate = \(2.0 \pm 0.1\) |
| Energy | Conservation | Change < 1% |
| Monk-S"uli | Published results | Rate matches paper |
Our test suite in tests/test_maxwell*.py implements all these checks, ensuring the solver remains correct as the code evolves.
8.8 Application: Dielectric Waveguides
Dielectric waveguides confine and guide electromagnetic waves through total internal reflection. They are fundamental to optical fibers, integrated photonics, and laser resonators. This section demonstrates how FDTD can simulate waveguide modes and compare with analytical solutions.
8.8.1 Waveguide Physics
A dielectric slab waveguide consists of a high-index core surrounded by lower-index cladding:
n_clad (cladding)
========================
n_core (core) d = core thickness
========================
n_clad (cladding)Guided modes exist when light undergoes total internal reflection at the core-cladding interface. This requires:
- \(n_{core} > n_{clad}\) (core has higher refractive index)
- Incidence angle exceeds the critical angle: \(\theta > \theta_c = \arcsin(n_{clad}/n_{core})\)
8.8.2 Analytical Mode Solutions
For TE modes in a symmetric slab waveguide, the eigenvalue equation is: \[ \tan\left(\frac{k_x d}{2}\right) = \frac{\gamma}{k_x} \quad \text{(symmetric modes)} \tag{8.44}\] \[ \cot\left(\frac{k_x d}{2}\right) = -\frac{\gamma}{k_x} \quad \text{(antisymmetric modes)} \]
where:
- \(k_x = k_0\sqrt{n_{core}^2 - n_{eff}^2}\) is the transverse wavenumber in core
- \(\gamma = k_0\sqrt{n_{eff}^2 - n_{clad}^2}\) is the decay constant in cladding
- \(k_0 = 2\pi/\lambda_0\) is the free-space wavenumber
- \(n_{eff}\) is the effective index of the mode
The effective index satisfies \(n_{clad} < n_{eff} < n_{core}\). Waves with \(n_{eff}\) below this range are not guided (radiation modes).
8.8.3 Using the Waveguide Module
The src.em.waveguide module provides tools for computing analytical modes:
from src.em import SlabWaveguide, cutoff_wavelength
import numpy as np
# Silicon nitride waveguide in silica
waveguide = SlabWaveguide(
n_core=2.0, # Silicon nitride
n_clad=1.45, # Silica cladding
thickness=0.4e-6, # 400 nm core
wavelength=1.55e-6 # Telecom wavelength
)
print(f"V-number: {waveguide.V:.2f}")
print(f"Max modes: ~{waveguide.max_modes}")
# Find guided modes
modes = waveguide.find_modes()
for mode in modes:
print(f"Mode {mode.mode_number}: n_eff = {mode.n_eff:.4f}, "
f"type = {mode.symmetry}")8.8.4 Mode Profiles
Each mode has a characteristic transverse field profile:
import matplotlib.pyplot as plt
import numpy as np
# Coordinate across the waveguide
x = np.linspace(-1e-6, 1e-6, 500) # +/-1 micron
fig, ax = plt.subplots(figsize=(8, 5))
for mode in modes[:3]: # First three modes
profile = waveguide.mode_profile(mode, x)
ax.plot(x * 1e6, profile, label=f"Mode {mode.mode_number}")
# Show core region
ax.axvline(-waveguide.thickness/2 * 1e6, color='gray', linestyle='--')
ax.axvline(waveguide.thickness/2 * 1e6, color='gray', linestyle='--')
ax.fill_betweenx([-1.5, 1.5],
-waveguide.thickness/2 * 1e6,
waveguide.thickness/2 * 1e6,
alpha=0.2, label='Core')
ax.set_xlabel('Position [um]')
ax.set_ylabel('E-field (normalized)')
ax.set_title('Waveguide Mode Profiles')
ax.legend()
ax.set_xlim(-1, 1)Key observations:
- Fundamental mode (mode 0): Symmetric, single peak in core
- Higher-order modes: More oscillations, wider extent
- Evanescent tails: Field decays exponentially in cladding
8.8.5 FDTD Simulation of Waveguide
We can verify the analytical solutions using FDTD. The approach:
- Create a 2D grid with the waveguide structure
- Launch a wave at one end
- Observe mode confinement and propagation
- Measure effective index from propagation phase
from src.em import solve_maxwell_2d
from src.em.materials import DielectricMaterial, create_layered_model
import numpy as np
# Waveguide parameters (scaled for FDTD)
n_core = 2.0
n_clad = 1.45
core_thickness = 0.4e-6
wavelength = 1.55e-6
# Grid setup
Lx = 10e-6 # 10 um propagation length
Ly = 3e-6 # 3 um transverse
Nx = 400
Ny = 120
# Create material model
eps_r = np.ones((Nx + 1, Ny + 1)) * n_clad**2
# Core region
y = np.linspace(0, Ly, Ny + 1)
core_mask = np.abs(y - Ly/2) < core_thickness/2
for i in range(Nx + 1):
eps_r[i, core_mask] = n_core**2
# Initial condition: Gaussian beam at input
def E_init(X, Y):
sigma_y = core_thickness / 2
return np.exp(-((Y - Ly/2)**2) / (2*sigma_y**2)) * \
np.exp(-((X - 0.5e-6)**2) / (0.2e-6)**2)
# Run simulation
result = solve_maxwell_2d(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny,
T=50e-15, # 50 fs
CFL=0.5,
eps_r=eps_r,
E_init=E_init,
pml_width=15,
save_history=True,
save_every=5,
)8.8.6 Measuring Effective Index
The effective index can be measured from the FDTD simulation by tracking the phase of the propagating mode:
# Extract field along propagation direction at core center
y_center = Ny // 2
E_along_x = result.E_z[:, y_center]
# Find phase velocity from interference pattern
# (requires CW excitation or Fourier analysis of pulse)
# Alternative: Compare with analytical mode profile
mode = modes[0] # Fundamental mode
print(f"Analytical n_eff: {mode.n_eff:.4f}")
print(f"Confinement factor: {waveguide.confinement_factor(mode):.2%}")8.8.7 Single-Mode Condition
For telecommunications, single-mode operation is often desired. The cutoff condition is given by the V-number: \[ V = \frac{\pi d}{\lambda_0} \sqrt{n_{core}^2 - n_{clad}^2} < \frac{\pi}{2} \tag{8.45}\]
from src.em import single_mode_condition
# Maximum core thickness for single-mode at 1.55 um
d_max = single_mode_condition(n_core=2.0, n_clad=1.45, wavelength=1.55e-6)
print(f"Max single-mode thickness: {d_max*1e9:.1f} nm")
# Cutoff wavelength for our 400 nm core
lambda_c = cutoff_wavelength(n_core=2.0, n_clad=1.45,
thickness=0.4e-6, mode_number=1)
print(f"First higher-order mode cutoff: {lambda_c*1e9:.1f} nm")8.8.8 Bending Loss
When a waveguide bends, the outer edge must travel faster than the speed of light in the cladding, causing radiation loss. FDTD naturally captures this effect by including curved waveguide geometries.
8.8.9 Summary
This application demonstrates:
- Analytical foundation: Eigenvalue equations for guided modes
- FDTD verification: Numerical simulation confirms analytical results
- Practical parameters: Effective index, confinement factor, cutoff
- Physical insight: Mode profiles show field confinement
The waveguide analysis techniques extend to more complex structures like photonic crystal waveguides, ring resonators, and directional couplers, all of which can be simulated with FDTD.
8.9 Application: Ground Penetrating Radar
Ground Penetrating Radar (GPR) is a geophysical method that uses electromagnetic pulses to image subsurface structures. FDTD is widely used for GPR simulation, enabling interpretation of field data and system design.
8.9.1 GPR Fundamentals
A GPR system consists of:
- Transmitter: Emits short EM pulses (typically 100 MHz - 2 GHz)
- Receiver: Records reflected signals
- Controller: Triggers pulses and records data
The transmitted pulse propagates into the ground and reflects from:
- Soil layer interfaces
- Buried objects (pipes, voids, artifacts)
- Water table
- Bedrock
8.9.2 EM Properties of Soil
Unlike free space, soil is a lossy dielectric with:
- Relative permittivity \(\varepsilon_r\): 3-40 (depends on water content)
- Conductivity \(\sigma\): 0.001-0.1 S/m (causes attenuation)
The wave velocity in soil is: \[ v = \frac{c}{\sqrt{\varepsilon_r}} \approx 0.1-0.15 \text{ m/ns (typical)} \]
Topp’s empirical equation relates permittivity to volumetric water content \(\theta\): \[ \varepsilon_r = 3.03 + 9.3\theta + 146\theta^2 - 76.7\theta^3 \tag{8.46}\]
from src.em.materials import topp_equation
import numpy as np
water_content = np.linspace(0, 0.5, 50)
eps_r = [topp_equation(w) for w in water_content]
# At 20% water content
print(f"eps_r at 20% water: {topp_equation(0.2):.1f}")8.9.3 Attenuation in Lossy Media
The electric field attenuates as: \[ E(z) = E_0 e^{-\alpha z} \]
where the attenuation coefficient is: \[ \alpha = \omega\sqrt{\frac{\mu\varepsilon}{2}\left(\sqrt{1 + \left(\frac{\sigma}{\omega\varepsilon}\right)^2} - 1\right)} \]
The skin depth \(\delta = 1/\alpha\) is the depth at which the field decays to \(1/e\) of its surface value.
from src.em.materials import DielectricMaterial
# Wet clay
wet_clay = DielectricMaterial(name="Wet clay", eps_r=25, sigma=0.05)
skin_depth = wet_clay.skin_depth(frequency=500e6)
print(f"Skin depth at 500 MHz: {skin_depth:.2f} m")8.9.4 GPR Source Wavelets
GPR transmitters emit short pulses. The Ricker wavelet (Mexican hat) is commonly used:
from src.em.gpr import ricker_wavelet, wavelet_spectrum
import numpy as np
t = np.linspace(0, 20e-9, 1000) # 20 ns window
f0 = 500e6 # 500 MHz center frequency
wavelet = ricker_wavelet(t, f0=f0)
# Analyze spectrum
freq, spectrum = wavelet_spectrum(wavelet, t[1] - t[0])Other common wavelets include:
- Gaussian derivative: Broader bandwidth
- Blackman-Harris: Lower sidelobes
8.9.5 1D GPR Simulation
A simple 1D simulation models vertical propagation:
from src.em.gpr import run_gpr_1d, two_way_travel_time
from src.em.materials import DRY_SAND
# Simulate GPR over dry sand with a buried reflector
result = run_gpr_1d(
depth=2.0, # 2 m total depth
eps_r_soil=4.0, # Dry sand
sigma_soil=0.001, # Low loss
frequency=500e6, # 500 MHz
target_depth=1.0, # Reflector at 1 m
target_eps_r=1.0, # Air void
)
# Expected travel time to 1 m depth
expected_twtt = two_way_travel_time(depth=1.0, eps_r=4.0)
print(f"Expected TWTT: {expected_twtt*1e9:.2f} ns")8.9.6 B-Scan Radargram
A B-scan is a 2D image showing reflections along a survey line:
- Horizontal axis: Antenna position
- Vertical axis: Two-way travel time (depth)
- Color/intensity: Reflection amplitude
Buried objects create characteristic hyperbolic diffraction patterns.
from src.em.gpr import run_gpr_bscan_2d, hyperbola_travel_time
import matplotlib.pyplot as plt
import numpy as np
# Simulate B-scan over buried pipe
result = run_gpr_bscan_2d(
Lx=2.0, Ly=1.5, # 2 m survey line, 1.5 m depth
eps_r_background=9.0, # Wet sand
sigma_background=0.02,
frequency=500e6,
n_traces=50,
target_center=(1.0, 0.5), # Pipe at center, 50 cm depth
target_radius=0.05, # 10 cm diameter pipe
)
if result.bscan is not None:
plt.figure(figsize=(10, 6))
plt.imshow(result.bscan, aspect='auto',
extent=[result.positions[0], result.positions[-1],
result.t[-1]*1e9, 0],
cmap='gray')
plt.xlabel('Position [m]')
plt.ylabel('Two-way travel time [ns]')
plt.title('GPR B-scan')
plt.colorbar(label='Amplitude')8.9.7 Hyperbolic Reflections
A point scatterer creates a hyperbolic pattern because the travel time is: \[ t(x) = \frac{2}{v}\sqrt{(x - x_0)^2 + z_0^2} \]
where \((x_0, z_0)\) is the scatterer position. The apex of the hyperbola indicates the target location; the shape can be used to estimate velocity.
from src.em.gpr import fit_hyperbola
import numpy as np
# If we extract the hyperbola apex times from the B-scan:
# x_positions = np.linspace(0.5, 1.5, 20)
# travel_times = [pick times from radargram]
# x0, z0, v = fit_hyperbola(x_positions, travel_times)8.9.8 Comparison with gprMax
gprMax (Warren, Giannopoulos, and Giannakis 2016) is a widely-used open-source GPR simulator. Our implementation can be validated against gprMax results:
# Key benchmark: Buried cylinder B-scan
# Compare:
# 1. Direct wave arrival time
# 2. Reflection hyperbola shape
# 3. Peak amplitude ratios
# Expected results from gprMax:
# - Direct wave at t = d_antenna / c
# - Hyperbola apex at t = 2*depth / v_soil
# - Amplitude ratio depends on reflection coefficient8.9.9 Layered Earth Model
Real soils often have distinct layers:
from src.em.materials import create_layered_model, DielectricMaterial
# Three-layer model: topsoil, clay, sand
layers = [
(0.3, DielectricMaterial("Topsoil", eps_r=12, sigma=0.01)),
(0.5, DielectricMaterial("Clay", eps_r=25, sigma=0.05)),
(1.2, DielectricMaterial("Sand", eps_r=5, sigma=0.001)),
]
eps_r, sigma = create_layered_model(layers, Nx=200, L=2.0)Each interface produces a reflection with amplitude determined by the Fresnel coefficient: \[ R = \frac{\sqrt{\varepsilon_1} - \sqrt{\varepsilon_2}}{\sqrt{\varepsilon_1} + \sqrt{\varepsilon_2}} \]
8.9.10 Practical Considerations
Frequency Selection:
| Frequency | Resolution | Depth Penetration |
|---|---|---|
| 100 MHz | ~1 m | 10-30 m |
| 500 MHz | ~0.2 m | 2-5 m |
| 1 GHz | ~0.1 m | 0.5-2 m |
| 2 GHz | ~0.05 m | <0.5 m |
Higher frequencies provide better resolution but attenuate faster.
Clutter and Noise:
Real GPR data contains:
- Surface reflections (antenna-ground interface)
- Multiples (repeated bounces)
- Lateral waves
- System ringing
FDTD simulations help distinguish target responses from clutter.
8.9.11 Summary
GPR simulation with FDTD enables:
- System design: Optimize frequency, antenna spacing
- Data interpretation: Understand complex radargrams
- Forward modeling: Create synthetic data for inversion
- Training: Generate labeled data for ML-based detection
The src.em.gpr module provides a starting point; production GPR simulation codes like gprMax (Warren, Giannopoulos, and Giannakis 2016) offer additional features like dispersive media, antenna models, and parallel computation.
8.10 Exercises
8.10.1 Derivation Exercises
Exercise 8.1 (Derive 1D FDTD Update Equations) Starting from the 1D Maxwell equations 8.9 and 8.10, derive the discrete update equations 8.16 and 8.18. Show that the staggered grid arrangement leads to centered differences in both space and time.
Exercise 8.2 (CFL Condition in 2D) Show that the 2D stability condition 8.27 reduces to the 1D condition when \(\Delta x = \Delta y\) and one of the wavenumber components is zero (wave traveling along a grid axis).
Exercise 8.3 (Numerical Phase Velocity) Starting from the dispersion relation 8.37, derive an expression for the relative phase velocity error as a function of the Courant number \(C\) and points per wavelength \(N_\lambda = \lambda/\Delta x\). Verify that the error vanishes when \(C = 1\).
8.10.2 Implementation Exercises
Exercise 8.4 (Soft Source Implementation) Modify the 1D solver to implement a soft source (additive source term) at an interior point. Verify that waves propagate in both directions from the source. Compare with a hard source (field value directly set).
Exercise 8.5 (Lossy Material Implementation) Add a lossy material (\(\sigma \neq 0\)) to the 1D solver using the update coefficients 8.21. Verify that:
- The wave amplitude decays exponentially with distance
- The decay rate matches the analytical skin depth
Exercise 8.6 (Mur’s First-Order ABC) Implement Mur’s first-order absorbing boundary condition for the 1D solver: \[ E^{n+1}_0 = E^n_1 + \frac{c\Delta t - \Delta x}{c\Delta t + \Delta x}(E^{n+1}_1 - E^n_0) \]
Compare its performance (reflection coefficient) against:
- Simple PEC boundary
- The simple ABC already in the solver
8.10.3 Analysis Exercises
Exercise 8.7 (Dispersion Error Visualization) Plot the phase velocity error as a function of:
- Points per wavelength \(N_\lambda\) (for fixed \(C = 0.9\))
- Courant number \(C\) (for fixed \(N_\lambda = 10\))
Use the src.em.analysis.dispersion_maxwell module. What practical guideline would you recommend for choosing \(N_\lambda\)?
Exercise 8.8 (Magic Time Step Investigation) Numerically verify the “magic time step” property by running a plane wave simulation with \(C = 1\) exactly. Measure the phase error after the wave has traveled 100 wavelengths. Compare with \(C = 0.9\) and \(C = 0.5\).
Exercise 8.9 (Convergence Rate Study) Run a grid refinement study for the 1D Maxwell solver with grid sizes \(N = 50, 100, 200, 400, 800\). Plot the L2 error versus \(\Delta x\) on a log-log scale and verify second-order convergence. What happens at very fine grids (floating-point precision limits)?
8.10.4 Application Exercises
Exercise 8.10 (1D Fabry-P'erot Resonator) Simulate a 1D Fabry-P'erot resonator: two PEC boundaries separated by distance \(L\). Excite the system with a broadband pulse and observe the resonant modes. Verify that the resonant frequencies are \(f_n = nc/(2L)\).
Exercise 8.11 (Reflection at Dielectric Interface) Simulate a plane wave normally incident on an interface between two dielectrics. Measure the reflection and transmission coefficients and compare with the analytical Fresnel equations: \[ R = \left(\frac{\sqrt{\varepsilon_1} - \sqrt{\varepsilon_2}}{\sqrt{\varepsilon_1} + \sqrt{\varepsilon_2}}\right)^2 \]
Test for \(\varepsilon_2/\varepsilon_1 = 2, 4, 9\).
Exercise 8.12 (GPR Depth Estimation) Using the src.em.gpr module, simulate a GPR survey over a buried pipe. From the resulting B-scan, extract the hyperbolic diffraction pattern and use the fit_hyperbola function to estimate:
- The pipe depth
- The soil velocity
Compare with the known values used in the simulation.
8.10.5 Advanced Exercises
Exercise 8.13 (2D Scattering from Cylinder) Simulate 2D scattering of a plane wave from a dielectric cylinder. For a cylinder much larger than the wavelength, the shadow region and diffraction fringes should be visible. Compare qualitatively with Mie theory predictions.
Exercise 8.14 (PML Parameter Optimization) Investigate the effect of PML parameters on reflection:
- Vary the PML width from 5 to 30 cells
- Vary the polynomial order from 2 to 5
- Vary \(\sigma_{\max}\) from \(0.5\times\) to \(2\times\) the optimal value
Plot the maximum reflection coefficient versus each parameter and determine the optimal configuration.
Exercise 8.15 (Waveguide Bend Loss) Using 2D FDTD, simulate a \(90^\circ\) bend in a dielectric waveguide. Measure the transmitted power as a function of bend radius. At what radius does the loss become acceptable (<3 dB)?
8.10.6 Project Exercises
Exercise 8.16 (Simple Dipole Antenna) Extend the 2D solver to model a simple dipole antenna:
- Use a line current source
- Compute the near-field pattern
- Use a large PML-bounded domain to approximate far-field conditions
- Extract the radiation pattern (field vs. angle)
Exercise 8.17 (Photonic Crystal Waveguide) Design a 2D photonic crystal by arranging dielectric cylinders in a triangular lattice. Introduce a line defect to create a waveguide. Demonstrate wave guiding through the defect channel.
Exercise 8.18 (GPR with Machine Learning) Generate a dataset of GPR B-scans using the simulator with random buried target positions and depths. Train a simple neural network to detect and locate targets from the B-scan images. Evaluate detection accuracy and localization error.