# The perfectly matched layer (PML) for hyperbolic wave propagation problems:
A review

Kenneth Duru Corresponding author. Mathematical Sciences Institute, The
Australian National University, Australia. Gunilla Kreiss Division of
Scientific Computing, Department of Information Technology, Uppsala
University, Sweden.

###### Abstract

It is well-known that reliable and efficient domain truncation is crucial to
accurate numerical solution of most wave propagation problems. The perfectly
matched layer (PML) is a method which, when stable, can provide a domain
truncation scheme which is convergent with increasing layer width/damping. The
difficulties in using the PML are primarily associated with stability, which
can be present at the continuous level or be triggered by numerical
approximations. The mathematical and numerical analysis of the PML for
hyperbolic wave propagation problems has been an area of active research. It
is now possible to construct stable and high order accurate numerical wave
solvers by augmenting wave equations with the PML and approximating the
equations using summation-by-parts finite difference methods, continuous and
discontinuous Galerkin finite element methods. In this review we summarise the
progress made, from mathematical, numerical and practical perspectives, point
out some open problems and set the stage for future work. We also present
numerical experiments of model problems corroborating the theoretical
analysis, and numerical simulations of real-world wave propagation
demonstrating impact. Stable and parallel implementations of the PML in the
high performance computing software packages WaveQLab3D and ExaHyPE allow to
sufficiently limit the computational domain of seismological problems with
only a few grid points/elements around the computational boundaries where the
PML is active, thus saving as much as $96\%$ of the required computational
resources for a three space dimensional seismological benchmark problem.

Keyword: hyperbolic wave propagation, perfectly matched layer, Laplace
transforms, stability, discontinuous Galerkin methods, summation-by-parts
finite difference methods, penalty boundary procedures.

##  1 Introduction

##### Background.

Real-world wave propagation problems are often formulated in unbounded or very
large domains. Because of limited computational resources numerical
simulations are performed in truncated (smaller) computational domains by
introducing artificial boundaries. However, one of the main features of
propagating waves is that they can travel long distances relative to their
characteristic length-scale, the wavelength. As an example, a strong ground
shaking resulting from an earthquake or a nuclear explosion that occurred in
one continent can be recorded far away in another continent sometime after it
has occurred. This innate capability of waves to travel long distances while
carrying important information of their sources and the medium of propagation
allows us to use them to probe the world around us and develop modern multi-
media technologies. In practice, computer simulations of such wave phenomena
are restricted to smaller computational domains and areas of interest, for
example, where the effects of the strong ground motion is significant. This is
typical of numerical simulations of most real-world wave propagation problems
including earthquake seismology, seismic imaging, aero-acoustics, sonar
applications, wireless communication, ground penetrating radar technologies,
to name only a few.

For numerical simulations, it is precisely this essential feature of waves,
the radiation of waves to far field, that leads to the greatest difficulties.
In truncated computational domains outgoing wave radiations will result to
spurious reflections at artificial boundaries, which will travel back into the
simulation domain and destroy the accuracy of numerical simulations
everywhere. Therefore, in order to ensure the accuracy of numerical
simulations, artificial boundaries introduced to limit the computational
domain must be closed with reliable and accurate boundary conditions.

##### Absorbing boundary conditions and absorbing layers.

The endeavour to design efficient domain truncation strategies for waves began
over forty five years ago [52], and has yielded two powerful and equally
appealing approaches for realising effective domain truncation schemes. The
first approach is the absorbing boundary condition (ABC) [52, 4, 69, 53, 29,
85, 11, 64, 63, 56, 66, 11, 65, 14, 57, 97, 28], which is a boundary condition
enforced at an artificial boundary such that unwanted reflections from the
boundary are significantly minimised. The second approach corresponds to
extending the domain to an absorbing layer of finite thickness where the
underlying equations are transformed such that waves decay rapidly in the
layer [83, 7, 23, 15, 102, 6, 36, 41, 94, 77, 99, 59, 1, 73, 90, 72]. In order
for an absorbing layer to be effective the equations must be perfectly matched
[35, 67, 23, 15]. The perfectly matched layer (PML), [23, 15, 102, 6, 36, 41,
35, 55], are absorbing layers with the desirable properties that all waves
enter the layer and are absorbed without any reflections, regardless of
frequency and angle of incidence. The PML was initially introduced over twenty
five years ago for electromagnetic waves [23, 100, 27, 82]. However, the
approach has since then been extended to other wave propagation problems such
as elastic wave propagation, see for example [6, 49, 42, 96, 76, 35, 87, 86],
and acoustic wave propagation, for examples [59, 1, 73, 45, 44]. The
popularity of the PML is mainly due to its effective absorption properties,
versatility, simplicity, ease of derivation and implementation using standard
numerical methods. For example, using the so-called complex coordinate
stretching technique [27, 109, 95] the derivation of the PML for many linear
hyperbolic systems is simplified and takes the following three straightforward
steps. Consider the Cartesian spatial coordinates $(x,y,z)$ and a PML in
$\xi$-direction with $\xi\in\\{x,y,z\\}$. The steps for a standard PML model
are as follows:
$\bullet$ Take the Laplace transform in time of the underlying equations,
$\partial/\partial t\to s$,
$\bullet$ Transform the partial derivatives using the PML complex metric
$S_{\xi}=1+d_{\xi}(\xi)/s$,
$\partial/{\partial\xi}\to\left(1/S_{\xi}\right)\partial/{\partial\xi}$,
$\bullet$ Introduce local auxiliary variables and invert the Laplace
transforms.
Here, $s$ is the dual time variable and $d_{\xi}\geq 0$ is the absorption
function. Note that the PML metric $S_{\xi}=1+d_{\xi}(\xi)/s$ is standard, and
can be enhanced by including more parameters [27, 109, 95, 35, 55, 10] (see
also Section 3 of this review) such that the PML is more robust. The PML
transformation has some important mathematical and numerical consequences
which are critical for effective use in many practical computations. For the
past twenty years, the mathematical analysis of well-posedness and stability,
and numerical analysis of error and convergence of the PML has attracted
substantial attention in the literature, see for examples [35, 15, 102, 6, 39,
67, 16, 18, 10, 3, 2, 68, 70, 71, 46, 105, 84, 1, 12, 19, 102, 26, 13, 21, 20,
17].

##### Aim of the paper.

In this paper, we shall review the state-of-the-art results for the PML, both
when it comes to theoretical results and practical issues. We have collected
some useful results concerning how to discretely handle PML boundaries and
interfaces for several types of wave propagation problems. We will consider
PML initial boundary value problems (IBVPs) for first order hyperbolic systems
and mention some extensions to second order hyperbolic systems, on domains
partially or completely surrounded by a PML. Throughout the review we shall
highlight important results and point out some open problems.

Mathematical well-posedness and stability are fundamental features of any
useful model for accurate numerical absorption of waves. Analysis of well-
posedness and stability of the PML has been of significant interests [35, 15,
102, 6, 39, 67, 16, 18, 3, 104, 91, 2, 68, 70, 71, 13, 21, 20, 17, 21, 22]. We
will try to summarise the results without going into too much details of the
proofs. For initial value problems (IVPs) the stability of the PML can be
predicted by the so–called geometric stability condition [15] using classical
Fourier analysis. This is corroborated by the analysis in [10] by deriving
decaying energy densities in the Fourier space and in [67] using geometric
optics, and verified through several numerical experiments published in the
literature. We note that for some simple models such as the acoustic wave
equation with constant damping and constant medium parameters, analytical
solution of the PML can be derived, for instance using the Cagniard-De Hoop
method [19, 33]. However, several numerical experiments presented in the
literature suggest that, even when the corresponding IVP is stable, the PML
can suffer from exponential energy growth, in particular when the PML
interacts with certain boundary conditions. Furthermore, even when there are
no physical boundaries, in practice the PML must be implemented in a bounded
computational domain, as an IBVP. Therefore we must prove well-posedness and
stability for PML IBVPs and derive stable and convergent numerical
approximations of the PML IBVPs.

The theory and numerical methods to solve IBVPs are more elaborate and
complicated than those of a corresponding IVP; for a well-posed partial
differential equation (PDE) which is stable in the absence of boundaries can
support unstable solutions, or become ill-posed, when boundaries are
introduced [35, 60, 102, 36, 79, 42]. Indeed, the analysis of the IVP is a
necessary (first) step towards the analysis of the IBVP [60, 36, 79]. In
numerical approximations of IBVPs, most of the difficulties arise from the
boundaries. Similarly, a numerical method which is stable in the absence of
boundaries can support growth when boundary conditions are imposed [60, 62,
93, 62, 61].

We note that the mode analysis for PMLs has been extended to PML IBVPs, to
investigate the effects of boundary conditions on the well-posedness and
stability of PML IBVPs, see [35, 42, 36, 40]. The overarching result is that
the PML IBVP will not support exponentially growing solutions as long as
neither the underlying undamped IBVP nor the PML IVPs support exponentially
growing solutions. However these results are too technical to be extended to
the analysis of numerical approximations. An approach that has been
demonstrated to provide a cure for this problem is the extension of the
theoretical stability results using the energy method in the Laplace space
[36, 40, 49, 39].

One aim of this review is to demonstrate the importance and strength of the
theoretical results in the form of energy estimates. We will discuss the
results in both continuous and discrete settings, and demonstrate that the
discrete energy estimates are essential for the usefulness of PMLs. We will
focus here on summation-by-part (SBP) finite difference methods and DG finite
element methods, where there are well developed techniques to impose boundary
and interface conditions weakly with penalty terms. Such a technique allows
for mimicking continuous energy estimates.

For IBVPs, a standard approach to achieve numerical stability is to first
derive an energy estimate for the continuous problem. The next step is to
derive a discrete energy estimate, which typically is done by mimicking the
procedure in the continuous setting. If such an estimate is possible,
numerical stability is ensured. For PMLs this approach has not worked so well.
The main issue is that the PML is generally asymmetric. It becomes extremely
difficult to derive energy estimates that can be used to design stable
numerical methods for the PML in truncated domains. In many settings, a
straight forward approach would yield an exponentially growing estimate in
physical space. This would indicate well-posedness, but exponential growth is
not optimal for an absorbing model. A well-functioning PML should not support
exponential growth. We note that at constant coefficient an energy estimate
for the acoustic wave equation in second order form was recently derived for a
specific time-domain formulation of the PML [26, 12]. However, we will review
the use of energy methods in the Laplace space [36, 40, 49, 39], in general
media, which is applicable to most PML models and useful for developing stable
and high order accurate numerical methods.

Most modern numerical methods such as the discontinuous Galerkin (DG)
finite/spectral element methods and multi-block SBP finite difference methods
require interface conditions such that locally adjacent elements/blocks can be
coupled together. Interface conditions connecting elements together are often
implemented through numerical fluxes. For symmetric or symmetrisable linear
hyperlolic systems, numerical fluxes are often designed such that numerical
scheme obeys a discrete energy estimate, thus ensuring robustness of numerical
approximations. However, when the PML is present, the development of accurate
and stable numerical fluxes has proven to be a nightmare for practitioners.
Exponential and/or linear growth is often seen in numerical simulations using
modern numerical methods. As above, numerical procedures based on the energy
method in the Laplace space has proven to yield robust numerical methods and
provide a possible cure for this problem [36, 40, 49, 39]. We will summarise
the energy-based method in the Laplace [49, 39] for stable DG implementations
of the PML in acoustics and linear elastodynamics.

Another aim is to demonstrate impact on practical real-world simulations. Our
algorithms and the PML have been implemented in two different freely open
source high performance computing (HPC) software packages, WaveQLab3D [37, 38]
and ExaHyPE [49, 50], for large-scale simulation of seismic waves in
geometrically complex 3D Earth models. The software package WaveQLab3D [37,
38] is a high order accurate SBP-SAT finite difference solver. ExaHyPE is a DG
solver of arbitrary accuracy for large-scale numerical simulation of
hyperbolic wave propagation problems on dynamically adaptive curvilinear
meshes. Stable and parallel implementations of the PML in the HPC software
packages WaveQLab3D and ExaHyPE [37, 38, 39, 49, 50] allow us to sufficiently
limit the computational domain of seismological problems with only a few grid
points/elements around the computational boundaries where the PML is active.
Thus saving as much as $96\%$ of the required computational resources for a 3D
seismological benchmark problem [103]. We also present a real-world wave
propagation propagation problem which involves the simulation of 3D seismic
waves in a section of European Alpine region, with strong non-planar free-
surface topography [30, 50].

##### Structure of the paper

The remaining parts of the review will proceed as follows. In section 2 we
introduce a general model of linear hyperbolic systems in first and second
order forms. We perform dispersion relation analysis, introduce the geometric
stability condition and formulate well-posed boundary and interface
conditions. In section 3 we review the derivation of the PML using complex
coordinate stretching. Mathematical analysis of the PML at the continuous
level is reviewed in section 4. In section 5 we summarise theoretical results
on the numerical analysis of the PML discrestised with the DG finite element
method and SBP finite difference methods, respectively. In section 6 we
present some numerical experiments verifying the theoretical results and
demonstrating the importance of the continuous and discrete energy estimates.
In section 7 we present the use of PML in practical seismological applications
problems, and demonstrate impact. We draw conclusions in section 8 and
speculate on the directions for future work.

##  2 Linear hyperbolic partial differential equations

In this section, we present model equations that are representative of linear
hyperbolic wave propagation problems encountered in different application
areas, such as acoustics, seismology, electrodynamics and elastodynamics. We
will consider the equations in both first order form and second order
formulations. Plane waves will be used to discuss the wave propagation
properties of the medium at constant coefficients. We will also introduce
boundary conditions, comment on well-posedness and derived energy estimates.

###  2.1 First order systems

Consider the linear first order hyperbolic partial differential equation (PDE)
in Cartesian coordinates and a source free heterogeneous medium:

$\begin{split}\mathbf{P}^{-1}\frac{\partial{\mathbf{U}}}{\partial
t}=\sum_{\xi=x,y,z}\mathbf{A}_{\xi}\frac{\partial{\mathbf{U}}}{\partial\xi},\quad(x,y,z)\in\Omega\subset\mathbb{R}^{d},\quad
d=3,\quad t\geq 0,\end{split}$ (1)

subject to the initial condition

$\displaystyle\mathbf{U}(x,y,z,0)=\mathbf{U}_{0}(x,y,z)\in\mathrm{L}^{2}\left(\Omega\right).$

Here, $\mathbf{U}\in\mathbb{R}^{m}$, $m\geq 1$ is the unknown vector field,
the coefficient matrices are symmetric,
$\mathbf{A}_{\xi}=\mathbf{A}_{\xi}^{T}\in\mathbb{R}^{m\times m}$ and
$\mathbf{P}=\mathbf{P}^{T}\in\mathbb{R}^{m\times m}$, with
$\mathbf{U}^{T}\mathbf{P}\mathbf{U}>0$. In general the matrix $\mathbf{P}$
depends on the spatial coordinates $x,y,z$, and encodes the material
parameters of the underlying medium. The constant coefficients and non-
dimensional matrices $\mathbf{A}_{\xi}$ encapsulate the underlying linear
conservation law and the corresponding linear constitutive relation.

As will be shown below, depending on the coefficient matrices
$\mathbf{P},\mathbf{A}_{\xi}$ the system (1) describes first order linear
hyperbolic wave propagation problems encountered in different application
areas such acoustics, elastodynamics, electromagnetics, etc. For example in an
acoustic medium we denote $p$ the pressure and
$\mathbf{v}=(v_{x},v_{y},v_{z})^{T}$ the particle velocity, $\rho$ the
density, and $\kappa>0$ the bulk modulus of the medium, we have

$\displaystyle\mathbf{U}=\begin{bmatrix}p\\\
\mathbf{v}\end{bmatrix},\quad\mathbf{P}=\begin{pmatrix}\kappa&\mathbf{0}^{T}\\\
\mathbf{0}&\rho^{-1}\mathbf{I}_{3}\end{pmatrix},\quad\mathbf{A}_{\xi}=\begin{pmatrix}{0}&-\mathbf{e}_{\xi}^{T}\\\
-\mathbf{e}_{\xi}&\mathbf{0}_{3}\end{pmatrix},\quad\mathbf{0}=\left(0,0,0\right)^{T},\quad\mathbf{e}_{\xi}=\left({e}_{\xi
x},{e}_{\xi y},{e}_{\xi z}\right)^{T},$

where $\mathbf{I}_{3}$ is a 3-by-3 identity matrix, $\mathbf{0}_{3}$ is a
3-by-3 zero matrix, and

$\begin{split}{e}_{\xi\eta}=\left\\{\begin{array}[]{rl}0\quad{}\quad\text{if}\quad\xi\neq\eta,\\\
1\quad{}\quad\text{if}\quad\xi=\eta.\end{array}\right.\end{split}$

To describe wave propagation in elastic solids, we introduce the unknown wave
fields

$\displaystyle\mathbf{U}\left(x,y,z,t\right)=\begin{bmatrix}\mathbf{v}(x,y,z,t)\\\
\bm{\sigma}(x,y,z,t)\end{bmatrix},$ (2)

with the particle velocity vector is
$\mathbf{v}=\left[v_{x},v_{y},v_{z}\right]^{T}$, and the stress vector is
given by
$\bm{\sigma}=\left[\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy},\sigma_{xz},\sigma_{yz}\right]^{T}$.
The symmetric constant coefficient matrices $\mathbf{A}_{\xi}$ describing the
conservation of momentum and the constitutive relation, defined by Hooke’s
law, are given by

$\displaystyle\mathbf{A}_{\xi}=\begin{pmatrix}\mathbf{0}_{3}&\mathbf{a}_{\xi}\\\
\mathbf{a}_{\xi}^{T}&\mathbf{0}_{6}\end{pmatrix},\quad\mathbf{a}_{x}=\begin{pmatrix}1&0&0&0&0&0\\\
0&0&0&1&0&0\\\ 0&0&0&0&1&0\end{pmatrix},$
$\displaystyle\mathbf{a}_{y}=\begin{pmatrix}0&0&0&1&0&0\\\ 0&1&0&0&0&0\\\
0&0&0&0&0&1\end{pmatrix},\quad\mathbf{a}_{z}=\begin{pmatrix}0&0&0&0&1&0\\\
0&0&0&0&0&1\\\ 0&0&1&0&0&0\\\ \end{pmatrix},$

where $\mathbf{0}_{3}$ and $\mathbf{0}_{6}$ are the $3$-by-$3$ and $6$-by-$6$
zero matrices. The symmetric positive definite material parameter matrix
$\mathbf{P}$ is defined by

$\displaystyle\mathbf{P}=\begin{pmatrix}\rho^{-1}\mathbf{I}_{3}&\mathbf{0}\\\
\mathbf{0}^{T}&\mathbf{C}\end{pmatrix},\quad\mathbf{0}=\begin{pmatrix}0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&0&0&0\end{pmatrix},$ (3)

where $\rho(x,y,z)>0$ is the density of the medium, and
$\mathbf{C}=\mathbf{C}^{T}>0$ is the symmetric positive definite matrix of
elastic constants.

For real functions, we introduce the weighted $L^{2}$ inner product and the
corresponding energy norm

$\displaystyle\left(\mathbf{U},\mathbf{F}\right)_{P}=\int_{\Omega}{\frac{1}{2}[\mathbf{U}^{T}\mathbf{P}^{-1}\mathbf{F}]dxdydz},\quad\|\mathbf{U}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}=\left(\mathbf{U},\mathbf{U}\right)_{P}$
(4)

The weighted $L^{2}$-norm
$\|\mathbf{U}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}$ is the physical
energy in the medium.

For the IVP (1), $\Omega=\mathbb{R}^{3}$, the energy (conservation) equation
follows

$\displaystyle\frac{d}{dt}\|\mathbf{U}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}=0.$
(5)

The conservation of energy (5) indicates that the IVP (1) is well-posed and
strongly stable.

###  2.2 Second order systems

Now consider the linear second order hyperbolic PDE in Cartesian coordinates
in a source free heterogeneous medium

$\begin{split}\bm{\Phi}^{-1}\frac{\partial^{2}{\mathbf{u}}}{\partial
t^{2}}=\sum_{\xi=x,y,z}\frac{\partial}{\partial\xi}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right),\quad(x,y,z)\in\Omega\subset\mathbb{R}^{3},\quad
t\geq 0,\end{split}$ (6)

subject to the initial condition

$\displaystyle\mathbf{u}(x,y,z,0)=\mathbf{u}_{0}(x,y,z)\in\mathrm{H}^{1}\left(\Omega\right),\quad\frac{\partial}{\partial
t}\mathbf{u}(x,y,z,0)=\mathbf{v}_{0}(x,y,z)\in\mathrm{L}^{2}\left(\Omega\right).$

Here, $\mathbf{u}(x,y,z,t)\in\mathbb{R}^{n}$ is the unknown vector field,
$\bm{\Phi}\in\mathbb{R}^{n\times n}$ is a symmetric and positive definite
matrix, and the matrices $\mathbf{A}_{\xi\eta}\in\mathbb{R}^{n\times n}$ are
such that the potential energy matrix $\mathscr{P}\in\mathbb{R}^{dn\times dn}$
defined below is symmetric semi-positive definite, that is

$\mathscr{P}=\begin{pmatrix}\mathbf{A}_{xx}&\mathbf{A}_{xy}&\mathbf{A}_{xz}\\\
\mathbf{A}_{xy}^{T}&\mathbf{A}_{yy}&\mathbf{A}_{yz}\\\
\mathbf{A}_{xz}^{T}&\mathbf{A}_{yz}^{T}&\mathbf{A}_{zz}\end{pmatrix},\quad\mathbf{U}^{T}\mathscr{P}\mathbf{U}\geq
0,\quad\forall\mathbf{U}\in\mathbb{R}^{dn},\quad d=3.$

Many linear hyperbolic PDEs originally arise as second order form (6) but can
be reduced to a first order system (1), by introducing new variables, see eg.
the classical acoustic wave equation and the elastic wave equation. Some
models also initially appear as first order systems (1) but they can be
rewritten as second order systems (6), eg linearised Euler equations of fluid
dynamics, Maxwell’s equation, and linearised MHD equations, to name a few.

For example, consider an acoustic medium where $p$ denote the pressure,
$\rho,\kappa>0$ denote the density and the bulk modulus of the medium with the
acoustic waves peed $c=\sqrt{\kappa/\rho}$, we have

$\displaystyle\mathbf{u}=p,\quad\bm{\Phi}={\kappa},\quad\mathbf{A}_{\xi\eta}=\left\\{\begin{array}[]{rl}0\quad{}\quad\text{if}\quad\xi\neq\eta,\\\
\frac{1}{\rho}\quad{}\quad\text{if}\quad\xi=\eta.\end{array}\right.$

Define the energy in the medium

$\displaystyle
E(t)=\frac{1}{2}\int_{\Omega}{\left(\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\bm{\Phi}^{-1}\frac{\partial{\mathbf{u}}}{\partial
t}+\begin{pmatrix}\frac{\partial{\mathbf{u}}}{\partial x}\\\
\frac{\partial{\mathbf{u}}}{\partial y}\\\
\frac{\partial{\mathbf{u}}}{\partial
z}\end{pmatrix}^{T}\mathscr{P}\begin{pmatrix}\frac{\partial{\mathbf{u}}}{\partial
x}\\\ \frac{\partial{\mathbf{u}}}{\partial y}\\\
\frac{\partial{\mathbf{u}}}{\partial z}\end{pmatrix}\right)}dxdydz>0.$ (7)

The energy $E(t)$ given in (7) defines a semi-norm. As above, for Cauchy
problems, $\Omega=\mathbb{R}^{3}$, the energy equation follows

$\displaystyle\frac{d}{dt}E(t)=0.$ (8)

Thus the energy is conserved, $E(t)=E(0)$, for all $t\geq 0$. Again, this
analysis indicates that the Cauchy problem for (6) is well-posed and
asymptotically stable.

###  2.3 Dispersion relation

To understand the wave propagation properties of the models (1) and (6) it is
useful to consider wave-like solutions

$\displaystyle\mathbf{U}\left(x,y,z,t\right)=\mathbf{U}_{0}e^{st+i\left(k_{x}x+k_{y}y+k_{z}z\right)},\quad\left(k_{x},k_{y},k_{z}\right)\in\mathbb{R}^{3}.$
(9)

In (9), $\mathbf{k}=\left(k_{x},k_{y},k_{z}\right)$ is the wave vector, and
$\mathbf{U}_{0}$ is a vector of constant amplitude called the polarization
vector. By inserting (9) in (1) we have the solvability condition, often
called the dispersion relation

$\displaystyle
F\left(s,\mathbf{k}\right):=\Big{|}{P}^{-1}s-\sum_{\xi=x,y,z}{ik_{\xi}}\mathbf{A}_{\xi}\Big{|}=0.$
(10)

The indeterminate $s$ is related to the temporal frequency and will be
determined by the dispersion relation (10). Since (1) conserve energy the
roots $s$ must be zero or purely imaginary, that is $s\in\mathbb{C}$ with
$\Re{s}=0$. Otherwise if the roots $s$ have non-zero real parts the energy
will grow or decay, contradicting the energy equations (5).

Similarly, inserting the simple wave solution (9) in the second order system
(6), we have

$\displaystyle
G\left(s,\mathbf{k}\right):=\Big{|}\bm{\Phi}^{-1}s^{2}+\sum_{\xi=x,y,z}\sum_{\eta=x,y,z}{k_{\xi}}k_{\eta}\mathbf{A}_{\xi,\eta}\Big{|}=0.$

If (1) and (6) model the same physical phenomena, for some natural number
$r\in\mathbb{N}$, we have

$\displaystyle
F\left(s,\mathbf{k}\right)=s^{r}G\left(s,\mathbf{k}\right),\quad{r}\geq 0,$

where

$\displaystyle s=\\{s_{0}\in\mathbb{C}:\quad
G\left(s_{0},\mathbf{k}\right)=0,\quad|s_{0}|\neq 0\\}.$

Note that $F\left(s,\mathbf{k}\right)$ will have $r$ zero roots $s=0$ and
$m-r$ nonzero roots $s=\\{s_{0}\\}$. The nonzero roots $s$ of
$G\left(s,\mathbf{k}\right)$ correspond to propagating physical modes.
Specifically in acoustics with the wave speed $c=\sqrt{\kappa/\rho}$ we have

$G\left(s,\mathbf{k}\right)=\frac{s^{2}}{c^{2}}+|\mathbf{k}|^{2},\quad
F\left(s,\mathbf{k}\right)=s^{2}G\left(s,\mathbf{k}\right),$

where $|\mathbf{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$ and the nonzero
roots are $s=\pm ic|\mathbf{k}|$. In isotropic linear elastic medium we have

$G\left(s,\mathbf{k}\right)=\left(\frac{s^{2}}{c_{s}^{2}}+|\mathbf{k}|^{2}\right)\left(\frac{s^{2}}{c_{s}^{2}}+|\mathbf{k}|^{2}\right)\left(\frac{s^{2}}{c_{p}^{2}}+|\mathbf{k}|^{2}\right),\quad
F\left(s,\mathbf{k}\right)=s^{3}G\left(s,\mathbf{k}\right),$

where $c_{s}$, $c_{p}$ are the S-wave speed and P-wave speed, respectively, of
the medium, and the nonzero roots are are given by the simple roots $s=\pm
ic_{s}|\mathbf{k}|$, $s=\pm ic_{s}|\mathbf{k}|$, and $s=\pm
ic_{p}|\mathbf{k}|$. In anisotropic media the dispersion relation are
complicated nonlinear algebraic expressions, and the roots are difficult to
determine explicitly. However, as above, one can show that the roots are
purely imaginary, that is $s\in\mathbb{C}$ with $\Re{s}=0$.

We write $s=i\omega$ where $\omega\in\mathbb{R}$ is called the temporal
frequency, and introduce

$\displaystyle\mathbf{K}=\left(\frac{k_{x}}{|\mathbf{k}|},\frac{k_{y}}{|\mathbf{k}|},\frac{k_{z}}{|\mathbf{k}|}\right),\quad\text{normalized
propagation direction},$ (11)
$\displaystyle\mathbf{V}_{p}=\left(\frac{\omega}{k_{x}},\frac{\omega}{k_{y}},\frac{\omega}{k_{z}}\right),\quad\text{phase
velocity},$ (12)
$\displaystyle\mathbf{S}=\left(\frac{k_{x}}{\omega},\frac{k_{y}}{\omega},\frac{k_{z}}{\omega}\right),\quad\text{slowness
vector},$ (13)
$\displaystyle\mathbf{V}_{g}=\left(\frac{\partial\omega}{\partial
k_{x}},\frac{\partial\omega}{\partial k_{y}},\frac{\partial\omega}{\partial
k_{z}}\right),\quad\text{group velocity}.$ (14)

For a constant coefficient medium with no boundaries, or with periodic
boundary conditions, the dispersion relation
$F\left(i\omega,k_{x},k_{y},k_{z}\right)=0$ and the quantities $\mathbf{K}$,
$\mathbf{V}_{p}$, $\mathbf{S}$, $\mathbf{V}_{g}$, defined above, provide a
rather comprehensive description of the wave propagation properties of the
medium. When physical boundaries are introduced the additional presence of
boundary wave modes, such as Rayleigh waves, introduces more complex and
interesting wave phenomena. More general media, such as complex elastic media,
can also support guided waves, such as Love waves, Stoneley waves, Scholte
waves, conical waves, and other interface waves.

###### Remark 1

We note however that the analysis of the PML in the presence of interface
waves is still an open problem and will not be considered in this review.

Since (1) contains only first derivative terms and no lower order terms then
$G\left(s,\mathbf{k}\right)$ is homogeneous in $\omega,k_{x},k_{y},k_{z}$, and
we can rewrite $G\left(s,\mathbf{k}\right)=0$ as
$G\left(i,\mathbf{S}\right)=0$.

######  Definition 1 (Geometric stability condition)

At any point on the slowness surface $\mathbf{S}=(S_{x},S_{y},S_{z})$ with
$G\left(i,\mathbf{S}\right)=0$ , the geometric stability condition in the
$\xi$-direction (for $\xi=x,y,z$) is defined by

${V}_{p\xi}{V}_{g\xi}\geq 0.$

As was shown in [35, 15, 6] the Geometric stability condition is necessary for
the stability of the PML IVP. This will be discussed in more detail in Section
4 where the mathematical analysis of the PML is reviewed.

###  2.4 Boundary conditions

Effective treatments of boundary conditions are critical for the stability and
efficacy of the PML [35, 42, 36, 40, 49, 39]. In this section, we formulate
well-posed and energy stable boundary conditions for the first order system
(1) and the second order system (6). The boundary conditions will be extended
when the PML is introduced.

Consider the 3D cuboidal domain

$\Omega=\\{(x,y,z):-1\leq x\leq 1,\quad-1\leq y\leq 1,\quad-1\leq z\leq 1\\}.$
(15)

Stable and well-posed boundary conditions are needed to close the rectangular
surfaces of the boundaries of the cuboidal domain. Define the reference
boundary surface

$\displaystyle\widetilde{\Gamma}=[-1,1]\times[-1,1].$ (16)

####  2.4.1 Boundary condition for first order systems

The application of the energy method to (1) yields

$\displaystyle\frac{d}{dt}\|\mathbf{U}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}=\mathrm{BT}\left(\mathbf{U}\right),$
(17)

where the boundary term is given by

$\displaystyle\mathrm{BT}\left(\mathbf{U}\right)$
$\displaystyle\equiv\frac{1}{2}\oint_{\partial\Omega}\left[\mathbf{U}^{T}\left(\sum_{\xi=x,y,z}n_{\xi}\mathbf{A}_{\xi}\right)\mathbf{U}\right]dS$
$\displaystyle=\sum_{\xi=x,y,z}\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\Big{|}_{\xi=1}-\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\Big{|}_{\xi=-1}\right]\frac{dxdydz}{d\xi}.$
(18)

For the acoustic wave equation we have

$\frac{1}{2}\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}=-v_{\xi}p.$

For linear elastodynamics, at the boundaries $\xi=\pm 1$ we introduce with the
tractions and particle velocity on the boundary,

$\mathbf{T}:=\mathbf{a}_{\xi}\bm{\sigma}=(T_{x},T_{y},T_{z})^{T},\quad\mathbf{v}=(v_{x},v_{y},v_{z})^{T}.$

We have

$\frac{1}{2}\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}=\mathbf{v}^{T}\mathbf{T}=\sum_{\xi=x,y,z}v_{\xi}T_{\xi}.$

See [50, 49, 39] for details.

Well-posed boundary conditions are imposed such that the boundary term
$\mathrm{BT}\left(\mathbf{U}\right)$ defined in (2.4.1) is negative semi-
definite. We consider linear boundary operators

$\displaystyle\mathbf{L}\left(\mathbf{U}\right)=0,\quad\xi=\pm 1,$ (19)

such that

$\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\geq 0,\quad\text{at
$\xi=-1$},\quad\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\leq 0,\quad\text{at
$\xi=1$}.$

Such boundary conditions are formulated by prescribing the amplitude of the
incoming characteristic variables, see [60, 43, 50, 49] for more details. Thus
the number of boundary conditions must be the same as the number of ingoing
characteristics at the boundary.

Specifically for acoustics we introduce the impedance $Z=\rho c$, the
reflection coefficients $r_{\xi}$, with $|r_{\xi}|\leq 1$ and the well-posed
boundary condition

$\mathbf{L}\left(\mathbf{U}\right):=\frac{1-r_{\xi}}{2}Zv_{\xi}\mp\frac{1+r_{\xi}}{2}p=0,\quad|r_{\xi}|\leq
1,\quad\text{at}\quad\xi=\pm 1.$ (20)

For linear elasticity, on a given boundary $\xi=\pm 1$, we introduce the
impedance $Z_{\eta}>0$ (for $\eta=x,y,z$) and the well-posed boundary
conditions

$[\mathbf{L}\left(\mathbf{U}\right)]_{\eta}:=\frac{1-r_{\xi}}{2}Z_{\eta}v_{\eta}\pm\frac{1+r_{\xi}}{2}T_{\eta}=0,\quad|r_{\xi}|\leq
1,\quad\text{at}\quad\xi=\pm 1.$ (21)

The reflection coefficient $r_{\xi}=-1$ will correspond to boundary conditions
for velocities, that is, soft-wall boundary conditions ($v_{\xi}=0$) for
acoustics and clamped-wall boundary conditions ($v_{\eta}=0$) for
elastodynamics. With $r_{\xi}=1$ we will have boundary conditions for pressure
and traction, that is, hard-wall boundary conditions ($p=0$) for acoustics and
free-surface boundary conditions ($\mathbf{T}=0$) for elastodynamics. The zero
reflection coefficients, $r_{\xi}=0$, yield the classical first order ABC [52,
85, 11, 64, 63, 56, 66, 11, 65, 28].

The boundary conditions (20) and (21) with $|r_{\xi}|\leq 1$ satisfy
$\mathrm{BT}\left(\mathbf{U}\right)\leq 0$, and ensure bounded solutions
through

$\displaystyle\frac{d}{dt}\|\mathbf{U}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}=\mathrm{BT}\left(\mathbf{U}\right)\leq
0.$ (22)

This energy loss through the boundary is what a stable numerical method should
emulate.

####  2.4.2 Boundary condition for second order systems

Consider the second order system in the bounded domain $(x,y,z)\in[-1,1]^{3}$,
the energy method gives

$\displaystyle{\frac{d}{dt}{E}(t)=\mathrm{BT}\left(\mathbf{u}\right)},$

where the boundary term is given by

$\begin{split}&\mathrm{BT}\left(\mathbf{u}\right)\equiv\frac{1}{2}\oint_{\partial\Omega}\left[\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right)\right]dS\\\
&=\sum_{\xi=x,y,z}\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right)\Big{|}_{\xi=1}-\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right)\Big{|}_{\xi=-1}\right]\frac{dxdydz}{d\xi}.\end{split}$
(23)

We consider linear boundary operators

$\displaystyle\mathbf{\mathscr{B}}\mathbf{u}=0,$ (24)

such that the boundary term is never positive
$\mathrm{BT}\left(\mathbf{u}\right)\leq 0$. As above, we will consider
specifically the well-posed linear operator

$\displaystyle\mathbf{\mathscr{B}}=\frac{1-r_{\xi}}{2}\mathbf{Z}\frac{\partial{}}{\partial
t}\pm\frac{1+r_{\xi}}{2}\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{}}{\partial\eta},\quad\mathbf{Z}=\mathbf{Z}^{T}\geq
0,\quad\xi=\pm 1,$ (25)

with reflection coefficients $|r_{\xi}|\leq 1$. For the acoustic wave equation
we have $\mathbf{Z}=\rho c>0$, and for linear elasticity the impedance is a
$3$-by-$3$ diagonal matrix with diagonal entries

$\begin{split}&Z_{\eta,\eta}=\left\\{\begin{array}[]{rl}Z_{p}=\rho
c_{p}>0,{}{}&\text{if}\quad\eta=\xi,\\\ Z_{s}=\rho
c_{s}>0,{}{}&\text{if}\quad\eta\neq\xi.\end{array}\right.\end{split}$

Again, the reflection coefficient $r_{\xi}=-1$ will correspond to Dirichlet
boundary conditions

$\displaystyle\mathbf{u}=0\implies\frac{\partial{\mathbf{u}}}{\partial
t}=0,\quad\xi=\pm 1,$ (26)

equivalent to the soft-wall boundary conditions for acoustics and clamped-wall
boundary conditions for elastodynamics. With $r_{\xi}=1$ we have the Neumann
boundary condition

$\displaystyle\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}=0,\quad\xi=\pm
1,$ (27)

equivalent to hard-wall boundary conditions for acoustics and free-surface
boundary conditions for elastodynamics. The zero reflection coefficients,
$r_{\xi}=0$, yield the classical first order ABC,

$\displaystyle\frac{1}{2}\mathbf{Z}\frac{\partial{\mathbf{u}}}{\partial
t}\pm\frac{1}{2}\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}=0,\quad\xi=\pm
1.$ (28)

Note that for $|r_{\xi}|=1$ we have

$\mathrm{BT}\left(\mathbf{u}\right)=0,$

and if $|r_{\xi}|<1$

$\mathrm{BT}\left(\mathbf{u}\right)=-\sum_{\xi=x,y,z}\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\frac{1-r_{\xi}}{1+r_{\xi}}\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\mathbf{Z}\frac{\partial{\mathbf{u}}}{\partial
t}\Big{|}_{\xi=1}+\frac{1-r_{\xi}}{1+r_{\xi}}\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\mathbf{Z}\frac{\partial{\mathbf{u}}}{\partial
t}\Big{|}_{\xi=-1}\right]\frac{dxdydz}{d\xi}\leq 0.$

We have

$\displaystyle\frac{d}{dt}{E}(t)=\mathrm{BT}\left(\mathbf{u}\right)\leq 0.$
(29)

###  2.5 Interface conditions

Some numerical methods such as the DG and multi-block SBP finite difference
methods require interface conditions such that locally adjacent
elements/blocks can be coupled together. This is also the situation in many
settings where material parameters are discontinuous. We introduce physical
interface conditions which we will use to couple local elements to the global
domain. Later we will summarise how the interface conditions will be perturbed
by the PML.

Consider the IVP with a planar interface at $\xi=0$ and denote the
corresponding fields and material parameters in the positive/negative sides of
the interface with the superscripts $+/-$. We define the jumps in scalar or
vector valued fields by $[\\![{{a}]\\!]}=a^{+}-a^{-}$.

####  2.5.1 Interface condition for the first order system

As before, using the energy method we have

$\displaystyle\frac{d}{dt}\left(\|\mathbf{U}^{-}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}+\|\mathbf{U}^{+}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}\right)=\mathrm{IT}\left(\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\right).$
(30)

Here, the interface term is defined by

$\displaystyle\mathrm{IT}\left(\mathbf{a}\right)=-\frac{1}{2}\int_{\widetilde{\Gamma}}[\\![\mathbf{a}]\\!]\frac{dxdydz}{d\xi}.$
(31)

At interfaces the physics determines the interface conditions. In a standard
model with conservation of the physical energy, interface conditions

$\displaystyle\mathbb{I}\mathbb{T}\left(\mathbf{A}_{\xi},\mathbf{U}\right)=0,$
(32)

must guarantee that the interface term vanishes identically
$\mathrm{IT}\left(\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\right)=0$. For
well-posedness, the number of interface conditions must correspond to the
number of in-going characteristics at the interface.

In acoustics we prescribe the continuity of the pressure and the normal
velocity, having

$\mathbb{I}\mathbb{T}\left(\mathbf{A}_{\xi},\mathbf{U}\right):=\\{[\\![{{p}}]\\!]=0,\quad[\\![{{v}}_{\xi}]\\!]=0,\quad\text{at}\quad\xi=0\\}.$
(33)

In linear elastodynamics, we impose force balance, no slip and no opening
conditions

$\mathbb{I}\mathbb{T}\left(\mathbf{A}_{\xi},\mathbf{U}\right):=\\{[\\![{\mathbf{T}}]\\!]=0,\quad[\\![{\mathbf{v}}]\\!]=0,\quad\text{at}\quad\xi=0\\}.$
(34)

Note that the interface conditions (33)–(34) ensure
$\mathrm{IT}\left(\mathbf{U}^{T}\mathbf{A}_{\xi}\mathbf{U}\right)=0$ and
conserve energy

$\displaystyle\frac{d}{dt}\left(\|\mathbf{U}^{-}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}+\|\mathbf{U}^{+}\left(\cdot,\cdot,\cdot,t\right)\|_{P}^{2}\right)=0.$

An accurate and stable numerical method should as much as possible mimic this
energy estimate.

####  2.5.2 Interface condition for the second order system

As before, using the energy method we have

$\displaystyle\frac{d}{dt}\left({E}^{-}(t)+{E}^{+}(t)\right)=\mathrm{IT}\left(\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right)\right),$
(35)

where ${E}^{-}(t),{E}^{+}(t)$ are the energy in the negative and positive sub-
domains, and the interface term $\mathrm{IT}$ is the surface integral defined
by (31). Here the the well-posed interface condition is

$[\\![\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}]\\!]=0,\quad[\\![\frac{\partial{\mathbf{u}}}{\partial
t}]\\!]=0.$ (36)

The interface condition (36) imposes the continuity of the flux and the time-
derivatives of the wave fields, and ensures that the interface term
$\mathrm{IT}$ vanishes

$\displaystyle\mathrm{IT}\left(\frac{\partial{\mathbf{u}}}{\partial
t}^{T}\left(\sum_{\eta=x,y,z}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\right)\right)=0.$

We again have

$\displaystyle\frac{d}{dt}\left({E}^{-}(t)+{E}^{+}(t)\right)=0,$

and the total energy is conserved.

##  3 Perfectly matched layers

There are two standard approaches, the complex coordinate stretching technique
[27, 109, 95] and Bérenger’s splitting method [23, 15]. The two approaches are
mathematically analogous but yield two different formulations of the PML, the
unsplit PML and the split-field PML, respectively. The difference lies in the
choice of variables in the time domain. As shown in [36], using standard PML
metric, by appropriate choice of variables the split-field PML [23] can be
formulated as the unsplit PML [3, 36, 46]. This shows that, for linear
problems, the general solutions of the two PML formulations are identical.
However, the complex coordinate stretching technique [27] simplifies the
construction of the PML for several hyperbolic PDEs [36, 40, 6]

Here, we will first demonstrate the construction of the PML for the first
order system (1) by using the complex coordinate stretching technique [27].
The techniques extends to linear second order hyperbolic systems in a
straightforward manner [42, 41, 35, 108, 92, 33, 74, 76].

Let the Laplace transform, in time, of $\mathbf{U}\left(x,y,z,t\right)$ be
defined by

$\widetilde{\mathbf{U}}(x,y,z,s)=\int_{0}^{\infty}e^{-st}{\mathbf{U}}\left(x,y,z,t\right)\text{dt},\quad
s=a+ib,\quad\Re{s}=a>0.$ (37)

We consider a setup where the PML is included in all spatial coordinates. Take
the Laplace transform, in time, of equation (1). The PML can be constructed in
each coordinate, $\xi=x,y,z$, using $\partial/\partial{\xi}\to
1/S_{\xi}\partial/\partial{\xi}$. Here

$\displaystyle\widetilde{\xi}=\int_{0}^{\xi}S_{\eta}d\eta,\quad
S_{\xi}=\gamma_{\xi}\left(\xi\right)\left(1+\frac{d_{\xi}\left(\xi\right)}{s+\alpha_{\xi}(\xi)}\right),$
(38)

are the complex cordinates and the PML metric, with $s$ denoting the Laplace
dual time variable. The PML damping functions $d_{\xi}(\xi)\geq 0$ are zero in
the interior, where we solve the PDE (1), and take positive values in the PML.
The grid stretching/compression parameter $\gamma_{\xi}>0$ [41, 35] are unity
$\gamma_{\xi}\equiv 1$ in the interior where $d_{\xi}=0$ but
increases/decreases inside the PML. The nonnegative function $\alpha_{\xi}\geq
0$ is called the complex frequency shift (CFS) [82]. The standard PML metric
parameters will correspond to $\gamma_{\xi}=1$ and $\alpha_{\xi}=0$. It is
possible to use even a more general and complicated PML metric with more
parameters, see for example [10].

###  3.1 First order systems

Now, take the Laplace transform in time of (1), (19) and (32), and introduce
the complex change of variable $\partial/\partial{\xi}\to
1/S_{\xi}\partial/\partial{\xi}$. We have the PML in the Laplace space

$\begin{split}\mathbf{P}^{-1}s\widetilde{\mathbf{U}}(x,y,z,s)&=\sum_{\xi=x,y,z}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}(x,y,z,s)}}{\partial\widetilde{\xi}}=\sum_{\xi=x,y,z}\mathbf{A}_{\xi}\frac{1}{S_{\xi}}\frac{\partial{\widetilde{\mathbf{U}}(x,y,z,s)}}{\partial\xi},\end{split}$
(39)

subject to

$\displaystyle\mathbf{L}\left(\widetilde{\mathbf{U}}\right)=0,\quad\xi=\pm 1,$
(40)

and

$\displaystyle\mathbb{I}\mathbb{T}\left(\mathbf{A}_{\xi},\widetilde{\mathbf{U}}\right)=0.$
(41)

Note that the for the first order system the boundary condition (40) and the
interface condition (41) are not modified by the PML. Also note the identity

$\displaystyle\frac{1}{S_{\xi}}=\frac{1}{\gamma_{\xi}}-\frac{1}{S_{\xi}}\frac{d_{\xi}}{s+\alpha_{\xi}}.$
(42)

In order to localize the PML in time, we introduce the auxiliary variable

$\begin{split}\widetilde{\mathbf{w}}_{\xi}=\frac{1}{\left(s+\alpha_{\xi}\right)S_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi},\end{split}$
(43)

having

$\mathbf{P}^{-1}s\widetilde{\mathbf{U}}=\sum_{\xi=x,y,z}\mathbf{A}_{\xi}\frac{1}{S_{\xi}}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}=\sum_{\xi=x,y,z}\left[\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}-d_{\xi}\widetilde{\mathbf{w}}_{\xi}\right].$
(44)

We then invert the Laplace transforms in (44) and (43), and we have the time-
dependent PML

$\begin{split}\mathbf{P}^{-1}\frac{\partial{\mathbf{U}}}{\partial
t}=\sum_{\xi=x,y,z}\left[\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\mathbf{U}}}{\partial\xi}-d_{\xi}\mathbf{w}_{\xi}\right],\quad(x,y,z)\in\Omega,\quad
t\geq 0,\end{split}$ (45)
$\begin{split}\frac{\partial{\mathbf{w}_{\xi}}}{\partial
t}=\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\mathbf{U}}}{\partial\xi}-\left(\alpha_{\xi}+d_{\xi}\right)\mathbf{w}_{\xi},\end{split}$
(46)

with the boundary condition (19) and the interface conditions (32). In the
interior, where the PML is not activated, we have $\gamma_{\xi}=1,d_{\xi}=0,$
which recovers the original equation (1).

###### Remark 2

Note that for the first order systems the boundary/interface conditions, (19)
and (32), are not modified by the PML. However, the numerical treatments of
the boundary/interface conditions will be critical for numerical stability
when the PML is active.

###  3.2 Second order systems

For second order systems the ideas are similar. As above, take the Laplace
transform of (6) in time, and introduce the complex change of variable
$\partial/\partial{\xi}\to 1/S_{\xi}\partial/\partial{\xi}$,
$\partial/\partial{\eta}\to 1/S_{\eta}\partial/\partial{\eta}$. We have the
PML in the Laplace space

$\begin{split}\bm{\Phi}^{-1}s^{2}\widetilde{\mathbf{u}}&=\sum_{\widetilde{\xi}=\widetilde{x},\widetilde{y},\widetilde{z}}\frac{\partial}{\partial\widetilde{\xi}}\sum_{\widetilde{\eta}=\widetilde{x},\widetilde{y},\widetilde{z}}\mathbf{A}_{\widetilde{\xi}\widetilde{\eta}}\frac{\partial{\widetilde{\mathbf{u}}}}{\partial\widetilde{\eta}},\\\
&=\sum_{\xi=x,y,z}\frac{1}{S_{\xi}}\frac{\partial}{\partial\xi}\sum_{\eta=x,y,z}\frac{1}{S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{\widetilde{\mathbf{u}}}}{\partial\eta},\end{split}$
(47)

subject to the boundary condition

$\begin{split}\widetilde{\mathbf{\mathscr{B}}}\widetilde{\mathbf{u}}=0,\quad\widetilde{\mathbf{\mathscr{B}}}=\frac{1-r_{\xi}}{2}\mathbf{Z}\mathbf{I}s\pm\frac{1+r_{\xi}}{2}\sum_{\eta=x,y,z}\frac{1}{S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{}}{\partial\eta},\quad\mathbf{Z}=\mathbf{Z}^{T}\geq
0,\quad\xi=\pm 1.\end{split}$ (48)

and the interface conditions,

$[\\![\sum_{\eta=x,y,z}\frac{1}{S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{\widetilde{\mathbf{u}}}}{\partial\eta}]\\!]=0,\quad[\\![s{\widetilde{\mathbf{u}}}]\\!]=0,$
(49)

which ensure that the interface term vanishes

$\mathrm{IT}\left(\left(s{\widetilde{\mathbf{u}}}\right)^{\dagger}\left(\sum_{\eta=x,y,z}\frac{1}{S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{\widetilde{\mathbf{u}}}}{\partial\eta}\right)\right)=0.$

For second order systems, if the underlying boundary condition involves
spatial derivatives the PML will transform the boundary conditions.

Localising the PML in time for the second order system may involve significant
algebraic manipulations. To do this, first, we multiply (47) with
$S_{x}S_{y}S_{z}$, and (48) and (49) with $S_{x}S_{y}S_{z}/S_{\xi}$. We have

$\begin{split}\bm{\Phi}^{-1}s^{2}S_{x}S_{y}S_{z}\widetilde{\mathbf{u}}(x,y,z,s)&=\sum_{\xi=x,y,z}\frac{\partial}{\partial\xi}\sum_{\eta=x,y,z}\frac{S_{x}S_{y}S_{z}}{S_{\xi}S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{\mathbf{u}}}{\partial\eta}\end{split}$

with

$\widetilde{\mathbf{\mathscr{B}}}\widetilde{\mathbf{u}}=0,\quad\widetilde{\mathbf{\mathscr{B}}}=\frac{1-r_{\xi}}{2}\mathbf{Z}\mathbf{I}s\frac{S_{x}S_{y}S_{z}}{S_{\xi}}\pm\frac{1+r_{\xi}}{2}\sum_{\eta=x,y,z}\frac{S_{x}S_{y}S_{z}}{S_{\xi}S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{}}{\partial\eta},\quad\xi=\pm
1,$

and

$[\\![\sum_{\eta=x,y,z}\frac{S_{x}S_{y}S_{z}}{S_{\xi}S_{\eta}}\mathbf{A}_{\xi\eta}\frac{\partial{\widetilde{\mathbf{u}}}}{\partial\eta}]\\!]=0,\quad[\\![s{\widetilde{\mathbf{u}}}]\\!]=0.$

Second, introduce the auxiliary variables

$\widetilde{\bm{w}}_{\eta}=\left(\frac{\gamma_{\eta}S_{x}S_{y}S_{z}}{S_{\eta}\gamma_{x}\gamma_{y}\gamma_{z}}-1\right)\widetilde{\mathbf{u}}(x,y,z,s),\quad\widetilde{\bm{v}}_{\eta\xi}=\frac{1}{\left(\alpha+s\right)S_{\xi}}\frac{\partial}{\partial\eta}\left(\widetilde{\mathbf{u}}(x,y,z,s)+\widetilde{\bm{w}}_{\eta}\right),$

$\widetilde{\bm{\psi}}=\frac{1}{s+\alpha}\widetilde{\mathbf{u}}(x,y,z,s),\quad\widetilde{\bm{\phi}}=\frac{s}{(s+\alpha)^{2}}\widetilde{\mathbf{u}}(x,y,z,s),\quad\widetilde{\bm{\theta}}=\frac{1}{s+\alpha}\widetilde{\bm{\phi}}.$

And invert the Laplace transforms, we have

$\begin{split}&\bm{\Phi}^{-1}\left(\frac{\partial^{2}\mathbf{u}}{\partial
t^{2}}+\sum_{\xi=x,y,z}d_{\xi}\frac{\partial}{\partial
t}\left(\mathbf{u}-\alpha{\bm{\psi}}\right)+\prod_{\eta=x,y,z}d_{\eta}\sum_{\xi=x,y,z}\frac{1}{d_{\xi}}\left(\mathbf{u}-\alpha\left({\bm{\psi}}+\bm{\phi}\right)\right)+\prod_{\xi}d_{\xi}\left(\bm{\phi}-\alpha\bm{\theta}\right)\right)\\\
&=\sum_{\xi=x,y,z}\frac{1}{\gamma_{\xi}}\frac{\partial}{\partial\xi}\sum_{\eta=x,y,z}\left(\frac{1}{\gamma_{\eta}}\mathbf{A}_{\xi\eta}\left(\frac{\partial}{\partial\eta}\left({\mathbf{u}}+{\bm{w}}_{\eta}\right)-d_{\xi}\bm{v}_{\eta\xi}\right)\right),\end{split}$
(50) $\begin{split}\frac{\partial}{\partial
t}\bm{\psi}=\mathbf{u}-\alpha\bm{\psi},\quad\frac{\partial}{\partial
t}\bm{\phi}=\mathbf{u}-\alpha\left(\bm{\psi}+\bm{\phi}\right),\quad\frac{\partial}{\partial
t}\bm{\theta}=\bm{\phi}-\alpha\bm{\theta},\end{split}$ (51)
$\begin{split}\frac{\partial}{\partial
t}\bm{w}_{\eta}=\sum_{\xi\neq\eta}d_{\xi}\mathbf{u}+\prod_{\xi\neq\eta}d_{\xi}\bm{\psi}-\alpha\bm{w}_{\eta},\quad\frac{\partial}{\partial
t}\bm{v}_{\eta\xi}=\frac{\partial}{\partial\eta}\left({\mathbf{u}}+{\bm{w}}_{\eta}\right)-\left(d_{\xi}+\alpha\right)\bm{v}_{\eta\xi},\end{split}$
(52)

subject to the boundary conditions,

$\begin{split}&\frac{1-r_{\xi}}{2}\mathbf{Z}\left(\frac{\partial\mathbf{u}}{\partial
t}+\sum_{\eta\neq\xi}d_{\eta}\left(\mathbf{u}-\alpha{\bm{\psi}}\right)+\prod_{\eta\neq\xi}d_{\eta}\bm{\phi}\right)\\\
&\pm\\\
&\frac{1+r_{\xi}}{2}\sum_{\eta=x,y,z}\left(\frac{1}{\gamma_{\eta}}\mathbf{A}_{\xi\eta}\left(\frac{\partial}{\partial\eta}\left({\mathbf{u}}+{\bm{w}}_{\eta}\right)-d_{\xi}\bm{v}_{\eta\xi}\right)\right)=0,\end{split}$
(53)

and the interface conditions

$[\\![\sum_{\eta=x,y,z}\frac{1}{\gamma_{\eta}}\mathbf{A}_{\xi\eta}\left(\frac{\partial}{\partial\eta}\left({\mathbf{u}}+{\bm{w}}_{\eta}\right)-d_{\xi}\bm{v}_{\eta\xi}\right)]\\!]=0,\quad[\\![\frac{\partial{\mathbf{u}}}{\partial
t}]\\!]=0.$ (54)

The boundary/interface conditions, (24)–(25) and (36) for the second order
wave equation are transformed to (53) and (54) for the second order PML
equations (50)–(52).

###### Remark 3

Note that for the second order systems the boundary/interface conditions,
(24)–(25) and (36), are are transformed in a nontrivial manner to the boundary
conditions (53) and (54) for the PML equations. This is in contrast to the PML
for first order systems where the boundary conditions remain unchanged for the
PML.

Since, we have used the well known PML metric (38), the PML models (45)-(46)
and (50)–(52) can be shown to be analogous to other PML models such as [36,
40, 6, 15, 41, 107, 102, 42]. We will initialize the PML with zero initial
data and terminate the PML (45)-(46) with the boundary conditions (19) and the
PML (50)–(52) with the boundary condition (25). We again note that the PML
absorption functions and auxiliary functions vanish almost everywhere except
in the layers defining the PML where $d_{\xi}\geq 0$.

##  4 Mathematical analysis of PMLs at the continuous PDE level

The PML transformations (38) and (39) have some important mathematical and
physical properties which can be revealed through rigorous mathematical
analysis of PMLs. For a first order system and a second order system modelling
the same physical phenomena, the second order PML (50)–(52) can be rewritten
as the first order PML (45)-(46) model by introducing suitable variables. This
implies that the mathematical properties of the systems (45)-(46) and
(50)–(52) are equivalent. In this review, we will focus on the first order
system, and discuss _perfect matching_ , _well-posedness_ and _stability_.
Analysis of well-posedness and stability of the PML has been an area of active
research [35, 15, 102, 6, 39, 67, 16, 18, 10]. We will try to summarise the
results without going into too much details of the proofs. To simplify the
discussions we will often consider the PML only in the $x$-direction, with
$d_{x}(x)\geq 0$ and $\gamma_{\xi}=1$, $d_{\xi}=0$ for $\xi=y,z$.

###  4.1 Perfect matching

We will now consider a whole space PML problem, with $\Omega=\mathbb{R}^{3}$,
and the PML in the $x$-direction. In the interior $x\leq 0$ we have
$d_{x}(x)=0$, $\gamma_{x}=1$ and in the PML $x>0$, $d_{x}(x)>0$,
$\gamma_{x}(x)>0$. Probably, the most important property of the PML (45)-(46)
is that the equations are perfectly matched to the underlying hyperbolic PDEs.
This means that the restriction of the general solutions to (45)-(46) in the
interior of the domain, $x\leq 0$, coincides with the general solutions to
(1). Therefore, there are no reflections as waves propagate from the interior,
$x\leq 0$ into the layer $x>0$.

We will formalise the discussion below.

###### Definition 2

Let $\mathbf{U}$ denote the general solution of the wave equation (1) in the
unbounded domain and $\mathbf{U}_{\mathrm{pml}}$ denote the general solution
of the PML (45)-(46). The equations (45)-(46) are perfectly matched to the
wave equation (1) if

$\displaystyle\mathbf{U}_{\mathrm{pml}}=\mathbf{U},\quad x\leq 0,\quad\forall
t\geq 0.$

There are two standard methods [25, 23, 10, 35] that have been used to study
the perfect matching property of the PML. The approach [25, 23] uses plane
wave analysis and only accounts for propagating modes. The technique [10, 35],
which is rooted in the construction of the general solutions to the wave
equations in the Laplace–Fourier space, is more general since it includes both
the propagating mode regime and the evanescent mode regime.

Consider the wave equation (1), and take the Laplace transform in time and
Fourier transforms in the tangential directions

$\begin{split}\mathbf{P}^{-1}s\widehat{\widetilde{\mathbf{U}}}\left(x,ik_{y},ik_{z},s\right)&=\mathbf{A}_{x}\frac{d}{d{x}}{\widehat{\widetilde{\mathbf{U}}}\left(x,ik_{y},ik_{z},s\right)}-\sum_{\xi=y,z}ik_{\xi}\mathbf{A}_{\xi}\widehat{\widetilde{\mathbf{U}}}\left(x,ik_{y},ik_{z},s\right)\\\
\end{split}$ (55)

For (55) we can construct modal solutions

$\displaystyle\widehat{\widetilde{\mathbf{U}}}\left(x,ik_{y},ik_{z},s\right)=\widehat{\widetilde{\mathbf{U}}}_{0}\left(ik_{y},ik_{z},s\right)e^{\lambda
x}$

where $\lambda$ satisfies

$\displaystyle\left|\Upsilon\left(\lambda,ik_{y},ik_{z},s\right)\right|=0,\quad\Upsilon\left(\lambda,ik_{y},ik_{z},s\right)=s\mathbf{P}^{-1}+\sum_{\xi=y,z}ik_{\xi}\mathbf{A}_{\xi}-\lambda\mathbf{A}_{x}.$
(56)

and $\widehat{\widetilde{\mathbf{U}}}_{0}\left(ik_{y},ik_{z},s\right)$ is the
eigenfunction of $\Upsilon\left(\lambda,ik_{y},ik_{z},s\right)$. Similarly,
for the PML (39) we can construct modal solutions

$\displaystyle\widehat{\widetilde{\mathbf{U}}}_{\mathrm{pml}}\left({x},ik_{y},ik_{z},s\right)=\widehat{\widetilde{\mathbf{U}}}_{0}\left(ik_{y},ik_{z},s\right)e^{\lambda\widetilde{x}},$

with $\lambda$ satisfying (56). Note that in the interior $x\leq 0$ we have
$d_{x}=0$, $\gamma_{x}=1$, yielding $x=\widetilde{x}$, and we have

$\displaystyle\widehat{\widetilde{\mathbf{U}}}_{\mathrm{pml}}\left({x},ik_{y},ik_{z},s\right)=\widehat{\widetilde{\mathbf{U}}}\left(x,ik_{y},ik_{z},s\right),\quad\forall
x\leq 0.$

The solutions
$\widehat{\widetilde{\mathbf{U}}}\left({x},ik_{y},ik_{z},s\right)$ and
$\widehat{\widetilde{\mathbf{U}}}_{\mathrm{pml}}\left({x},ik_{y},ik_{z},s\right)=\widehat{\widetilde{\mathbf{U}}}\left(\widetilde{x},ik_{y},ik_{z},s\right)$
are identical in the interior $x\leq 0$. The solutions are perfectly matched
by construction.

###### Remark 4

Note the similarity between (39) and (55). Since the eigenfunction
$\widehat{\widetilde{\mathbf{U}}}_{0}\left(ik_{y},ik_{z},s\right)$ is
unchanged the PML (45)-(46) corresponds to (39), a complex change of
coordinates in the Laplace space. Any wave propagating to the right in the
interior corresponds to a wave propagating to the right in the PML. It then
follows by construction that (45)-(46) is perfectly matched to the underlying
hyperbolic PDE (1), see [10, 35].

###  4.2 Well-posedness

We begin by introducing the notions of well-posedness and stability. By a
well-posed problem, we mean that there is a unique solution which depends
continuously on the data of the problem. To be precise, consider the Cauchy
problem

$\displaystyle\frac{\partial\mathbf{Q}}{\partial
t}=\mathscr{D}\left({\partial}/{\partial\xi},d_{\xi},\gamma_{\xi},\alpha\right)\mathbf{Q},\quad\mathbf{Q}(\bar{\xi},0)=\mathbf{Q}_{0}\left(\bar{\xi}\right),\quad\bar{\xi}=(x,y,z)$
(57)

Here, $\mathbf{Q}=\left(\mathbf{U},\mathbf{w}_{\xi}\right)^{T}$, the
differential symbol
$\mathscr{D}\left({\partial}/{\partial\xi},d_{\xi},\gamma_{\xi},\alpha\right)$
denotes the spatial operator, and for the PML (45)-(46) it is given by

$\begin{split}\mathscr{D}\left({\partial}/{\partial\xi},d_{\xi},\gamma_{\xi},\alpha\right)=\underbrace{\begin{pmatrix}\sum_{\xi=x,y,z}\frac{1}{\gamma_{\xi}}\mathbf{P}\mathbf{A}_{\xi}\frac{\partial{}}{\partial\xi},&\mathbf{0}\\\
\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{}}{\partial\xi},&\mathbf{0}\end{pmatrix}}_{\text{principal
part}}-\begin{pmatrix}\mathbf{0},&\sum_{\xi=x,y,z}d_{\xi}\mathbf{P}\\\
\mathbf{0},&\left(d_{\xi}+\alpha\right)\mathbf{I}\end{pmatrix}.\end{split}$
(58)

###### Definition 3

The Cauchy problem (57) is weakly (resp. strongly) well-posed if for every
$t\geq 0$

$\|\mathbf{Q}\|\leq Ke^{\kappa t}\|\mathbf{Q}_{0}\|_{H^{\omega}},$

for $\mathbf{Q}_{0}$ giving in the Sobolev space $H^{\omega}$, ${\omega}>0$
(resp. ${\omega}=0$), where $\kappa$ and $K$ are constants independent of
$\mathbf{Q}_{0}$.

For the PML (45)-(46), we consider the principal part and the coefficient
matrices

$\widehat{\mathbf{A}}_{\xi}=\begin{pmatrix}\frac{1}{\gamma_{\xi}}\mathbf{P}\mathbf{A}_{\xi},&\mathbf{0}\\\
\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi},&\mathbf{0}\end{pmatrix}\sim\begin{pmatrix}\frac{1}{\gamma_{\xi}}\mathbf{P}\mathbf{A}_{\xi},&\mathbf{0}\\\
\frac{1}{\gamma_{\xi}}\mathbf{P}\mathbf{A}_{\xi},&\mathbf{0}\end{pmatrix},$

where $\sim$ denotes the similarity of matrices. Since the underlying system
(1), without the PML, is strongly hyperbolic the coefficient matrices
$\mathbf{P}\mathbf{A}_{\xi}$ have real eigenvalues and a complete system of
eigenvectors. Thus, by inspection the matrices $\widehat{\mathbf{A}}_{\xi}$
have real eigenvalues with multiple eigenvectors. Therefore, by standard
definitions [78], the PML (45)-(46) is weakly hyperbolic, and thus weakly
well-posed. For some systems such as the acoustic wave equation [36] or the
Maxwell’s equations [3, 36], it is possible to rewrite the equations, by
introducing new variables such that the PML is strongly hyperbolic, and thus
strongly well-posed. However, the resulting strongly hyperbolic system is
linearly equivalent to the original weakly hyperbolic systems, and their
general solutions are identical [3, 36, 46].

For time-dependent problems, it is not sufficient that the PML is well-posed,
it must also be stable in the sense that the solution remains bounded after a
sufficiently long time. By Definition 3, a well-posed problem can support
exponentially growing solutions. This is of course undesirable of an absorbing
model. Any growth in the layer can propagate into the computational domain and
pollute the solutions everywhere. In order for the PML to be useful it must
therefore be stable.

###  4.3 Stability

The stability of the PML has attracted considerable attention in the
literature, [35, 15, 102, 6, 39]. The perfectly matched layer is indeed a
variable coefficient problem, but many classical techniques require constant
coefficients, and therefore one often considers corresponding frozen
coefficient problems. Further simplifications involve considering sub-problems
where only one PML damping is nonzero, eg. the $x$-dependent PML strip with
$d_{x}>0$ and $d_{\xi}=0$ for $\xi\neq x$.

####  4.3.1 Fourier analysis for the PML Cauchy problem

Consider the Cauchy problem for Equations (45)-(46), with constant
coefficients, $d_{x}=d>0$, $d_{y}=d_{z}=0$, $\alpha_{x}=\alpha\geq 0$,
$\alpha_{y}=\alpha_{z}=0$, and $\gamma_{x}=\gamma>0$. The main idea is to look
for wave-like solutions of the form

$\displaystyle\mathbf{Q}=\mathbf{Q}_{0}e^{ik_{x}x+ik_{y}y+ik_{z}z+st},\quad
s=\Re{s}+i\Im{s},$ (59)

where $s\in\mathbb{C}$ is the unknown eigenvalue and will be determined. If
there is an eigenvalue with a positive real part the PML is unstable, since
this will correspond to a plane wave solution with exponentially growing
amplitude. We introduce
$|\mathbf{k}|=\sqrt{\left(k_{x}/\gamma_{x}\right)^{2}+k_{y}^{2}+k_{z}^{2}}$
and the scaling

$\lambda=\frac{s}{|\mathbf{k}|},\quad
k_{1}=\frac{k_{x}}{\gamma_{x}|\mathbf{k}|},\quad
k_{2}=\frac{k_{y}}{|\mathbf{k}|},\quad
k_{3}=\frac{k_{z}}{|\mathbf{k}|},\quad\epsilon=\frac{d}{|\mathbf{k}|},\quad\nu=\frac{\alpha}{|\mathbf{k}|},\quad
S_{x}\left(\lambda,\epsilon,\nu\right)=1+\frac{\epsilon}{\lambda+\nu}.$

Correspondingly, if there are $\Re{\lambda}>0$, the PML is unstable.

We insert the plane wave solution (59) in the PML and obtain the dispersion
relation

$\displaystyle
F\left(-i\lambda,\frac{1}{S_{x}\left(\lambda,\epsilon,\nu\right)}k_{1},k_{2},k_{3}\right):=\Big{|}-i\lambda\mathbf{I}+\frac{1}{S_{x}\left(\lambda,\epsilon,\nu\right)}k_{x}\mathbf{P}\mathbf{A}_{x}+\sum_{\xi=y,z}{k_{\xi}}\mathbf{P}\mathbf{A}_{\xi}\Big{|}=0,$
(60)

with $S_{x}\left(\lambda,\epsilon,\nu\right)=1+\frac{\epsilon}{\lambda+\nu}$.
The scaled eigenvalue $\lambda$ is the root of the complicated nonlinear
dispersion relation (60) for the PML. The roots $\lambda$ can be difficult to
determine. However standard perturbation arguments yield the following well-
known result, the so-called _geometric stability condition_ , see Definition 1
and [35, 15, 6].

######  Theorem 1 (Necessary condition for stability)

Consider the constant coefficient PML, with $d_{x}=d>0$,
$\alpha_{x}=\alpha\geq 0$ and $\gamma_{x}=\gamma>0$. If the geometric
stability condition in the $x$-direction is violated, then there are growing
modes with $\Re{\lambda}>0$ at all sufficiently high frequencies.

The theorem states that if the underlying medium violates the geometric
stability condition, then there are high frequency wave modes with
exponentially growing amplitude. See also Figure 1 for a schematic explanation
of Theorem 1.

###### Remark 5

At intermediate frequencies the parameter $\gamma_{x}$ can be used to
stabilise unstable modes, by compressing the grid using $\gamma_{x}<1$.
However at sufficiently high frequencies the instability persists. See [35,
41] for details.

$S$$V_{p}$$V_{g}$$k_{y}/\omega$$k_{x}/\omega$$\mathbf{Unstable}$$S$$V_{p}$$V_{g}$$k_{y}/\omega$$k_{x}/\omega$$\mathbf{Stable}$
Figure 1: The geometric stability condition of the PML

In Figure 2 we display the slowness diagrams of different wave media. Note
that the geometric stability condition is satisfied for the first three
diagrams in Figure 2(a)–(c). The last diagram Figure 2(d) which is an
anisotropic elastic medium violates the geometric stability condition in
certain parts of the slowness curve.

(a) Acoustic medium

(b) Isotropic elastic wave solid

(c) AM1: Stable anisotropic elastic solid

(d) AM2: Unstable anisotropic elastic solid

Figure 2: Slowness diagrams showing the solutions of the dispersion relation
for different wave media in two space dimensions.

A typical example where the geometric stability condition is violated is in
advective acoustics or linearised Euler equations with non-vanishing mean flow
[70, 105, 84, 1, 71, 63, 33]. Other examples often include certain crystals
and anisotropic elastic media [75, 35, 41, 15]. See for example the
anisotropic material AM2 displayed in Figure 2 d). For advective acoustics or
linearised Euler equations with non-vanishing mean flow, a space-time linear
transformation can allow for the derivation of stable PML, avoiding the
catastrophic exponential growth for the Cauchy PML [1, 71, 63, 33]. There is
however no cure yet for general systems such as linear anisotropic
elastodynamic equation.

In isotropic media and many anisotropic media the geometric stability
condition is satisfied. In particular, it can be shown for isotropic elastic
materials that there are no growing modes for any frequency. However, in
numerical simulations instability may persist even when the mode analysis
indicates that the Cauchy problem is stable. One possibility is that the
introduction of boundary conditions can ruin the stability or well-posedness
of the PML [102]. A second possible source of numerical growth can result from
unstable numerical approximation of the PML IBVP. Next we will review the
theory of well-posedness and stability for PML IBVP.

####  4.3.2 Stability analysis of PML IBVP

The mode analysis for PMLs has been extended to include boundary conditions,
see [35, 42, 36, 40]. The central result is that the PML IBVP will not support
exponentially growing solutions as long as neither the underlying undamped
IBVP nor the PML Cauchy problem support exponentially growing solutions. We
will summarise the analysis below. For more elaborate discussion we refer the
reader to the references above.

As before we assume constant coefficients and consider a PML with damping only
in the horizontal $(x)$ direction. We analyze a left half-space problem with a
boundary at $x=0$, and a lower half-space problem with a boundary at $y=0$,
seperatly. In both cases the boundary condition is of the same form as in the
original equation. We will look for simple wave like solutions

$\displaystyle\mathbf{U}=\bm{\phi}_{0}(x)e^{ik_{y}y+ik_{z}z+st},\quad\int_{-\infty}^{0}|\bm{\phi}_{0}(x)|^{2}dx<\infty,$
(61)

and

$\displaystyle\mathbf{U}=\bm{\phi}_{0}(y)e^{ik_{x}x+ik_{z}z+st},\quad\int_{-\infty}^{0}|\bm{\phi}_{0}(y)|^{2}dy<\infty,$
(62)

respectively, where $s=\Re{s}+i\Im{s}$.

We consider only physical media where the _geometric stability condition_ is
satisfied.

##### PML IBVPs.

We consider the PML in a domain with boundaries and boundary conditions, that
is 2a) left half-space problem, and 2b) lower half-space problem.

As above the PML is not symmetric, we consider first the boundary terminating
the PML, to be precise
2a) The left half-spce PML problem:
$(x,y,z)\in(-\infty,0)\times(-\infty,\infty)\times(-\infty,\infty)$, with
$\mathbf{L}\mathbf{U}=0$, $x=0$. We look for simple wave-like solutions of the
form

Next, we consider
2b) The lower half-space PML problem:
$(x,y,z)\in(-\infty,\infty)\times(-\infty,0)\times(-\infty,\infty)$, with
$\mathbf{L}\mathbf{U}=0$, $y=0$. We also look for simple wave solutions of the
form

###### Definition 4

The IBVP PML 2a) or 2b) is not stable if there are nontrivial solutions
$\mathbf{U}$ of the form (61) or (62) with $\Re{s}>0$.

We make the following technical assumption

###### Assumption 1

In the absence of the PML, $d_{x}=0$, the IBVPs are stable, that is _all
nontrivial eigen-pairs are in the stable half of the complex plane,
$\Re{s}\leq 0$_.

Note that the underlying IBVPs without the PML, $d_{x}=0$, satisfy strictly
non-growing energy estimates, (22) and (29). Therefore, Assumption 1 is
satisfied, otherwise the energy in the medium will grow, thereby contradicting
the energy estimate. The results below would determine how the eigenvalues
will be perturbed when the PML is present, $d_{x}>0$.

###### Theorem 2

Consider the constant coefficient PML IBVP 2a) with $d_{x}=d\geq 0$ subject to
the boundary condition $\mathbf{L}\mathbf{U}=0$ at $x=0$. If Assumption 1
holds and the geometric stability condition is satisfied, then the IBVP has no
nontrivial solution of the form (61) with $\Re{s}>0$.

###### Theorem 3

Consider the constant constant coefficient lPML IBVP 2b) with $d_{x}=d\geq 0$,
subject to the boundary condition $\mathbf{L}\mathbf{U}=0$ at $y=0$. If
Assumption 1 holds and the geometric stability condition is satisfied then the
PML damping $d_{x}>0$ will move all eigenvalues further into the stable
complex plane, $\Re{s}<0$.

The proofs of Theorems 2 and 3 can be adapted from [36, 40, 42].

Theorems 1–3 together prove that if the geometric stability condition is
satisfied the constant coefficient PML IBVP does not support exponentially
growing solutions. However these results are too technical to be extended to
the analysis of numerical approximations of the PML IBVP.

Below, we will review extensions of the results using the energy method. The
energy estimates enable the development of stable and accurate numerical
methods for the PML using SBP finite difference methods and DG methods. More
detailed discussions can be found in [36, 40, 49, 39].

###### Remark 6

We remark that the analysis above and Theorems 1–3 hold for both first order
systems and second order hyperbolic systems. In particular some of these
results, for example Theorem 3, were initially formulated for the 2D elastic
wave equation in second order form [42], however the results hold for first
order systems as well [40] and extend to 3D [49, 39].

###### Remark 7

We also note that Theorems 1–3 proves that the PML is stable in the presence
of boundary waves such as Rayleigh waves propagating on the surface of elastic
solids. The proof of stability of the PML in the presence of interface waves
such as Stoneley waves and Scholte waves is still an open problem. However, we
believe that a similar technique for boundary waves can be used to establish
stability for interface wave modes.

###  4.4 Energy methods

The PML is generally asymmetric. It becomes extremely difficult to derive
energy estimates that can be used to design stable numerical methods for the
PML in truncated domains. In many settings, a straight forward approach would
yield an exponentially growing estimate in physical space. This would indicate
well-posedness, but exponential growth is not optimal for an absorbing model.
A well-functioning PML should not support exponential growth.

We note that at constant coefficient an energy estimate for the acoustic wave
equation in second order form became possible [12] for a specific time-domain
formulation of the PML. However, we will briefly review the energy method in
the Laplace space in general media, which is applicable to most PML models and
useful for developing stable and accurate numerical methods.

To begin with we reformulate the IBVP by introducing new variables by
subtracting the product of the initial function and $e^{-t}$ from each of the
unknown functions. The new unknown functions satisfy the same system as the
original unknowns but with homogeneous initial data and inhomogeneous source
terms in the equations. Denote the source terms by
$\mathbf{F}_{Q}(x,y,z,t)=(\mathbf{F}_{U}(x,y,z,t),\mathbf{F}_{w_{\xi}}(x,y,z,t))^{T}.$
Laplace transformation in time of the PML equations (45)-(46), yields

$\begin{split}\mathbf{P}^{-1}s\widetilde{\mathbf{U}}=\sum_{\xi=x,y,z}\left[\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}-d_{\xi}\left(\xi\right)\widetilde{\mathbf{w}}_{\xi}\right]+\widetilde{\mathbf{F}}_{U},\end{split}$
(63)
$\begin{split}s\widetilde{\mathbf{w}}_{\xi}=\frac{1}{\gamma_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}-\left(\alpha_{\xi}\left(\xi\right)+d_{\xi}\left(\xi\right)\right)\widetilde{\mathbf{w}}_{\xi}+\widetilde{\mathbf{F}}_{w_{\xi}},\end{split}$
(64)

Further, we use (64) and eliminate the auxiliary variables
$\widetilde{\mathbf{w}}_{\xi}$ in (63). We obtain

$\begin{split}\mathbf{P}^{-1}s\widetilde{\mathbf{U}}=\sum_{\xi=x,y,z}\frac{1}{S_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}+\mathbf{P}^{-1}\widetilde{\mathbf{F}},\quad\widetilde{\mathbf{F}}=\mathbf{P}\left(\widetilde{\mathbf{F}}_{U}-\sum_{\xi=x,y,z}\frac{d_{\xi}}{s+\alpha_{\xi}+d_{\xi}}\widetilde{\mathbf{F}}_{w_{\xi}}\right).\end{split}$
(65)

Transforming the boundary conditions (19) and the interface conditions (32)
gives

$\begin{split}L\widetilde{\mathbf{U}}=0,\quad\xi=\pm 1,\end{split}$ (66)

such that

$\frac{1}{2}\widetilde{\mathbf{U}}^{\dagger}\mathbf{A}_{\xi}\widetilde{\mathbf{U}}\leq
0,$

and

$\displaystyle\mathbb{I}\mathbb{T}\left(\mathbf{A}_{\xi},\widetilde{\mathbf{U}}\right)=0,\quad\xi=0,$
(67)

with
$\mathrm{IT}\left(\widetilde{\mathbf{U}}^{\dagger}\mathbf{A}_{\xi}\widetilde{\mathbf{U}}\right)=0$,
and $\widetilde{\mathbf{U}}^{\dagger}$ denotes the conjugate transpose of
$\widetilde{\mathbf{U}}$.

For Laplace transformed systems, it is necessary to extend the weighted scalar
product (4) and the corresponding norm to complex functions. For complex
functions, we define the weighted scalar product

$\displaystyle\left(\widetilde{\mathbf{U}},\widetilde{\mathbf{F}}\right)_{{P}}=\int_{\Omega}\frac{1}{2}\left[\widetilde{\mathbf{U}}^{\dagger}\mathbf{P}^{-1}\widetilde{\mathbf{F}}\right]dxdydz,\quad\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}^{2}=\left(\widetilde{\mathbf{U}},\widetilde{\mathbf{U}}\right)_{{P}}.$
(68)

Again, $\widetilde{\mathbf{U}}^{\dagger}$ denotes the complex conjugate
transpose of $\widetilde{\mathbf{U}}$.

As before we will consider the acoustic wave equation and proceed later to the
general linear hyperbolic system system. We will now review the results
presented in [49, 39]

####  4.4.1 Energy estimate of the PML for the acoustic wave equation

In [39], the PML for the acoustics wave equation was considered with the
standard metric parameters, $\gamma_{\xi}=1$, $\alpha=0$ and $d_{\xi}(\xi)\geq
0$. The estimate can be extended to the general metric parameters (38) used in
this review. Consider specifically the acoustic wave equation with
$\mathbf{U}=\left(p,\mathbf{v}\right)^{T}$, having

$\begin{split}\frac{1}{\kappa}{s\widetilde{p}}+{\gradient}_{d}\cdot\widetilde{\mathbf{v}}=\frac{1}{\kappa}\widehat{F}_{p},\quad\rho{s\widetilde{\mathbf{v}}}+{\gradient}_{d}\widetilde{p}=\rho\mathbf{\widehat{f}},\end{split}$
(69)

where
${\gradient}_{d}=\left(1/S_{x}\partial/\partial{{x}},1/S_{y}\partial/\partial{{y}},1/S_{z}\partial/\partial{{z}}\right)^{T}$
where $\widehat{F}_{p}$ and $\mathbf{\widehat{f}}$ are source terms that
depend on the initial condition. Further, we eliminate the velocity fields, in
(69), having

$\begin{split}\frac{1}{\kappa}{s^{2}\widetilde{p}}-{\gradient}_{d}\cdot\left(\frac{1}{\rho}{\gradient}_{d}\widetilde{p}\right)&=\frac{s}{\kappa}\widehat{F}_{p}-{\gradient}_{d}\cdot\mathbf{\widehat{f}},\\\
\end{split}$ (70)

with the boundary conditions at ${\xi}=\pm 1$

$\frac{1+r_{\xi}}{2Z}s\widetilde{p}\pm\frac{1-r_{\xi}}{2}\frac{1}{S_{\xi}}\frac{\partial\widetilde{p}}{\partial\xi}=0,\quad\text{with}\quad|r_{\xi}|\leq
1.$ (71)

The jump condition (33) and (41) translate to

$\displaystyle\widetilde{p}^{-}=\widetilde{p}^{+}=\widetilde{p},\quad[\\![\frac{1}{S_{\xi}\rho}\frac{\partial\widetilde{p}}{\partial\xi}]\\!]=0.$
(72)

For $s\neq 0$ and $d_{\xi}\geq 0$, denote $s=a+ib$ and

$\displaystyle\real\left(\frac{\left(sS_{\xi}\right)^{*}}{S_{\xi}}\right)=a+\epsilon_{\xi}(s,d_{\xi}),\quad\epsilon_{\xi}(s,d_{\xi}):=d_{\xi}b^{2}\frac{(a+\alpha)^{2}+b^{2}+(a+\alpha)a+\alpha
d_{\xi}}{\left((a+\alpha)(\alpha+a+d_{\xi})\right)^{2}+(d_{\xi}b)^{2}}\geq 0,$
(73)

with

$\displaystyle\real\left(\frac{1}{S_{\xi}}\right)=\frac{(a+\alpha)(a+\alpha+d_{\xi})+b^{2}}{\gamma_{\xi}\left(\left(a+\alpha+d_{\xi}\right)^{2}+b^{2}\right)}>0.$
(74)

Define the energies

$\begin{split}&\widetilde{E}^{2}_{p}\left(s,d_{\xi}\right)=\Big{\|}{s}\widetilde{p}\Big{\|}^{2}_{1/\kappa}+\sum_{\xi=x,y,z}\Big{\|}\frac{1}{S_{\xi}}\frac{\partial\widetilde{p}}{\partial{\xi}}\Big{\|}^{2}_{1/\rho}>0,\\\
&\widetilde{E}_{f}^{2}\left(s,d_{\xi}\right)=\Big{\|}{s}\widehat{F}_{p}\Big{\|}^{2}_{1/\kappa}+\sum_{\xi=x,y,z}\Big{\|}\frac{1}{S_{\xi}}\frac{\partial\widehat{f}_{\xi}}{\partial{\xi}}\Big{\|}_{\kappa}^{2}>0.\end{split}$
(75)

The following theorem and corollary were proven in [39]

###### Theorem 4

Consider the PML equation in the Laplace space (70) with constant damping
$d_{\xi}\geq 0$, and source terms $\widetilde{\mathbf{F}}(x,y,z,s)$,
$\widetilde{\mathbf{f}}$, subject to homogeneous initial data, the boundary
conditions (71), and the jump condition (72) at discontinuities in material
parameters. Let the energy norms $\widetilde{E}_{p}\left(s,d_{\xi}\right)>0$,
$\widetilde{E}_{f}\left(s,d_{\xi}\right)>0$ in the Laplace space be defined in
(75). For any $\Re{s}=a>0$ we have

$\begin{split}&a\widetilde{E}_{p}^{2}\left(s,d_{\xi}\right)+\sum_{\xi=x,y,z}\Big{\|}\frac{1}{S_{\xi}}\frac{\partial\widetilde{p}}{\partial{\xi}}\Big{\|}^{2}_{\epsilon_{\xi}/\rho}+\widetilde{\mathrm{BT}}\left(s,d_{\xi}\right)\leq
2\widetilde{E}_{p}\left(s,d_{\xi}\right)\widetilde{E}_{f}\left(s,d_{\xi}\right),\\\
&\widetilde{\mathrm{BT}}=\int_{-1}^{1}\int_{-1}^{1}\sum_{\xi=x,y,z}\real\left(\frac{1}{S_{\xi}}\right)\left(\frac{s^{*}}{S_{\xi}\rho}\widetilde{p}^{*}\frac{\partial\widetilde{p}}{\partial\xi}\Big{|}_{\xi=-1}-\frac{s^{*}}{S_{\xi}\rho}\widetilde{p}^{*}\frac{\partial\widetilde{p}}{\partial\xi}\Big{|}_{\xi=1}\right)\frac{dxdydz}{d\xi}\geq
0.\end{split}$ (76)

###### Corollary 1

Consider the same PML IBVP, as in Theorem 4. The IBVP, with $d_{\xi}\geq 0$,
is asymptotically stable in the sense that no exponentially growing solutions
are supported.

####  4.4.2 Energy estimate of the PML for general first order systems

The discussions here are the contributions from [49]. For general systems we
consider the scaled equation

$\begin{split}&\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=\sum_{\xi=x,y,z}\frac{S_{x}}{S_{\xi}}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}+\mathbf{P}^{-1}\widetilde{\mathbf{F}},\\\
&\widetilde{\mathbf{F}}=\mathbf{P}\left(S_{x}\widetilde{\mathbf{F}}_{U}-\sum_{\xi=x,y,z}\frac{d_{\xi}S_{x}}{s+\alpha_{\xi}+d_{\xi}}\widetilde{\mathbf{F}}_{w_{\xi}}\right),\end{split}$
(77)

with the boundary conditions (66). The form (77) is chosen so that with
suitable assumptions the spatial operator remains in its original form, and
allows using the standard integration by parts methodology. For an energy
estimate we also need the coefficient $\left(sS_{x}\right)$ in the left hand
side to have a positive real part. Not that

$\displaystyle\Re{\left(sS_{x}\right)}={\gamma_{x}}\left(a+\left(\frac{a\left(a+\alpha\right)+b^{2}}{{|s+\alpha|^{2}}}\right)d_{x}\right)\geq
a\gamma_{x},$ $\displaystyle
0<\frac{1}{\Re{\left(sS_{x}\right)}}=\frac{1}{\gamma_{x}\left(a+\left(\frac{a\left(a+\alpha\right)+b^{2}}{{|s+\alpha|^{2}}}\right)d_{x}\right)}\leq\frac{1}{a\gamma_{x}}.$

##### The PML strip problem.

Consider specifically the $x$-dependent PML strip problem, that is (45)-(46)
with $d_{x}\geq 0$, $\gamma_{x}=\gamma>0$ and $\partial/\partial\xi=0$,
$d_{\xi}=0$, $\gamma_{\xi}=1$ for $\xi=y,z$. This simplification results in
the 1D PML problem

$\begin{split}\left(sS_{x}\right)\mathbf{P}^{-1}\widetilde{\mathbf{U}}={A}_{x}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial
x}+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$ (78)

Equation (78) lives in a 3D domain defined in (15), but has been simplified by
restricting the initial data and the forcing $\widetilde{\mathbf{F}}$ to
functions that vary only in 1D, the $x$-axis. The following theorem was proven
in [49].

###### Theorem 5

Consider the 1D PML equation (78) in the Laplace space, with piecewise
constant $d_{x}(x)=d\geq 0$, constant grid stretching $\gamma_{x}=\gamma>0$
and $\alpha_{x}=\alpha\geq 0$, subject to the boundary conditions (66), and
the interface condition (67) at discontinuities of $d_{x}(x)$. For all $s$
such that $\real(s)\geq a>0$ we have

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}^{2}\leq\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}\|\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{P}+\widetilde{\mathrm{BT}},\quad\real(sS_{x})\geq
a\gamma>0,\\\
&\widetilde{\mathbf{F}}=\mathbf{P}\left(S_{x}\widetilde{\mathbf{F}}_{U}-\frac{d}{s+\alpha}\widetilde{\mathbf{F}}_{w_{x}}\right),\quad\widetilde{\mathrm{BT}}=\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\widetilde{\mathbf{U}}^{\dagger}{A}_{x}\widetilde{\mathbf{U}}\right]\Big{|}_{x=-1}^{x=1}dydz\leq
0.\end{split}$ (79)

###### Remark 8

It is important to note that similar energy estimates (85) are also valid for
the PML strip problems in the $y$-axis or $z$-axis, that is when
$d_{y}=d>0,d_{x}=d_{z}=0$ or $d_{z}=d>0,d_{x}=d_{y}=0$.

##### The PML edge problem.

Next, we consider specifically the $xy$–edge PML problem, that is (45)-(46)
with $d_{x}=d_{y}=d>0$, $\gamma_{x}=\gamma_{y}=\gamma>0$,
$\alpha_{x}=\alpha_{y}=\alpha\geq 0$ and $d_{z}=0$, $\partial/\partial z=0$.
This simplification results in the 2D $xy$-edge PML problem,

$\begin{split}\left(sS_{x}\right)\mathbf{P}^{-1}\widetilde{\mathbf{U}}=\sum_{\xi=x,y}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$
(80)

Note also that equation (80) lives in a 3D domain defined in (15), but has
been simplified by restricting the initial data and the forcing
$\widetilde{\mathbf{F}}$ to functions that vary only in 2D, the $xy$-plane. We
have the result which was proven in [49].

###### Theorem 6

Consider the 2D PML equation (80) in the Laplace space, with piecewise
constant $d_{x}(x)=d_{y}(y)=d\geq 0$, constant grid stretching
$\gamma_{x}=\gamma_{y}=\gamma>0$ and $\alpha_{x}=\alpha_{y}=\alpha\geq 0$,
subject to the boundary conditions (66), and the interface condition (67) at
discontinuities of $d_{x}(x)$, $d_{y}(y)$. For all $s$ such that $\real(s)\geq
a>0$ we have

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}\leq\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}\|\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{P}+\widetilde{\mathrm{BT}},\quad\real(sS_{x})\geq
a\gamma>0,\\\
&\widetilde{\mathbf{F}}=\mathbf{P}\left(S_{x}\widetilde{\mathbf{F}}_{U}-\sum_{\xi=x,y}\frac{d}{s+\alpha}\widetilde{\mathbf{F}}_{w_{\xi}}\right),\quad\widetilde{\mathrm{BT}}=\sum_{\xi=x,y}\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\widetilde{\mathbf{U}}^{\dagger}\mathbf{A}_{\xi}\widetilde{\mathbf{U}}\right]\Big{|}_{-1}^{1}\frac{dxdydz}{d\xi}\leq
0.\end{split}$ (81)

###### Remark 9

It is also noteworthy that similar energy estimates (81) are valid for the PML
edge problems in the $xz$-edge or $yz$-edge, that is when
$d_{x}=d_{z}=d>0,d_{y}=0$ or $d_{y}=d_{z}=d>0,d_{x}=0$.

##### The PML corner problem.

Consider the corresponding PML corner problem, (45)-(46) with
$d_{x}=d_{y}=d_{z}=d>0$, $\gamma_{x}=\gamma_{y}=\gamma_{z}=\gamma>0$, and
$\alpha_{x}=\alpha_{y}=\alpha_{z}=\alpha>0$. Then, all PML metrics are
identical $S_{y}=S_{x}=S_{z}$. We have $S_{x}/S_{\xi}=1$ and

$\begin{split}\left(sS_{x}\right)\mathbf{P}^{-1}\widetilde{\mathbf{U}}=\sum_{\xi=x,y,z}\mathbf{A}_{\xi}\frac{\partial{\widetilde{\mathbf{U}}}}{\partial\xi}+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$
(82)

###### Theorem 7

Consider the 3D PML equation (82) in the Laplace space, with piecewise
constant $d_{x}(x)=d_{y}(y)=d_{z}(z)=d\geq 0$, constant grid stretching
$\gamma_{x}=\gamma_{y}=\gamma_{z}=\gamma>0$ and
$\alpha_{x}=\alpha_{y}=\alpha_{z}=\alpha\geq 0$, subject to the boundary
conditions (66), and the interface condition (67) at discontinuities of
$d_{x}(x)$, $d_{y}(y)$, $d_{z}(z)$. For all $s$ such that $\real(s)\geq a>0$
we have

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}^{2}\leq\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}\|\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{P}+\widetilde{\mathrm{BT}},\quad\real(sS_{x})\geq
a\gamma>0,\\\
&\widetilde{\mathbf{F}}=\mathbf{P}\left(S_{x}\widetilde{\mathbf{F}}_{U}-\sum_{\xi=x,y,z}\frac{d}{s+\alpha}\widetilde{\mathbf{F}}_{w_{\xi}}\right),\quad\widetilde{\mathrm{BT}}=\sum_{\xi=x,y,z}\int_{\widetilde{\Gamma}}\frac{1}{2}\left[\widetilde{\mathbf{U}}^{\dagger}\mathbf{A}_{\xi}\widetilde{\mathbf{U}}\right]\Big{|}_{-1}^{1}\frac{dxdydz}{d\xi}\leq
0.\end{split}$ (83)

We introduce the energy norms in the physical space

$\begin{split}\|{\mathbf{U}}(\cdot,\cdot,\cdot,t)\|_{P}^{2}=\|\mathscr{L}^{-1}\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}^{2},\quad\|{\mathbf{F}}(\cdot,\cdot,\cdot,t)\|_{P}^{2}=\|\mathscr{L}^{-1}\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{P}^{2},\end{split}$
(84)

###### Theorem 8

Consider the energy estimate in the Laplace space

$\begin{split}&\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}^{2}\leq\frac{1}{a\gamma}\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{P}\|\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{P},\quad\real(s)\geq
a>0.\end{split}$ (85)

For any $a>0$, $\gamma>0$ and $T>0$ we have

$\displaystyle\int_{0}^{T}e^{-2at}\|{\mathbf{U}}(\cdot,\cdot,\cdot,t)\|_{P}^{2}dt\leq\frac{1}{\left(a\gamma\right)^{2}}\int_{0}^{T}e^{-2at}\left(\|{\mathbf{F}}(\cdot,\cdot,\cdot,t)\|_{P}^{2}\right)dt.$
(86)

Since the boundary terms $\widetilde{\mathrm{BT}}$ in (76), (83), (81) and
(79) are never positive, we can use Theorem 8 to invert the the estimates in
Theorem 4, 5, 6 and 7, and get an energy estimate in the physical space. The
estimate seemingly allows the energy of the solution to grow exponentially
with time for general data. However, by choosing $a>0$ in relation to the
length of the time interval of interest, a bound involving only algebraic
growth in time follows.

###### Remark 10

For general systems the energy estimates for the first order formulation,
Theorems 5–7, are not directly applicable to general second order formulation,
in particular when mixed spatial derivatives are present. For general second
order systems, such as displacement formulation of linear elastodynamics, we
may need other simplifying assumptions and model reductions in order to
formulate results that will enable the development of accurate and robust
numerical methods.

##  5 Numerical analysis of the discrete PML

In this section, we review the derivation of stable numerical methods and
numerical analysis of the PML. We will formulate the ideas in a general DG
framework. However, the results extend to finite difference approximation
based on the SBP-SAT method. As shown in [36, 40, 49, 39], the goal is to
design the numerical method for the PML IBVP in a manner that allows the
derivation of discrete energy estimates analogous to the continuous energy
estimates presented in Theorems 4–8.

As presented in [49, 39] we will use the physically motivated numerical fluxes
developed in [43, 50] to patch DG elements together into the global domain.
For the undamped problem (1), the physically motivated numerical flux is
upwind by construction and gives an energy estimate analogous to (22). We
believe that the general technique can be extended other DG methods with a
different but stable numerical flux.

###  5.1 Integral formulation of the boundary and inter-element procedures

Consider the spatial domain $(x,y,z)\in\Omega=[-1,1]\times[-1,1]\times[-1,1]$
and discretise it into $K\times L\times M$ elements, where the $klm$-th
element is denoted by
$\Omega_{klm}=[x_{k},x_{k+1}]\times[y_{l},y_{l+1}]\times[z_{m},z_{m+1}]$, for
$k=1,2,\dots,K$, $l=1,2,\dots,L$, $m=1,2,\dots,M$ with $x_{1}=-1$, $y_{1}=-1$,
$z_{1}=-1$ and $x_{K+1}=1$, $y_{L+1}=1$, $z_{M+1}=1$. The physical element
$\Omega_{lmn}=[x_{k},x_{k+1}]\times[y_{l},y_{l+1}]\times[z_{m},z_{m+1}]$ is
mapped to a reference element $(q,r,s)\in\widetilde{\Omega}=[-1,1]^{3}$ using
the linear transformation

$\displaystyle x=x_{k}+\frac{\Delta{x}_{k}}{2}\left(1+q\right),\quad
y=y_{l}+\frac{\Delta{y}_{l}}{2}\left(1+r\right),\quad
z=z_{m}+\frac{\Delta{z}_{m}}{2}\left(1+s\right),$ (87)

with
$\Delta{x}_{k}=x_{k+1}-x_{k},\quad\Delta{y}_{l}=y_{l+1}-y_{l},\quad\Delta{z}_{m}=y_{m+1}-y_{m}.$
The nonzero metric derivatives and the Jacobian of the transformation are

$q_{x}=\frac{2}{\Delta{x}_{k}},\quad r_{y}=\frac{2}{\Delta{y}_{l}},\quad
s_{z}=\frac{2}{\Delta{z}_{m}},\qquad
J=\frac{\Delta{x}_{k}}{2}\frac{\Delta{y}_{l}}{2}\frac{\Delta{z}_{m}}{2}>0.$

The elemental volume integral yields

$\displaystyle\int_{\Omega_{klm}}f(x,y,z)dxdydz=\int_{\widetilde{\Omega}}f(q,r,s)Jdqdrds,$
(88)

We introduce the scalar products defined by the volume and surface integrals

$\displaystyle\left(\bm{\phi},\bm{\psi}\right)=\int_{\widetilde{\Omega}}\bm{\phi}^{T}\bm{\psi}Jdqdrds,\quad\left(\bm{\phi},\bm{\psi}\right)_{\Gamma_{\xi}}=\int_{\widetilde{\Gamma}}\bm{\phi}^{T}\bm{\psi}J\frac{dqdrds}{d\xi},$
(89)

and adopt the notation

$\displaystyle\mathbf{w}_{q}=\mathbf{w}_{x},\quad\mathbf{w}_{r}=\mathbf{w}_{y},\quad\mathbf{w}_{s}=\mathbf{w}_{z},\quad\mathbf{A}_{q}=\frac{q_{x}}{\gamma_{x}}\mathbf{A}_{x},\quad\mathbf{A}_{r}=\frac{r_{y}}{\gamma_{y}}\mathbf{A}_{y},\quad\mathbf{A}_{s}=\frac{s_{z}}{\gamma_{z}}\mathbf{A}_{z},$
$\displaystyle{d}_{q}={d}_{x},\quad{d}_{r}={d}_{y},\quad{d}_{s}={d}_{z},\quad{\alpha}_{q}={\alpha}_{x},\quad{\alpha}_{r}={\alpha}_{y},\quad{\alpha}_{s}={\alpha}_{z}.$

The elemental integral formulation reads

$\displaystyle\left(\bm{\phi},\mathbf{P}^{-1}\frac{\partial\mathbf{U}}{\partial
t}\right)=\sum_{\xi=q,r,s}\left[\left(\bm{\phi},\mathbf{A}_{\xi}\frac{\partial{\mathbf{U}}}{\partial\xi}-d_{\xi}\left(\xi\right)\mathbf{w}_{\xi}\right)-\left(\left(\bm{\phi}\Big{|}_{\xi=-1},\mathbf{FL}_{\xi}\right)_{\Gamma_{\xi}}+\left(\bm{\phi}\Big{|}_{\xi=1},\mathbf{FR}_{\xi}\right)_{\Gamma_{\xi}}\right)\right],$
(90) $\displaystyle\left(\bm{\phi},\frac{\partial{\mathbf{w}_{\xi}}}{\partial
t}\right)=\left(\bm{\phi},\mathbf{A}_{\xi}\frac{\partial{\mathbf{U}}}{\partial\xi}-\left(\alpha_{\xi}(\xi)+d_{\xi}\left(\xi\right)\right)\mathbf{w}_{\xi}\right)-\underbrace{\theta_{\xi}\left(\left(\bm{\phi}\Big{|}_{\xi=-1},\mathbf{FL}_{\xi}\right)_{\Gamma_{\xi}}+\left(\bm{\phi}\Big{|}_{\xi=1},\mathbf{FR}_{\xi}\right)_{\Gamma_{\xi}}\right)}_{\text{PML
stabilizing flux fluctuations}},$ (91)

where $\mathbf{FL}_{\xi}$ and $\mathbf{FR}_{\xi}$ are physics based flux
fluctuation vectors designed to enforce the boundary and inter-element
conditions. Please see [43, 50, 49, 39] for more details. We note however for
exact solutions that satisfy the IBVP the flux fluctuations vanish identially
$\mathbf{FL}_{\xi}=\mathbf{FR}_{\xi}=0$. The stabilising PML parameters
$\theta_{\xi}$ to be determined by requiring stability of the discrete PML, in
the sense corresponding to Theorems 5, 6, 7 and 8.

###### Remark 11

In [42, 40, 49, 39] we have demonstrated for DG/SBP approximation that
appropriate choice of the stabilising PML parameters $\theta_{\xi}$ is
necessary for the stability of the discrete PML for several systems. These
important results will be summarised in the next sections.

###  5.2 The DG/SBP approximation of the PML

We follow a standard finite element procedure and approximate the elemental
solution by a polynomial interpolant

$\mathbf{U}\left(q,r,s,t\right)=\sum_{i=1}^{N+1}\sum_{j=1}^{N+1}\sum_{k=1}^{N+1}\mathbf{U}_{ijk}(t)\bm{\phi}_{ijk}\left(q,r,s\right),$
(92)

where $\mathbf{U}_{ijk}(t)$, are the elemental degrees of freedom to be
determined, and $\bm{\phi}_{ijk}(q,r,s)$ are the $ijk$-th interpolating
polynomials. We consider tensor products of nodal basis with
$\bm{\phi}_{ijk}(q,r,s)=\mathscr{L}_{i}(q)\mathscr{L}_{j}(r)\mathscr{L}_{k}(s)$,
where $\mathscr{L}_{i}(q)$, $\mathscr{L}_{j}(r)$, $\mathscr{L}_{k}(s)$, are
one dimensional nodal interpolating Lagrange polynomials of degree $N$. We
will only use Gauss-type quadrature rules such that for all polynomial
integrand $f(\xi)$ of degree $\leq 2N-1$, the corresponding one dimensional
rule is exact, $\sum_{m=1}^{N+1}f(\xi_{m})h_{m}=\int_{-1}^{1}f(\xi)d\xi.$

We now make a classical Galerkin approximation by choosing test functions in
the same space as the basis functions, so that the residual is orthogonal to
the space of test functions. We rearrange the elemental degrees of freedom
$[\mathbf{U}_{ijk}(t)]$ row-wise as a vector, $\mathbf{U}(t)$, of length
$n_{f}(N+1)^{d}$ where $d=3$ is the number of space dimensions, $n_{f}$
($n_{f}=4$ for acoustics and $n_{f}=9$ for elasticity) is the number of fields
in the physical variables. We have the semi-discrete approximation

$\begin{split}\mathbf{P}^{-1}\frac{d{\mathbf{U}}}{dt}=\sum_{\xi=x,y,z}\left[\mathbf{A}_{\xi}\mathbf{D}_{\xi}\mathbf{U}-\mathbf{d}_{\xi}\mathbf{w}_{\xi}-\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\mathbf{FL}_{\xi}+{\mathbf{e}_{\xi}(1)}\mathbf{FR}_{\xi}\right)}\right],\end{split}$
(93)
$\begin{split}\frac{d{\mathbf{w}_{\xi}}}{dt}=\mathbf{A}_{\xi}\mathbf{D}_{\xi}\mathbf{U}-\left(\mathbf{d}_{\xi}+\bm{\alpha}_{\xi}\right)\mathbf{w}_{\xi}-\underbrace{\theta_{\xi}\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\mathbf{FL}_{\xi}+{\mathbf{e}_{\xi}(1)}\mathbf{FR}_{\xi}\right)}}_{\text{PML
stabilizing flux fluctuations}}.\end{split}$ (94)

Here the discrete operators are defined by

$\mathbf{D}_{x}=\frac{2}{\Delta{x}_{k}}\left(I_{f}\otimes D\otimes I\otimes
I\right),\quad\mathbf{D}_{y}=\frac{2}{\Delta{y}_{l}}\left(I_{f}\otimes
I\otimes D\otimes
I\right),\quad\mathbf{D}_{z}=\frac{2}{\Delta{z}_{m}}\left(I_{f}\otimes
I\otimes I\otimes D\right),$ (95)

and the discrete norms are given by

$\displaystyle\mathbf{H}_{x}=\frac{\Delta{x}_{k}}{2}\left(I_{f}\otimes
H\otimes I\otimes
I\right),\quad\mathbf{H}_{y}=\frac{\Delta{y}_{l}}{2}\left(I_{f}\otimes
I\otimes H\otimes
I\right),\quad\mathbf{H}_{z}=\frac{\Delta{z}_{m}}{2}\left(I_{f}\otimes
I\otimes I\otimes H\right),$

where

$D=H^{-1}Q,\quad H=\mathrm{diag}[h_{1},h_{2},\cdots,h_{P+1}],\quad
Q_{ij}=\sum_{m=1}^{P+1}h_{m}\mathscr{L}_{i}(q_{m}){\mathscr{L}_{j}^{\prime}(q_{m})},$
(96)

is a spectral difference approximation of the first derivative, in one space
dimension. We also have

$\mathbf{P}=\left({P}\otimes I\otimes I\otimes
I\right),\quad\mathbf{A}_{\xi}=\left({A}_{\xi}\otimes I\otimes I\otimes
I\right),\quad\mathbf{H}=\mathbf{H}_{x}\mathbf{H}_{y}\mathbf{H}_{z},$

where $I$ is the $(N+1)\times(N+1)$ identity matrix, $I_{f}$ is the
$n_{f}\times n_{f}$ identity matrix, and $\otimes$ denotes the Kronecker
product.

We also introduce the projection matrices

$\displaystyle\mathbf{e}_{x}(\eta)=\left(I_{n_{f}}\otimes\bm{e}(\eta)\otimes
I\otimes I\right),\quad\mathbf{e}_{y}(\eta)=\left(I_{n_{f}}\otimes
I\otimes\bm{e}(\eta)\otimes I\right),$
$\displaystyle\mathbf{e}_{z}(\eta)=\left(I_{n_{f}}\otimes I\otimes
I\otimes\bm{e}(\eta)\right),\quad\mathbf{B}_{\eta}(\psi,\xi)=\mathbf{e}_{\eta}(\psi)\mathbf{e}_{\eta}^{T}(\xi),$

where

$\bm{e}(\eta)=[\mathscr{L}_{i}(\eta),\mathscr{L}_{i}(\eta),\cdots,\mathscr{L}_{P+1}(\eta)]^{T}.$

The derivative oopertaors (95) satisfy the discrete SBP property

$\displaystyle\mathbf{D}_{\xi}=-\mathbf{H}_{\xi}^{-1}\mathbf{D}_{\xi}^{T}\mathbf{H}_{\xi}+\mathbf{H}_{\xi}^{-1}\left(\mathbf{B}_{\xi}\left(1,1\right)-\mathbf{B}_{\xi}\left(-1,-1\right)\right),\quad\xi=x,y,z.$
(97)

###### Remark 12

Note that the spatial discrete operators (96), (95) satisfying the SBP
property (97) have been derived using a spectral approach and have full
accuracy, (that is $(P+1)$th accuracy) within the element. Standard SBP finite
difference operators on equidistant grids have also been used with reduced
accuracy close to the boundaries. However, the stability results are
identical.

Introduce the discrete inner products and norm

$\displaystyle\Big{\langle}{\mathbf{U}},{\mathbf{V}}\Big{\rangle}_{\mathbf{H}}={\mathbf{V}}^{T}\mathbf{H}{\mathbf{U}},\quad\|{\mathbf{U}}\|_{\mathbf{H}}^{2}=\Big{\langle}{\mathbf{U}},{\mathbf{U}}\Big{\rangle}_{\mathbf{H}}>0,$

and the discrete energy norm

$\displaystyle\|\mathbf{U}\|^{2}_{hP}:=\Big{\langle}{\mathbf{U}},P^{-1}{\mathbf{U}}\Big{\rangle}_{\mathbf{H}}>0.$

We have

###### Theorem 9

Consider the semi-discrete approximation (93)-(94). When all PML absorption
functions vanish, $d_{\xi}=0$, the solution of the semi-discrete approximation
satisfies the energy identity

$\frac{d}{dt}\|\mathbf{U}\|^{2}_{hP}\leq 0.$

The proof of Theorem 9 can be adapted from [43, 50]. We will not repeat it
here.

Note that Theorem 9 is not applicable when the PML is present, $d_{\xi}\neq 0$
for any $\xi=x,y,z$.

###  5.3 Discrete energy estimate for the PML in Laplace space

Taking the Laplace transform, in time, of the semi-discrete problem (93)-(94)
gives

$\begin{split}\mathbf{P}^{-1}s{\widetilde{\mathbf{U}}}=\sum_{\xi=x,y,z}\left[\mathbf{A}_{\xi}\mathbf{D}_{\xi}\widetilde{\mathbf{U}}-\mathbf{d}_{\xi}\mathbf{w}_{\xi}-\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}\right]+\widetilde{\mathbf{F}}_{U},\end{split}$
(98)
$\begin{split}s\widetilde{\mathbf{w}}_{\xi}=\mathbf{A}_{\xi}\mathbf{D}_{\xi}\widetilde{\mathbf{U}}-\left(\mathbf{d}_{\xi}+\bm{\alpha}_{\xi}\right)\widetilde{\mathbf{w}}_{\xi}-\underbrace{\theta_{\xi}\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}}_{\text{PML
stabilizing flux fluctuations}}+\widetilde{\mathbf{F}}_{w_{\xi}}.\end{split}$
(99)

Next, we use (99) and eliminate the auxiliary variable
$\widetilde{\mathbf{w}}_{\xi}$ from (98). We have

$\begin{split}\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=&\sum_{\xi=x,y,z}\frac{S_{x}}{S_{\xi}}\left[\frac{1}{2}\left(\mathbf{A}_{\xi}\mathbf{D}_{\xi}-\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\mathbf{D}_{\xi}^{T}\mathbf{H}_{\xi}+\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\left(\mathbf{B}_{\xi}\left(1,1\right)\right)-\mathbf{B}_{\xi}\left(-1,-1\right)\right)\widetilde{\mathbf{U}}\right]\\\
-&\sum_{\xi=x,y,z}\frac{S_{x}}{S_{\xi}}\left[\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}\right]+\mathbf{P}^{-1}\widetilde{\mathbf{F}}\\\
+&\underbrace{\sum_{\xi=x,y,z}\frac{\mathbf{d}_{\xi}S_{x}\left(1-\theta_{\xi}\right)}{S_{\xi}\left(s+\alpha_{\xi}\right)}\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}}_{\text{destabilizing
PML flux term}}.\end{split}$ (100)

Note that in (100), we have made use of the discrete integration-by-parts
property (97) of the spatial derivative operator $\mathbf{D}_{\xi}$. The last
term in the right hand side of (100) is a destabilising PML flux term at
element faces. For a DG method in one element or a standard SBP finite
difference approximation in one block the destabilising terms will appear at
external boundaries. And for a multi-element/block DG/SBP discretisation the
destabilising PML flux terms will also appear every internal element
boundaries in the PML. In particular, the DG is an element based method, these
destabilizing flux terms appear almost everywhere in the PML, including at
internal and external element faces. The destabilizing flux terms can be
eliminated using the parameter $\theta_{\xi}$. We set the penalty parameter
$\theta_{\xi}=1$, extending the numerical implementation of the boundary
conditions and inter-element conditions to the auxiliary differential
equations. The last term in the right hand side of (100) vanishes, having

$\begin{split}\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=&\sum_{\xi=x,y,z}\frac{S_{x}}{S_{\xi}}\left[\frac{1}{2}\left(\mathbf{A}_{\xi}\mathbf{D}_{\xi}-\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\mathbf{D}_{\xi}^{T}\mathbf{H}_{\xi}+\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\left(\mathbf{B}_{\xi}\left(1,1\right)\right)-\mathbf{B}_{\xi}\left(-1,-1\right)\right)\widetilde{\mathbf{U}}\right]\\\
-&\sum_{\xi=x,y,z}\frac{S_{x}}{S_{\xi}}\left[\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}\right]+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$
(101)

We will now make the review specific for the acoustic wave equation and the
elastic wave equations.

####  5.3.1 Discrete energy estimate of the PML for the acoustic wave equation

As before we will consider first the acoustic wave equation, where
$\widetilde{\mathbf{U}}=(\mathbf{\widetilde{p}},\widetilde{\mathbf{v}})^{T}$,
and summarise the discrete stability results. More elaborate discussions can
be found in [36, 39]. For simplicity we consider a single element DG
approximation, by restricting (101) to one element, and divide through by
$S_{x}$, we have

$\begin{split}\frac{1}{\kappa}{s\mathbf{\widetilde{p}}}+\widetilde{{\gradient}}_{D}\cdot\widetilde{\mathbf{v}}=&\frac{1}{\kappa}\mathbf{\widehat{f}}_{p}-\sum_{\xi}\frac{1}{S_{\xi}}\mathbf{H}_{\xi}^{-1}{\left(\frac{\mathbf{B}_{\xi}(-1,-1)}{Z}\left(\frac{1-r_{\xi}}{2}Z\mathbf{\widetilde{v}}_{\xi}+\frac{1+r_{\xi}}{2}\mathbf{\widetilde{p}}\right)-\frac{\mathbf{B}_{\xi}(1,1)}{Z}\left(\frac{1-r_{\xi}}{2}Z\mathbf{\widetilde{v}}_{\xi}-\frac{1+r_{\xi}}{2}\mathbf{\widetilde{p}}\right)\right)},\end{split}$
(102)
$\begin{split}\rho{s\widetilde{\mathbf{v}}}+\widetilde{{\gradient}}_{D}\mathbf{\widetilde{p}}&=\rho\mathbf{\widehat{f}}-\sum_{\xi}\frac{\mathbf{n}}{S_{\xi}}\otimes\mathbf{H}_{\xi}^{-1}{\left({\mathbf{B}_{\xi}(-1,-1)}\left(\frac{1-r_{\xi}}{2}Z\mathbf{\widetilde{v}}_{\xi}+\frac{1+r_{\xi}}{2}\mathbf{\widetilde{p}}\right)+{\mathbf{B}_{\xi}(1,1)}\left(\frac{1-r_{\xi}}{2}Z\mathbf{\widetilde{v}}_{\xi}-\frac{1+r_{\xi}}{2}\mathbf{\widetilde{p}}\right)\right)}.\end{split}$
(103)

In (102)–(103), only external boundaries are present. Introducing the modified
discrete operators

$\begin{split}&\widetilde{\mathbf{D}}_{\xi}=\left(\mathbf{D}_{\xi}+\frac{1+r_{\xi}}{2}\mathbf{H}_{\xi}^{-1}\left(\mathbf{B}_{\xi}\left(-1,-1\right)-\mathbf{B}_{\xi}\left(1,1\right)\right)\right),\quad\widetilde{\mathbf{D}}_{0\xi}=\left(\mathbf{D}_{\xi}+\frac{1-r_{\xi}}{2}\mathbf{H}_{\xi}^{-1}\left(\mathbf{B}_{\xi}(-1,-1)-\mathbf{B}_{\xi}(1,1)\right)\right),\\\
&\widetilde{\mathbf{H}}_{\xi}=\mathbf{H}\widetilde{\mathbf{I}}_{\xi},\quad\widetilde{\mathbf{I}}_{\xi}=\left(\mathbf{I}+\frac{c\left(1-r_{\xi}\right)}{2sS_{\xi}}\mathbf{H}_{\xi}^{-1}\left(\mathbf{B}_{\xi}(-1,-1)+\mathbf{B}_{\xi}(1,1)\right)\right)^{-1},\end{split}$
(104)

with $c=Z/\rho$, and eliminating the velocity fields having

$\begin{split}&s^{*}\mathbf{H}s{\bm{\kappa}}^{-1}s\widetilde{\bm{p}}+\sum_{\xi}\left(\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{\xi}\right)^{\dagger}\left(\frac{(s^{*}S_{\xi}^{*})}{\bm{\rho}S_{\xi}}\widetilde{\mathbf{H}}_{\xi}\right)\left(\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{\xi}\right)\widetilde{\bm{p}}\\\
&=|s|^{2}\sum_{\xi}\frac{1+r_{\xi}}{2ZS_{\xi}}{\mathbf{H}}{\mathbf{H}_{\xi}}^{-1}\left(\mathbf{B}_{\xi}{(-1,-1)}+\mathbf{B}_{\xi}{(1,1)}\right)\widetilde{\bm{p}}s^{*}\mathbf{H}{\bm{\kappa}}^{-1}s\widehat{\mathbf{F}}_{p}-\sum_{\xi}s^{*}{\mathbf{H}}\left(\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{0\xi}\right)\widehat{\mathbf{\widetilde{f}}}_{\xi},\end{split}$

where
$\widehat{\mathbf{\widetilde{f}}}_{\xi}=\widetilde{\mathbf{I}}_{\xi}\widehat{\mathbf{{f}}}_{\xi}$.

Note that the discrete operators $\widetilde{\mathbf{D}}_{0\xi}$,
$\widetilde{\mathbf{D}}_{\xi}$, $\widetilde{\mathbf{H}}_{\xi}$ are modified by
the boundary conditions. On GLL collocation nodes
$\widetilde{\mathbf{H}}_{\xi}$ is diagonal. In particular if we consider a
hard wall $r_{\eta}=1$, the boundary terms will vanish, we recover
$\widetilde{\mathbf{H}}_{\xi}=\mathbf{H}=\mathbf{H}^{T}>0$. Consider
$\mathbf{H}=\mathbf{H}^{T}>0$, introduce the discrete scalar product and the
corresponding norm

$\displaystyle\Big{\langle}\widetilde{\mathbf{u}},\widetilde{\mathbf{v}}\Big{\rangle}_{\mathbf{H}}=\widetilde{\mathbf{v}}^{\dagger}\mathbf{H}\widetilde{\mathbf{u}},\quad\|\widetilde{\mathbf{v}}\|_{\mathbf{H}}^{2}=\Big{\langle}\widetilde{\mathbf{v}},\widetilde{\mathbf{v}}\Big{\rangle}_{\mathbf{H}}.$
(105)

For the modified discrete operator
$\Re{\left(\frac{(s^{*}S_{\xi}^{*})}{S_{\xi}}\widetilde{\mathbf{H}}_{\xi}\right)}$
the scalar product can be decomposed into

$\displaystyle\mathbf{v}^{\dagger}\Re{\left(\frac{(s^{*}S_{\xi}^{*})}{S_{\xi}}\widetilde{\mathbf{H}}_{\xi}\right)}\mathbf{v}=a\left\langle\mathbf{v},\mathbf{v}\right\rangle_{\widehat{\mathbf{H}}_{\xi}}+\epsilon_{\xi}(s,d_{\xi})\left\langle\mathbf{v},\mathbf{v}\right\rangle_{\widehat{\mathbf{H}}_{\xi}}+\left(1-r_{\xi}\right)\mathbf{BT}_{\mathrm{num}}^{(\xi)},$
(106)

with
$\left\langle\mathbf{v},\mathbf{v}\right\rangle_{\widehat{\mathbf{H}}_{\xi}}>0,$
and $\mathbf{BT}_{\mathrm{num}}^{(\xi)}>0.$

Introduce the discrete energy norms

$\begin{split}&\mathscr{\widetilde{E}}^{2}_{p}(s)=\Big{\langle}{s}\widetilde{\mathbf{p}},{s}\widetilde{\mathbf{p}}\Big{\rangle}_{\mathbf{H}/\kappa}+\sum_{\xi=x,y,z}\Big{\langle}\frac{1}{S_{\eta}}\widetilde{\mathbf{D}}_{\xi}\widetilde{\mathbf{p}},\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{\xi}\widetilde{\mathbf{p}}\Big{\rangle}_{\widehat{\mathbf{H}}_{\xi}/\rho}>0,\end{split}$
$\begin{split}&\mathscr{\widetilde{E}}_{f}^{2}(s)=\Big{\langle}{s}\widehat{\mathbf{F}}_{p},{s}\widehat{\mathbf{F}}_{p}\Big{\rangle}_{\mathbf{H}/\kappa}+\sum_{\xi=x,y,z}\Big{\langle}\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{0\xi}\widehat{\mathbf{\widetilde{f}}}_{\xi},\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{0\xi}\widehat{\mathbf{\widetilde{f}}}_{\xi}\Big{\rangle}_{\kappa\mathbf{H}}>0.\end{split}$

Note that

$\displaystyle\real\left(\frac{1}{S_{\xi}}\right)=\frac{(a+\alpha_{\xi})(a+\alpha_{\xi}+d_{\xi})+b^{2}}{\gamma_{\xi}\left(\left(a+\alpha_{\xi}+d_{\xi}\right)^{2}+b^{2}\right)}>0.$

If $\Re{s}=a>0$ then $\real\left(\frac{1}{S_{\xi}}\right)>0$ for any
$d_{\xi}\geq 0$. The following theorem was proven in [39]

###### Theorem 10

Consider the one element DGSEM approximation of the PML in the Laplace space,
(104), with constant damping $d_{\xi}\geq 0$, $\alpha_{\xi}=\alpha\geq 0$,
$\gamma_{\xi}=\gamma>0$ and $\Re{s}=a>0$. If $\theta_{\xi}=1$, then we have

$\begin{split}&(\gamma
a)\mathscr{\widetilde{E}}_{p}^{2}(s)+\sum_{\xi=x,y,z}\Big{\langle}\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{\xi}\widetilde{\mathbf{p}},\frac{1}{S_{\xi}}\widetilde{\mathbf{D}}_{\xi}\widetilde{\mathbf{p}}\Big{\rangle}_{\mathscr{\epsilon}_{\xi}\widehat{\mathbf{H}}_{\xi}/\rho}+\sum_{\xi}\frac{1-r_{\xi}}{\rho}\mathbf{BT}_{\mathrm{num}}^{(\xi)}+\mathrm{BT}(s)\leq
2\mathscr{\widetilde{E}}_{p}(s)\mathscr{\widetilde{E}}_{f}(s),\\\
&\mathrm{BT}(s)=|s|^{2}\sum_{\xi}\Re{\left(\frac{1}{S_{\xi}}\right)}\frac{1+r_{\xi}}{2Z}\widetilde{\bm{p}}^{\dagger}\left[{\mathbf{H}}{\mathbf{H}_{\xi}^{-1}}\left(\mathbf{B}_{\xi}{(-1,-1)}+\mathbf{B}_{\xi}{(1,1)}\right)\right]\widetilde{\bm{p}}\geq
0.\end{split}$

Theorem 10 is completely analogous to the continuous counterpart Theorem 4.
The result can be extended to multiple DG elements, see [39] for details.

We will now review the discrete energy estimates for general systems, such as
the linear elastic wave equation, as presented in [49].

####  5.3.2 Discrete energy estimate of the PML for general systems

We introduce the discrete weighted inner product and norm

$\displaystyle\Big{\langle}{\widetilde{\mathbf{U}},\widetilde{\mathbf{F}}}\Big{\rangle}_{h{P}}=\Big{\langle}\widetilde{\mathbf{U}},{\mathbf{P}}^{-1}\widetilde{\mathbf{F}}\Big{\rangle}_{\mathbf{H}},\quad\|\widetilde{\mathbf{U}}\|^{2}_{hP}=\Big{\langle}{\widetilde{\mathbf{U}},\widetilde{\mathbf{U}}}\Big{\rangle}_{h{P}}.$

Following [49], we consider specifically the elastic wave equation in general
media, and for the two elements model we introduce the surface terms,
including the external boundary terms

$\displaystyle\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{-},\widehat{\widetilde{T}}_{\eta}^{-}\right)=\frac{\Delta{y}}{2}\frac{\Delta{z}}{2}\sum_{i=1}^{P+1}\sum_{k=1}^{P+1}\sum_{\eta=x,y,z}\left(\left(\widehat{\widetilde{T}}_{\eta}^{-*}\widehat{\widetilde{v}}_{\eta}^{-}\right)\Big{|}_{-1}\right)_{ik}h_{i}h_{k}\leq
0,$
$\displaystyle\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{+},\widehat{\widetilde{T}}_{\eta}^{+}\right)=\frac{\Delta{y}}{2}\frac{\Delta{z}}{2}\sum_{i=1}^{P+1}\sum_{k=1}^{P+1}\sum_{\eta=x,y,z}\left(\left(\widehat{\widetilde{T}}_{\eta}^{+*}\widehat{\widetilde{v}}_{\eta}^{+}\right)\Big{|}_{1}\right)_{ik}h_{i}h_{k}\leq
0,$

the interface term

$\displaystyle\mathscr{IT}_{s}\left(\widehat{\widetilde{v}}^{\pm},\widehat{\widetilde{T}}^{\pm}\right)=-\frac{\Delta{y}}{2}\frac{\Delta{z}}{2}\sum_{i=1}^{P+1}\sum_{k=1}^{P+1}\sum_{\eta=x,y,z}\left({\widehat{\widetilde{T}}}_{\eta}^{*}[\\![{{\widehat{\widetilde{v}}}_{\eta}]\\!]}\right)_{ik}h_{i}h_{k}\equiv
0,$

and the fluctuation term

$\displaystyle\mathscr{F}_{luc}\left({\widetilde{G}},Z\right)=-\frac{\Delta{y}}{2}\frac{\Delta{z}}{2}\sum_{\eta=x,y,z}\sum_{i=1}^{P+1}\sum_{k=1}^{P+1}\left(\left(\frac{1}{Z_{\eta}}|{\widetilde{G}}_{\eta}|^{2}\right)\Big{|}_{-1}+\left(\frac{1}{Z_{\eta}}|{\widetilde{G}}_{\eta}|^{2}\right)\Big{|}_{1}\right)_{i,k}h_{i}h_{k}\leq
0.$

Here $v^{*}$ denotes the complex conjugate of $v$, and the surface integrals
have been approximated by quadrature rules.

As in the continuous setting, we will consider first the 1D PML strip problem,
and proceed later to the edge and the corner regions.

##### The discrete PML strip problem.

In (101), let $d_{x}=d\geq 0$, $d_{\xi}\equiv 0$,
$\mathbf{D}_{\xi}=\widetilde{\mathbf{FL}}_{\xi}=\widetilde{\mathbf{FR}}_{\xi}=0$
for $\xi=y,z$. We arrive at

$\begin{split}\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=&\left[\frac{1}{2}\left(\mathbf{A}_{x}\mathbf{D}_{x}-\mathbf{H}_{x}^{-1}\mathbf{A}_{x}\mathbf{D}_{x}^{T}\mathbf{H}_{x}+\mathbf{H}_{x}^{-1}\mathbf{A}_{x}\left(\mathbf{B}_{\xi}\left(1,1\right)\right)-\mathbf{B}_{x}\left(-1,-1\right)\right)\widetilde{\mathbf{U}}\right]\\\
-&\left[\mathbf{H}_{x}^{-1}{\left({\mathbf{e}_{x}(-1)}\widetilde{\mathbf{FL}}_{x}+{\mathbf{e}_{x}(1)}\widetilde{\mathbf{FR}}_{x}\right)}\right]+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$
(107)

We can prove the discrete equivalence of Theorem 5.

###### Theorem 11

Consider the semi-discrete PML equation in the Laplace space (107), with
$d_{x}=d\geq 0$, $\alpha_{x}=\alpha\geq 0$, $\gamma_{x}=\gamma>0$ and
$\real(s)\geq a>0$. We have $\real(sS_{x})\geq\gamma a>0$ and

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{-}\left(s\right)\|^{2}_{hP}+\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{+}\left(s\right)\|^{2}_{hP}\leq\|\widetilde{\mathbf{U}}^{-}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{-}\left(s\right)\|_{hP}+\|\widetilde{\mathbf{U}}^{+}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{+}\left(s\right)\|_{hP}+\widetilde{\mathrm{BT}}_{h},\\\
&\widetilde{\mathrm{BT}}_{h}=\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{-},\widehat{\widetilde{T}}_{\eta}^{-}\right)+\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{+},\widehat{\widetilde{T}}_{\eta}^{+}\right)+\mathscr{IT}_{s}\left(\widehat{\widetilde{v}}^{\pm},\widehat{\widetilde{T}}^{\pm}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{-},Z^{-}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{+},Z^{+}\right)\leq
0,\end{split}$

where the surface terms $\widetilde{\mathrm{BT}}_{h}$ are negative semi-
definite.

##### The discrete PML edge problem.

Consider now the discrete PML (101) in the $xy$-edge region. In (101), let
$d_{x}=d_{y}=d\geq 0$, $d_{z}\equiv 0$,
$\mathbf{D}_{z}=\widetilde{\mathbf{FL}}_{z}=0$.

$\begin{split}\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=&\sum_{\xi=x,y}\left[\frac{1}{2}\left(\mathbf{A}_{\xi}\mathbf{D}_{\xi}-\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\mathbf{D}_{\xi}^{T}\mathbf{H}_{\xi}+\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\left(\mathbf{B}_{\xi}\left(1,1\right)\right)-\mathbf{B}_{\xi}\left(-1,-1\right)\right)\widetilde{\mathbf{U}}\right]\\\
-&\sum_{\xi=x,y}\left[\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}\right]+\mathbf{P}^{-1}\widetilde{\mathbf{F}}\end{split}$
(108)

We formulate the discrete equivalence of Theorem 6.

###### Theorem 12

Consider the semi-discrete PML equation in the Laplace space (108), with
$d_{x}=d_{y}=d\geq 0$, $\alpha_{x}=\alpha_{y}=\alpha\geq 0$,
$\gamma_{x}=\gamma_{y}=\gamma>0$ and $\real(s)\geq a>0$. We have
$\real(sS_{x})\geq\gamma a>0$ and

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{-}\left(s\right)\|^{2}_{hP}+\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{+}\left(s\right)\|^{2}_{hP}\leq\|\widetilde{\mathbf{U}}^{-}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{-}\left(s\right)\|_{hP}+\|\widetilde{\mathbf{U}}^{+}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{+}\left(s\right)\|_{hP}+\widetilde{\mathrm{BT}}_{h},\\\
&\widetilde{\mathrm{BT}}_{h}=\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{-},\widehat{\widetilde{T}}_{\eta}^{-}\right)+\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{+},\widehat{\widetilde{T}}_{\eta}^{+}\right)+\mathscr{IT}_{s}\left(\widehat{\widetilde{v}}^{\pm},\widehat{\widetilde{T}}^{\pm}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{-},Z^{-}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{+},Z^{+}\right)\leq
0,\end{split}$

where the surface terms $\widetilde{\mathrm{BT}}_{h}$ are negative semi-
definite.

##### The discrete PML corner problem.

Consider now a the discrete PML (101) in the corner region, where all damping
functions are nonzero, $d_{\xi}=d\geq 0$ for all $\xi=x,y,z$. For this case
the PML complex metrics are identical $S_{y}=S_{x}=S_{z}$. Thus
$S_{x}/S_{\xi}=1$, and from (101) we have

$\begin{split}\mathbf{P}^{-1}\left(sS_{x}\right)\widetilde{\mathbf{U}}=&\sum_{\xi=x,y,z}\left[\frac{1}{2}\left(\mathbf{A}_{\xi}\mathbf{D}_{\xi}-\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\mathbf{D}_{\xi}^{T}\mathbf{H}_{\xi}+\mathbf{H}_{\xi}^{-1}\mathbf{A}_{\xi}\left(\mathbf{B}_{\xi}\left(1,1\right)\right)-\mathbf{B}_{\xi}\left(-1,-1\right)\right)\widetilde{\mathbf{U}}\right]\\\
-&\sum_{\xi=x,y,z}\left[\mathbf{H}_{\xi}^{-1}{\left({\mathbf{e}_{\xi}(-1)}\widetilde{\mathbf{FL}}_{\xi}+{\mathbf{e}_{\xi}(1)}\widetilde{\mathbf{FR}}_{\xi}\right)}\right]+\mathbf{P}^{-1}\widetilde{\mathbf{F}}.\end{split}$
(109)

We will now state the discrete equivalence of Theorem 7.

###### Theorem 13

Consider the semi-discrete PML equation in the Laplace space (109), with
$d_{\xi}=d\geq 0$, $\alpha_{\xi}=\alpha\geq 0$, $\gamma_{\xi}=\gamma>0$ and
$\real(s)\geq a>0$. We have $\real(sS_{x})\geq\gamma a>0$ and

$\begin{split}&\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{-}\left(s\right)\|^{2}_{hP}+\|\sqrt{\real(sS_{x})}\widetilde{\mathbf{U}}^{+}\left(s\right)\|^{2}_{hP}\leq\|\widetilde{\mathbf{U}}^{-}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{-}\left(s\right)\|_{hP}+\|\widetilde{\mathbf{U}}^{+}\left(s\right)\|_{hP}\|\widetilde{\mathbf{F}}^{+}\left(s\right)\|_{hP}+\widetilde{\mathrm{BT}}_{h},\\\
&\widetilde{\mathrm{BT}}_{h}=\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{-},\widehat{\widetilde{T}}_{\eta}^{-}\right)+\mathscr{BT}_{s}\left(\widehat{\widetilde{v}}_{\eta}^{+},\widehat{\widetilde{T}}_{\eta}^{+}\right)+\mathscr{IT}_{s}\left(\widehat{\widetilde{v}}^{\pm},\widehat{\widetilde{T}}^{\pm}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{-},Z^{-}\right)+\mathscr{F}_{luc}\left({\widetilde{G}}^{+},Z^{+}\right)\leq
0,\end{split}$

where the surface terms are negative semi-definite.

A discrete equivalence of Theorem 8 also holds:

###### Corollary 2

Consider the semi-discrete energy estimate in the Laplace space

$\begin{split}&\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{hP}^{2}\leq\frac{1}{\gamma
a}\|\widetilde{\mathbf{U}}(\cdot,\cdot,\cdot,s)\|_{hP}\|\widetilde{\mathbf{F}}(\cdot,\cdot,\cdot,s)\|_{hP},\quad\real(s)\geq
a>0.\end{split}$

For any $a>0$, $\gamma>0$ and $T>0$ we have

$\displaystyle\int_{0}^{T}e^{-2at}\|{\mathbf{U}}(\cdot,\cdot,\cdot,t)\|_{hP}^{2}dt\leq\frac{1}{\left(\gamma
a\right)^{2}}\int_{0}^{T}e^{-2at}\left(\|{\mathbf{F}}(\cdot,\cdot,\cdot,t)\|_{hP}^{2}\right)dt.$

We can use Theorem 2 to invert the the estimates in Theorems 10, 11, 12 and 13
to get a discrete energy estimate in the physical space. As before the
estimate seemingly allows the energy of the solution to grow exponentially
with time for general data. However, by choosing $a>0$ in relation to the
length of the time interval of interest, a bound involving only algebraic
growth in time follows.

The following remarks are of important note.

###### Remark 13

For general first order systems the continuous energy estimates, Theorems 5–7,
and the discrete energy estimates, Theorems 11–13, are valid for all media
parameters including those that violate the so-called geometric stability
condition. These results may not eliminate instabilities that exist in the
continuous problems but they ensure robustness and guarantee longtime
numerical stability for the discrete PML, in particular for physical models
satisfying the geometric stability conditions.

###### Remark 14

We also remark that the discrete energy estimates, Theorems 11–13, are valid
for all discrete derivative operators that satisfy the SBP property (97).
These include, but not limited to, DG spectral difference operators defined on
Gaussian-type quadrature nodes and standard SBP finite difference operator on
equidistant grids.

###### Remark 15

For the general second order systems, (50)–(52) with the boundary condition
(25) and the interface conditions (54), such as displacement formulation of
linear elastodynamics, there is no general result yet to guide the derivation
of provably stable numerical methods for the PML. The development of such a
general procedure will be useful in designing reliable and efficient numerical
methods for second order systems, such as multi-block SBP finite difference
methods [51, 48, 89, 88], and energy-based and symmetric interior penalty DG
methods [58, 8, 9, 5, 31].

##  6 Numerical problems, What can go wrong? How it can be fixed?

We will perform some numerical experiments to review and demonstrate what can
go wrong with a numerical PML and show how it can be fixed. To do this we
consider the computational setup, Figure 3, a semi-infinite waveguide

Figure 3: A semi-infinite waveguide truncated by the PML of finite width
$\delta_{x}>0$.

We truncate the waveguide by the PML of finite width $\delta_{x}>0$, having
$\left(x,y\right)\in\Omega=[-L_{x},L_{x}+\delta_{x}]\times[0,L_{y}]$. The
domain $\Omega$ can be an electromagnetic medium, an acoustic medium or an
elastic medium, including isotropic and anisotropic elastic properties, inside
the waveguide. Throughout we will consider medium parameters for which the
geometric stability condition is satisfied and the unbounded constant
coefficient PML is stable.

We will consider DG discretisations but the results extend to SBP finite
difference methods. The discrete PMLs for different wave media exhibit very
similar numerical instabilities. As will be shown these numerical difficulties
can be cured by the theory developed in [40, 49, 36, 39], and summarised in
previous sections of this paper.

We consider a 2D problem with $\partial/\partial z=0$, $d_{z}(z)=0$, the
velocity component $v_{z}$ and the auxiliary variables $\mathbf{w}_{z}$ drop
out. For linear elasticity the stress field components
$\sigma_{xz},\sigma_{yz},\sigma_{zz}$ also drop out. We will also consider the
vertical PML strip problem (with $d_{x}(x)\geq 0$, $d_{y}(y)=0$), that is a
PML in the $x$-direction truncating the right boundary. In all experimental
setups, we will review the effectiveness and robustness of the PML stabilising
parameter $\theta_{x}=1$, as presented [40, 49, 36, 39].

The damping profile is a cubic monomial

$\begin{split}&d_{x}\left(x\right)=\left\\{\begin{array}[]{rl}0{}{}&\text{if}\quad
x\leq L_{x},\\\
d_{0}\Big{(}\frac{x-L_{x}}{\delta_{x}}\Big{)}^{3}&\text{if}\quad x\geq
L_{x},\end{array}\right.\end{split}$ (110)

where $d_{0}>0$ is the damping strength. We set the CFS $\alpha=0.15$, the
damping strength

$d_{0}=\frac{4c_{p}}{2\delta_{x}}\ln{\frac{1}{{tol}}},$ (111)

where $c_{p}$ denotes the P-wave speed, and ${tol}=10^{-3}$ is the magnitude
of the relative PML error [101]. We note that evanescent waves can be
problematic for standard PML models [24, 32, 54], unless the width of the
layer is expanded. Since the derivation of the formula (111) is based on
damping of propagating waves, and it must be considered as a heuristic formula
when evanescent waves are present. However, it has been demonstrated in [100]
that the inclusion of the complex frequency shift, CFS, in the PML metric will
enhance the absorption of evanescent waves. Also numerical simulations
presented in [40, 49, 39, 38, 98] suggest that the formula (111) is a good way
to choose the PML damping coefficient for practical simulations.

We will evolve the Gaussian initial wave profile

$\displaystyle f(x,y)=e^{-\log\left({2}\right)\frac{x^{2}+(y-Ly/2)^{2}}{9}},$
(112)

At the top and bottom boundaries, $y=0,L_{y}$, we set a free-surface boundary
condition with $r_{y}=1$ in (20)–(21), the left boundary at $x=-L_{x}$ is a
soft/clamped wall with $r_{x}=-1$ in (20)–(21) and the PML boundary at
$x=L_{x}+\delta_{x}$ is terminated with the classical first order ABC, by
setting $r_{x}=0$ in (20)–(21). When the PML damping vanishes, the set-up will
correspond to a semi-infinite waveguide terminated by the classical first
order ABC, and will reflect all waves with non-normal incident angle. The
first order ABC is exact for waves impinging at normal incidence on the
boundary. Specifically, for linear elasticity this is the first order Lysmer-
Kuhlemeyer ABC which is exact for plane shear and pressure waves with normal
incidence angles at the boundary [85, 11]. We will compare the accuracy of the
PML and the ABC.

We set the domain width $L_{\xi}=50$ km, the PML width $\delta_{x}=10$ km. We
discretise the domain with a uniform element size $\Delta{x}=\Delta{y}=5$ km,
spanning the PML with two DG elements, and approximate the solution by a
polynomial of degree $N=4$, with the sub-cell resolution $h=\Delta{x}/(N+1)=1$
km. Time-integration is performed using the high order ADER scheme [106, 34]
of the same order of accuracy with the spatial discretization. We use the time
step

$\displaystyle\Delta{t}=\frac{CFL}{\sqrt{d}\left(2N+1\right)c_{p}}\min{\left(\Delta{x},\Delta{y}\right)},$
(113)

with the $CFL=0.9$ number and the P-wave speed $c_{p}$ [km/s] and $d=2$ is the
spatial dimension. The final time is $t=200$ s.

In order to access numerical PML errors we compute a reference solution in a
larger domain $\left(x,y\right)\in\Omega_{ref}=[-L_{x},3L_{x}]\times[0,L_{y}]$
and for a limited time such that reflection from the reference domain do not
re-enter the interior computational domain. By comparing the reference
solution and the PML solution in the interior,
$\left(x,y\right)\in[-L_{x},L_{x}]\times[0,L_{y}]$, we get an accurate measure
of the PML error.

###  6.1 Acoustic waveguide

We consider first a 2D acoustic medium with constant acoustic wave speed
$c_{p}=1.484$ km/s and constant medium density $\rho=1$ g/cm3. We set the
initial condition (112), $p(x,y,t=0)=f(x,y)$, for the pressure field, and zero
initial condition for the velocity fields and the auxiliary variables. The
snapshots of the absolute pressure fields are plotted in Figure 4, for
$\theta_{x}=0,1$, and the time history L∞-norm of the pressure field plotted
in Figure 5(a).

[5pt][5pt][5pt]$\theta_{x}=0$.[5pt]$\theta_{x}=1$.

Figure 4: Absolute pressure field $|p(x,y,t)|$ at $t=25,40,200$ s for
$\theta_{x}=0,1$. The penalty parameter $\theta_{x}=1$ stabilises the PML
solution at long times, $t=200$ s.

(a) L∞ norm of $p$ for the PML with $\theta_{x}=0,1$

(b) PML error versus ABC error for $p$

Figure 5: Time series of the L∞ norm of the pressure field for the PML with
$\theta_{x}=0,1$.

We can see from Figure 4 that without the stabilisation parameter, that is
with $\theta_{x}=0$, the solution in the PML explodes after some time steps.
See also Figure 5(a). As shown in [39], the initiation time of the explosive
numerical mode depends on the quadrature rule used and the mesh resolution. On
a finer mesh the growth persists, but it starts at a much later time. The PML
stabilisation parameter $\theta_{x}=1$ guarantees both accuracy and numerical
stability of the PML at long times. See also Figure 5(a) and 5(b). While the
ABC is numerically stable, it however generates unacceptable large spurious
reflections errors. The errors can be reduced by moving the artificial
boundary further away. This however will lead to a significant increase in the
computational cost.

These experiments demonstrate the significant importance of proving continuous
and discrete energy estimates for the PML in the Laplace space, and its impact
in ensuring long-time stability of the PML. The results extend to 2D Maxwell’s
equations of electrodynamics and to SBP-SAT finite difference method [36]. The
initial ideas of proving continuous and discrete energy estimates for the PML
in the Laplace space were first derived for the 2D Maxwell’s equation using
the SBP-SAT finite difference method [36]. In [39], these were recently
extended to acoustic waves, and to 3D and DG methods. Although the proofs are
derived for constant and piece-wise damping coefficients, however, numerical
experiments demonstrate that the analysis may be valid for continuously
varying damping profiles.

###  6.2 Isotropic and anisotropic elastic waveguides

Next we consider the PML truncating a linear elastic waveguide. Both isotropic
and anisotropic media properties will be considered. The density of the media
and elastic constants are given in Table 1. In both isotropic and anisotropic
elastic media we have chosen the media properties such that the maximum P-wave
speed is $c_{p}=6$ km/s. The dispersion relation of the isotropic and
anisotropic elastic media considered here are shown in Figure 2(b)–2(c), and
indicate that the Cauchy PML problem is stable.

Table 1: Elastic media parameters. Parameters | Isotropic medium | Anisotropic medium (AM1)
---|---|---
$c_{11}$ [GPa] | 97.20 | 20
$c_{12}$ [GPa] | 36.85 | 3.8
$c_{22}$ [GPa] | 97.20 | 4
$c_{33}$ [GPa] | 30.17 | 2
$\rho$ [gm/cm3] | $2.7$ | $20/36$

In elastic media, we set the initial condition (112),
$v_{x}(x,y,t=0)=v_{y}(x,y,t=0)=f(x,y)$, for the velocity field, and zero
initial conditions for the stress fields and the auxiliary variables.

The snapshots of the absolute velocity are plotted in Figure 6, for
$\theta_{x}=0,1$, and the time history L∞-norm of the absolute velocity and
the numerical PML errors are plotted in Figure 7.

[5pt][5pt][5pt]$\theta_{x}=0$.[5pt]$\theta_{x}=1$.

(a) Isotropic medium

[5pt][5pt][5pt]$\theta_{x}=0$.[5pt]$\theta_{x}=1$.

(b) Anisotropic medium

Figure 6: Snapshots of the absolute particle velocity
$\sqrt{v_{x}^{2}+v_{y}^{2}}$ in 2D elastic solids.

(a) Isotropic medium: L∞ norm of $\mathbf{v}$ for the PML.

(b) Isotropic medium: PML error versus ABC error.

(c) Anisotropic medium: L∞ norm of $\mathbf{v}$ for the PML.

(d) Anisotropic medium: PML error versus ABC error.