# Integral Equations and Iterative Schemes for Acoustic Scattering Problems

Xavier Antoine, Marion Darbas

► To cite this version:

Xavier Antoine, Marion Darbas. Integral Equations and Iterative Schemes for Acoustic Scattering Problems. F. Magoulès. Numerical Methods for Acoustics Problems, Saxe-Coburg Editors, 2016, 978-1-874672-46-3. hal-00591456

HAL Id: hal-00591456

https://hal.science/hal-00591456v1

Submitted on 9 May 2011

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Integral Equations and Iterative Schemes for Acoustic Scattering Problems

Xavier ANTOINE and Marion DARBAS

## Abstract

The aim of this chapter is to provide an introduction to the iterative Krylov solution of integral equations for time-harmonic acoustic scattering. From the point of view of computational methods, considering large frequencies in acoustics is challenging since it leads to solving highly indefinite large scale complex linear systems which generally implies a convergence breakdown of iterative methods. Most specifically, we develop the problematic and explain some partial solutions through analytical preconditioning for high frequency scattering and the introduction of new combined field integral equations.

**Keywords:** acoustic scattering, high-frequency, integral equation, GMRES, Krylov solver, preconditioners, well-conditioned integral equation

## 1 Introduction

The numerical solution of time-harmonic (exterior) scattering problems in the high frequency regime remains a challenging problem due to its specific computational
bottlenecks. A few possible directions can be considered to formulate the problem and try to solve it [15]. Here, we choose the point of view of integral equations formulations [33, 47, 55]. In the most recent industrial code developments, the integral equation formulation is approximated by a boundary element method. The resulting linear system is next solved by an iterative Krylov solver [58] coupled with a Fast Multilevel Multipole method [36, 37, 57, 62]. However, the Helmholtz operator for scattering problems is a highly indefinite complex-valued linear operator. As a consequence, the associated matrix resulting from the boundary element discretization is also highly indefinite and complex. This results in a breakdown of the iterative Krylov solver in many applications or to an extremely slow convergence in best cases. The aim of this chapter is to explain the specific difficulties linked to this problem and to introduce some partial solutions. In particular, the convergence of the solver is closely related to the spectral distribution of the integral operators. We propose here to modify the spectrum in such a way that we get a fast convergence. This is related to the idea of preconditioning but not only. Furthermore, the originality of our approach is that we work at the continuous level, meaning that we directly manipulate the Helmholtz operator and not its matrix representation. This point of view gives rise to an *analytical* operator theory for building preconditioners which must be compared to the standard purely *algebraic* approach. Finally, this chapter can be seen as a companion paper to the chapter [66] by E. Turkel on iterative scheme for acoustic scattering.

The plan of the chapter is the following. In Section 2, we give some basics about numerical linear algebra and Krylov solvers. In particular, we describe on a few explicit examples the problems of convergence that can arise when solving indefinite linear systems. This leads us to analyze and understand why high frequency exterior Helmholtz problems also suffer from a lack of convergence. A brief review of some notions and recent developments on *algebraic* preconditioners is also proposed for convergence improvement. As already said, one of the originalities of our approach for preconditioning integral operators is to build continuous *analytical* preconditioners. To this aim, Section 3 gives an introduction to the theory of pseudodifferential operators and associated symbolic calculus which will be used in the sequel. In particular, we give a first simple application of this theory for proposing an original interpretation of standard *algebraic* Laplace shifted preconditioners. Section 4 focuses on the theory of potential which is used for writing integral equation representations for acoustics scattering. We develop the elementary notions for obtaining integral equations and review standard direct and indirect, first- and second-kind Fredholm integral operator formulations. This allows us to insist on the spectral properties of each representation and to precise the well-posedness results. Section 5 gives a short presentation of the way actual integral equations solvers are built. This mainly highlights the fact that the construction of a preconditioner should be based on the only assumption that an integral operator is never available through its matrix but rather through a "black-box" which is able to compute fast matrix-vector products. Therefore, a suitable preconditioner should be matrix-free which is a strong restriction to usual algebraic preconditioners. We analyze in Section 6 the spectral distribution of standard integral operators in the case of a sphere where integral operators can be
diagonalized. This gives a thorough understanding of the eigenvalue distribution of each operator and how it is related to the frequency parameter $k$ and density $n_λ$ of discretization points per wavelength. In Section 7, we describe two possible solutions for obtaining some efficient and robust analytical preconditioners. The first strategy uses integral operators through the so-called Calderón relations while the second one is based on the theory of pseudodifferential operators. Some examples show that these two directions are promising. In Section 8, we propose to build some new alternative well-posed integral equations with a spectral distribution well-suited to get a fast convergence of an iterative solver. Examples show that these new integral formulations are well-adapted for solving high frequency problems. Finally, Section 9 gives a conclusion.

# 2 Difficulties for the iterative solution of scattering problems

## 2.1 Notations and background

Let us begin by some background in linear algebra. A vector norm on a vector space $\mathbb{Y}$ is a real-valued function $y \mapsto \|y\|$ on $\mathbb{Y}$ which satisfies:

* $\|y\| \ge 0, \forall y \in \mathbb{Y}$, and $\|y\| = 0$ if and only if $y = 0$ (positivity),

* $\|\alpha y\| = |\alpha|\|y\|, \forall y \in \mathbb{Y}, \forall \alpha \in \mathbb{C}$ (scaling),

* $\|y + x\| \le \|y\| + \|x\|, \forall y, x \in \mathbb{Y}$ (triangular inequality).

For the particular case when $\mathbb{Y} = \mathbb{C}^n$, the most commonly used vector norms are the Hölder norms ($p$-norms)

$$ \|y\|_p = (\sum_{i=1}^{n} |y_i|^p)^{1/p}, $$

with $y = (y_i)_{i=1,...,n}$ and $p \ge 1$. The most useful norms in practice are

$$ \|y\|_1 = \sum_{i=1}^{n} |y_i|, \quad \|y\|_2 = (\sum_{i=1}^{n} |y_i|^2)^{1/2}, \quad \text{and} \quad \|y\|_{\infty} = \max_{i=1,...,n} |y_i|. $$

We can use norms to measure the magnitude of a matrix. The same axioms stated above for vector norms apply here

* $\|A\| \ge 0, \forall A \in \mathbb{C}^{n \times n}$, and $\|A\| = 0$ if and only if $A = 0$ (positivity),

* $\|\alpha A\| = |\alpha|\|A\|, \forall A \in \mathbb{C}^{n \times n}, \forall \alpha \in \mathbb{C}$ (scaling),

* $\|A + B\| \le \|A\| + \|B\|, \forall A, B \in \mathbb{C}^{n \times n}$ (triangular inequality).
The most important class of norms are the induced matrix p-norms

$$ \|A\|_p = \max_{\mathbf{x} \in \mathbb{C}^n, \mathbf{x} \neq \mathbf{0}} \frac{\|A\mathbf{x}\|_p}{\|\mathbf{x}\|_p} $$

A fundamental property is that

$$ \|AB\|_p \le \|A\|_p \|B\|_p, \forall A, B \in \mathbb{C}^{n \times n}. $$

When $p = 1$ and $p = \infty$, we have the simple formulas $(A = (a_{ij})_{1 \le i,j \le n})$

$$ \|A\|_1 = \max_{j=1, \dots, n} \sum_{i=1}^{n} |a_{ij}| \quad \text{and} \quad \|A\|_{\infty} = \max_{i=1, \dots, n} \sum_{j=1}^{n} |a_{ij}|. $$

The induced 2-norm is also called the spectral norm and is given by

$$ \|A\|_2 = \rho(A^* A)^{1/2}, $$

with $A^*$ the transpose conjugate matrix (= $\bar{A}^T$) and $\rho(A)$ the maximum modulus of the eigenvalues of A (spectral radius).

To end this section, let us give some other useful definitions.

**Definition 1.** Let us introduce the following definitions

* A non Hermitian matrix A of size $n \times n$ is said to be positive-definite (respectively negative-definite) if and only if $\Re(\mathbf{x}^* A \mathbf{x}) > 0$ (respectively $\Re(\mathbf{x}^* A \mathbf{x}) < 0$) for all non-zero $\mathbf{x} \in \mathbb{C}^n$.

* A Hermitian matrix A ($A = A^*$) of size $n \times n$ is said to be positive-definite (respectively negative-definite) if and only if $\mathbf{x}^* A \mathbf{x} > 0$ (respectively $\mathbf{x}^* A \mathbf{x} < 0$) for all non-zero $\mathbf{x} \in \mathbb{C}^n$.

* A Hermitian matrix A of size $n \times n$ is said to be positive-semidefinite (respectively negative-semidefinite) if and only if $\mathbf{x}^* A \mathbf{x} \ge 0$ (respectively $\mathbf{x}^* A \mathbf{x} \le 0$) for all $\mathbf{x} \in \mathbb{C}^n$.

* A Hermitian matrix which is neither positive- or negative-semidefinite is called indefinite.

We also have the following Proposition.

**Proposition 1.** A Hermitian matrix is positive-definite if and only if all its eigenvalues are real and strictly positive.
## 2.2 Iterative algorithms

We consider a linear system

$$Ax = b, \tag{1}$$

where A is an invertible complex-valued matrix of size $n \times n$ and $b \in \mathbb{C}^n$. Methods for the solution of linear systems fall into two classes: direct methods and iterative methods.

Direct methods [31] produce an exact solution in a predictable finite number of elementary operations if no rounding errors are present. When the matrix is symmetric positive-definite, a Cholesky algorithm is applied. For a nonsymmetric matrix, a gaussian elimination solver is used. In this case, the memory storage and computational times costs scale as $O(n^2)$ and $O(n^3)$ respectively. Then, a direct solution is clearly out of reach when A is large and dense, requiring hence too large memory and prohibitive computational times.

Iterative schemes can be considered as an alternative to direct methods for the solution of large linear systems. Stationary relaxation-type methods [31] (Jacobi, Gauss-Seidel...) have the disadvantage of slow convergence and concern certain classes of matrices only. Projection methods are more general and robust. This is most particularly the case of Krylov subspace methods [58]. These techniques are motivated by the Cayley-Hamilton theorem [50] which allows to construct the inverse of a matrix as a polynomial according to A. Krylov subspace methods consist in seeking an approximate solution $x^{(m)}$ to (1) from an affine subspace

$$x^{(0)} + \mathcal{K}_m(A, \mathbf{r}^{(0)}), \tag{2}$$

with

* $x^{(0)}$ the initial guess,

* $\mathbf{r}^{(0)} = \mathbf{b} - A\mathbf{x}^{(0)}$ the initial residual vector,

* $\mathcal{K}_m(A, \mathbf{r}^{(0)}) = \text{span}\{\mathbf{r}^{(0)}, A\mathbf{r}^{(0)}, A^2\mathbf{r}^{(0)}, \dots, A^{m-1}\mathbf{r}^{(0)}\}$ the Krylov subspace of dimension $m$,

such that the orthogonality condition (Petrov-Galerkin condition)

$$(B(\mathbf{x} - \mathbf{x}^{(m)}), \mathbf{v}) = 0, \quad \forall \mathbf{v} \in \mathcal{K}_m(A, \mathbf{r}^{(0)}), \tag{3}$$

is fulfilled. The successive approximations are clearly expressed by $x^{(m)} = x^{(0)} + P_{m-1}(A)\mathbf{r}^{(0)}$, where $P_{m-1}$ is a polynomial of degree $m-1$. The different Krylov methods therefore differ in the choice of the matrix B (generally symmetric positive-definite). If B defines a scalar product, then the orthogonality condition (3) becomes equivalent to the following minimization problem (optimality condition)

$$\|\mathbf{x} - \mathbf{x}^{(m)}\|_B = \min_{\mathbf{y} \in \mathcal{K}_m(A, \mathbf{r}^{(0)})} \|\mathbf{x} - \mathbf{y}\|_B.$$

Two standard examples of such methods are
* Conjugate gradient methods: If A is symmetric positive-definite, then $B = A$ and (3) is equivalent to

$$ \|\mathbf{x} - \mathbf{x}^{(m)}\|_A = \min_{\mathbf{y} \in \mathcal{K}_m(A, \mathbf{r}^{(0)})} \|\mathbf{x} - \mathbf{y}\|_A, $$

* Generalized Minimum RESidual method (GMRES) : For a general matrix A, we choose $B = A^T A$ and get

$$ \|\mathbf{b} - A\mathbf{x}^{(m)}\|_2 = \min_{\mathbf{x} \in \mathcal{K}_m(A, \mathbf{r}^{(0)})} \|\mathbf{b} - A\mathbf{x}\|_2. $$

A Krylov subspace method is then represented by a construction algorithm of a basis of the affine subspace $\mathcal{K}_m(A, \mathbf{r}^{(0)})$ and by an optimality criterion to determine the approximate solution $\mathbf{x}^{(m)}$ of (1).

Let us now focus our attention on the GMRES which constitutes the reference algorithm in this chapter. This algorithm, introduced by Saad and Schultze [58, 59], is well-adapted to solve large nonsymmetric linear systems. It corresponds to the choice of the Krylov subspace $\mathcal{K}_m(A, \mathbf{v}_1)$ with $\mathbf{v}_1 = \mathbf{r}^{(0)}/\|\mathbf{r}^{(0)}\|_2$. At each iteration, the Krylov subspace $\mathcal{K}_m(A, \mathbf{v}_1)$ has to be constructed. In practice, we have to generate a set of basis of this subspace. The natural basis $(\mathbf{v}_1, A\mathbf{v}_1, A^2\mathbf{v}_1, \dots, A^{m-1}\mathbf{v}_1)$ cannot be used because of its numerical degeneracy. A solution is to construct an orthonormal basis $(\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m)$ of $\mathcal{K}_m(A, \mathbf{v}_1)$ via the Arnoldi-Modified Gramm-Schmidt algorithm. We denote by $V_m = (\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_m)$ the $n \times m$ matrix with columns vectors $\mathbf{v}_i, i = 1, \dots, n$, and by $\overline{H}_m$ the $(m+1) \times m$ Hessenberg matrix where nonzero entries $h_{ij} = (A\mathbf{v}_j, \mathbf{v}_i)$ are defined by the Arnoldi-Modified Gramm-Schmidt algorithm. Any vector $\mathbf{x}$ in $\mathbf{x}^{(0)} + \mathcal{K}_m(A, \mathbf{v}_1)$ can thus be written as $\mathbf{x} = \mathbf{x}^{(0)} + V_m\mathbf{y}$, where $\mathbf{y} \in \mathbb{C}^m$. Moreover, the relation $AV_m = V_{m+1}\overline{H}_m$ holds. It results in

$$ \mathbf{b} - A\mathbf{x} = \mathbf{b} - A(\mathbf{x}^{(0)} + V_m\mathbf{y}) = V_{m+1}(\beta\mathbf{e}_1 - \overline{H}_m\mathbf{y}), $$

by setting $\beta = \|\mathbf{r}^{(0)}\|_2$ and $\mathbf{e}_1 = (1, 0, \dots, 0)^T \in \mathbb{R}^{m+1}$. Then, we have the optimality criterion

$$ \min_{\mathbf{x} \in \mathcal{K}_m(A, \mathbf{r}^{(0)})} \|\mathbf{b} - A\mathbf{x}\|_2 = \min_{\mathbf{y} \in \mathbb{R}^m} \|\beta\mathbf{e}_1 - \overline{H}_m\mathbf{y}\|_2, $$

exploiting the fact that the column-vectors of $V_{m+1}$ are orthonormal. The approximate solution $\mathbf{x}^{(m)}$ can be obtained as $\mathbf{x}^{(m)} = \mathbf{x}^{(0)} + V_m\mathbf{y}^{(m)}$, where $\mathbf{y}^{(m)}$ minimizes the functional $J(\mathbf{y}) = \|\beta\mathbf{e}_1 - \overline{H}_m\mathbf{y}\|_2$. The pseudocode for basic form of the GMRES algorithm can now be given as

1. Initialization

$$ \text{Compute } \mathbf{r}^{(0)} = \mathbf{b} - A\mathbf{x}^{(0)}, \beta := \|\mathbf{r}^{(0)}\|_2 \text{ and } \mathbf{v}_1 := \mathbf{r}^{(0)} / \beta $$

2. Define the $(m+1) \times m$ matrix $\overline{H}_m = \{h_{ij}\}_{1 \le i \le m+1, 1 \le j \le m}$. Set $\overline{H}_m = 0$.

3. Construction of the Arnoldi's basis (steps 3-11)

for $j = 1, 2, \dots, m$, do
4. compute $\mathbf{w}_j := A\mathbf{v}_j$

5. for $i = 1, \dots, j$, do

6. $h_{ij} := (\mathbf{w}_j, \mathbf{v}_i)$

7. $\mathbf{w}_j := \mathbf{w}_j - h_{ij}\mathbf{v}_i$

8. enddo

9. $h_{j+1,j} = \|\mathbf{w}_j\|_2$. If $h_{j+1,j} = 0$ set $m := j$ and go to 12

10. $\mathbf{v}_{j+1} = \mathbf{w}_j/h_{j+1,j}$

11. enddo

12. Minimization problem

Compute $\mathbf{y}^{(m)}$ as the minimizer of $\|\beta\mathbf{e}_1 - \bar{\mathbf{H}}_m\mathbf{y}\|_2$ and $\mathbf{x}^{(m)} = \mathbf{x}^{(0)} + V_m\mathbf{y}^{(m)}$.

In order to solve the least-squares problem (step 12), the most adapted technique is to employ a QR-decomposition of the Hessenberg matrix $\bar{\mathbf{H}}_m$ (see [58] for details).

**Remark 1.** The parameter $m$ is not determined a priori. In general, a maximum number $m_{\max}$ is fixed and is typically dictated by the computational ressources. If the maximum number $m_{\max}$ of iterations has been reached without triggering the convergence test, then a restarting is done, i.e. GMRES is started afresh with the last approximation $\mathbf{x}^{(m_{\max})}$ as the initial guess. This method is called restarted GMRES($m_{\max}$). The residual $\|\mathbf{b} - A\mathbf{x}^{(m)}\|_2$ is generally used as a stopping criterion.

Essentially, the computational cost of GMRES is related to

i) the total number of iterations $N^{\text{iter}}$ required to reach an *a priori* fixed tolerance $\varepsilon$ on the residual norm,

ii) the cost of one iteration which is mainly the Matrix-Vector Product (MVP) $\mathbf{w}_j = A\mathbf{v}_j$ involved at step 4.

If this algorithm is directly used for a full complex-valued matrix $A$, then the total cost is $O(N^{\text{iter}}n^2)$. In terms of memory storage, the algorithm still needs $O(n^2)$ entries for $A$. In the background of integral equations, efficient compression algorithms have been proposed during the last two decades. For example, the multilevel Fast Multipole Method (FMM) [32, 36, 37, 40, 62] both computes the MVP in $O(n \log n)$ operations and requires $O(n)$ entries. Essentially, the FMM proposes to compute only the near-field entries exactly and the far-field entries approximately (but with controllable error). Other fast algorithms exist. We do not detail all these techniques and rather refer to [23] for technical points. This drastic reduction both in computational cost and memory storage for one iteration of the GMRES gives expectations for solving high frequency problems.
Let us now briefly review the most important results on the convergence behaviour of GMRES (for proof and further details see [39, 53, 58]). If $A \in \mathbb{C}^{n \times n}$ is non singular, we define the quantity

$$ \kappa_p(A) = \|A\|_p \|A^{-1}\|_p, $$

which is called the condition number of the linear system (1) with respect to the induced matrix p-norm (cf. section 2.1). We begin by giving a global convergence result.

**Theorem 1.** Let $A \in \mathbb{C}^{n \times n}$ be a non singular matrix. The full GMRES algorithm is guaranteed to converge in at most $n$ iterations.

This is true in the absence of arithmetic errors. Unfortunately, computer arithmetic is not exact. Moreover, this would be impractical for large linear systems if there were many steps required to reach convergence. In these situations, a preconditioner can be used to reduce the number of iterations (see e.g. Section 2.5). In order to predict the behaviour of GMRES, the convergence analysis is concerned with the derivation of upper bounds on the residual norms.

**Proposition 2.** Let $A \in \mathbb{C}^{n \times n}$ be a non singular matrix and $x^{(0)} \in \mathbb{C}^n$ an initial guess. If we can diagonalize $A$ in the form $A = U\Lambda U^{-1}$ where $\Lambda$ is the diagonal matrix of eigenvalues $(\lambda_i)_i$ of $A$ corresponding to the appropriate eigenvectors in $U$, then a bound on the residual norm at iteration $m$ is expressed by

$$ \|\mathbf{r}^{(m)}\|_2 \le \kappa_2(U) \min_{q \in \mathbb{P}_m, q(0)=1} \max_{i=1,...,n} q(\lambda_i) \|\mathbf{r}^{(0)}\|_2, \quad (4) $$

with $\kappa_2(U) = \|U\|_2 \|U^{-1}\|_2$ the condition number of $U$ in the 2-norm.

This bound was the first convergence result for GMRES [59]. However, even if $A$ is normal ($\kappa_2(U) = 1$), the bound (4) may fail to provide any reasonable information about the rate of reduction of the GMRES residual norms. For diagonalizable but non normal $A$ (far from unitary), $\kappa_2(U)$ might be very large, and the bound in (4) might be a large overestimate of the residual norm. Moreover, it is not clear that only the conditioning of the eigenvectors of $A$ should influence the convergence behaviour of GMRES (see Section 2.3). To obtain computable bounds on the residual norm that generate a good prescribed convergence curve for a general matrix is a difficult challenge. Theoretical results (min-max approximations on matrix eigenvalues) in this field are still partial. We see in the sequel of the chapter that the analysis of the distribution of the eigenvalues in the complex plane gives a very useful approach for predicting the convergence of the GMRES. Concerning the restarted GMRES, the following proposition holds.

**Theorem 2.** If $A$ is symmetric positive-definite, then GMRES($m$) converges at any $m \ge 1$.

However, the restarted GMRES algorithm can stagnate when the matrix is not definite positive. We will observe this difficulty in the following section.
## 2.3 Convergence problems for indefinite linear systems

In iterative methods, a common belief is that the condition number $\kappa_2(A)$ (in the 2-norm) of A is a good measure of the convergence rate of the algorithm. This is generally not true for a complex-valued matrix A where indeed the distribution of the eigenvalues in the complex plane is the most crucial point to observe. In many applications, the eigenvalues of the matrix A are computed in order to examine if the whole spectrum is included or not in a given part of the complex plane. Unfortunately, the matrix is often defined with a given precision and the computed eigenvalues may differ from the real ones, especially in highly non normal cases. To answer to this problem, the difficult notion of $\epsilon$-pseudospectrum of a matrix was introduced [65]. In a few words, the idea is to compute the set of the eigenvalues of perturbed matrices $A + E$ for some $E$, with $\|E\| \le \epsilon$. Generally, this approach is costly in terms of computation.

To illustrate the fact that the distribution of the eigenvalues plays a crucial role in the convergence of the GMRES, let us consider the seven simple complex-valued diagonal matrices $A_j$ defined by

$$ 
\begin{aligned}
A_1 &= \operatorname{diag}(1/\ell^2)_{1\le\ell\le n}, & A_2 &= \operatorname{diag}(e^{i\ell}/\ell^2)_{1\le\ell\le n}, & A_3 &= \operatorname{diag}(e^{i\ell})_{1\le\ell\le n}, \\
A_4 &= \operatorname{diag}(1+e^{i\ell})_{1\le\ell\le n}, & A_5 &= \operatorname{diag}(1.5+e^{i\ell})_{1\le\ell\le n}, & \\
A_6 &= \operatorname{diag}(1+e^{i\ell}/\ell)_{1\le\ell\le n}, & A_7 &= \operatorname{diag}(1+e^{i\ell}/\ell^2)_{1\le\ell\le n}.
\end{aligned} 
\quad (5) 
$$

Concerning $A_1$, the matrix is real-valued and the eigenvalues tend towards zero as $\ell \to +\infty$ (Figure 1(a)). In terms of operators (see section 3), it corresponds to the notion of real-valued elliptic positive pseudodifferential operator of order $-2$. More or less, the underlying operator is the inverse of the one-dimensional Laplacian operator $(\partial_x^2)^{-1}$. In particular, the matrix is symmetric and *positive-definite*. We can observe that there is an eigenvalue clustering around zero. The condition number $\kappa_2(A_1) = n^2$ becomes large as $n$ grows and so the convergence of GMRES(50) takes more iterations (Figure 1(b)). However, the convergence is observed. The tolerance of the iterative solver is fixed equal to $\epsilon = 10^{-14}$ in all the examples.

The second test-case is related to a complex-valued matrix $A_2$. The case is close to the previous one but the distribution of the eigenvalues in the complex plane is completely different even if their modulus tends to zero (Figure 1(c)). The condition number is again $n^2$ but we can quickly observe the divergence of GMRES(50) e.g. for $n = 100$ (Figure 1(d)). Here, the matrix is indefinite (eigenvalues with negative or positive real parts). In terms of operators, this corresponds to an indefinite complex-valued integral (or pseudodifferential) operator of order $-2$. This is also a first-kind Fredholm integral operator (see Section 4.3.5). These two situations show how the convergence of the GMRES strongly depends on the distribution of the eigenvalues in the complex plane.

The third situation, with $A_3$, is also very interesting. The complex-valued matrix has a condition number equals to 1 but the convergence of GMRES(50) takes many iterations to converge for $n = 100$ and even diverges for $n = 1000$. This example
illustrates clearly that considering indefinite matrices leads to instabilities of the GMRES and eventually to its divergence. This is still true if one translates the spectrum of $A_3$ from 1 to the right and get the matrix $A_4$ (Figure 2(a)). Now, if one translates $A_3$ from 1.5 and obtain $A_5$, then the convergence holds and is pretty fast since all the eigenvalues have a large positive real part and are sufficiently far away from the origin (Figure 2(b)).

Let us consider the two matrices $A_6$ and $A_7$. They look like $A_4$ but the perturbative term to 1 has a modulus which tends towards zero linearly or quadratically (Figures 2(c) and 2(e)). This rate of convergence of the sequence is related to the order of the underlying operator: first- or second-order. This kind of matrices (close to positive-definite matrices for large values of $n$) corresponds to an integral operator called a second-kind Fredholm operator (see Section 4.3.5). Their spectrum clusters around a complex value $\alpha$ (= 1 here) up to a sequence converging to zero ($e^{i\ell}/\ell$ or $e^{i\ell}/\ell^2$ here). These configurations lead to converging iterative schemes with a rate depending on the decay of the sequence to zero (Figures 2(d) and 2(f)).

## 2.4 Why this happens in acoustic scattering: a simple example

We have just seen that an indefinite matrix without any robust eigenvalue clustering implies bad convergence properties of the GMRES. Such a situation naturally arises in acoustic scattering.

Let us look at the following simple scattering problem. We consider a plane wave $u^{\text{inc}}(x) = e^{-ikx}$ coming from $+\infty$ and illuminating the left positive domain $]-\infty; 0]$. The real-valued constant wavenumber is denoted by $k$. The sound-hard scattering problem reads: find the scattered field $u$ in the unbounded domain $[0; +\infty[$ such that

$$ \left\{ \begin{array}{ll} (\partial_x^2 + k^2)u = 0, & \text{in } [0; +\infty[, \\ \partial_x u(0) = ik, & \\ u \text{ travels to the right.} & \end{array} \right. \qquad (6) $$

Of course, the solution is trivial: $u(x) = e^{ikx}$. Let us note that other boundary conditions could be used as for example the sound-soft one: $u(0) = -1$ or an impedance boundary condition. Since problem (6) is set in an unbounded domain, one usually introduces an Artificial Boundary Condition (ABC) at the fictitious boundary $\Sigma = \{1\}$ (another point could be chosen). Its aim is to replace the sentence "u travels to the right" to get a bounded domain boundary value problem. Here, the boundary condition is trivial and is $(\partial_n - ik)u = 0$, where $n$ is the outwardly directed unit normal vector to $\Omega$ at $\Sigma$. Finally, the problem which is solved by a finite element method is

$$ \left\{ \begin{array}{ll} (\partial_x^2 + k^2)u = 0, & \text{in } \Omega, \\ \partial_x u(0) = ik, & \\ (\partial_n - ik)u = 0, & \text{on } \Sigma, \end{array} \right. \qquad (7) $$

with $\Omega = ]0; 1[$ and $\Sigma = \{1\}$. Let us consider that we use a linear continuous Galerkin
<figure><img src="image_2.png" /><figcaption>(a) $A_1$: Eigenvalue distribution ($n = 100$)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) $A_2$: Eigenvalue distribution ($n = 100$)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(e) $A_3$: Eigenvalue distribution ($n = 1000$)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) $A_1$: Convergence history</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) $A_2$: Convergence history</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(f) $A_3$: Convergence history</figcaption></figure>

<caption>Figure 1: Matrices $A_1$, $A_2$, and $A_3$.</caption>
<figure><img src="image_8.png" /><figcaption>(a) $A_4$: Convergence history</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) $A_5$: Convergence history</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) $A_6$: Eigenvalue distribution ($n = 100$)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) $A_6$: Convergence history ($n = 100$)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) $A_7$: Eigenvalue distribution ($n = 100$)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) $A_7$: Convergence history</figcaption></figure>

Figure 2: Matrices $A_4$, $A_5$, $A_6$ and $A_7$.
approximation of (7). Then, the variational formulation writes down: find $u$ such that

$$ \int_{\Omega} \partial_x u \partial_x v - k^2 u v dx - iku(1)v(1) = -ikv(0), \qquad (8) $$

for all well-chosen test-functions $v$. Let us denote by $M_h$ and $S_h$ respectively the mass and stiffness matrices associated with the linear Finite Element Method (FEM), for a uniform discretization of $\Omega$. The length of one element is $h$ and the total number of degrees of freedom of the FEM is $n_h$ (which is equal to the number of segments plus one). We denote by $u_h \in C^{n_h+1}$ the linear approximation of $u$ solution to (7) or (8). The term $-iku(1)v(1)$ related to the ABC contributes in the linear system by an additional term on the last row and column of the linear system. We denote this matrix term by $-ikB_h$. The right-hand side is related to $-ikv(0)$ and gives birth to a right-hand side vector $b_h \in C^{n_h+1}$. Finally, the linear finite element approximation of (8) leads to the solution of the system

$$ (S_h - k^2 M_h - ikB_h)u_h = b_h. \qquad (9) $$

The sparse matrix involved in system (9) is $i$) complex-valued because of the boundary term $-ikB_h$ and $ii$) non positive-definite since we have the contribution $S_h - k^2 M_h$. This is most particularly penalizing when solving high frequency problems ($k$ large). We report on Figure 3(a) the behaviour of GMRES(50) applied to solving (9) for $k = 60$ and $n_h = 100$. The convergence is extremely slow. Let us note that $n_h$ must be quite large for this value of $k$ because of the pollution error [15, 64] into the FEM. Figure 3(b) shows the distribution of the eigenvalues of $(S_h - k^2 M_h - ikB_h)$ as well as $(S_h - k^2 M_h)$. We can clearly see that the problem has many eigenvalues lying in the left half-plane with null imaginary part for $(S_h - k^2 M_h)$. As a consequence, the matrix is non positive-definite. Furthermore, the reason why all the eigenvalues are real is that we are rather solving an interior Helmholtz problem with Neumann boundary condition. In the case of an (exterior) scattering problem, adding the ABC (which means that we consider the additional term $-ikB_h$) leads to a complex spectrum. This can be observed on Figure 3(b). Then, this is worst for having a converging iterative scheme.

## 2.5 How to improve the convergence: preconditioners

For practical purpose, large scale sparse linear systems are solved by an iterative Krylov solver. This is the case in acoustic scattering where we must be able to design 1) a convergent iterative solver and 2) methods that converge fastly if convergence occurs. One very convenient way to do this is to preconditioning the linear system by a sparse matrix $P$ called the preconditioner. The idea is the following. Let us assume that we are solving a linear system $Ax = b$ and that we are able to build a matrix $P$ such that $P$ is close to $A^{-1}$. This means that in some sense, for a suitable matrix norm $\|\cdot\|$, the quantity $\|PA - I\|$ or $\|AP - I\|$ is small. Here, $I$ is the identity matrix. We say that $P$ is a left preconditioner if $PA \approx I$ or a right preconditioner if $AP \approx I$. Since
<figure><img src="image_14.png" /><figcaption>Figure 3: Matrix $(S_h - k^2 M_h - ikB_h)$.</figcaption></figure>

the condition number $\kappa(PA)$ is close to $\kappa(I) = 1$, we can think that the solution to the new linear system $PAx = Pb$ by an iterative method is convergent and fast. This can be expected for positive-definite matrices $A$ but not necessarily for indefinite matrices as we previously noticed. If applying $P$ only requires a Matrix-Vector Product (MVP) in the GMRES, we say that $P$ is an *explicit* preconditioner. If its application needs the solution of a sparse linear system, $P$ is said to be *implicit*.

Since we work at a matrix level, then the construction of the preconditioner $P$ can only be based on *algebraic* considerations (the entries of $A$). Many algebraic preconditioners have been proposed during the last decades. Let us mention e.g. Incomplete LU preconditioners (ILU), SParse Approximate Inverses (SPAI) or Algebraic Recursive Multilevel Solvers (ARMS) and all their variants. We refer to [17] for some of these techniques for general linear systems. Concerning the solution of scattering problems, these preconditioners have been tested and improved for instance in [18, 26, 27, 41, 44, 45, 56, 67]. However, even if these preconditioners provide an improved convergence, convergence breakdown and slow convergence still arise when medium and high wavenumbers $k$ are considered. This challenging problem is so still open.

An alternative to the algebraic preconditioners is developed in the sequel. Essentially, we can summarize the derivation of the algebraic preconditioners to: 1) first, take a continuous (partial differential or integrodifferential) operator, 2) discretize it to get a matrix $A$ and 3) build an approximate algebraic inverse $P$. Our solution leading to what we call *analytical* preconditioners is based on the following three points: 1) take a partial differential operator or an integral operator, 2) build an approximate inverse at the continuous level and 3) discretize this operator for example by a finite element method to get a preconditioner $P$. The interesting aspect of this approach is that point 2) keeps the information related to the underlying structure of the operator to approximate, information which is *a priori* lost by considering the algebraic viewpoint. To attain this goal, we need a mathematical theory which allows to compute
practically an approximation of the inverse of a general operator at the continuous level. This is the aim of the next section where we introduce the tools from pseudodifferential operators theory and associated microlocal calculus.

# 3 Elements on pseudodifferential operator theory

We introduce in this section the basics of the theory of pseudodifferential operators. The aim of pseudodifferential operators and associated microlocal symbolic calculus is to naturally generalize partial differential and integral operators through the notion of symbol. After the definitions, we give some practical rules for symbolic calculus which are the keystone for building approximations of pseudodifferential operators, with future applications to analytical preconditioners (Sections 3.3, 7 and 8). For further reading, we refer e.g. to [63] where the theory of pseudodifferential operators is presented with more details that we cannot address in this short introduction.

## 3.1 Definitions: pseudodifferential operator and symbol

Let $\Omega \subset \mathbb{R}^d$ be an open set and $D'(\Omega)$ the space of distributions on $\Omega$ which is the dual space of $D(\Omega) := C_0^\infty(\Omega)$ [60, 63]. We introduce the vectorial differential operator $D = (D_1, \dots, D_d)$, setting $D_j := -i\partial_j = -i\partial_{x_j}$. A variable coefficients partial differential operator $P(x, D)$ of order $m$ has the general form

$$P(x, D) = \sum_{|\alpha| \le m} a_{\alpha}(x) D^{\alpha}, \qquad (10)$$

where $x := (x_1, \dots, x_d) \in \Omega$, $\alpha = (\alpha_1, \dots, \alpha_d)$ is a multi-index in $\mathbb{N}^d$ and the coefficients $a_\alpha$ are $C^\infty(\Omega)$ smooth functions. The operator $D^\alpha$ is: $D^\alpha = D_1^{\alpha_1} \dots D_d^{\alpha_d}$. For the sake of clarity, we will sometimes precise the derivation variable used for $D^\alpha$ like for example $D_x^\alpha$ or $\partial_\xi^\alpha$. The polynomial $p: \Omega \times \mathbb{R}^d \to \mathbb{C}$

$$p(x, \xi) = \sum_{|\alpha| \le m} a_{\alpha}(x) \xi^{\alpha},$$

is called the symbol of the operator $P$ (in the sequel $\xi = (\xi_1, \dots, \xi_d)$ is the Fourier covariable). The *principal symbol* of order $m$ of $P$, denoted by $\sigma_p(P)$, represents the homogeneous part of degree $m$ in $\xi$ of $p$,

$$\sigma_p(P)(x, \xi) = \sum_{|\alpha|=m} a_{\alpha}(x) \xi^{\alpha}.$$

For example, for the operator $P(x, D) := \text{div}(A_2(x)\nabla\cdot) + a_0(x)$, we have $p(x, \xi) = -\xi^T A_2(x)\xi + a_0(x)$ and $\sigma_p(P)(x, \xi) = -\xi^T A_2(x)\xi$. Here, $\text{div } v := \partial_1 v_1 + \dots + \partial_d v_d$ is the usual divergence operator of a vector field $v$.
Let $f$ be a function in the Schwartz space $S$ of $C^\infty$ functions that rapidly decay at infinity. Then, its Fourier transform $\hat{f} \in S$ is defined by

$$ \hat{f}(\xi) = \int_{\mathbb{R}^d} e^{-i\boldsymbol{x} \cdot \boldsymbol{\xi}} f(\boldsymbol{x}) d\boldsymbol{x}, $$

and we have

$$ f(\boldsymbol{x}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{i\boldsymbol{x} \cdot \boldsymbol{\xi}} \hat{f}(\boldsymbol{\xi}) d\boldsymbol{\xi}. $$

It is well-known that the Fourier transform is an isomorphism from $S$. Moreover we have the property $\tilde{D}^\alpha u(\boldsymbol{\xi}) = \boldsymbol{\xi}^\alpha \hat{u}(\boldsymbol{\xi})$. The inverse Fourier transform is used to rewrite the partial differential operator $P: S \to S$ through its symbol $p$

$$ P(\boldsymbol{x}, D)u(\boldsymbol{x}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{i\boldsymbol{x} \cdot \boldsymbol{\xi}} p(\boldsymbol{x}, \boldsymbol{\xi}) \hat{u}(\boldsymbol{\xi}) d\boldsymbol{\xi}, $$

with $u \in S$. Pseudodifferential operators are a generalization of differential operators. The motivation is to replace polynomial symbols $p$ in $\boldsymbol{\xi}$ with more general symbols. Working not only locally in space for $\boldsymbol{x}$ but also in the cotangent space with respect to $\boldsymbol{\xi}$ is known as *microlocal analysis*. Let us introduce the space of admissible symbols that define a pseudodifferential operator.

**Definition 2.** Let $m \in \mathbb{R}$. We denote by $S^m(\Omega)$ the space of functions $a \in C^\infty(\Omega \times \mathbb{R}^d)$ such that for every compact subset $K \subset \Omega$ and every $\alpha \in \mathbb{N}^d$, $\beta \in \mathbb{N}^d$, there exists a constant $C = C(K, \alpha, \beta) \in \mathbb{R}$ such that

$$ |\partial_{\boldsymbol{x}}^{\beta} \partial_{\boldsymbol{\xi}}^{\alpha} a(\boldsymbol{x}, \boldsymbol{\xi})| \leq C(1 + |\boldsymbol{\xi}|)^{m-|\alpha|}, \quad \forall (\boldsymbol{x}, \boldsymbol{\xi}) \in K \times \mathbb{R}^d. $$

The notation $|\boldsymbol{\xi}|$ designates the euclidian norm of vector $\boldsymbol{\xi} \in \mathbb{R}^d$, i.e. $|\boldsymbol{\xi}| = \sqrt{\boldsymbol{\xi} \cdot \boldsymbol{\xi}}$, and, for a multi-index, we set $|\alpha| = \alpha_1 + \dots + \alpha_d$. Elements of $S^m(\Omega)$ are called symbols of order $m$ and we write $a \in S^m$.

Let us define now a pseudodifferential operator $A$ of order $m$ through its symbol $a \in S^m$ and the inverse Fourier transform.

**Definition 3.** A symbol $a \in S^m(\Omega)$ defines a continuous linear operator $A = \text{Op } a : C_0^\infty(\Omega) \to C_0^\infty(\Omega)$ by

$$ A(\boldsymbol{x}, D)u(\boldsymbol{x}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{i\boldsymbol{x} \cdot \boldsymbol{\xi}} a(\boldsymbol{x}, \boldsymbol{\xi}) \hat{u}(\boldsymbol{\xi}) d\boldsymbol{\xi}. $$

Operators of this type are called *pseudodifferential operators of symbol a and of order m*. We set then $A \in \psi^m(\Omega)$.

Pseudodifferential operators of arbitrarily small order $m$ are called regularizing and we have $\psi^{-\infty}(\Omega) = \bigcap_{m \in \mathbb{R}} \psi^m(\Omega)$. The theory of pseudodifferential operators offers an interesting and useful property through the symbolic calculus. This leads to practical computations on operators working at the algebraic level with their symbols (see section 3.2). To complete some of the definitions below, we introduce the notion of elliptic pseudodifferential operators.
**Definition 4.** A symbol $a \in S^m(\Omega)$ is called elliptic (of degree $m$) if for every $K \subset \Omega$ compact, there exists a constant $C$ such that

$$|a(\boldsymbol{x}, \boldsymbol{\xi})| \geq C(1 + |\boldsymbol{\xi}|)^m, \quad \forall \boldsymbol{x} \in K, |\boldsymbol{\xi}| \geq \frac{1}{C}.$$

An elliptic pseudodifferential operator $A \in \psi^m(\Omega)$ is an operator with an elliptic symbol $a \in S^m(\Omega)$.

Pseudodifferential operators have regularity properties which are related to their order $m$. For example, if we denote by $H^s(\Omega)$ the Sobolev space of order $s \in \mathbb{R}$ of distributions $u$ defined on $\Omega$ [63], then it can be proved that $A$ is a continuous linear operator acting from $H^s(\Omega)$ onto $H^{s-m}(\Omega)$, for any $s \in \mathbb{R}$.

Finally, pseudodifferential operators have the pseudolocal property. We say that an operator $A$ acting on a distribution $u$ is local if $Au$ is smooth in the same set as $u$. Pseudolocal means that the set where $A$ is smooth includes the set where $u$ is smooth. This implies that $A$ could smooth out a nonsmoothness of $u$ (a more rigorous mathematical definition uses the notion of support and singular support of a distribution [63]). Partial differential operators with smooth coefficients are local operators and every local operator is a differential operator. Examples of pseudolocal operators include integral operators $A$ of the form

$$Au(\boldsymbol{x}) = \int_{\Omega} G(\boldsymbol{x}, \boldsymbol{y})u(\boldsymbol{y})d\boldsymbol{y},$$

where $G$ is a smooth kernel.

## 3.2 Practical symbolic calculus rules

The most important property of pseudodifferential operators for practical applications in scientific computing is the fact that all the computations like the composition of two operators or the transposition of an operator can be performed algebraically at the symbol level. We give here the main results useful for the sequel and most particularly introduce the idea of asymptotic expansion of a symbol.

A symbol $a \in S^m(\Omega)$ is said to be homogeneous of degree $m$ in $\boldsymbol{\xi}$ if the following relation holds

$$\forall (\boldsymbol{x}, \boldsymbol{\xi}) \in \Omega \times (\mathbb{R}^d \setminus 0), \forall \lambda > 0, a(\boldsymbol{x}, \lambda\boldsymbol{\xi}) = \lambda^m a(\boldsymbol{x}, \boldsymbol{\xi}).$$

We can then give the following definition.

**Definition 5.** Let $a \in S^m(\Omega)$. Consider a decreasing sequence of real constants $(m_j)_{j=1}^\infty$ with $\lim_{j \to \infty} m_j = -\infty$. Let $(a_j)_{j=1}^\infty$ a sequence of homogeneous symbols $a_j \in S^{m_j}(\Omega)$ such that

$$a - \sum_{j=1}^{k} a_j \in S^{m_{k+1}}(\Omega), \qquad (11)$$
for every $k \in \mathbb{N}$, $k \neq 0$. Then, we say that $(a_j)_{j=1}^{\infty}$ is an asymptotic expansion of a. In this case, we write $a \sim \sum_j a_j$. The first term $a_1$ of order $m_1$ in the expansion is called the principal symbol.

Not every symbol $a \in S^m(\Omega)$ has an asymptotic expansion. The set of symbols of the form (11) is called classical. It is denoted by $S_{cl}^m(\Omega)$ and the corresponding operators belong to $\psi_{cl}^m(\Omega)$.

**Remark 2.** If two symbols $a$ and $b$ have the same asymptotic expansion, then they differ from a smoothing pseudodifferential operator

$$a - b = (a - \sum_{j=1}^{k} a_j) - (b - \sum_{j=1}^{k} a_j) \in S^{m_{k+1}}(\Omega),$$

for all $k$ and $\lim_{j \to +\infty} m_j = -\infty$, so $(a-b) \in S^{-\infty}(\Omega)$, with $S^{-\infty}(\Omega) = \bigcap_{m \in \mathbb{R}} S^m(\Omega)$.

As said above, one of the crucial points of pseudodifferential operators is that we have algebraic rules for computing some operations on pseudodifferential operators. Two extremely important properties are the following.

**Proposition 3.** Let $A \in \psi_{cl}^{m_1}(\Omega)$, $B \in \psi_{cl}^{m_2}(\Omega)$ with symbols $a \in S_{cl}^{m_1}(\Omega)$, $b \in S_{cl}^{m_2}(\Omega)$ respectively. Then

1. The transpose $A^t$ of A is a pseudodifferential operator. The symbol $a^t \in S_{cl}^{m_1}(\Omega)$ of $A^t \in \psi_{cl}^{m_1}(\Omega)$ is given by the following expansion

$$a^t(\mathbf{x}, \boldsymbol{\xi}) \sim \sum_{\alpha \in \mathbb{N}^d} \frac{1}{\alpha!} \partial_{\boldsymbol{\xi}}^{\alpha} D_{\mathbf{x}}^{\alpha} a(\mathbf{x}, -\boldsymbol{\xi}), \quad (12)$$

where $\alpha \in \mathbb{N}^d$ is a multi-index.

2. The composition of two operators $A$ and $B$, denoted by $AB$, is a pseudodifferential operator. The symbol $a \sharp b \in S_{cl}^{m_1+m_2}(\Omega)$ of $AB \in \psi_{cl}^{m_1+m_2}(\Omega)$ is

$$a^{\sharp} b(\mathbf{x}, \boldsymbol{\xi}) \sim \sum_{\alpha \in \mathbb{N}^d} \frac{1}{\alpha!} \partial_{\boldsymbol{\xi}}^{\alpha} a(\mathbf{x}, \boldsymbol{\xi}) D_{\mathbf{x}}^{\alpha} b(\mathbf{x}, \boldsymbol{\xi}). \quad (13)$$

In particular, this shows that $a(\mathbf{x}, -\boldsymbol{\xi})^t$ is the principal symbol of $A^t$ and $\sigma_p(A)\sigma_p(B)$ of $AB$.

In addition to Proposition 3, let us remark that $D_{\mathbf{x}}^{\alpha} a \in S_{cl}^{m_1}$ and $\partial_{\boldsymbol{\xi}}^{\alpha} a \in S_{cl}^{m_1 - |\alpha|}$ if $a \in S_{cl}^{m_1}$. Furthermore, we have the following Theorem.

**Theorem 3.** Let $A$ be an elliptic pseudodifferential operator of order $m$. There exists a pseudodifferential operator $B$ (inverse of $A$) with order $-m$ such that $AB - I$ (right inverse of $A$) and $BA - I \in \psi^{-\infty}(\Omega)$ (left inverse of $A$). The operator $I$ is the identity operator (with symbol 1).
We can see that combining formula (13) and Theorem 3 gives a practical way of computing an approximate inverse B of a given pseudodifferential operator A if we know its symbol a or at least the first terms of its asymptotic expansion. In practice, this can often be done. This point of view is at the basis of what we call *analytical preconditioners* which is an alternative approach to the purely *algebraic* methods described briefly in section 2.5.

### 3.3 A first and simple application to preconditioning of the Helmholtz equation

An example of application of pseudodifferential operators theory and symbolic calculus rules to preconditioning is the following. Let us consider the Helmholtz operator $L = \Delta + k^2$ of symbol $\sigma_L(x, \xi) = \sigma_L(\xi) = k^2 - |\xi|^2$ (the dependence in $x$ really occurs for an inhomogeneous media). Then, an exact analytical preconditioner would be an operator $A$ such that $AL = I$. Since $L$ is a pseudodifferential (in fact, partial differential) operator of order 2, $A$ is a pseudodifferential operator of order $-2$. Formally, $A$ is equal to the nonlocal operator $(\Delta + k^2)^{-1}$ with symbol $(k^2 - |\xi|^2)^{-1}$. This point of view is purely theoretical since the practical computation of $A$ is exactly what we wish to obtain by solving the initial scattering problem. A first approximation is to consider a static approximation of the symbol of $A$ as: $\sigma_A \approx (-|\xi|^2)^{-1}$. This means that the corresponding analytical preconditioner is $A_0 = Op((-|\xi|^2)^{-1}) = \Delta^{-1}$, implying that the associated preconditioned operator to solve is: $A_0L = \Delta^{-1}(\Delta + k^2) = I + k^2\Delta^{-2}$. Hence $A_0L$ is a second-kind integral operator (see Section 4.3.5) with eigenvalue clustering around $(1, 0)$ in the complex plane for large frequencies $|\xi|$. This idea was introduced differently by Bayliss, Goldberg and Turkel in [16]. It can be shown (see also Figure 4(a) for a one-dimensional example) that this clearly improves the convergence of an iterative solver when considering low-frequencies $k$ (close to the static problem). However, for larger values of $k$ (medium and high frequency regimes), the convergence may fail or strongly depends on $k$. As an example, let us consider $k = 60$ and $n_h = 80$ for the one-dimensional case. We represent on Figure 4(a) the eigenvalue distribution in the complex plane of the corresponding discrete preconditioned matrix $\mathbb{I} - k^{-2}\mathbb{S}_h^{-1}\mathbb{M}_h$. If we would zoom around the origin, then one would observe a clustering of the eigenvalues around $(1, 0)$. However, as we can see, many eigenvalues remain in the left half-plane leading to an indefinite matrix when $k$ is sufficiently large. An improved solution is to consider the smoothing of $A$ by a complex parameter $\alpha \in \mathbb{C}$: $A_\alpha = (\Delta + k^2\alpha)^{-1}$. Adding this parameter leads to the solution of a dissipative Helmholtz equation. This approach is called Laplace shifted preconditioning approach [38]. In terms of eigenvalues, the symbol of the preconditioned operator $A_\alpha L$ is

$$ \sigma_{A_\alpha L}(x, z) = \frac{1 + z}{\alpha + z}, \qquad (14) $$
setting $z = -|\xi|^2/k^2 \in \mathbb{R}^-$. The eigenvalues are then distributed along a circular arc as a consequence of the Poincaré map (14). The center of the circle is

$$ (\frac{1}{2}, \frac{(a^2 + b^2 - a)}{2b}) $$

and its radius is

$$ \sqrt{\frac{1}{4} + \frac{(a - (a^2 + b^2))^2}{4b^2}}, $$

setting $\alpha = a+ib$. The endpoints are $\alpha^{-1}$ for the low-frequency spatial values $|\xi| \ll k$ ($z \approx 0$ for physical propagative modes) and $(1,0)$ for large frequencies $|\xi| \gg k$ ($z \approx -\infty$ for evanescent waves). For the "grazing" waves related to $|\xi| \approx k$ ($z \approx -1$), one gets positive eigenvalues close to the origin. These remarks can be observed on Figures 4(b) and 4(c) for respectively $\alpha = 1 + 0.5i$ and $1 + i$ ("analytical arc"). From a practical point of view, $A_\alpha$ is computed by a few step of a multigrid algorithm or by an ILU factorization and not a direct solver (see [38, 66]). These preconditioners a priori lead to preconditioned matrices with eigenvalues having positive real parts (and so positive-definite matrices correspond). However, for a scattering problem, a boundary contribution related to the ABC must be considered in the global system matrix and must be included into the preconditioner. On figures 4(b) and 4(c), we draw the numerical eigenvalues of the preconditioned matrices $A_{\alpha,h}^{ik} \mathbb{L}_h$ with ABCs, setting $\mathbb{L}_h = S_h - k^2 M_h - ik B_h$ and $A_{\alpha,h}^{ik} = (S_h - k^2 \alpha M_h - ik B_h)^{-1}$. We observe that the introduction of the boundary term modifies the circular arc. Most particularly, the associated spectrum has eigenvalues with negative real parts, meaning that $A_{\alpha,h}^{ik} \mathbb{L}_h$ is an indefinite matrix now. The residual history of the GMRES(50) for computing the solution to the scattering problem with different preconditioners are reported on Figure 4(d). This shows the divergence of the restarted GMRES without preconditioner and its convergence for $\alpha = 0, \alpha = 1 + i$ and $\alpha = 1 + 0.5i$, which is compatible with the literature on the subject [38, 66]. The fastest convergence is obtained with $\alpha = 1 + 0.5i$, which can be expected from Figure 4(c) since most eigenvalues have significative positive real parts.

This first simple application shows that the pseudodifferential operator theory and symbolic calculus lead to some possible constructions of robust and efficient preconditioners for PDEs. The aim of the next section is to provide other applications for preconditioning well-known integral equation formulations used in acoustic scattering problems.

# 4 Potential theory - integral equations

In this section, we present how to solve an exterior boundary value problem set in an unbounded computational domain via the integral equations method. First, we recall the basic acoustic scattering problem and some notations. Next, we give elements of potential theory that are crucial for the integral equations method. Finally, we discuss the derivation of the classical direct and indirect integral equations for both Dirichlet
<figure><img src="image_15.png" /><figcaption>(a) $\alpha = 0$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) $\alpha = 1 + i$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) $\alpha = 1 + 0.5i$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Convergence history</figcaption></figure>

<caption>Figure 4: Laplace shifted preconditioning.</caption>
and Neumann boundary conditions. We also describe their properties in view of a numerical solution by a Krylov iterative solver (like the GMRES).

## 4.1 Acoustic scattering problems

Let us define a $d$-dimensional bounded domain $\Omega^- \subset \mathbb{R}^d$ representing an impenetrable body with boundary $\Gamma := \partial\Omega^-$. We denote by $\Omega^+ := \mathbb{R}^d \setminus \bar{\Omega}^-$ the associated homogeneous exterior domain of propagation. Consider the scattering of an incident time-harmonic acoustic wave $u^{\text{inc}}$ by the obstacle $\Omega^-$. The scattered field $u^+$ satisfies the following exterior boundary-value problem

$$ \left\{ \begin{array}{ll} \Delta u^+ + k^2 u^+ = 0, & \text{in } \Omega^+, \\ u^+|_{\Gamma} = -u^{\text{inc}}|_{\Gamma} \text{ or } \partial_n u^+|_{\Gamma} = -\partial_n u^{\text{inc}}|_{\Gamma}, & \text{on } \Gamma, \\ \lim_{\|x\| \to +\infty} \|x\|^{(d-1)/2} (\nabla u^+ \cdot \frac{x}{\|x\|} - iku^+) = 0, & \end{array} \right. \quad (15) $$

setting $\Delta := \sum_{i=1}^d \partial_{x_j}^2$. We consider an incident time-harmonic plane wave of the form:

$$ u^{\text{inc}}(x) = e^{-ik\theta^{\text{inc}} \cdot x}. $$

This wavefield is characterized by the wavenumber $k := 2\pi/\lambda$, setting $\lambda$ as the wavelength of the signal. In the two-dimensional case ($d = 2$), the direction of incidence $\theta^{\text{inc}}$ is given by the relation $\theta^{\text{inc}} = (\cos(\theta^{\text{inc}}), \sin(\theta^{\text{inc}}))^T$, where $\theta^{\text{inc}}$ is the scattering angle in the polar coordinates system. In the three-dimensional case ($d = 3$), we have $\theta^{\text{inc}} = (\cos(\theta^{\text{inc}}) \sin(\phi^{\text{inc}}), \sin(\theta^{\text{inc}}) \sin(\phi^{\text{inc}}), \cos(\phi^{\text{inc}}))^T$. The scattering angles $(\theta^{\text{inc}}, \phi^{\text{inc}})$ are expressed in the spherical coordinates system. We define by $n$ the outwardly directed unit normal to $\Omega^-$ at the boundary $\Gamma$. The boundary condition on $\Gamma$ (second equation of (15)) depends on the physical problem under study. The sound-soft or Dirichlet (respectively sound-hard or Neumann) boundary condition on $\Gamma$ corresponds to the first (respectively second) boundary condition in (15). Finally, the last equation is the well-known Sommerfeld radiation condition or the outgoing wave condition. This condition completely characterizes the behaviour of the solution to the Helmholtz equation at infinity and guarantees the uniqueness of the solution to the exterior problem (15).

Let us introduce the functional spaces [55]

$$ \begin{aligned} H_{\text{loc}}^{s}(\overline{\Omega^{+}}) &:= \{v \in \mathcal{D}'(\Omega^{+}) / \psi v \in H^{s}(\Omega^{+}), \forall \psi \in \mathcal{D}(\mathbb{R}^{d})\}, \quad s \ge 1, \\ H_{-}^{1}(\Delta) &:= H^{1}(\Delta, \overline{\Omega^{-}}) := \{u \in H^{1}(\overline{\Omega^{-}}); \Delta u \in L^{2}(\overline{\Omega^{-}})\}, \\ H_{+}^{1}(\Delta) &:= H_{\text{loc}}^{1}(\Delta, \overline{\Omega^{+}}) := \{u \in H_{\text{loc}}^{1}(\overline{\Omega^{+}}); \Delta u \in L_{\text{loc}}^{2}(\overline{\Omega^{+}})\}. \end{aligned} $$

For $u \in H_{\pm}^{1}(\Delta)$, the exterior (+) and interior (-) trace operators of order $j$ ($j = 0$ or 1) can be defined by

$$ \begin{array}{rcl} \gamma_j^{\pm} : H_{\pm}^{1}(\Delta) & \to & H^{1/2-j}(\Gamma) \\ u & \mapsto & \gamma_j^{\pm} u^{\pm} = \partial_n^j u^{\pm}|_{\Gamma}. \end{array} \quad (16) $$
In this functional setting, the existence and uniqueness of the solution to the scattering problem

$$
\left\{
\begin{array}{l}
\text{Find } u^+ \in H_{\text{loc}}^1(\bar{\Omega}^+) \text{ such that} \\
\Delta u^+ + k^2 u^+ = 0, \quad \text{in } \mathcal{D}'(\Omega^+), \\
\gamma_j^+ u^+ = g := -\gamma_j^+ u^{\text{inc}}, \quad \text{in } H^{1/2-j}(\Gamma), \quad \text{for } j = 0 \text{ or } 1, \\
\lim_{\|x\| \to +\infty} \|x\|^{(d-1)/2} (\nabla u^+ \cdot \frac{x}{\|x\|} - iku^+) = 0,
\end{array}
\right.
\quad (17)
$$

can be proved [33].

The first main difficulty arising in the numerical solution of the exterior boundary-value problem (17) is related to the unboundedness of the computational domain $\Omega^+$. A solution is to apply the integral equations method [29, 33]. This approach allows to reformulate the initial boundary-value problem as an integral equation defined on the boundary $\Gamma$ of the scattering obstacle $\Omega^-$. Then, this method reduces the dimension of the problem to $d-1$. Boundary integral equations are derived from potential theory. Let us give in the following Section some elements of this theory.

## 4.2 Potential theory: basic relations - properties

The essential property is that any solution to the Helmholtz equation can be represented as the linear combination of a single- and a double-layer potentials. The following proposition holds.

**Proposition 4.** Let us define the outgoing Green's function $G$ associated with the Helmholtz operator in $\mathbb{R}^d$ by

$$ G(\mathbf{x}, \mathbf{y}) = \begin{cases} \frac{i}{4} H_0^{(1)}(k \|\mathbf{x} - \mathbf{y}\|), & \text{for } d=2, \\ \frac{1}{4\pi} \frac{e^{ik\|\mathbf{x} - \mathbf{y}\|}}{\|\mathbf{x} - \mathbf{y}\|}, & \text{for } d=3, \end{cases} \quad (18) $$

where $H_0^{(1)}$ designates the first-kind Hankel function of order zero. Let $(v^{-}, v^{+}) \in H^{1}(\Omega^{-}) \times H_{\text{loc}}^{1}(\bar{\Omega}^{+})$ satisfying

$$ \Delta v^{-} + k^{2}v^{-}, \text{ in } \Omega^{-}, $$

and

$$ \left\{ \begin{array}{ll} \Delta v^{+} + k^2 v^{+}, & \text{in } \Omega^{+}, \\ v^{+} \text{ outgoing wave.} & \end{array} \right. $$

Then, we have

$$ L([\partial_n v(\mathbf{y})]_{\Gamma})(\mathbf{x}) - D([v(\mathbf{y})]_{\Gamma})(\mathbf{x}) = \begin{cases} v^{-}(\mathbf{x}), & \text{for } \mathbf{x} \in \Omega^{-}, \\ v^{+}(\mathbf{x}), & \text{for } \mathbf{x} \in \Omega^{+}, \end{cases} \quad (19) $$

where

$$ [v]_{\Gamma} := \gamma_{0}^{-} v^{-} - \gamma_{0}^{+} v^{+}, \quad [\partial_{n} v]_{\Gamma} := \gamma_{1}^{-} v^{-} - \gamma_{1}^{+} v^{+}, $$
and

$$Lp(\boldsymbol{x}) := \int_{\Gamma} G(\boldsymbol{x}, \boldsymbol{y})p(\boldsymbol{y})d\Gamma(\boldsymbol{y}), \quad \boldsymbol{x} \notin \Gamma, \qquad (20)$$

$$D\phi(\boldsymbol{x}) := \int_{\Gamma} \partial_{n(\boldsymbol{y})} G(\boldsymbol{x}, \boldsymbol{y})\phi(\boldsymbol{y})d\Gamma(\boldsymbol{y}), \quad \boldsymbol{x} \notin \Gamma, \qquad (21)$$

$$\text{for } (p, \phi) \in H^{1/2}(\Gamma) \times H^{-1/2}(\Gamma).$$

The operators $L$ and $D$ defined above are called the single- and the double-layer potentials respectively. To obtain integral equations set on the boundary $\Gamma$, we need the trace formulae for these two potentials (see for instance [55]).

**Proposition 5.** The first and second traces on $\Gamma$ of the single-layer and the double-layer potentials are given by

$$\left\{ \begin{array}{l} \gamma_0^- \circ L = \gamma_0^+ \circ L = \mathcal{L} \\ \gamma_1^- \circ L = \left(\frac{\mathcal{I}}{2} + \mathcal{N}\right) \\ \gamma_1^+ \circ L = \left(-\frac{\mathcal{I}}{2} + \mathcal{N}\right) \end{array} \right. \qquad (22)$$

$$\left\{ \begin{array}{l} \gamma_0^- \circ D = \left(-\frac{\mathcal{I}}{2} + \mathcal{D}\right) \\ \gamma_0^+ \circ D = \left(\frac{\mathcal{I}}{2} + \mathcal{D}\right) \\ \gamma_1^- \circ D = \gamma_1^+ \circ D = \mathcal{S} \end{array} \right. \qquad (23)$$

where $\mathcal{I}$ is the identity operator and $\mathcal{L}, \mathcal{N}, \mathcal{D}$ and $\mathcal{S}$ the four elementary boundary integral operators expressed, for all $\boldsymbol{x} \in \Gamma$, by

$$\begin{align}
\mathcal{L}p(\boldsymbol{x}) &:= \int_{\Gamma} G(\boldsymbol{x}, \boldsymbol{y})p(\boldsymbol{y})d\Gamma(\boldsymbol{y}) \nonumber \\
\mathcal{N}p(\boldsymbol{x}) &:= \int_{\Gamma} \partial_{n(\boldsymbol{x})} G(\boldsymbol{x}, \boldsymbol{y})p(\boldsymbol{y})d\Gamma(\boldsymbol{y}) \nonumber \\
\mathcal{D}\phi(\boldsymbol{x}) &:= \int_{\Gamma} \partial_{n(\boldsymbol{y})} G(\boldsymbol{x}, \boldsymbol{y})\phi(\boldsymbol{y})d\Gamma(\boldsymbol{y}) \tag{24} \\
\mathcal{S}\phi(\boldsymbol{x}) &:= \oint_{\Gamma} \frac{\partial^2 G}{\partial \boldsymbol{n}(\boldsymbol{x}) \partial \boldsymbol{n}(\boldsymbol{y})} (\boldsymbol{x}, \boldsymbol{y})\phi(\boldsymbol{y})d\Gamma(\boldsymbol{y}). \nonumber
\end{align}$$

Note that the expression defining $\mathcal{S}$ is not an integral (its singularity is not integrable) but a finite part expression associated with a hypersingular kernel. We preferred to keep formally the integral expression for the sake of clarity. Let us now summarize the continuity properties of the elementary boundary integral operators (see for instance [55, Theorem 4.4.1] or Theorems 7.1 and 7.2 in [51]).

**Proposition 6.** For a smooth boundary $\Gamma$, the boundary integral operators given in Proposition 5 define the following continuous mappings

$$\begin{align*}
\mathcal{L} : H^s(\Gamma) &\longrightarrow H^{s+1}(\Gamma), \\
\mathcal{N} : H^s(\Gamma) &\longrightarrow H^s(\Gamma), \\
\mathcal{D} : H^s(\Gamma) &\longrightarrow H^s(\Gamma), \\
\mathcal{S} : H^s(\Gamma) &\longrightarrow H^{s-1}(\Gamma),
\end{align*} \qquad (25)$$
for all $s \in \mathbb{R}$. Moreover, the operators $\mathcal{N}$ and $\mathcal{D}$ are compact from $H^s(\Gamma)$ onto itself for all $s \in \mathbb{R}$.

In the case of a Lipschitz boundary [34, 51], the above continuity properties still hold for $-1 \le s \le 0$ (respectively for $0 \le s \le 1$) for operators $\mathcal{L}$ and $\mathcal{N}$ (respectively $\mathcal{D}$ and $\mathcal{S}$), while the compactness properties of $\mathcal{N}$ and $\mathcal{D}$ fail. A possible approach to rigorously extend the following developments is to use e.g. some regularizing techniques [24].

The representation (19) allows to determine the near-field around the scatterer. Another physical quantity of interest is the scattering amplitude (or the far-field pattern). For instance, in the two-dimensional case, we have

$$a_0(\boldsymbol{\theta}) = \frac{i}{4} \sqrt{\frac{2}{i\pi k}} \int_{\Gamma} e^{-ik\boldsymbol{y}\cdot\boldsymbol{\theta}} (\gamma_1^+ u^+(\boldsymbol{y}) + ik\boldsymbol{\theta} \cdot \boldsymbol{n}(\boldsymbol{y})\gamma_0^+ u^+(\boldsymbol{y})) d\Gamma(\boldsymbol{y}).$$

## 4.3 Standard integral equations formulations

The Helmholtz representation formula (19) leads to the construction of an infinite number of integral equations (equivalent if invertible) in the case of a closed surface. In the case of an open surface, only one integral equation can be written. The aim of this part is to introduce the most standard integral equations for both Dirichlet and Neumann boundary conditions. We usually distinguish between direct and indirect integral equations, each of them having their own mathematical properties.

Let us introduce the following notations

* $K_D(\Omega^-) = \{k_m^D, m \in \mathbb{N}\}$, the set of Dirichlet irregular frequencies (interior Dirichlet eigenvalues), is the set of values of $k$ such that the boundary value problem

$$\begin{cases} -\Delta v = k^2 v, & \text{in } \Omega^-, \\ \gamma_0^- v = 0, & \text{on } \Gamma, \end{cases}$$

admits a non vanishing solution.

* $K_N(\Omega^-) = \{k_m^N, m \in \mathbb{N}\}$, the set of Neumann irregular frequencies (interior Neumann eigenvalues), is the set of values of $k$ such that the boundary value problem

$$\begin{cases} -\Delta v = k^2 v, & \text{in } \Omega^-, \\ \gamma_1^- v = 0, & \text{on } \Gamma, \end{cases}$$

admits a non vanishing solution.

### 4.3.1 The Dirichlet problem: direct integral formulations

The total field $w$ is expressed by $w := u^+ + u^{\text{inc}}$. The direct formulations consist in seeking the total field under the form

$$w(x) = Lp(x) + u^{\text{inc}}(x), \quad x \in \Omega^{+}. \qquad (26)$$
The integral representation (26) ensures that $w$ is solution to the Helmholtz equation in $\Omega^- \cup \Omega^+$, and satisfies the Sommerfeld radiation condition. Then, we have to determine the unknown $p$ such that $w$ satisfies also the Dirichlet boundary condition $(\gamma_0^+ w = 0)$. More precisely, the representation (26) corresponds to the particular choice of solutions $(v^-, v^+) := (-u^{\text{inc}}, u^+)$ in proposition 4, i.e.

$$[v]_{\Gamma} = 0 \qquad [\partial_n v]_{\Gamma} = -\gamma_1^+ w|_{\Gamma} := p$$

and

$$Lp(x) = \begin{cases} -u^{\text{inc}}(x) & \text{for } x \in \Omega^{-} \\ u^{+}(x) & \text{for } x \in \Omega^{+}. \end{cases}$$

Then, we get the following single-layer potential representation of the total field:

$$Lp(x) + u^{\text{inc}}(x) = \begin{cases} 0 & \text{for } x \in \Omega^{-} \\ w(x) & \text{for } x \in \Omega^{+}. \end{cases}$$

Clearly, this formulation is completely equivalent to extending artificially the total field by zero inside $\Omega^-$, which explains that this approach is also referred sometimes to as the *null field method*.

The next step is to obtain an integral equation for the physical unknown $p = -\gamma_1^+ w|_{\Gamma} \in H^{-1/2}(\Gamma)$. To achieve this, the idea is to apply a trace operator to the relation

$$Lp(x) + u^{\text{inc}}(x) = 0, \quad \forall x \in \Omega^{-}. \qquad (27)$$

At this point, many choices are available. Let us cite three of them leading to classical integral equations of potential theory.

* **EFIE**: This equation is obtained by applying the trace operator $\gamma_0^-$ to (27). Thanks to the trace relations of Proposition 5, this leads to the well-known Electric Field Integral Equation (EFIE):

$$\mathcal{L}p = -\gamma_0^+ u^{\text{inc}}, \quad \text{on } \Gamma. \qquad (28)$$

* **MFIE**: This equation is obtained by applying the normal trace operator $\gamma_1^-$ to (27). Thanks to the trace relations of Proposition 5, this leads to the so-called Magnetic Field Integral Equation (MFIE):

$$(\frac{\mathcal{I}}{2} + \mathcal{N})p = -\gamma_1^+ u^{\text{inc}}, \quad \text{on } \Gamma. \qquad (29)$$

* **CFIE**: This equation is obtained by applying to (27) the Fourier-Robin (impedance) trace operator $\gamma_1^- + \eta\gamma_0^-$, with $\eta \neq 0$. Once again, the trace relations of Proposition 5 give the Combined Field Integral Equation (CFIE):

$$\left\{ \left( \frac{\mathcal{I}}{2} + \mathcal{N} \right) + \eta \mathcal{L} \right\} p = -(\gamma_1^+ u^{\text{inc}} + \eta \gamma_0^+ u^{\text{inc}}), \quad \text{on } \Gamma. \qquad (30)$$
The existence and uniqueness results for the above direct integral equations ((28), (29) or (30)) are given in the following theorem.

**Theorem 4.** The following properties hold.

1. The operator $\mathcal{L}$ defines an isomorphism from $H^{-1/2}(\Gamma)$ onto $H^{1/2}(\Gamma)$ if and only if $k \notin K_D(\Omega^-)$. Under this condition, the EFIE (28) is uniquely solvable in $H^{-1/2}(\Gamma)$.

2. The operator

$$ (\frac{\mathcal{I}}{2} + \mathcal{N}) $$

defines an isomorphism from $H^{-1/2}(\Gamma)$ onto $H^{1/2}(\Gamma)$ if and only if $k \notin K_N(\Omega^-)$. Under this condition, the MFIE (29) is uniquely solvable in $H^{-1/2}(\Gamma)$.

3. The operator

$$ (\frac{\mathcal{I}}{2} + \mathcal{N}) + \eta \mathcal{L} $$

defines an isomorphism from $H^{-1/2}(\Gamma)$ onto $H^{1/2}(\Gamma)$ for all $k > 0$ provided $\mathfrak{S}(\eta) \neq 0$ (imaginary part of $\eta$). Under this condition, the CFIE (30) is uniquely solvable in $H^{-1/2}(\Gamma)$ for any frequency $k > 0$.

In the case where $k$ in an irregular frequency, the integral equations EFIE and MFIE have non zero kernels. Nevertheless, it can be shown that the spurious modes of the EFIE will not radiate in the exterior. Thus, the field is not corrupted outside the object: $\mathcal{L}p = 0$ on $\Gamma \implies Lp = 0$ in $\Omega^+$. Then, the EFIE provides accurate computations and often represents a reference solution. Unlike the EFIE, the spurious solutions of the MFIE do radiate in the exterior domain, leading hence to a wrong solution. Finally, by its construction itself, the CFIE is free of the internal-resonance problem. We consider in the sequel $\eta = -ik\alpha/(1 - \alpha)$, $\alpha \in ]0, 1[$,

$$ (1 - \alpha) \frac{i}{k} (\frac{\mathcal{I}}{2} + \mathcal{N}) + \alpha \mathcal{L} = -((1 - \alpha) \frac{i}{k} \gamma_1^+ u^{\text{inc}} + \alpha \gamma_0^+ u^{\text{inc}}), \quad \text{on } \Gamma. \quad (31) $$

A common choice of $\alpha$ for engineering computations is $\alpha = 0.2$ which gives an almost minimal condition number for the CFIE.

### 4.3.2 The Dirichlet problem: indirect integral formulations

The indirect formulations are based on the assumption that the solution can be expressed in terms of a source density function defined on the boundary. The unknowns are then generally non-physical quantities. The physical variables are solved afterwards in terms of these source densities. Here, we focus on the most commonly used indirect integral formulation independently proposed by Burton-Miller [25] and Brakhage-Werner [19]. The idea is to seek the exterior field as a superposition of the single- and double-layer potentials acting on a fictitious surface density $\psi$:

$$ u^{+}(\boldsymbol{x}) = (D + \eta L)\psi(\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \Omega^{+}, \quad (32) $$
where $\eta$ is a complex-valued coupling parameter to choose. The above expression leads, thanks to the trace relations (5), to the following integral equation :

$$ \left\{ \left( \frac{\mathcal{I}}{2} + \mathcal{D} \right) + \eta \mathcal{L} \right\} \psi = -\gamma_0^+ u^{\text{inc}}, \quad \text{on } \Gamma. \qquad (33) $$

We consider the above integral equation in the space $H^{1/2}(\Gamma)$ and we can prove the following result.

**Theorem 5.** The operator

$$ \left( \frac{\mathcal{I}}{2} + \mathcal{D} \right) + \eta \mathcal{L} $$

defines an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{1/2}(\Gamma)$ for all $k > 0$ provided $\Im(\eta) \neq 0$. Under this condition, (33) is uniquely solvable in $H^{1/2}(\Gamma)$ for all frequency $k > 0$.

This integral equation is uniquely solvable if and only if $\Im(\eta) > 0$. An almost optimal value of $\eta$ has been obtained in [3, 46, 48] as: $\eta = ik$. We will see in Section 8 that other more subtle choices can be made in view of an iterative Krylov solution.

### 4.3.3 The Neumann problem: direct integral formulations

Let us now briefly discuss the derivation of direct integral equations in the case of a Neumann boundary condition. The total field $w := u^+ + u^{\text{inc}}$ is sought under the form

$$ w(x) = D\phi(x) + u^{\text{inc}}(x), \quad x \in \Omega^{+}. \qquad (34) $$

Proposition 4 for $(v^{-}, v^{+}) := (-u^{\text{inc}}, u^{+})$ leads to:

$$ [v]_{\Gamma} = -\gamma_{0}^{+} w := \phi, $$

$$ [\partial_{\mathbf{n}} v]_{\Gamma} = 0, $$

and

$$ D\phi(x) = \begin{cases} -u^{\text{inc}}(x) & \text{for } x \in \Omega^{-} \\ u^{+}(x) & \text{for } x \in \Omega^{+}. \end{cases} $$

Then, we get

$$ D\phi(x) + u^{\text{inc}}(x) = \begin{cases} 0 & \text{for } x \in \Omega^{-} \\ w(x) & \text{for } x \in \Omega^{+}. \end{cases} $$

Applying a trace operator to the relation

$$ D\phi(x) + u^{\text{inc}}(x) = 0, \quad \forall x \in \Omega^{-}, \qquad (35) $$

the physical unknown $\phi = -\gamma_0^+ w \in H^{1/2}(\Gamma)$ is solution to the following direct integral equations
• **EFIE :**
$$S\phi = -\gamma_1^+ u^{\text{inc}}, \quad \text{on } \Gamma. \qquad (36)$$

• **MFIE :**
$$(-\frac{\mathcal{I}}{2} + \mathcal{D})\phi = -\gamma_0^+ u^{\text{inc}}, \quad \text{on } \Gamma. \qquad (37)$$

• **CFIE :**
$$\{(-\frac{\mathcal{I}}{2} + \mathcal{D}) + \eta S\} \phi = -(\eta\gamma_1^+ u_{|\Gamma}^{\text{inc}} + \gamma_0^+ u^{\text{inc}}), \quad \text{on } \Gamma. \qquad (38)$$

The existence and uniqueness results for the above integral equations are summarized in the next result.

**Theorem 6.** *The following properties hold*

1. The operator $S$ defines an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{-1/2}(\Gamma)$ if and only if $k \notin K_N(\Omega^-)$. Under this condition, the EFIE (36) is uniquely solvable in $H^{1/2}(\Gamma)$.

2. The operator
$$(-\frac{\mathcal{I}}{2} + \mathcal{D})$$
defines an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{1/2}(\Gamma)$ if and only if $k \notin K_D(\Omega^-)$. Under this condition, the MFIE (37) is uniquely solvable in $H^{1/2}(\Gamma)$.

3. The operator
$$(-\frac{\mathcal{I}}{2} + \mathcal{D}) + \eta S$$
defines an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{-1/2}(\Gamma)$ for all $k > 0$ provided $\mathfrak{S}(\eta) \neq 0$. Under this condition, the CFIE (38) is uniquely solvable in $H^{1/2}(\Gamma)$ for all frequency $k > 0$.

The reference CFIE in this chapter is

$$(1 - \alpha) \frac{i}{k} \left(-\frac{\mathcal{I}}{2} + \mathcal{D}\right) - \frac{\alpha}{k^2} S = -\left((1 - \alpha) \frac{i}{k} \gamma_0^+ u^{\text{inc}} - \frac{\alpha}{k^2} \gamma_1^+ u^{\text{inc}}\right), \quad \text{on } \Gamma. \qquad (39)$$

### 4.3.4 The Neumann problem: indirect integral formulations

The Burton-Miller (or Brakhage-Werner) integral representation of the exterior field is expressed by

$$u^{+}(x) = (L + \eta D)\varphi(x), \quad \forall x \in \Omega^{+}, \qquad (40)$$

where $\eta$ is a complex-valued coupling parameter to determine. Then, the field (40) solves the exterior boundary-value problem (17) if the surface density $\varphi$ is solution to the following integral equation

$$\left\{ \left(-\frac{\mathcal{I}}{2} + \mathcal{N}\right) + \eta S \right\} \varphi = -\gamma_1^+ u^{\text{inc}}, \qquad (41)$$
called the Burton-Miller or Brakhage-Werner (BW) integral equation. We have the following existence and uniqueness result.

**Theorem 7.** *The operator*

$$(-\frac{\mathcal{I}}{2} + \mathcal{N}) + \eta\mathcal{S}$$

defines an isomorphism from $H^{1/2}(\Gamma)$ onto $H^{-1/2}(\Gamma)$ for all $k > 0$ provided $\mathfrak{S}(\eta) \neq 0$. Under this condition, (41) is uniquely solvable in $H^{1/2}(\Gamma)$ for all frequency $k > 0$.

An almost optimal value of $\eta$ has been numerically discussed in [3, 46, 48] as: $\eta = 1/ik$. Again, we will see in Section 8 that better choices can be considered.

### 4.3.5 First- vs second-kind fredholm integral equation

For surface formulations of acoustic scattering problems, we have seen in the previous sections that various integral equations, direct or indirect, combined or not, can be employed. Let us underline that all these integral equations are applicable to closed geometries, leaving EFIE ((28) and (36) in the Dirichlet and Neumann cases respectively) as the only choice for the solutions of scattering problems by open surfaces.

In view of an iterative solution of an integral equation, two properties are essential to achieve a fast convergence rate: on the one hand the existence and uniqueness of the solution, on the other hand the clustering of the eigenvalues of the underlying integral operator (and hence a well-conditioning). The first property, as we have seen, is not ensured for each integral equation. Only combined field integral equations (CFIE, Brakhage-Werner) provide a well-possedness for all frequencies. To observe the second property, let us recall some definitions. Given an integral operator $\mathcal{A} \in \mathcal{L}(X)$ on a Hilbert space $X$, an integral equation is called of first-kind if it is of the form

$$\mathcal{A}\varrho = f$$

of second-kind if it is of the form

$$(\mathcal{I} + \mathcal{A})\varrho = f. \qquad (42)$$

Moreover, if $\mathcal{A}: X \to X$ is compact, the above equations are respectively called Fredholm integral equations of the first-kind and second-kind. As classically known [47], the spectrum of compact operators is composed in the infinite dimensional case of 0 and a sequence of discrete eigenvalues possibly accumulating at the origin. Therefore, second-kind Fredholm integral equation have large clusters of eigenvalues accumulating at the real value point 1 (cf. Section 2.3). This is a very interesting spectral property in view of an iterative solution using a Krylov solver. Indeed, eigenvalue clustering in the complex plane generally implies fast convergence.

The integral operator $\mathcal{L}$ is bounded from $H^s(\Gamma)$ onto $H^{s+1}(\Gamma)$, and it is compact from $L^2(\Gamma)$ onto itself. Therefore, the EFIE (28) is a Fredholm first-kind integral equation on $L^2(\Gamma)$ in the case of a Dirichlet condition. On the contrary, the MFIE
(29), the CFIE (30) and the Burton-Miller integral equation (33) are all second-kind integral equations. Then, the combined field integral equations have all the properties required for an efficient iterative solving.

In the case of a Neumann boundary condition, the situation is more complex. In fact, only the MFIE is a second-kind Fredholm integral equation but it is unfortunately ill-posed. The EFIE, CFIE and Burton-Miller integral equations involve the operator $S$ which is a first-order, strongly singular and non-compact operator. Therefore, these equations are first-kind integral equations. There are two alternatives to expect an eigenvalue clustering. The first one is to precondition the integral operator (see Section 7). The second possibility is to incorporate a suitable operator of order $-1$ which has a regularizing effect on $S$ and leads to a second-kind integral equation (see Section 8). Let us note that in the case of the numerical solution of scattering problems by an open surface, the preconditioning represents the only way to improve convergence. We come back to this point in section 7.

# 5 Developing fast converging solutions: strategy and problematics

Let us consider now that we have chosen an integral equation representation

$$A\varrho = f$$

that we want to numerically solve. Here, $A$ is one of the previous first- or second-kind integral operators, $f$ is a right-hand side given by the incident field and $\varrho$ is the density that we would like to compute. For a numerical calculation, the surface of the scatterer needs to be discretized by using for example $N_K$ triangles $K$ in three-dimensions, resulting in a discrete surface $\Gamma_h = \cup_{j=1}^{N_K} K$, where $h$ is the meshsize (see Section 8.3 for more details). Let us consider a regular triangulation $\mathcal{T}_h$ based on triangles. The linear Galerkin boundary element method is based on the approximation space

$$V_h = \{v_h \in C^0(\Gamma_h); v_h|_{\Gamma} \in \mathbb{P}_1, \forall T \in \mathcal{T}_h\},$$

of dimension $N_V$ (equal to the number of degrees of freedom). The density of discretization points per wavelength $n_\lambda$ is given by $n_\lambda = \lambda/h_{max}$, where $h_{max}$ is the maximal length of the edges of the triangles. Then, the operator $A$ is also approximated by a matrix $[A]$ and the right hand-side $f$ by a complex-valued vector $f_h$. The sizes of the vector and the matrix $[A]$ are equal to the number of degrees of freedom of the boundary element method used to approximate the density $\varrho$ by a vector $\varrho_h$. For example, using linear boundary elements leads to $N_V$ degrees of freedom, where $N_V$ is the number of vertices of the triangular mesh. One problem that we do not address here but which is a hard task in integral equations consists in integrating the kernel singularities. However, since it is out of our goal, we do not develop this point here.

Consider now that we want to solve the linear system

$$[A]\varrho_h = f_h.$$
The matrix [A] is complex-valued, dense and highly non-definite positive. It requires a storage of the order of $O(N_V^2)$ (if a linear boundary element method is used) and its solution by a direct gaussian solver requires $O(N_V^3)$ operations. For high frequency problems, then $N_V$ becomes extremely large and makes the direct approach not applicable. Therefore, since the introduction, about twenty years, of the Fast Multilevel Multipole (FMM) method [32, 36, 37, 40, 62], the strategy for solving an integral equation has fundamentally changed. The FMM method (as well as other similar techniques, see [21, 22, 23]) allows to compute with a low-storage $O(N_V)$ and in a fast way $O(N_V \log N_V)$ the Matrix-Vector Product (MVP): $x_h \to [A]x_h$ ("black box"). This matrix-free approach allowed to solve higher frequency problems if it is coupled with an iterative Krylov solver like the GMRES. As previously noticed, the global algorithm is efficient and robust if the GMRES converges with a few iterations. This is closely related to the spectral properties of A as well as the way a preconditioner is built. Since we have a non positive-definite matrix, the iterative algorithm may diverge. If not, its convergence can be extremely slow. Furthermore, since we do not have the matrix at hand (because of the matrix-free FMM approach), building the preconditioner must be done without having access to the whole matrix.

The aim of the next sections are 1) to propose a spectral analysis of integral operators in some simplified situations to understand their properties (Section 6) and 2) to explain two recent matrix-free analytical preconditioning techniques (Sections 7 and 8) for solving integral equation formulations.

# 6 Spectral analysis of integral operators for the sphere

Let $\Gamma = S_1$ be the unit sphere centered at the origin. Let us introduce the spherical harmonics $Y_m^n$ as the functions of order $m$ for $n = -m \cdots m$, with $m \in \mathbb{N}$ [33], given by

$$Y_m^n(\theta_1, \theta_2) = \sqrt{\frac{2m+1}{4\pi} \frac{(m-|n|)!}{(m+|n|)!}} P_m^n(\cos\theta_1) e^{in\theta_2}, \quad (\theta_1, \theta_2) \in S_1,$$

where $P_m^n$ are the Legendre polynomials. The functions $Y_m^n$ form a complete orthonormal system of $L^2(S_1)$. Furthermore, they also constitute a basis of eigenvectors for the four elementary integral operators $\mathcal{L}$, $\mathcal{N}$, $\mathcal{D}$ and $\mathcal{S}$ (cf. proposition 5). More precisely, we have the following proposition [4, 46].

**Proposition 7.** The eigenvalues $\mathcal{L}_m$, $\mathcal{N}_m$, $\mathcal{D}_m$ and $\mathcal{S}_m$ of multiplicity $(2m + 1)$ respectively associated with the elementary integral operators $\mathcal{L}$, $\mathcal{N}$, $\mathcal{D}$ and $\mathcal{S}$ are given by

$$\mathcal{L}Y_m^n = \{ikj_m(k)h_m^{(1)}(k)\}Y_m^n = \mathcal{L}_m Y_m^n,$$

$$\mathcal{N}Y_m^n = \{-\frac{1}{2} + ik^2 j_m'(k)h_m^{(1)}(k)\}Y_m^n = \mathcal{N}_m Y_m^n,$$

$$\mathcal{D}Y_m^n = \mathcal{N}^t Y_m^n = \mathcal{N}_m Y_m^n = \mathcal{D}_m Y_m^n,$$

$$\mathcal{S}Y_m^n = \{ik^3 j_m'(k)(h_m^{(1)})'(k)\}Y_m^n = \mathcal{S}_m Y_m^n,$$
where $j_m$ denotes the spherical Bessel function and $h_m^{(1)}$ the spherical Hankel function of the first-kind. Their derivatives are specified with a prime'.

Let us begin with the case of a sound-soft sphere. We consider the EFIE operator $\mathcal{L}$ (28), the CFIE operator (31) and the usual BW operator (33) with $\eta = ik$. A direct computation gives the eigenvalues $C_m^D$ and $B_m^D$ of the combined operators (31) and (33) respectively

$$C_m^D = ik^2 h_m^{(1)}(k)(\alpha j_m'(k) + i(1-\alpha)j_m(k)), \\ B_m^D = ik^2 h_m^{(1)}(k)(j_m'(k) - ij_m(k)).$$

We draw on Fig. 5 the spectrum of the usual integral operators taking $k = 30$ and a maximal number of modes $m_{max} = 12k$. To understand the spectral distribution, we need to introduce three zones: the elliptic zone of evanescent modes ($|m| \gg k$), the hyperbolic region of propagative modes ($|m| \ll k$) and the transition (hyperbolic-elliptic) zone of physical surface modes satisfying $|m| \approx k$. To make the connection with pseudodifferential operators, think for easiness that $\xi = m$. For the EFIE (Figure 5(a)), we can see that some eigenvalues spread out the complex plane and are related to the finite number of propagative modes. Next, a loop begins to appear for modes close to the transition zone. Finally, we observe an accumulation point at $(0, 0)$ for large modes $|m|$. Since $\mathcal{L}$ is a first-kind operator of order $-1$, its eigenvalues asymptotically behave like $|m|^{-1}$ for large values of $|m|$ (elliptic zone). This implies that the smallest eigenvalue is related to the largest mode $m_{max}$ and makes appear a dependence of the condition number according to the density of discretization points per wavelength $n_\lambda$ for integral equations [29]. The largest eigenvalue is related to a propagative mode and leads to a dependence of the condition number according to $k$. Therefore, the EFIE has a spectrum which is not really appropriate to an iterative solver and preconditioning will have to be considered, most particularly for large wavenumbers $k$ and large densities of discretization points. Two large clusters of eigenvalues at points $(1/2, 0)$ and $(1, 0)$ can be observed for the BW operator (Figure 5(b)). The first appears in the elliptic zone and the second one in the hyperbolic zone. The few eigenvalues that form the loop between low-order and high-order modes $|m|$ correspond to surface modes. The CFIE operator also offers an interesting and quite similar distribution of eigenvalues (Figure 5(c)). These good spectral properties are linked to the Fredholm second-kind character of these operators. As a consequence, the standard CFIE and BW formulations are well-adapted for the iterative solution of the acoustic sound-soft problem.

Let us now consider the hard sphere. The eigenvalues $C_m^N$ and $B_m^N$ of the CFIE (39) and the BW ((33), $\eta = i/k$) operators are expressed respectively by

$$C_m^N = k(h_m^{(1)})'(k)((1-\alpha)j_m(k) + i\alpha j_m'(k)), \\ B_m^N = -ik^2(h_m^{(1)})'(k)(j_m(k) + ij_m'(k)).$$

We draw on Fig. 5 the distribution of the eigenvalues for the EFIE, BW and CFIE o-perators. We fix $k = 30$ and $m_{max} = 12k$ modes. The EFIE operator is of order 1 and its eigenvalues associated with high order modes become large since they
<figure><img src="image_19.png" /><figcaption>(a) Sound-soft: EFIE</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Sound-soft: Brakhage-Werner (BW)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) Sound-soft: CFIE</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) Sound-hard: EFIE</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(e) Sound-hard: Brakhage-Werner (BW)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(f) Sound-hard: CFIE</figcaption></figure>

Figure 5: Sound-soft and sound-hard spheres: distribution of the eigenvalues.
behave like $|m|$ (Figure 5(d)). The smallest eigenvalue is associated with the prop- agative modes. For the usual BW integral operator (Figure 5(e)), we observe a cluster of eigenvalues linked to the low-order modes. A large number oonf eigenvalues cor- responding to the evanescent modes are on the line $x = 1/2$ and do not cluster like for the sound-soft case. This behaviour penalizes the convergence rate of the GM- RES. Similar conclusions arise for the CFIE (Figure 5(f)). This is due to the fact that these three integral equations are Fredholm first-kind integral equations and so are ill-conditionned. In particular, the convergence rate of a Krylov iterative solver will depend on the wavenumber $k$ as well as the mesh density. Using some asymptotics of the special functions appearing above and some formal replacements, some similar conclusions can be drawn for these integral equations and a general convex surface $\Gamma$. Finally, other geometrical configurations lead to spectral estimates of integral opera- tors (see [28, 29]).

From this analysis, we clearly see that some efforts on analytical (or algebraic) preconditioning must be directed towards the EFIE for the sound-soft problem and EFIE, BW and CFIE for the sound-hard case.

# 7 A first direction to improve the convergence: algebraic/analytical preconditioning techniques for the EFIE

As seen above, the EFIE has a condition number which depends on both the wavenum- ber $k$ and the density of discretization points per wavelength $n_\lambda$ for both the sound- soft and sound-hard scattering problems. Furthermore, the eigenvalues spread out the complex plane which is very penalizing in view of an iterative Krylov solution. One way to improve the convergence properties of the EFIE is to precondition it. Algebraic preconditioners have been proposed over the years but may fail, most particularly for large wavenumbers. Furthermore, they need to handle the full matrix which is in- compatible with the idea of a matrix-free solver imposed by the FMM method. We propose here two possible directions for building analytical preconditioners: the first one uses the Calderón integral relations and the second one elements from the pseu- dodifferential operator theory.

## 7.1 Integral operator preconditioning

The idea of integral operator preconditioning uses the following Calderón relations [55]

$$ -\mathcal{L}S = \frac{\mathcal{I}}{4} - \mathcal{N}^2, \quad -\mathcal{S}\mathcal{L} = \frac{\mathcal{I}}{4} - (\mathcal{N}^T)^2. \qquad (43) $$

From these two equations, we can see that, since $\mathcal{N}$ is a compact operator, then $\mathcal{N}^2$ and $(\mathcal{N}^T)^2$ are compact and their eigenvalues tend towards zero. This also means that
the eigenvalues of $-\mathcal{L}S$ and $-\mathcal{SL}$ cluster around $(1/4, 0)$ in the complex plane since they are equal to the operator $\mathcal{I}/4$ up to a compact perturbation. As a consequence, $-\mathcal{L}$ (respectively $-\mathcal{S}$) is a pseudo inverse operator of $S$ (respectively $-\mathcal{L}$) and can serve as a preconditioner. This is an interesting property since the application of this preconditioner in a FMM environment only involves some evaluations of standard integral operators. This first application of Calderón relations to preconditioning has been developed by Steinbach and Wendland in [61]. It has been next applied to the EFIE for acoustics and electromagnetism by Christiansen and Nédélec in [30]. Extensions and other studies related to this approach are available in [9, 10, 11].

To show the improvement induced by a preconditioner based on Calderón relations, we consider the transmission scattering problem which consists in solving an exterior Helmholtz problem with wavenumber $k_2$ coupled with an interior Helmholtz problem for a wavenumber $k_2\sqrt{N}$ with transmission boundary conditions. This physically corresponds to scattering by a penetrable homogeneous isotropic scatterer. Following [11], we propose a Calderón preconditioner for the two-field integral equation solution to the transmission problem. We refer to [11] for more details. The Calderón preconditioner for the two-field integral equation is denoted by $C$. We consider on Figure 6(a) the convergence history of GMRES(50) for the preconditioned linear system in the case of the unit square cylinder (centered at the origin and with sidelength 2). We represent the dependence of the preconditioned algorithm according to both the wavenumber $k_2$ and the density $n_\lambda$ of discretization points per wavelength. The index $\sqrt{N}$ for the interior problem is $\sqrt{N} = 2 + i$. We can observe that the convergence is independent of the density $n_\lambda$ since it can be shown that the preconditioned integral equation is Fredholm second-kind. However, the convergence depends moderately on $k_2$, which is one of the limitations for Calderón preconditioners. Another example is given by the scattering problem by a penetrable kite-shaped object (see Figure 6(b)) for $\sqrt{N} = 1.55 + 0.64i$. The same conclusions can be made from Figure 6(c). We can see on Figure 6(d) (for $k_2 = 40$ and $n_\lambda = 10$) that even if some eigenvalues spread out the complex plane and imply a $k_2$-dependence of the iterative solver, a large cluster of eigenvalues characterizes the preconditioned matrix and results in the $n_\lambda$ independence. This first example shows the impact that the point of view of analytical preconditioning can have on solving a scattering problem. Let us finally remark that examples in [11] show that ILU preconditioners fail in general leading to a breakdown of the GMRES.

## 7.2 Pseudodifferential operator preconditioning

A second approach for preconditioning the EFIE uses the following result which precises the principal symbol of the single-layer $\mathcal{L}$ and normal derivative trace $\mathcal{S}$ of the double-layer potentials.

**Proposition 8.** Let $\mathcal{L}$ and $\mathcal{S}$ be respectively the single-layer and normal derivative trace of the double-layer potentials defined by the expressions (24). Let $\xi$ be the dual variable of $x$ by Fourier transform for $x$ restricted to $\Gamma$. Then, the principal symbols
<figure><img src="image_25.png" /><figcaption>(a) Convergence history of GMRES(50) for the unit square cylinder.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) The kite-shaped scatterer.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c) Convergence history of GMRES(50) for the kite-shaped scatterer.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(d) Eigenvalue distribution.</figcaption></figure>

<caption>Figure 6: Calderón preconditioner for transmission problems.</caption>
of $\mathcal{L}$ and $S$, denoted by $\sigma_p(\mathcal{L})$ and $\sigma_p(S)$, are given by

$$ \sigma_p(\mathcal{L}) = \frac{i}{2\sqrt{k^2 - |\xi|^2}} \quad \text{and} \quad \sigma_p(S) = -\frac{\sqrt{k^2 - |\xi|^2}}{2i}. \qquad (44) $$

Following an approach similar to the one developed in Section 3.3 but for integral equations, an appropriate preconditioner for the sound-soft scattering problem based on the EFIE representation is $Op(\sigma_p(\mathcal{L})^{-1}) = -4Op(\sigma_p(S))$ since then

$$ -Op(\sigma_p(S))\mathcal{L} = \frac{\mathcal{I}}{4} + R, \qquad (45) $$

where $R$ is a pseudodifferential operator of order less or equal to $-1$. One of the crucial points in this approach is that $k$ is also considered as a symbol (associated with the time derivative operator $\partial_t$ if we come back to a wave equation) of order 1. Then, the symbol $\sigma_p(S)$ is homogeneous of order 1 and defines a corresponding classical pseudodifferential operator of order 1. We see here that the equation (45) is typically a pseudodifferential operator version of the second Calderón relation (43). The main difference is that the operator $Op(\sigma_p(S))$ is a correct approximation of the inverse of $\mathcal{L}$ in both its hyperbolic and elliptic zones. This is not the case of $S$ which is only a good approximation of the inverse in the elliptic part since the error term $-(\mathcal{N}^T)^2$ is compact (for large values of $|\xi|$ compared to $k$). Consequently, the convergence rate of the iterative Krylov solver is independent of the mesh refinement parameter $n_\lambda$ and the dependence according to $k$ is weakened. For the sound-hard scattering problem, a suitable choice is $-4Op(\sigma_p(\mathcal{L}))$ for preconditioning $S$ and the same conclusions arise.

Concerning the implementation in a FMM code, including the operator $Op(\sigma_p(S))$ can be made independently of the integral equations. For effectiveness, a complex Padé approximation of the square-root must be considered (see Section 8.2), resulting in the numerical solution of a coupled system of surface dissipative Helmholtz equations on $\Gamma$. Their approximate solution can be considered by using an ILU preconditioner and an iterative scheme at the cost $O(N_V)$, which is much less than applying an integral operator. We refer to [10] for the interested reader where much more details about the implementation are given. Some partial elements are also developed in Section 8.3.

To show the robustness of such an approach, we report on Figures 7 the number of MVPs required for solving the sound-hard scattering problem by a strip of length 2 with GMRES(50) for a tolerance $\epsilon = 10^{-6}$. We use the EFIE and represent the number of MVPs with respect to the density $n_\lambda$ (for $k = 15$, figure 7(a)) and the wavenumber $k$ (for $n_\lambda = 20$, figure 7(b)). The EFIE is referred to as "DBIE" here, the pseudodifferential preconditioner based on $-4Op(\sigma_p(\mathcal{L}))$ and a Padé approximation by "Padé-type", the Calderón integral preconditioner $[\mathcal{L}]$ by "Calderon". Furthermore, an algebraic SPAI (SParse Approximate Inverse) preconditioner with sparsity level $N_v = 5$ is also used for comparison. The choice $N_v = 5$ makes the preconditioning matrix already quite dense (about 20% of nonzero coefficients). Another drawback is that it is built on the full initial EFIE matrix $[S]$ which is not available when considering a FMM solution. On this simple example, we can see that the number of iterations
<figure><img src="image_29.png" /><figcaption>(a) MVPs vs. $n_\lambda$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) MVPs vs. $k$.</figcaption></figure>

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

<figure><img src="image_32.png" /><figcaption>(c) Convergence history of GMRES(50) for the (d) Convergence history of GMRES(50) for the cobra-like scatterer.</figcaption></figure>

<caption>Figure 7: Pseudodifferential preconditioner for scattering problems.</caption>

is independent of both $n_\lambda$ and $k$. We report on Figure 7(c) the convergence history for all methods. We can see the superiority here of the two analytical preconditioners over the algebraic SPAI. A second example is given on Figure 7(d) for a cobra-like shaped scatterer which is an open curved resonator modeling an inlet. We can again observe the very good behaviour of the iterative solvers compared with analytical preconditioners.

In references [10], other numerical examples show that some problems can arise for the Calderón preconditioner if the scatterer is composed of two parallel strips or screens due to resonance phenomena. It seems that, for this kind of situations, the pseudodifferential preconditioners are more robust. However, for a single three-dimensional open surface, some problems are still open like for instance how to correctly handle the edges effects in a pseudodifferential approach. In the case of Calderón operators, this point is quite naturally treated. A hybrid Calderón-pseudodifferential preconditioner is proposed in [35].
# 8 A second direction to improve the convergence: generalized Combined Field Integral Equations

We restrict our presentation to a Neumann boundary condition. Indeed, we have seen in Section 6 that the standard combined field integral equations (both BW and CFIE) are second-kind Fredholm integral equations for a sound-soft obstacle. As a consequence, the convergence rate of an iterative solver applied to these equations is independent of the density $n_λ$ but still depends slightly on $k$. In the same spirit as what is proposed below, this dependence can be weakened by developing generalized (direct and indirect) combined field integral equations. In the case of the sound-hard obstacle, much work is required since the standard integral equations are first-kind Fredholm integral equations. Furthermore, our construction of new second-kind integral equations preserves the well-posedness of these new integral equations, meaning that they are all free of any spurious resonance.

## 8.1 Construction of generalized combined field integral equations

The starting point of the method proposed in [12, 13, 35] is based on the exterior Helmholtz integral formulation

$$u^+ = -L\gamma_1^+ u^+ + D\gamma_0^+ u^+, \text{ in } \Omega^+. \quad (46)$$

Applying the exterior normal derivative trace to (46), we get thanks to Proposition 5

$$\gamma_1^+ u^+ = \left(-\frac{\mathcal{I}}{2} + \mathcal{N}\right) \gamma_1^+ u^+ + \mathcal{S}\gamma_0^+ u^+, \text{ in } H^{1/2}(\Gamma). \quad (47)$$

Let us assume that the exact exterior Neumann-to-Dirichlet (NtD) map $V^{\text{ex}}$ is known

$$V^{\text{ex}}: \begin{array}{ccc} H^{-1/2}(\Gamma) & \to & H^{1/2}(\Gamma) \\ \gamma_1^+ u^+ & \mapsto & \gamma_0^+ u^+ = V^{\text{ex}}\gamma_1^+ u^+. \end{array} \quad (48)$$

Then, from (47) and (48), the following identity holds

$$-\frac{\mathcal{I}}{2} + \mathcal{N} + \mathcal{S}V^{\text{ex}} = \mathcal{I}. \quad (49)$$

The NtD operator $V^{\text{ex}}$ is a nonlocal pseudodifferential operator of order $-1$. The equation (49) means that plugging the NtD operator into the integral formulation leads to the identity operator which can be solved trivially in one iteration (direct solver). Unfortunately, the NtD operator cannot be computed explicitly for a general surface, otherwise, this would mean that an explicit integral equation can be written for a scattering problem and solved in one step.

An alternative is to consider an accurate approximation $\tilde{V}$ of the exact NtD operator (48) such that

$$\varphi = \tilde{V}\psi,$$
with $(\varphi, \psi)$ an approximation of the Cauchy data $(\gamma_0^+ u^+, \gamma_1^+ u^+)$. In the same spirit as for the exact NtD operator, we propose to solve the following integral equation

$$ \left(-\frac{\mathcal{I}}{2} + \mathcal{N} + S\tilde{V}\right)\psi = -\gamma_1^+ u^{\text{inc}}, \quad \text{in } H^{1/2}(\Gamma), \qquad (50) $$

based on the exterior integral representation

$$ u^+ = -L\psi + D\varphi, \quad \text{in } \Omega^+. \qquad (51) $$

Let us for example consider the lowest order approximation $\tilde{V} = -i/k$ (Sommerfeld radiation condition). This choice leads to the standard BW integral equation (41) with the optimal parameter of Kress [46]. Hence, the integral equation (50) can be seen as a generalization of the standard BW integral equation. A similar approach can also be developed for constructing a generalized CFIE

$$ \left( (1-\alpha)\frac{i}{k}\left(-\frac{\mathcal{I}}{2} + \mathcal{D}\right) - \frac{\alpha}{k^2}\tilde{V}\mathcal{S} \right)\phi = -\left( (1-\alpha)\frac{i}{k}\gamma_0^+u^+ + \tilde{V}\gamma_1^+u^{\text{inc}} \right), \quad \text{in } H^{1/2}(\Gamma), \quad (52) $$

by composition of the operator $\mathcal{S}$ (EFIE) with $-\tilde{V}$ and then adding the contribution of the MFIE by a suitable linear combination. One of the aims of the pseudodifferential operator $\tilde{V}$ is to regularize the operator $\mathcal{S}$ of order 1. It should therefore be of order $-1$ to get the identity operator. This is not possible by considering the lowest order approximation which is of order zero. For this reason, high-order approximations of the NtD operator must be considered. This can be done for example by adapting techniques related to the Beam Propagation Method [42] or On-Surface Radiation Conditions (OSRCs) methods [6, 7, 8, 14, 43, 49]. In particular, we refer to the recent review chapter book which gives an overview of OSRCs approaches. From [14], an accurate approximation of the NtD map is given by the square-root operator

$$ \tilde{V} = \frac{1}{ik}\left(1 + \frac{\Delta_{\Gamma}}{k_{\delta}^{2}}\right)^{-1/2}, \qquad (53) $$

where the operator $\Delta_{\Gamma}$ is the usual Laplace-Beltrami operator over the surface $\Gamma$ and the parameter $k_{\delta} = k + i\delta$ is a complex wavenumber. The aim of the damping parameter $\delta > 0$ is to regularize the square-root operator in the transition zone of grazing rays [14]. An optimal value of this parameter, related to the mean curvature of $\Gamma$ and to $k$ is given in [14]. Incorporating the operator $\tilde{V}$, the generalized combined integral equations (50) and (52) are uniquely solvable [35] for any frequency $k > 0$ and damping parameter $\delta > 0$. Moreover, they are second-kind Fredholm integral equations unlike the standard formulations (50) and (39). This is a positive point in view of an iterative Krylov solution.

## 8.2 The special case of the unit sphere

The eigenvalues of the operator $-\Delta_{\Gamma}$ are expressed by $\mu_m = m(m+1)$. A direct computation gives the eigenvalues $B_m^{N,\delta}$ and $C_m^{N,\delta}$ of respectively the generalized BW
<figure><img src="image_33.png" /><figcaption>Figure 8: Sound-hard sphere: distribution of the eigenvalues for the generalized combined field integral equations.</figcaption></figure>

and the CFIE operators as

$$ \mathcal{B}_m^{N,\delta} = -ik^2(h_m^{(1)})'(k)[ij'_m(k)(1 - \frac{\mu_m}{k_\delta^2})^{-1/2} + j_m(k)], $$

and

$$ \mathcal{C}_m^{N,\delta} = k^2(h_m^{(1)})'(k)[-i(1-\alpha)j_m(k) + \alpha(1-\frac{\mu_m}{k_\delta^2})^{-1/2}j'_m(k)]. $$

We draw on Fig. 8 the eigenvalues of the generalized integral operators by taking $k = 30$ and $m_{max} = 12k$ modes. The eigenvalues of the generalized BW and CFIE operators cluster around the points $(1, 0)$ and $(1/2, 0)$ respectively. This observation is partially related to the property that these integral equations are Fredholm second-kind but also to the fact that an accurate approximation of the NtD operator is injected. By using suitable asymptotic expansions of the Bessel and Hankel functions [1], we justify this remark by the following approximations of the eigenvalues [35]

* in the elliptic zone (evanescent modes) for large values of $m$

$$ \mathcal{B}_m^{N,\delta} = \frac{1}{2} + \frac{k_\delta}{2k} + \mathcal{O}(\frac{1}{m^3}), \quad \text{and} \quad \mathcal{C}_m^{N,\delta} = \frac{1}{2} + \alpha(\frac{k}{2k_\delta} - \frac{1}{2}) + \mathcal{O}(\frac{1}{m^3}), \quad \alpha \in [0, 1]. $$

* in the hyperbolic zone (propagative modes) for large wavenumbers $k$

$$ \mathcal{B}_m^{N,\delta} = 1 + \mathcal{O}(\frac{1}{k}), \quad \text{and} \quad \mathcal{C}_m^{N,\delta} = \frac{1}{2} + \mathcal{O}(\frac{1}{k}). $$

The coupling between low- (propagative) and high-order (evanescent) modes is observable as a loop around the accumulation points $(1, 0)$ and $(1/2, 0)$ for the generalized BW and CFIE respectively. From these observations, we can expect the fast convergence of a Krylov solver for computing the solution to (50) or (52).
## 8.3 Numerical issues: implementation and examples

Let us denote by $[P]$ the matrix associated with the linear discretization of a given integral operator $P$. The discretization of (50) leads to the $N_V \times N_V$ dense complex linear system

$$ \left(-\frac{[\mathcal{I}]}{2} + [\mathcal{N}] + [\mathcal{S}][\tilde{V}]\right)\psi_h = g_h. \qquad (54) $$

The complex valued vectors $\psi_h$ and $g_h$ of $\mathbb{C}^{N_V}$ are respectively the $\mathbb{P}_1$ interpolated unknown density $\psi$ and the right-hand side of (50). The matrix representation $[\mathcal{I}]$ of the identity operator $\mathcal{I}$ is the surface mass matrix. We refer for example to [5, 36] for the direct approximation of the integral operators $\mathcal{N}$ and $\mathcal{S}$. Of course, a FMM implementation can also be derived in the background of Galerkin methods for a fast evaluation of MVP. Concerning the discretization of the pseudodifferential operator $\tilde{V}$, a rational Padé approximant of the square-root operator is used. In [14], the authors use the rotating branch-cut technique introduced by Milinazzo al. [54]. It consists in using the rational approximation defined by

$$ \sqrt{1+X} \approx e^{i\theta_p/2} R_{N_p}((1+X)e^{-i\theta_p}-1) = C_0 + \sum_{j=1}^{N_p} \frac{A_j X}{1+B_j X}, $$

where the parameter $\theta_p$ is the rotation angle of the usual branch-cut $\{z \in \mathbb{C}, z < -1\}$ of the square-root operator, and $R_{N_p}$ denotes the usual Padé approximant of order $N_p$. The complex-valued coefficients $C_0, (A_j, B_j)_{j=1,...,N_p}$ are related to the real Padé coefficients and the angle $\theta_p$ [14]. Since the square-root operator is implemented in an iterative solver, we must be able to compute $\mathbf{y} = [\tilde{V}]\mathbf{x}$ efficiently, for a given vector $\mathbf{x} \in \mathbb{C}^{N_V}$. This can be achieved by writing down

$$ (1 + \frac{\Delta_{\Gamma}}{k_{\delta}^{2}})\mathbf{y} = \frac{1}{ik}(1 + \frac{\Delta_{\Gamma}}{k_{\delta}^{2}})^{1/2}\mathbf{x}. $$

Then, the MVP $\mathbf{y} = [\tilde{V}]\mathbf{x}$ is computed by solving $N_p$ dissipative surface Helmholtz equations

$$ \left(\frac{B_j}{k_{\delta}^2}[\Delta_{\Gamma}] + [\mathcal{I}]\right)\mathbf{x}_j = [\mathcal{I}]\mathbf{x}, \quad j = 1, \dots, N_p, \qquad (55) $$

using a variational formulation and boundary elements and then by computing $\mathbf{y}$ solution to

$$ ([\mathcal{I}] + \frac{[\Delta_{\Gamma}]}{k_{\delta}^2})\mathbf{y} = \frac{1}{ik}(C_0[\mathcal{I}]\mathbf{x} + \sum_{j=1}^{N_p} \frac{A_j}{k_{\delta}^2}[\Delta_{\Gamma}]\mathbf{x}_j). \qquad (56) $$

The matrix $[\Delta_{\Gamma}]$ is the surface stiffness matrix on $\Gamma_h$. Let us mention that the solution of the $(N_p + 1)$ linear systems (55) and (56) is computed by an ILUT preconditioned GMRES and only requires 2 or 3 iterations to get the solution [35]. Therefore, applying the square-root operator results in a cost of the order $\mathcal{O}(N_V)$ which is less than applying an elementary integral operator (even with a FMM implementation).
<figure><img src="image_34.png" /><figcaption>Figure 9: Sound-hard unit square cylinder: number of MVPs for the four integral formulations with respect to the wavenumber *k* (left: $n_{\lambda} = 10$) and the density $n_{\lambda}$ of discretization points per wavelength (right: $k = 10$).</figcaption></figure>

We present now some numerical results to compare the efficiency of the GMRES without restart and for a tolerance equal to $\varepsilon = 10^{-6}$ when using the usual and generalized BW and CFIE. We report on Fig. 9 an example of computation for the scattering of an incident plane wave of incidence 45 degrees by the unit square cylinder. For the usual formulations, we can see that the number of MVPs grows linearly according to the density $n_{\lambda}$. Indeed, these two formulations are first-kind Fredholm. Furthermore, the number of iterations also depends on *k* (see the next example for the standard BW formulation). The problem is related to the fact that injecting the Sommerfeld approximation is not enough to handle the *k* dependence of the convergence. These two points can be corrected by considering the generalized formulations. These conclusions are confirmed for the three-dimensional case. In Fig. 10, we present the case of the scattering of an incident plane wave of incidence zero degree by an ellipsoidal scatterer of semi-axis *a* = 1 and *b* = *c* = 0.5 respectively along the $x_1$-, $x_2$- and $x_3$-directions and centered at the origin. Here again, the convergence of the GMRES is independent of *k* and $n_{\lambda}$ while this is not the case for the two standard formulations.

## 8.4 Other recent directions

Other recent and closely related approaches have been derived e.g. by Levadoux and al [2, 52], Bruno and al [20]. We refer to their papers for more details. Let us note that generalized combined field integral equations have also been proposed recently in the framework of electromagnetism for Maxwell's equations.

## 9 Conclusion

The aim of this chapter was to introduce some recent achievements into the topic of preconditioning of integral equations of acoustics. Most particularly, we discussed the
<figure><img src="image_35.png" /><figcaption>Figure 10: Sound-hard ellipsoidal scatterer: number of MVPs for the four integral formulations with respect to the wavenumber *k* (left: *n*<sub>λ</sub> = 10) and the density *n*<sub>λ</sub> of discretization points per wavelength (right: *k* = 5).</figcaption></figure>

high frequency regime which is a very active area where much work remains to do.
A selected review of some advances on analytical preconditioning has been proposed.
Different points of view can be adopted like for example the integral or pseudodiffer-
ential approaches. Many open questions are still to address on this problem and can
have implications on other close situations like electromagnetism or elasticity where
the numerical solution by means of integral equations is also extremely challenging.

Other future exciting directions of research include: understanding of the effect of non convexity and geometrical singularities on an iterative solver, coupling of analytical and algebraic preconditioners, application of pseudodifferential operator theory to the design of hybrid analytical-algebraic preconditioners for the finite element solution of PDEs...

# 10 Acknowledgements.

The authors are supported partially by the French ANR fundings under the project
MicroWave NT09_460489 (http://microwave.math.cnrs.fr/).

References

[1] M. Abramowitz and I.A. Stegun, *Handbook of Mathematical Functions*, Dover, New York, (1974).

[2] F. Alouges, S. Borel and D. Levadoux, *A stable well-conditioned integral equation for electromagnetism scattering*, J. Comput. Appl. Math. **204** (2) (2007), pp.440-451.
[3] S. Amini, On the choice of the coupling parameter in boundary integral formulations of the exterior acoustic problem, Appl. Anal. **35** (1-4), (1990), pp.75-92.

[4] S. Amini and S.M. Kirkup, Solution of Helmholtz equation in exterior domain by elementary boundary integral equations, J. Comput. Phys. **118** (1995), pp. 208-221.

[5] X. Antoine, Conditions de Radiation sur le Bord, Ph.D. Thesis, University of Pau, 1997.

[6] X. Antoine, Fast approximate computation of a time-harmonic scattered field using the On-Surface Radiation Condition method, IMA J. Appl. Math. **66** (2001), pp. 83-110.

[7] X. Antoine, Advances in the On-Surface Radiation Condition Method: Theory, Numerics and Applications, Book Chapter in Computational Methods for Acoustics Problems, Editor F. Magoulès, Saxe-Coburg Publications, 2008, pp.169-194 (ISBN: 978-1-874672-30-2).

[8] X. Antoine, H. Barucq and A. Bendali, Bayliss-Turkel-like radiation condition on surfaces of arbitrary shape, J. Math. Anal. Appl. **229** (1999), pp. 184-211.

[9] X. Antoine, A. Bendali and M. Darbas, Analytic Preconditioners for the Electric Field Integral Equation, Int. J. Numer. Meth. Eng. **61** (2004), pp. 1310-1331.

[10] X. Antoine, A. Bendali and M. Darbas, Analytic Preconditioners for the Boundary Integral Solution of the Scattering of Acoustic Waves by Open Surfaces, J. Comput. Acoust. **13** (3), (2005), pp. 477-498.

[11] X. Antoine and Y. Boubendir, An integral preconditioner for solving the two-dimensional scattering transmission problem using integral equations, Int. J. Comput. Math., **85** (10), (2008), pp.1473-1490.

[12] X. Antoine and M. Darbas, Alternative integral equations for the iterative solution of acoustic scattering problems, Quarterly J. Mech. Appl. Math. **58** (2005), pp. 107-128.

[13] X. Antoine and M. Darbas, Generalized combined field integral equations for the iterative solution of the three-dimensional Helmholtz equation, Math. Model. Numer. Anal., **41** (1), (2007), pp.147-167.

[14] X. Antoine, M. Darbas and Y.Y. Lu, An Improved Surface Radiation Condition for High-Frequency Acoustics Scattering Problems, Comput. Meth. Appl. Mech. Eng. **195** (33-36) (2006), pp. 4060-4074.

[15] X. Antoine, C. Geuzaine and K. Ramdani, Computational Methods for Multiple Scattering at High Frequency with Applications to Periodic Structures Calculations, M. Ehrhardt (Ed.), Wave Propagation in Periodic Media - Analysis, Numerical Techniques and Practical Applications, Progress in Computational Physics, Vol. 1, Bentham Science Publishers Ltd., 2009.

[16] A. Bayliss, C.I. Goldstein and E. Turkel, An iterative method for the Helmholtz equation, J. Comput. Phys., **49**, (1983), pp.443-457.

[17] M. Benzi, Preconditioning techniques for large linear systems: A survey, J. Comput. Phys. **182**(1), (2002), p.418-477.

[18] M. Bollhoefer, M.J. Grote and O. Schenk, Algebraic multilevel preconditioner
for the Helmholtz equation in heterogeneous media, SIAM J. Sc. Comput. 31
(5), (2009), pp.3781-3805.

[19] H. Brakhage and P. Werner, Über das Dirichletsche Aussenraumproblem für die Helmholtzsche Schwingungsgleichung, Numer. Math., **16**, (1965), pp.325–329.

[20] O.P. Bruno, T. Elling, and C. Turc, Well conditioned high-order algorithms for the solution of three-dimensional surface acoustic scattering problems with Neumann boundary conditions, submitted, 2009.

[21] O. Bruno and C. Geuzaine, An $O(1)$ integration scheme for three-dimensional surface scattering problems, J. Comput. Appl. Math., **204**, (2007), pp.463–476.

[22] O. Bruno, C. Geuzaine, J. Monro, Jr., and F. Reitich, Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case, Philosophical Transactions of the Royal Society (Series A: Mathematical, Physical and Engineering Sciences), **362** (1816), (2004), pp.629-645.

[23] O.P. Bruno and L.A. Kunyansky, High-order algorithm for the solution of surface scattering problems: Basic implementation, tests, and applications, J. Comput. Phys. **169**, (2001), pp.80–110.

[24] A. Buffa and R. Hiptmair, Regularized combined field integral equations, Numer. Math., **100** (1), (2005), pp.1-19.

[25] A. J. Burton and G. F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proc. Roy. Soc. London. Ser. A, **323**, (1971), pp.201-210.

[26] B. Carpentieri, I. S. Duff, and L. Giraud, Sparse pattern selection strategies for robust Frobenius-norm minimization preconditioners in electromagnetics, Numer. Linear Algebra Appl. 7 (2000), pp. 667-685.

[27] B. Carpentieri, I. S. Duff and L. Giraud, Experiments with sparse approximate preconditioning of dense linear problems from electromagnetic applications, Tech. Rep. TR/PA/00/04, CERFACS, France, (2000).

[28] S.N. Chandler-Wilde, I.G. Graham, S. Langdon and M. Lindner, Condition number estimates for combined potential boundary integral operators in acoustic scattering, J. Int. Eq. Appl. **21** (2) (2009), pp.229-279.

[29] W.C. Chew, J.M. Jin, E. Michielssen and J. Song, "Fast and Efficient Algorithms in Computational Electromagnetics", Artech House Antennas and Propagation Library, Norwood, 2001.

[30] S.H. Christiansen, Résolution des équations intégrales pour la diffraction d'ondes acoustiques et électromagnétiques. Stabilisation d'algorithmes itératifs et aspects de l'analyse numérique, PhD Thesis, Ecole Polytechnique, Palaiseau, France, (2001).

[31] P.G. Ciarlet, "Introduction à l'analyse numérique matricielle et à l'optimisation", Masson, Paris, 1982.

[32] R. Coifman, V. Rokhlin and S. Wandzura, "The Fast Multipole Method for the wave equation: A pedestrian description", IEEE Trans. on Ant. and Prop., **35**(3), (1993), pp.7-12.

[33] D. L. Colton and R. Kress, "Integral equation methods in scattering theory", Pure
and Applied Mathematics, John Wiley and Sons Inc., 1983.

[34] M. Costabel, *Integral equation methods in scattering theory*, SIAM J. Math. Anal. **19**, (1988), pp.613-626.

[35] M. Darbas, *Préconditionneurs analytiques de type Calderon pour les formulations intégrales des problèmes de diffraction d'ondes*, PhD Thesis, Toulouse, 2004.

[36] E. Darrigrand, *Couplage Méthodes Multipôles Rapides et Discrétisation Microlocale pour les Equations Intégrales de l'Electromagnétisme*, PhD Thesis, Bordeaux, 2002.

[37] E. Darve, *Méthodes multipôles rapides: Résolution des équations de Maxwell par formulations intégrales*, Ph.D. Thesis, Paris 6, 1999.

[38] Y. A. Erlangga, C. Vuik and C.W. Oosterlee, *On a class of preconditioners for the Helmholtz equation*, Appl. Numer. Math., **50**, (2004), pp.409-425.

[39] A. Greenbaum, *Iterative Methods for Solving Linear Systems*, Frontiers in Applied Mathematics, **27**, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.

[40] L. Greengard and V. Rokhlin, *A fast algorithm for particle simulations*, J. Comput. Phys., **73** (2), (1987), pp.325-348.

[41] M.J. Grote and T. Huckle, *Parallel preconditioning with sparse approximate inverses*, SIAM J. Sc. Comput. **18** (3), (1997), pp.838-853.

[42] P.L. Ho and Y.Y. Lu, *Improving the beam propagation method for TM polarization*, Opt. Quantum Electronics, **35**(4), (2003), pp. 507-519.

[43] D.S. Jones, *An approximate boundary condition in acoustics*, J. Sound Vibr. **121** (1) (1988), pp. 37-45.

[44] R. Kerchroud, A. Soulaimani and X. Antoine, *Performance study of plane wave finite element methods with a Padé-type artificial boundary condition in acoustic scattering*, Adv. Eng. Soft. **40** (2009), pp. 738-750.

[45] R. Kerchroud, A. Soulaimani and Y. Saad, *Preconditioning techniques for the solution of the Helmholtz equation by the finite element method*, Math. Comput. Simul., **65** (4-5), (2004), pp.303-321.

[46] R. Kress, *Minimizing the condition number of boundary integral operators in acoustic and electromagnetic scattering*, Quart. J. Mech. Appl. Math., **38** (2), (1985), pp.323-341.

[47] R. Kress, *Linear Integral Equations*, Second Ed., Applied Mathematical Sciences (82), New York, Springer-Verlag, 1999.

[48] R. Kress and W. T. Spassov, *On the condition number of boundary integral operators for the exterior Dirichlet problem for the Helmholtz equation*, Numer. Math., **42** (1), (1983), pp.77-95.

[49] G.A. Kriegsmann, A. Taflove and K.R. Umashankar, *A new formulation of electromagnetic wave scattering using the on-surface radiation condition method*, IEEE Trans. Antennas Propag. **35** (1987), pp. 153-161.

[50] S. Lang, *Linear Algebra*, Springer, 3rd edition, 2004.

[51] W. McLean, *Strongly Elliptic Systems and Boundary Integral Equations*, Cambridge University Press, Cambridge, 2000.
[52] D. P. Levadoux, *Etude d'une Equation Intégrale Adaptée à la Résolution Hautes Fréquences de l'Equation de Helmholtz*, PhD thesis, Université Paris VI, 2001.

[53] J. Liesen and P. Tichy, *Convergence analysis of Krylov subspace methods*, GAMM-Mitteilungen, **27**(2), (2004), pp. 153-173.

[54] F.A. Milinazzo, C.A. Zala et G.H. Brooke, *Rational square-root approximations for parabolic equation algorithms*, J. Acoust. Soc. Am. **101**(2), (1997), pp. 760-766.

[55] J.-C. Nédélec, *Acoustic and Electromagnetic Equations*, Applied Mathematical Sciences (144), New York, Springer-Verlag, 2001.

[56] D. Osei-Kuffuor and Y. Saad, *Preconditioning Helmholtz linear systems*, Applied Numerical Mathematics (in press, 2010).

[57] V. Rokhlin, *Rapid solution of integral equations of scattering theory in two dimensions*, J. Comput. Phys. **86**(2) (1990), pp. 414-439.

[58] Y. Saad, *Iterative Methods for Sparse Linear Systems*, PWS Publishing Company, Boston, 1996.

[59] Y. Saad and M.H. Schultz, *A generalized minimal residual algorithm for solving nonsymmetric linear systems*, SIAM J. Sci. Comput., **7**(3), (1986), pp.856-869.

[60] L. Schwartz, *Théorie des Distributions*, Broché, 1997.

[61] O. Steinbach and W. L. Wendland, *The construction of some efficient preconditioners in the boundary element method*, Adv. Comput. Math. **9**(1-2) (1998), pp. 191-216.

[62] G. Sylvand, *La Méthode Multipôle Rapide en Electromagnétisme: Performances, Parallélisation, Applications*, Thèse de Doctorat, ENPC, 2002.

[63] M.E. Taylor, *Pseudodifferential Operators*, Princeton Mathematical Series, 34, Princeton University Press, 1981, xi+452.

[64] L.L. Thompson, *A review of finite-element methods for time-harmonic acoustics*, J. Acoust. Soc. Amer., **119**(3), (2006), pp. 2272-2293.

[65] L. N. Trefethen, *Pseudospectra of Matrices*, editors D. F. Griffiths, G. A. Watson, in "Numerical Analysis", Longman Scientific and Technical, 1991.

[66] E. Turkel, *Boundary Conditions and Iterative Schemes for the Helmholtz Equation in Unbounded Regions*, Book Chapter in Computational Methods for Acoustics Problems, Editor F. Magoulès, Saxe-Coburg Publications, 2008, pp. 127-159 (ISBN: 978-1-874672-30-2).

[67] Y. Wang, J. Lee and J. Zhang, *A short survey on preconditioning techniques for large-scale dense complex linear systems in electromagnetics*, Int. J. Comp. Math. **84**(8) (2007), pp. 1211-1223.