# Perfectly matched layers for time-harmonic acoustics in the presence of a uniform flow

Eliane Bécache, Anne-Sophie Bonnet-Ben Dhia, Guillaume Legendre

► To cite this version:

Eliane Bécache, Anne-Sophie Bonnet-Ben Dhia, Guillaume Legendre. Perfectly matched layers for time-harmonic acoustics in the presence of a uniform flow. [Research Report] RR-5486, INRIA. 2005, pp.35. inria-00070521

HAL Id: inria-00070521

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

Submitted on 19 May 2006

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.
<figure><img src="image_1.png" /></figure>

INRIA

INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE

*Perfectly matched layers for time-harmonic acoustics*
*in the presence of a uniform flow*

Eliane Bécache — Anne-Sophie Bonnet-Ben Dhia — Guillaume Legendre

N° 5486

Fevrier 2005

THÈME 4

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>
# Perfectly matched layers for time-harmonic acoustics in the presence of a uniform flow

Eliane Bécache*, Anne-Sophie Bonnet-Ben Dhia†, Guillaume Legendre‡

Thème 4 — Simulation et optimisation
de systèmes complexes
Projet ONDES (POems)

Rapport de recherche n° 5486 — Fevrier 2005 — 35 pages

**Abstract:** This paper is devoted to the resolution of the time-harmonic linearized Galbrun's equation, which models, via a mixed Lagrangian-Eulerian representation, the propagation of acoustic and hydrodynamic perturbations in a given flow of a compressible fluid. We consider here the case of a uniform subsonic flow in an infinite, two-dimensional duct. Using a limiting amplitude process, we characterize the outgoing solution radiated by a compactly supported source. Then, we propose a Fredholm formulation with perfectly matched absorbing layers for approximating this outgoing solution. The convergence of the approximated solution to the exact one is proved, and error estimates with respect to the parameters of the absorbing layers are derived. Several significant numerical examples are included.

**Key-words:** aeroacoustics, Galbrun's equation, limiting absorption principle, Perfectly Matched Layers, acoustic waveguide, modal decomposition

<footnote>* Laboratoire POems, UMR 2706 CNRS/INRIA/ENSTA, INRIA, Domaine de Voluceau-Rocquencourt, BP 105, 78153 Le Chesnay cedex, France (eliane.becache@inria.fr).</footnote>

<footnote>† Laboratoire POems, UMR 2706 CNRS/INRIA/ENSTA, ENSTA, 32 boulevard Victor, 75739 Paris cedex 15, France (bonnet@ensta.fr).</footnote>

<footnote>‡ Laboratoire de mathématiques appliquées, UMR 7641 du CNRS, Université de Versailles Saint-Quentin-en-Yvelines, 45 avenue des États-Unis, 78035 Versailles, France (legendre@math.uvsq.fr).</footnote>
Couches parfaitement adaptées pour l'aéro-acoustique en régime
harmonique en présence d'un écoulement uniforme

Résumé : Ce papier est consacré à la résolution des équations de Galbrun linéarisées en régime har-
monique, qui modélisent, via une représentation mixte Lagrangienne-Eulérienne, la propagation de per-
turbations acoustiques et hydrodynamiques dans un écoulement donné de fluide compressible. Nous
considérons ici le cas d'un écoulement uniforme subsonique dans un conduit infini bidimensionnel. En
utilisant un procédé d'amplitude limite, nous caractérisons la solution sortante rayonnée par une source
à support compact. Nous proposons ensuite, pour approcher cette solution sortante, une formulation de
Fredholm avec des couches absorbantes parfaitement adaptées. Nous démontrons la convergence de la
solution approchée vers la solution exacte et donnons des estimations d'erreur, dépendant des paramètres
des couches absorbantes. Finalement, nous présentons plusieurs expériences numériques significatives.

Mots-clés : aéroacoustique, équations de Galbrun, principe d'absorption limite, Couches parfaitement adaptées, PML, guide d'onde acoustique, décomposition modale
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>The physical problem posed in an infinite waveguide</td><td>4</td></tr><tr><td>3</td><td>Well-posedness – the limiting absorption principle</td><td>5</td></tr><tr><td>3.1</td><td>The dissipative problem</td><td>5</td></tr><tr><td>3.2</td><td>Study of the dissipative problem</td><td>6</td></tr><tr><td>3.3</td><td>Convergence of the dissipative problem</td><td>8</td></tr><tr><td>3.3.1</td><td>Use of potentials</td><td>8</td></tr><tr><td>3.3.2</td><td>Limit and convergence of the acoustic problem</td><td>8</td></tr><tr><td>3.3.3</td><td>Limit and convergence of the hydrodynamic problem</td><td>10</td></tr><tr><td>3.3.4</td><td>Conclusion</td><td>11</td></tr><tr><td>3.3.5</td><td>Another characterization of ϕ<sub>h</sub></td><td>11</td></tr><tr><td>4</td><td>Setting of the problem with perfectly matched layers</td><td>14</td></tr><tr><td>4.1</td><td>The PML formulation</td><td>14</td></tr><tr><td>4.2</td><td>Well-posedness</td><td>15</td></tr><tr><td>5</td><td>Convergence results of PMLs for Galbrun's equation</td><td>18</td></tr><tr><td>5.1</td><td>Convergence results of PMLs for scalar problems</td><td>18</td></tr><tr><td>5.1.1</td><td>Problem A</td><td>19</td></tr><tr><td>5.1.2</td><td>Problem B</td><td>21</td></tr><tr><td>5.1.3</td><td>Problem C</td><td>21</td></tr><tr><td>5.2</td><td>Well-posedness and convergence analysis for ϕ<sub>h</sub><sup>L</sup></td><td>22</td></tr><tr><td>5.3</td><td>Well-posedness and convergence analysis for ϕ<sub>a</sub><sup>L</sup></td><td>24</td></tr><tr><td>5.4</td><td>Conclusion</td><td>26</td></tr><tr><td>5.5</td><td>Remark on the use of PMLs without regularization</td><td>26</td></tr><tr><td>6</td><td>Numerical applications</td><td>27</td></tr><tr><td>6.1</td><td>Mode propagation in a rigid duct</td><td>27</td></tr><tr><td>6.1.1</td><td>Acoustic and vortical modes: some definitions</td><td>27</td></tr><tr><td>6.1.2</td><td>Description of the simulations</td><td>28</td></tr><tr><td>6.1.3</td><td>Numerical results for acoustic modes</td><td>28</td></tr><tr><td>6.1.4</td><td>Numerical results for vortical modes</td><td>28</td></tr><tr><td>6.2</td><td>Radiation of compactly supported sources</td><td>31</td></tr><tr><td>6.2.1</td><td>Acoustic source</td><td>31</td></tr><tr><td>6.2.2</td><td>Rotational source</td><td>32</td></tr><tr><td>6.3</td><td>A concluding remark on the non-regularized PML problem</td><td>32</td></tr><tr><td>A</td><td>Study of a second order transport equation</td><td>32</td></tr></table>

# 1 Introduction

Numerical simulations of wave phenomena in unbounded domains have always been concerned with artificial boundary conditions. The computational domain being necessarily finite, one has to impose boundary conditions such that the outgoing waves are not reflected. In computational aeroacoustics,
problems involve wave radiation and convection, adding to the complexity of defining an adequate condition. Many studies have been conducted in this field over the last thirty years, and various methods have been developed. Some are based on the localization of the Dirichlet-to-Neumann map [12], or on asymptotic expansions of the far-field solution [2], or on the use of characteristics [25]. Alternatively, buffer or sponge layer techniques can also be used, in which outgoing waves are damped through artificial dissipation [13] or grid stretching and filtering [9] in a region adjacent to the artificial boundary.

In this work, we make use of absorbing layers in the manner of perfectly matched layers (PML), which were introduced by Bérenger [5] in the context of computational electromagnetics. A major difference between the PML and other buffer zones is that they are designed in such a way that, at any angle of incidence, waves are transmitted with (theoretically) no reflection. While PMLs have already been used in aeroacoustics (see, for instance, [19, 18, 24, 1, 20]), the applications were concerned with the linearized Euler equations solved in the time domain. In a previous paper [3], we dealt with the PML for the convected Helmholtz equation (that is, in the frequency domain) in a waveguide and proved the well-posedness and convergence of the method. However, that model being scalar, acoustic waves were the only ones taken into account, whereas the main difficulty when applying the PML in aeroacoustics lies in the appropriate treatment in the layers of the entropy and vorticity waves, which are both convected downstream of the mean flow.

The model used in this paper is the one introduced by Galbrun [14, 23], which is based on a Lagrangian-Eulerian description of the perturbations, in the sense that the Lagrangian perturbations of the quantities are expressed in terms of Eulerian variables with respect to the mean flow, and it assumes small perturbations of an isentropic flow of a perfect fluid. The model consists of a linear partial differential equation involving the Lagrangian displacement perturbation. This equation, being of second order in time and space, is amenable to variational methods, which is not the case for the hyperbolic system of linearized conservation laws usually solved in aeroacoustics simulations. However, its resolution by standard (i.e., nodal) finite element methods is subject to difficulties quite similar to those observed for Maxwell's equations in electromagnetism [7]. We previously proposed a regularized formulation of the Galbrun's equation in a uniform flow and time-harmonic dependence [8] that allowed the use of nodal finite elements for the discretization of the problem.

In the present article, we discuss the use of the PML method to solve the radiation problem of a perturbation source placed in a two-dimensional guide in the presence of a uniform mean flow. The outline is the following. The problem to be solved is introduced in section 2. In section 3, we characterize through a limiting absorption principle the unique solution which represents the periodic regime. The problem with PMLs is posed and analyzed in section 4 and its convergence is subsequently proved in section 5 via the combined use of vector potentials and scalar modal analysis. Finally, numerical applications are presented in section 6.

## 2 The physical problem posed in an infinite waveguide

We consider a displacement formulation for the propagation of acoustic waves in an infinite, rigid duct and in the presence of a uniform mean flow of subsonic speed $v_0$. A time-harmonic dependence of the form $\exp(-i\omega t)$, $\omega > 0$ being the pulsation, is assumed throughout the paper. The problem is then modeled by the following equation and boundary condition

$$ (1) \qquad \mathrm{D}^2\boldsymbol{u} - \nabla (\mathrm{div} \, \boldsymbol{u}) = \boldsymbol{f} \quad \text{in } \Omega, $$

$$ (2) \qquad \boldsymbol{u} \cdot \boldsymbol{n} = 0 \quad \text{on } \partial\Omega, $$

where $\Omega$ and $\partial\Omega$ denote, respectively, the infinite duct of height $l$ and its rigid walls (i.e., $\Omega = \{(x_1, x_2); x_1 \in \mathbb{R} \text{ and } 0 < x_2 < l\}$). In equation (1), the letter D stands for the material derivative in the uniform flow
and the time-harmonic dependence, given by

$$D\boldsymbol{u} = -ik\boldsymbol{u} + M \frac{\partial\boldsymbol{u}}{\partial x_1},$$

where $k = \frac{\omega}{c_0}$ is the wave number, and $M = \frac{v_0}{c_0}$ is the Mach number ($0 < M < 1$). Thus, in extended form, equation (1) reads

$$-k^2 \boldsymbol{u} - 2ikM \frac{\partial \boldsymbol{u}}{\partial x_1} + M^2 \frac{\partial^2 \boldsymbol{u}}{\partial x_1^2} - \nabla (\operatorname{div} \boldsymbol{u}) = \boldsymbol{f} \text{ in } \Omega.$$

An additional hypothesis is made on the compactly supported source $f$, which is assumed to admit the Helmholtz decomposition

$$f = \nabla g_a + \operatorname{curl} g_h,$$

where $g_a$ and $g_h$ are also compactly supported. From a physical point of view, the source term $f$ is meant to contain an “acoustic” part $g_a$, which generates irrotational perturbations (i.e., pressure fluctuations), and a vortical part $g_h$, which creates hydrodynamic perturbations. Note here that $\operatorname{curl} f = \frac{\partial f}{\partial x_2} e_1 - \frac{\partial f}{\partial x_1} e_2$, where $e_1$ and $e_2$ are the vectors of the canonical basis of $\mathbb{R}^2$, is the vectorial form of the curl operator when applied to scalar functions. We denote by $\operatorname{curl} v = \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2}$ the dual form of this operator when applied to vector fields.

The source $f$ is also assumed to belong to the space $H(\operatorname{curl}; \Omega)$. Of course, this assumption implies some regularity on $g_a$ and $g_h$. In our case, sufficient conditions of regularity are

$$g_a \in H^1(\Omega) \text{ and } g_h \in H^2(\Omega).$$

The problem (1)-(2) admits an infinite number of solutions as long as an additional condition at infinity is not given. We are interested in the unique solution associated with the time-harmonic problem. In the next section, we characterize this solution through the study of a dissipative problem.

## 3 Well-posedness – the limiting absorption principle

The problem described in the previous section is not completely defined before the behavior of the solution at infinity is prescribed and the notion of an outgoing solution of this problem is specified. This will be done in a very natural manner, using the limiting absorption principle [11]. The underlying idea consists of considering the physical solution as the limit solution of a dissipative problem, which is well-posed due to the absorption.

### 3.1 The dissipative problem

A dissipative problem associated with (1)-(2) is readily obtained by replacing the real wave number $k$ by a complex number $k_\varepsilon$ such that

$$k_\varepsilon = k + i\varepsilon,$$

where $\varepsilon$ is a positive real number. The physical case then becomes the limiting case in which $\varepsilon$ tends to 0. In what follows, we prove that the unique solution of “finite energy” (i.e., which belongs to the space $H^1(\Omega)^2$) of the dissipative problem converges as $\varepsilon$ tends to zero in the $H^1_{\text{loc}}(\Omega)^2$ sense to a limit, that will be called the “outgoing” solution of (1)-(2).
## 3.2 Study of the dissipative problem

We now consider the dissipative problem. We seek a function $u^\epsilon$ in $H^1(\Omega)^2$ which solves the following equation with boundary condition

$$ (3) \qquad D_{\epsilon}^{2} u^{\epsilon} - \nabla (\operatorname{div} u^{\epsilon}) = f \quad \text{in } \Omega, $$

$$ (4) \qquad u^{\epsilon} \cdot n = 0 \text{ on } \partial\Omega, $$

where $D_\epsilon = -ik_\epsilon + M \frac{\partial}{\partial x_1}$. To be able to prove the well-posedness of this problem, it must be regularized, as proposed in [8]. To this end, we introduce the function $\psi^\epsilon = \operatorname{curl} u^\epsilon$, belonging to $L^2(\Omega)$, which is a solution of the following ordinary differential equation with constant coefficients (obtained by taking the curl of equation (3))

$$ (5) \qquad D_{\epsilon}^{2} \psi^{\epsilon} = \operatorname{curl} f \quad \text{in } \Omega. $$

First, we state a preliminary result.

**Lemma 1** Equation (5) has a unique solution $\psi^\epsilon$ in $L^2(\Omega)$. This solution vanishes upstream of the support of the source $f$.

*Proof.* Introducing the causal Green's function of the differential operator $D_\epsilon^2$

$$ G_{\epsilon}(x_1) = \frac{x_1}{M^2} H(x_1) e^{i \frac{k_{\epsilon}}{M} x_1}, \quad \forall x_1 \in \mathbb{R}, $$

where $H$ denotes the Heaviside function, one can derive the following particular solution of (5)

$$ \psi^{\epsilon}(x_1, x_2) = G_{\epsilon} * \operatorname{curl} f(\cdot, x_2)(x_1) = \frac{1}{M^2} \int_{-\infty}^{x_1} (x_1 - z) e^{i \frac{k_{\epsilon}}{M} (x_1 - z)} \operatorname{curl} f(z, x_2) dz. $$

One can easily verify that this function vanishes upstream of the source and belongs to $L^2(\Omega)$ (see the appendix). Uniqueness of the solution follows from the fact that the solutions of the homogeneous equation $D_\epsilon^2 \psi^\epsilon = 0$ are of the following form

$$ \psi^{\epsilon}(x_1, x_2) = (a_{\epsilon}(x_2) + x_1 b_{\epsilon}(x_2)) e^{i \frac{k_{\epsilon}}{M} x_1}, \quad \forall (x_1, x_2) \in \Omega, $$

and therefore do not belong to $L^2(\Omega)$, except for the trivial solution $\psi^\epsilon \equiv 0$. $\square$

Now, if $u^\epsilon$ is a solution of (3)-(4), it clearly satisfies the so-called regularized or augmented problem

$$ D_{\epsilon}^{2} u^{\epsilon} - \nabla (\operatorname{div} u^{\epsilon}) + \operatorname{curl} (\operatorname{curl} u^{\epsilon} - \psi^{\epsilon}) = f \quad \text{in } \Omega, $$

$$ (6) \qquad u^{\epsilon} \cdot n = 0 \text{ on } \partial\Omega, $$

$$ \operatorname{curl} u^{\epsilon} = \psi^{\epsilon} \text{ on } \partial\Omega. $$

Setting $V(\Omega) = \{v \in H^1(\Omega)^2 \mid v \cdot n = 0 \text{ on } \partial\Omega\}$, a variational formulation of this last problem reads: find $u^\epsilon \in V(\Omega)$ such that

$$ (7) \qquad a_{\Omega}(k_{\epsilon}; u^{\epsilon}, v) = \int_{\Omega} (f \cdot \bar{v} + \psi^{\epsilon}(\operatorname{curl} \bar{v})) \, d\mathbf{x}, \quad \forall v \in V(\Omega), $$

where the sesquilinear form $a_\Omega(k_\epsilon; \cdot, \cdot)$ is defined by

$$ a_{\Omega}(k_{\epsilon}; \mathbf{u}, \mathbf{v}) = \int_{\Omega} \left( (\operatorname{div} \mathbf{u})(\operatorname{div} \bar{\mathbf{v}}) + (\operatorname{curl} \mathbf{u})(\operatorname{curl} \bar{\mathbf{v}}) - M^2 \frac{\partial \mathbf{u}}{\partial x_1} \cdot \frac{\partial \bar{\mathbf{v}}}{\partial x_1} \right) d\mathbf{x} \\ + \int_{\Omega} \left( -k_{\epsilon}^2 \mathbf{u} \cdot \bar{\mathbf{v}} - 2ik_{\epsilon}M \frac{\partial \mathbf{u}}{\partial x_1} \cdot \bar{\mathbf{v}} \right) d\mathbf{x}. $$
**Theorem 1** The variational problem (7) is well-posed.

*Proof.* Integrating by parts gives

$$ \int_{\Omega} \frac{\partial \mathbf{u}}{\partial x_1} \cdot \overline{\mathbf{u}} \, d\mathbf{x} = - \int_{\Omega} \mathbf{u} \cdot \frac{\partial \overline{\mathbf{u}}}{\partial x_1} \, d\mathbf{x} = - \overline{\int_{\Omega} \frac{\partial \mathbf{u}}{\partial x_1} \cdot \overline{\mathbf{u}} \, d\mathbf{x}}, \quad \forall \mathbf{u} \in H^1(\Omega)^2, $$

hence $\int_{\Omega} \frac{\partial \mathbf{u}}{\partial x_1} \cdot \overline{\mathbf{u}} \, d\mathbf{x} \in i\mathbb{R}$. We then have

$$ \operatorname{Im} \left( -\frac{1}{k_{\varepsilon}} a_{\Omega}(k_{\varepsilon}; \mathbf{u}, \mathbf{u}) \right) = \int_{\Omega} \frac{\operatorname{Im}(k_{\varepsilon})}{|k_{\varepsilon}|^2} \left( |\operatorname{div} \mathbf{u}|^2 + |\operatorname{curl} \mathbf{u}|^2 - M^2 \left| \frac{\partial \mathbf{u}}{\partial x_1} \right|^2 \right) d\mathbf{x} \\ + \int_{\Omega} \operatorname{Im}(k_{\varepsilon}) |\mathbf{u}|^2 d\mathbf{x}. $$

Since $M^2 < 1$ and $\operatorname{Im}(k_{\varepsilon}) > 0$, the sesquilinear form $a_{\Omega}(k_{\varepsilon}; \cdot, \cdot)$ is coercive on $H^1(\Omega)^2$, due to theorem 4.1 of [10]. It is also clear that this form is continuous on the same space.

Moreover, estimate 49 (see the appendix) allows one to establish the continuity of the antilinear form simply by using the Cauchy-Schwarz inequality. The well-posedness of problem (7) is then a consequence of the Lax-Milgram lemma. □

By construction, every solution of (3)-(4) belonging to $H^1(\Omega)$ is a solution of (7). The converse statement results from the following theorem.

**Theorem 2** The solution $\mathbf{u}^{\varepsilon}$ of problem (7) is such that $\operatorname{curl} \mathbf{u}^{\varepsilon} = \psi^{\varepsilon}$.

*Proof.* Taking as test functions $\mathbf{v} = \operatorname{curl} \varphi$ with $\varphi \in \{\phi \in H^3(\Omega) \mid \phi|_{\partial\Omega} = 0\}$, we obtain, after some integrations by parts, and the use of boundary conditions of problem (7), the following orthogonality relation

$$ \int_{\Omega} (\operatorname{curl} \mathbf{u}^{\varepsilon} - \psi^{\varepsilon}) (\mathcal{H}_{k_{\varepsilon}, M} \overline{\varphi}) \, d\mathbf{x} = 0. $$

Here $\mathcal{H}_{k_{\varepsilon}, M}$ denotes the operator $D_{\varepsilon}^2 - \Delta$. By a density result (theorem 1.6.2 of [17]), this relation holds for any function $\varphi$ of $D(\mathcal{H}_{k_{\varepsilon}, M}) = H^2(\Omega) \cap H_0^1(\Omega)$. To conclude that $\operatorname{curl} \mathbf{u}^{\varepsilon} = \psi^{\varepsilon}$ in $L^2(\Omega)$, it suffices to show that $\mathcal{H}_{k_{\varepsilon}, M}$ is surjective from $D(\mathcal{H}_{k_{\varepsilon}, M})$ to $L^2(\Omega)$. For all $\varphi$ in $D(\mathcal{H}_{k_{\varepsilon}, M})$, we have

$$ (\mathcal{H}_{k_{\varepsilon}, M} \varphi, \varphi) = \int_{\Omega} \left( -k_{\varepsilon}^2 |\varphi|^2 + 2ik_{\varepsilon}M \frac{\partial \varphi}{\partial x_1} \overline{\varphi} + |\nabla \varphi|^2 - M^2 \left| \frac{\partial \varphi}{\partial x_1} \right|^2 \right) d\mathbf{x}. $$

As in the proof of the preceding theorem, we have $\int_{\Omega} \frac{\partial \varphi}{\partial x_1} \overline{\varphi} \, d\mathbf{x} \in i\mathbb{R}$, and we deduce that

$$ \operatorname{Im} \left( -\frac{1}{k_{\varepsilon}} (\mathcal{H}_{k_{\varepsilon}, M} \varphi, \varphi) \right) = \int_{\Omega} \left( \frac{\operatorname{Im}(k_{\varepsilon})}{|k_{\varepsilon}|^2} \left( |\nabla \varphi|^2 - M^2 \left| \frac{\partial \varphi}{\partial x_1} \right|^2 \right) + \operatorname{Im}(k_{\varepsilon}) |\varphi|^2 \right) d\mathbf{x}. $$

Again, the surjectivity of the operator is a consequence of the Lax-Milgram lemma applied to the sesquilinear form $(\mathcal{H}_{k_{\varepsilon}, M}, \cdot, \cdot)$. □

**Corollary 1** Problem (3)-(4) has a unique solution in $H^1(\Omega)^2$ which is the solution of problem (7).

*Proof.* We choose $\mathbf{v} \in \mathcal{D}(\Omega)^2 \subset V(\Omega)$ in the variational formulation (7). Using integration by parts and the previous theorem, we obtain that the unique solution $\mathbf{u}^{\varepsilon}$ of (7) verifies equation (3) in the distributional sense. The boundary condition (4) is also satisfied since $\mathbf{u}^{\varepsilon} \in V(\Omega)$. □
## 3.3 Convergence of the dissipative problem

We will prove in this subsection that, if $k$ is not a cut-off wave number for acoustic modes, the solution $u_\varepsilon$ of problem (7) converges in $H_{\text{loc}}^1(\Omega)^2$ to a limit $u$ as $\varepsilon$ tends to 0. This limit is clearly a solution of (1)-(2) and, contrary to $u_\varepsilon$, does not belong to $H^1(\Omega)^2$, since it does not decrease a infinity. The proof of convergence is based on a Helmholtz decomposition of the field $u_\varepsilon$ and the use of convergence results for scalar problems.

### 3.3.1 Use of potentials

Let us consider the following problems: find $\varphi_a^\varepsilon \in H^1(\Omega)$ such that

$$ (8) \qquad \begin{array}{l} D_\varepsilon^2 \varphi_a^\varepsilon - \Delta \varphi_a^\varepsilon = g_a \quad \text{in } \Omega, \\ \frac{\partial \varphi_a^\varepsilon}{\partial \boldsymbol{n}} = 0 \quad \text{on } \partial \Omega, \end{array} $$

and: find $\varphi_h^\varepsilon \in L^2(\Omega)$ such that

$$ (9) \qquad D_\varepsilon^2 \varphi_h^\varepsilon = g_h \quad \text{in } \Omega. $$

Both problems are well-posed. Indeed, problem (8) has been studied in [6] (theorem 1) and we obtain, using the regularity of $g_a$ and domain $\Omega$, that its solution $\varphi_a^\varepsilon$ belongs to the space $H^2(\Omega)$. Problem (9) was dealt with in Lemma 1, the regularity of $g_h$ implying that $\varphi_h^\varepsilon$ is in $H^2(\Omega)$ (see the appendix). It then clearly follows that the function $\nabla \varphi_a^\varepsilon + \text{curl} \varphi_h^\varepsilon$ is a solution of (3)-(4) (or equivalently of problem (7)) since

$$ \nabla (D_\varepsilon^2 \varphi_a^\varepsilon - \Delta \varphi_a^\varepsilon) + \text{curl} (D_\varepsilon^2 \varphi_h^\varepsilon) = \nabla g_a + \text{curl} g_h = f \quad \text{in } \Omega, $$

and

$$ \frac{\partial \varphi_a^\varepsilon}{\partial \boldsymbol{n}} + \text{curl} \varphi_h^\varepsilon \cdot \boldsymbol{n} = 0 \quad \text{on } \partial \Omega. $$

Indeed, since the function $g_h$ is compactly supported, $\varphi_h^\varepsilon$ vanishes on the boundary $\partial\Omega$, which implies that $\text{curl} \varphi_h^\varepsilon \cdot \boldsymbol{n} = 0$ on $\partial\Omega$. Hence, the uniqueness of the solution to (7) implies that

$$ (10) \qquad \boldsymbol{u}^\varepsilon = \nabla \varphi_a^\varepsilon + \text{curl} \varphi_h^\varepsilon. $$

We now prove the convergence of the respective solutions of problems (8) and (9) as $\varepsilon$ tends to zero.

### 3.3.2 Limit and convergence of the acoustic problem

In order to get the limit in $H_{\text{loc}}^2(\Omega)$ of $\varphi_a^\varepsilon$ as $\varepsilon$ tends to zero, we use some theoretical results previously established in [6] for scalar problems of the same type. First, problem (8) is equivalently set in a bounded domain $\Omega_b$ (see figure 1), which contains the supports of $g_a$ and $g_h$ and is situated in between two vertical boundaries $\Sigma_\pm$, respectively located at $x_1 = x_\pm$. To this end, we make use of the Dirichlet-to-Neumann (DtN) operators $T_\pm^{N,\varepsilon}$, defined as follows (the superscript $N$ refers here to the Neumann boundary condition in problem (14))

$$ \begin{array}{rcl} T_{\pm}^{N,\varepsilon} : & H^{1/2}(\Sigma_{\pm}) & \rightarrow H^{-1/2}(\Sigma_{\pm}) \\ & \phi & \mapsto \mp i \sum_{n=0}^{+\infty} \beta_n^{\varepsilon \pm} (\phi, C_n)_{L^2(\Sigma_{\pm})} C_n(x_2), \end{array} $$
<figure><img src="image_4.png" /><figcaption>Figure 1: The partition of the domain $\Omega$.</figcaption></figure>

where

$$ (11) \qquad \beta_n^{\varepsilon \pm} = \frac{-k_{\varepsilon} M \pm \sqrt{k_{\varepsilon}^2 - \frac{n^2 \pi^2}{l^2} (1 - M^2)}}{1 - M^2}, $$

with the following definition of the complex square root

$$ (12) \qquad \sqrt{z} = \sqrt{|z|} e^{i \frac{\arg(z)}{2}}, \quad 0 \le \arg(z) < 2\pi, $$

and where

$$ (13) \qquad C_0(x_2) = \sqrt{\frac{1}{l}} \quad \text{and} \quad C_n(x_2) = \sqrt{\frac{2}{l}} \cos\left(\frac{n\pi}{l}x_2\right), \quad \forall n \in \mathbb{N}^*. $$

An equivalent formulation of problem (8) is then: find $\varphi_a^\varepsilon \in H^1(\Omega_b)$ such that

$$ (14) \qquad \begin{aligned} D_{\varepsilon}^2 \varphi_a^\varepsilon - \Delta \varphi_a^\varepsilon &= g_a \quad \text{in } \Omega_b, \\ \frac{\partial \varphi_a^\varepsilon}{\partial \boldsymbol{n}} &= 0 \quad \text{on } \partial\Omega \cap \partial\Omega_b, \\ \frac{\partial \varphi_a^\varepsilon}{\partial \boldsymbol{n}} &= -T_{\pm}^{N,\varepsilon} \varphi_a^\varepsilon \quad \text{on } \Sigma_{\pm}. \end{aligned} $$

We now can formally derive a limit problem for (14). Observe that, because of definition (12), one has

$$ \lim_{\varepsilon \to 0} \sqrt{k_{\varepsilon}^2 - \frac{n^2 \pi^2}{l^2} (1 - M^2)} = \begin{cases} \sqrt{k^2 - \frac{n^2 \pi^2}{l^2} (1 - M^2)} & \in \mathbb{R}_+ \text{ if } k \ge \frac{n\pi}{l} \sqrt{1 - M^2}, \\ i\sqrt{\frac{n^2 \pi^2}{l^2} (1 - M^2) - k^2} & \in i\mathbb{R}_+ \text{ if } k < \frac{n\pi}{l} \sqrt{1 - M^2}. \end{cases} $$

Hence, the respective limits $\beta_n^{\pm}, \forall n \in \mathbb{N}$, of axial wave numbers $\beta_n^{\pm}$ are simply defined by

$$ (15) \qquad \beta_n^{\pm} = \begin{cases} \frac{-kM \pm \sqrt{k^2 - \frac{n^2\pi^2}{l^2}(1-M^2)}}{1-M^2} & \text{if } k \ge \frac{n\pi}{l}\sqrt{1-M^2}, \\ \frac{-kM \pm i\sqrt{\frac{n^2\pi^2}{l^2}(1-M^2) - k^2}}{1-M^2} & \text{if } k < \frac{n\pi}{l}\sqrt{1-M^2}. \end{cases} $$

The limit problem to be considered is then: find $\varphi_a \in H^1(\Omega_b)$ such that

$$ (16) \qquad \begin{aligned} D^2\varphi_a - \Delta\varphi_a &= g_a \quad \text{in } \Omega_b, \\ \frac{\partial\varphi_a}{\partial \boldsymbol{n}} &= 0 \quad \text{on } \partial\Omega \cap \partial\Omega_b, \\ \frac{\partial\varphi_a}{\partial \boldsymbol{n}} &= -T_{\pm}^N \varphi_a \quad \text{on } \Sigma_{\pm}, \end{aligned} $$
with the following obvious definition for the DtN operators $T_{\pm}^N$

$$ T_{\pm}^N : H^{1/2}(\Sigma_{\pm}) \rightarrow H^{-1/2}(\Sigma_{\pm}) \\ \phi \mapsto \mp i \sum_{n=0}^{+\infty} \beta_n^{\pm} (\phi, C_n)_{L^2(\Sigma_{\pm})} C_n(x_2), $$

**Theorem 3** Problem (16) is well-posed, except if $k = k_n$, for $n \in \mathbb{N}$, with $k_n = \sqrt{1 - M^2 \frac{n\pi}{l}}$.

A proof of this theorem is available in [3] (theorem 2.2). We are now in a position to prove the convergence for problem (8). The scalars $k_n, \forall n \in \mathbb{N}$, are the cut-off wave numbers of the modes.

**Theorem 4** If $k \neq k_n, \forall n \in \mathbb{N}$, the solution $\varphi_a^\varepsilon$ of problem (8) tends in $H^2(\Omega_b)$ to $\varphi_a$ as $\varepsilon$ tends to zero, $\varphi_a$ being the solution of (16).

*Proof.* Theorem 4 of [6] gives the convergence of $\varphi_a^\varepsilon$ to $\varphi_a$ in $H^1(\Omega_b)$. We then have

$$ (1 - M^2) \frac{\partial^2 \varphi_a^\varepsilon}{\partial x_1^2} + \frac{\partial^2 \varphi_a^\varepsilon}{\partial x_2^2} \underset{\varepsilon \to 0}{\longrightarrow} (1 - M^2) \frac{\partial^2 \varphi_a}{\partial x_1^2} + \frac{\partial^2 \varphi_a}{\partial x_2^2} \text{ in } L^2(\Omega_b). $$

The domain $\Omega_b$ being convex, we deduce $\varphi_a^\varepsilon \underset{\varepsilon \to 0}{\longrightarrow} \varphi_a$ in $H^2(\Omega_b)$ (see [16]). □

### 3.3.3 Limit and convergence of the hydrodynamic problem

The solution of problem (9) is explicitely given by the convolution product $\varphi_h^\varepsilon(x_1, x_2) = G_\varepsilon * g_h(., x_2)(x_1)$, where the kernel $G_\varepsilon$ denotes the causal Green's function of the differential operator $D_\varepsilon^2$. Introducing

$$ G(x_1) = \frac{x_1}{M^2} H(x_1) e^{i \frac{k}{M} x_1} $$

as the formal limit of $G_\varepsilon$ as $\varepsilon$ tends to 0, one can show that $G_\varepsilon$ converges to $G$ in $L_{\text{loc}}^2(\mathbb{R})$. Let $\varphi_h(x_1, x_2) = G * g_h(., x_2)(x_1)$ be a solution of the limit problem: find $\varphi_h \in L_{\text{loc}}^2(\Omega)$ such that

$$ (17) \qquad D^2 \varphi_h = g_h \text{ in } \Omega. $$

One has

$$ |\varphi_h^\varepsilon - \varphi_h| = |(G_\varepsilon - G) * g_h(., x_2)|, \forall x_2 \in [0, l], $$

and, using the Cauchy-Schwarz inequality,

$$ \|\varphi_h^\varepsilon - \varphi_h\|_{L^2(\Omega_b)} \le \left( \int_{x_-}^{x_+} |G_\varepsilon(z) - G(z)|^2 dz \right)^{1/2} \|g_h\|_{L^2(\Omega)}. $$

Using the convergence of $G_\varepsilon$ to $G$ in $L_{\text{loc}}^2(\Omega)$, we obtain that $\varphi_h^\varepsilon$ converges to $\varphi_h$ in $L^2(\Omega_b)$ as $\varepsilon$ tends to 0. Since $g_h \in H^2(\Omega)$, we obtain as well, using classical properties of the convolution of distributions,

$$ \|\varphi_h^\varepsilon - \varphi_h\|_{H^2(\Omega_b)} \le \left( \int_{x_-}^{x_+} |G_\varepsilon(z) - G(z)|^2 dz \right)^{1/2} \|g_h\|_{H^2(\Omega)}. $$

and deduce the following result.

**Theorem 5** The solution $\varphi_h^\varepsilon$ of (9) tends to $\varphi_h$ in $H^2(\Omega_b)$ as $\varepsilon$ tends to zero.
### 3.3.4 Conclusion

We finally infer from theorems 4 and 5 the following result.

**Theorem 6** If $k$ is not a cut-off wave number, the solution $u^\epsilon$ of problem (3)-(4) tends to $u = \nabla \varphi_a + \text{curl} \varphi_h$ in $H^1(\Omega_b)^2$ as $\epsilon$ tends to zero, where $\varphi_a$ is the unique solution of (16) and $\varphi_h(x_1, x_2) = G * g_h(., x_2)(x_1)$.

The potential $\varphi_a$ can be extended via its modal expansion to the whole domain $\Omega$. The field $u$ is therefore also defined in the whole duct $\Omega$, where it obviously satisfies the equations (1) and (2). This field $u$ will be refered to in the sequel as the outgoing solution of (1) and (2).

**Corollary 2** The function $\psi^\epsilon = \text{curl } u^\epsilon$ tends to $\text{curl } u$ in $L^2(\Omega_b)$. We set $\psi = \text{curl } u$.

In the remainder of the paper, we assume that $k$ is not a cut-off wave number, i.e.,

$$ (18) \qquad k \neq k_n, \forall n \in \mathbb{N}. $$

### 3.3.5 Another characterization of $\varphi_h$

We note that using the Helmholtz decomposition (10) of $u^\epsilon$ in the regularized problem (6) would lead, as $\text{curl} (\text{curl}) = -\Delta$, to the following problem for the hydrodynamic potential: find $\varphi_h^\epsilon \in H^1(\Omega)$ such that

$$ (19) \qquad \begin{aligned} D_\epsilon^2 \varphi_h^\epsilon - \Delta \varphi_h^\epsilon &= g_h + \psi^\epsilon \quad \text{in } \Omega, \\ \varphi_h^\epsilon &= 0 \quad \text{on } \partial\Omega, \end{aligned} $$

in place of (9). Nevertheless, problems (9) and (19) are equivalent, since $-\Delta \varphi_h^\epsilon = \psi^\epsilon$, but using the latter to find the limit of $\varphi_h^\epsilon$ as $\epsilon$ tends to zero is less natural and more delicate than what has been proposed in subsection 3.3.3. Still, we detail this alternative approach here, since it provides another characterization of the potential $\varphi_h$, which will prove to be useful in what follows.

Notice that problem (19) is very similar to (8), except for the homogeneous Dirichlet boundary condition, which replaces the homogeneous Neumann boundary condition of (8), and for the right hand side term, which does not have compact support. To prove the convergence of this problem, we define a similar problem with compactly supported data, which then fits into the previous framework.

We first introduce a function $\psi^{\epsilon, \infty}$, which is the sum of two functions with separated variables and coincides with $\psi^\epsilon$ downstream of the support of $\text{curl } f$. More precisely, if we set

$$ d_- = \min_{x_2 \in [0, l]} \{x_1 \in \mathbb{R} \mid (x_1, x_2) \in \text{supp}(\text{curl } f)\} $$

and

$$ d_+ = \max_{x_2 \in [0, l]} \{x_1 \in \mathbb{R} \mid (x_1, x_2) \in \text{supp}(\text{curl } f)\}, $$

the function $\psi^{\epsilon, \infty}$ is such that

$$ \psi^{\epsilon, \infty}(x_1, x_2) = (a_\epsilon(x_2) + x_1 b_\epsilon(x_2)) e^{i \frac{k_\epsilon}{M} x_1}, \forall (x_1, x_2) \in \Omega, $$

and

$$ \psi^{\epsilon, \infty}(x_1, x_2) = \psi^{\epsilon}(x_1, x_2), \forall (x_1, x_2) \in ]d_+, +\infty[ \times [0, l]. $$

Notice that $\psi^{\epsilon, \infty}$ does not generally vanish upstream of the support of the source term $g_a$, contrary to $\psi^\epsilon$.
Taking advantage of the particular form of $\psi^{\epsilon,\infty}$, one can explicitly determine a function $\zeta^{\epsilon}$, which is a solution of the following problem: find $\zeta^{\epsilon} \in H^1(\Omega)$ such that

$$D_{\epsilon}^2 \zeta^{\epsilon} - \Delta \zeta^{\epsilon} = \psi^{\epsilon,\infty} \text{ in } \Omega,$$
$$\zeta^{\epsilon} = 0 \text{ on } \partial \Omega,$$

which we seek of the same form, i.e.

$$\zeta^{\epsilon}(x_1, x_2) = (A_{\epsilon}(x_2) + x_1 B_{\epsilon}(x_2)) e^{i \frac{k_{\epsilon}}{M} x_1}, \quad \forall (x_1, x_2) \in \Omega.$$

This leads to the resolution of the two following problems

$$(20) \qquad \begin{aligned} -B_{\epsilon}''(x_2) + \frac{k_{\epsilon}^2}{M^2} B_{\epsilon}(x_2) &= b_{\epsilon}(x_2), \\ B_{\epsilon}(0) &= B_{\epsilon}(l) = 0, \end{aligned}$$

$$(21) \qquad \begin{aligned} -A_{\epsilon}''(x_2) + \frac{k_{\epsilon}^2}{M^2} A_{\epsilon}(x_2) &= 2i \frac{k_{\epsilon}}{M} B_{\epsilon}(x_2) + a_{\epsilon}(x_2), \\ A_{\epsilon}(0) &= A_{\epsilon}(l) = 0. \end{aligned}$$

One can easily prove that problem (20) is well-posed in the space $H_0^1([0, l])$ and consequently compute the function $B_{\epsilon}$, which allows in turn to determine $A_{\epsilon}$ (using the same argument) and to finally find the function $\zeta^{\epsilon}$.

Now using their explicit formulas (see the appendix), one can show $a_{\epsilon}$ and $b_{\epsilon}$ respectively tend to some functions $a$ and $b$ as $\epsilon$ tends to zero, which are related to $\psi = \text{curl } u$ through the relation

$$\psi(x_1, x_2) = (a(x_2) + x_1 b(x_2)) e^{i \frac{k}{M} x_1}, \quad \forall (x_1, x_2) \in ]d_+, +\infty[ \times [0, l].$$

We are then able to compute the respective limits of $A_{\epsilon}$ and $B_{\epsilon}$, $A$ and $B$, by simply solving similar problems. This way, we define the function

$$\zeta(x_1, x_2) = (A(x_2) + x_1 B(x_2)) e^{i \frac{k}{M} x_1}, \quad \forall (x_1, x_2) \in \Omega,$$

which verifies

$$D^2 \zeta - \Delta \zeta = \psi, \quad \forall (x_1, x_2) \in ]d_+, +\infty[ \times [0, l].$$

Now, consider a cut-off function $\chi \in C^{\infty}(\mathbb{R})$ such that

$$\chi(x_1) = \begin{cases} 1 & \text{if } x_1 > d_+, \\ 0 & \text{if } x_1 < d_-. \end{cases}$$

It is easy to see that $\tilde{\varphi}_h^\epsilon = \varphi_h^\epsilon - \chi \zeta^\epsilon$ satisfies the following problem: find $\tilde{\varphi}_h^\epsilon \in H^1(\Omega)$ such that

$$D_{\epsilon}^2 \tilde{\varphi}_h^\epsilon - \Delta \tilde{\varphi}_h^\epsilon = \tilde{g}_h^\epsilon \text{ in } \Omega,$$
$$\tilde{\varphi}_h^\epsilon = 0 \text{ on } \partial \Omega.$$

where $\tilde{g}_h^\epsilon = g_h + \psi^\epsilon - (D_{\epsilon}^2(\chi \zeta^\epsilon) - \Delta(\chi \zeta^\epsilon))$.
By construction, $\psi^\epsilon - (D_\epsilon^2(\chi\zeta^\epsilon) - \Delta(\chi\zeta^\epsilon))$ has a compact support contained in $\Omega_b$ and therefore $\tilde{g}_h^\epsilon$ is compactly supported in $\Omega_b$. Using the DtN operators, this last problem may equivalently be rewritten as a problem set in $\Omega_b$: find $\tilde{\varphi}_h^\epsilon \in H^1(\Omega_b)$ such that

$$ (22) \qquad \begin{aligned} D_\epsilon^2 \tilde{\varphi}_h^\epsilon - \Delta \tilde{\varphi}_h^\epsilon &= \tilde{g}_h^\epsilon \quad \text{in } \Omega_b, \\ \tilde{\varphi}_h^\epsilon &= 0 \quad \text{on } \partial\Omega \cap \partial\Omega_b, \\ \frac{\partial \tilde{\varphi}_h^\epsilon}{\partial \boldsymbol{n}} &= -T_\pm^{D,\epsilon} \tilde{\varphi}_h^\epsilon \quad \text{on } \Sigma_\pm, \end{aligned} $$

where the operators $T_\pm^{D,\epsilon}$ are defined as (the superscript D refers here to the Dirichlet boundary condition in problem (22))

$$ (23) \qquad T_\pm^{D,\epsilon}: H^{1/2}(\Sigma_\pm) \to H^{-1/2}(\Sigma_\pm) \\ \phi \mapsto \mp i \sum_{n=1}^{+\infty} \beta_n^{\epsilon\pm} (\phi, S_n)_{L^2(\Sigma_\pm)} S_n(x_2), $$

the numbers $\beta_n^{\epsilon\pm}$ being defined in (11) and

$$ (24) \qquad S_n(x_2) = \sqrt{\frac{2}{l}} \sin\left(\frac{n\pi}{l}x_2\right), \quad \forall n \in \mathbb{N}^*. $$

One can prove that this problem is well-posed, due to hypothesis (18). It is now possible to pass to the limit as $\epsilon$ tends to zero in the same way as done for the acoustic potential $\varphi_a^\epsilon$ and finally show that

$$ \lim_{\epsilon \to 0} \tilde{\varphi}_h^\epsilon = \tilde{\varphi}_h \text{ in } H^1(\Omega_b), $$

where $\tilde{\varphi}_h$ is the solution of: find $\tilde{\varphi}_h \in H^1(\Omega_b)$ such that

$$ (24) \qquad \begin{aligned} D^2 \tilde{\varphi}_h - \Delta \tilde{\varphi}_h &= \tilde{g}_h \quad \text{in } \Omega_b, \\ \tilde{\varphi}_h &= 0 \quad \text{on } \partial\Omega \cap \partial\Omega_b, \\ \frac{\partial \tilde{\varphi}_h}{\partial \boldsymbol{n}} &= -T_\pm^D \tilde{\varphi}_h \quad \text{on } \Sigma_\pm, \end{aligned} $$

where $\tilde{g}_h = g_h + \psi - (D^2(\chi\zeta) - \Delta(\chi\zeta))$ and with the following definition of the DtN operators $T_\pm^D$

$$ (25) \qquad T_\pm^D: H^{1/2}(\Sigma_\pm) \to H^{-1/2}(\Sigma_\pm) \\ \phi \mapsto \mp i \sum_{n=0}^{+\infty} \beta_n^{\pm} (\phi, S_n)_{L^2(\Sigma_\pm)} S_n(x_2). $$

Moreover, the function $\zeta^\epsilon$ obviously tends to $\zeta$ as $\epsilon$ tends to zero, and one has

$$ \lim_{\epsilon \to 0} \varphi_h^\epsilon = \lim_{\epsilon \to 0} (\tilde{\varphi}_h^\epsilon + \chi\zeta^\epsilon) = \tilde{\varphi}_h + \chi\zeta. $$

By uniqueness of the limit, we conclude that

$$ \varphi_h = \tilde{\varphi}_h + \chi\zeta. $$

In conclusion, we derived two ways of characterizing the limit $\varphi_h$ of the potential $\varphi_h^\epsilon$. We have shown in subsection 3.3.3 that the function $\varphi_h$ is the unique solution of (17) which vanishes upstream of the source. Besides, we have just demonstrated that $\varphi_h$ can also be defined as the sum $\varphi_h = \tilde{\varphi}_h + \chi\zeta$, where $\tilde{\varphi}_h$ is the unique (outgoing) solution of problem (24).
# 4 Setting of the problem with perfectly matched layers

Our goal in this section is to develop a finite element method to compute an approximation of the outgoing solution $u$ of equations (1) and (2). To do so, we must address two main difficulties. First, this problem is set in an unbounded domain. Second, the operator in Galbrun's equation is not coercive, making the finite element method unstable.

As already seen during the study of the dissipative problem, the coerciveness can be restored by applying a regularization technique. On the other hand, we plan to use PMLs (see for instance [3] and references therein for a presentation of this methodology) in order to truncate the computational domain. A posteriori, the regularization will prove to be necessary, not only for the finite element method, but also for the PML method (see subsection 5.5).

In a previous paper [3], we proved the convergence of the solutions of PML formulations of the scalar problem (16). Two different models of PMLs were considered: a “classical” one, derived directly from Bérenger's original model, and a modified one, designed to avoid a possible growing of the solution in the downstream layer (due to the presence of the so-called inverse upstream modes). This last property will be useful for the vectorial problem at hand and the modified model (recalled in the next subsection) will be the only one studied.

## 4.1 The PML formulation

We introduce the bounded domain $\Omega^L = \Omega_b \cup \Omega_\pm^L$ (see figure 2), composed of domain $\Omega_b$ and surrounding layers $\Omega_\pm^L$ defined by

$$ \Omega_-^L = \{(x_1, x_2) \in \Omega^L, x_- - L < x_1 < x_-\} $$

and

$$ \Omega_+^L = \{(x_1, x_2) \in \Omega^L, x_+ < x_1 < x_+ + L\}. $$

The external boundaries of the layers $\Omega_\pm^L$ are respectively denoted by $\Sigma_\pm$.

<figure><img src="image_5.png" /><figcaption>Figure 2: The bounded domain $\Omega^L$.</figcaption></figure>

The so-called modified PML model consists of the transformation of the differential operator

$$ (25) \qquad \frac{\partial}{\partial x_1} \rightarrow \alpha(x_1) \frac{\partial}{\partial x_1} + i\lambda(x_1) $$

in the governing equations of the problem. The complex function $\alpha$ is assumed to be unity in $\Omega_b$ and, in order to simplify the subsequent study, constant and equal to the complex scalar $\alpha^*$, satisfying the following hypotheses

$$ (26) \qquad \mathrm{Re}(\alpha^*) > 0, \mathrm{Im}(\alpha^*) < 0, $$

in $\Omega \setminus \Omega_b$ (see [3] for a justification). More generally, the function $\alpha$ depends on the variable $x_1$ in the layers. For instance, in the original Bérenger's model, one has

$$ \alpha(x_1) = \frac{-i\omega}{-i\omega + \sigma(x_1)}, $$
where the real, positive function $\sigma$ vanishes in $\Omega_b$. The function $\lambda$ is assumed to be zero in $\Omega_b$ and to be constant and equal to

$$ (27) \qquad \lambda^* = -\frac{kM}{1 - M^2} $$

in $\Omega \setminus \Omega_b$. In Bérenger's model, one has $\lambda \equiv 0$.

As a consequence of the transformation (25), the various modified operators will now be indexed by $\alpha$ and $\lambda$. For instance, we have the following definitions

$$ \nabla_{\alpha, \lambda} \mathbf{v} = \begin{pmatrix} \alpha \frac{\partial v_1}{\partial x_1} + i\lambda v_1 & \alpha \frac{\partial v_2}{\partial x_1} + i\lambda v_2 \\ \frac{\partial v_1}{\partial x_2} & \frac{\partial v_2}{\partial x_2} \end{pmatrix}, \quad D_{\alpha, \lambda} \mathbf{v} = -i\omega \mathbf{v} + M \left( \alpha \frac{\partial \mathbf{v}}{\partial x_1} + i\lambda \mathbf{v} \right), $$

and so on.

As seen in section 3, Galbrun's equation must be regularized in order to be numerically solved in a stable fashion by a nodal (i.e., $H^1$-conforming) finite element method. The PML formulation of the problem is no exception to this rule. Consequently, we introduce a function $\psi_{\alpha,\lambda}$, defined as the unique solution of

$$ D_{\alpha, \lambda}^2 \psi_{\alpha, \lambda} = \mathrm{curl} \, \mathbf{f} \, \text{in} \, \Omega, $$

which vanishes upstream of the source. One can easily verify that $\psi_{\alpha,\lambda} = \psi$ in $\Omega_b$ and that, $\forall (x_1, x_2) \in \Omega_+$,

$$ \psi_{\alpha,\lambda}(x_1, x_2) = e^{i\left(\frac{k}{M}x_++\left(\frac{k}{M}-\lambda^*\right)\frac{x_1-x_+}{\alpha^*}\right)} \left(a(x_2) + \left(x_+ + \frac{x_1-x_+}{\alpha^*}\right)b(x_2)\right). $$

We then pose the following problem for the approximated displacement field: find $\mathbf{u}^L \in H^1(\Omega^L)^2$ such that

$$ (28) \qquad \begin{array}{l} D_{\alpha, \lambda}^2 \mathbf{u}^L - \nabla_{\alpha, \lambda} (\mathrm{div}_{\alpha, \lambda} \mathbf{u}^L) + \mathrm{curl}_{\alpha, \lambda} (\mathrm{curl}_{\alpha, \lambda} \mathbf{u}^L - \psi_{\alpha, \lambda}) = \mathbf{f} \quad \text{in } \Omega^L, \\ \mathbf{u}^L \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega^L, \\ \mathrm{curl}_{\alpha, \lambda} \mathbf{u}^L = \psi_{\alpha, \lambda} \quad \text{on } \partial\Omega^L. \end{array} $$

## 4.2 Well-posedness

We establish a variational formulation of problem (28): find $\mathbf{u}^L \in V(\Omega^L) = \{\mathbf{v} \in H^1(\Omega^L)^2 \mid \mathbf{v} \cdot \mathbf{n} = 0 \text{ on } \partial\Omega^L\}$ such that

$$ (29) \qquad a_{\Omega^L}(\mathbf{u}^L, \mathbf{v}) + b_{\Omega^L}(\mathbf{u}^L, \mathbf{v}) = l_{\Omega^L}(\mathbf{v}), \quad \forall \mathbf{v} \in V(\Omega^L), $$

where

$$ \begin{align*} a_{\Omega^L}(\mathbf{u}, \mathbf{v}) &= \int_{\Omega^L} \left( \mathbf{u} \cdot \overline{\mathbf{v}} + \frac{1}{\alpha} (\mathrm{div}_{\alpha, \lambda} \mathbf{u} \cdot \mathrm{div}_{\alpha, -\lambda} \overline{\mathbf{v}} + \mathrm{curl}_{\alpha, \lambda} \mathbf{u} \cdot \mathrm{curl}_{\alpha, -\lambda} \overline{\mathbf{v}}) - \alpha M^2 \frac{\partial \mathbf{u}}{\partial x_1} \cdot \frac{\partial \overline{\mathbf{v}}}{\partial x_1} \right) \mathrm{d}\mathbf{x}, \\ b_{\Omega^L}(\mathbf{u}, \mathbf{v}) &= \int_{\Omega^L} \left( \frac{1}{\alpha} (-k^2 + 2kM\lambda - M^2\lambda^2 - \alpha) \mathbf{u} \cdot \overline{\mathbf{v}} + (i\lambda M^2 - 2ikM) \frac{\partial \mathbf{u}}{\partial x_1} \cdot \overline{\mathbf{v}} - i\lambda M^2 \mathbf{u} \cdot \frac{\partial \overline{\mathbf{v}}}{\partial x_1} \right) \mathrm{d}\mathbf{x}, \\ l_{\Omega^L}(\mathbf{v}) &= \int_{\Omega^L} \frac{1}{\alpha} (\mathbf{f} \cdot \overline{\mathbf{v}} + \psi_{\alpha, \lambda} \mathrm{curl}_{\alpha, -\lambda} \overline{\mathbf{v}}) \mathrm{d}\mathbf{x} - \int_{\Sigma_{\pm}^L} M^2 \psi_{\alpha, \lambda} \overline{v_2}(\mathbf{n} \cdot \mathbf{e}_1) \mathrm{d}\sigma. \end{align*} $$

**Theorem 7** If the assumptions (26) is satisfied, the variational problem (29) is of Fredholm type.
*Proof.* We prove the sesquilinear form $a_{\Omega^L}(\cdot, \cdot)$ defines, via the Riesz representation theorem, an operator which is the sum of an isomorphism and a compact operator on $V(\Omega^L)$. We write $a_{\Omega^L}(\cdot, \cdot)$ as

$$a_{\Omega^L}(\mathbf{u}, \mathbf{v}) = a_{\Omega^L}^0(\mathbf{u}, \mathbf{v}) + \lambda a_{\Omega^L}^1(\mathbf{u}, \mathbf{v}) + \lambda^2 a_{\Omega^L}^2(\mathbf{u}, \mathbf{v}),$$

where the sesquilinear forms $a_{\Omega^L}^i(\cdot, \cdot)$, $i = 0, 1, 2$, which are independent of $\lambda$, are defined by

$$a_{\Omega^L}^0(\mathbf{u}, \mathbf{v}) = \int_{\Omega^L} \left( \mathbf{u} \cdot \bar{\mathbf{v}} + \frac{1}{\alpha} (\operatorname{div}_{\alpha,0} \mathbf{u} \operatorname{div}_{\alpha,0} \bar{\mathbf{v}} + \operatorname{curl}_{\alpha,0} \mathbf{u} \operatorname{curl}_{\alpha,0} \bar{\mathbf{v}}) - \alpha M^2 \frac{\partial \mathbf{u}}{\partial x_1} \cdot \frac{\partial \bar{\mathbf{v}}}{\partial x_1} \right) \mathrm{d}\mathbf{x},$$

$$a_{\Omega^L}^1(\mathbf{u}, \mathbf{v}) = \int_{\Omega^L} \frac{i}{\alpha} (u_1 \operatorname{div}_{\alpha,0} \bar{\mathbf{v}} - \bar{v}_1 \operatorname{div}_{\alpha,0} \mathbf{u} + u_2 \operatorname{curl}_{\alpha,0} \bar{\mathbf{v}} - \bar{v}_2 \operatorname{curl}_{\alpha,0} \mathbf{u}) \mathrm{d}\mathbf{x},$$

$$a_{\Omega^L}^2(\mathbf{u}, \mathbf{v}) = \int_{\Omega^L} \frac{1}{\alpha} \mathbf{u} \cdot \bar{\mathbf{v}} \mathrm{d}\mathbf{x}.$$

The sesquilinear form $a_{\Omega^L}^0(\cdot, \cdot)$ is coercive on $V(\Omega^L)$. Indeed, for all $\mathbf{u}, \mathbf{v}$ in $V(\Omega^L)$, we have

$$\int_{\Omega^L} (\operatorname{curl}_{\alpha,0} \mathbf{u} \operatorname{curl}_{\alpha,0} \bar{\mathbf{v}} + \operatorname{div}_{\alpha,0} \mathbf{u} \operatorname{div}_{\alpha,0} \bar{\mathbf{v}}) \mathrm{d}\mathbf{x} = \int_{\Omega^L} \nabla_{\alpha,0} \mathbf{u} : \nabla_{\alpha,0} \bar{\mathbf{v}} \mathrm{d}\mathbf{x}.$$

and due to hypotheses (26), we have for all $\mathbf{u}$ of $V(\Omega^L)$

$$\begin{aligned} \operatorname{Re} (a_{\Omega^L}^0(\mathbf{u}, \mathbf{u})) &= \int_{\Omega^L} \left( |\mathbf{u}|^2 + \operatorname{Re}(\alpha)(1 - M^2) \left| \frac{\partial \mathbf{u}}{\partial x_1} \right|^2 + \operatorname{Re} \left( \frac{1}{\alpha} \right) \left| \frac{\partial \mathbf{u}}{\partial x_2} \right|^2 \right) \mathrm{d}\mathbf{x} \\ &\geq C_{\alpha} \|\mathbf{u}\|_{H^1(\Omega^L)^2}, \end{aligned}$$

with $C_{\alpha} = \min(1 - M^2, \operatorname{Re}(\alpha^*)(1 - M^2), \operatorname{Re}(\frac{1}{\alpha^*}))$.

On the other hand, the forms $a_{\Omega^L}^1(\cdot, \cdot)$ and $a_{\Omega^L}^2(\cdot, \cdot)$ both define a compact operator on $V(\Omega^L)$, due to the compactness of the embedding of $H^1(\Omega^L)$ into $L^2(\Omega^L)$. The same argument is used to show that the bounded operator defined on $H^1(\Omega^L)^2$ by the sesquilinear form $b_{\Omega^L}(\cdot, \cdot)$ is compact. Finally, the forms $a_{\Omega^L}(\cdot, \cdot)$, $b_{\Omega^L}(\cdot, \cdot)$ and $l_{\Omega^L}(\cdot)$ are clearly continuous on their respective spaces, which ends the proof. □

We now prove the equivalence between the variational formulation (29) and the following (strong) problem: find $\mathbf{u}^L \in H^1(\Omega^L)^2$ such that

$$\begin{aligned} & \mathrm{D}^2_{\alpha, \lambda} \mathbf{u}^L - \nabla_{\alpha, \lambda} (\operatorname{div}_{\alpha, \lambda} \mathbf{u}^L) = \mathbf{f} \quad \text{in } \Omega^L, \\ & \operatorname{curl}_{\alpha, \lambda} \mathbf{u}^L = \psi_{\alpha, \lambda} \quad \text{in } \Omega^L, \\ & \mathbf{u}^L \cdot \mathbf{n} = 0 \quad \text{on } \partial \Omega^L. \end{aligned}$$

**Lemma 2** There exists a strictly positive constant $\delta$, depending on $k$, $M$ and $\theta = \arg(\alpha^*)$, such that if $L/|\alpha^*| \ge \delta$, any solution of variational problem (29) is such that $\operatorname{curl}_{\alpha, \lambda} \mathbf{u}^L = \psi_{\alpha, \lambda}$ in $\Omega^L$.

*Proof.* Considering a test function of the form $\mathbf{v} = \operatorname{curl}_{\alpha, -\lambda} \phi$ with $\phi \in H^3(\Omega^L) \cap H_0^1(\Omega^L)$ in (29), we obtain, after integration by parts,

$$\int_{\Omega^L} \frac{1}{\alpha} \operatorname{curl}_{\alpha, \lambda} \mathbf{u}^L (\mathrm{D}^2_{\alpha, -\lambda} \bar{\phi} - \Delta_{\alpha, -\lambda} \bar{\phi}) \mathrm{d}\mathbf{x} = \int_{\Omega^L} \frac{1}{\alpha} (\bar{\phi} \operatorname{curl}_{\alpha, \lambda} \mathbf{f} - \psi_{\alpha, \lambda} \Delta \bar{\phi}) \mathrm{d}\mathbf{x} \\ + \int_{\Sigma^L \cup \Sigma^L_+} \alpha M^2 \psi_{\alpha, \lambda} \frac{\partial \bar{\phi}}{\partial x_1} (\mathbf{n} \cdot \mathbf{e}_1) \mathrm{d}\sigma.$$

INRIA
Knowing that the function $\psi_{\alpha,\lambda}$ satisfies

$$D_{\alpha,\lambda}^2 \psi_{\alpha,\lambda} = \mathrm{curl}_{\alpha,\lambda} \boldsymbol{f} \text{ in } \Omega^L,$$

we finally have

$$\int_{\Omega^L} \frac{1}{\alpha} (\mathrm{curl}_{\alpha,\lambda} \boldsymbol{u}^L - \psi_{\alpha,\lambda}) (D_{\alpha, -\lambda}^2 \bar{\phi} - \Delta_{\alpha, -\lambda} \bar{\phi}) \, d\boldsymbol{x} = 0.$$

By a density argument, this equality is valid for any function $\phi$ in $H^2(\Omega^L) \cap H_0^1(\Omega^L)$. Then, if $L/|\alpha^*|$ is large enough, the operator $D_{\alpha, -\lambda}^2 - \Delta_{\alpha, -\lambda}$ is surjective from $H^2(\Omega^L) \cap H_0^1(\Omega^L)$ onto $L^2(\Omega^L)$ (see subsection 5.1) and we deduce that $\mathrm{curl}_{\alpha,\lambda} \boldsymbol{u}^L = \psi_{\alpha,\lambda}$ almost everywhere in $\Omega^L$. □

We finish this study by an existence and uniqueness result. In addition to hypothesis (18) on the wave number $k$, we now assume that

$$(30) \qquad k \neq \frac{n\pi}{l}, \quad \forall n \in \mathbb{N}.$$

The need for this second assumption will be made clear in the next section.

**Theorem 8** Assume that hypotheses (26) are satisfied and choice (27) is made. Then, there exists a strictly positive constant $\delta$ such that problem (29) admits a unique solution if $L/|\alpha^*| \ge \delta$.

*Proof.* The problem being of Fredholm type, proving the existence of its solution amounts to proving uniqueness. Let us consider two solutions of problem (29) and let $w$ denote their difference. By the previous lemma, there exists $\delta > 0$ such that if $L/|\alpha^*| \ge \delta$ then $w$ a solution to the problem: find $w \in H^1(\Omega^L)^2$ such that

$$D_{\alpha,\lambda}^2 w - \nabla_{\alpha,\lambda}(\mathrm{div}_{\alpha,\lambda} w) = 0 \text{ in } \Omega^L,$$

$$\mathrm{curl}_{\alpha,\lambda} w = 0 \text{ in } \Omega^L,$$

$$w \cdot n = 0 \text{ on } \partial\Omega^L.$$

We then consider the function $\tilde{w}(x_1, x_2) = w(x_1, x_2) e^{i\lambda x_1}$, so that $\mathrm{curl}_{\alpha,0}\tilde{w} = \mathrm{curl}_{\alpha,\lambda}w e^{i\lambda x_1}$, and use the following result.

**Theorem 9** A function $v$ in $L^2(\Omega^L)^2$ satisfies $\mathrm{curl}_{\alpha,0}v = 0$ in $\Omega^L$ if and only if there exists a function $\phi$ in $H^1(\Omega^L)$ such that $v = \nabla_{\alpha,0}\phi$. This function is unique up to an additive constant.

*Proof.* The proof is simply adapted from the proof of theorem 2.9 of [15]. □

The field $\tilde{w}$ then derives from a scalar potential $\phi$, which is a solution of

$$\nabla_{\alpha,0} (D_{\alpha,0}^2 \phi - \Delta_{\alpha,0} \phi) = 0 \text{ in } \Omega^L,$$

$$\nabla_{\alpha,0} \phi \cdot n = 0 \text{ on } \partial \Omega^L,$$

hence

$$(31) \qquad D_{\alpha,0}^2 \phi - \Delta_{\alpha,0} \phi = C \text{ in } \Omega^L,$$

$$\nabla_{\alpha,0} \phi \cdot n = 0 \text{ on } \partial \Omega^L,$$

where $C$ is a complex constant. This last problem is well-posed if $L/|\alpha^*|$ is large enough (to be proved later, see subsection 5.1). Besides, one can easily verify that $\phi = -C/k^2$ is a solution of system (31) and, thus, the unique solution of the problem. As a consequence, one has $w = 0$ in $\Omega^L$. □
# 5 Convergence results of PMLs for Galbrun's equation

Our objective is to prove that $u^L$, the solution to (29), tends to $u = \nabla\varphi_a + \text{curl}\varphi_h$ in $\Omega_b$ when the ratio $L/|\alpha^*|$ tends to infinity. A natural idea is to introduce two approximate potentials $\varphi_a^L$ and $\varphi_h^L$, which converge respectively to $\varphi_a$ and $\varphi_h$ and are such that $u^L = \nabla\varphi_a^L + \text{curl}\varphi_h^L$.

We now briefly present a sketch of the convergence proof and point out several difficulties that need to be addressed in the analysis as well. Indeed, it appears there is not a unique Helmholtz decomposition for $u^L$, which leaves us with the delicate task of choosing the adequate approximate potentials. Secondly, the "natural candidates" for these potentials satisfy scalar problems with boundary conditions that are, as we shall see, a priori unusual for PML problems.

Let us characterize the potentials $\varphi_h^L$ and $\varphi_a^L$ as the respective solutions of the following scalar problems: find $\varphi_h^L \in H^1(\Omega^L)$ such that

$$ (32a) \qquad D_{\alpha,\lambda}^2 \varphi_h^L - \Delta_{\alpha,\lambda} \varphi_h^L = g_h + \psi_{\alpha,\lambda} \quad \text{in } \Omega^L, $$

$$ (32b) \qquad \varphi_h^L = 0 \quad \text{on } \partial\Omega^L \cap \partial\Omega, $$

$$ (32c) \qquad -\Delta_{\alpha,\lambda} \varphi_h^L = \psi_{\alpha,\lambda} \quad \text{on } \Sigma_{\pm}^L, $$

and: find $\varphi_a^L \in H^1(\Omega^L)$ such that

$$ (33a) \qquad D_{\alpha,\lambda}^2 \varphi_a^L - \Delta_{\alpha,\lambda} \varphi_a^L = g_a \quad \text{in } \Omega^L, $$

$$ (33b) \qquad \nabla_{\alpha,\lambda} \varphi_a^L \cdot n = 0 \quad \text{on } \partial\Omega^L \cap \partial\Omega, $$

$$ (33c) \qquad \nabla_{\alpha,\lambda} \varphi_a^L \cdot n = -\text{curl}_{\alpha,\lambda} \varphi_h^L \cdot n \quad \text{on } \Sigma_{\pm}^L. $$

The first step of the proof consists of showing that these two problems are well-posed. The field $\nabla\varphi_a^L + \text{curl}\varphi_h^L$ is then clearly a solution of problem (28), which establishes that $u^L = \nabla\varphi_a^L + \text{curl}\varphi_h^L$. It remains to prove the convergence of the potentials $\varphi_h^L$ and $\varphi_a^L$ to their counterparts $\varphi_h$ and $\varphi_a$, using the convergence analysis previously done in [3]. However, problems (32) and (33) do not enter exactly the framework considered in this reference, for three main reasons. First, the boundary condition (32c) appears to be non-standard and possesses a non-homogeneous datum. Second, the right-hand side term in equation (32a) has a part of its support contained in the PMLs (we nevertheless observe both data are exponentially decreasing in the layers). Third, the boundary condition (33c) is not homogeneous. Notice this last condition links the potential $\varphi_a^L$ to the tangential trace of $\text{curl}\varphi_h^L$ on the boundaries $\Sigma_{\pm}^L$. Since $\varphi_a$ satisfies an analogous problem with homogeneous boundary conditions, we expect this trace to tend to zero in order to be in a position to prove the convergence of the method. As a consequence, we will first investigate the convergence of $\varphi_h^L$ and then that of $\varphi_a^L$.

Before doing so, we first recall and give convergence results for scalar problems of the type of (32) and (33), but with compactly supported data and homogeneous boundary conditions.

## 5.1 Convergence results of PMLs for scalar problems

Some straightforward extensions of the results of [3] will be needed in the sequel. We do not detail their proofs here, as they are only slightly modified from the ones in the above reference.

We study the following model problem. Suppose $g \in L^2(\Omega_b)$ is a source with compact support and consider the scalar field $\varphi$, which satisfies

$$ (34) \qquad \begin{aligned} D^2\varphi - \Delta\varphi &= g \quad \text{in } \Omega, \\ \varphi &= 0 \quad \text{on } \partial\Omega, \end{aligned} $$
and an additional radiation condition at infinity, which selects the outgoing solution. As already seen in subsection 3.3.5, this (nonlocal) condition may be expressed through the DtN operators $T_{\pm}^D$ on $\Sigma_{\pm}$.

In the following subsections, we give three PML formulations of this model problem, each one of them using a different boundary condition at the end of the layers. One should note that the results obtained here are still valid if a homogeneous Neumann boundary condition is applied on $\partial\Omega$ and, for the sake of brevity, we do not duplicate the statements. Indeed, such a change induces only a modification of the modal basis that appears in the demonstrations, the sine functions $(S_n)_{n \in \mathbb{N}^*}$, introduced in (23), being replaced by the cosine functions $(C_n)_{n \in \mathbb{N}}$ from (13).

### 5.1.1 Problem A

We consider a PML formulation of problem (34) with a homogeneous Dirichlet boundary condition on $\Sigma_{\pm}^L$

$$ (35a) \qquad D_{\alpha, \lambda}^2 \varphi^L - \Delta_{\alpha, \lambda} \varphi^L = g \quad \text{in } \Omega^L, $$

$$ (35b) \qquad \varphi^L = 0 \quad \text{on } \partial\Omega^L. $$

Equation (35a) has to be understood in the distributional sense, so that it implies the following transmission conditions at the interfaces between $\Omega_b$ and $\Omega_{\pm}^L$

$$ (36) \qquad [\varphi^L]_{\Sigma_{\pm}} = 0 \quad \text{and} \quad [\nabla_{\alpha, \lambda} \varphi^L]_{\Sigma_{\pm}} = 0. $$

In $\Omega_b$, the function $\varphi^L$, solution to (35), is meant to be an approximation of $\varphi$, solution to (34).

Adapting the proofs of [3], one can easily show the

**Lemma 3** Assume problem (35) has a solution. Then, this solution can be written as

$$ \varphi^L(x_1, x_2) = \sum_{n=1}^{+\infty} (\varphi^L(x_{\pm}, .), S_n)_{L^2(\Sigma_{\pm})} (A_n^{+}(\pm L) e^{i\gamma_n^{+}(x_1-x_{\pm})} + A_n^{-}(\pm L) e^{i\gamma_n^{-}(x_1-x_{\pm})}) S_n(x_2) $$

in the layers $\Omega_{\pm}^L$, where $\gamma_n^{\pm} = \frac{\beta_n^{\pm} - \lambda^*}{\alpha^*}$ and $A_n^{\pm}(L) = \mp \frac{e^{i\beta_n^{\pm}L/\alpha^*}}{e^{i\beta_n^{\pm}L/\alpha^*} - e^{i\beta_n^{-}L/\alpha^*}}$.

We easily check that the scalars $(A_n^{\pm}(\pm L))_{n \in \mathbb{N}}$ are always defined. Actually, the denominator would vanish if $\frac{(\beta_n^{+} - \beta_n^{-})L}{\alpha^*} \in 2\pi\mathbb{Z}$, but, due to assumptions (18) and (26), it is never zero and always has non-zero imaginary part. We are then able to write exact boundary conditions satisfied by $\varphi^L|_{\Omega_{\pm}^L} = \varphi_{\pm}^L$ on $\Sigma_{\pm}$

$$ \left(\frac{\partial \varphi_{\pm}^{L}}{\partial x_{1}}\right)_{|\Sigma_{\pm}} = i \sum_{n=1}^{+\infty} (\varphi_{\pm}^{L}(x_{\pm}, .), S_{n})_{L^{2}(\Sigma_{\pm})} (A_{n}^{+}(\pm L)\gamma_{n}^{+} + A_{n}^{-}(\pm L)\gamma_{n}^{-}) S_{n}, $$

which in turn yield, using the transmission conditions (36), exact boundary conditions satisfied by $\varphi^L|_{\Omega_b} = \varphi_b^L$ on $\Sigma_{\pm}$

$$ \left(\frac{\partial \varphi_b^L}{\partial x_1}\right)_{|\Sigma_{\pm}} = i \sum_{n=1}^{+\infty} (\varphi_b^L(x_{\pm}, .), S_n)_{L^2(\Sigma_{\pm})} \nu_n(\pm L) S_n, $$

where we have set $\nu_n(L) = A_n^+(L)\beta_n^+ + A_n^-(L)\beta_n^-$, which is equal to

$$ \nu_n(L) = \beta_n^+ + \frac{\beta_n^- - \beta_n^+}{1 - e^{i(\beta_n^- - \beta_n^+)L/\alpha^*}}. $$
Observe this formula does not depend on the value of $\lambda^*$.

Having in mind the comparison between $\varphi^L$ and $\varphi$ in $\Omega_b$, we reformulate (35) as an equivalent problem posed solely in this domain: find $\varphi^L \in H^1(\Omega_b)$ such that

$$ (37) \qquad \begin{aligned} D\varphi^L - \Delta\varphi^L &= g \quad \text{in } \Omega_b, \\ \varphi^L &= 0 \quad \text{on } \partial\Omega_b \cap \partial\Omega, \\ \frac{\partial\varphi^L}{\partial n} &= -T_{\pm}^L \varphi^L \quad \text{on } \Sigma_{\pm}, \end{aligned} $$

with

$$ T_{\pm}^{L}: H^{1/2}(\Sigma_{\pm}) \to H^{-1/2}(\Sigma_{\pm}) \\ \phi \mapsto \mp i \sum_{n=1}^{+\infty} \nu_n(\pm L) (\phi, S_n)_{L^2(\Sigma_{\pm})} S_n(x_2). $$

On the other hand, problem (34) has the following equivalent formulation: find $\varphi \in H^1(\Omega_b)$ such that

$$ \begin{aligned} D\varphi - \Delta\varphi &= g \quad \text{in } \Omega_b, \\ \varphi &= 0 \quad \text{on } \partial\Omega_b \cap \partial\Omega, \\ \frac{\partial\varphi}{\partial n} &= -T_{\pm}^D \varphi \quad \text{on } \Sigma_{\pm}. \end{aligned} $$

Since, the scalars $\nu_n(\pm L)$ tend to $\beta_n^\pm$ as $L/|\alpha^*|$ tends to $+\infty$ for any integer $n$, the operators $T_{\pm}^L$ converge in some sense to the operators $T_{\pm}^D$. As a consequence, one can prove that problem (37) is well-posed if the ratio $L/|\alpha^*|$ is large enough, as done in [3]. Moreover, we have the following convergence result (see [3] for a proof).

**Theorem 10** Suppose that assumptions (18) and (26) are verified and let $g \in L^2(\Omega^L)$ with compact support in $\Omega_b$. Then, there exists a strictly positive constant $\delta$, depending on $k$, $M$ and $\theta = \arg(\alpha^*)$, such that if $L/|\alpha^*| \ge \delta$, problem (35) is well-posed. Furthermore, the restriction to $\Omega_b$ of the solution $\varphi^L$ of problem (35) converges to the restriction of the solution $\varphi$ of problem (34) as $L/|\alpha^*|$ tends to $+\infty$. There also exists a constant $C$, depending on $M$ and $k$, such that

$$ \|\varphi - \varphi^L\|_{H^2(\Omega_b)} \le C e^{-\eta|\alpha^*|} \|\varphi\|_{H^2(\Omega_b)}, $$

where

$$ (38) \qquad \eta = \frac{2k}{1 - M^2} \min \left( -\sin(\theta) \sqrt{1 - \frac{N_0^2}{K_0^2}}, \cos(\theta) \sqrt{\frac{(N_0 + 1)^2}{K_0^2} - 1} \right), $$

where $K_0 = \frac{kl}{\pi\sqrt{1 - M^2}}$, $N_0$ is the floor of $K_0$ (i.e., the closest integer to $K_0$ which does not exceed it) and $\theta = \arg(\alpha^*)$.

More generally, if the datum $g$ is not compactly supported and/or if the Dirichlet boundary condition considered at the end of the layers is not homogeneous, it is still possible to prove the well-posedness of the approximated problem by simply using the Fredholm alternative and then deducing the uniqueness from the well-posedness of problem (37).
### 5.1.2 Problem B

We now choose to impose the homogeneous Neumann boundary condition

$$ (39) \qquad \nabla_{\alpha,\lambda} \varphi^L \cdot n = 0 \text{ on } \Sigma_{\pm}^{L}, $$

instead of the previous Dirichlet boundary condition. In this case, the claim of Theorem 10 remains true if we furthermore assume that $k \neq \frac{n\pi}{L}$, $\forall n \in \mathbb{N}$. Indeed, the sketch of the proof for the problem at hand is nearly identical to the preceding one, with instead the following values

$$ A_{n}^{\pm}(L) = \mp \frac{\beta_{n}^{\mp} e^{i\beta_{n}^{\mp} L/\alpha^{*}}}{\beta_{n}^{+} e^{i\beta_{n}^{+} L/\alpha^{*}} - \beta_{n}^{-} e^{i\beta_{n}^{-} L/\alpha^{*}}} $$

and

$$ \nu_n(L) = \beta_n^+ + \frac{\beta_n^+(\beta_n^- - \beta_n^+)}{\beta_n^+ - \beta_n^- e^{i(\beta_n^- - \beta_n^+)L/\alpha^*}}. $$

Again, one can verify that the scalars $(A_n^\pm(L))_{n \in \mathbb{N}}$ are always defined if the assumptions (18) and (26) are satisfied. Obviously, the scalar $\nu_n(L)$ tends to $\beta_n^+$ as $L/|\alpha^*|$ tends to $+\infty$ for any integer $n$. On the other hand, if there exists an integer $j$ such that $k = \frac{j\pi}{L}$, then the corresponding axial wave number $\beta_j^+$ vanishes and the scalar

$$ \nu_j(-L) = \beta_j^- + \frac{\beta_j^-(\beta_j^+ - \beta_j^-)}{\beta_j^- - \beta_j^+ e^{i(\beta_j^- - \beta_j^+)L/\alpha^*}} = 0 $$

cannot converge to $\beta_j^- \neq 0$ when $L/|\alpha^*|$ tends to $+\infty$.

Additionally, the well-posedness of problem (35a)-(39) in the more general case of a source $g$ with non-compact support and/or if the boundary condition (39) is not homogeneous can be proved by means of the Fredholm alternative.

### 5.1.3 Problem C

We finally consider the following homogeneous condition

$$ (40) \qquad \Delta_{\alpha,\lambda} \varphi^L = 0 \text{ on } \Sigma_{\pm}^{L}. $$

Theorem 10 is still valid. This time, in the proof, we have

$$ (41) \qquad A_n^\pm(L) = \mp \frac{(M\beta_n^\mp - k)^2 e^{i\beta_n^\mp L/\alpha^*}}{(M\beta_n^+ - k)^2 e^{i\beta_n^+ L/\alpha^*} - (M\beta_n^- - k)^2 e^{i\beta_n^- L/\alpha^*}} $$

and

$$ \nu_n(L) = \beta_n^+ + \frac{(M\beta_n^+ - k)^2(\beta_n^- - \beta_n^+)}{(M\beta_n^+ - k)^2 - (M\beta_n^- - k)^2 e^{i(\beta_n^- - \beta_n^+)L/\alpha^*}}. $$

Assuming hypotheses (18) and (26), one can check the scalars $(A_n^\pm(L))_{n \in \mathbb{N}}$ are always defined in this case as well.

Contrary to both previous problems, the well-posedness of problem (35a)-(40) cannot be extended to the case of an arbitrary source term $g$, since the boundary condition (40) does not allow to write a variational formulation of the problem.
## 5.2 Well-posedness and convergence analysis for $\varphi_h^L$

Once again, the idea is to deal with an equivalent problem, whose source term is compactly supported and boundary conditions are homogeneous, and the approach previously used in subsection 3.3.5 is followed closely. This new problem then permits the construction of the solution via a modal decomposition.

We first introduce the functions $\psi_{\alpha,\lambda}^{\infty}$, such that, $\forall (x_1, x_2) \in \Omega_+^L$,

$$ \psi_{\alpha,\lambda}^{\infty}(x_1, x_2) = \begin{cases} (a(x_2) + b(x_2)x_1) e^{\frac{i}{M}x_1}, & \forall (x_1, x_2) \in ]d_+, x_+[\times]0, l[, \\ \left(a(x_2) + b(x_2)\left(x_+ + \frac{x_1 - x_+}{\alpha^*}\right)\right) e^{\frac{i}{M}\left(x_++\left(\frac{k}{M}-\lambda^*\right)\frac{x_1-x_+}{\alpha^*}\right)}, & \end{cases} $$

and $\zeta_{\alpha,\lambda}$, such that, $\forall (x_1, x_2) \in \Omega_+^L$,

$$ \zeta_{\alpha,\lambda}(x_1, x_2) = \begin{cases} (A(x_2) + B(x_2)x_1) e^{\frac{i}{M}x_1}, & \forall (x_1, x_2) \in ]d_+, x_+[\times]0, l[, \\ \left(A(x_2) + B(x_2)\left(x_+ + \frac{x_1 - x_+}{\alpha^*}\right)\right) e^{\frac{i}{M}\left(x_++\left(\frac{k}{M}-\lambda^*\right)\frac{x_1-x_+}{\alpha^*}\right)}, & \end{cases} $$

the functions $a, b, A$ and $B$ having been characterized in the subsection 3.3.5. We then set

$$ \varphi_h^L = \tilde{\varphi}_h^L + \chi\zeta_{\alpha,\lambda}, $$

where the function $\tilde{\varphi}_h^L$ satisfies

$$ \begin{aligned} D_{\alpha,\lambda}^2 \tilde{\varphi}_h^L - \Delta_{\alpha,\lambda} \tilde{\varphi}_h^L &= \tilde{g}_h \text{ in } \Omega^L, \\ \tilde{\varphi}_h^L &= 0 \text{ on } \partial\Omega^L \cap \partial\Omega, \\ -\Delta_{\alpha,\lambda} \tilde{\varphi}_h^L &= 0 \text{ on } \Sigma_{\pm}^L. \end{aligned} $$

Indeed, the quantity $g_h + \psi_{\alpha,\lambda} - D_{\alpha,\lambda}^2(\chi\zeta_{\alpha,\lambda}) - \Delta_{\alpha,\lambda}(\chi\zeta_{\alpha,\lambda})$ coincides with $\tilde{g}_h$, since $\psi_{\alpha,\lambda}$ and $\zeta_{\alpha,\lambda}$ respectively coincide with $\psi$ and $\zeta$ in $\Omega_b$. This last problem is well-posed, due to the results of subsection 5.1 and we easily deduce the following result.

**Theorem 11** If the ratio $L/|\alpha^*|$ is large enough, problem (32) has a unique solution which is $\varphi_h^L = \tilde{\varphi}_h^L + \chi\zeta_{\alpha,\lambda}$. Moreover, the function $\varphi_h^L$ converges to $\tilde{\varphi}_h + \chi\zeta = \varphi_h$ in $H^2(\Omega_b)$ as $L/|\alpha^*|$ tends to $+\infty$, and one has the following estimate

$$ \|\varphi_h - \varphi_h^L\|_{H^2(\Omega_b)} \le C e^{-\frac{\eta}{2}|\alpha^*|} \|\varphi_h\|_{H^2(\Omega_b)}, $$

where the constant $C$ depends on $k$ and $M$ and $\eta$ is defined in (38).

*Proof.* We remark that, in the domain $\Omega_b$, one has $\varphi_h^L - \varphi_h = \tilde{\varphi}_h^L + \chi\zeta_{\alpha,\lambda} - \tilde{\varphi}_h - \chi\zeta = \tilde{\varphi}_h^L - \tilde{\varphi}_h$. The convergence result then directly follows from the subsection 5.1. □

**Corollary 3** Suppose assumptions (26) and (27) hold. The traces $(\mathrm{curl}_{\alpha,\lambda}\varphi_h^L \cdot \mathbf{n})|_{\Sigma_{\pm}^L}$ tend to zero in $H^{1/2}(\Sigma_{\pm}^L)$ as $L/|\alpha^*|$ tends to $+\infty$. More precisely, for $L/|\alpha^*|$ large enough, one has the estimate

$$ \|\mathrm{curl}_{\alpha,\lambda} \varphi_h^L \cdot \mathbf{n}\|_{H^{1/2}(\Sigma_{\pm}^L)} \le C \left( e^{-\frac{\eta}{2}|\alpha^*|} \|\varphi_h\|_{H^2(\Omega_b)} + e^{\left(\frac{k}{M}-\lambda^*\right)\sin(\theta)} \frac{L}{|\alpha^*|} \|g_h\|_{L^2(\Omega_b)} \right), $$

where the constant $C$ depends on $k$ and $M$, $\eta$ is defined in (38) and $\theta$ denotes the argument of $\alpha^*$. □
Proof. Using a modal decomposition for $\tilde{\varphi}_h^L$ in the layers $\Omega_{\pm}^L$, similar to the ones used in the subsection 5.1, one has, on $\Sigma_+^L$ for instance,

$$ \frac{\partial \tilde{\varphi}_h^L}{\partial x_2}(x_++L, x_2) = \sum_{n=1}^{+\infty} (\tilde{\varphi}_h^L(x_+, \cdot), S_n)_{L^2([0,l])} (A_n^+(L)e^{i\gamma_n^+ L} + A_n^-(L)e^{i\gamma_n^- L}) \frac{n\pi}{l} C_n(x_2), \forall x_2 \in [0, l], $$

the scalars $(A_n^{\pm}(L))_{n \in \mathbb{N}}$ being defined in (41).

Setting $\tau_n = A_n^+(L) e^{i\gamma_n^+ L} + A_n^-(L) e^{i\gamma_n^- L}$, we obtain

$$ \tau_n = \frac{M(\beta_n^- - \beta_n^+)(2k + M(\beta_n^+ + \beta_n^-)) e^{i\gamma_n^+ L}}{(M\beta_n^+ - k)^2 e^{i(\beta_n^+ - \beta_n^-)\frac{L}{\alpha^*}} - (M\beta_n^- - k)^2}, \forall n \in \mathbb{N}, $$

and, if $L/|\alpha^*|$ is large enough and after some majorizations,

$$ |\tau_n| \le C e^{-\frac{\eta}{2} \frac{L}{|\alpha^*|}}, \forall n \in \mathbb{N}, $$

the constant $\eta$ being defined in (38). Thus, we have

$$ \begin{aligned} \left\| \frac{\partial \tilde{\varphi}_h^L}{\partial x_2} \right\|_{H^{1/2}(\Sigma_+^L)}^2 &= \sum_{n=1}^{+\infty} \left(1 + \frac{n^2 \pi^2}{l^2}\right)^{1/2} \left| \left( \frac{\partial \tilde{\varphi}_h^L}{\partial x_2}(x_+ + L, \cdot), C_n \right)_{L^2([0,l])} \right|^2 \\ &\le \sum_{n=1}^{+\infty} \left(1 + \frac{n^2 \pi^2}{l^2}\right)^{3/2} |\tau_n|^2 \left| (\tilde{\varphi}_h^L(x_+, \cdot), S_n)_{L^2([0,l])} \right|^2 \\ &\le C e^{-\eta \frac{L}{|\alpha^*|}} \left\| \tilde{\varphi}_h^L \right\|_{H^{3/2}(\Sigma_+)}^2. \end{aligned} $$

We then deduce an estimate of this quantity using the convergence of $\tilde{\varphi}_h^L$ to $\tilde{\varphi}_h$ in $H^2(\Omega_b)$ and a trace theorem.

On the other hand, one has

$$ \frac{\partial \zeta_{\alpha, \lambda}}{\partial x_2}(x_+ + L, x_2) = \left( A'(x_2) + B'(x_2) \left( x_+ + \frac{L}{\alpha^*} \right) \right) e^{i \left( \frac{k}{M} x_+ + \left( \frac{k}{M} - \lambda^* \right) \frac{L}{\alpha^*} \right)}, \forall x_2 \in [0, l]. $$

Hence,

$$ \begin{aligned} \left\| \frac{\partial \zeta_{\alpha, \lambda}}{\partial x_2} \right\|_{H^{1/2}(\Sigma_+^L)}^2 &\le \sum_{n=0}^{+\infty} \left(1 + \frac{n^2 \pi^2}{l^2}\right)^{1/2} \left| (A', C_n)_{L^2([0,l])} + (B', C_n)_{L^2([0,l])} \left(x_+ + \frac{L}{\alpha^*}\right) \right|^2 e^{2\left(\frac{k}{M}-\lambda^*\right)\sin(\theta)\frac{L}{|\alpha^*|}} \\ &\le \left( \left\| A \right\|_{H^2([0,l])}^2 + \left\| B \right\|_{H^2([0,l])}^2 \left| x_+ + \frac{L}{\alpha^*} \right|^2 \right) e^{2\left(\frac{k}{M}-\lambda^*\right)\sin(\theta)\frac{L}{|\alpha^*|}} \\ &\le C \left(1 + \frac{L}{|\alpha^*|}\right)^2 e^{2\left(\frac{k}{M}-\lambda^*\right)\sin(\theta)\frac{L}{|\alpha^*|}} \left\| g_h \right\|_{L^2(\Omega)}^2. \end{aligned} $$

One can obviously obtain similar estimates on $\Sigma_-^L$ and finally obtain the announced result. □

We emphasize here that this corollary is only valid for the PML model corresponding to the transformation (25) with the fixed choice (27) for $\lambda^*$. Indeed, every propagative mode, and more particularly every inverse upstream mode, has to be damped in the layers in order to prove the claim. This appears to be quite different from the convergence results previously given in [3] and subsequently extended in the subsection 5.1, which are also true for the original Bérenger's model (i.e., for $\lambda \equiv 0$).
## 5.3 Well-posedness and convergence analysis for $\varphi_a^L$

Since the analysis done in section 5.1 was concerned with PML problems with homogeneous boundary conditions on $\Sigma_{\pm}^L$, the only difficulty in proving the convergence of problem (33) comes from the non-homogeneous boundary condition (33c). Consequently, the field $\varphi_a^L$ is split in the following manner

$$ \varphi_a^L = \varphi_a^{L,1} + \varphi_a^{L,2}, $$

with

$$ (42) \qquad \begin{aligned} D_{\alpha,\lambda}^2 \varphi_a^{L,1} - \Delta_{\alpha,\lambda} \varphi_a^{L,1} &= g_a \quad \text{in } \Omega^L, \\ \nabla_{\alpha,\lambda} \varphi_a^{L,1} \cdot \mathbf{n} &= 0 \quad \text{on } \partial\Omega^L, \end{aligned} $$

and

$$ (43) \qquad \begin{aligned} D_{\alpha,\lambda}^2 \varphi_a^{L,2} - \Delta_{\alpha,\lambda} \varphi_a^{L,2} &= 0 \quad \text{in } \Omega^L, \\ \nabla_{\alpha,\lambda} \varphi_a^{L,2} \cdot \mathbf{n} &= 0 \quad \text{on } \partial\Omega^L \cap \partial\Omega, \\ \nabla_{\alpha,\lambda} \varphi_a^{L,2} \cdot \mathbf{n} &= -\mathbf{curl}_{\alpha,\lambda} \varphi_h^L \cdot \mathbf{n} \quad \text{on } \Sigma_{\pm}^L. \end{aligned} $$

As seen in subsection 5.1, both of these problems are well-posed. Moreover, the results of subsection 5.1.2 readily give the convergence of the solution $\varphi_a^{L,1}$ of problem (42) to $\varphi_a$ in $H^2(\Omega_b)$ as the ratio $L/|\alpha^*|$ tends to infinity.

We next prove the following lemma.

**Lemma 4** If the ratio $L/|\alpha^*|$ is large enough, the solution $\varphi_a^{L,2}$ of problem (43) satisfies the following estimate

$$ \|\varphi_a^{L,2}\|_{H^2(\Omega_b)} \le C e^{-\frac{\eta}{2}|\alpha^*|} \|\mathbf{curl}_{\alpha,\lambda} \varphi_h^L \cdot \mathbf{n}\|_{H^{1/2}(\Sigma_{\pm}^L)}, $$

where the constant $C$ depends on $k$ and $M$ and the constant $\eta$ is defined in (38).

*Proof.* We first set $q^{\pm} = - (\mathbf{curl}_{\alpha,\lambda} \varphi_h^L \cdot \mathbf{n})|_{\Sigma_{\pm}^L}$. The main tool of the proof is again a modal decomposition. Since the functions $q^{\pm}$ respectively belong to $H^{1/2}(\Sigma_{\pm}^L)$, they can be written as

$$ q^{\pm}(x_2) = \sum_{n=0}^{+\infty} q_n^{\pm} C_n(x_2) \quad \text{on } \Sigma_{\pm}^{L}, $$

and we have

$$ \|q^{\pm}\|_{H^{1/2}(\Sigma_{\pm}^L)}^2 = \sum_{n=0}^{+\infty} \left(1 + \frac{n^2\pi^2}{l^2}\right)^{1/2} |q_n^{\pm}|^2. $$

We look for a solution of problem (43) of the same form

$$ (44) \qquad \varphi_a^{L,2}(x_1, x_2) = \sum_{n=0}^{+\infty} \phi_n(x_1) C_n(x_2) \quad \text{in } \Omega^L, $$

which yields the following ODE's

$$ D_{\alpha,\lambda}^2 \phi_n - \left(\alpha \frac{d}{dx_1} + i\lambda\right)^2 \phi_n + \frac{n^2\pi^2}{l^2} \phi_n = 0, $$
with the boundary conditions

$$ \pm \left( \alpha \frac{d}{dx_1} + i\lambda \right) \phi_n = q_n^{\pm} \quad \text{on } \Sigma_{\pm}^{L}, $$

and the transmission conditions between $\Omega_b$ and the PMLs

$$ [\phi_n]_{\Sigma_{\pm}} = 0 \quad \text{and} \quad \left[ \left( \alpha \frac{d}{dx_1} + i\lambda \right) \phi_n \right]_{\Sigma_{\pm}} = 0. $$

There are three different zones

$$ \phi_n(x_1) = \begin{cases} A_n^- e^{i\gamma_n^- x_1} + A_n^+ e^{i\gamma_n^+ x_1} & \text{if } x_1 < x_-, \\ B_n^- e^{i\beta_n^- (x_1-x_+)} + B_n^+ e^{i\beta_n^+ (x_1-x_-)} & \text{if } x_- < x_1 < x_+, \\ C_n^- e^{i\gamma_n^- x_1} + C_n^+ e^{i\gamma_n^+ x_1} & \text{if } x_1 > x_+. \end{cases} $$

Expressing the boundary and transmission conditions gives a 6 by 6 linear system to solve in order to obtain the coefficients ($A_n^{\pm}$, $B_n^{\pm}$, $C_n^{\pm}$). After some manipulations, we finally obtain

$$ B_n^- = -i \frac{q_n^-}{\beta_n^-} e^{-i\gamma_n^- L} \frac{1+z_n^+}{1-z_n^+ z_n^-}, \quad \text{and} \quad B_n^+ = i \frac{q_n^+}{\beta_n^+} e^{i\gamma_n^+ L} \frac{1+z_n^-}{1-z_n^+ z_n^-}, $$

where

$$ z_n^+ = e^{i\beta_n^+ (x_+-x_-)} e^{2i\gamma_n^+ L} \quad \text{and} \quad z_n^- = e^{-i\beta_n^- (x_+-x_-)} e^{-2i\gamma_n^- L}. $$

Note that, due to assumption (30), the scalars $B_n^{\pm}$ are well defined.

Using classical properties of propagative and evanescent modes, it is easy to show, for $L/|\alpha^*|$ sufficiently large, the following estimates

$$ |B_n^{\pm}| \le C \left| \frac{q_n^{\pm}}{\beta_n^{\pm}} \right| e^{\mp \operatorname{Im}(\gamma_n^{\pm}) L}, $$

which then yield

$$ |B_n^{\pm}| \le C \left| \frac{q_n^{\pm}}{\beta_n^{\pm}} \right| e^{-\frac{\eta}{2} \frac{L}{|\alpha^*|}}, $$

the constant $\eta$ being the one defined in (38).

Successively integrating (44) with respect to $x_2$ and $x_1$ and using the estimates above, we conclude that

$$ \|\varphi_a^{L,2}\|_{H^2(\Omega_b)} \le C e^{-\frac{\eta}{2} \frac{L}{|\alpha^*|}} \|q^{\pm}\|_{H^{1/2}(\Sigma_{\pm}^L)}, $$

which ends the proof. □

Using Corollary 3, we finally deduce the following

**Corollary 4** The function $\varphi_a^{L,2}$ converges to 0 as $L/|\alpha^*|$ tends to $+\infty$. More precisely, for $L/|\alpha^*|$ large enough, one has the estimate

$$ \|\varphi_a^{L,2}\|_{H^2(\Omega_b)} \le C \left( e^{-\eta \frac{L}{|\alpha^*|}} \|\varphi_h\|_{H^2(\Omega_b)} + e^{\left(\left(\frac{k}{M}-\lambda^*\right)\sin(\theta)-\frac{\eta}{2}\right)\frac{L}{|\alpha^*|}} \|g_h\|_{L^2(\Omega_b)} \right), $$

where the constant $C$ depends on $k$ and $M$, $\eta$ is defined in (38) and $\theta$ denotes the argument of $\alpha^*$. □
## 5.4 Conclusion

The gathering of the preceding results yields the following inequality

$$ \|\mathbf{u}^L - \mathbf{u}\|_{H^1(\Omega_b)^2} \le C e^{-\frac{\eta}{2} \frac{L}{|\alpha^*|}} \left( \|\nabla \varphi_a\|_{H^1(\Omega_b)^2} + \|\mathbf{curl} \varphi_h\|_{H^1(\Omega_b)^2} + \left(1 + \frac{L}{|\alpha^*|}\right) e^{\left(\frac{k}{M} - \lambda^*\right) \sin(\theta) \frac{L}{|\alpha^*|}} \|g_h\|_{L^2(\Omega)} \right) $$

and consequently allows us to state one final theorem, relative to the convergence of the solution of problem (28) when the ratio $L/|\alpha^*|$ tends to infinity.

**Theorem 12** Suppose assumptions (26) and (27) hold. Then, the field $\mathbf{u}^L$ tends to $\mathbf{u}$ in $H^1(\Omega_b)^2$ as $L/|\alpha^*|$ tends to $+\infty$. Furthermore, for $L/|\alpha^*|$ large enough, one has the estimate

$$ \|\mathbf{u}^L - \mathbf{u}\|_{H^1(\Omega_b)^2} \le C e^{-\frac{\eta}{2} \frac{L}{|\alpha^*|}}, $$

where the constant $C$ depends on $k$, $M$ and the solution $\mathbf{u}$, the constant $\eta$ being defined in (38).

## 5.5 Remark on the use of PMLs without regularization

Let us finally tackle the claim made in the introduction of section 4, namely that PMLs do not work without regularization. Indeed, assume that the approximated displacement field $\mathbf{u}^L$ is computed as the solution of the following problem

$$ \begin{array}{l} D_{\alpha, \lambda}^2 \mathbf{u}^L - \nabla_{\alpha, \lambda} (\mathrm{div}_{\alpha, \lambda} \mathbf{u}^L) = \mathbf{f} \quad \text{in } \Omega^L, \\ \mathbf{u}^L \cdot \mathbf{n} = 0 \quad \text{on } \partial \Omega^L, \end{array} \qquad (45) $$

completed by an additional boundary condition at the end of the layers, which can be chosen as

$$ \mathrm{curl}_{\alpha, \lambda} \mathbf{u}^L = 0 \quad \text{on } \partial \Omega^L. $$

Other choices will produce the same type of results.

Then, the function

$$ \psi^L = \mathrm{curl}_{\alpha, \lambda} \mathbf{u}^L $$

is a solution of the following problem

$$ \begin{array}{l} D_{\alpha, \lambda}^2 \psi^L = \mathrm{curl} \, \mathbf{f} \quad \text{in } \Omega^L, \\ \psi^L = 0 \quad \text{on } \Sigma_{\pm}^L. \end{array} \qquad (46) $$

If the PML model works, the field $\mathbf{u}^L$ must converge to $\mathbf{u}$ in $\Omega_b$ as the ratio $L/|\alpha^*|$ tends to $+\infty$, and, consequently, $\psi^L$ must converge to $\psi = \mathrm{curl} \, \mathbf{u}$ in $\Omega_b$. We now show that this convergence does not hold. Indeed, the solution $\psi^L$ of problem (46) can be sought of the form

$$ \psi^L = \psi_{\alpha, \lambda} + \tilde{\psi}^L, $$

where $\tilde{\psi}^L$ solves the following problem

$$ \begin{array}{l} D_{\alpha, \lambda}^2 \tilde{\psi}^L = 0 \quad \text{in } \Omega^L, \\ \tilde{\psi}^L = 0 \quad \text{on } \Sigma_{-}^L, \\ \tilde{\psi}^L = -\psi_{\alpha, \lambda} \quad \text{on } \Sigma_{+}^L. \end{array} \qquad (47) $$

INRIA
Using the expression of $\psi_{\alpha,\lambda}$ derived above, we find that

$$ \psi_{\alpha,\lambda}(x_+ + L, x_2) = \left( a(x_2) + b(x_2) \left( x_+ + \frac{L}{\alpha^*} \right) \right) e^{i \left( \frac{k}{M} x_+ + \left( \frac{k}{M} - \lambda^* \right) \frac{L}{\alpha^*} \right)} , \quad \forall x_2 \in ]0, l[. $$

Then, the expression of $\tilde{\psi}^L$ can be easily derived, by seeking a solution of the form

$$ \tilde{\psi}^L(x_1, x_2) = \begin{cases} \tilde{a}(x_2) \left( 1 + t \frac{x_1 - x_-}{\alpha^*} \right) e^{i \left( \frac{k}{M} x_- + \left( \frac{k}{M} - \lambda^* \right) \frac{x_1 - x_-}{\alpha^*} \right)} & \text{if } x_1 < x^-, \\ \tilde{a}(x_2) (1 + t(x_1 - x_-)) e^{i \frac{k}{M} x_1} & \text{if } x^- \le x_1 \le x^+, \\ \tilde{a}(x_2) \left( 1 + t \left( x_+ - x_- + \frac{x_1 - x_+}{\alpha^*} \right) \right) e^{i \left( \frac{k}{M} x_+ + \left( \frac{k}{M} - \lambda^* \right) \frac{x_1 - x_+}{\alpha^*} \right)} & \text{if } x_1 > x^+, \end{cases} $$

where the function $\tilde{a}$ and the real number $t$ need to be determined. By construction, the transmission conditions on $\Sigma_{\pm}$ are automatically satisfied and we simply use the boundary conditions imposed on $\Sigma_{\pm}^L$, which give us the two following identities

$$ \begin{gathered} 1 - t \frac{L}{\alpha^*} = 0 \\ \tilde{a}(x_2) \left( 1 + t \left( x_+ - x_- + \frac{L}{\alpha^*} \right) \right) = - \left( a(x_2) + b(x_2) \left( x_+ + \frac{L}{\alpha^*} \right) \right) \end{gathered} $$

so that, finally, $t = \alpha^*/L$ and

$$ \tilde{a}(x_2) = - \frac{a(x_2) + b(x_2) \left( x_+ + \frac{L}{\alpha^*} \right)}{2 + \frac{\alpha^*}{L} (x_+ - x_-)} . $$

Clearly, we see that $\tilde{\psi}^L$ does not converge to zero as $L/|\alpha^*|$ tends to $+\infty$.

# 6 Numerical applications

We are interested in simulating the radiation of a compactly supported source situated in a two-dimensional rigid duct, problem for which no explicit reference solution is available. Nonetheless, we consider as a preliminary study the propagation of acoustic and vortical modes in order to validate the method.

## 6.1 Mode propagation in a rigid duct

### 6.1.1 Acoustic and vortical modes: some definitions

The so-called modes are solutions with separated variables of the homogeneous, non-regularized Galbrun's equation

$$ D^2 u - \nabla (\operatorname{div} u) = 0 \text{ in } \Omega, $$

with the rigid wall boundary condition

$$ u \cdot n = 0 \text{ on } \partial \Omega. $$

These solutions are of two distinct kinds, called, respectively, acoustic and vortical modes. The acoustic modes are of the form

$$ u(x_1, x_2) = \begin{cases} C e^{i \beta_n^{\pm} x_1} e_1 & \text{if } n = 0, \\ C e^{i \beta_n^{\pm} x_1} \left( -\frac{i \beta_n^{\pm} l}{n\pi} \cos\left(\frac{n\pi}{l}x_2\right) e_1 + \sin\left(\frac{n\pi}{l}x_2\right) e_2 \right) & \text{if } n \in \mathbb{N}^*, \end{cases} $$
where C is a complex constant and the axial wave numbers $\beta_n^{\pm}$ are given by (15). One can see that these fields are irrotational, hence their name. The acoustic modes associated to real valued axial wave numbers are called propagative and they are called evanescent otherwise. Propagative modes with positive (resp. negative) group velocity $\frac{\partial\omega}{\partial\beta}$ are called downstream (resp. upstream) modes, since their energy propagates downstream (resp. upstream) of the mean flow. Finally, we have seen that there may exist modes with a positive group velocity and a negative phase velocity $\frac{\omega}{\beta}$, which are called inverse upstream modes. These modes are known to cause instabilities in PMLs in time domain applications (see, for instance, [24, 4]).

The second “family” of solutions of Galbrun’s equation consists of a continuum of fields such that

$$u(x_1, x_2) = C e^{i \frac{k}{M} x_1} \left( \frac{ik}{M} \varphi'(x_2) e_1 + \varphi(x_2) e_2 \right),$$

where C is a complex constant and $\varphi$ denotes a scalar function belonging to $H_0^1([0, l])$. These solutions only propagate downstream and are called vortical modes, since they are divergence-free.

## 6.1.2 Description of the simulations

The following numerical simulations consist of solving problems similar to (28), with only one absorbing layer downstream (denoted by $\Omega_+^L$ on figure 2). Each problem was designed in such a way that one of the previously introduced modes is its exact solution. For the propagation of an acoustic mode, we have $f \equiv 0$ and $\psi_{\alpha,\lambda} \equiv 0$, the mode being imposed via a non-homogeneous condition on the boundary $\Sigma_-$ for the normal displacement $u \cdot n$. In the case of a vortical mode, we still have $f \equiv 0$ and a non-homogeneous boundary condition on $\Sigma_-$, but the field $\psi_{\alpha,\lambda}$ has to be computed a priori as the curl of the considered mode.

All computations were done with the finite element library MÉLINA [22]. We used $P_2$ Lagrange finite elements on a non-structured mesh and the length of the PML was equal to 10% of the length of the domain $\Omega_b$. As in the theoretical work already presented, the function $\alpha$ is constant in the layer and the argument of the complex number $\alpha^*$ is fixed and equal<fnref n="1" /> to $-\frac{\pi}{4}$, its modulus $|\alpha^*|$ being a parameter in the simulations.

## 6.1.3 Numerical results for acoustic modes

In the chosen configuration, characterized by the values $l = 1$, $k = 8$ and $M = 0.4$, six (i.e., three upstream and three downstream) acoustic modes are propagative. The curves plotting the relative error in the $H^1(\Omega_b)$ norm for the computed displacement versus the modulus of $\alpha^*$ for the propagative downstream modes are shown in Figure 3. We observe that each curve contains a minimum plateau where the relative error is below a few percent. For large values of $|\alpha^*|$, the error increases due to the reflection at the end of the layer and behaves as theoretically predicted. For small values of $|\alpha^*|$, the method diverges, the mesh resolution being too coarse to adequately represent the modes in the PML medium, thus producing spurious numerical errors. Similar results were obtained for the propagative upstream modes.

We show in Figures 4 to 6 the contours of the components of the computed displacement for a value of $|\alpha^*|$ such that the error of the method is below one percent.

## 6.1.4 Numerical results for vortical modes

For the study of the method on vortical perturbations, the transverse dependence of the modes is arbitrarily chosen as $\varphi(x_2) = \sin(\frac{m\pi}{l} x_2)$, where m is a given nonzero integer. In Figure 7 we show the

<footnote n="1">¹While apparently arbitrary, this choice simply makes the quantities $-\sin(\arg(\alpha^*))$ and $\cos(\arg(\alpha^*))$, which appear in the definition (38) of the coefficient $\eta$, equal.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3: Relative error in $H^1(\Omega_b)^2$ norm as a function of $|\alpha^*|$ for the computation of propagative downstream acoustic modes, $k = 8$ and $M = 0.4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Contours of the real part of the component $u_1$ of the computed displacement field for the propagative downstream mode $n = 0, k = 8, M = 0.4, \alpha^* = 0.25(1 - i)$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Contours of the real part of the components of the computed displacement field for the propagative downstream mode $n = 1, k = 8, M = 0.4, \alpha^* = 0.25(1 - i)$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Contours of the real part of the components of the computed displacement field for the propagative downstream mode $n = 2, k = 8, M = 0.4, \alpha^* = 0.25(1 - i)$.</figcaption></figure>

relative error of the method plotted versus the modulus of the coefficient $\alpha^*$, for values of the integer $m$ equal to 1, 2 and 3. The contours of the components of the corresponding computed solutions when the relative error is below one percent are presented in Figures 8 to 10.

The convergence of the PML method is also obtained in this case for appropriate values of $|\alpha^*|$. However, by comparing Figures 3 and 7, one can already point out a potential difficulty that may be encountered in practice when both irrotational and vortical perturbations are present, as the values of
<figure><img src="image_10.png" /><figcaption>Figure 7: Relative error in $H^1(\Omega_b)^2$ norm as a function of $|\alpha^*|$ for the vortical modes, $k = 8$ et $M = 0.4$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Contours of the real part of the components of the computed displacement field for a vortical mode $m = 1, k = 8, M = 0.4, \alpha^* = 0.65(1 - i)$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Contours of the real part of the components of the computed displacement field for a vortical mode $m = 2, k = 8, M = 0.4, \alpha^* = 0.65(1 - i)$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Contours of the real part of the components of the computed displacement field for a vortical mode $m = 3, k = 8, M = 0.4, \alpha^* = 0.65(1 - i)$.</figcaption></figure>

$|\alpha^*|$ that allow a good agreement between the exact and the computed solutions in the two cases are quite different. Indeed, the propagation constant of the vortical modes, which is equal to $\frac{k}{M}$, leads to more significant damping in the layer for these modes than for the acoustic ones. On the other hand, the finite element error for vortical modes is higher, since their wavelength is generally much shorter
than that of their acoustic counterparts. This fact may cause some discretization issues as $|\alpha^*|$ becomes small<fnref n="2" />.

However, a compromise can be found by using a thicker layer. For instance, we have repeated the previous simulations with a layer of length equal to 25% of the length of the domain $\Omega_b$. The obtained results are presented in Figure 11 and one can now observe a partial match between the respective ranges of values of $|\alpha^*|$ for which the relative errors are below a few percent.

<figure><img src="image_14.png" /><figcaption>Figure 11: Relative error in $H^1(\Omega_b)^2$ norm as a function of $|\alpha^*|$ for several modes, $k = 8$ et $M = 0.4$.</figcaption></figure>

## 6.2 Radiation of compactly supported sources

### 6.2.1 Acoustic source

We next simulated the radiation of an irrotational, compactly supported source $f$, placed in the duct. This case happens to be more complex than the previous one because of the absence of any reference solutions, which would permit the measuring of the precision of the method and the choosing of an adequate value for the parameter $|\alpha^*|$. The real part of the components of the computed displacement field are shown in Figure 12. The acoustic mode mainly radiated by the source is the one of index $n = 2$. The convective effect of the uniform flow can clearly be seen, as the wavelength of the computed solution is shorter upstream of the source than downstream.

<figure><img src="image_15.png" /><figcaption>Figure 12: Contours of the real part of the components of the computed displacement field for the radiation of an acoustic source ($k = 8, M = 0.4, \alpha^* = 0.5(1-i)$).</figcaption></figure>

<footnote n="2">²Note that these issues, related to the mesh size, appear as well without PMLs, when the speed of the fbw is slow (i.e., when the Mach number $M$ is close to 0) or at high frequencies (i.e., when the wave number $k$ is large).</footnote>
### 6.2.2 Rotational source

The simulation of the radiation of a compactly supported source which curl is nonzero is presented in Figure 13. It is located slightly upstream of the center of the domain $\Omega_b$. One can observe the hydrodynamic wake generated by the source and convected by the flow, whose amplitude increases linearly with respect to the coordinate $x_1$. The acoustic perturbations, whereas present, have a negligible amplitude compared to the vortical ones.

<figure><img src="image_16.png" /><figcaption>Figure 13: Contours of the real part of the components of the computed displacement field for the radiation of a vortical source ($k = 8, M = 0.4, \alpha^* = 0.5(1 - i)$).</figcaption></figure>

### 6.3 A concluding remark on the non-regularized PML problem

We did some numerical experiments and confirmed that PMLs produce very bad results if used with a non-regularized formulation of Galbrun's equation. Although we previously stressed that the regularization step was necessary in order to obtain a reliable numerical solution, free of spurious modes when solving the problem with the usual nodal finite elements, it has been observed that approximations using quadrilateral-based Lagrange elements on structured meshes are not subject to spectral pollution and appear to be stable [21]. Therefore, we can assert that the large errors observed in the simulations involving the non-regularized equation are due to the PML method and not to the finite element method.

## A Study of a second order transport equation

This appendix is devoted to the second order transport equation (5), which is written here as

$$ (48) \qquad D_{\varepsilon}^2 \psi^{\varepsilon} = g, $$

with $g$ a compactly supported source belonging to $L^2(\Omega)$. We know from Lemma 1 that equation (48) has a unique solution in $L^2(\Omega)$. We also have the following result.

**Lemma 5** *The solution in $L^2(\Omega)$ of equation (48) vanishes upstream of the support of the source $g$ and satisfies the estimate*

$$ (49) \qquad \|\psi^j\|_{L^2(\Omega)} \le C_{\varepsilon} \|g\|_{L^2(\Omega)}, $$

where $C_{\varepsilon}$ denotes a positive constant depending on the parameter $\varepsilon$.
*Proof.* From the proof of Lemma 1, we know that

$$ \psi^{\epsilon}(x_1, x_2) = G_{\epsilon} * g(\cdot, x_2)(x_1) = \frac{1}{M^2} \int_{-\infty}^{x_1} (x_1 - z) e^{i \frac{k_{\epsilon}}{M} (x_1 - z)} g(z, x_2) dz, $$

where the kernel $G_{\epsilon}$ denotes the causal Green's function of the differential operator $D_{\epsilon}^2$. Setting

$$ d_{-} = \min_{x_2 \in [0, h]} \{x_1 \in \mathbb{R} \mid (x_1, x_2) \in \text{supp } g\} $$

and

$$ d_{+} = \max_{x_2 \in [0, h]} \{x_1 \in \mathbb{R} \mid (x_1, x_2) \in \text{supp } g\}, $$

we consider the following cases.

If $x_1 < d_-$, one has $]-\infty, x_1] \cap [d_-, d_+] = \emptyset$, thus $\psi^{\epsilon}(x_1, x_2) \equiv 0$ for all $x_2 \in [0, l]$. The solution then vanishes upstream the support of the source.

If $d_- \le x_1 \le d_+$, one has

$$ \psi^{\epsilon}(x_1, x_2) = \frac{1}{M^2} \int_{d_{-}}^{x_1} (x_1 - z) e^{i \frac{k_{\epsilon}}{M} (x_1 - z)} g(z, x_2) dz. $$

The Cauchy-Schwarz inequality then yields

$$ |\psi^{\epsilon}(x_1, x_2)|^2 \le \frac{1}{M^4} \left( \int_{d_{-}}^{x_1} (x_1 - z)^2 e^{-\frac{2\epsilon}{M}(x_1 - z)} dz \right) \left( \int_{d_{-}}^{x_1} |g(z, x_2)|^2 dz \right), $$

and one obtains

$$ \int_{d_{-}}^{d_{+}} \int_{0}^{l} |\psi^{\epsilon}(x_1, x_2)|^2 dx_1 dx_2 \le C_1 \|g\|_{L^2(\Omega)}^2, $$

with $C_{1\epsilon} = \frac{(d_+ - d_-)^4}{M^4} e^{-\frac{2\epsilon}{M}(d_+ - d_-)}$.

Finally, if $x_1 > d_+$, we have

$$ \begin{aligned} \psi^{\epsilon}(x_1, x_2) &= \frac{1}{M^2} \int_{d_{-}}^{d_{+}} (x_1 - z) e^{i \frac{k_{\epsilon}}{M} (x_1 - z)} g(z, x_2) dz \\ &= \left( -\frac{1}{M^2} \int_{d_{-}}^{d_{+}} z e^{-i \frac{k_{\epsilon}}{M} z} g(z, x_2) dz + \frac{x_1}{M^2} \int_{d_{-}}^{d_{+}} e^{-i \frac{k_{\epsilon}}{M} z} g(z, x_2) dz \right) e^{i \frac{k_{\epsilon}}{M} x_1}. \end{aligned} $$

Observing that the variables in $\psi^{\epsilon}$ can be separated, so that

$$ \psi^{\epsilon}(x_1, x_2) = (a_{\epsilon}(x_2) + x_1 b_{\epsilon}(x_2)) e^{i \frac{k_{\epsilon}}{M} x_1}, $$

one arrives at

$$ \int_{d_{+}}^{+\infty} \int_{0}^{l} |\psi^{\epsilon}(x_1, x_2)|^2 dx_1 dx_2 \le \|a_{\epsilon}\|_{L^2(0,l)} \int_{d_{+}}^{+\infty} e^{-\frac{2\epsilon}{M}x_1} dx_1 + \|b_{\epsilon}\|_{L^2(0,l)} \int_{d_{+}}^{+\infty} |x_1| e^{-\frac{2\epsilon}{M}x_1} dx_1. $$

We have

$$ |a_{\epsilon}(x_2)|^2 \le \frac{1}{M^4} \int_{d_{-}}^{d_{+}} z^2 e^{-\frac{2\epsilon}{M}z} dz \int_{d_{-}}^{d_{+}} |g(z, x_2)|^2 dz $$

and

$$ |b_{\epsilon}(x_2)|^2 \le \frac{1}{M^4} \int_{d_{-}}^{d_{+}} e^{-\frac{2\epsilon}{M}z} dz \int_{d_{-}}^{d_{+}} |g(z, x_2)|^2 dz, $$
hence

$$ \int_{d_{+}}^{+\infty} \int_{0}^{l} |\psi^{\epsilon}(x_1, x_2)|^2 dx_1 dx_2 \leq C_{2\epsilon} \|g\|_{L^2(\Omega)}^2, $$

with $C_{2\epsilon} = \frac{1}{M^4} \max \left( \int_{d_{-}}^{d_{+}} z^2 e^{-\frac{2\epsilon}{M}z} dz \int_{d_{+}}^{+\infty} e^{-\frac{2\epsilon}{M}x_1} dx_1, \int_{d_{-}}^{d_{+}} e^{-\frac{2\epsilon}{M}z} dz \int_{d_{+}}^{+\infty} |x_1| e^{-\frac{2\epsilon}{M}x_1} dx_1 \right)$.

We then finally set $C_{\epsilon} = \sqrt{\max(C_{1\epsilon}, C_{2\epsilon})}$. □

## References

[1] S. ABARBANEL, D. GOTTLIEB, AND J. S. HESTHAVEN, Well-posed perfectly matched layers for advective acoustics, J. Comput. Phys., 154 (1999), pp. 266–283.

[2] A. BAYLISS AND E. TURKEL, Radiation boundary conditions for wave-like equations, Comm. Pure Appl. Math., 33 (1980), pp. 707–725.

[3] E. BÉCACHE, A.-S. BONNET-BEN DHIA, AND G. LEGENDRE, Perfectly matched layers for the convected Helmholtz equation, SIAM J. Numer. Anal., 42 (2004), pp. 409–433.

[4] E. BÉCACHE, S. FAUQUEUX, AND P. JOLY, Stability of perfectly matched layers, group velocities and anisotropic waves, J. Comput. Phys., 188 (2003), pp. 399–433.

[5] J.-P. BÉRENGER, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys., 114 (1994), pp. 185–200.

[6] A.-S. BONNET-BEN DHIA, L. DAHI, E. LUNÉVILLE, AND V. PAGNEUX, Acoustic diffraction by a plate in a uniform flow, Math. Models Methods Appl. Sci., 12 (2002), pp. 625–647.

[7] A.-S. BONNET-BEN DHIA, C. HAZARD, AND S. LOHRENGEL, A singular field method for the solution of Maxwell's equations in polyhedral domains, SIAM J. Appl. Math., 59 (1999), pp. 2028–2044.

[8] A.-S. BONNET-BEN DHIA, G. LEGENDRE, AND E. LUNÉVILLE, Analyse mathématique de l'équation de Galbrun en écoulement uniforme, C. R. Acad. Sci. Paris Sér. IIb Méc., 329 (2001), pp. 601–606.

[9] T. COLONIUS, S. K. LELE, AND P. MOIN, Boundary conditions for direct computation of aerodynamic sound generation, AIAA J., 31 (1993), pp. 1574–1582.

[10] M. COSTABEL, A coercive bilinear form for Maxwell's equations, J. Math. Anal. Appl., 157 (1991), pp. 527–541.

[11] D. M. EIDUS, The principle of limiting absorption, Amer. Math. Soc Transl., 47 (1965), pp. 157–191.

[12] B. ENGQUIST AND A. MAJDA, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp., 31 (1977), pp. 629–651.

[13] J. B. FREUND, Proposed inflow/outflow boundary condition for direct computation of aerodynamic sound, AIAA J., 35 (1997), pp. 740–742.

[14] H. GALBRUN, Propagation d'une onde sonore dans l'atmosphère terrestre et théorie des zones de silence, Gauthier-Villars, Paris, France, 1931.
[15] V. GIRAULT AND P.-A. RAVIART, *Finite element methods for Navier-Stokes equations, theory and algorithms*, vol. 5 of Springer series in computational mathematics, Springer-Verlag, Berlin, Germany, 1986.

[16] P. GRISVARD, *Elliptic problems in nonsmooth domains*, vol. 24 of Monographs and studies in mathematics, Pitman, London, Great Britain, 1985.

[17] —, * Singularities in boundary value problems*, vol. 22 of Research notes in applied mathematics, Masson, Paris, France, 1992.

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

[19] F. Q. HU, *On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer*, J. Comput. Phys., 129 (1996), pp. 201–219.

[20] —, *A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables*, J. Comput. Phys., 173 (2001), pp. 455–480.

[21] G. LEGENDRE, *Rayonnement acoustique dans un fluide en écoulement : analyse mathématique et numérique de l'équation de Galbrun*, PhD thesis, Université Paris VI, Paris, France, 2003.

[22] D. MARTIN, *On line documentation of MÉLINA*. http://perso.univ-rennes1.fr/daniel.martin/melina/www/homepage.html.

[23] B. POIRÉE, *Les équations de l'acoustique linéaire et non linéaire dans un écoulement de fluide parfait*, Acustica, 57 (1985), pp. 5–25.

[24] C. K. W. TAM, L. AURIAULT, AND F. CAMBULI, *Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains*, J. Comput. Phys., 144 (1998), pp. 213–234.

[25] K. W. THOMPSON, *Time dependent boundary conditions for hyperbolic systems*, I, J. Comput. Phys., 68 (1987), pp. 1–24.
<figure><img src="image_17.png" /></figure>

Unité de recherche INRIA Rocquencourt

Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex (France)

Unité de recherche INRIA Lorraine : LORIA, Technopôle de Nancy-Brabois - Campus scientifique 615, rue du Jardin Botanique - BP 101 - 54602 Villers-lès-Nancy Cedex (France)

Unité de recherche INRIA Rennes : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex (France)

Unité de recherche INRIA Rhône-Alpes : 655, avenue de l'Europe - 38330 Montbonnot-St-Martin (France)

Unité de recherche INRIA Sophia Antipolis : 2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France)

Éditeur

INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France)

http://www.inria.fr

ISSN 0249-6399