PML

In numerical electromagnetic simulations, one of the most significant challenges is simulating an infinite or unbounded domain within a finite computational grid. When electromagnetic waves reach the edge of a truncated domain, they typically undergo artificial reflections, which contaminate the solution and degrade accuracy. The Perfectly Matched Layer (PML) is a sophisticated absorbing boundary condition designed to mitigate this issue.

By constructing an artificial, lossy medium surrounding the primary computational domain, the PML allows waves to enter the layer and undergo exponential decay without significant reflection at the interface.

1. Complex Coordinate Stretching

The mathematical core of the PML lies in the concept of complex coordinate stretching. Instead of using real-valued spatial coordinates, the PML effectively transforms the wave vector $\mathbf{k}$ into a complex-valued vector. For each spatial dimension $i \in {x, y, z}$, the transformation is expressed as:

$$\tilde{k}_i = \frac{k_i}{1 + j\sigma_i / \omega}$$

where $\sigma_i$ represents the absorption coefficient in the $i$-th direction and $\omega$ is the angular frequency. This transformation ensures that the wave impedance remains perfectly matched to the adjacent medium at the interface, theoretically eliminating reflections regardless of the angle of incidence or frequency.

2. Exponential Attenuation

As the electromagnetic wave propagates through the PML, its amplitude decays according to the relationship $\exp(-\alpha d)$, where $d$ is the depth into the layer and $\alpha$ is proportional to the absorption coefficients $\sigma_i$. If the absorption profile is optimized, the wave energy can be attenuated below the machine's floating-point precision within a very small number of grid cells.

3. Equivalence to Anisotropic Media

From a physical perspective, the coordinate stretching can be mathematically reformulated as the introduction of an anisotropic, complex conductivity tensor within the PML region. This allows the standard Maxwell’s equations to maintain their original form, provided that the material parameters (permittivity $\varepsilon$ and permeability $\mu$) are replaced by complex, direction-dependent tensors.

Comparative Analysis of PML Implementations

Different numerical solvers require different flavors of PML to balance accuracy, stability, and computational overhead.

For the purpose of this technical guide, we will focus on the UPML implementation within a 2D TM-mode FDTD framework, as it offers an optimal balance between numerical stability and implementation simplicity.

Implementation Workflow for UPML in FDTD

Implementing a robust UPML requires a systematic approach to parameter definition and field updates.

Step 1: Determining Layer Thickness

The thickness of the PML (number of grid cells) is a critical trade-off. A common practice is to use 8 to 12 cells. While increasing the thickness improves absorption, it also increases the total memory footprint and computational time.

Step 2: Constructing the Absorption Profile

To minimize reflections caused by the sudden change in material properties at the interface, the absorption coefficient $\sigma$ should not be a constant. Instead, it should follow a polynomial growth profile:

$$\sigma_i (n) = \sigma_{\max}\left(\frac{n}{N_{\text{PML}}}\right)^m$$

where $n$ is the cell index within the PML, $N_{\text{PML}}$ is the total thickness, and $m$ is the polynomial order (typically $m=3$ or $4$). The maximum conductivity $\sigma_{\max}$ can be estimated using the empirical formula:
$$\sigma_{\max} = \frac{m+1}{150\pi\Delta}$$
where $\Delta$ is the grid spacing.

Step 3: Calculating Update Coefficients

In the UPML framework, we introduce auxiliary parameters $\kappa$ and $\alpha$ to modify the field updates. For a given direction $i$:

• $\kappa_i$ controls the scaling of the field components.
• $\alpha_i$ manages the damping.

These parameters are pre-calculated based on the desired frequency range and the polynomial profile of $\sigma$.

Step 4: Modifying the Update Equations

In a standard FDTD Yee cell, the $E_z$ field is updated using the curl of $H$. In the presence of UPML, we introduce convolutional auxiliary variables ($\psi$) to account for the lossy, anisotropic nature of the layer. The modified update for $E_z$ becomes:

$$\psi_{E_x}^{n+1} = b_{E_x}\psi_{E_x}^{n} + c_{E_x}(\text{curl } H)$$
$$E_z^{n+1} = a_{E_z}E_z^{n} + d_{E_z}(\psi_{E_x}^{n+1} - \psi_{E_y}^{n+1})$$

The coefficients $a, b, c,$ and $d$ are derived from the $\sigma, \kappa,$ and $\alpha$ values.

MATLAB Implementation Example (2D TM-FDTD + UPML)

The following code demonstrates a simplified 2D TM-mode simulation using a UPML-inspired approach.

% -------------------------------------------------
% 2D TM FDTD with UPML-inspired Absorption
% -------------------------------------------------
clear; clc;

% ---------- Physical Constants ----------
c0   = 299792458;               
eps0 = 8.854187817e-12;
mu0  = 4*pi*1e-7;
f0   = 3e9;                     
lambda0 = c0/f0;
dx = dy = lambda0/40;          
dt = 0.99/(c0*sqrt(1/dx^2+1/dy^2)); 

% ---------- Grid Setup ----------
Nx = 200; Ny = 200;            
npml = 12;                     

% ---------- Material Initialization ----------
eps = eps0 * ones(Nx, Ny);
mu  = mu0 * ones(Nx, Ny);

% ---------- PML Parameter Setup ----------
m = 4;                         
sigma_max = (m+1)/(150*pi*dx);
sigma_x = zeros(Nx,1);
sigma_y = zeros(Ny,1);

for i=1:npml
    % Apply polynomial profile to boundaries
    val = sigma_max * ((npml-i+1)/npml)^m;
    sigma_x(i) = val;
    sigma_x(Nx-i+1) = val;
    sigma_y(i) = val;
    sigma_y(Ny-i+1) = val;
end

% Expand sigma to 2D matrices
Sigma_x = repmat(sigma_x, 1, Ny);
Sigma_y = repmat(sigma_y.', Nx, 1);

% ---------- Pre-calculate Update Coefficients (Ez Example) ----------
% These coefficients incorporate the damping effect
aEz = (1 - dt*Sigma_x./(2*eps))./(1 + dt*Sigma_x./(2*eps));
bEz = dt./(eps.*dx)./(1 + dt*Sigma_x./(2*eps));
cEz = dt./(eps.*dy)./(1 + dt*Sigma_y./(2*eps));

% ---------- Field Initialization ----------
Ez = zeros(Nx, Ny);
Hx = zeros(Nx, Ny-1);
Hy = zeros(Nx-1, Ny);

% ---------- Time Stepping Loop ----------
nSteps = 800;
for n = 1:nSteps
    % 1. Update Magnetic Fields (H)
    Hx(:,1:end) = Hx(:,1:end) - (dt/mu(:,1:end)).*(Ez(:,2:end)-Ez(:,1:end-1))/dy;
    Hy(1:end,:) = Hy(1:end,:) + (dt/mu(1:end,:)).*(Ez(2:end,:)-Ez(1:end-1,:))/dx;
    
    % 2. Update Electric Field (Ez) with PML coefficients
    curlH = (Hy(2:end,:)-Hy(1:end-1,:))/dx - (Hx(:,2:end)-Hx(:,1:end-1))/dy;
    % Apply the modified update equation to the interior and PML regions
    Ez(2:end-1,2:end-1) = aEz(2:end-1,2:end-1).*Ez(2:end-1,2:end-1) ...
                         + bEz(2:end-1,2:end-1).*curlH;
    
    % 3. Soft Source (Gaussian Pulse)
    t0 = 40; spread = 15;
    Ez(Nx/2, Ny/2) = Ez(Nx/2, Ny/2) + exp(-((n-t0)/spread)^2);
    
    % 4. Visualization
    if mod(n,20)==0
        imagesc(Ez.');
        colorbar; axis equal tight;
        title(['Ez Field Distribution at Step ', num2str(n)]);
        drawnow;
    end
end

Parameter Optimization and Tuning Strategies

To achieve high-fidelity results, one must fine-tune several parameters:

1. PML Thickness vs. Frequency: For low-frequency problems (long wavelengths), a thicker PML (15+ layers) is recommended to ensure sufficient attenuation. For high-frequency problems, 8–10 layers are usually sufficient.
2. Polynomial Order ($m$):
   • A lower order ($m=3$) provides a smoother transition, which is better for broadband stability.
   • A higher order ($m=4$ or $5$) allows for more aggressive absorption in a shorter distance but can introduce numerical reflections if the grid is not fine enough.
3. Frequency Considerations: The $\sigma_{\max}$ must be calculated based on the highest frequency component in the simulation. If the absorption is too low for high-frequency waves, they will "leak" through the boundary.
4. Numerical Stability: Extreme absorption coefficients can lead to numerical instability. If the simulation crashes, consider reducing the time step dt (e.g., using a 0.95 CFL factor instead of 0.99).
5. Validation: Always validate your PML by simulating a plane wave incident on the boundary. The goal is a reflection coefficient $R = |E_{\text{ref}}/E_{\text{inc}}|$ of less than $10^{-3}$ (approximately -60 dB).

Summary

The Perfectly Matched Layer is an indispensable tool in modern computational electromagnetics. By leveraging the principle of complex coordinate stretching, it creates a boundary that is theoretically reflectionless and highly absorptive. Whether using the classic Berenger method or the more modern UPML and CPML approaches, the key to success lies in the careful selection of the absorption profile, layer thickness, and frequency-dependent coefficients. Mastering these nuances allows researchers to simulate complex, open-space electromagnetic phenomena with unprecedented accuracy and efficiency.