<figure><img src="image_1.png" /></figure>

HAL
open science

# A time domain analysis of PML models in acoustics.

Julien Diaz, Patrick Joly

## ► To cite this version:

Julien Diaz, Patrick Joly. A time domain analysis of PML models in acoustics. Computer Methods in Applied Mechanics and Engineering, 2006, 195 (29-32), pp.3820-3853. 10.1016/j.cma.2005.02.031. inria-00410313

HAL Id: inria-00410313

https://inria.hal.science/inria-00410313

Submitted on 19 Aug 2009

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A time domain analysis of PML models in
acoustics.

J. Diaz, P. Joly

INRIA Rocquencourt - Domaine de Voluceau - 78153 Le Chesnay - France.

Projet POEMS - UMR 2706 - CNRS - ENSTA - INRIA

Abstract

We present in this paper a time-domain analysis of PML's for non-adptive and
adptive acoustics. We focus our attention on time-stability and error estimates
(with respect to the parameters of the layers). The main new technical tool is the
Cagniard-de Hoop method. Our theoretical results are validated and illustrated by
various numerical results.

Key words: Perfectly Matched Layers; Acoustics; Stability; Error Estimates; Cagniard-de Hoop.

PACS:

The PML technique for the numerical absorption of waves, initially introduced
about ten years ago by Bérenger [1] in electromagnetism, is now widely used
for simulating the propagation of waves in unbounded domains, in particular
in time domain acoustics [2–8].

From the mathematical point of view, much work has been devoted to the
stability analysis of the PML model [9,10]. The question of the accuracy of
PML's in time domain is much less discussed. Existing results mainly con-
cern the plane wave analysis [4,11] (reflection of plane waves at the artificial
boundary) or convergence results in the time harmonic regime [12–14].

The object of the present paper is to regroup some results obtained by the two
authors concerning the time domain analysis for acoustic wave propagation
models.

Email addresses: julien.diaz@inria.fr (J. Diaz), patrick.joly@inria.fr
(P. Joly).
In section 1, we first derive an error analysis for the PML for non advective acoustics. This work is a continuation of [15] that concerned absorbing boundary conditions. The main point is the derivation of an explicit solution using the Cagniard-De Hoop technique [16–18]. Our results show the exponential accuracy of PML’s with respect to the damping factor $\sigma$ and the width $L$ of the layer.

Section 2 is devoted to classical PML’s for advective acoustics: the propagation of waves in the presence of a uniform mean flow is modeled by linearized Euler equations. We revisit the well-known instability of standard PML’s by the Cagniard-De Hoop technique.

In section 3, we are concerned by the question of the stabilization of PML’s for advective acoustics when the flow is orthogonal to the artificial boundary. We proposed a stable PML (already introduced in [19]) which consists, in a first step, of sticking to the physical domain a non advective propagation model designed in such a way that no reflection is produced at the interface between the physical domain and the artificial one (this is the ”PML property”), and in a second step, of using standard PML’s for the artificial domain. We derive the stability and accuracy analysis for this new PML.

Our theoretical results are illustrated by various numerical results. Of course the present work has to be compared with some recent ones in the literature, in particular by Hagström and Nazarov [20], Abarbanel, Gottlieb and Hesthaven [21] and more specially by Hu [22] whose work on stabilized PML for advective acoustics appears to be close to what we are doing in this paper as it will be briefly discussed in the last section of this paper.

An appendix is devoted to a partial solution to the problem of constructing stable PML’s in the case of an oblique (with respect to the boundaries of the rectangular computational domain) mean flow.

# 1 PML for acoustics

Let us consider the two dimensional wave equation

$$ \frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - \Delta u = f(x, t), \quad x = (x_1, x_2) \in \mathbb{R}^2, \quad t > 0. \quad (1) $$

If the initial data of the problem and the domain of interest are supported in the left half-space

$$ \mathbb{R}^2 = \mathbb{R}_- \times \mathbb{R} $$

it is natural to try to reduce the computation to this half-space. There are two main classes of method to do so. The first ones are the so-called absorbing
boundary condition, whose principle consists in imposing a specific condition on the boundary designed so that the waves leave the computational domain without being (too much) reflected. Another way to bound the computational domain consists in adding a vertical layer where (1) is modified by an absorption term (see figure 1). In [1,6], Bérenger proposed to introduce a very special

<figure><img src="image_2.png" /><figcaption>Fig. 1. Computational domain and Perfectly Matched Layer</figcaption></figure>

anisotropic damping term acting only in the direction which is orthogonal to the interface. An alternative to the original split form of its model (see section 2) consists in replacing the derivation with respect to $x_1$ by the operator $D_{x_1}^\sigma$ formally defined by :

$$D_{x_1}^\sigma = \left( \frac{\partial}{\partial t} + \sigma(x_1) \right)^{-1} \frac{\partial}{\partial t} \frac{\partial}{\partial x_1},$$

where $\sigma$ is a non-negative function of $x_1$.

More rigorously $\psi = D_{x_1}^\sigma \phi$ means :

$$\left(\frac{\partial}{\partial t} + \sigma(x_1)\right) \psi = \frac{\partial}{\partial t} \frac{\partial}{\partial x_1} \phi \quad \text{and} \quad \psi(t=0) = 0.$$

The equation in the PML then becomes :

$$\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - (D_{x_1}^\sigma)^2 u - \frac{\partial^2 u}{\partial x_2^2} = 0, \quad \mathbf{x} = (x_1, x_2) \in [0; L] \times \mathbb{R}, \quad t > 0. \quad (2)$$

The most interesting property of the PML is that a wave propagating in the computational domain is transmitted to the absorbing layer without being reflected. As this method is moreover really easy to implement (even in the corners of the domain), it has rapidly attracted a lot of people in different fields of application.
## 1.1 The Main results

**The case of an infinite PML.** We are first interested in the explicit expression of the fundamental solution $u^{\sigma, \infty}$ of the wave equation in the left half space $\mathbb{R}^2_-$ with an infinite PML in the right half-space. $u^{\sigma, \infty}$ satisfies :

$$ egin{cases} 
\text{Find } u : \mathbb{R}^2 \times \mathbb{R}_+ \mapsto \mathbb{R}, \text{ null for } t < 0 & \text{such that} \\[1em] 
\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - (D_{x_1}^{\sigma})^2 u - \frac{\partial^2 u}{\partial x_2^2} = \delta(\mathbf{x} - \mathbf{x}_S) \times \delta(t), & \text{in } \mathbb{R}^2 \times \mathbb{R}_+ 
\end{cases} 
\quad (3) $$

where $\mathbf{x}_S = (-h, 0)$, $\sigma(x_1) = 0$ if $x_1 \le 0$ and $\sigma(x_1) > 0$ if $x_1 > 0$.

Let us define $r(\mathbf{x}) = |\mathbf{x} - \mathbf{x}_S|$ and the function $\theta(\mathbf{x})$ by :

$$ \theta(\mathbf{x}) \in [0, 2\pi], \quad \mathbf{x} - \mathbf{x}_S = (r(\mathbf{x}) \cos \theta(\mathbf{x}), r(\mathbf{x}) \sin \theta(\mathbf{x}))^t, $$

and the function $\Sigma(x_1)$ by : $\Sigma(x_1) = \int_0^{x_1} \sigma(x) dx$. The functions $A(\mathbf{x}, t)$ and $B(\mathbf{x}, t)$ are defined by ($r = r(\mathbf{x})$ and $\theta = \theta(\mathbf{x})$) :

$$ A(\mathbf{x}, t) = |\cos \theta| \frac{\Sigma(x_1) t}{r} > 0 \quad \text{and} \quad B(\mathbf{x}, t) = |\sin \theta| \Sigma(x_1) \sqrt{\frac{t^2}{r^2} - \frac{1}{c^2}}. $$

**Theorem 1** The expression of the solution $u^{\sigma, \infty}(\mathbf{x}, t) = G_i^{\sigma, \infty}(\mathbf{x}, t)$ of problem (3) is given by :

$$ egin{cases} 
G_i^{\sigma, \infty}(\mathbf{x}, t) = \frac{H(c t - r(\mathbf{x}))}{2\pi\sqrt{t^2 - \frac{r(\mathbf{x})^2}{c^2}}}, & x_1 \le 0 \\[1em] 
G_i^{\sigma, \infty}(\mathbf{x}, t) = \frac{H(c t - r(\mathbf{x}))}{2\pi\sqrt{t^2 - \frac{r(\mathbf{x})^2}{c^2}}} e^{-A(\mathbf{x}, t)} \cos[B(\mathbf{x}, t)], & x_1 > 0. 
\end{cases} 
\quad (4) $$

where $H$ denotes the Heaviside function

**Remark 2** The restriction of $G_i^{\sigma, \infty}(\mathbf{x}, t)$ to the left half-space $\mathbb{R}^2_-$ does not depend on $\sigma$ and is nothing more than the restriction to the left half-space of the fundamental solution of the wave equation in the whole plane : the infinite PML does not produce any parasitical reflection.

**The case of a finite PML.** We now consider a PML of finite width $L$ with
a Neumann condition at its bottom end. Problem (3) then becomes :

$$
\left\{
\begin{array}{ll}
\text{Find } u : \mathbb{R}^2 \times \mathbb{R}_+ \to \mathbb{R}, \text{ null for } t < 0 & \text{such that} \\[1em]
\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - (D_{x_1}^\sigma)^2 u - \frac{\partial^2 u}{\partial x_2^2} = \delta(\mathbf{x} - \mathbf{x}_S) \times \delta(t), & \text{in } [-\infty; L] \times \mathbb{R} \times \mathbb{R}_+, \\[1em]
\frac{\partial u}{\partial x_1}(L, x_2, t) = 0.
\end{array}
\right.
\tag{5}
$$

We define the image point source, symmetrical to the point source with respect to the line $x_1 = L$ by :

$$
\boldsymbol{x}_S^* = (h + 2L, 0),
$$

we denote $r^*(\boldsymbol{x}) = |\boldsymbol{x} - \boldsymbol{x}_S^*|$ and we define the function $\theta^*(\boldsymbol{x})$ by :

$$
\theta^*(\boldsymbol{x}) \in ]\frac{\pi}{2}, \frac{3\pi}{2}] , \quad \boldsymbol{x} - \boldsymbol{x}_S^* = ( r^*(\boldsymbol{x}) \cos \theta^*(\boldsymbol{x}), r^*(\boldsymbol{x}) \sin \theta^*(\boldsymbol{x}) )^t.
$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of the notation</figcaption></figure>

We also define the function $\Sigma^*(x_1)$ by :

$$
\Sigma^*(x_1) = \int_0^L \sigma(x) \, dx + \int_{x_1}^L \sigma(x) \, dx \quad (= 2\Sigma(L) \text{ for } x_1 < 0),
$$

and the functions $A^*(\mathbf{x}, t)$ and $B^*(\mathbf{x}, t)$ by ($r^* = r^*(\mathbf{x})$ and $\theta^* = \theta^*(\mathbf{x})$):

$$
A^*(\mathbf{x}, t) = |\cos \theta^*| \Sigma^*(x_1) \frac{t}{r^*} > 0 \quad \text{and} \quad B^*(\mathbf{x}, t) = |\sin \theta^*| \Sigma^*(x_1) \sqrt{\frac{t^2}{r^{*2}} - \frac{1}{c^2}}.
$$

**Theorem 3** The expression of the solution $u^{\sigma,L}(\mathbf{x},t) = G^{\sigma,L}(\mathbf{x},t)$ of problem (5) is given by :

$$
G^{\sigma,L}(\mathbf{x}, t) = G_{i}^{\sigma,\infty}(\mathbf{x}, t) + G_{r}^{\sigma,L}(\mathbf{x}, t)
$$

with

$$
G_r^{\sigma,L}(\mathbf{x},t) = \frac{H(c t - r^*(\mathbf{x}))}{2\pi\sqrt{t^2 - \frac{r^*(\mathbf{x})^2}{c^2}}} e^{-A^*(\mathbf{x},t)} \cos[B^*(\mathbf{x},t)]. \quad (6)
$$
This theorem can easily be deduced from theorem 1 thanks to the image principle. Therefore, only theorem 1 will be proved in this paper.

**Error Estimates.** Let us consider the approximation, in the upper half-space, of the solution $u$ of problem (1) (the wave equation with a “regular” point source), by the function $u^{\sigma,L}$, solution of problem :

$$ egin{cases} 
\text{Find } u^{\sigma,L}: \mathbb{R}^2 \times \mathbb{R}_+ \mapsto \mathbb{R}, \text{ null for } t < 0 & \text{such that} \\ 
\frac{1}{c^2} \frac{\partial^2 u^{\sigma,L}}{\partial t^2} - (D_{x_1}^{\sigma})^2 u^{\sigma,L} - \frac{\partial^2 u^{\sigma,L}}{\partial x_2^2} = \delta(x - x_S) \times f(t), \text{ in } [-\infty; L] \times \mathbb{R} \times \mathbb{R}_+, \\ 
\frac{\partial u^{\sigma,L}}{\partial x_1}(L, x_2, t) = 0. 
\end{cases} 
\tag{7} $$

After having defined the function $\Phi(t)$, null for $t < (h + 2L)/c$, by :

$$ egin{cases} 
\Phi(t) = \text{Log} \left( \frac{t + \sqrt{t^2 - ((2L + h)/c)^2}}{((2L + h)/c)} \right), & \text{if } \frac{h + 2L}{c} \le t \le \frac{h + 2L}{c} + T, \\ 
\Phi(t) = \text{Log} \left( \frac{t + \sqrt{t^2 - (t - T)^2}}{t - T} \right), & \text{if } t > \frac{h + 2L}{c} + T. 
\end{cases} $$

we can state the

**Theorem 4** One has the uniform estimates :

$$ \|(u - u^{\sigma,L})(\cdot, t)\|_{L^{\infty}(\mathbb{R}^2)} \le \frac{e^{-2\frac{2L+h}{c^2t}L\bar{\sigma}}}{2\pi} \Phi(t) \|f\|_{L^{\infty}}. \tag{8} $$

where $\bar{\sigma} = \frac{1}{L} \int_0^L \sigma(x)dx$ is the mean value of $\sigma$ in the PML.

This result appeals the following comments :

* The error converges spectrally to 0 (in the uniform norm) when $\bar{\sigma}$ (or $L$) goes to infinity.

* For given $t$, the upper bound in the estimate diminishes when the distance $h$ from the source to the absorbing layer increases. This is coherent with the physical intuition and numerical observations.

* Concerning the behavior of the error for large $t$, if we assume that $T < +\infty$, we observe that the right hand side in the estimate behaves for large $t$ as

$$ \frac{1}{2\pi} \sqrt{\frac{2T}{t}}, $$

which shows that, for all $\bar{\sigma}$, $L$ and $h$, the error converges uniformly to 0 when $t$ tends to $+\infty$.
* As the plane wave analysis might suggest, the shape of the function $\sigma(x_1)$ does not influence the estimates, only its mean value between 0 and L has an importance.

* We could have imagined, once again from the plane wave analysis, that $\bar{\sigma}$ and L did play the same role (i.e. increasing L is equivalent to increasing $\bar{\sigma}$), but the uniform estimates show that increasing L allows us, not only to increase the absorption thanks to the term $L \bar{\sigma}$, but also to move away the point source from the absorbing layer thanks to the term $2L + h$.

**The case of a rectangular domain surrounded by PML.** One of the advantages of PML when compared to the absorbing boundary conditions (ABC's) is the simplicity of their implementation in the corner. By the principle of images and Cagniard-de Hoop method, it is also easy to compute the analytical solution in the case of a rectangular domain surrounded by PML's (whereas it is not possible with ABC's).

We now suppose that the computational domain $\Omega$ is the square $[-h; h]^2$ and is surrounded by PML's of width L. We impose a Neumann condition at the end of the layers. We then consider the problem :

$$ 
\left\{
\begin{array}{ll}
\text{Find } u : \mathbb{R}^2 \times \mathbb{R}_+ \mapsto \mathbb{R} & \text{such that} \\
\\
\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} - (D_{x_1}^\sigma)^2 u - (D_{x_2}^\sigma)^2 u = \delta(\mathbf{x}) \times \delta(t), \quad \text{in } [-L; L]^2 \times \mathbb{R}_+, \\
\\
\frac{\partial u}{\partial x_1}(\mathbf{x}, t) = 0, & \text{for } x_1 = \pm L \\
\\
\frac{\partial u}{\partial x_2}(\mathbf{x}, t) = 0, & \text{for } x_2 = \pm L, \\
\\
u(\mathbf{x}, t) = 0, & \text{for } t < 0.
\end{array}
\right.
\quad (9) $$

Let us define the image points source $(\mathbf{x}_{ij})_{(i,j) \in \mathbb{Z}^2}$ as

$$ \mathbf{x}_{ij} = (2i(h + L), 2j(h + L)). $$

We set $r_{ij}(\mathbf{x}) = |\mathbf{x} - \mathbf{x}_{ij}|$ and we define the functions $\theta_{ij}(\mathbf{x})$ by :

$$ \theta_{ij}(\mathbf{x}) \in ] -\pi, \pi], \quad \mathbf{x} = (r_{ij}(\mathbf{x}) \sin \theta_{ij}(\mathbf{x}), r_{ij}(\mathbf{x}) \cos \theta_{ij}(\mathbf{x}))^t, $$
and the non-negative functions $(\Sigma_i(x))_{i \in \mathbb{Z}}$ by :

$$ \left\{ \begin{array}{ll} \Sigma_i(x) = (2i-1) \int_0^L \sigma(x) \, dx + \int_{-L}^x \sigma(x) \, dx & \text{if } i < 0, \\ \\ \Sigma_i(x) = (2i-1) \int_0^L \sigma(x) \, dx + \int_x^L \sigma(x) \, dx & \text{if } i > 0, \\ \\ \Sigma_0(x) = \left| \int_0^x \sigma(x) \, dx \right| & \end{array} \right. $$

Finally we define the functions $A_{ij}(\mathbf{x}, t)$ and $B_{ij}(\mathbf{x}, t)$ by :

$$ A_{ij}(\mathbf{x}, t) = (|\cos \theta_{ij}(\mathbf{x})| \Sigma_j(x_2) + |\sin \theta_{ij}(\mathbf{x})| \Sigma_i(x_1)) \frac{t}{r_{ij}(\mathbf{x})} > 0 $$

and

$$ B_{ij}(\mathbf{x}, t) = (|\cos \theta_{ij}(\mathbf{x})| \Sigma_i(x_1) + |\sin \theta_{ij}(\mathbf{x})| \Sigma_j(x_2)) \sqrt{\frac{t^2}{r_{ij}(\mathbf{x})^2} - \frac{1}{c^2}} $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. The image points source</figcaption></figure>

**Theorem 1.1** The solution $u^{\sigma,h,L}(x,t) = G^{\sigma,h,L}(x,t)$ of problem (9) is given by :

$$ G^{\sigma,h,L}(\mathbf{x}, t) = \sum_{(i,j) \in \mathbb{Z}^2} G_{i,j}^{\sigma,h,L}(\mathbf{x}, t) $$
where

$$G_{ij}^{\sigma,h,L}(\boldsymbol{x}, t) = \frac{H(c t - r_{ij}(\boldsymbol{x}))}{2\pi\sqrt{t^2 - \frac{r_{ij}(\boldsymbol{x})^2}{c^2}}} e^{-A_{ij}(\boldsymbol{x},t)} \cos[B_{ij}(\boldsymbol{x},t)]. \quad (10)$$

**Remark 5** For given $t$ and $\boldsymbol{x}$, $G^{\sigma,h,L}(\boldsymbol{x}, t)$ is a finite sum. Indeed if $i$ and $j$ are such that $r_{ij}(\boldsymbol{x}) < ct$, then $G_{ij}^{\sigma,h,L}(\boldsymbol{x}, t)$ is null.

We do not present here the details of the proof which is very similar to the half-plane case, we refer the reader to [23] for details.

## 1.2 Proof of the theorems

**Expression of the Green's function.** Let $G_i^{\sigma,\infty}(\boldsymbol{x}, t)$ be the solution of problem (3). The first step of the Cagniard-de Hoop's method consists in applying successively to $G_i^{\sigma,\infty}(\boldsymbol{x}, t)$ the Laplace transform in time and the Fourier transform in the space variable $x_2$ (we denote by $s$ (resp. $k$) the dual variable of $t$ (resp. $x_2$)). The resulting function $x_1 \mapsto \hat{G}_i^{\sigma,\infty}(x_1, k, s)$ satisfies the ordinary differential equation :

$$-\frac{s}{s + \sigma(x_1)} \frac{d}{dx_1} \left( \frac{s}{s + \sigma(x_1)} \frac{d\hat{G}_i^{\sigma,\infty}}{dx_1} \right) + \left(k^2 + \frac{s^2}{c^2}\right) \hat{G}_i^{\sigma,\infty} = \delta(x_1 + h),$$

Introducing $\hat{G}(X_1(x_1, s) = x_1 + \frac{1}{s}\int_0^{x_1} \sigma(x)dx, k, s) = \hat{G}_i^{\sigma,\infty}(x_1, k, s)$ we get :

$$-\frac{d^2 \hat{G}}{dX_1^2} + \left(k^2 + \frac{s^2}{c^2}\right) \hat{G} = \delta(X_1 + h).$$

which leads to the expressions :

$$\hat{G}(X_1, k, s) = \frac{e^{-(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}} |X_1 + h|}}{2(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}} \implies \hat{G}_i^{\sigma,\infty}(x_1, k, s) = \frac{e^{-(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}} |x_1 + h + \frac{\Sigma(x_1)}{s}|}}{2(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}}.$$

Let us remark that, if $x_1$ is negative, $\Sigma(x_1) = 0$ and

$$\hat{G}_i^{\sigma,\infty}(x_1, k, s) = \frac{e^{-(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}} |x_1 + h|}}{2(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}},$$

which is nothing but the Fourier-Laplace transform of the fundamental solution of the wave equation :

$$\hat{G}_i^{\sigma,\infty}(x_1, x_2, t) = \frac{H(c t - r(\boldsymbol{x}))}{2\pi\sqrt{t^2 - \frac{r(\boldsymbol{x})^2}{c^2}}}, \quad x_1 < 0.$$
Let us now consider $x_1 > 0$ so that $x_1 + h + \frac{\Sigma(x_1)}{s} > 0$ and

$$ \hat{G}_{i}^{\sigma, \infty}(x_1, k, s) = \frac{e^{-(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}(x_1 + h + \frac{\Sigma(x_1)}{s})}}{2(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}}. $$

Applying the $x_2$ inverse Fourier transform to $\hat{G}_{i}^{\sigma, \infty}$ we get:

$$ \tilde{G}_{i}^{\sigma, \infty}(x_1, x_2, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{e^{-(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}(x_1 + h + \frac{\Sigma(x_1)}{s}) - ikx_2}}{(k^2 + \frac{s^2}{c^2})^{\frac{1}{2}}} dk. $$

Contrary to what a reader familiar with the Cagniard-de Hoop could expect, the term in the exponential is not an homogeneous function in $(k, s)$. Applying the usual change of variable $k = ps/c$ allows us however to obtain:

$$ \tilde{G}_{i}^{\sigma, \infty}(x_1, x_2, s) = \int \frac{e^{-\frac{\Sigma(x_1)}{c}(1+p^2)^{\frac{1}{2}}} e^{-\frac{s}{c}[(1+p^2)^{\frac{1}{2}}(x_1+h)+ipx_2]}}{4\pi(1+p^2)^{\frac{1}{2}}} dp, \quad (\equiv \int_{-\infty}^{+\infty} \Psi(p) dp). $$

We then have to deal with the product of two exponential functions. The first one does not depend on $s$ and the other one is the one that usually appears in most applications of the Cagniard-de Hoop method.

Let us introduce the so-called Cagniard-de Hoop contour $\Gamma$ defined as:

$$ \Gamma = \Gamma^{+} \cup \Gamma^{-} \text{ with } \Gamma^{\pm} = \left\{ p = \gamma^{\pm}(t) \equiv -i \frac{ct}{r} \sin\theta \pm \cos\theta \sqrt{\frac{c^2t^2}{r^2} - 1}, t \ge \frac{r}{c} \right\}. $$

so that $(1 + p^2)^{\frac{1}{2}}(x_1 + h) + ipx_2 = ct \in \mathbb{R}_+$ for all $p \in \Gamma$.

It is clear that the two curves $\Gamma^{\pm}$ are symmetric the one from the other with respect to the imaginary axis and meet at point $-i \sin\theta$ (for $t = r/c$). Actually $\Gamma$ is nothing but the branch of the hyperbola of equation:

$$ \frac{Y^2}{\sin^2 \theta} - \frac{X^2}{\cos^2 \theta} = 1, \quad (p = X + iY, (X, Y) \in \mathbb{R}^2), $$

which is located in the upper half-space $Y = \mathfrak{Im} p > 0$ if $\sin\theta < 0$. Note that this hyperbola does not intersect the two branch cuts of $\Psi$. All this information is summarized in figure 4.

Let us denote by $D$ the real line and by $\Omega$ the connected part of the complex plane delimited by $D$ and $\Gamma$. Let $\rho > 0$ a parameter devoted to tend to $+\infty$. We set:

$$ D_{\rho} = \{p \in D / |p| \le \rho\}, \quad \Gamma_{\rho} = \{p \in \Gamma / |p| \le \rho\} \text{ and } C_{\rho} = \{p \in \Omega / |p| = \rho\}. $$

Note that $D_{\rho} \cup C_{\rho} \cup \Gamma_{\rho}$ is a closed curve (see figure 5). Since $\Psi(p)$ is analytic
<figure><img src="image_5.png" /><figcaption>Fig. 4. The contours $\Gamma$ and $D$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The closed contour $D_{\rho} \cup C_{\rho} \cup \Gamma_{\rho}$</figcaption></figure>

in $\Omega$, the integral of $\Psi$ along $D_{\rho} \cup C_{\rho} \cup \Gamma_{\rho}$, that we orient according to figure 5, is identically 0:

$$ \int_{D_{\rho}} \Psi(p) dp + \int_{C_{\rho}} \Psi(p) dp + \int_{\Gamma_{\rho}} \Psi(p) dp = 0. $$

Thanks to the choice of the square root and since $x_1 + h > 0$, the function $\Psi(p)$ decays exponentially to 0 when $\Im p$ goes to $+\infty$. As a consequence, it is easy to show that (Jordan's lemma):

$$ \lim_{\rho \to +\infty} \int_{C_{\rho}} \Psi(p) dp = 0, \text{ which implies } \int_{-\infty}^{+\infty} \Psi(p) dp = - \int_{\Gamma} \Psi(p) dp \quad \text{and} $$

$$ \tilde{G}_{i}^{\sigma, \infty}(x_1, x_2, s) = -\frac{1}{4\pi} \int_{\Gamma} \frac{e^{-\frac{\Sigma(x_1)}{c}(1+p^2)^{\frac{1}{2}}} e^{-s \left[ (1+p^2)^{\frac{1}{2}} \frac{x_1+h}{c} + ip \frac{x_2}{c} \right]}}{(1+p^2)^{\frac{1}{2}}} dp. $$

We use the parametrisations $p = \gamma^{+}(t)$ and $p = \gamma^{-}(t)$, for $t \ge r/c$, respectively along $\Gamma^{+}$ and $\Gamma^{-}$ and remark that:

$$ \bullet (1 + p^2)^{\frac{1}{2}} \left( \frac{x_1 + h}{c} \right) + ip \frac{x_2}{c} = t, \quad \text{(by construction)}, $$

$$ \bullet \frac{dp}{(1 + p^2)^{\frac{1}{2}}} = \pm \frac{dt}{(t^2 - \frac{r^2}{c^2})^{\frac{1}{2}}}, \quad \text{on } \Gamma^{\pm}, $$

$$ \bullet \frac{(1 + p^2)^{\frac{1}{2}}}{c} = \frac{t}{r^2}(x_1 + h) \mp i \frac{x_2}{r^2} \left(t^2 - \frac{r^2}{c^2}\right)^{\frac{1}{2}}, \quad \text{on } \Gamma^{\pm}. $$

Therefore, since $t$ goes from $+\infty$ to $\frac{r}{c}$ on $\Gamma^{+}$ and from $\frac{r}{c}$ to $+\infty$ on $\Gamma^{-}$:

$$ \tilde{G}_{i}^{\sigma, \infty}(x_1, x_2, s) = \frac{1}{2\pi} \int_{\frac{r}{c}}^{+\infty} \frac{e^{-(x_1+h)\Sigma(x_1)\frac{t}{r^2}}}{(t^2 - \frac{r^2}{c^2})^{\frac{1}{2}}} \cos \left[ \frac{x_2}{r^2} \Sigma(x_1) (t^2 - \frac{r^2}{c^2})^{\frac{1}{2}} \right] e^{-st} dt. $$

We can now conclude, using the injectivity of the Laplace-Fourier transform.

**Error Estimates** Let $u$ and $u^{\sigma,L}$ be the respective solutions of (1) and (7). We introduce the error (or reflected field) defined as $e^{\sigma,L} = u^{\sigma,L} - u$. Using the fact that $f$ is supported in $[0, T]$ and $G_{r}^{\sigma,L}$ in $[r^*/c, \infty]$, we obviously have
$e^{\sigma,L}(\boldsymbol{x},t) = 0$ if $t \le r^*/c$ and :

$$e^{\sigma,L}(\boldsymbol{x},t) = \int_{\max(r^*/c, t-T)}^{+\infty} G_r^{\sigma,L}(\boldsymbol{x}, \tau) f(t-\tau) d\tau \quad \text{if } t > \frac{r^*}{c}. \quad (11)$$

We deduce that :

$$\left\{ \begin{array}{l} |e^{\sigma,L}(\boldsymbol{x},t)| \le \|f\|_{L^{\infty}(0,t)} \cdot \|G_r^{\sigma,L}(\boldsymbol{x}, \cdot)\|_{L^1(\frac{x^*}{c},t)} \quad \text{if } \frac{r^*}{c} \le t \le \frac{r^*}{c} + T, \\ \\ |e^{\sigma,L}(\boldsymbol{x},t)| \le \|f\|_{L^{\infty}(0,T)} \cdot \|G_r^{\sigma,L}(\boldsymbol{x}, \cdot)\|_{L^1(t-T,t)} \quad \text{if } t > \frac{r^*}{c} + T. \end{array} \right. \quad (12)$$

Let us now introduce the two disjoint sets $\Omega_1(t)$ and $\Omega_2(t)$ :

$$\left| \begin{array}{l} \Omega_1(t) = \{ \boldsymbol{x} \in \mathbb{R}_+^2 / c(t-T) < r^*(\boldsymbol{x}) \le ct \} \\ \\ \Omega_2(t) = \{ \boldsymbol{x} \in \mathbb{R}_+^2 / r^*(\boldsymbol{x}) \le c(t-T) \} \end{array} \right. \quad (13)$$

These two sets are represented on figures 6 and 7 for two values of $t$. Note that $\Omega_1(t)$ is not empty as soon as $t > (h + 2L)/c$ while $\Omega_2(t)$ is not empty as soon as $t > (h + 2L)/c + T$.

<figure><img src="image_7.png" /><figcaption>Fig. 6. The sets $\Omega_1(t)$ and $\Omega_2(t)$, if $(h + 2L)/c \le t < T + (h + 2L)/c$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. The sets $\Omega_1(t)$ and $\Omega_2(t)$, if $t \ge T + (h + 2L)/c$.</figcaption></figure>

According to (12), in order to derive an $L^\infty$ estimate of $e^{\sigma,L}(.,t)$, we need an upper bound for the quantity:

$$\sup_{x \in \Omega_1(t)} \|G_r^N(x, \cdot)\|_{L^1(\frac{x^*}{c}, t)} \quad \text{when } t > \frac{h+2L}{c},$$

and for the quantity:

$$\sup_{x \in \Omega_2(t)} \|G_r^N(x, \cdot)\|_{L^1(t-T, t)} \quad \text{when } t > \frac{h+2L}{c} + T.$$
We thus first have to estimate $L^1$-norms in time of $G_r^{\sigma,L}(\mathbf{x}, \cdot)$. Estimating the cosine function by 1, we get:

$$|G_r^{\sigma,L}(\mathbf{x}, t)| < \frac{e^{-A^*(\mathbf{x},t)}}{2\pi\sqrt{t^2 - \frac{r^*(\mathbf{x})^2}{c^2}}}.$$

Let us remark that the function $t \mapsto e^{-A^*(\mathbf{x},t)}$ is decreasing for $t > 0$, thus, recalling that, for $x_1 < 0$, $\Sigma^*(x_1) = 2\Sigma(L) = 2L\bar{\sigma}$:

$$ \left\| G_r^{\sigma,L}(\mathbf{x}, \cdot) \right\|_{L^1(\frac{r^*}{c}, t)} \le \frac{1}{2\pi} e^{-A^*(\mathbf{x}, \frac{r^*}{c})} \int_{\frac{r^*}{c}}^{t} \frac{d\tau}{\sqrt{t^2 - \frac{r^{*2}}{c^2}}}, \quad (14) \\ = \frac{1}{2\pi} e^{-\frac{2L}{c}\bar{\sigma} |\cos \theta^*|} \log \left( \frac{ct + \sqrt{c^2 t^2 - r^{*2}}}{r^*} \right) $$

in the same way:

$$ \|G_r^{\sigma,L}(\mathbf{x}, \cdot)\|_{L^1(t-T,t)} \le \frac{e^{-2L\bar{\sigma}|\cos \theta^*|\frac{t-T}{r^*}}}{2\pi} \log \left( \frac{ct + \sqrt{c^2 t^2 - r^{*2}}}{c(t-T) + \sqrt{c^2(t-T)^2 - r^{*2}}} \right). \quad (15) $$

Since $\cos \theta^*(\mathbf{x}) = (x_1 - 2L - h)/r^*(\mathbf{x})$ we easily verify that

$$ \inf_{\mathbf{x} \in \Omega_1(t)} |\cos \theta^*(\mathbf{x})| = \frac{2L+h}{ct} \quad \text{and} \quad \inf_{\mathbf{x} \in \Omega_2(t)} \frac{|\cos \theta^*(\mathbf{x})|}{r^*(\mathbf{x})} = \frac{2L+h}{c^2(t-T)^2}, $$

hence

$$ \sup_{\mathbf{x} \in \Omega_1(t)} e^{-\frac{2L\bar{\sigma}}{c} |\cos \theta^*(\mathbf{x})|} = e^{-\frac{2L(2L+h)\bar{\sigma}}{c^2 t}} \quad \text{and} \quad \sup_{\mathbf{x} \in \Omega_2(t)} e^{-2L\bar{\sigma} \sin \theta^* \frac{t-T}{r^*}} = e^{-\frac{2L(2L+h)\bar{\sigma}}{c^2 t}}. $$

On the other hand, using the fact that the two functions:

$$ \begin{cases} r \mapsto \frac{ct}{r} + \sqrt{\frac{c^2 t^2}{r^2} - 1}, & r \in [0, ct], \\ r \mapsto \frac{ct + \sqrt{c^2 t^2 - r^2}}{c(t-T) + \sqrt{c^2(t-T)^2 - r^2}}, & r \in [0, c(t-T)], \quad (t > T) \end{cases} $$

are respectively decreasing and increasing, we deduce that:

$$ \sup_{x \in \Omega_1(t)} \mathrm{Log} \left| \frac{ct}{r^*} + \sqrt{\frac{c^2 t^2}{r^{*2}} - 1} \right| = \Phi(t) $$

$$ \text{and} \quad \sup_{x \in \Omega_2(t)} \mathrm{Log} \left| \frac{ct + \sqrt{c^2 t^2 - r^{*2}}}{c(t-T) + \sqrt{c^2(t-T)^2 - r^{*2}}} \right| = \Phi(t). $$
It is now obvious to obtain :

$$
\left\{
\begin{array}{ll}
\displaystyle \sup_{x \in \Omega_1(t)} \|G_r^{\sigma,L}(x, \cdot)\|_{L^1(\frac{r^*}{c}, t)} & \le \frac{1}{2\pi} e^{-\frac{2L(2L+h)\bar{\sigma}}{c^2t}} \Phi(t) \\[1.5ex]
\displaystyle \sup_{x \in \Omega_2(t)} \|G_r^{\sigma,L}(x, \cdot)\|_{L^1(t-T, t)} & \le \frac{1}{2\pi} e^{-\frac{2L(2L+h)\bar{\sigma}}{c^2(t-T)}} \Phi(t)
\end{array}
\right.
\qquad (16)
$$

Since, for $t > \frac{2L+h}{c} + T$, $\exp(-1/t) > \exp(-1/(t-T))$, it is easy to obtain the uniform estimates from the first inequality.

1.3 Numerical Results

**The Green's function.** The main difficulty, if one wants to represent numerically the function $G_r^{\sigma,L}(x, \cdot)$ is the singularity on the circle $r^*(x) = ct$. To overcome this problem we introduce the relative error defined by :

$$
\gamma_r^{\sigma,L}(\mathbf{x}, t) = \frac{G_r^{\sigma,L}(\mathbf{x}, \cdot)}{G_r^{0,L}(\mathbf{x}, \cdot)} = e^{-A(\mathbf{x},t)} \cos[B(\mathbf{x},t)]. \quad (17)
$$

Note that $G_r^{0,L}(x, \cdot)$ represents the field we would have obtained by imposing a Neumann condition on the line $x_1 = L$ without adding a damping term. That is why $\gamma^{\sigma,L}$ is called a relative error. We choose $h = 1$ and $c = 1$ The figure 8 represents the level curves of $\gamma^{\sigma,L}$ for three different values of $\bar{\sigma}$ ($\bar{\sigma}$=10, 20 and 50 from the left to the right). We remark that the amplitude of the error decrease strongly with $\bar{\sigma}$ : the error level is 0.6 for $\bar{\sigma}$ = 10, 0.3 for $\bar{\sigma}$ = 20 and 0.06 for $\bar{\sigma}$ = 50. When $\bar{\sigma}$ increases, the amplitude of the error concentrates more and more at the neighborhood of the layer. Moreover, its dependence with respect to the space variable is much more complicated because of the influence of the cosine function.

<figure><img src="image_9.png" /><figcaption>Fig. 8. $x \to \gamma_{r}^{\sigma,0.1}(x)$, $\bar{\sigma} = 10, 20, 50$.</figcaption></figure>

**The case of a source term.** We have implemented a Matlab code to com-
pute the convolution integral (11). To validate our solution we have compared it to the one obtained by a finite difference code (based on a standard central finite difference scheme of leap-frog type, second order accurate in space and time).

In our experiment, the source function is a truncated first derivative of a Gaussian:

$$f(t) = \frac{d}{dt} \{e^{-2\pi f_0 (t-t_0)^2}\} H(2t_0 - t), \quad f_0 = 10, \ t_0 = 1/f_0 \quad (18)$$

On figures 9 and 10 we have compared the “analytical” solution (top picture in each figure) to the numerical one (bottom picture in each figure) for two values of $\bar{\sigma}$: $\bar{\sigma} = 1.5$ and 3. In each picture we represent the level lines of the solution at time $t = 0.4$. The left pictures represent the total field while the right pictures represent the reflected field (the error). For the representation, the reflected field has been amplified by a factor which depends on $\bar{\sigma}$: 4 for $\bar{\sigma} = 1.5$ and 10 for $\bar{\sigma} = 3$. In each case, the results reveal a very good agreement between the two solutions.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between analytical (top picture) and numerical (bottom picture) solution $\bar{\sigma} = 1.5$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison between analytical (top picture) and numerical (bottom picture) solution $\bar{\sigma} = 3$.</figcaption></figure>

On figures 11 and 12 we have compared both solutions at point $(0.9, 0.1)$ as functions of time. The solid curves represent the “analytical” solution and the dashed curves the numerical one for two values of $\bar{\sigma}$: $\bar{\sigma} = 1.5$ and 3. As before the left pictures represent the total field while the right pictures represent the reflected field.

**$L^\infty$ error estimates.** On figure 13 we have compared the $L_\infty$ norm of the reflected field (the solid curves) to the uniform estimates (8) given by Theorem 4 (the dashed curves) for $\bar{\sigma} = 1, 10, 100$. The source is a step function in time: $f(t) = 1$ if $0 \le t \le 2$ and $f(t) = 0$ otherwise. Our estimate appears to be
<figure><img src="image_12.png" /><figcaption>Fig. 11.a. The total field.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11.b. The reflected field.</figcaption></figure>

Fig. 11. Comparison between analytical (top picture) and numerical (bottom picture) solution, $\sigma = 1.5$

<figure><img src="image_14.png" /><figcaption>Fig. 12.a. The total field.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12.b. The reflected field.</figcaption></figure>

Fig. 12. Comparison between analytical (top picture) and numerical (bottom picture) solution, $\sigma = 3$

very sharp for $\bar{\sigma} = 1$ and becomes less accurate (although quite acceptable) when $\bar{\sigma}$ increases.

<figure><img src="image_16.png" /><figcaption>Fig. 13.a. $\bar{\sigma} = 1$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13.b. $\bar{\sigma} = 10$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 13.c. $\bar{\sigma} = 100$.</figcaption></figure>

Fig. 13. Error estimates (dashed curve) and effective error (solid curve)

**The case of a rectangular domain surrounded by PML.** We consider a domain $[0; 1]^2$ surrounded by a PML of width $L = 0.1$. We use a point source in space located at the center of the domain, the expression of the source in time is given by (18). As for the case of the half-plane, we have compared our results to a numerical solution obtained by the same code as for the previous paragraph. These results are represented in figures 14, 15 and 16, respectively at time 7, 12 and 17 (take care to the fact that the color
scale differs from one picture to another). Here again the results show a very good agreement between the analytical and numerical solution, even for long time, after several reflections.

<figure><img src="image_19.png" /><figcaption>Fig. 14. The analytical (left picture) and numerical (right picture) solutions at time 7</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 15. The analytical (left picture) and numerical (right picture) solutions at time 12</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 16. The analytical (left picture) and numerical (right picture) solutions at time 17</figcaption></figure>

## 2 Classical PML for advective acoustics

The Bérenger's PML were adapted to the aeroacoustics equations for instance by Hu [22] in 1996. Let us first recall his model, considering an horizontal
uniform mean flow M. We have to rewrite the aeroacoustics equations whose unknowns are the pressure and the velocity $(p, U = (u, v))$,

$$ egin{cases} 
\frac{\partial p}{\partial t} + M \frac{\partial p}{\partial x_1} + \frac{\partial u}{\partial x_1} + \frac{\partial v}{\partial x_2} = 0, \\ 
\frac{\partial u}{\partial t} + M \frac{\partial u}{\partial x_1} + \frac{\partial p}{\partial x_1} = 0, \\ 
\frac{\partial v}{\partial t} + M \frac{\partial v}{\partial x_1} + \frac{\partial p}{\partial x_2} = 0, 
\end{cases} 
\quad (19) $$

in a so-called “split form” with the unknowns $(p_1, p_2, u, v_1, v_2)$ where $p_1$ and $p_2$ (resp. $v_1$ and $v_2$) are non physical variables, whose sum is equal to the pressure (resp. the second component of the velocity). One then obtains the PML model by adding a zero-order absorption term, for instance in the equations containing the derivatives with respect to $x_1$ for a vertical layer :

$$ egin{cases} 
\frac{\partial p_1}{\partial t} + \sigma(x_1)p_1 + M \frac{\partial p}{\partial x_1} + \frac{\partial u}{\partial x_1} = 0, \\ 
\frac{\partial p_2}{\partial t} + \frac{\partial v}{\partial x_2} = 0, \\ 
\frac{\partial u}{\partial t} + \sigma(x_1)u + M \frac{\partial u}{\partial x_1} + \frac{\partial p}{\partial x_1} = 0, \\ 
\frac{\partial v_1}{\partial t} + \sigma(x_1)v_1 + M \frac{\partial v}{\partial x_1} = 0, \\ 
\frac{\partial v_2}{\partial t} + \frac{\partial p}{\partial x_2} = 0. 
\end{cases} 
\quad (20) $$

Though this problem is (weakly) well-posed (see [24]), it is unstable as soon as $|M| > 0$ : it may admit solution growing exponentially with the time variable.

## 2.1 Analysis of the instabilities via the slowness curve

The instabilities can be analyzed thanks to a plane wave analysis and to the slowness curves : we search for $(\omega, k_1, k_2)$ such that

$$ [p_1, p_2, u, v_1, v_2]^t(x_1, x_2, t) = [p_1^0, p_2^0, u^0, v_1^0, v_2^0]^t e^{i(\omega t + k_1 x_1 + k_2 x_2)} $$

It is well-known that, if $\sigma = 0$ (i.e for the aeroacoustics equations), $(\omega, k_1, k_2)$ is solution of one of the two following dispersion equations :

* the dispersion equation of the advective waves :

$$ (\omega + Mk_1)^2 - k_1^2 - k_2^2 = 0; \quad (21) $$
• the dispersion equation of the vorticity waves :

$$ \omega + Mk_1 = 0. \qquad (22) $$

The slowness curves are the representation in the so-called slowness plane $(k_1/\omega, k_2/\omega)$ of the solutions of these equations. One can easily show that the slowness curve associated to the advective waves is an ellipse whose center is $(M/(1 - M^2), 0)$ and that the slowness curve associated to the vorticity is a line of equation $k_1/\omega = -1/M$. Let us also define the slowness vector and the group velocity respectively by : $l = [k_1/\omega, k_2/\omega]^t$ and $v_g = \nabla_k \omega$. It is immediate to deduce the slowness vector from the slowness curve. Moreover the group velocity vector is orthogonal to the slowness curves (see figure 17). It is well-known that, as soon as $|M| > 0$, it exists waves whose first component of slowness and of group velocity vector have opposite signs (in red in figure 17).

It has been shown [25,2,10] that the high-frequency instabilities are due to these so-called back-propagating modes. It is thus easy to deduce the set of unstable waves : the advective waves such that $0 < k_1/\omega < M/(1 - M^2)$ (see figure 18). An illustration of this instabilities is given in figure 19 for $M = 0.5$.

<figure><img src="image_22.png" /><figcaption>Fig. 17. The slowness and group velocity vectors</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 18. The set of unstable points</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 19. The instabilities in the PML, t=1, 3, 8, 30 and 50</figcaption></figure>
## 2.2 The main results

As for the acoustics equation, it is possible to obtain an analytic expression of the fundamental solution of problem 20. Let us first remark that this problem can be formally rewritten as (the calculations are left to the reader) :

$$ 
\left\{
\begin{aligned}
\frac{\partial p}{\partial t} + MD_{x_1}^{\sigma} p + D_{x_1}^{\sigma} u + \frac{\partial v}{\partial x_2} &= 0, \\
\frac{\partial u}{\partial t} + MD_{x_1}^{\sigma} u + D_{x_1}^{\sigma} p &= 0, \\
\frac{\partial v}{\partial t} + MD_{x_1}^{\sigma} v + \frac{\partial p}{\partial x_2} &= 0.
\end{aligned}
\right.
\qquad (23)
$$

Then, straightforward calculations lead us to the following decoupled system :

$$ 
\left\{
\begin{aligned}
\frac{\partial^2 p}{\partial t^2} + 2MD_{x_1}^{\sigma} \frac{\partial p}{\partial t} - (1 - M^2)(D_{x_1}^{\sigma})^2 p - \frac{\partial^2 p}{\partial x_2^2} &= 0, \\
\frac{\partial u}{\partial t} + MD_{x_1}^{\sigma} u + D_{x_1}^{\sigma} p &= 0, \\
\frac{\partial v}{\partial t} + MD_{x_1}^{\sigma} v + \frac{\partial p}{\partial x_2} &= 0.
\end{aligned}
\right.
\qquad (24)
$$

The first equation is nothing but the advective wave equation where the derivation with respect to $x_1$ has been replaced by the operator $D_{x_1}^{\sigma}$. This equation, which we choose to study here, is the cause of the instabilities.

As for the acoustics case we consider the problem

$$ 
\left\{
\begin{array}{ll}
\text{Find } u : \mathbb{R}^2 \times \mathbb{R}_+ \mapsto \mathbb{R}, \text{ null for } t < 0 & \text{such that} \\
\frac{\partial^2 p}{\partial t^2} + 2M \frac{\partial}{\partial t} D_{x_1}^{\sigma} p - \left[ (1 - M^2)(D_{x_1}^{\sigma})^2 p + \frac{\partial^2}{\partial x_2^2} p \right] = \delta(\mathbf{x} - \mathbf{x}_S)\delta(t) \text{ in } \mathbb{R}^2 \times \mathbb{R}_+ & \\
\text{where } \mathbf{x}_S = (-h, 0), \sigma(x_1) = 0 \text{ if } x_1 \le 0 \text{ and } \sigma(x_1) = 0 \text{ if } x_1 > 0. & (25)
\end{array}
\right.
$$

We define the functions $r(\mathbf{x})$ and $\theta(\mathbf{x})$ ("generalized polar coordinates") as :

$$ r^2(\mathbf{x}) = \frac{(x_1 + h)^2}{(1 - M^2)^2} + \frac{x_2^2}{1 - M^2}, \cos \theta(\mathbf{x}) = \frac{x_1 + h}{r(1 - M^2)} \text{ and } \sin \theta(\mathbf{x}) = \frac{x_2}{r\sqrt{1 - M^2}}. $$

Note that the level lines of $r(\mathbf{x})$ and $\theta(\mathbf{x})$ are ellipses and half-lines (see Fig.20).
<figure><img src="image_25.png" /><figcaption>Fig. 20. The level lines of $r(\boldsymbol{x})$ and $\theta(\boldsymbol{x})$</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 21. $r \mapsto \gamma(r, \theta_0, t)$</figcaption></figure>

We also introduce the two functions $A(\boldsymbol{x}, t)$ and $B(\boldsymbol{x}, t)$ defined by

$$A(\boldsymbol{x}, t) = \Sigma(x_1) \left( \frac{t}{r} \cos \theta - M \sin^2 \theta \right)$$

$$\text{and } B(\boldsymbol{x}, t) = \Sigma(x_1) \sin \theta \sqrt{\left(\frac{t}{r} + M \cos \theta\right)^2 - 1}.$$

**Theorem 6** The expression of the solution $p(\boldsymbol{x}, t) = G_i^\sigma(\boldsymbol{x}, t)$ of problem (25) is given by :

$$\left| \begin{array}{l} G_i^\sigma(\boldsymbol{x}, t) = \frac{H(t - r(1 - M \cos \theta))}{2\pi\sqrt{1 - M^2}\sqrt{(t + M r \cos \theta)^2 - r^2}}, \quad x_1 \le 0 \\ \\ G_i^\sigma(\boldsymbol{x}, t) = \frac{H(t - r(1 - M \cos \theta))}{2\pi\sqrt{1 - M^2}\sqrt{(t + M r \cos \theta)^2 - r^2}} e^{-A(\boldsymbol{x}, t)} \cos [B(\boldsymbol{x}, t)], \quad x_1 > 0. \end{array} \right. \qquad (26)$$

For $x_1 > 0$ and $t > t_0 = r(1 - M \cos \theta)$ we define the function

$$\gamma(\boldsymbol{x}, t) = \frac{G_i^\sigma(\boldsymbol{x}, t)}{G_i^0(\boldsymbol{x}, t)} = e^{-A(\boldsymbol{x}, t)} \cos [B(\boldsymbol{x}, t)].$$

**Theorem 7** For fixed $\boldsymbol{x}$, the function $t \to \gamma(\boldsymbol{x}, t)$ is bounded : $\gamma(\boldsymbol{x}, t) \le C(\boldsymbol{x})$. However, given a direction $\theta_0$, let us denote

$$r^{\max}(t) = \frac{t}{(1 - M \cos \theta_0)},$$

then, $\gamma(r^{\max}(t), \theta_0, t) > e^{\alpha(\theta_0)t}$ with $\alpha(\theta_0) > 0$ if $\cos \theta_0 < M$ and $\alpha(\theta_0) \le 0$ if $\cos \theta_0 \ge M$. Consequently the $L^\infty(\mathbb{R}^+)$-norm of $\gamma(\boldsymbol{x}, t)$ grows exponentially
with t :

$$ \|\gamma(\cdot, t)\|_{L^{\infty}} \geq C e^{\alpha t}, \quad \alpha > 0. $$

The kind of instability expressed by this theorem is often called “convective instability” because the instability is “moving with the wave front”. To illustrate this point we have represented the variations of $\gamma$ with respect to $r$ for a given unstable direction $\theta_0$ and for four values of $t$ in figure 21.

We have also implemented a Matlab code to compute the solution of problem (25) where the source function is once again a truncated first derivate of a Gaussian. In the following experiment the frequency of the source is 1 and the source is located at point $(0, -2)$. Figures (22) and (23) represents the propagation of an advective wave, respectively with $M = 0.5$ and $M = 0.9$. We have represented in red the half lines $\theta(x) = \arccos(M)$ and $\theta(x) = -\arccos(M)$ that define the set of unstable points.

<figure><img src="image_27.png" /><figcaption>Fig. 22. Instabilities in the PML with $M = 0.5$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 23. Instabilities in the PML with $M = 0.9$</figcaption></figure>

## 2.3 Proof of the theorems

**Expression of the analytic solution.** Let us first apply to $G_i^\sigma$ the Laplace transform in $t$ and the Fourier transform in $x_2$. The resulting function $x_1 \mapsto$
$\hat{G}_i^\sigma(x_1, k, s)$ satisfies the ordinary differential equation :

$$ -(1 - M^2) \frac{s}{s + \sigma} \frac{d}{dx_1} \left( \frac{s}{s + \sigma} \frac{d\hat{G}_i^\sigma}{dx_1} \right) + 2Ms \frac{s}{s + \sigma} \frac{\partial \hat{G}_i^\sigma}{\partial x_1} + (k^2 + s^2) \hat{G}_i^\sigma = \delta(x_1 + h) \quad (27) $$

We denote

$$ \lambda(x_1) = \begin{cases} M + \sqrt{1 + (1 - M^2) \frac{k^2}{s^2}}, & \text{if } x_1 < -h, \\ M - \sqrt{1 + (1 - M^2) \frac{k^2}{s^2}}, & \text{if } x_1 > -h, \end{cases} $$

$\hat{G}_i^\sigma(x_1, k, s)$ has the following form

$$ \hat{G}_i^\sigma(x_1, k, s) = \mathcal{A}(k, s) e^{\frac{s(x_1+h)+\Sigma(x_1)}{1-M^2}\lambda(x_1)} $$

Taking into account the jump condition

$$ \left[ -(1 - M^2) \left( \frac{s}{s + \sigma} \right)^2 \frac{\partial \hat{G}_i^\sigma}{\partial x_1} + 2M \frac{s}{s + \sigma} \hat{G}_i^\sigma \right]_{x_1 = -h} = 1, \quad (28) $$

one obtains, if $\sigma(-h) = 0$, $\mathcal{A}(k, s) = \frac{1}{2\sqrt{s^2 + (1 - M^2)k^2}}$.

If $x_1$ is negative, $\Sigma(x_1) = 0$ and $\hat{G}_i^\sigma(x_1, k, s) = \frac{e^{\frac{s(x_1+h)}{1-M^2}\lambda(x_1)}}{2\sqrt{1+(1-M^2)\frac{k^2}{s^2}}}$.

This is nothing but the Fourier-Laplace transform of the fundamental solution of the advective wave equation :

$$ G_i^\sigma(\boldsymbol{x}, t) = \frac{H(c t - r(1 - M \cos \theta))}{2\pi\sqrt{(t + M r \cos \theta)^2 - r^2}}, \quad x_1 \le 0. $$

Let us now consider $x_1 > 0$ so that $\lambda(x_1) = M - \sqrt{1 + (1 - M^2)\frac{k^2}{s^2}}$ and

$$ \hat{G}_i^\sigma(x_1, k, s) = \frac{e^{\frac{s(x_1+h)+\Sigma(x_1)}{1-M^2}\left(M - \sqrt{1+(1-M^2)\frac{k^2}{s^2}}\right)}}{2\sqrt{s^2 + (1-M^2)k^2}} $$

Applying the $x_2$ inverse Fourier transform to $\hat{G}_i^\sigma(x_1, k, s)$, we get

$$ \tilde{G}_i^\sigma(x_1, x_2, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{e^{\frac{s(x_1+h)+\Sigma(x_1)}{1-M^2}\left(M - \sqrt{1+(1-M^2)\frac{k^2}{s^2}}\right)}}{\sqrt{s^2 + (1-M^2)k^2}} e^{-ikx_2} dk, \quad (29) $$
We apply the change of variable $k = \frac{ps}{\sqrt{1-M^2}}$ and obtain:

$$ \tilde{G}_i^\sigma(x_1, x_2, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{e^{\Sigma(x_1) \frac{M-\sqrt{1+p^2}}{1-M^2}} e^{-s \left[ \frac{(M-\sqrt{1+p^2})}{1-M^2} (x_1+h) - ip \frac{x_2}{\sqrt{1-M^2}} \right]}}{\sqrt{(1-M^2)(1+p^2)}} dp. \quad (30) $$

We now have to search a path in the complex plane such that

$$ \frac{M - \sqrt{1 + p^2}}{1 - M^2} (x_1 + h) - ip \frac{x_2}{\sqrt{1 - M^2}} = -t, \quad \text{for } t > 0, \quad (31) $$

For $t \ge t_0 = r(1 - M \cos \theta)$ we define the functions $\gamma^{\pm}(t)$ by

$$ \gamma^{\pm}(t) = -i \left( \frac{t}{r} + M \cos \theta \right) \sin \theta \pm \cos \theta \sqrt{\left( \frac{t}{r} + M \cos \theta \right)^2 - 1}, $$

and we introduce the Cagniard-de Hoop contour $\Gamma$ defined as:

$$ \Gamma = \Gamma^+ \cup \Gamma^-, \quad \Gamma^{\pm} = \{ p = \gamma^{\pm}(t) \}. \quad (32) $$

The contour $\Gamma$ is very similar to the one exposed at the previous section (consider denoting $T = t + Mr \cos \theta$) : in particular it does not intercept the branch cut of the function $(1+p^2)^{-1/2}$. Then, as for the acoustics case, we easily deduce :

$$ \tilde{G}_i^\sigma(x_1, x_2, s) = -\frac{1}{4\pi} \int_{\Gamma} \frac{e^{\Sigma(x_1) \frac{M-\sqrt{1+p^2}}{1-M^2}} e^{-s \left[ \frac{(M-\sqrt{1+p^2})}{1-M^2} (x_1+h) - ip \frac{x_2}{\sqrt{1-M^2}} \right]}}{\sqrt{1+p^2}} dp. $$

We use the parameterizations $p = \gamma^+(t)$ and $p = \gamma^-(t)$, for $t \ge t_0$, respectively along $\Gamma^+$ and $\Gamma^-$ and remark that:

$$ \bullet \quad \frac{(M - \sqrt{1 + p^2})}{1 - M^2} (x_1 + h) - ip \frac{x_2}{\sqrt{1 - M^2}} = -t, \quad (\text{by construction}), $$

$$ \bullet \quad \frac{dp}{(1 + p^2)^{\frac{1}{2}}} = \pm \frac{dt}{\left[ (t + M r \cos \theta)^2 - r^2 \right]^{\frac{1}{2}}}, \quad \text{on } \Gamma^{\pm}, $$

$$ \bullet \quad \Sigma(x_1) \frac{M - \sqrt{1 + p^2}}{1 - M^2} = -A(\boldsymbol{x}, t) \mp iB(\boldsymbol{x}, t), \quad \text{on } \Gamma^{\pm} $$

Therefore, since $t$ goes from $+\infty$ to $t_0$ on $\Gamma^+$ and from $t_0$ to $+\infty$ on $\Gamma^-$:

$$ \tilde{G}_i^\sigma(x_1, x_2, s) = \frac{1}{4\pi} \int_{t_0}^{+\infty} \frac{e^{-A(\boldsymbol{x}, t)} \cos(B(\boldsymbol{x}, t))}{\sqrt{1 - M^2} \sqrt{(t + M r \sin \theta)^2 - r^2}} e^{-st} dt. $$
It is then easy to conclude.

**Properties of the function $\gamma$.** For fixed $x$, the function $B(x, t)$ is bounded by 1 in time. Moreover the function $t \to A(x, t)$ is obviously increasing with $t$ on $[t_0; +\infty]$ and its minimum is then $A(x, t_0) = \Sigma(x_1)(\cos \theta - M)$. Then

$$ \max_{t \ge t_0} \gamma(x, t) \le e^{-\Sigma(x_1)(\cos \theta - M)}. $$

For given $t$ we denote by $\Omega(t)$ the support of $\gamma(., t)$, it is a disk of center $(Mt - h, 0)$ and of radius $t$. Let us remark that the point of generalized polar coordinates $(r^{\max}(t), \theta_0)$ belongs to $\partial\Omega_t$. Denoting

$$ x_1^{\max}(t) = (1 - M^2)r^{\max}(t) \cos \theta_0 - h = \max_{\Delta(\theta_0, t)} x_1, $$

it is easy to check that

$$ \gamma(r^{\max}(t), \theta_0, t) = e^{-\Sigma(x_1^{\max}(t))(\cos \theta_0 - M)}, $$

since $B(r^{\max}(t), \theta_0, t) = 0$, then, since $x_1^{\max}(t)$ grows linearly with $t$ and $\Sigma$ grows with $x_1$ it is easy to conclude the proof

# 3 Stabilized PML for advective acoustics

Various solutions have been proposed in the literature to prevent instabilities in the PML. Hu [5] used a low-pass filter inside the absorbing layer. Tam, Auriault and Cambuli [2] proposed to use selective damping coefficients. Hesthaven [25] decelerated progressively the flow in the layer. Lions, Metral and Vacus [26] defined a new layer by regularizing the damping terms.

However all these methods may be not perfectly matched and, although these problems are well-posed, their stability has not been proven. Recently a significative advance has been achieved by Abarbanel, Gottlieb and Hesthaven [21] and by Hu [22] who has presented a stable and perfectly matched model, close to the one we are going to present (though the two models have been obtained independently). Hagström [20], using a different method, has also constructed a similar layer.

## 3.1 The advective wave equation

We are first concerned by stabilizing the PML for the advective equation whose associated slowness curve is the ellipse represented in the figure 18. As we said before, the instabilities produced by the PML can be explained by the
fact that the center of this ellipse is not the origin of the axes. Therefore an idea for stabilizing the PML consists in using a change of variable to translate the ellipse by the vector $(-M/(1 - M^2), 0)$ for its center to be located at the origin (see figure 24). Therefore we introduce the new unknowns $k_1^*$ and $\omega^*$

<figure><img src="image_29.png" /><figcaption>Fig. 24. The slowness curve before the change of variable (in blue/solid) and after (in magenta/dashed)</figcaption></figure>

such that

$$ \frac{k_1^*}{\omega^*} = \frac{k_1}{\omega} - \frac{M}{1 - M^2}. $$

If we set $\omega^* = \omega$ this relation becomes

$$ k_1^* = k_1 - \frac{M}{1 - M^2} \omega. $$

This change of variable can also be expressed through the time-space variables :

$$ \frac{\partial}{\partial x_1} - \frac{M}{1 - M^2} \frac{\partial}{\partial t} \Rightarrow \frac{\partial}{\partial x_1^*}, $$

which is obtained by using the new variables :

$$ x_1^* = x_1, \quad x_2^* = x_2, \quad \text{and} \quad t^* = t - \frac{M}{1 - M^2} x_1. \qquad (33) $$

Note that this change of variable does not affect the spatial coordinates and consequently the interface between the computational domain and the layer. The choice $\omega^* = \omega$ has been made to obtain this property.

**Construction of the stabilized PML model.** Let us first show how to obtain an infinite PML in $\mathbb{R}_+^2$ ($x_1 > 0$). The steps of the methods are :

(1) Rewrite the problem as a transmission problem between $\mathbb{R}_-$ and $\mathbb{R}_+^2$ :
$$
\left\{
\begin{aligned}
& \frac{\partial^2 p}{\partial t^2} + 2M \frac{\partial p}{\partial t \partial x_1} - (1 - M^2) \frac{\partial^2 p}{\partial x_1^2} - \frac{\partial^2 p}{\partial x_2^2} = \delta(\mathbf{x} - \mathbf{x}_s)\delta(t), \quad x_1 < 0, \quad (34a) \\
& \frac{\partial^2 p}{\partial t^2} + 2M \frac{\partial p}{\partial t \partial x_1} - (1 - M^2) \frac{\partial^2 p}{\partial x_1^2} - \frac{\partial^2 p}{\partial x_2^2} = 0, \quad x_1 > 0, \quad (34b) \\
& p(0^{-}, x_2, t) = p(0^{+}, x_2, t), \quad (34c) \\
& \frac{\partial p}{\partial x_1}(0^{-}, x_2, t) = \frac{\partial p}{\partial x_1}(0^{+}, x_2, t). \quad (34d)
\end{aligned}
\right.
$$

(2) Apply the change of variable

$$p^*(x_1, x_2, t) = p(x_1, x_2, t - \frac{M}{1 - M^2}x_1),$$

in the PML to obtain the new transmission problem :

$$
\left\{
\begin{aligned}
& \frac{\partial^2 p}{\partial t^2} + 2M \frac{\partial p}{\partial t \partial x_1} - (1 - M^2) \frac{\partial^2 p}{\partial x_1^2} - \frac{\partial^2 p}{\partial x_2^2} = \delta(\mathbf{x} - \mathbf{x}_s)\delta(t), \quad x_1 < 0, \quad (35a) \\
& \frac{1}{1 - M^2} \frac{\partial^2 p^*}{\partial t^2} - (1 - M^2) \frac{\partial^2 p^*}{\partial x_1^2} - \frac{\partial^2 p^*}{\partial x_2^2} = 0 \quad x_1 > 0, \quad (35b) \\
& p(0, x_2, t) = p^*(0, x_2, t), \quad (35c) \\
& \frac{\partial p}{\partial x_1}(0, x_2, t) = \frac{\partial p^*}{\partial x_1}(0, x_2, t) + \frac{M}{1 - M^2} \frac{\partial p^*}{\partial t}(0, x_2, t). \quad (35d)
\end{aligned}
\right.
$$

Note that the advective wave equation has been turned into an equation very close to the classical wave equation by this change of variable.

Let us now define the function $\rho(\mathbf{x})$ by $\rho(\mathbf{x}) = 1$ if $x_1 < 0$ and $\rho(\mathbf{x}) = 1/(1 - M^2)$ otherwise.

**Theorem 8** *The following energy is conserved :*

$$E(t) = \frac{1}{2} \int_{\mathbb{R}^2} \left[ \rho(\mathbf{x}) \left| \frac{\partial p}{\partial t} \right|^2 + (1 - M^2) \left| \frac{\partial p}{\partial x_1} \right|^2 + \left| \frac{\partial p}{\partial x_2} \right|^2 \right].$$

**Proof.** Let us respectively multiply (35a) and (35b) by $\partial p/\partial t$ and $\partial p^*/\partial t$. After having integrated both this equations over their respective defini-
tion domain ($\mathbb{R}^2_-$ and $\mathbb{R}^2_+$) and added them, we check that :

$$ \begin{aligned} \frac{dE}{dt}(t) + \int_{\mathbb{R}^2_-} 2M \frac{\partial p}{\partial t \partial x_1} \frac{\partial p}{\partial t} + \int_{\Gamma} (1-M^2) \left( \frac{\partial p^*}{\partial x_1} \frac{\partial p^*}{\partial t} - \frac{\partial p}{\partial x_1} \frac{\partial p}{\partial t} \right) \\ + \int_{\Gamma} \left( \frac{\partial p^*}{\partial x_2} \frac{\partial p^*}{\partial t} - \frac{\partial p}{\partial x_2} \frac{\partial p}{\partial t} \right) = 0 \end{aligned} $$

where $\Gamma$ is the line of equation $x_1 = 0$. Using the continuity of $p$ (and consequently of its derivatives with respect to $x_2$ ant $t$) it is obvious that the last integral in this equation is zero. Furthermore

$$ \int_{\mathbb{R}^2_-} M \frac{\partial p}{\partial t \partial x_1} \frac{\partial p}{\partial t} = \int_{\Gamma} M \frac{\partial p}{\partial t} \frac{\partial p}{\partial t} - \int_{\mathbb{R}^2_-} M \frac{\partial p}{\partial t} \frac{\partial p}{\partial t \partial x_1} $$

Then

$$ \begin{aligned} & \int_{\mathbb{R}^2_-} 2M \frac{\partial p}{\partial t \partial x_1} \frac{\partial p}{\partial t} + \int_{\Gamma} (1-M^2) \left( \frac{\partial p^*}{\partial x_1} \frac{\partial p^*}{\partial t} - \frac{\partial p}{\partial x_1} \frac{\partial p}{\partial t} \right) \\ & = \int_{\Gamma} \frac{\partial p}{\partial t} \left[ M \frac{\partial p}{\partial t} + (1-M^2) \left( \frac{\partial p^*}{\partial x_1} - \frac{\partial p}{\partial x_1} \right) \right] \end{aligned} $$

We deduce from (35d) that this expression is zero :

$$ \frac{dE}{dt}(t) = 0. $$

(3) It only remains to replace formally the derivatives with respect to $x_1$ by the operator $D_{x_1}^\sigma$ :

$$ \frac{\partial^2 p}{\partial t^2} + 2M \frac{\partial p}{\partial t \partial x_1} - (1 - M^2) \frac{\partial^2 p}{\partial x_1^2} - \frac{\partial^2 p}{\partial x_2^2} = \delta(\boldsymbol{x} - \boldsymbol{x}_s), \quad x_1 < 0, (36a) $$

$$ \frac{1}{1 - M^2} \frac{\partial^2 p^*}{\partial t^2} - (1 - M^2) (D_{x_1}^\sigma)^2 p^* - \frac{\partial^2 p^*}{\partial x_2^2} = 0 \quad x_1 > 0, (36b) $$

$$ p(0, x_2, t) = p^*(0, x_2, t), \qquad (36c) $$

$$ \frac{\partial p}{\partial x_1} = D_{x_1}^\sigma p^*(0, x_2, t) + \frac{M}{1 - M^2} \frac{\partial p^*}{\partial t}. \qquad (36d) $$

**Analysis of the new model.** Let us define the functions $r(\boldsymbol{x})$ and $\theta(\boldsymbol{x})$ by

$$ r^2(\boldsymbol{x}) = \frac{(x_1 + h)^2}{(1 - M^2)^2} + \frac{x_2^2}{1 - M^2}, \cos \theta(\boldsymbol{x}) = \frac{x_1 + h}{r(1 - M^2)}, \text{ and } \sin \theta(\boldsymbol{x}) = \frac{x_2}{r\sqrt{1 - M^2}}. $$
Let $A(\boldsymbol{x}, t)$ and $B(\boldsymbol{x}, t)$ be defined by

$$A(\boldsymbol{x}, t) = \Sigma(x_1) \left( t + \frac{Mh}{1 - M^2} \right) \frac{\cos \theta}{r} > 0$$

$$\text{and } B(\boldsymbol{x}, t) = \Sigma(x_1) \frac{\sin \theta}{r} \sqrt{\left(t + \frac{Mh}{1 - M^2}\right)^2 - 1}.$$

**Theorem 9** The expression of the solution $(p, p^*)(\boldsymbol{x}, t) = (G_i, G_i^*)(\boldsymbol{x}, t)$ of problem (36) is given by:

$$\left\{
\begin{array}{ll}
G_i(\boldsymbol{x}, t) = \frac{H(t - r(1 - M \cos \theta))}{2\pi\sqrt{1 - M^2}\sqrt{(t + M r \cos \theta)^2 - r^2}}, & x_1 \le 0 \\
\\
G_i^*(\boldsymbol{x}, t) = \frac{H(t + \frac{Mh}{1 - M^2} - r)}{2\pi\sqrt{1 - M^2}\sqrt{\left(t + \frac{Mh}{1 - M^2}\right)^2 - r^2}} e^{-A(\boldsymbol{x}, t)} \cos[B(\boldsymbol{x}, t)], & x_1 > 0.
\end{array}
\right.
\qquad (37)$$

This expression is similar to one exposed at the previous section. The difference lies in the fact that the function $A$ is always positive which guarantee the stability of the layer.

**Proof via the Cagniard-de Hoop technique.** After applying a Laplace transform in $t$ and a Fourier transform in $x_2$ we obtain the following differential system:

$$\left\{
\begin{array}{ll}
(s^2 + k^2) \hat{G}_i + 2sM \frac{d\hat{G}_i}{dx_1} - (1 - M^2) \frac{d^2 \hat{G}_i}{dx_1} = \delta(x_1 + h), & x_1 < 0, (38a) \\
\\
\left(\frac{s^2}{1 - M^2} + k^2\right) \hat{G}_i^* - (1 - M^2) \frac{s}{s + \sigma} \frac{d}{dx_1} \left(\frac{s}{s + \sigma} \frac{d\hat{G}_i}{dx_1}\right) = 0 & x_1 > 0, (38b) \\
\\
\hat{G}_i(0) = \hat{G}_i^*(0), & \\
\\
\frac{d\hat{G}_i}{dx_1}(0) = \frac{s}{s + \sigma} \frac{d\hat{G}_i^*}{dx_1}(0) + \frac{Ms}{1 - M^2} \hat{G}_i^*(0). & (38d)
\end{array}
\right.$$

After calculations (similar to the ones of section 2) we obtain:

$$\hat{G}_i(x_1, k, s) = \frac{e^{\frac{s(x_1 + h)}{1 - M^2} \lambda(x_1)}}{2\sqrt{s^2 + (1 - M^2)k^2}},$$
$$ \hat{G}_{i}^{*}(x_{1}, k, s) = \frac{e^{\frac{1}{1-M^{2}}\left(shM - \left(x_{1} + h + \frac{\Sigma(x_{1})}{s}\right)\sqrt{s^{2} + (1-M^{2})k^{2}}\right)}}{2\sqrt{s^{2} + (1-M^{2})k^{2}}}, $$

where $\lambda$ has been defined at the section 2. As $\hat{G}_i$ is once again the Fourier-Laplace transform of the classical advective wave equation it is obvious that the new PML does not produce parasitical reflection. Applying the $x_2$ inverse Fourier transform to $\hat{G}_i^*(x_1, k, s)$, we get

$$ \tilde{G}_{i}^{*}(x_{1}, x_{2}, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{e^{\frac{1}{1-M^{2}}\left(shM - \left(x_{1} + h + \frac{\Sigma(x_{1})}{s}\right)\sqrt{s^{2} + (1-M^{2})k^{2}}\right)}}{2\sqrt{s^{2} + (1-M^{2})k^{2}}} e^{-ikx_{2}} dk, \quad (39) $$

We apply the change of variable $k = \frac{ps}{\sqrt{1-M^2}}$ and obtain:

$$ \tilde{G}_{i}^{*}(x_{1}, x_{2}, s) = \frac{1}{4\pi} \int_{-\infty}^{+\infty} \frac{e^{-\Sigma(x_{1})\frac{\sqrt{1+p^{2}}}{1-M^{2}}} e^{s\left[\frac{Mh-(x_{1}+h)\sqrt{1+p^{2}}}{1-M^{2}} - ip\frac{x_{2}}{\sqrt{1-M^{2}}}\right]}}{\sqrt{(1-M^{2})(1+p^{2})}} dp. \quad (40) $$

We now have to search a path in the complex plane such that

$$ \frac{Mh - (x_1 + h)\sqrt{1 + p^2}}{1 - M^2} - ip \frac{x_2}{\sqrt{1 - M^2}} = -t, \quad \text{for } t > 0, \quad (41) $$

For $t \ge t_0 = r - \frac{Mh}{1-M^2} > 0$ we define the functions $\gamma^{\pm}(t)$ by

$$ \gamma^{\pm}(t) = -i \left( \frac{t}{r} + \frac{Mh}{r(1-M^2)} \right) \sin\theta \pm \cos\theta \sqrt{\left( \frac{t}{r} + \frac{Mh}{r(1-M^2)} \right)^2 - 1}, $$

and we introduce the Cagniard-de Hoop contour $\Gamma$ defined as:

$$ \Gamma = \Gamma^{+} \cup \Gamma^{-}, \quad \Gamma^{\pm} = \{ p = \gamma^{\pm}(t) \}. \quad (42) $$

As for the acoustics case we easily deduce:

$$ \tilde{G}_{i}^{*}(x_{1}, x_{2}, s) = -\frac{1}{4\pi} \int_{\Gamma} \frac{e^{-\Sigma(x_{1})\frac{\sqrt{1+p^{2}}}{1-M^{2}}} e^{s\left[\frac{Mh-(x_{1}+h)\sqrt{1+p^{2}}}{1-M^{2}} - ip\frac{x_{2}}{\sqrt{1-M^{2}}}\right]}}{\sqrt{(1-M^{2})(1+p^{2})}} dp. $$
We use the parameterizations $p = \gamma^+(t)$ and $p = \gamma^-(t)$, for $t \ge t_0$, respectively along $\Gamma^+$ and $\Gamma^-$ and remark that:

$$ \bullet \quad \frac{Mh - (x_1 + h)\sqrt{1 + p^2}}{1 - M^2} - ip \frac{x_2}{\sqrt{1 - M^2}} = -t, \quad (\text{by construction}), $$

$$ \bullet \quad \frac{dp}{(1 + p^2)^{\frac{1}{2}}} = \pm \frac{dt}{\left[ \left( t + \frac{Mh}{1 - M^2} \right)^2 - r^2 \right]^{\frac{1}{2}}}, \quad \text{on } \Gamma^{\pm}, $$

$$ \bullet \quad - \Sigma(x_1) \frac{\sqrt{1 + p^2}}{1 - M^2} = -A(\boldsymbol{x}, t) \mp iB(\boldsymbol{x}, t), \quad \text{on } \Gamma^{\pm} $$

Therefore, since $t$ goes from $+\infty$ to $t_0$ on $\Gamma^+$ and from $t_0$ to $+\infty$ on $\Gamma^-$:

$$ \widetilde{G}_i(x_1, x_2, s) = \frac{1}{4\pi} \int_{t_0}^{+\infty} \frac{e^{-A(\boldsymbol{x}, t)} \cos(B(\boldsymbol{x}, t))}{\sqrt{1 - M^2} \sqrt{\left(t + \frac{Mh}{1 - M^2}\right)^2 - r^2}} e^{-st} dt. $$

It is then easy to conclude.

**Numerical illustration.** In figure 25 we compare three numerical experiments:

(1) (top picture) the propagation of an advective wave in a large computational domain:

(2) (middle picture) the same experiment but the change of variable (33) has been applied in the left and right hand layers;

(3) (bottom picture) in the last one a damping term has been added in the layers.

<figure><img src="image_30.png" /><figcaption>Fig. 25. Comparison between three experiments</figcaption></figure>

If we compare the two first experiments we remark that the change of variable (33) slows down the wave in the right hand layer whereas it accelerates in the left hand one (actually, this change of variable “makes the flow vanish”). The last experiment shows that the waves are properly damped in the PML.
## 3.2 Extension to the linearized Euler equations

Looking at the shape of slowness curve associated to the vorticity waves (a line orthogonal to the axis $k_1/\omega$), we can reasonably assume that the change of variable (33), which translate the slowness curves of $M/(1 - M^2)$ to the right (see figure 26), does not generate instabilities.

<figure><img src="image_31.png" /><figcaption>Fig. 26. The slowness curve before (in blue/solid) and after (in magenta/dashed) the change of variable</figcaption></figure>

It is useful to introduce the new unknowns $s = p + u$ and $d = p - u$, and to rewrite the aeroacoustics equation as :

$$ \left\{ \frac{\partial s}{\partial t} + (1 + M) \frac{\partial s}{\partial x_1} + \frac{\partial v}{\partial x_2} = 0, \right. \qquad (43a) $$

$$ \left\{ \frac{\partial d}{\partial t} - (1 - M) \frac{\partial d}{\partial x_1} + \frac{\partial v}{\partial x_2} = 0, \right. \qquad (43b) $$

$$ \left. \frac{\partial v}{\partial t} + M \frac{\partial v}{\partial x_1} + \frac{1}{2} \left( \frac{\partial s}{\partial x_2} + \frac{\partial d}{\partial x_2} \right) = 0. \right. \qquad (43c) $$

Then the three steps of the construction of a stable PML are :

(1) Write a transmission problem between $x_1 < 0$ and $x_1 > 0$.

(2) Apply the change of variable (33) to obtain the following new equations in the domain $x_1 > 0$:

$$ \left\{ \frac{1}{1 - M} \frac{\partial s^*}{\partial t} + (1 + M) \frac{\partial s^*}{\partial x_1} + \frac{\partial v^*}{\partial x_2} = 0, \right. \qquad (44a) $$

$$ \left\{ \frac{1}{1 + M} \frac{\partial d^*}{\partial t} - (1 - M) \frac{\partial d^*}{\partial x_1} + \frac{\partial v^*}{\partial x_2} = 0, \right. \qquad (44b) $$

$$ \left. \frac{1}{1 - M^2} \frac{\partial v^*}{\partial t} + M \frac{\partial v^*}{\partial x_1} + \frac{1}{2} \left( \frac{\partial s^*}{\partial x_2} + \frac{\partial d^*}{\partial x_2} \right) = 0, \right. \qquad (44c) $$
that are coupled to (43) by the continuity conditions :

$$s^*(0, x_2, t) = s(0, x_2, t), \quad d^*(0, x_2, t) = d(0, x_2, t), \quad v^*(0, x_2, t) = v(0, x_2, t).$$

(3) Apply to (44) the usual PML model.

**Remark 10** The well-posedness of the transmission problem obtained after step 2 is guaranteed by the following theorem, whose proof, omitted here, is similar to the one of theorem (8).

**Theorem 11** The following energy is conserved :

$$E_2(t) = \frac{1}{2} \int_{\mathbb{R}^{-}} [s^2 + d^2 + 2v^2] + \frac{1}{2} \int_{\mathbb{R}^{+}} \left[ \frac{s^{*2}}{1-M} + \frac{d^{*2}}{1+M} + 2 \frac{v^{*2}}{1-M^2} \right]. \quad (45)$$

**Comparison with the Hu's model.** The PML recently proposed by Hu are close to the one we have presented here. Both are based on the same change of variable, the difference lies in the fact that Hu applies this change of variable in the whole domain whereas we only apply it for $x_1 > 0$. Therefore he does not need to impose transmission condition at the interface. Then he replace formally the derivatives with respect to $x_1$ by the operator $D_{x_1}^\sigma$. Finally he use the inverse change of variable in the whole domain to come back to the original variables. This inverse change of variable makes the operator $D_{x_1}^\sigma$ slightly more complicated.

## 3.3 Numerical Results

In figure (27) we represent the propagation of an advective wave in an uniform mean flow with $M = 0.5$ in open domain. The domain of interest is $[-5; 5]^2$ and we impose an initial condition $p = \exp(-4(x_1^2+x_2^2)\ln(2))$. The width of the PML is $L = 1$ and the profile of the damping coefficient is $\sigma(x) = \sigma_0 x^2/D$. We have represented the level lines of the pressure at time $t = 1, 3, 8, 30$ and 50 We have repeated the same experiment but with an initial condition

<figure><img src="image_32.png" /><figcaption>Fig. 27. Propagation of an advective wave $t = 1, 3, 8, 30$ and 50</figcaption></figure>

$u = x_1 \exp(-(x_1^2+x_2^2)\ln(2))$ and $v = -x_2 \exp(-(x_1^2+x_2^2)\ln(2))$ designed to simulate the propagation of a vorticity wave. Figure 28 represent the level lines of the modulus of the velocity at time $t = 1, 3, 7, 11$ and 50.
<figure><img src="image_33.png" /><figcaption>Fig. 28. Propagation of a vorticity wave t= 1, 3, 7, 11 and 50</figcaption></figure>

A The case of a non-parallel mean-flow

The case of a non-parallel mean flow (with respect to the boundaries) is of interest for practical applications. Here we restrict ourselves to the case of the advective wave equation. Surprisingly, linearized Euler equations still arised open questions that are under study.

A.1 Construction of a vertical PML for the advective equation

A non-parallel mean flow is characterized by the Mach vector $\mathbf{M} = [M_1 ; M_2]^t$ and we set $M = \|\mathbf{M}\|$. The advective wave equation is then :

$$ \frac{\partial^2 p}{\partial t^2} + 2M_1 \frac{\partial^2 p}{\partial t \partial x_1} + 2M_2 \frac{\partial^2 p}{\partial t \partial x_2} + 2M_1 M_2 \frac{\partial^2 p}{\partial x_1 \partial x_2} - (1-M_1^2) \frac{\partial^2 p}{\partial x_1^2} (1-M_2^2) \frac{\partial^2 p}{\partial x_2^2} = 0. $$

The slowness curve associated to this equation is represented in figure A.1 (the set of unstable points for a vertical layer has been highlighted in red), it is an ellipse of center $C = (M_1/(1 - M^2), M_2/(1 - M^2))$. The construction of stabilized PML's is inspired by what has been done in the case of an horizontal mean flow. The novelty lies in the apparition of a new step (step (2.2))

(1) Write a transmission problem between $x_1 < 0$ and $x_1 > 0$.

(2.1) In the region $x_1 > 0$, get rid of the convective term by applying a change of variable :

$$ p^*(x_1, x_2, t) = p \left( x_1, x_2, t - \frac{M_1}{1 - M^2} x_1 - \frac{M_2}{1 - M^2} x_2 \right). $$

The new slowness curve (see figure A.2) still gives rise to instability if one applies classical PML. This is due to the presence of the cross-derivative with respect to $x_1$ and $x_2$.

(2.2) In the region $x_1 > 0$ we apply the second change of variable (which let invariant the axis $x_1 = 0$)

$$ \tilde{x}_1 = \alpha x_1 \text{ and } \tilde{x}_2 = x_2 + \beta x_1. $$
Writing the equations in the new variables, one observes that the cross derivatives disappear if $\alpha$ and $\beta$ are related by:

$$ \beta = -\frac{M_1 M_2}{1 - M_1^2} \alpha. $$

Choosing $\alpha = 1$, we see that the new function

$$ \tilde{p}(x_1, x_2, t) = p^* \left( x_1, x_2 - \frac{M_1 M_2}{1 - M_1^2} x_1, t \right) $$

satisfies the following equation, whose slowness curve is shown in Fig. A.3:

$$ \frac{1}{1 - M^2} \frac{\partial^2 \tilde{p}}{\partial t^2} - (1 - M_1^2) \frac{\partial^2 \tilde{p}}{\partial x_1^2} - (1 - \frac{M_2^2}{1 - M_1^2}) \frac{\partial^2 \tilde{p}}{\partial x_2^2} = 0. \quad (A.1) $$

(3) Apply the standard PML technique to (A.1).

<figure><img src="image_34.png" /><figcaption>Fig. A.1. The slowness curve associated to the advective waves</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. A.2. The slowness curve after the first change of variable</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. A.3. The slowness curve after the second change of variable</figcaption></figure>

## A.2 The corner problem

The solution we have adopted for corners is much less justified from a numerical point of view but appears to give satisfactory results (see Fig. A.4). Let us for instance suppose that the computational domain is the quarter-plane $\mathbb{R}^- \times \mathbb{R}^-$. In the horizontal and vertical layers, according to the previous section, we use the following models

$$ \left| \frac{1}{1 - M^2} \frac{\partial^2 \tilde{p}}{\partial t^2} - (1 - M_1^2)(D_{x_1}^\sigma)^2 \tilde{p} - (1 - \frac{M_2^2}{1 - M_1^2}) \frac{\partial^2 \tilde{p}}{\partial x_2^2} = 0 \quad (\text{vertical}) \right. \\ \left. \frac{1}{1 - M^2} \frac{\partial^2 \tilde{p}}{\partial t^2} - (1 - \frac{M_1^2}{1 - M_2^2}) \frac{\partial^2 \tilde{p}}{\partial x_1^2} - (1 - M_2^2)(D_{x_2}^\sigma)^2 \tilde{p} = 0 \quad (\text{horizontal}). \right. $$
In the corner we have chosen to use the following equation in which variables $x_1$ and $x_2$ are not treated in a symmetric way!

$$ \frac{1}{1 - M^2} \frac{\partial^2 \tilde{p}}{\partial t^2} - (1 - M_1^2)(D_{x_1}^\sigma)^2 \tilde{p} - (1 - \frac{M_2^2}{1 - M_1^2})(D_{x_2}^\sigma)^2 \tilde{p} = 0 $$

We refer to [23] and present in figure A.4 a numerical result obtained with this method.

<figure><img src="image_37.png" /><figcaption>Fig. A.4. The case of a non-parallel mean flow, t= 1 3 7 15 and 50</figcaption></figure>

## References

[1] J. P. Bérenger, A Perfectly Matched Layer for the Absorption of Electromagnetic Waves, J. of Comp. Phys. 114 (1994) 185-200.

[2] C. K. W. Tam, L. Auriault, F. Cambuli, Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains, J. Comput. Phys. 144 (1) (1998) 213-234.

[3] L. Zhao, A. C. Cangellaris, A General Approach for the Development of Unsplit-Field Time-Domain Implementations of Perfectly Matched Layers for FDTD Grid Truncation, IEEE Microwave and Guided Letters 6 (5).

[4] F. Collino, C. Tsogka, Application of the pml absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics 66 (1) (2001) 294-307.

[5] F. Q. Hu, On absorbing boundary conditions for linearized euler equations by a perfectly matched layer, J. Comp. Phys. (129) (1996) 201.

[6] J. P. Bérenger, Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. of Comp. Phys., 127 (1996) 363-379.

[7] S. Fauqueux, Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en régime transitoire, Ph.D. thesis, Université Paris IX (2003).

[8] F. Collino, P. Monk, The Perfectly Matched Layer in curvilinear coordinates, SIAM J. Scient. Comp. 164 (1998) 157-171.
[9] S. Abarbanel, D. Gottlieb, A mathematical analysis of the PML method, J. Comput. Phys. 134 (2) (1997) 357-363.

[10] E. Bécache, S. Fauqueux, P. Joly, Stability of perfectly matched layers, group velocities and anisotropic waves, J. Comput. Phys. 188 (2) (2003) 399-433.

[11] F. Collino, P. Monk, Optimizing the Perfectly Matched Layer, Comp. Meth. Appl. Mech. Engrg 164 (1998) 157-171.

[12] M. Lassas, E. Somersalo, On the existence and convergence of the solution of pml equations, Computing - Wien 60 (3) (1998) 229-242.

[13] T. Hohage, F. Schmidt, L. Zschiedrich, Solving time-harmonic scattering problems based on the pole condition. II. Convergence of the PML method, SIAM J. Math. Anal. 35 (3) (2003) 547-560 (electronic).

[14] E. Bécache, A.-S. Bonnet-Ben Dhia, G. Legendre, Perfectly matched layers for the convected Helmholtz equation, SIAM J. Numer. Anal. 42 (1) (2004) 409-433 (electronic).

[15] J. Diaz, P. Joly, An analysis of higher boundary conditions for the wave equation, to be published in SIAP.

[16] L. Cagniard, Reflection and refraction of progressive seismic waves, McGraw-Hill, 1962.

[17] A. T. de Hoop, The surface line source problem, Appl. Sci. Res. B 8 (1959) 349-356.

[18] A. T. de Hoop, P. van den Berg, F. Remis, Analytic time-domain performance analysis of absorbing boundary conditions and perfectly matched layers, in: Proc. IEEE Antennas and Propagation Society International Symposium, Vol. 4, 2001, pp. 502-505.

[19] J. Diaz, P. Joly, Stabilized perfectly matched layer for advective acoustics, in: The Sixth International Conference on Mathematical and Numerical Aspects of Wave Propagation (Waves 2003), 2003, pp. 115-119.

[20] T. Hagström, I. Nazarov, Absorbing layers and radiation boundary condition for jet flow simulations, Tech. Rep. Technical Report AIAA 2002-2606, AIAA (2002).

[21] S. Abarbanel, D. Gottlieb, J. S. Hesthaven, Well-posed perfectly matched layers for advective acoustics, J. Comput. Phys. 154 (2) (1999) 266-283.

[22] F. Q. Hu, A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables, J. Comp. Phys. 173 (2) (2001) 455-480.

[23] J. Diaz, Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. application au couplage fluide-structure et aux méthodes de couches parfaitement adaptées, Ph.D. thesis, Université Paris 6 (2005).
[24] E. Bécache, P. Joly, On the analysis of Bérenger's perfectly matched layers for Maxwell's equations., M2AN Math. Model. Numer. Anal. 36 (1) (2002) 87–119.

[25] J. S. Hesthaven, On the analysis and construction of perfectly matched layers for the linearized Euler equations, J. Comput. Phys. 142 (1) (1998) 129–147.

[26] J.-L. Lions, J. Métral, O. Vacus, Well-posed absorbing layer for hyperbolic problems, Numer. Math. 92 (3) (2002) 535–562.