For the transverse electric (TE) case, the original Berenger’s Perfectly Matched Layer (PML) paper [1] writes
Propagation of a Plane Wave in an APML Medium
We consider a plane wave of magnitude (\(E_{0},H_{zx0},H_{zy0}\))
and pulsation \(\omega\) propagating in the APML medium with an
angle \(\varphi\) relative to the x axis
(36)\[E_{x} = -E_{0}\sin \varphi \: e^{i\omega \left( t-\alpha x-\beta y\right) }\]
(37)\[E_{y} = E_{0}\cos \varphi \: e^{i\omega \left( t-\alpha x-\beta y\right) }\]
(38)\[H_{zx} = H_{zx0} \: e^{i\omega \left( t-\alpha x-\beta y\right) }\]
(39)\[H_{zy} = H_{zy0} \: e^{i\omega \left( t-\alpha x-\beta y\right) }\]
where \(\alpha\) and \(\beta\) are two complex constants to
be determined.
Introducing Eqs. (36), (37),
(38) and (39)
into Eqs. (31), (32), (33)
and (34) gives
(40)\[\varepsilon _{0}E_{0}\sin \varphi -i\frac{\sigma _{y}}{\omega }E_{0}\sin \varphi = \beta \frac{c_{y}}{c}\left( H_{zx0}+H_{zy0}\right) +i\frac{\overline{\sigma }_{y}}{\omega }\left( H_{zx0}+H_{zy0}\right)\]
(41)\[\varepsilon _{0}E_{0}\cos \varphi -i\frac{\sigma _{x}}{\omega }E_{0}\cos \varphi = \alpha \frac{c_{x}}{c}\left( H_{zx0}+H_{zy0}\right) -i\frac{\overline{\sigma }_{x}}{\omega }\left( H_{zx0}+H_{zy0}\right)\]
(42)\[\mu _{0}H_{zx0}-i\frac{\sigma ^{*}_{x}}{\omega }H_{zx0} = \alpha \frac{c^{*}_{x}}{c}E_{0}\cos \varphi -i\frac{\overline{\sigma }^{*}_{x}}{\omega }E_{0}\cos \varphi\]
(43)\[\mu _{0}H_{zy0}-i\frac{\sigma ^{*}_{y}}{\omega }H_{zy0} = \beta \frac{c^{*}_{y}}{c}E_{0}\sin \varphi +i\frac{\overline{\sigma }^{*}_{y}}{\omega }E_{0}\sin \varphi\]
Defining \(Z=E_{0}/\left( H_{zx0}+H_{zy0}\right)\) and using Eqs. (40)
and (41), we get
(44)\[\beta = \left[ Z\left( \varepsilon _{0}-i\frac{\sigma _{y}}{\omega }\right) \sin \varphi -i\frac{\overline{\sigma }_{y}}{\omega }\right] \frac{c}{c_{y}}\]
(45)\[\alpha = \left[ Z\left( \varepsilon _{0}-i\frac{\sigma _{x}}{\omega }\right) \cos \varphi +i\frac{\overline{\sigma }_{x}}{\omega }\right] \frac{c}{c_{x}}\]
Adding \(H_{zx0}\) and \(H_{zy0}\) from Eqs. (42)
and (43) and substituting the expressions
for \(\alpha\) and \(\beta\) from Eqs. (44)
and (45) yields
\[\begin{split}\begin{aligned}
\frac{1}{Z} & = \frac{Z\left( \varepsilon _{0}-i\frac{\sigma _{x}}{\omega }\right) \cos \varphi \frac{c^{*}_{x}}{c_{x}}+i\frac{\overline{\sigma }_{x}}{\omega }\frac{c^{*}_{x}}{c_{x}}-i\frac{\overline{\sigma }^{*}_{x}}{\omega }}{\mu _{0}-i\frac{\sigma ^{*}_{x}}{\omega }}\cos \varphi \nonumber
\\
& + \frac{Z\left( \varepsilon _{0}-i\frac{\sigma _{y}}{\omega }\right) \sin \varphi \frac{c^{*}_{y}}{c_{y}}-i\frac{\overline{\sigma }_{y}}{\omega }\frac{c^{*}_{y}}{c_{y}}+i\frac{\overline{\sigma }^{*}_{y}}{\omega }}{\mu _{0}-i\frac{\sigma ^{*}_{y}}{\omega }}\sin \varphi
\end{aligned}\end{split}\]
If \(c_{x}=c^{*}_{x}\), \(c_{y}=c^{*}_{y}\), \(\overline{\sigma }_{x}=\overline{\sigma }^{*}_{x}\), \(\overline{\sigma }_{y}=\overline{\sigma }^{*}_{y}\), \(\frac{\sigma _{x}}{\varepsilon _{0}}=\frac{\sigma ^{*}_{x}}{\mu _{0}}\) and \(\frac{\sigma _{y}}{\varepsilon _{0}}=\frac{\sigma ^{*}_{y}}{\mu _{0}}\) then
(46)\[Z = \pm \sqrt{\frac{\mu _{0}}{\varepsilon _{0}}}\]
which is the impedance of vacuum. Hence, like the PML, given some
restrictions on the parameters, the APML does not generate any reflection
at any angle and any frequency. As for the PML, this property is not
retained after discretization, as shown subsequently.
Calling \(\psi\) any component of the field and \(\psi _{0}\)
its magnitude, we get from Eqs. (36), (44),
(45) and (46) that
(47)\[\psi =\psi _{0} \: e^{i\omega \left( t\mp x\cos \varphi /c_{x}\mp y\sin \varphi /c_{y}\right) }e^{-\left( \pm \frac{\sigma _{x}\cos \varphi }{\varepsilon _{0}c_{x}}+\overline{\sigma }_{x}\frac{c}{c_{x}}\right) x} e^{-\left( \pm \frac{\sigma _{y}\sin \varphi }{\varepsilon _{0}c_{y}}+\overline{\sigma }_{y}\frac{c}{c_{y}}\right) y}.\]
We assume that we have an APML layer of thickness \(\delta\) (measured
along \(x\)) and that \(\sigma _{y}=\overline{\sigma }_{y}=0\)
and \(c_{y}=c.\) Using (47), we determine
that the coefficient of reflection given by this layer is
\[\begin{split}\begin{aligned}
R_{\mathrm{APML}}\left( \theta \right) & = e^{-\left( \sigma _{x}\cos \varphi /\varepsilon _{0}c_{x}+\overline{\sigma }_{x}c/c_{x}\right) \delta }e^{-\left( \sigma _{x}\cos \varphi /\varepsilon _{0}c_{x}-\overline{\sigma }_{x}c/c_{x}\right) \delta },\nonumber
\\
& = e^{-2\left( \sigma _{x}\cos \varphi /\varepsilon _{0}c_{x}\right) \delta },
\end{aligned}\end{split}\]
which happens to be the same as the PML theoretical coefficient of
reflection if we assume \(c_{x}=c\). Hence, it follows that for
the purpose of wave absorption, the term \(\overline{\sigma }_{x}\)
seems to be of no interest. However, although this conclusion is true
at the infinitesimal limit, it does not hold for the discretized counterpart.
Discretization
In the following we set \(\varepsilon_0 = \mu_0 = 1\). We discretize Eqs. (26), (27), (28), and (29) to obtain
\[\frac{E_x|^{n+1}_{j+1/2,k,l}-E_x|^{n}_{j+1/2,k,l}}{\Delta t} + \sigma_y \frac{E_x|^{n+1}_{j+1/2,k,l}+E_x|^{n}_{j+1/2,k,l}}{2} = \frac{H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}}{\Delta y}\]
\[\frac{E_y|^{n+1}_{j,k+1/2,l}-E_y|^{n}_{j,k+1/2,l}}{\Delta t} + \sigma_x \frac{E_y|^{n+1}_{j,k+1/2,l}+E_y|^{n}_{j,k+1/2,l}}{2} = - \frac{H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}}{\Delta x}\]
\[\frac{H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l}-H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l}}{\Delta t} + \sigma^*_x \frac{H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l}+H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l}}{2} = - \frac{E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}}{\Delta x}\]
\[\frac{H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l}-H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l}}{\Delta t} + \sigma^*_y \frac{H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l}+H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l}}{2} = \frac{E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}}{\Delta y}\]
and this can be solved to obtain the following leapfrog integration equations
\[\begin{split}\begin{aligned}
E_x|^{n+1}_{j+1/2,k,l} & = \left(\frac{1-\sigma_y \Delta t/2}{1+\sigma_y \Delta t/2}\right) E_x|^{n}_{j+1/2,k,l} + \frac{\Delta t/\Delta y}{1+\sigma_y \Delta t/2} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}\right)
\\
E_y|^{n+1}_{j,k+1/2,l} & = \left(\frac{1-\sigma_x \Delta t/2}{1+\sigma_x \Delta t/2}\right) E_y|^{n}_{j,k+1/2,l} - \frac{\Delta t/\Delta x}{1+\sigma_x \Delta t/2} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}\right)
\\
H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l} & = \left(\frac{1-\sigma^*_x \Delta t/2}{1+\sigma^*_x \Delta t/2}\right) H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l} - \frac{\Delta t/\Delta x}{1+\sigma^*_x \Delta t/2} \left(E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}\right)
\\
H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l} & = \left(\frac{1-\sigma^*_y \Delta t/2}{1+\sigma^*_y \Delta t/2}\right) H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l} + \frac{\Delta t/\Delta y}{1+\sigma^*_y \Delta t/2} \left(E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}\right)
\end{aligned}\end{split}\]
If we account for higher order \(\Delta t\) terms, a better approximation is given by
\[\begin{split}\begin{aligned}
E_x|^{n+1}_{j+1/2,k,l} & = e^{-\sigma_y\Delta t} E_x|^{n}_{j+1/2,k,l} + \frac{1-e^{-\sigma_y\Delta t}}{\sigma_y \Delta y} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}\right)
\\
E_y|^{n+1}_{j,k+1/2,l} & = e^{-\sigma_x\Delta t} E_y|^{n}_{j,k+1/2,l} - \frac{1-e^{-\sigma_x\Delta t}}{\sigma_x \Delta x} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}\right)
\\
H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_x\Delta t} H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l} - \frac{1-e^{-\sigma^*_x\Delta t}}{\sigma^*_x \Delta x} \left(E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}\right)
\\
H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_y\Delta t} H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l} + \frac{1-e^{-\sigma^*_y\Delta t}}{\sigma^*_y \Delta y} \left(E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}\right)
\end{aligned}\end{split}\]
More generally, this becomes
\[\begin{split}\begin{aligned}
E_x|^{n+1}_{j+1/2,k,l} & = e^{-\sigma_y\Delta t} E_x|^{n}_{j+1/2,k,l} + \frac{1-e^{-\sigma_y\Delta t}}{\sigma_y \Delta y}\frac{c_y}{c} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}\right)
\\
E_y|^{n+1}_{j,k+1/2,l} & = e^{-\sigma_x\Delta t} E_y|^{n}_{j,k+1/2,l} - \frac{1-e^{-\sigma_x\Delta t}}{\sigma_x \Delta x}\frac{c_x}{c} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}\right)
\\
H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_x\Delta t} H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l} - \frac{1-e^{-\sigma^*_x\Delta t}}{\sigma^*_x \Delta x}\frac{c^*_x}{c} \left(E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}\right)
\\
H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_y\Delta t} H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l} + \frac{1-e^{-\sigma^*_y\Delta t}}{\sigma^*_y \Delta y}\frac{c^*_y}{c} \left(E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}\right)
\end{aligned}\end{split}\]
If we set
\[\begin{split}\begin{aligned}
c_x & = c \: e^{-\sigma_x\Delta t} \frac{\sigma_x \Delta t}{1-e^{-\sigma_x\Delta t}} \\
c_y & = c \: e^{-\sigma_y\Delta t} \frac{\sigma_y \Delta t}{1-e^{-\sigma_y\Delta t}} \\
c^*_x & = c \: e^{-\sigma^*_x\Delta t} \frac{\sigma^*_x \Delta t}{1-e^{-\sigma^*_x\Delta t}} \\
c^*_y & = c \: e^{-\sigma^*_y\Delta t} \frac{\sigma^*_y \Delta t}{1-e^{-\sigma^*_y\Delta t}}\end{aligned}\end{split}\]
then this becomes
\[\begin{split}\begin{aligned}
E_x|^{n+1}_{j+1/2,k,l} & = e^{-\sigma_y\Delta t} \left[ E_x|^{n}_{j+1/2,k,l} + \frac{\Delta t}{\Delta y} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}\right) \right]
\\
E_y|^{n+1}_{j,k+1/2,l} & = e^{-\sigma_x\Delta t} \left[ E_y|^{n}_{j,k+1/2,l} - \frac{\Delta t}{\Delta x} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}\right) \right]
\\
H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_x\Delta t} \left[ H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l} - \frac{\Delta t}{\Delta x} \left(E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}\right) \right]
\\
H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l} & = e^{-\sigma^*_y\Delta t} \left[ H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l} + \frac{\Delta t}{\Delta y} \left(E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}\right) \right]
\end{aligned}\end{split}\]
When the generalized conductivities are zero, the update equations are
\[\begin{split}\begin{aligned}
E_x|^{n+1}_{j+1/2,k,l} & = E_x|^{n}_{j+1/2,k,l} + \frac{\Delta t}{\Delta y} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j+1/2,k-1/2,l}\right)
\\
E_y|^{n+1}_{j,k+1/2,l} & = E_y|^{n}_{j,k+1/2,l} - \frac{\Delta t}{\Delta x} \left(H_z|^{n+1/2}_{j+1/2,k+1/2,l}-H_z|^{n+1/2}_{j-1/2,k+1/2,l}\right)
\\
H_{zx}|^{n+3/2}_{j+1/2,k+1/2,l} & = H_{zx}|^{n+1/2}_{j+1/2,k+1/2,l} - \frac{\Delta t}{\Delta x} \left(E_y|^{n+1}_{j+1,k+1/2,l}-E_y|^{n+1}_{j,k+1/2,l}\right)
\\
H_{zy}|^{n+3/2}_{j+1/2,k+1/2,l} & = H_{zy}|^{n+1/2}_{j+1/2,k+1/2,l} + \frac{\Delta t}{\Delta y} \left(E_x|^{n+1}_{j+1/2,k+1,l}-E_x|^{n+1}_{j+1/2,k,l}\right)
\end{aligned}\end{split}\]
as expected.
