Review: Perfectly Matched Layers Implementation & Stability

Jul 12, 2026

1. Brief Intro: PML

(1) Bérenger’s PML with split fields

Wave Impedance of Lossy Materials

When the magnetic loss exists, the characteristic impedance can be expressed as $\newcommand{\dif}{\mathrm{d}} \eta = \sqrt{\frac{\tilde{\mu}}{\tilde{\epsilon}}}=\sqrt{\frac{\epsilon-i\sigma_e/\omega}{\mu-i\sigma_m/\omega}}$ which degenerates to the lossless expression when $\sigma_m/\mu=\sigma_e/\epsilon$ .

Proof. From the frequency domain Maxwell equation $\mathbf{\nabla} \times \mathbf{E} = -i\omega\mu \mathbf{H} -\sigma_m\mathbf{H}=-i\omega\tilde{\mu}\mathbf{H},\quad \mathbf{\nabla}\times \mathbf{H}=i\omega\epsilon \mathbf{E}+\sigma_e\mathbf{E}=i\omega\tilde{\epsilon}\mathbf{E}$ The characteristic impedance is defined as the magnitude ratio of the electric and magnetic fields of the wave travelling through free space $\eta = \frac{|\mathbf{E}|}{|\mathbf{H}|} = \sqrt{\frac{\tilde{\mu}}{\tilde{\epsilon}}}$

For simplicity, we only deal with 2D transverse electric ( $$E_z =0$$ ) Maxwell equation. The simplified dimensionless Maxwell equations in a lossy medium are

\frac{\partial E_y}{\partial t} + \sigma_1 E_y = - \frac{\partial H_z}{\partial x},\quad \frac{\partial E_x}{\partial t} + \sigma_2 E_x = \frac{\partial H_z}{\partial y},\quad \frac{\partial H_z}{\partial t} + \rho^* H_z = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x}

The Bérenger’s idea is that the orthogonal magnetic field $$H_z$$ can be split into nonphysical components $$H_x+H_y$$ . The system becomes

\frac{\partial E_y}{\partial t} + \sigma_1 E_y =-\frac{\partial (H_x+H_y)}{\partial x},\quad \frac{\partial E_x}{\partial t}+\sigma_2 E_x = \frac{\partial (H_x+H_y)}{\partial y}, \quad \frac{\partial H_x}{\partial t}+\sigma_1^* H_x = -\frac{\partial E_y}{\partial x},\quad \frac{\partial H_y}{\partial t}+\sigma_2^* H_y = \frac{\partial E_x}{\partial y}

Notice here instead of using only one scalar $\rho^*$ , we use anisotropic $\mathbf{\rho}^* = \mathrm{diag}(\sigma_1^*, \sigma_2^*)$ . Let’s then consider a wave whose electric field of magnitude $$E_0$$ forms an angle $\varphi$ with the y-axis. If a plane wave propagates in the PML, then the 4-component field can be expressed as

E_x = -E_0\sin\varphi \,e^{i\omega(t-\alpha x-\beta y)},\;E_y = E_0\cos\varphi\, e^{i\omega(t-\alpha x-\beta y)},\, H_x = H_{x0}\,e^{i\omega (t-\alpha x-\beta y)},\, H_y = H_{y0}\,e^{i\omega(t-\alpha x-\beta y)}

By substituting them back into the equations, we get

E_0\cos\varphi-i\frac{\sigma_1}{\omega} E_0\cos\varphi=\alpha( H_{x0} + H_{y0})

E_0\sin\varphi -i \frac{\sigma_2}{\omega}E_0\sin\varphi = \beta(H_{x0}+H_{y0})

H_{x0} - i\frac{\sigma_1^*}{\omega}H_{x0} = \alpha E_0 \cos\varphi

H_{y0}-i\frac{\sigma_2^*}{\omega}H_{y0} = \beta E_0 \sin\varphi

By substituting $H_{x0},\,H_{y0}$ back into the first two equations, we get

(1-i\frac{\sigma_2}{\omega})\sin\varphi = \beta \left[ \frac{\alpha \cos\varphi}{1-i\sigma_1^*/\omega} + \frac{\beta\sin\varphi}{1-i\sigma_2^*/\omega} \right]

(1-i\frac{\sigma_1}{\omega})\cos\varphi = \alpha \left[ \frac{\alpha \cos\varphi}{1-i\sigma_1^*/\omega} +\frac{\beta\sin\varphi}{1-i\sigma_2^*/\omega} \right]

Therefore, $\alpha,\beta$ satisfy the following identity

\frac{\beta}{\alpha} = \frac{\sin\varphi}{\cos\varphi} \frac{1-i\sigma_2/\omega}{1-i\sigma_1/\omega}

There are two sets of $(\alpha,\beta)$ with opposite signs, which stand for two opposite propagation directions. Only keeping the positive sign, we have

\alpha = \frac{1}{G}\left(1-i\frac{\sigma_1}{\omega}\right)\cos\varphi,\quad \beta = \frac{1}{G}\left(1-i\frac{\sigma_2}{\omega}\right)\sin\varphi

G=\sqrt{w_x\cos^2\varphi + w_y \sin^2\varphi},\quad w_1 = \frac{1-i\sigma_1/\omega}{1-i\sigma_1^*/\omega},\quad w_2=\frac{1-i\sigma_2/\omega}{1-i\sigma_2^*/\omega}

The field can be written as

\psi = \psi_0 \exp i\omega [t-(x\cos\varphi+y\sin\varphi)/G]\cdot\exp[-\sigma_1 x\cos\varphi / G]\cdot\exp[-\sigma_2y\sin\varphi/G]

in which $\psi = \{E_x,E_y,H_x,H_y\}$ . The amplitudes satisfy

H_{x0} = \frac{E_0}{G}w_x\cos^2\varphi,\; H_{y0} = \frac{E_0}{G}w_y\sin^2\varphi\Rightarrow H_0=GE_0

Thus, the impedance is

Z= \frac{E_0}{H_0}=\frac{1}{G}\quad(=1\text{ if } \sigma_{1,2}^*=\sigma_{1,2})

When $$G=1$$ we get the field solutions

\psi = \psi_0\exp i\omega[t-(x\cos\varphi+y\sin\varphi)]\cdot \exp[-\sigma_1 x\cos\varphi]\cdot \exp[-\sigma_2 y \sin\varphi]

which decays exponentially along $$x$$ and $$y$$ . Moreover, as the impedance matches the vacuum impedance, there should be no reflection from the interface.

To better understand the reflection at PML-PML interfaces, we study the inclined incidence. Now we first consider the interface normal to $$x$$ along $$y$$ . Taking two points, A, B, along the interface, the field components satisfy

\frac{\psi_t(B)}{\psi_i(B)}=\frac{\psi_t(A)}{\psi_i(A)},\quad \psi_{i,t}(B)=\psi_{i,t}(A)\exp[-i\omega\,d\sin\theta_{i,t}/G_{i,t}-(\sigma_2)_{i,t}d\sin\theta_{i,t}/G_{i,t}]

As $$d$$ is random, the angles $\theta_{i,t}$ should satisfy

\left(1-i\frac{\sigma_{2,i}}{\omega}\right)\frac{\sin \theta_i}{G_i}=\left(1-i\frac{\sigma_{2,t}}{\omega}\right)\frac{\sin \theta_t}{G_t}

where

G_k = \sqrt{w_{1k}\cos^2\theta_k+w_{2k}\sin^2\theta_k},\quad k = i,t

Along with $\theta_r = \pi-\theta_i$ , it forms the Snell-Descartes law at the PML-PML interface. Furthermore, for the TE case, the reflection is

r_p = \frac{Z_t\cos\theta_t-Z_i\cos\theta_i}{Z_t\cos\theta_t+Z_i\cos\theta_i}=\frac{G_i\cos\theta_t-G_t\cos\theta_i}{G_i\cos\theta_t+G_t\cos\theta_i}

when the two PMLs are impedance matched, i.e. $\sigma_{1,t} = \sigma_{1,t}^*,\,\sigma_{1,i}=\sigma_{1,i}^*,\,\sigma_2=\sigma_2^*\Rightarrow G_i=G_t$ , we finally have $$r_p=0$$ . Similarly, if an interface is normal to $$y$$ along $$x$$ , the reflection is null if two PMLs matched.

Now we can move on to the setup of PMLs. In simulation, a PML is first used to absorb the outgoing wave energy and then the tail is reflected with a PEC. Theoretically, a PML with uniform $\sigma$ can be applied without causing arbitrary small reflection, e.g. when thickness is $\delta$ and the incidence angle is $\theta$

\psi(\rho) = \psi(0) e^{-\sigma\cos\theta\,\rho},\;R(\theta)= e^{-2\sigma\cos\theta\,\delta}

However, numerical reflections will happen due to discretization in practice. Graded conductivity is then used to mitigate numerical reflection, e.g. $\sigma(\rho) = \sigma_m (\rho/\delta)^n$

R(\theta) = \exp\left[-2\cos\theta \int_0^\delta \sigma(\rho)\,\dif\rho\right]=\exp\left[-\frac{2}{n+1}\sigma_m\delta \cos\theta\right]

Now we can move on to the stability analysis. Define $$U = (H_z,E_y,E_x)^T$$ and

\mathrm{E}(t) = \int_{\mathbb{R}^2} \left(\left|\frac{\partial U}{\partial t}\right|^2+\left|\frac{\partial U}{\partial x}\right|^2+\left|\frac{\partial U}{\partial y}\right|^2\right)\dif x\dif y,\quad ||H_x(t)||^2_0=\int_{\mathbb{R}^2}H_x^2\dif x\dif y

The two conditions are called Weak Stability, since the $$L^2$$ -norms of $H_{x,y}$ are not only bounded by the $$L^2$$ -norms of the initial data, but also bounded by that of the derivatives. This property may cause instabilities in numerical computations.

(2) Unsplit-field UPML

Using the modified Maxwell’s equations in $$k$$ -space

\mathbf{k}_e \times \mathbf{E} = \omega \mu\mathbf{ H},\quad \mathbf{k}_h \times \mathbf{H}=-\omega \epsilon\mathbf{ E}

\mathbf{k}_e = \hat{x}\frac{k_x}{e_x}+\hat{y}\frac{k_y}{e_y}+\hat{z}\frac{k_z}{e_z},\quad \mathbf{k}_h = \hat{x}\frac{k_x}{h_x}+\hat{y}\frac{k_y}{h_y}+\hat{z}\frac{k_z}{h_z}

Combining the equation above, we have

-\omega^2\mu\epsilon \mathbf{H} = \mathbf{k}_e\times \mathbf{k}_h \times \mathbf{H}=\mathbf{k}_h(\mathbf{k}_e\cdot\mathbf{H}) - (\mathbf{k}_e\cdot\mathbf{k}_h)\mathbf{H}

As the magnetic divergence equation leads to $\mathbf{k}_e\cdot\mathbf{H}=0$ for homogeneous medium, we finally get the dispersion relation

\omega^2\mu\epsilon = \mathbf{k}_e\cdot \mathbf{k}_h = k^2 = \sum_i\frac{1}{e_i h_i} k_i^2

When perfectly matched (Condition 1 - Wave Impedance Match), i.e., $e_i=h_i\equiv s_i$ , the wave impedance is therefore

\eta = \frac{|\mathbf{E}|}{|\mathbf{H}|} = \frac{|\mathbf{k}_h|}{\omega\epsilon} = \sqrt{\frac{\mu}{\epsilon}}

However, this condition doesn’t necessarily mean no reflection at any angles. Further conditions should be satisfied. Now consider an $$xy$$ -interface separating two PML layers denoted as $$i,t$$ , both of which holds a dispersion relation

\omega^2\mu\epsilon=\frac{k_{ix}^2}{s_{ix}^2}+\frac{k_{iy}^2}{s_{iy}^2}+\frac{k_{iz}^2}{s_{iz}^2}=\frac{k_{tx}^2}{s_{tx}^2}+\frac{k_{ty}^2}{s_{ty}^2}+\frac{k_{tz}^2}{s_{tz}^2}

The phase matching requires the tangential wave vectors should be continuous, i.e. $k_{ix} = k_{tx},\, k_{iy}=k_{ty}$ . In addition, the tangential fields should be continuous as well, i.e.,

\text{TE mode: }\qquad\qquad(\mathrm{E}_i)_{xy} = (\mathrm{E}_t)_{xy},\quad (\mathrm{H}_i)_{xy} = (\mathrm{H}_t)_{xy} \; \Leftrightarrow\; (\mathbf{k}_{e,i}\times \mathbf{E}_i)_{xy} = (\mathbf{k}_{e,t}\times \mathbf{E}_t)_{xy} \;\Rightarrow\; \frac{s_{iz}\mu_i}{k_{iz}} = \frac{s_{tz}\mu_t}{k_{tz}}

\text{TM mode: }\qquad\qquad(\mathrm{H}_i)_{xy} = (\mathrm{H}_t)_{xy},\quad (\mathrm{E}_i)_{xy} = (\mathrm{E}_t)_{xy} \; \Leftrightarrow\; -(\mathbf{k}_{h,i}\times \mathbf{H}_i)_{xy} = -(\mathbf{k}_{h,t}\times \mathbf{H}_t)_{xy} \;\Rightarrow\; \frac{s_{iz}\epsilon_i}{k_{iz}} = \frac{s_{tz}\epsilon_t}{k_{tz}}

Assuming both sides have the same vacuum $\epsilon_i = \epsilon_t = \epsilon_0,\,\mu_i=\mu_t=\mu_0$ , we conclude the previous conditions as the Condition 2 - Arbitrary Angle Match

\begin{cases} k_{ix} = k_{tx},\, k_{iy}=k_{ty},\,\frac{k_{iz}}{s_{iz}} = \frac{k_{tz}}{s_{tz}},\ \frac{k_{ix}^2}{s_{ix}^2}+\frac{k_{iy}^2}{s_{iy}^2}+\frac{k_{iz}^2}{s_{iz}^2}=\frac{k_{tx}^2}{s_{tx}^2}+\frac{k_{ty}^2}{s_{ty}^2}+\frac{k_{tz}^2}{s_{tz}^2} \end{cases} \quad \Rightarrow \quad s_{ix}=s_{tx},\,s_{iy}=s_{ty},\, s_{iz},\,s_{tz}

This indicates that a PML should be stretched along the normal direction, while the other two directions remain unchanged. Around edges and corners where the normal directions are ill-defined, several strategies can be applied

Prioritized Bypass: Simply assign the edges/corners to one PML class with respect to priority e.g. x-PML > y-PML > z-PML
Multi-direction Stretching: For example, at the edge of x-PML $$(s_x,1,1,s_x,1,1)$$ and y-PML $$(1,s_y,1,1,s_y,1)$$ , we can use $$(s_x,s_y,1,s_x,s_y,1)$$ which supports modes decaying along both $$x$$ and $$y$$ direction.

Moreover, the coordinate stretching is equivalent to a transformation of $\epsilon, \mu$ according to the transformation optics.

Transformation Optics Fundamental Principle

The transformation optics describes how Maxwell equations changed under an arbitrary differentiable coordinate transformation. Assuming we employ a coordinate transformation $x'^i = f^{ij}(x_j),\quad \Lambda^i_j = \frac{\partial x'^i}{\partial x^j},\quad J = \det(\Lambda),\quad g_{ij} = \Lambda_i^k\Lambda_j^l g'_{kl},\quad \sqrt{g} = \sqrt{\det(g_{ij})}$ We have the following transformation

Field Components:
- Covariant vectors: $E'_i=(\Lambda^{-1})_i^jE_j,\,H'_i=(\Lambda^{-1})^j_i H_j$
- Contravariant densities: $D'^i = (\det \Lambda)^{-1} \Lambda^i_j D^j,\,B'^i = (\det\Lambda)^{-1} \Lambda^i_j B^j$
Material Properties:
- Permittivity/Permeability: $\epsilon'^{ij}=(\det \Lambda)^{-1}\Lambda^i_k\Lambda^j_q\epsilon^{kq},\,\mu'^{ij}=(\det \Lambda)^{-1}\Lambda^i_k\Lambda^j_q\mu^{kq}$
- Electric/Magnetic conductivity: $\sigma'^{ij}=(\det \Lambda)^{-1}\Lambda^i_k\Lambda^j_q\sigma^{kq},\,\kappa'^{ij}=(\det \Lambda)^{-1}\Lambda^i_k\Lambda^j_q\kappa^{kq}$
Source:
- Current density: $J'^i = (\det \Lambda)^{-1} \Lambda^i_j J^j$
- Charge density: $\rho' = (\det \Lambda)^{-1}\rho$

For the coordinate scaling, i.e. $\mathbf{\Lambda}=\mathrm{diag}(s_x,s_y,s_z)$ , used in UPMLs, we can therefore obtain equivalent expressions in conventional Euclidean space with

\mathbf{\epsilon} = \det(\mathbf{\Lambda})\mathbf{\Lambda}^{-1} \epsilon_0 \mathbf{\Lambda}^{-T} = \epsilon_0 \,\mathrm{diag}(\frac{s_ys_z}{s_x}, \frac{s_xs_z}{s_y}, \frac{s_xs_y}{s_z}),\quad \mathbf{\mu} = \mu_0 \,\mathrm{diag}(\frac{s_ys_z}{s_x}, \frac{s_xs_z}{s_y}, \frac{s_xs_y}{s_z})