Keywords: multiple scattering, acoustic scattering, small obstacles acoustic scattering, single layer methods, integral equation methods, preconditioning in multiple scattering

We investigate efficient numerical methods for the problem of multiple-scattering of obstacles in homogeneous media. This is a first step towards the more general problem in strongly inhomogeneous media. The inhomegeneity for the multiple-scattering problem is caused by the presence of obstacles. For formations composed of a small number of medium-sized obstacles, satisfactory results can be obtained with optimized softwares based on standard discretization technique such as Finite Element Method (FEM). However, constraint by its need for meshing, a FEM loses its robustness, as the number of obstacles increases, or when their size decreases. As an alternative, we work with a Galerkin Integral Equation method, which we call Fourier Series -Single Layer (FS-SL) method, which describes the scattered wave as a superposition of single layer potentials and uses truncation of Fourier series to discretize the continuous problem. We describe in details the systems generated by the method, accompanied by a well-posedness study, for penetrable and impenetrable obstacles, the later involving Dirichlet, Neuman and Impedance boundary conditions. To study the numerical performance of the method, we limit ourselves to the case of disc -shaped obstacles. We first compare the results of our mesh-free method with Montjoie (a FE-based software) to validate the robustness for problems with a large number of small obstacles. We then investigate then efficiency of different solver types in the resolution of the dense linear system generated by FS-SL method. The study is done for Direct Solvers (Mumps, Lapack and Scalapack) and iterative GMRES-type Solvers with various preconditioners. We show that the optimal choice depends on the distance between obstacles, their size and number.

Une étude sur l'efficacité numérique de la méthode de potentiel de simple couche avec une base de Fourier pour la diffraction multiple par des obstacles dans un milieu homogène.

Résumé :

Nous nous intéressons aux méthodes numérique pour simuler avec efficacité la diffraction multiple d'une onde acoustique par des obstacles dans un milieu homogène. Ce project constitute un première pas vers un problème plus géneral incluant un milieu fortement héterogène. L'inhomogénéité dans le cas présent se caractérise par la présence des obstacles. Lorsque le milieu contient peu d'obstacles et qu'ils sont de taille moyenne, les méthodes numériques basées sur les techniques de discrétization telles que les élements finis (EF) sont efficaces. En revanche, lorsque le nombre d'obstacles augmente, ou que leur taille diminue, de telles méthodes, de par leur besoin en maillage, perdent en performance. Comme alternative, nous travaillons sur une méthode de type Galerkin Équation Intégrale, que nous appellons 'Fourier Series -Single Layer' (FS-SL). La méthode décrit l'onde diffractée comme une superposition des potentiels de simple couche, et utilise la troncature des Séries de Fourier pour discrétiser le problème continu. Dans ce rapport, nous donnons les systèmes ainsi engendrés par la méthode, accompagnés d'une étude détaillée sur leur comportement du type Fredholm, dans le cas d'obstacles pénétrables et impénétrables. Nous proposons ensuite une étude des performances numériques de notre méthode dans le cas d'obstacles circulaires. Nous comparons dans un premier temps notre méthode, qui ne nécessite pas de discrétisation spatiale (maillage), avec une méthode d'élements finis implémenteée dans Montjoie. Nous étudions aussi l'éfficacité de plusieurs types de solveurs, pour la résolution du système linéaire plein généré par la méthode FS-SL. Nous comparons les solveurs directs (Mumps, Lapack et Scalapack) et de type GMRES avec plusieurs préconditionneurs. Nous montrons que le choix optimal de solveur dépend de la distance entre les obstacles, leur taille et leur nombre.

Mots-clés : diffraction multiples, diffraction multiple par des disques, diffraction des petites obstacles, potentiel de simple couche, solveur itératifs.

Introduction

This project is part of a program which explores efficient numerical methods to solve direct and inverse problems for the propagation of acoustic wave in strongly inhomogeneous media in low-frequency regime. We work in close collaboration with the acoustic research lab I2M (at l'Université de Bordeaux), aiming at comparing data predicted by numerical simulations with those obtained from physical experiments (carried out by I2M). This goal necessitates the development of robust 'in-house' codes that can be evolved according to our needs, especially for applications to solving inverse problems, e.g. to detect defects in materials.

We start our investigation with inhomogeneities created by compactly-supported and nonoverlapping obstacles (also called scatterers); the overall problem is thus called the multiple scattering of obstacles in homogeneous media. The usual tools on small domains with a small number of sizable obstacles are highly optimized softwares based on Finite Element (FE) Method, e.g. Montjoie (montjoie.gforge.inria.fr). However, when the size of the domain of interest and number of obstacles increase, or when the size of the obstacles decreases, a FE method generates large linear systems due to its need of extremely refined mesh. For example, when the obstacles are very small compared to the wavelength of the incident waves, in order to capture the interaction between them, the mesh has to be refined, at least in the region around the obstacles. In addition, the meshing of the obstacles may not be straightforward and requires some effort, in particular when there is a large number of them. Finally, any FE method requires finite domain of calculation and thus a numerical method for domain truncation (e.g. by techniques of absorbing boundary conditions or perfectly matched layer), which increases the level of technicality in its implementation.

To overcome these challenges, we study a method in the family of Galerkin Boundary Integral Equations (BIE) methods, which we will call Fourier Series -Single Layer method (FS-SL). The method does not require a mesh, and thus allows to study the scattering problem on infinite domains and with very small obstacles. The single layer part comes from viewing the scattered response as a superposition of waves, which are scattered by each of the obstacles and expressed as acoustic Single-Layers (SL). In particular, for a formation of N Obs non-overlapping obstacles with boundary denoted by Γ J , 1 ≤ J ≤ N Obs , we write

u scatt = N Obs J=1 SΓ J V J , V J ∈ C(Γ J ).
Here, the acoustic single-layer SΓ I with density φ ∈ C(Γ I ) are defined as

SΓ I φ (x) := Γ I φ(y) G κ (x, y) dσ(y) , x ∈ R 2 \ Γ I ,
with G κ being the fundamental solution of the Helmholtz equation at wavenumber κ 2 , G κ (x, y) := i 4 H

(1) 0 (κ |x -y|) , x = y .

With this formulation, the unknowns of the problem are now the family of SL densities {V J }.

The Fourier Series part (in the name of the method) indicates the choice of Galerkin spaces, for which the continuous densities V J -s are approximated by V J,h -s which are obtained from generalized Fourier series approximation on simply-closed continuous curves,

V J,h = m k=-m v J,k SΓ J w J,k .

Inria

Here, m is the order of approximation, and the unknowns of the discrete problem are {v J,k } the generalized Fourier Series coefficients of SL densities V J,h , c.f. Subsection 3.6 for more details. We choose to work with the single-layer operator S, since it gives rise to the simplest1 integral equations. Moreover, the intrinsic problem of S regarding invertibility (and hence instability) does not arise with small obstacles2 , e.g. the formation with the exterior wavenumber κ 2 e satisfying κ e r circumscribed circle of Ω I < 2 , 1 ≤ I ≤ N Obs .

Together with the single layer Ansatz, the choice of Fourier Series basis is natural in the following sense. When there is only one obstacle which is disc-shaped, we obtain the well-known exact solution in separated variables (polar coordinates); when there are more than one obstacles, each disc-shaped, we obtain the same solution given by Multipole theory, c.f. Proposition 7 and Remark 5.

We make a brief digression to compare our (FS-SL) method with other approaches usually found in the vast literature of multiple obstacle scattering, e.g. various techniques of asymptotics, Foldy-Lax, Born-approximation, Fast Multipole, etc.; we also refer the readers to [START_REF] Martin | Multiple Scattering : Interaction of Time-Harmonic Waves with N Obstacles[END_REF] for an extensive exposition on the history and reviews of the subject. The Foldy-Lax model [START_REF] Foldy | The multiple scattering of waves[END_REF], which approximates each small obstacle by an isotropic point scatterer, corresponds to our method (and the Multipole Method) at first order approximation, i.e with m = 0. In the coefficient matrix of the linear system, c.f. Subsection 3.2, the off-diagonal blocks describe the interaction between different obstacles, while the diagonal blocks the self-interaction within an obstacle. If we ignore the off-diagonal blocks, we obtain the Born-series approximation, c.f. [START_REF] Marengo | Subspace-based localization and inverse scattering of multiply scattering point targets[END_REF]6]); in other words, we are treating the problem as 'single-scattering', an approximation which does not describe adequately the corresponding physical phenomena, unless the distances between the obstacles are much larger than their sizes and the wavelength of the incident wave. An improved version of the single-scattering approximation is given by Generalized Born Series, c.f. [START_REF] Martin | Multiple scattering: an invitation[END_REF], which is equivalent to using GMRES-solver with block Jacobi preconditioner, the efficiency of which will be discussed below. There is also matched asymptotic expansion technique, which gives a comprehensive interaction between the obstacles, c.f. [4] in dimension 3; however, a similar attempt in dimension 2 would be technically cumbersome, since one will have to deal with power series in the radius variable and its logarithm. The use of integral equation technique in multiplescattering is not new, several variants can be found in the literature. Towards the end of this project, we learnt of a series of work done by [3], [START_REF] Thierry | µ-diff: an open-source matlab toolbox for computing multiple scattering problems by disks[END_REF], also using single-layer technique to study the multiple-scattering problem at low and high frequencies. Our work can be considered as a complement to their work; in addition, we offer further extensive numerical experiments, whose parameters are more relevant to our research questions. Finally, as mentioned at the beginning, the need for an in-house code that can be adapted at will and ease to our research needs rest unchanged, despite numerous ideas that might already exist in literature.

We next discuss the range of applicability of our method. The FS-SL method ( in fact layer operator technique in general) is applicable to obstacles of arbitrary (C 2 and convex) shapes in arbitrary formation. In Section 3, we list the detailed linear systems satisfied by impenetrable obstacles (with Dirichlet, Neumann and Impedance boundary conditions) and penetrable ones, together with a study/review of the Fredholmness of these linear systems. On the other hand, for immediate and practical need to compare with physical experiments, we devote our numerical study to periodic formations of disc-shaped obstacles; however, the method and the codes do not make use of the periodicity in the formation. As mentioned previously, in the case of disc geometry, the FS-SL method coincides with the Multipole method; as a result, one has analytical expressions for linear system, giving the method an efficiency boost, c.f. Subsection 4.1. In order to maintain the robustness observed for disc geometry, an efficient integration quadrature rule will be needed to handle the weak singularity of the Green kernel, c.f. (44). For the case of one obstacle, modified integral equation and Nystrom quadrature have been used, c.f. [START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF] and the references therein.

The code is written in Fortran 90 in double precision, and uses a parallel architecture 3 . Users have the option to choose a solver type among the Direct Solvers (Mumps, Lapack and Scalapack), or the iterative GMRES-type Solvers [START_REF] Frayssé | A set of gmres routines for real and complex arithmetics on high performance computers[END_REF] with various preconditioners, and between modes of post-processing (exact or interpolation). Below we discuss the comparative reviews of these solvers. In order to discuss in details the efficiency of the FS-SL method, we distinguish between pre-processing and post-processing time. By pre-processing time, we mean the CPU time needed for the construction of the linear system, whose unknowns are the generalized Fourier Series coefficients {v J,k }, and its resolution. Unlike the linear system for FEM-s, those for BIE methods are generally dense. However, this is not a problem, since the size in IE methods are generally much smaller; for an approximation order of m, the size of the FS-SL system is

(Number of obstacles N Obs ) × (2 × m + 1) .
Hence, the denseness does not present difficulty, especially when the problem can be handled with efficient solvers like Lapack and Scalapack; e.g. , with parallelization and parallel solvers Scalapack, we have tested up to 10 4 small and closely-spaced disc-shaped obstacles, c.f. Subsection 6.6. In addition, even with low order of approximation (< 4), one is already in an acceptable precision range (around 10 -7 ), c.f. the numerical convergence in Subsection 4.3, Section 5 and 6.

The post-processing time depends on the purpose of usage. The cost depends on the number of points of evaluation, since one needs to evaluate the single layers for each of these points. Most expensive for a BIE-type method is a visualisation of solution on a two-dimensional grid (with the cost increasing with the degree of visualisation resolution). On the other hand, for the purpose of inverse problems, for which one is mostly concerned with the far field pattern, or the value on an at most one-dimensional curve, the evaluation cost is very low. It should be noted that the two processes are decoupled, giving the flexibility in the sense that one is not constrained to a fixed domain of interest. More specifically, for a fixed formation of obstacles, once the SL densities {V h,J } are obtained and saved, a user can obtain the value of the field at whichever point in the infinite domain, and thus can either zoom or enlarge the domain of visualisation, with varying degree of resolution. In the case of disc geometry, the post-processing time involves the evaluation of Hankel functions on the visualization grid. The drawback of the FS-SL method can be greatly diminuished, if we use interpolation, e.g. Hermite cubic spline c.f. [START_REF] Hairer | Lecture notes introduction à l'analyse numérique[END_REF][p.48-50], together with a parallelization of the post-processing codes, c.f. the numerical tests in Section 5 and 6.

The numerical comparisons are separated in two groups. For the first one, we validate the observations made above regarding the shortcomings of a FE-based method for our settings, and compare the performance between our method (with direct solver Mumps) and the optimized software Montjoie, c.f. Section 5. For the second group, we address the questions whether a direct or an iterative solver is better, and within the iterative family, whether preconditioning is needed, and if this is the case, which kind of preconditioner and position of preconditioning will give the fastest convergence, c.f. Section 6. We will see that different formations of obstacles require different solvers to obtain optimality or even just convergence (for the iterative ones).

The comparison is done for the list of preconditioners listed in Appendix C.4, and thus offers an extensive complement to the study done in [3] which mentioned two preconditioners, the first one being the block Jacobi and the second one comparable to our 2nd-order Jacobi4 . Although the usage in [3] is aimed for high frequency, one encounters the same problem in low frequency, regarding the convergence of the iterative solvers, when there are large number of obstacles which are closed together, as is noted in [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circuler cylinders[END_REF].

For closely-spaced obstacles, we observe that the direct solvers outperform the iteratives ones, with Lapack and its parallel version Scalapack leading in efficiency. For our setting, we also observe a need to precondition the system. Among the preconditioners, the Jacobi family of preconditioners, which contains only the information of diagonal blocks, dooes not perform as well as Gauss-Seidel type, which contain informations of off-diagonal blocks. The first group has difficulty in attaining convergence, even for 200 obstacles, c.f. Subsection 6.1-6.4. Intuitively, with the strong interaction among closely-spaced obstacles, the coefficient matrix ceases to be diagonally dominant, with the off-diagonal blocks (describing interaction between different obstacles) being comparable in size to the diagonal ones (describing self-interaction). We note the prominent robustness of the Lower-upper Symmetric Gauss-Seidel (LU-SGS) and Symmetric Gauss-Seidel (SGS) preconditioners above the rest. When the obstacles are further apart, the numbers of iterations needed for the GMRES solvers drastically drop. We observe that LU-SGS and SGS take less time than Lapack, and are almost comparable in performance to Scalapack, c.f. the experiments in Subsection 6.5. It should be noted that the current codes for the GMRES solvers are still sequential. This means that any parallelization and optimization can reduce the time cost even further, making the iterative solvers a promising candidate for cases where Scalapack or Lapack will fail, e.g. beyond 10 4 sparsely-spaced obstacles.

Mathematical Statement and Well-posedness of the continuous problems

The propagation of a time-harmonic acoustic wave of frequency ω 2π in a homogeneous medium with wavespeed c is described by u inc (x) e iωt , where u inc satisfies the Helmholtz equation

(-∆ -κ 2 e ) u inc = 0 in R 2 ,
with the wavenumber κ e which satisfies the dispersion relation κ e = ω c . In the presence of obstacles, the above incident wave is scattered and becomes a superposition of the incident wave and a scattered response u, see Figure 1,

u total = u inc + u .
We require both the total field u total and scattered one u to solve the Helmholtz equation (with the wavenumber κ e ) in the region exterior to the obstacles. The scattered field u also needs to satisfy the κ e -outgoing condition at infinity,

lim r→∞ √ r (∂ r u -i κ e u) = 0 ; r = |x| . (1)
This condition assures that the scattered wave does not re-enter the domain of interest. How the obstacles reflect the incident wave is prescribed as boundary conditions posed along their boundary. To describe concretely these boundary conditions, we first specify the geometry We assume that there are no overlapping between them, i.e Ω I ∩ Ω J = ∅ for I = J. The region outside the obstacles is denoted by Ω ext , the region within Ω int , and their common boundary Γ Obs ,

Ω ext u inc us c a t t
Ω ext = R 2 \ Ω int ; Ω int = N Obs I=1 Ω I ; Γ Obs = Ω int ∩ Ω ext = N Obs I=1 Γ I .
The boundary conditions, imposed on Γ Obs , reflect how the obstacles interact with an incident wave. For impenetrable cases, the waves are fully reflected, e.g. when the obstacles are holes. The boundary value problems in this case are further classified as follows. For sound-soft obstacles, the normal velocity of the total wave vanishes on the boundary, which leads to the Dirichlet condition,

u total Γ Obs = 0 . (2)
For sound-hard obstacles, we impose the Neumann condition,

∂ n u total Γ Obs = 0 . (3)
The more general setting is when the normal velocity is proportional to the wave on the boundary, which leads to an impedance condition,

(∂ n + i λ ) u total Γ Obs = 0 . (4)
For the penetrable cases, waves are both reflected and transmitted inside the obstacles. Denote by u i the wave field in Ω int . In addition to the conditions imposed on the scattered wave u, the transmission problem comprises in addition of

(-∆ 2 -κ 2 i,J ) u i = 0 , x ∈ Ω int u total = u i ; ∂ n u total = µ ∂ n u i , x ∈ Γ Obs . (5)
In short, the multiple scattering problem by N Obs obstacles of a time-harmonic incident wave is formulated as an exterior boundary-value problem (BVP) for impenetrable obstacles where the wave is fully reflected, and as a transmission one for penetrable obstacles where reflection is only partial. For the sake of clarity, we restate the above problems in terms of the scattered wave u and the traces operators. We first cite list the definitions and sign convention for the trace operators. The normal vector n(x) is chosen to point outward. For f ∈ H 2 , define the normal derivative associated to normal vector n,

∂ ∂n f = lim h→+0 n(x) • ∇f (x -h n(x)).

Inria

In terms of γ 0,int and γ 0,int , this can be written as,

γ 1,int f := (γ 0,int ∇f ) • n , γ 1,ext f := (γ 0,ext ∇f ) • n.
The jump at an interface is defined as,

f := γ 0,ext f -γ 0,int f ; ∂ n f := γ 1,ext f -γ 1,int f = ∇f • n.
Boundary value problems :

     (-∆ -κ 2 e ) u = 0 , in Ω ext γ 0,I,ext (u + u inc ) = 0 , 1 ≤ I ≤ N Obs lim r→∞ √ r (∂ r u -i κ e u) = 0 , r = |x| Exterior Diriclet Problem (EDP) (6)
We have denoted by λ the impedance parameter

     (-∆ -κ 2 e ) u = 0 , in Ω ext (γ 1,I,ext + iλγ 0,I,ext ) (u + u inc ) = 0 , 1 ≤ I ≤ N Obs lim r→∞ √ r (∂ r u -i κ e u) = 0 , r = |x| Exterior Impedance Problem (ENP) (7)
When λ = 0, we have the Neumann problem (ENP). . For the transmission problem, we look for the scattered field u in Ω ext and the interior transmitted ones u I in Ω I for 1 ≤ I ≤ N Obs . We have denoted by µ the transmission parameter.

               (-∆ -κ 2 e ) u = 0 , in Ω ext (-∆ -κ 2 I ) u I = 0 , in Ω I , 1 ≤ I ≤ N Obs γ 0,I,ext (u + u inc ) = γ 0,I,int u I , 1 ≤ I ≤ N Obs γ 1,I,ext (u + u inc ) = µ γ 1,I,int u I , 1 ≤ I ≤ N Obs lim r→∞ √ r (∂ r u -i κ e u) = 0 , r = |x| Transmission Problem (TP) (8)
Well-posedness of the boundary value problems : We summarize the uniqueness results from [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF] and [START_REF] Kleinman | On single integral equations for the transmission problem of the acoustics[END_REF]. We note that the well-posedness results in [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF] allow for multi-component interior domains, as in our setting, c.f. [8][Section 2.1].

Theorem 1 (Uniqueness). Hypothesis for the boundary :

Γ Obs ∈ C 2 .
1. The Exterior Dirichlet (EDP) (6) and Neumann problems (ENP) have at most one solution.

The Exterior Impedance Problem (EIP)(7) has at most one solution if

Im (κ e i λ) ≥ 0 , on Γ Obs .

3.

Let κ e be such that κ e ∈ R + or Im κ e > 0. Let κ I and µ be such that

µ = 0 , Im (µ κ e ) ≥ 0 , Im (µ κ e κ 2 I ) ≥ 0 . ( 9
)
Then the transmission problem (TP) (8) has at most one solution.

Proof.

1. For EDP and ENP, we refer [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF][Thm 3.13].

2. For EIP, we refer to [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF][Thm 3.37].

3. The statement is verbatim the Uniqueness Theorem from [START_REF] Kleinman | On single integral equations for the transmission problem of the acoustics[END_REF].

Remark 1. As noted in Example 1 [START_REF] Kleinman | On single integral equations for the transmission problem of the acoustics[END_REF], the conditions for uniqueness are satisfied when all of the variables are real. Specifically, we will be working the following situation

0 < κ e < ∞ , 0 ≤ κ I < ∞ , 0 < µ < ∞ ,
and in this case, condition (9) is satisfied.

We can show the well-posedness (existence and uniqueness) of the above problems by first reducing to ones on a bounded domain. This is achieved by replacing the outgoing radiation condition with an exact boundary condition

∂ ∂ν u = T u , on ∂ B R , (10)
placed at the boundary of a disc B R of size R, with R large enough so that B R contains all of the obstacles. Here, we have used the Dirichlet-to -Neunmann map (DtN) T, defined as,

T : H 1/2 (∂B R ) → H -1/2 (∂B R ) bounded g → ∂ ∂ν w ,
where w is the unique solution the exterior Dirichlet problem in R 2 \ B R with boundary data w| ∂B R = g. Note that this solution is obtained in the form of a series expansion involving Hankel functions, by taking advantage of the disk geometry of B R . For proof of the boundedness of the mapping, see e.g.. and do the same for the first order traces and the interior versions. Since we assume that the incidence wave u inc is smooth, its zero-th and first order trace along a C 2 curve Γ is also smooth from both sides; as a result, we will drop the distinction 'int' and 'ext' from the notation of its traces, and simply note γ 0,Γ u inc and γ 1,Γ u inc .

Inria

We list the operators from potential theory that we will use. We following the notation of Colton Kress theorem 3.1 , p 39 [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF], see also [5][Section 7.1]. The acoustic single-layer and double-layer potential with density φ ∈ C(Γ I ) are defined correspondingly as,

SΓ I ,κ φ (x) := Γ I φ(y) G κ (x, y) dσ(y) , x ∈ R 2 \ Γ I ; DΓ I ,κ φ (x) := Γ I φ ∂ ∂n(y) G κ (x, y) dσ(y) , x ∈ R 2 \ Γ I ,
where G κ is the fundamental solution of the Helmholtz equation at wavenumber κ 2

G κ (x, y) := i 4 H (1) 0 (κ |x -y|) , x = y .
To describe the traces of the layer potentials (along the interface where it is defined), we will need the following surface operators.

(S Γ I ,κ φ) (x) := Γ I φ(y) G κ (x, y) ds(y) , x ∈ Γ I ; (D Γ I ,κ φ) (x) := Γ I φ ∂ ∂n(y) G κ (x -y) dσ(y) , x ∈ Γ I ; D Γ I ,κ φ (x) := Γ I φ ∂ ∂n(x) G κ (x -y) dσ(y) , x ∈ Γ I ; (T Γ I ,κ φ) (x) := ∂ ∂n(x) Γ I φ ∂ ∂n(y) G κ (x -y) dσ(y) , x ∈ Γ I .
To simplify the notation, we write SI,κ := SΓ I ,κ ; S I,κ = S Γ I ,κ , and do the same for surface operators D, D and T.

We will need to extend the notations to the multiple obstacles setting. Since the kernel of SI,κ is smooth for x ∈ R 2 \ Γ I , for φ ∈ H -1/2 (Γ I ), we have SI,κ φ ∈ H 1 (R 2 \ Γ I ), c.f. Appendix A. As a result, the function and its normal derivative are continuous across Γ J for J = I, γ 0,Γ J Sκ,I φ = 0 ;

γ 1,Γ J Sκ,I φ = 0 . (11)
For J = I, and φ ∈ H -1/2 (Γ J ), we define

S IJ,κ φ := γ 0,I Sκ,J φ ; D IJ φ := γ 1,I Sκ,J . (12)

Linear systems

As announced in the introduction, we use a single-layer potential Ansatz to describe the scattered response u = u scatt ,

u(x) := N Obs J=1 SΓ J ,κe ṽJ . ( 13
)
We refer the readers to Appendix A for a short summary of the layer theory as well as the definitions for operators involved.

Using the trace identities [START_REF]3 GMRES with right conditioning[END_REF], we obtain the exterior zero-th order trace along Γ I , γ 0,I,ext u =

Similarly, from [START_REF]3 GMRES with right conditioning[END_REF], we obtain the exterior first order trace along Γ I ,

γ 1,I,ext u = N Obs J=1 γ 1,I,ext SJ,κe ṽJ = D I,κe - 1 2 Id ṽI + N Obs J=1 J =I D IJ,κe ṽJ . ( 15
)
Dirichlet problem : u of the form (13) solves the EDP (6), if and only if {ṽ I } satisfies correspondingly of following equations, for each I with 1 ≤ I ≤ N Obs , γ 0,I,ext u = -γ 0,I,ext u inc .

Applying [START_REF] Kress | Linear Integral Equation[END_REF], we obtain

N Obs J=1 S IJ,κe ṽJ = -γ 0,I u inc . ( 16
)
Impedance problem : u of the form (13) solves the EIP [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], if and only if {ṽ I } satisfies correspondingly of following equations, for each I with

1 ≤ I ≤ N Obs , (γ 1,I,ext + i λ γ 0,I,ext ) u = -(γ 1,I,ext + i λ γ 0,I,ext ) u inc .
Applying ( 14) and [START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF], we obtain,

D I,κe - 1 2 Id ṽI + N Obs J=1 J =I D IJ,κe ṽJ + i λ N Obs

J=1

S IJ,κe ṽJ = -(γ 1,I,ext + i λ γ 0,I,ext ) u inc .

After rearrangement, the above equality becomes,

D I,κe - 1 2 Id + i λ S I,κe ṽI + J=1 J =I D IJ,κe + i λ S IJ,κe ṽJ = -(γ 1,I + i λ γ 0,I ) u inc . ( 17
)
Neumann problem : When µ = 0 in [START_REF] Liesen | Krylov subspace methods: principles and analysis[END_REF], we obtain the Neumann problem,

D I,κe - 1 2 Id ṽI + N Obs J=1 J =I D IJ,κe ṽJ = -γ 1,I u inc . ( 18
)
Transmission problem : To describe the wave inside obstacle I, we use the Green's representation c.f. [7, Thm 2.1]. We first rewrite the Green's representation in terms of the surface operators and the trace operators along Γ I ,

u I := u i | Ω I = SI,κ I • γ 1,I,int u i -DI,κ I • γ 0,I,int u i .
We next use the transmission conditions to express the traces of u I in terms of those of the exterior wave u total , γ 1,I,int u I = µ -1 γ 1,I,ext u total ; γ 0,I,int u I = γ 0,I,ext u total .

Inria

Substituting the above expressions into the Green's representation for u I , we get

u I = µ -1 SI,κ I • γ 1,I,ext u total -DI,κ I • γ 0,I,ext u total .
Next we replace u total by u total = u inc + u,

u I = µ -1 SI,κ I • γ 1,I,ext u -DI,κ I • γ 0,I,ext γ 0,e u + f . (19)
Here, we have denoted by f ,

f = µ -1 SI,κ I • γ 1,I,ext u inc -DI,κ I • γ 0,I,ext u inc = µ -1 SI,κ I • γ 1,I u inc -DI,κ I • γ 0,I u inc ⇒ γ 0,I,int f = µ -1 (S I,κ I • γ 1,I ) u inc -D I,κ I -1 2 Id (γ 0,I u inc ) . (20)
Taking the trace from the interior along Γ I on both sides of [START_REF] Martin | Multiple scattering: an invitation[END_REF], we obtain

γ 0,I,int u I = µ -1 γ 0,I,int • SI,κ I (γ 1,I,ext u) -γ 0,I,int • DI,κ I (γ 0,I,ext u) + γ 0,I,int f = µ -1 S I,κ I (γ 1,I,ext u) -D I,κ I -1 2
Id (γ 0,I,ext u) + γ 0,I,int f . Now replace the traces of u by that of its single layer Ansatz, ( 14) and ( 15), we get

γ 0,I,int u I = µ -1 S I,κ I D I,κe - 1 2 Id ṽI + µ -1 S I,κ I N Obs J=1,J =I D IJ,κe ṽJ -D I,κ I -1 2 Id N Obs

J=1

S IJ,κe ṽJ + γ 0,I,int f .

On the other hand, by the first transmission condition, the LHS can be written in terms of u total , γ 0,I,int u I = γ 0,I,ext u total = γ 0,I,ext u + γ 0,I,ext u inc = N Obs J=1 S IJ,κe ṽJ + γ 0,I u inc .

Thus replace the LHS of ( 21) by [START_REF] Thierry | µ-diff: an open-source matlab toolbox for computing multiple scattering problems by disks[END_REF], and after some simplification, we get

µ -1 S I,κ I • D I,κe - 1 2 Id ṽI - D I,κ I + 1 2 Id • N Obs J=1 S IJ,κe ṽJ + µ -1 S I,κ I • N Obs J=1,J =I D IJ,κe ṽJ = -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I + 1 2 Id (γ 0,I u inc ) .
This is rewriten as,

µ -1 S I,κ I • D I,κe - 1 2 Id - D I,κ I + 1 2 Id • S I,κe ṽI + N Obs J=1,J =I µ -1 S I,κ I • D IJ,κe -D I,κ I + 1 2 Id • S IJ,κe ṽJ = -γ 0,I,int f + γ 0,I u inc . (23)

RR n°8988

We further rewrite the RHS, using expression (20)

-γ 0,I,int f + γ 0,I u inc = -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I -1 2 Id (γ 0,I u inc ) + γ 0,I u inc = -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I + 1
2 Id (γ 0,I u inc ) .

We finally arrive at the linear equation for the transmission problem along Γ I ,

µ -1 S I,κ I • D I,κe - 1 2 Id - D I,κ I + 1 2 Id • S I,κe ṽI + N Obs J=1,J =I µ -1 S I,κ I • D IJ,κe -D I,κ I + 1 2 Id • S IJ,κe ṽJ = -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I + 1 2 Id (γ 0,I u inc ) . (24)
When N Obs = 1, this corresponds to Equation (5.25) (referred to as the MVM) in [START_REF] Kleinman | On single integral equations for the transmission problem of the acoustics[END_REF]. We have shown that if {u, u 1 , . . . , u N Obs } with u defined by

u = N Obs J=1 SJ,κe v J ; u I = µ -1 SI,κ I • γ 1,I,ext (u + u inc ) -DI,κ I • γ 0,1,ext (u + u inc ) , 1 ≤ I ≤ N Obs , (25)
solves the transmission problem [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF] with incident source u inc , then {ṽ I } has to satisfy linear problem (24). The following proposition states the reverse direction. In other words, the transmission problem, using single-layer Ansatz for the external field and Green's representation theorem for the internal ones, is equivalent to linear problem (24). We rearrange the equation satisfied by (v I ) 1≤I≤N Obs to obtain that

0 = µ -1 S I,κ I • D I,κe -1 2 Id - D I,κ I + 1 2 Id • S I,κe ṽI + N Obs J=1,J =I µ -1 S I,κ I • D IJ,κe -D I,κ I + 1 2 Id • S IJ,κe ṽJ -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I + 1 2 Id (γ 0,I u inc ) .
The RHS in the above equality is exactly γ 0,I,int u I -γ 0,I,ext (u + u inc ), which can be seen as Inria follows,

γ 0,I,int u I -γ 0,I,ext (u + u inc ) = µ -1 S I,κ I γ 1,I,ext u -(D I,κ I -1 2 Id)γ 0,I,ext u + µ -1 S I,κ I γ 1,I u inc -(D I,κ I -1 2 Id)γ 0,I u inc -γ 0,I,ext (u + u inc ) = µ -1 S I,κ I γ 1,I,ext u -(D I,κ I + 1 2 Id)γ 0,I,ext u + µ -1 S I,κ I γ 1,I u inc -(D I,κ I + 1 2 Id)γ 0,I u inc = µ -1 S I,κ I D I,κe -1 2 Id ψ I + µ -1 S I,κ I N Obs J=1,J =I D IJ,κe ψ J -D I,κ I + 1 2 Id N Obs J=1 S IJ,κe ψ J + µ -1 (S I,κ I • γ 1,I ) u inc -D I,κ I + 1 2 Id (γ 0,I u inc ) .
(26) Hence, we have shown that the first transmission condition is satisfied, i.e,

γ 0,I,int u I = γ 0,I,ext u , 1 ≤ I ≤ N Obs .
We next verify the second transmission condition. For 1 ≤ I ≤ N Obs , extend u I to outside of Ω I and denote this function by w I . Being a linear combination of a single layer and a double layer at wavenumber κ I , we have that w I solves the Helmholtz equation and satisfies the outgoing radiation condition (1), at wavenumber κ I . We next consider its external traces at Γ I .

γ 0,I,ext w I = µ -1 S I,κ I γ 1,I,ext u -D I,κ I + 1 2 Id • γ 0,I,ext u + µ -1 S I,κ I γ 1,I,ext u inc -D I,κ I + 1 2 Id • γ 0,I,ext u inc .
Comparing with the first equality of (26), we obtain that γ 0,I,ext w I = γ 0,I,int u I -γ 0,I,ext (u + u inc ) = 0 . This means that w I satisfies the exterior Dirichlet boundary value problem for the Helmholtz equation with zero Dirichlet boundary value. As a result,

w I ≡ 0 in Ω ext ⇒ γ 1,I,ext w I = 0 .
On other hand, from its definition, we have

γ 1,I,ext w I = µ -1 D I,κ I -1 2 Id γ 1,I,ext u -T I,κ I • γ 0,I,ext u + µ -1 D I,κ I -1 2 Id γ 1,I u inc -T I,κ I • γ 0,I u inc . The RHS is exactly γ 1,I,int u I -µ -1 γ 1,I,ext (u + u inc ), which can be seen as follows, µ γ 1,I,int u I -γ 1,I,ext (u + u inc ) = D I,κ I + 1 2 Id γ 1,I,ext u -µT I,κ I • γ 0,I,ext u + D I,κ I + 1 2 Id γ 1,I u inc -µT I,κ I • γ 0,I u inc -γ 1,I,ext (u + u inc ) = D I,κ I -1 2 Id γ 1,I,ext u -µT I,κ I • γ 0,I,ext u + D I,κ I -1 2 Id γ 1,I u inc -µT I,κ I • γ 0,I u inc .
As a result, we have shown that the second transmission condition is satisfied,

0 = µ γ 1,I,ext w I = µ γ 1,I,int u I -γ 1,I,ext (u + u inc ) . RR n°8988

Operator-valued matrix form

We summarize the results of previous subsection in a more compact form. We write α = D (Dirichlet) , N (Neumann) , I (Impedance), and T (Transmission). Denote by

H s (Γ Obs ) := H s (Γ 1 ) × . . . × H s (Γ N Obs ) . ( 27
)
The single-layer formulation of the boundary problems ( 6) -( 8) can be written as

A α V = F α . (28)
Below, we give the definition and meaning of each object. The unknown V is the product of single-layer densities, appearing in the definition of u scatt ; specifically,

u scatt = N Obs J=1 SJ V J , with V = (V 1 , . . . , V N Obs ) ∈ H 1/2 (Γ Obs ) for α = I, N H -1/2 (Γ Obs ) for α = D, T
.

Next, we give the definition of the matrix A α , an operator-valued matrix

A α =        A α ; 1 A α ; 12 . . . A α ; 1 (N Obs -1) A α ; 1 N Obs A α ; 21 A α ; 2 . . . A α ; 2 (N Obs-1 ) A α ; 2 N Obs . . . . . . . . . . . . . . . A α ; (N Obs -1) 1 A α ; (N Obs -1) 2 . . . A α ; N Obs -1 A α ; (N Obs -1)N Obs A α ; N Obs 1 A α ; N Obs 2 . . . A α ; N Obs (N Obs -1) A α ; N Obs       
.

with mapping property 5 ,

A α : H 1/2 (Γ Obs ) -→ H 1/2 (Γ Obs ) , α = I, N ; A α : H -1/2 (Γ Obs ) -→ H 1/2 (Γ Obs ) , α = D, T .
The operators on the diagonal are given by

A α ; I =          S I,κe , α = D D I,κe -1 2 Id + i λ S I,κe , α = I D I,κe -1 2 Id , α = N µ -1 S I,κ I • D I,κe -1 2 Id - D I,κ I + 1 2 Id • S I,κe , α = T . ( 29
)
The off-diagonal operators are given by

A α ; IJ =          S IJ,κe , α = D D IJ,κe + i λ S IJ,κe , α = I D IJ,κe , α = N µ -1 S I,κ I • D IJ,κe -D I,κ I + 1 2 Id • S IJ,κe , α = T . ( 30
)
Finally, we describe the right hand side (RHS) corresponding to incident wave u inc ,

F α ; I =          γ 0,I u inc , α = D (γ 1,I + i λ γ 0,I ) u inc , α = I γ 1,I u inc , α = N -µ -1 S I,κ I (γ 1,I u inc ) + D I,κ I + 1 2 Id (γ 0,I u inc ) , α = T . ( 31
)

Fredholmness and Invertibility

While the single-layer surface operator S is compact and hence is not invertible, we will show that apart from the set of Dirichlet eigenvalues, which is a discrete and infinite subset of the positive real line, the linear problems obtained from using single-layer Ansatz are invertible. In order to do this, we will need to show that A α -s are Fredholm and explore under which condition they are injective.

Proposition 3. For α = D, I and N , the operator A α is Fredholm. For α = T , if the transmission coefficient µ satisfies µ + 1 = 0, then A T is Fredholm.

Proof. For each type of boundary condition, we will show that A α can be written as a product of an invertible operator and one that is a compact perturbation of the identity map, i.e.,

A α = B α ( Id + K α ) ,
where B α is invertible and K α compact. For this, we will use extensively mapping properties of the surface potentials in the following proof listed in Appendix Section A.

Dirichlet : Decompose A D as

A D = B + C D where B :=       S 1,i 0 . . . 0 0 S 2,i . . . 0 . . . . . . . . . . . . 0 . . . 0 S N Obs ,i       . (32)
We have S I,i :

H -1/2 (Γ I ) → H 1/2 (Γ I ) invertible , with bounded inverse S -1 I,i : H 1/2 (Γ I ) → H -1/2 (Γ I ) bounded . Thus B is invertible with bounded inverse B -1 =       S -1 1,i 0 . . . 0 0 S -1 2,i . . . 0 . . . . . . . . . . . . 0 . . . 0 S -1 N Obs ,i       : H 1/2 (Γ Obs ) -→ H -1/2 (Γ Obs ) . (33)
The components of C are given by

C IJ = S IJ,κe , I = J S I,κe -S I,i , I = J .
We have S IJ,κe :

H -1/2 (Γ J ) → H 1/2 (Γ I ) is bounded and compact S I,κe -S I,i : H -1/2 (Γ J ) → H 1/2 (Γ I ) is bounded and compact . As a result C IJ : H -1/2 (Γ Obs ) -→ H 1/2 (Γ Obs ) is compact and bounded .
This means, for K defined as K = B -1 C, K has following mapping property

K : H -1/2 (Γ Obs ) -→ H -1/2 (Γ Obs ) is compact .
As a result,

(A D ) H -1/2 →H 1/2 = B H -1/2 →H 1/2 Id H -1/2 →H -1/2 + K H -1/2 →H -1/2 (34)
is a Fredholm operator.

Impedance and Neuman : For α = I and N , we can readily decompose

A α = - 1 2 (Id -K α ) (35)
where

K IJ = 2 D IJ,κe + 2 i λS IJ,κe .
We have

S IJ,κe , D IJ,κe : H -1/2 (Γ J ) → H -1/2 (Γ I ) is bounded and compact . As a result K α : H -1/2 (Γ Obs ) -→ H -1/2 (Γ Obs ) is compact.
Transmission : With B defined in (32), we decompose

A T = - µ + 1 2 B + C T .
The component of C T are

(C T ) IJ = µ -1 S I,κ I • D I,κe -1 2 µ -1 (S I,κ I -S I,i ) -D I,κ I • S I,κe -1 2 (S I,κe -S I,i ) , I = J µ -1 S I,κ I • D IJ,κe -D I,κ I + 1 2 Id • S IJ,κe , I = J .
We use the following mapping properties,

S I,κe , S I,κ I : H -1/2 (Γ I ) → H 1/2 (Γ I ) is bounded S I,κe -S I,i , S I,κ I -S I,i : H -1/2 (Γ I ) → H 1/2 (Γ I ) is bounded and compact S IJ,κe , S IJ,κ I : H -1/2 (Γ J ) → H 1/2 (Γ I ) is bounded and compact , I = J D IJ,κe : H s (Γ J ) → H s (Γ I ) is bounded and compact , s = -1 2 , 1 2
to show that

(C T ) IJ : H -1/2 (Γ J ) → H 1/2 (Γ I ) is compact and bounded .
As a result,

C T : H -1/2 (Γ Obs ) -→ H 1/2 (Γ Obs ) is compact and bounded .
Using the mapping property of B -1 , c.f. (33), we obtain

K T := B -1 • C T : H -1/2 (Γ Obs ) -→ H -1/2 (Γ Obs ) is compact and bounded. Inria As a result, for µ + 1 = 0, A T A T = - µ + 1 2 B H -1/2 →H 1/2 Id H -1/2 →H -1/2 - 2 µ + 1 K T ( 36
)
is Fredholm .

By Fredholmness, one has invertibility as soon as injectivity is guaranteed. In the following lemma, we investigate the conditions under which, the operators A α are injective. For the proof, we use the well-posednesses of the exterior boundary values problems, c.f. Section 2.

Lemma 4. If κ 2

e is not a Dirichlet eigenvalue for -∆ for any Ω I with 1 ≤ 1 ≤ N Obs , then A α is injective, with α =D, N, I, and T.

Putting together the above Lemma and Prop. 3, we readily have the following result.

Proposition 5. If κ 2
e is not a Dirichlet eigenvalue for -∆ for any Ω I with 1 ≤ 1 ≤ N Obs , then A α is invertible, with α =D, N and I. For α = T , under the additional assumption that the transmission coefficient µ satisfies µ + 1 = 0, then A T is invertible.

Proof for Lemma 4.

For ψ I ∈ H -1/2 (Γ I ) such that A α ψ = 0 , ψ = (ψ 1 , . . . , ψ N Obs ) ,
we want to show that ψ = 0. This would be the case, if u defined as

u := S1,κe ψ 1 + . . . + S2,κe ψ N ,
is zero in Ω ext . This can be seen as follows. Since SI extends to all R 2 , u as defined solves the Helmholtz equation (-∆ -κ 2 e ) u = 0 both in Ω ext and Ω I , 1 ≤ I ≤ N Obs . Under the assumption that u = 0 in Ω ext , we have

γ 0,I,ext u = 0 , 1 ≤ I ≤ N Obs .
By the continuity of zero-th order jump of the single layers, we thus have γ 0,I,int u = γ 0,I,ext u = 0 .

In another word, u solves the Helmholtz equation at wavenumber κ e with homogeneous Dirichlet boundary value in each domain Ω I . Under the assumption that κ 2 e is not a Dirichlet eigenvalue on Ω I , this means u = 0 in Ω I . As a result,

γ 1,I,int u = 0 , 1 ≤ I ≤ N Obs . Since u = 0 (as currently assumed) in Ω ext , [γ 1,Γ I u] = 0 , I = 1, . . . N Obs .
On the other hand, using the jump identity for single layer, c.f. (68), we have

[γ 1,Γ I SI,κe ψ I ] = ψ I ; [γ 1,Γ I SJ,κe ψ J ] = 0 , I = J . ⇒ [γ 1,Γ I u] = [γ 1,Γ I S 1 ψ 1 ] + . . . + [γ 1,Γ I S 2 ψ N ] = ψ I .

RR n°8988

As a result, under the assumption that u = 0 in Ω ext , we have shown the jump

ψ I = [γ 1,Γ I u] = 0 ⇒ ψ = 0 .
To finish the proof of the lemma, it remains to prove that under non-Dirichlet eigenvalue assumption on κ Transmission : For the interior domain Ω I , we consider u I defined as

u I := µ -1 SI,κ I • γ 1,I,ext u -DI,κ I • γ 0,I,ext u .
By Prop. 2, that ψ is in the kernel of A T means that {u, u 1 , . . . , u N Obs } solves the transmission problem ( 5) with u inc = 0. The well-posedness of the transmission for Helmholtz gives that u = 0 in Ω ext .

As a result of the above Proposition, as long as we stay away from the set of Dirichlet eigenvalues (DEV), which is a discrete and infinite subset of the positive real line, the linear problem obtained from using single-layer Ansatz is invertible. For our application, we focus on small obstacles, and thus it would be useful to know numerically how small the obstacles should be in order to be in the region of invertibility. This is determined from a positive lower bound of the first eigenvalue. The following remark cites a lower bound for a general domain.

Remark 2.

1. We cite the Isoperimetric inequality (Rayleigh-Faber-Kahn) for Dirichlet eigenvalues of the Laplacian

λ 1 (Ω) = inf u∈H 1 0 (Ω)\{0} Ω |∇u| 2 dx Ω u 2 dx ≥ π |Ω| j 2 0,1
, where j 0,1 first zeros of the Bessel function J 0 , and |Ω| denotes area of the membrane. The equality is attained if and only if the membrane is circular.

2. If κ satisfies κ r circumscribed circle of Ω < 2 , ( 37
)
then κ 2 is not a Dirichlet eigenvalue on Ω. This is due to j 0,1 ∼ 2.40 > 2, and

κ 2 < π 2 2 π r 2 circumvent circle < π 2 2 |Ω| < π |Ω| j 2 0,1 < λ 1 (Ω) .

Inria

For a disc domain, we have explicit expressions for all of the DEV-s.

Remark 3. The Dirichlet eigenvalues on a disk of radius R are

λ n,m = j n,m R 2 , m ≥ 0 , m ≥ 1
where j n,m is the m-th positive root of J n (r) = 0. We have λ n,m has multiplicity 2 when n ≥ 1.

The corresponding eigenfunctions

1. To λ 0,m m ≥ 1, dimension eigenspace is 1 generated by J 0 j 0,m R r . 2. To λ n,m n, m ≥ 1, dimension eigenspace is 2 generated by J n j n,m R r cos(nθ) , J n j n,m R r sin(nθ) .
3. The first 4 roots : j 0,1 ∼ 2.40 ; j 1,1 ∼ 3.83 ; j 2,1 ∼ 5.13 ; j 1,2 ∼ 5.52 .

Variational forms

For α = D and T : From the proof of Proposition 3, we have the decomposition,

A α = Ãα + K α ,
where Ãα is coercive with respect to H -1/2 (Γ Obs ), and

Ãα : H -1/2 (Γ Obs ) -→ H 1/2 (Γ Obs
) is bounded and invertible .

On the other hand,

Kα : H -1/2 (Γ Obs ) -→ H -1/2 (Γ Obs ) is bounded and compact .
Each of the operators listed above corresponds to a sesquilinear form, denoted respectively by a α , ãα , and k α , with

a α , ãα , k α : H -1/2 (Γ Obs ) × H -1/2 (Γ Obs ) -→ C a α (ψ, φ) := A α ψ , φ H 1/2 , H -1/2 ãα (ψ, φ) := Ãα ψ , φ H 1/2 , H -1/2 k α (ψ, φ) := K α ψ , φ H 1/2 , H -1/2
Here, ψ = (ψ

I ) 1≤J≤N Obs , φ = (φ J ) 1≤J≤N Obs ∈ H -1/2 (Γ Obs ) .
We also have a α = ãα + k α .

For α = I and N : From the proof of Proposition 3, we have the decomposition,

A α = 1 2 Id + K α ,
There correspond sesquilinear forms, a α , k α , with a α , k α :

H 1/2 (Γ Obs ) × H -1/2 (Γ Obs ) -→ C a α (ψ, φ) := 1 2 ψ , φ H 1/2 , H -1/2 + k α (ψ, φ) k α (ψ, φ) := K α ψ , φ H 1/2 , H -1/2 .
Also define the linear function corresponds to the RHS F α ,

α : H -1/2 (Γ Obs ) -→ C α (ϕ) = F α , ϕ H 1/2 ,H -1/2 = N Obs I=1 F α,I , ϕ I H 1/2 (Γ I ),H -1/2 (Γ I ) . ( 38
)
Variational form : The variational form of the problems in (28) are written as :

For F α ∈ H 1/2 (Γ Obs ) , find ψ α ∈ H -1/2 (Γ Obs ) for α = D, T H 1/2 (Γ Obs ) for α = I, N such that a α (ψ, ϕ) = α (ϕ) , ∀ ϕ ∈ H -1/2 (Γ Obs ).
(39)

General comments on discretization and error analysis

The discretization of the variational problem (39) is via a dense sequence of finite-dimensional spaces in

H -1/2 (Γ Obs ) (for α = D, T ) or H 1/2 (Γ Obs ) (for α = I, N ). Let {V l } l∈N be a dense sequence of finite-dimensional subspaces in H -1/2 (Γ) for α = D, T , and H 1/2 (Γ) for α = I, N . Define V m := V m (Γ 1 ) × . . . × V m (Γ N Obs ) .
Discretized problem : With α defined in (38), the discretized version of the problem (28) for approximation of order m is written as :

For F α ∈ H 1/2 (Γ Obs ) , find ψ h ∈ V m such that a α (ψ h , ϕ) = α (ϕ) , ∀ ϕ ∈ V m . (40)
Well-posedness of the discrete problem and Error Analysis : In the previous section, we have obtained a decomposition of a α to show that it satisfies the hypothesis of [START_REF] Sauter | Boundary Element Methods[END_REF][Thm 4.2.9,p.229]. As a result, we have the well-posedness of the approximate problem, when injectivity requirement is satisfied, i.e κ 2 e is not an interior Dirichlet eigenvalue for Ω I , 1 ≤ I ≤ N Obs , c.f. Lemma 4. Note that we already established this for the continuous problem, under this condition, c.f. Prop 5. In addition, we have a Céa-type estimate for the approximation error. For convenience, we cite Thm 4.2.9 from [START_REF] Sauter | Boundary Element Methods[END_REF].

Theorem 6 ([21]

). Let H be a Hilbert space and {V k } k a dense sequence of finite-dimensional subspaces in H. We assume 1. sesquilinear form a is elliptic :

|a(u, u)| ≥ C u 2
H with α > 0 ; 2. the operator T ∈ L(H, H ) associated with sesquilinear form t is compact ; 3. we have injectivity : a(u, v) + t(u, v) = 0, ∀v ∈ H \ {0} ⇒ u = 0 ; Inria Then there exists a constant k 0 > 0 such that for all k ≥ k 0 , the Galerkin equations

For F ∈ H , find u k ∈ V k such that a(u k , v k ) + t(u k , v k ) = F (v k ) , ∀ v k ∈ V k have unique solution u k in V k .
Moreover, u k converges to u, and , for k ≥ k 0 , we have the quasi-optimal error estimate

u -u k H ≤ C inf v k ∈V k u -v k H
Céa estimate .

Fourier Series Galerkin Basis

We introduce the basis coming from the Fourier series on the closed curve Γ I . Denote by T m the set of trigonometric polynomials on [0, 2π],

T m ([0, 2π]) = m k=-m a k e ikθ , a k ∈ C . This set is dense in H r (0, 2π) for all r ∈ R. Fix a C 2 parametrization of Γ J , γ J : [0, 2π] -→ Γ J ⊂ R 2 .
When Ω J is disc-shaped, we work with the following parametrization,

θ -→ x = x J + r J (cos θ , sin θ) . (41)
The basis on Γ I is defined such that, when pulled back to [0, 2π] via γ Γ J gives ∪ m ∈N T m . More specifically, we define the basis functions w I,m on Γ I , by

γ * I w I,m = w J,k • γ Γ J = e -i k θ . ( 42
)
For this type of basis, the approximation is a special case of a projection method, thus general error analysis can also be obtained from the theory of projection theory, c.f. [START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF][Section 5]. In particular, as m → ∞, the discrete solution u scatt,h in V m is of the form

u scatt,h = N Obs J=1 SJ m k=-m v J,k w J,k ,
and converges in the mean square to

u scatt = N Obs J=1 SJ ∞ k=-∞ v J,k w J,k .
Substituting the above form of u scatt,h in (40), we obtain the linear system in the unknowns {v J,m } have to satisfy,

N Obs J=1 m k=-m v J,k a α SJ w J,k , w I,l = α (w I,l ) , ∀ 1 ≤ I ≤ N Obs , l ∈ Z . ( 43
)
with α defined in (38), corresponding to F α ∈ H 1/2 (Γ Obs ), which is defined in terms of u inc , c.f. (31). In fact, the above expression gives an infinite matrix representation of A α,IJ , with components given by

(A α,IJ ) kl := a α SJ w J,k , w I,l , k, l ∈ Z . ( 44
)
Note that the RHS is a double line integral on Γ J × Γ I . For example, when α = D, this double line integral defining (A D,IJ ) kl is

(A D,IJ ) kl = Γ I w Γ I ;k Γ J G κe (x, y) w Γj ;l ds(Γ J ) ds(Γ I ) = 2π 0 2π 0 e ikθ G κe γ I (θ), γ J ( θ) e il θ |γ I (θ)| γ J ( θ) d θ dθ . ( 45
)
Infinite matrix+vector form : We can further write (28) in infinite matrix + vector form

A α V = F α , V = (V J,l ) 1≤J≤N Obs , l∈Z , (46)
where A α is composed of N Obs × N Obs block matrix, each of which is an infinite matrix, with components defined by (44), and the RHS

F α = (F 1 , . . . , F N Obs ) ; F Il = α (w I,l ) , l ∈ Z , 1 ≤ I ≤ N Obs .
Here we use the same notation for the operator A IJ and its infinite matrix representation, similarly for F α . The discretized/approximating problem (43) at order m is a linear system of size 2m + 1, and is written as

A α,h V h = F α,h , V h = (V J,l ) 1≤J≤N Obs , -m≤l≤m , (47)
where A α,h is a truncated version of A α , i.e composed of N Obs × N Obs block matrix, each of which is a matrix of size (2m + 1) × (2m + 1), with components defined by (44), i.e

(A α,h,IJ ) lk := a α ( SJ w J,k , w I,l ) , -m ≤ k, l ≤ m , 1 ≤ I, J ≤ N Obs . Similarly, F α,h is a truncated version of F α F α,h,Il = α (w I,l ) , -m ≤ l ≤ m , 1 ≤ I ≤ N Obs .
We call this approximating method Fourier Series -Single Layer (FS-SL).

Multiple scattering for Disc-shaped obstacles

Notations on the formation of obstacles : We will consider a formation of N Obs nonoverlapping disc-shaped obstacles, each represented by the disc B(x I , r I ) centered at x I ∈ R 2 of radius r I . For 1 ≤ I ≤ N Obs , denote the polar coordinates relative to x I ∈ R 2 by (r I (•), θ I (•)). Specifically,

r I (x) := |x -x I | ; θ I (x) ∈ [0, 2π) with cos θ I (x), sin θ I (x) = x -x I |x -x I | ; Inria and thus x = x I + r I (x) cos θ I (x) , sin θ I (x) . ( 48
)
Denote by d IJ the distance between the centers of obstacle I and J, and by θ IJ and θ JI their relative polar coordinates. In the above notation (48), their definitions are given by

d IJ = |x I -x J | ; θ IJ = θ I (x J ) ; θ JI = θ J (x I ) ; with x I = x J + d IJ (cos θ JI , sin θ JI ) ; x J = x I + d IJ (cos θ IJ , sin θ IJ ) .
The non-overlapping assumption is given by

d IJ > r I + r J . (49)

Linear systems

We recall the total wave in the exterior of the obstacles, denoted by u total , is a superposition of the unknown scattered wave u and the incidence wave u inc . In addition, we have chosen the single layer Ansatz for u, and decompose it further as linear combination of single-layer with Fourier series basis w J,l as potential,

u total = u + u inc ; u inc = N Obs J=1 m∈Z v J,m SJ,κe w J,m .
Thus, our unknown is the infinite scalar vector v J,m indexed by the degree of the Fourier node m ∈ Z, and the index of the obstacle J with 1 ≤ J ≤ N Obs .

In order to arrive at the linear equations satisfied by v J,m , we will need the multipole expansion of the incident plane wave and the single layer SJ,κe w J,l . We rewrite these results in the notation (48) from Appendix B.3 and B.4.

• For I = J : the single layer for points x such that r I (x) < d IJ , in terms of the relative polar coordinates with respect to x I , is given by

SJ w J,l (rI (x), θI (x)) = iπ rJ 2 J l (κe rJ ) ∞ m=-∞ H (1)
l-m (κe dIJ ) e i(l-m)θ JI Jm(κe rI (x)) e imθ I (x) .

(50)

Taking the exterior zero-th and first order traces along Γ I ,

γ 0,I,e SJ w J,l (θ I (x)) = iπ r J 2 J l (κ e r J ) ∞ m=-∞ H (1) l-m (κ e d IJ ) e i(l-m)θ JI J m (κ e r I ) e imθ I (x) ; γ 1,I,e SJ w J,l (θ I (x)) = iπ r J κ e 2 J l (κ e r J ) ∞ m=-∞ H (1) l-m (κ e d IJ ) e i(l-m)θ JI J m (κ e r I ) e imθ I (x) . (51)
• For I = J, the single layer for points x such that x ∈ Ω i,I , in terms of the relative polar coordinates with respect to x I , is given by,

SJ w J,l (r J (x), θ J (x)) = iπ r J 2 e i l θ J (x) J l (κ e r J ) H (1)
l (κ e r J (x)) .

RR n°8988

Taking the exterior zero-th and first order traces along Γ I , γ 0,I,e SJ w J,l (θ J (x)) = i π r J 2 e i l θ J (x) J l (κ e r J ) H

l (κ e r J ) ; γ 1,I,e SJ w J,l (θ J (x)) = i π r J κ e 2 e i l θ J (x) J l (κ e r J ) H

(1) l

(κ e r J ) .

(52)

We cite the Jacobi-Anger expansion for the (time-harmonic) acoustic plane wave incident at angle α inc u inc (rθ) = e i κx •(cos αinc,sin αinc) = e i κ r cos(θ-αinc) ,

given by, c.f. Appendix B.3

u inc (x) = u inc (x I ) ∞ l=-∞ i l J l (κ e r I (x)) e i l (θ I (x)-αinc) . ( 54
)
Taking the exterior zero-th and first order traces along Γ I ,

(γ 0,I,e u inc ) (θ I (x)) = u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) ; (γ 1,I,e u inc ) (θ I (x)) = κ e u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) . (55)
Linear systems : We can now obtain the linear equations satisfied by unknowns {v I,l } 1≤I≤N Obs , l∈Z .

This can be obtained either by taking variational pairing of the boundary conditions against test functions w J,l -s, or by simply matching the coefficient of the Fourier series. Both processes in essence are the same, and use the orthogonality of the basis {w I,l } l∈Z with respect to L 2 (Γ I ); however for a compact exposition, we choose the second one.

For Dirichlet problem : The coefficient of e ilθ I (x) in Dirichlet condition γ 0,I,ext u total = 0 along Γ I gives,

u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) + ∞ l=-∞ v I,l iπ r I 2 e i l θ I (x) J l (κ e r I ) H (1)
l (κ e r I ) = 0 .

(

) 56
For Impedance problem : The coefficient of e ilθ I (x) in the Impedance condition

( i λ γ 0,I,ext + γ 1,I,ext ) u total = 0 along Γ I gives, iλ u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) + i λ ∞ l=-∞ v I,l iπ r I 2 e i l θ I (x) J l (κ e r I ) H (1)
l (κ e r I ) +

Inria + i λ ∞ l=-∞ N Obs J=1 J =I v J,l iπ r J 2 J l (κ e r J ) ∞ m=-∞ H (1)
l-m (κ e d IJ ) e i(l-m)θ JI J m (κ e r I ) e imθ I (x)

+ κ e u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) + iπ r I κ e 2 ∞ l=-∞ v I,l e i l θ I (x) J l (κ e r I ) H (1) l
(κ e r I )

+ iπ κ e 2 ∞ l=-∞ N Obs J=1 J =I v J,l r J J l (κ e r J ) J l (κ e r J ) ∞ m=-∞ H (1)
l-m (κ e d IJ ) e i(l-m)θ JI J m (κ e r I ) e imθ I (x) = 0 .

We arrive at the simplified equation satisfied by v I,l for 1 ≤ I ≤ N Obs and l ∈ Z,

i λ v I,l iπ r I 2 J l (κ e r I ) H (1)
l (κ e r I ) +

- π λ 2 N Obs J=1 J =I ∞ m=-∞ v J,m r J J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I )

+ iπ r I κ e 2 v I,l J l (κ e r I ) H (1) l
(κ e r I )

+ iπ κ e 2 N Obs J=1 J =I ∞ m=-∞ v J,m r J J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I ) = -i λ u inc (x I )i l J l (κ e r I ) e -i l αinc -κ e u inc (x I ) i l J l (κ e r I ) e -i l αinc .
This further simplifies to

v I,l π r I J l (κ e r I ) 2 -λ H (1)
l (κ e r I ) + i κ e H

(1) l

(κ e r I )

+ π 2 N Obs J=1 J =I -λ J l (κ e r I ) + i κ e J l (κ e r I ) ∞ m=-∞ v J,m r J J m (κ e r J ) H (1)
m-l (κ e d IJ ) e i(m-l)θ JI = -u inc (x I ) i l e -i l αinc i λ J l (κ e r I ) + κ e J l (κ e r I ) .

(57)

For Neumann problem : The linear equations for v J,m is obtained by setting λ = 0 in (57),

-u inc (x I ) i l e -i l αinc κ e J l (κ e r I ) = v I,l π r I J l (κ e r I ) 2 i κ e H (1) l
(κ e r I )

+ π 2 N Obs J=1 J =I i κ e J l (κ e r I ) ∞ m=-∞ v J,m r J J m (κ e r J ) H (1)
m-l (κ e d IJ ) e i(m-l)θ JI .

(

) 58
For Transmission problem : Instead of using linear system (24), i.e using Green's representation for the interior wave inside of each obstacle, we will re-derive the linear system equation, taking advantage of the disc geometry, which gives a solution of the Helmholtz equation inside a disc. Explicitly, the interior wave u i inside Ω I has the form,

u i (r I (x), θ I (x)) = ∞ l=-∞ c I,l J l (κ I r I (x)) e i l θ I (x) , x ∈ Ω I .
Taking the interior zero-th and first order trace along Γ I gives,

(γ 0,I,i u i ) (θ I (x)) = ∞ l=-∞ c I,l J l (κ I r I ) e i l θ I (x) , x ∈ Γ I ; (γ 1,I,i u i ) (θ I (x)) = κ I ∞ l=-∞ c I,l J l (κ I r I ) e i l θ I (x) , x ∈ Γ I . ( 59
)
For the moment, we have the second set of unknowns c I,l associated with the interior wave; however, after some elimination, we will end up with a linear system only in terms of v J,m . We remark that this approach does not give a different linear system, but simply the (simplified) form of (24) in the setting of disc geometry.

The first transmission condition γ 0,I,e u total = γ 0,I,i u i along Γ I gives

u inc (x I ) ∞ l=-∞ i l J l (κ e r I ) e i l (θ I (x)-αinc) + ∞ l=-∞ v I,l iπ r I 2 e i l θ I (x) J l (κ e r I ) H (1) l (κ e r I ) + + ∞ l=-∞ N Obs J=1 J =I v J,l iπ r J 2 J l (κ e r J ) ∞ m=-∞ H (1)
l-m (κ e d IJ ) e i(l-m)θ JI J m (κ e r I ) e imθ I (x) = ∞ l=-∞ c I,l J l (κ I r I ) e i l θ I (x) , x ∈ Γ I .
Matching the coefficient of e ilθ I (x) we obtain the first linear system for v J,m and c J,m , c I,l J l (κ I r I ) = u inc (x I ) i l J l (κ e r I ) e -i l αinc + v I,l iπ r I 2 J l (κ e r I ) H

l (κ e r I )

+ N Obs J=1 J =I ∞ m=-∞ v J,m iπ r J 2 J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I ).

(

) 60
The coefficient of e ilθ I (x) in the second transmission condition γ 1,I,e u total = µ γ 1,I,i u i along Γ I gives

µ κ I κ e c I,l J l (κ I r I ) = u inc (x I ) i l J l (κ e r I ) e -i l αinc + v I,l iπ r I 2 J l (κ e r I ) H (1) l
(κ e r I )

+ N Obs J=1 J =I ∞ m=-∞ v J,m iπ r J 2 J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I ) .

(

) 61
Next we will eliminate c I,l -s from ( 60) and (61). Multiply both sides of (60) by µ κ I κ e J l (κ I r I ).

and of (61) by J l (κ I r I ), then subtract the resulting linear equations, we then obtain ones involving only the variables v J,l -s,

µ κ I κ e J l (κ I r I ) u inc (x I ) i l J l (κ e r I ) e -i l αinc + v I,l iπ r I 2 J l (κ e r I ) H (1)
l (κ e r I )

Inria + N Obs J=1 J =I ∞ m=-∞ v J,m iπ r J 2 J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I )

= J l (κ I r I ) u inc (x I ) i l J l (κ e r I ) e -i l αinc + v I,l iπ r I 2 J l (κ e r I ) H (1) l
(κ e r I )

+ N Obs J=1 J =I ∞ m=-∞ v J,m iπ r J 2 J m (κ e r J )H (1)
m-l (κ e d IJ ) e i(m-l)θ JI J l (κ e r I ) .

This simplifies to

u inc (x I ) i l e -i l αinc -µ κ I κ e J l (κ e r I ) J l (κ I r I ) + J l (κ e r I ) J l (κ I r I ) = v I,l iπ r I 2 J l (κ e r I ) µ κ I κ e H (1)
l (κ e r I ) J l (κ I r I ) -H

(κ e r I ) J l (κ I r I )

+ N Obs J=1 J =I iπ r J 2 ∞ m=-∞ v J,m J m (κ e r J ) H (1)
m-l (κ e d IJ ) e i(m-l)θ JI µ κ I κ e J l (κ e r I ) J l (κ I r I ) -J l (κ e r I ) J l (κ I r I ) .

Remark 4. The well-posedness of the transmission problem is not affected, even when κ I is an interior Dirichlet or Neumann eigenvalue of Ω I . The eigenvalue condition means that, for some l ∈ Z, we have J l (κ I r I ) = 0 for Dirichlet, and J l (κ I r I ) = 0 for Neumann. In fact, in these cases, the transmission problem can be written in terms of a mixed boundary values problem as follows.

First, we note that a wavenumber cannot both be a Dirichlet and Neumann eigenvalues. In the case of disc geometry, this can also be seen from the fact that the zeros of J l and J For I ∈ N D ∪ N N , let n(I) be such that such that

J n(I) (κ I r I ) = 0 or J n(I) (κ I r I ) = 0 .
The solution to the transmission problem {u total , u 1 , . . . , u N Obs } is given by the unique solution to the mixed problem as follows :

• First, we solve for {u total , (u I ) I∈N T } in the mixed boundary value problem,

                   (-∆ -κ 2 e ) u total = 0 x ∈ Ω ext (-∆ -κ 2 I ) u I = 0 x ∈ Ω I , I ∈ N T γ 0,I,e u total = 0 x ∈ Γ I , I ∈ N D γ 1,I,e u total = 0 x ∈ Γ I , I ∈ N N γ 0,I,e u total -γ 0,I,i u I = 0 ; γ 1,I,e u total -µγ 1,I,i u I = 0 x ∈ Γ I , I ∈ N T lim r→∞ √ r (∂ r -i κ e ) (u total -u inc ) = 0 ; r = |x| . • Then, for I ∈ N D ∪ N N , u I is obtained by u I = c I J n(I) (κ I r I (x)) e in(I) θ I (x) .
For I ∈ N D , the factor c I is uniquely determined by (61), and for I ∈ N N , by (60).

Matrix forms of the continuous problem and the discrete one

We summarize the results of subsection 4.1 in matrix form as follows. Note that the following formulas are the simplification of the general ones ( 46) and (47) in the disc geometry. For α = D, N, I and T ,

Continuous problem A α V = F α ; Approximating problem A α,h V h = F α,h . with V = (V J,l ) 1≤J≤N Obs , l∈Z ; V h = (V J,l ) 1≤J≤N Obs , -m≤l≤m . ( 62
)
Here m is the order of approximation. Below, we give the definition and meaning of each object. The coefficient matrix A α is described in terms of its diagonal and off-diagonal blocks,

A α =        A α ; 1 A α ; 12 . . . A α ; 1 (N Obs -1) A α ; 1 N Obs A α ; 21 A α ; 2 . . . A α ; 2 (N Obs-1 ) A α ; 2 N Obs . . . . . . . . . . . . . . . A α ; (N Obs -1) 1 A α ; (N Obs -1) 2 . . . A α ; N Obs -1 A α ; (N Obs -1)N Obs A α ; N Obs 1 A α ; N Obs 2 . . . A α ; N Obs (N Obs -1) A α ; N Obs       
.

We recall from Subsection 3.6 that A α,h is a truncated version of A α , i.e composed of N Obs ×N Obs block matrix, each of which is a matrix of size (2m + 1) × (2m + 1). The diagonal blocks A α,I are diagonal infinite matrices, with diagonal components given by

(A α,I ) ll = i π r I J l (κ e r I ) 2 ×            H (1)
l (κ e r I ) , α = D κ e H (1) l (κ e r I )
, α = N i λ H

l (κ e r I ) J l (κ I r I ) -H

(κ e r I ) J l (κ I r I ) , α = T .

(63) For I = J, the components of the off-diagonal block A α,IJ are given by,

(A α,IJ ) lm = i π r J 2 J m (κ e r J ) H (1) m-l (κ e d IJ ) e i(m-l)θ JI ×            J l (κ e r I ) , α = D κ e J l (κ e r I ) , α = N i λ J l (κ e r I ) + κ e J l (κ e r I ) , α = I µ κ I κ e J l (κ e r I ) J l (κ I r I ) -J l (κ e r I ) J l (κ I r I ) , α = T . (64)
They are the simplification of the general ones (44) for disc geometry.

For the scattering of the plave wave (54), the components of the RHS are given by

F α;I,l = -u inc (x I ) i l e -i l αinc ×            J l (κ e r I ) , α = D κ e J l (κ e r I )
, α = N i λ J l (κ e r I ) + κ e J l (κ e r I )

, α = I µ κ I κ e J l (κ e r I ) J l (κ I r I ) -J l (κ e r I ) J l (κ I r I ) , α = T .

(65) After solving linear system (62) to obtain V h , we obtain the numerical solution u scatt , expressed below in multipole expansion, c.f. Appendix B.4.

Inria

Proposition 7. The approximate solution given by method FS-SL at approximation order m is given by u total,h (x) = u planewave (x) + u h,scatt (x) ,

with u h,scatt (x) := N J=1 m l=-m v J,l SΓ J w J,l ;
where {v J,l } solves the linear system (62), and

SΓ J w J,l (r J (x), θ J (x)) = iπ r J 2 e i l θ J (x) J l (κ r J ) H (1)
l (κ r J (x)) .

Remark 5. The above result can be rewritten in the language of Multipole Theory; in the end, what is at work is the separable geometry which enables the use of separation of variables to solve the PDE. To describe a wave scattered by an obstacle, we use as basis functions the multipoles ψ n -s defined as ψ n (r, θ) := H (1) n (κ r) e inθ , n ∈ Z . Here r and θ are the relative polar coordinates with respect to the obstacle. In another word, the wave emitted by obstacle J, denoted by u scatt,J , is a superposition of multipoles of all orders placed at the center of the scatterer,

u scatt,J (x) = k∈Z c Jk ψ k (r J (x), θ J (x)) .
On the other hand, to describe a wave incident on an obstacle, we use the basis functions ψ n -s, defined as ψ n (r, θ) := J n (κ r) e inθ , n ∈ Z .

We next use the boundary conditions which give the interaction between an obstacle and a waving acting on it: e.g.. in the case of soft-scattering with the incident planewave u pw (54),

k∈Z c Ik , ψ k (r I , θ I ) + J =I 1≤J≤N Obs c J;k ψ k (r J , θ J ) = -u pw (x) , x ∈ Γ I .
As before, we use Graf 's addition theorem Appendix B.2, to expand ψ k (r J , θ J ) in the coordinates (r I , θ I ), and the Jacobi-Anger expansion for plane wave, c.f. Appendix B.3, in the basis ψ n (r I , θ I ). Finally, to obtain an infinite system for unknown c I,k , we either use the orthogonality of {e i k θ } or just simply match the coefficients.

It should be noted that although they give the same solutions, FS-SL method gives more explicit description of the coefficients, in the sense that the unknowns of FS-SL method are {v J,k } while those in Multipole method are {c J,k } with c J,k playing the role of 1 2 v J,k iπr J J k (κr J ).

Numerical Convergence

General theoretical convergence is guaranteed either by Céa Lemma, c.f. the discussion in Section 3.5, or by projection operator theory, c.f. Subsection 3.6, when the wavenumber κ 2 e is away from the set of Interior Dirichlet Eigenvalues. We validate this fact with the following experiments, showing the numerical convergence for a formation of 10 arbitrarily spaced obstacles. We observe that the rate of convergence improves as the order of the approximation increases and the size of the obstacles decreases. Remark 6. For the disc geometry, the rate of convergence can be described explicitly in the form of an upper bound for the approximation error, in terms of the size of the obstacles, distance between the obstacles, and number of the obstacles. The technical ingredient needed for this is an upper and lower bound in terms of these parameters for Bessel functions J n , Y n and H n . For example, when the obstacles are small enough, one can use the fact that J n (t) is an alternating series for 0 ≤ |t| < 2, and thus can be bounded from below and above by appropriate truncated sums from its defining series. The final bound is then obtained by using Hilbert-Schmidt theory. For a fixed wavelength and small enough obstacles, it should be expected that the rate of convergence increases as the size of the obstacles decreases, as the distance between the obstacles decreases, and decreases as the number of the obstacles increases.

We recall the scattered field is described as a sum of single layer potentials

u scatt = N Obs J=1 SJ V h,J ; V J single layer density . ( 66
)
The norms used in the comparison are the Sobolev norm H s , c.f. ( 27), of the single-layer densities {V J,h }, and the L ∞ of u scatt on the boundary of the domain in interest.

For the results in Figure 2 (in Sobolev norm) and 4 (in L ∞ ), for each stated ratio Obstacle Radius Wavelength , we investigate the rate of convergence with respect to the order of approximation. To obtain the data, we calculate • V ref has Order of Approximation = 16 ;

• then V test is calculated for Order of Approx ranging in 1 : 10 ;

• and obtain the datum : the error V h,test -V h,ref . This error in Sobolev norm H s is calculated as,

V h,test -V h,ref Hs = N Obs J=1 m ref k=-mtarget |v h,test;J,k -v h,ref;J,k | 2 (1 + k 2 ) s . ( 67
)
For the results in Figure 3 (in Sobolev norm) and 5, for each stated order of approximation OrdApp test , we investigate the rate of convergence with respect to the size of the obstacles. Specifically, we obtain results for following grid size. • and obtain the datum : the error V h,test -V h,ref .

At each ratio

The following figure shows that, for a fixed precision tolerance, smaller obstacles need lower approximation order. In particular, for precision 10 -7 , with the range of obstacle shown, the approximating order needed decreases from 8,7,4 to 3. The following figure shows that rate of convergence increases with the order of approximation. More specifically, in the asymptotic region, the error ∼ (Obs Radius) slope of curve , and the slope increases with the order of approximation.

Radius Obstacle Wavelength

Error in H 1/2

-norm

Figure 3: Convergence curve with respect to size of obstacles (10 randomly distributed obstacles). OrderApproximation = 2 ; = 3 ; = 4 ; = 6

The following figures show the same analysis but in the L ∞ norm of the scattered field on the boundary of a rectangle away from the obstacles. The rate of convergence is much bigger than that corresponding to the Sobolev norm of the density; in other words, the second norm is more 'pessimistic'. For example, for precision 10 -7 , the largest obstacles require approximating order 4, and the smallest one 2, which are much lower that required when one uses Sobolev norm to measure error, as shown in Figure 2. In the following figure, the convergence is so strong that, for the approximation order 6, we are mostly in the precision machine region 10 -17 . In working with 'far-field' wave, for small enough obstacles, one can use approximating order as low as 2.

Radius Obstacle Wavelength

Error in L ∞ -norm on boundary In this section, we compare the performance of the optimized software Montjoie, as a representative of the FEM family, with that of our method (FS-SL) using solver Mumps. Being mesh-free, a IE-based method generates a smaller linear system than a FE-based method, especially on large domain with obstacles that are small compared to the wavelength of the incident wave. On the other hand, the post-processing time for the FS-SL method, which involves evaluation of Hankel functions on a structured grid, can be expensive. With that said, even with exact evaluation, the FS-SL method still takes less time than Montjoie. The drawback of single-layer method can be diminuished, with interpolation, for e.g Hermite cubic spline, c.f. [START_REF] Hairer | Lecture notes introduction à l'analyse numérique[END_REF][p.48-50]. In this case, we see a great improvement for FS-SL:

• for 6 obstacles, FS-SL + Mumps takes 0.06 secs, compared to 14 secs for Montjoie, for a precision of 1.e-64, c.f. subsection 5.1;

• and for 200 obstacles, SL + Mumps takes 4.8 secs, compared to 33.8 secs for Montjoie, for a precision of 1.e-3, c.f. subsection 5.2.

Last but not least, one can switch more readily to higher precision tolerance with FS-SL, and hence obtain the convergence curve more easily. This is because the size of linear system for FS-SL method is N Obs × (2 × Order of Approx + 1) .

For Montjoie, the size of the linear system (already large) increases sharply, when one goes from Qn to Qn+1 or refine the mesh. In addition, the precision of Montjoe is also constrained by the order of PML. For 200 obstacles, it takes a lot of memory and time to obtain the convergence curve for Montjoie; in fact we stopped at order 16 at error 1.e-7 for this case, c.f. Figure 11.

Test 6 holes

Parameters for Montjoie

In geo file, lc = 1.0 ; AddPML = YES PML XY 3.0 AUTO ; DampingPML = 4.0

Reference solutions and numerical convergence : The numerical convergence for both methods are shown in Figure 7. We choose MJ Ref (Q15) and FS-SL Ref [START_REF] Hettlich | Frechet derivatives in inverse obstacle scattering[END_REF] as reference solutions (at precision tolerance 1.e-12), see also Figure 8.

Inria

Reference solutions and numerical convergence : The numerical convergence for both methods are shown in Figure 11. For 200 obstacles, it takes a lot of memory and time to obtain the convergence curve for Montjoie (in fact we stopped at order 16 at error 1.e-7). We choose MJ Ref (Q17) and FS-SL Ref ( 14) as reference solutions (at precision tolerance 1.e-8), see also Figure 12. In this section, we restrict ourselves to the FS-SL method, and do a pre-processing time-cost comparison among the solvers, which include direct solvers (Mumps, Lapack and Scalapack) and GMRES iterative solvers (with various choice of preconditioners). We refer the readers to Appendix C, for the definition of the pre-conditioners in the following investigation. The GMRES iterative solvers are obtained by modifying the codes given by [START_REF] Frayssé | A set of gmres routines for real and complex arithmetics on high performance computers[END_REF], to which we added the codes for preconditioners. The 'quality' of the solution is based on the relative difference in product Sobolev norm H 1/2 (Γ Obs ) of SL densities {v h,I } 1≤I≤N Obs , c.f. (67). Since this criteria can be quite pessimistic, for certain tests, we also use the L 2 norm (of the domain of visualization) of the final scattered field. To reduce the post-processing time, we also use interpolation with Hermite cubic spline, c.f. [11][p.48-50].

Observations:

Due to the denseness of the coefficient matrix A α , Mumps, even when run in parallel, loses its efficiency compared to Lapack and Scalapack. This behavior becomes more pronounced, as the number of the obstacles increases, c.f. subsection 6.4. This is an instrinsic property of general integral equation (IE) methods. For closed together obstacles, the direct solvers outperform the GMRES family. The second one has difficulty in converging; even in the case of convergence, they give 'lower quality' solutions. For very large number of obstacles (around 10 4 ), the only solver that can still run and give result is Scalapack, the pre-processing time of which is 22 mins 20 secs on 48 cores, c.f. Subsection 6.6.

Observation among the GMRES solvers :

1. The preconditioners derived from Gauss-Seidel outperform those from Jacobi or no preconditioning, with the best ones being the family of SGS and LUSGS. For 200 obstaces, there is a big difference in the number of iterations required, with the SGS and LUSGS preconditioners requiring 70 iterations, compared with 650 for Jacobi.

2. For obstacles that are spaced closely, convergence is not obtained except for family LUSGS and SGS, c.f. subsection 6.2 (for 1616 obstacles). For this last family, the number of iterations increases drastically, from 70 (for 200 obstacles) to close to 900 (for 1616 obstacles), c.f. subsections 6.1, 6.2, 6.4. Lowering the (GMRES) precision tolerance has little effect on the number of iterations required, c.f. Subsection 6.3.

3. However, when the obstacles are further apart, we observe a 'come-back' of iterative solvers. For 2000 obstacles, for a precision at 1.e-7 (in L 2 of domain of visualization), SGS and LUSGS only need around 70 iterations. This places them ahead of Lapack, and head-tohead with Scalapack, c.f. Subsection 6.5.

Conclusions:

The factors that affect the choice of solver types : number of obstacles, the size of obstacles, and the distance between two adjacent obstacles. For obstacles that are spaced closely, the preferable choice is Lapack (in sequential codes) and its parallel version Scalapack. For 2000 obstacles, for a visualization on a 800 × 800 grid 6 , Scalapack (on 16 cores) and postprocessing with Hermit interpolation (on 16 cores) require 1 min 10 secs, c.f. subsection 6.4. When the obstacles are further apart, the preferable choice is GMRES with either LUSGS or SGS preconditioners, which outperform Lapack, and are at head-to-head with Scalapack (run on several processors). As noted in the introduction, the codes for the iterative solvers are currently sequential; hence, there is a possibility to improve further their performance with parallelization. Currently, for the most difficult case with very large number of closed together obstacles, Scalapack is the only7 choice that can stilll run and gives result; e.g. for 10 4 obstacles, the pre-processing time on 48 cores is ∼ 24 mins , c.f. Subsection 6.6.

Ratio Obstacle Radius

Obs Distance = 0.03 0.30 = 0.1. In the above comparison, we obtain convergence for all options (with or without preconditioning) of GMRES. In general, the direct solvers take less time than iterative ones, with the exception of the SGS and LUSGS family which are head-to-head with Mumps. However, the iterative solvers give lower quality density compared to that produced by direct solvers (with final error at 10 -12 ). The error tolerance for GMRES (currently set at 10 -6) is not a good indication for the final error (10 -3 ) of the density. The various choices of preconditioners fall into 4 distinct performance groups, with big differences in terms of number of iterations between them:

Ratio

• no pre-conditioning takes 820 iterations;

• preconditining with Jacobi takes ∼ 500 iterations: Left Jacobi, Left2ndJacobi, Right Jacobi, Right2ndJacobi;

• preconditioning with one 'levels' of Gauss-Seidel takes ∼ 199 iterations: LeftFGS, LeftBGS, RightFGS, RightBGS;

• the most efficient group is preconditioning with two 'levels' of Gauss-Seidel which takes ∼ 76 iterations: LeftSGS, LeftLUSGS, RightSGS, RightLUSGS.

As a result, the most appropriate choice of preconditioners is the group with two 'levels' of Gauss-Seidel.

In the next numerical experiments, we investigate whether increasing the quality of the density produced by the iterative solvers will be too costly, this is done by increasing the GMRES error tolerance from 10 -6 to 10 -9 .

GMRES parameters : Precision tol 1.e-9 ; Maximum number of iterations = 2000 ; Size of Krylov space = restart = 100 . In this case, while increasing the GMRES tolerance requires more iterations, the time increase is small and is compensated by a great increase in the quality of density, having relative error at 10 -6 , compared to 10 -3 for the last setting.

Name

Case 2000 obstacles

For the following experiments, we study

• whether varying the stop criteria (precision tolerance, Niter Max, and number of restart) can help improve the convergence of GMRES with the two 'levels' of Gauss-Seidel preconditioners, for cases with large number of obstacles ;

• whether the performance (time cost) of Mumps can be improved, if it is run in parallel;

• for certain comparisons, we also use a different criteria, L 2 of the scattered field on a bounded domain (containing the obstacles), to see whether the Sobolev norm of the density is too pessimistic as an error criteria.

Parameters

Angle of incidence of plane wave = 90.0. FS-SL method order =2.

Domain of visualization

Reference solutions : The reference solution is obtained by Mumps (for linear solver) and Exact evaluation (for post-processing), see Figure 14. For parallel execution, the number after -n indicates the number of processors being used. We note that parallel execution of Mumps does not improve the pre-processing time. This is perhaps due to the fact that Mumps is not suitable for dense matrices. We also see that for a precision tolerance of 10 -6 , Hermite interpolation cuts down drastically (by 10) the postprocessing time. With the current choice of GMRES, SGS family is head-to-head with Mumps.

Method

Comparison 2 :

The following comparison is in relative difference in Sobolev norm of (single layer) density. GMRES parameters : Precision tol = 1.e-8; Maximum number of iterations = 5000; Size of Krylov space = restart = 400 . Here, we see that the relative error of the density is much more pessimistic than the relative L 2 error of the final scattered field on a bounded domain. In these experiments, we introduce the solver Scalapack, which greatly enhances the efficiency of the method.

Name

Parameters

Angle of incidence of plane wave = 90.0.

Ratio Obstacle Radius

Obs Distance = 0.03 0.30 = 0.1.

Ratio Obstacle Distance Wavelength = 0.3 0.63 ∼ 0.48.
Visualization is carried out on the structured grid 800 × 800.

For following comparison, the reference solution is produced with Lapack + Exact Eval. The parameters accompanying GMRES methods are (Error Tol, No of Iter Max, No Restart). For parallel execution, the number after -n indicates the number of processors being used. In order to simulate multiple scattering by 2000 small and closed-together obstacles, the optimal choice is using Scalapack and Hermite interpolation, which gives a visualization of the solution on a 800 × 800 grid, taking a total of 1 min and 10 secs.

Case 2000 obstacles -Large distance and small obstacles

From previous tests, we see that the direct solvers outperformed the iteratives ones. We next investigate a setting where one can see the advantages of iterative ones. It turns out that when the obstacles are further distanced than in the previous cases, the number of iterations needed drastically dropped, which places iterative solvers (written in sequential codes) at the same level of performance as Scalapack, the parallel version of Lapack, run on 16 cores.

Parameters

Angle of incidence of plane wave = 90.0.

Value for b and a (in the notation of the codes) : b = 1.00, a = 0.01

Domain

Observations

• The current Hermite interpolation which reduces the post-processing time by a third, is satisfactory for a precision toleration of 10 -5 .

• For the current configuration in which the obstacles are further distanced, the number of iterations for GMRES drastically dropped, taking only 56 iterations. This is quite impressive, since this is what is needed for 200 closed-together obstacles. In fact, GMRES with LUSGS preconditioner is faster than Lapack, and is head-to-head with Scalapack run on 16 processors.

• The distance between obstacles has small impact on the direct solvers. Their performance time stays approximately the same compared to previous case : 1 min 12 secs .

As a result, for the case where the obstacles are far-way, under low precision, one can have the option of using GMRES with either LUSGS or SGS as preconditioners. However, as soon as the obstacles are close together, direct solvers are more efficient.

Inria

Case 10000 obstacles

For very large number of obstacles starting 10000, the only solver that can still run and give result is Scalapack, see Figures 19 and20 for a visualization of the solution with the following parameters.

Parameters

Angle of incidence of plane wave = 90.0. We see that for a low definition visualization for 1000 obstacles, the method implemented on 48 processors takes 24 mins 40 secs. This means that with a very large number of obstacles, unless GMRES is optimized further to handle large-sized linear system, the only current option is the robust Direct solver Scalapack. As observed from previous experiment, the direct solver is less sensitive to the distance between obstacles.

Domain of visualization

Conclusions and Future Problems

We have seen that the FS-SL method is very robust in solving the multiple scattering problem for small circular obstacles in homogeneous media. Being mesh-free, it can handle efficiently a large number of obstacles on an infinite homogeneous domain, while allows the size of the obstacles to be very small compared to the incident wave. In addition, the linear systems generated by the method have simple definition, and thus enable easy coding and implementation. We also observe that, for optimality, different formations of obstacles require different solver types. Direct Solvers (Lapack and Scalapack) are more efficient than iterative solvers in treating the cases where the obstacles are closed together. The current codes using Scalapack can handle up to 10 4 small and closely spaced obstacles. On the other hand, iterative solvers regain in performance when the obstacles are far apart; in particular, GMRES methods with preconditioners LUSGS and SGS outperform Lapack and are head-to-head with parallel Scalapack. Below we discuss ideas that would widen the range of application and increase the efficiency of FS-SL method even further.

Improvement for codes written for disc-geometry We list ideas that can readily improve the efficiency of the method in general, and in particular that of iterative solvers. Firstly, the codes for the iterative solvers are sequential; a parallelization would allow to launch this family of solvers on several processors. Secondly, we have not taken advantage of the special structure of the off-diagonal blocks of the coefficient matrix, which can be written as the product of two diagonal matrices and a Toeplitz one, as used in [START_REF] Antoine | On the numerical approximation of high-frequency acoustic multiple scattering problems by circuler cylinders[END_REF] and [3]. This idea would lower the CPU storage. Thirdly, we could sparsify the blocks (of the preconditioners) associated with the interaction between obstacles that are far-apart. These improvements would allow the method to handle even larger-sized problems, especially for cases where direct solvers will fail.

Future questions

Immediate extension and application

The current transmission problem can be extended to circular solid inclusions. In its current form, the codes can be incorporated to solve inverse problems, e.g. in the reconstruction of the position, size and number of defects (in the form of holes) in a homogeneous material.

General geometry

In order to maintain the robustness observed for disc geometry (which also includes solid inclusions), an efficient integration quadrature rule will be needed to handle the weak singularity of the Green kernel.

Toward inhomogeneous media In the transmission problem we have considered, the wave speed is a step function, and is thus a special case of variable wave speeds having sharp discontinuities at the boundaries of the obstacles. We have seen that it can be treated efficiently with the FS-SL method, which is in the family of boundary integral equation. For wave speed that is constant everywhere except on a compact domain where it varies continuously, the general approach is to use volume integral equation methods, e.g. the Lippmann-Schwinger equation approach [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF].

On the other hand, despite its shortcomings in the current context, the general consensus is that FEM provides more flexibility in treating media of complex heterogeneities than Integral Equations Method (IEM); e.g. for strongly (continuous varying) inhomogeneities, the most common technique found in literature is homogenization based on Finite Element, called Finite Element Heterogeneous Multiscale Method. Since our future goal is to study material whose wave speed Inria is homogeneous everywhere except on a compact region where it varies, it would be promising to combine FEM with IEM to take advantage of the strength of both methods.

A Layer Potential Theory

Trace operaters : By our convention, the normal vector n(x) points outward. For f ∈ H 2 , define the normal derivative associated to normal vector n,

∂ ∂n f = lim h→+0 n(x) • ∇f (x -h n(x)).
In terms of γ 0,int and γ 0,int , this can be written as,

γ 1,int f := (γ 0,int ∇f ) • n , γ 1,ext f := (γ 0,ext ∇f ) • n.
Definitions of jumps:

f := γ 0,ext f -γ 0,int f ; ∂ n f := γ 1,ext f -γ 1,int f = ∇f • n.
We follow the convention of signs from [START_REF] Kress | Linear Integral Equation[END_REF] and [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF].

Layer potentials : Let Γ be a simple closed curve ∈ C 2 , i.e. Γ forms the boundary of a simply connected region Ω.

Define acoustic single-layer potential with density φ ∈ C(Γ)

Sκ,Γ φ (x) := Γ φ(y) G κ (x, y) dσ(y) , x ∈ R 2 \ Γ .
where G κ is the fundamental solution of the Helmholtz equation at wavenumber κ 2 G κ (x, y) := i 4 H

(1) 0 (κ |x -y|) , x = y .

Define acoustic double-layer potential with density φ

DΓ φ (x) := Γ φ ∂ ∂n(y) G κ (x, y) dσ(y) , x ∈ R 2 \ Γ .
Surface potentials :

(S Γ,κ φ) (x) := Γ φ(y) G κ (x, y) dσ(y) , x ∈ Γ ; (D Γ,κ φ) (x) := Γ φ ∂ ∂n(y) G κ (x -y) dσ(y) , x ∈ Γ ; D Γ,κ φ (x) := Γ φ ∂ ∂n(x) G κ (x -y) dσ(y) , x ∈ Γ ; (T Γ,κ φ) (x) := ∂ ∂n(x) Γ φ ∂ ∂n(y) G κ (x -y) dσ(y) , x ∈ Γ .

Properties of layer potentials

For -1 ≤ s ≤ 1, S : H s-1/2 (Γ) -→ H s+1 loc (R 2 ) bounded ; D : H s+1/2 (Γ) -→ H s+1 loc (R 2 \ Ω) bounded .
They both give solution to the Helmholtz equation on R 2 \ Γ, which satisfy the outgoing Sommerfeld radiation condition.

Trace identities : Following the notation of Colton Kress theorem 3.1 , p 39 [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF], see also [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF]

[Section 7.1]. For -1 ≤ s ≤ 1, γ 0,Γ,ext SΓ,κ φ = γ 0,Γ,ext SΓ,κ φ = S Γ,κ φ , φ ∈ H -1/2+s (Γ) ; γ 1,int D = γ 1,ext D = T , φ ∈ H 1/2+s (Γ) . γ 1,int S = D + 1 2 Id , γ 1,ext S = D -1 2 Id , φ ∈ H -1/2+s (Γ) . γ 0,int D = D -1 2 Id , γ 0,ext D = D + 1 2 Id , φ ∈ H 1/2+s (Γ) . (68)
Thus we obtain the jump identities,

[ Sψ] = 0 in H 1/2 (Γ) ; [γ 1 Sψ] = -ψ in H -1/2 (Γ) . ( 69
)
Mapping properties of surface potentials : We cite the following useful theorem from

[5][Thm 7.3], 1. S i is H -1/2 (Γ) -coercive, i.e S i φ , φ H 1/2 (Γ) , H -1/2 (Γ) ≥ C φ 2 H -1/2 (Γ) .
2. If κ 2 is not an eigenvalue of -∆ in Ω, then S κ : H -1/2 (Γ) → H 1/2 (Γ) is an isomorphism, i.e with bounded inverse.

For

-1 ≤ s ≤ 1, S κ : H s-1/2 (Γ) -→ H s+1/2 loc (Γ) bounded S κ : H s-1/2 (Γ) -→ H s-1/2 loc (Γ) bounded and compact S κ -S i : H -1/2 (Γ) -→ H 1/2 (Γ) bounded and compact D κ : H s+1/2 (Γ) -→ H s+1/2 loc (Γ) bounded and compact D κ : H s-1/2 (Γ) -→ H s-1/2 loc (Γ) bounded and compact T κ : H s+1/2 (Γ) -→ H s-1/2 loc (Γ) bounded . (70)

B Multipole Expansions

Notations of polar coordinates : Polar coordinates of x relative to the origin 0 R 2 is given by |x| and θ

(x) = θ 0 R 2 (x) x = |x| cos θ(x) , sin θ(x) . (71)
Polar coordinates relative to x ∈ R 2 ,

x = x + r x (x) cos θ x (x) , sin θ x (x) ; r x (x) = |x -x|. (72)

Inria

B.1 Multipole expansion for solutions of Helmholtz equation

First we recall the series reprentation of solutions of Helmholtz equation (-∆ -κ 2 )u = 0.

Proposition 8.

1. If u solves the Helmholtz equation in the ball |x| < R 0 for some R 0 , then u is of the form, for some coeffcients c nm (κ), for r with 0

< r < R 0 u(rθ) = ∞ n=-∞ c n (κ) J n (κr) e i n θ . 2. If u solves the Helmholtz in the annulus R 0 < |x| < R 1 with 0 ≤ R 0 < R 1 then u is of the form for r with R 0 < r < R 1 u(rθ) = ∞ n=-∞ a n (κ) H (1) n (κr) + b n , H (2)
n (κr) e i n θ .

In the

case R 1 = ∞, if u is κ-outgoing then b n (κ) = 0, while if u is κ-incoming then a n (κ) = 0 .

B.2 Graf 's addition theorem

We next consider two points x 1 , x 2 ∈ R 2 . Relative to x 2 , x 1 is given by

x 1 = (|x 1 -x 2 | , θ x2 (x 1 )) .
Relative to x 1 , x 2 is given by

x 2 = (|x 1 -x 2 | , θ x1 (x 2 )
) . We simplify the notations, by writing, for i = 1, 2

r i = r i (x) = r xi (x) ; θ i = θ i (x) = θ xi (x) .
We cite Graf's addition in [START_REF] Martin | Multiple Scattering : Interaction of Time-Harmonic Waves with N Obstacles[END_REF][Thm 2.12]. We would like to express the m-th order multipole centered at x 2 , in terms of the polar coordinates relative to x 1 , i.e, in terms of r 1 (x) and θ 1 (x). Theorem 9.

1. For r 1 (x) < x 1 -x 2 ,

H (1) m (κ r 2 ) e imθ2 = ∞ n=-∞ H (1) m-n (κ x 1 -x 2 ) e i(m-n)θx 2 (x1) J n (κ r 1 )e inθ1 . 2. For r 1 > x 1 -x 2 , H (1) m (κ r 2 ) e imθ2 = ∞ n=-∞ J m-n (κ x 1 -x 2 ) e i(m-n)θx 2 (x1) H (1) n (κ r 1 )e inθ1 .
The above identities can also be written in terms θ x1 (x 2 ) using the following relation between the two relative polar angles.

Remark 7. The relation between the two angles is given by,

|θ x2 (x 1 ) -θ x1 (x 2 )| = π.
Thus we have e inθx 2 (x1) = e inθx 1 (x2) e ±inπ = (-1) ±n e inθx 1 (x2) . As a result, we have the identity e inθx 2 (x1) = (-1) ±n e inθx 1 (x2) .

(73)

B.3 Multipole expansion for planewave

We have the Jacobi-Anger expansion, see for e.g [20, eqn (2.17)],

e i t cos ϕ = ∞ k=-∞ i k J k (t) e ikϕ . (74)
We assume that the incident field is given by a (time-harmonic) acoustic plane wave incident at angle α inc , i.e u inc (rθ) = e i κx •(cos αinc,sin αinc) = e i κ r cos(θ-αinc) .

In the current notation,

x = x + r x (x) cos θ x (x) , sin θ x (x) ; r x (x) = |x -x| ,
we can rewrite u inc in polar coordinate with respect to an arbitrary x ∈ R 2 , u inc (r θ) = u inc (x) e i κ rx(x) cos(θx(x)-αinc) .

The multipole expansion of the plane wave relative to the origin 0 R 2 and to x ∈ R 2 are given correspondingly by

u inc (x) = ∞ m=-∞ i m J m (κ r) e i m(θ-αinc) = u inc (x) ∞ m=-∞ i m J m (κ r x (x)) e i m(θx(x)-αinc) . (76)
As a result, we obtain Lemma 10. With u inc given by (54), for l ∈ Z, we have

γ 0,Γ J u inc = u inc (x J ) ∞ m=-∞ i m J m (κ r J ) e i m(θx J (x)-αinc) .

B.4 Multipole expansions for single layer with potential w I,l

For Γ = ∂B(x, r), the acoustic single layer potential with continuous density w m is defined as

u m = SΓ w m = Γ G κ (|x -y|)w m ds(y).
Lemma 11. In terms of the polar coordinates (r, θ) centered at x, with r = r x (x) and θ = θ x (x), 1. The acoustic single layer potential u m is given by

u m (r, θ) = iπ r 2 e i m θ × H (1) m (κ r) J m (κr) , r ≤ r J m (κ r) H (1) m (κr) , r > r . ( 77
)
2. The trace of u m on Γ is then given by

(γ 0,Γ u m ) (θ(x)) = S Γ w m = iπ r 2 e i m θ H (1) m (κ r) J m (κ r). (78)
Inria 3. The exterior first order trace of u m on Γ is then given by (γ 1,e u l ) (θ(x)) = iπ r κ 2 e i l θ(x) J l (κ r) H

(1) l (κr) .

Proof. Being a single layer, u m satisfies the Helmhotz in R n \ Γ, with u m κ-going at infinity and regular at x, c.f. Appendix A. By Proposition 8, in terms of the polar coordinates centered at x, u m is of the form

u m (r, θ) = ∞ n=-∞ α n (κ) J n (κr) e inθ , r < r ; ∞ n=-∞ β n (κ) H (1)
n (κr) e inθ , r > r ;

and

∂ r u(r, θ) = ∞ n=-∞ κ α n (κ) J n (κr) e imθ , r < r ; ∞ n=-∞ κ β n (κ) H (1)
n (κr) e inθ , r > r .

Secondly, using the formula for the jumps of single layer operator, c.f. ( 68), u m satisfies the following transmission conditions along Γ,

u Γ = 0 ; ∂ n u Γ = -w m . Thus β n (κ) H (1)
n (κ r) -α n (κ) J n (κ r) = 0; κ β n (κ) H (1) n (κr) -κα n (κ) J n (κr) = -δ nm . (81)
With W (f, g) denoting the Wronskian of two functions (f, g), by [START_REF] Lebedev | Special Functions and Their Applications[END_REF], we have

W (J n , H (1) n )(x) = J n H (1) n -J m H (1) n = 2i πz .
This means that the system (81) is invertible. As a result of this, the sum in (79) only contains the term at level m; for levels n with n = m,

α n (κ) = β n (κ) = 0 , n = m.
At the level m, α m (κ) and β m (κ) solve the linear problem

J n (κ r) -H (1)
m (κ r) J m (κ r) -H (1) m (κ r) α m (κ) β m (κ) = 0 1/κ .
As a result, we have

α = 1 W (J m , H (1)
m )(κ r)

det 0 H (1)
n (κ r)

1 κ H (1) n (κ r) = - π r 2i H (1) m (κ r) = iπ r 2 H (1) m (κ r); β = -1 W (J m , H (1)
m )(κ r)

det J m (κ r) 0 J m,(1) (κ r) 1 κ = - π r 2i J m (κ r) = iπ r 2 J m (κ r).
Plugging in the above form of α and β into (79), we obtain the formula (77) for u m .

To obtain the trace of u m at Γ, we let r = r in (77). Lastly, since the convergence of the series is uniform on compact subsets, this allows us to take derivative in the radial variable to obtain the normal derivative trace. Denote by Γ x = ∂B(x, r) .

For f = f • ϕ J f (θ) = ∞ m=-∞ a m e imθ
We would like to calculate the coefficient of the Fourier expansion of G f • ϕ Γx , i.e

(G f • ϕ Γx ) (θ) = ∞ m=-∞ a m e imθ
where a m is given by

1 2π 2π 0 i 4 R 2 y H (1)
0 (κr x(r,θx) (y)) f (y) dy ds(Γ x ) .

(85) note that since x ∈ Γ x ,

x = x(r, θ x ) = x + r(cos θ x (x), sin θ x (x)) .

Lemma 13. For F where Supp F ∩ B(x, r) = ∅, we have

(G F ) • ϕ Γx (θ) = ∞ m=-∞ a m e imθ
where a m is given by a m = i 4 (-1) l J l (κ r)

Supp F H (1)
-l (κ r x (y)) e -ilθx(y) F (y) dy .

Proof. We consider the double integral in (85) without the constant.

C GMRES Solvers

C.1 GMRES with no conditioning

We also refer reader to [START_REF] Frayssé | A set of gmres routines for real and complex arithmetics on high performance computers[END_REF] for a brief general description of the GMRES algorithm.

Intuition : Denote by q(•) the minimal polynomial for A. By definition, we have 0 = q(A) = α 0 I + . . . α m A m .

Note that α 0 = 0 if and only if A is nonsingular, and under this assumption, we can write its inverse as,

A -1 = - 1 α 0 m-1 j=0 α j+1 A j .
As a result, the unique solution to Ax = b, denoted by x = A -1 b, lies in the Krylov space K m (A, b).

Approach 1 : The above calculation shows that we should look for x h which approximates x in the Krylov spaces K n (A, b). Specifically, we look for x h of the form x h = K n c for some vector c such that the following residue is minimized,

r h 2 = Ax h -b 2 = AK n c -b 2 .
For this, one could use a QR-factorization of A; however, this can be unstable and expensive. More efficiently, one uses the Arnoldi iteration with matrix A and initial vector b to produce an orthonormal basis for the Krylov space K n , c.f. Appendix C.5. We denote this basis by Q n . Moreover, by (95), we have

A Q n = Q n+1 Hn . ( 88
)
In terms of the orthonomal basis Q n , we look for x h of the form

x h = Q n y ,
for some vector y so that the following residue is minimized,

Ax h -b 2 = AQ n y -b 2 (88) = Q n+1 Hn y -b 2 .
Since Q n+1 is unitary and multiplication by a unitary matrix does not change the norm • 2 ,

Ax h -b = Q * n+1 Q n+1 Hn y -Q * n+1 b 2 = Hn y -Q * n+1 b 2 . ( 89
)
This is an improvement from optimizing with AQ since AQ n is of size m × n while Hn is only (n+1)×n. Moreover, the right hand side Q * n+1 b has a much simpler expression. By the definition of Q n , we have,

Q * n+1 b =    q 1 b . . . q * n+1 b   
Since q j -s form an orthonormal basis, and q 1 = b b 2 , we can simplify Q * n+1 b to

Q * n+1 b = b e 1 .
With this, (89) simplifies to and retrieve the approximating solution x h to x by

x h = Q n y .
The error (residue) associated with x h is given by

Ax h -b 2 = Hn y -b e 1 2 .
Approach 2 : Instead of looking for x , we can start with an initial guess x 0 and approximate the correction, i.e x = x 0 + p .

The correction p satisfies

p = A -1 b -A -1 Ax 0 = A -1 r 0 ; r 0 = b -Ax 0 .
This means that we look for p which solves A p = r 0 ; r 0 = b -Ax 0 .

In another word, we approximate the exact solution x * = A -1 b by x h of the form

x h = x 0 + p h .
Look for p h in K n (A, r 0 ) so that the following residue is minimized, i.e

Ap h -r 0 2 = min p∈Kn(A,r0)

Ap -r 0 2 .

For the minimization problem, we use the Arnoldi iteration, as in Approach 1. Specifically, we look for p h of the form p h = Q n y. We have Ap h -r 0 2 = AQ n y -r 0 2 = Q n+1 Hn y -r 0 2 .

Inria

The last equality comes from the Arnoldi process AQ n = Q n+1 Hn . Since Q n+1 is unitary and multiplication by a unitary matrix does not change the norm • 2 , we have

Ap h -r 0 2 = Q * n+1 Q n+1 Hn y -Q * n+1 r 0 2 = Hn y -Q * n+1 r 0 2 .
This is an improvement compared to optimisation with AQ n since AQ n is of size m × n while Hn is only (n + 1) × n. Moreover, the right hand side Q * n+1 r 0 has a much simpler expression. By the definition of Q n , we have,

Q * n+1 r 0 =    q 1 r 0 . . . q * n+1 r 0    .
Since q j -s form an orthonormal basis, and q 1 = r 0 r 0 2 , we can simplify Q * n+1 r 0 to

Q * n+1 r 0 = r 0 e 1 .
As a result, we are left with the minimization problem Hn y -r 0 e 1 and retrieve the approximating solution x h to x by

x h = x 0 + p h ; p h = Q n y .
The error (residue) associated with x h is given by Ax h -b 2 = Hn y -r 0 e 1 2 .

Algorithm for approach 2 : We work with Krylov spaces of dimension m.

1. Start : Initial guess x 0 . Compute r 0 = b -Ax 0 . We work with the K m (A, r 0 ). Define q 1 := r 0 r 0 .

2. Use Arnoldi iteration process to construct Q m and H m for K m (A, r 0 ). For j = 1, . . . , m (a) Compute w := Aq j .

(b) For k = 1, . . . , j : do h ij = (w, q i ) ; w := w -h ij q i .

(c) h j+1,j = w 2 and q j+1 = w h (j+1)j . 4. Restart : if x h satisfies the stop criteria then stop, else set x 0 as x h and go to step 2.

RR n°8988

C.2 GMRES with left conditioning

Instead of solving A p m = r 0 , we solve P -1 Ap m = P -1 r 0 ; r 0 = b -Ax 0 .

In other words, we approximate x = A -1 f by an element of the form

x h = x 0 + K m (P -1 A , P -1 r 0 ).

This means, that we look for solution in the Krylov space K m (P -1 A , P -1 r 0 ) that minimizes p m = argmin p∈Km(P -1 A , P -1 r0)

P -1 Ap -P -1 r 0 2 .

We next carry out the Arnoldi process for P -1 A with initial vector P -1 r 0 . We will only calculate its action when needed (in two steps), and not the matrix P -1 A itself.

Algorithm : We work with Krylov spaces of dimension m.

1. Start : Start with Initial guess x 0 . Compute r 0 = b -Ax 0 . Since we work with Krylov space K m (P -1 A , P -1 r 0 ), we define q 1 = P -1 r 0 P -1 r 0 2 .

2. We use Arnoldi iterations to obtain orthogonal space Q m = (q 1 , . . . , q m ) and Hessenberg matrix H m for K m (P -1 A , P -1 r 0 ). (c) For k = 1, . . . , j : do h ij = (w, q i ) ; w := w -h ij q i .

(d) h j+1,j = w 2 and q j+1 = w h (j+1)j .

3. Solving the least square problem for y m which minimizes r 0 2 e 1 -Hm y m 2 = min y ∈ Km(P -1 A , P -1 r0) r 0 2 e 1 -Hm y 2 .

Compute the approximate solution

x m = x 0 + p m ; p m = Q m y m .

While the actual residue/error is Ax -b = Ap m -r 0 , with left preconditioning, we only have the preconditioned residue at our disposal, P -1 (Ap m -r 0 ) 2 = r 0 2 e 1 -Hm y m 2 .

5.

Restart : if the preconditioned residue satisfies the stop criteria, then stop, else set x 0 as x m and go to step 1.

Inria

C.3 GMRES with right conditioning

Instead of solving A p m = r 0 , we solve A P -1 P p m = r 0 .

In another word, we approximate x = A -1 f by an element of the form

x m = x 0 + P -1 K m (AP -1 , r 0 ) .
This means, that we look for pm in the Krylov space K m (AP -1 , r 0 ) that minimizes pm = argmin p∈Km(AP -1 , r0)

A P -1 p -r 0 2 .

Algorithm : We work with Krylov spaces of dimension m 1. Start : Start with Initial guess x 0 .Compute r 0 = b -Ax 0 . Since we work with the Krylov space K m (AP -1 , r 0 ), define q 1 = r 0 r 0 2 .

2. We use Arnoldi iterations to construct orthogonal space Q m = (q 1 , . . . , q m ) and Hessenberg matrix H m for K m (AP -1 , r 0 ).

For j = 1, . . . , m : do (a) Compute z j = P -1 q j ;

(b) Compute w := A z j ;

(c) For k = 1, . . . , j : do h ij = (w, q i ) ; w := w -h ij q i .

(d) h j+1,j = w 2 and q j+1 = w h (j+1)j .

5.

Restart : if the residue satisfies the stop criteria then stop, else set x 0 as x m and go to step 1.

In particular, the n-th column of AQ n satisfies Aq n = h 1n q 1 + h 2n q 2 + . . . + h nn q n + h n+1,n q n+1 ; ⇒ q n+1 = Aq n -n j=1 h jn q i h (n+1)n ;

⇒ q * k q n+1 = q * k Aq n -n j=1 h jn q * k q i h (n+1)n , 1 ≤ k ≤ n .

Since we require Q n to be orthonormal, we have for all

0 = q * k Aq n -h kn h (n+1)n ⇒ h kn = q * k Aq n , 1 ≤ k ≤ n.
In particular for k = n, this is the Rayleigh quotient. For n = 1, we have

0 = q * 1 Aq 1 -h 11 q * 1 q 1 ⇒ h 11 = q * 1 Aq 1 q * 1 q 1 .
Algorithm : This is the modified Gram-Schmidt implementation of the Arnoldi algorithm, c.f. [START_REF] Liesen | Krylov subspace methods: principles and analysis[END_REF][Algorithm 2.4.2]. The algorithm produces orthonormal vector q 1 , . . . , q d such that span {q 1 , . . . , q n } = K n (A, y) , 1 ≤ n ≤ d .

1. Start : Define q 1 := y y 2 .

2. For n = 1, 2, . . . w := A q n ; For i = 1, . . . , n, do : h i,n := q * j w ; w := w -h i,n q i . h n+1,n = w 2 if h n+1,n = 0 then stop. Else q n+1 := w h n+1,n .

Remark 10. The most expensive operation in the algorithm is matrix vector multiplications. As a result, the robustness of the method relies on an efficient implementation of the matrix-vector product which should be tailored to the problem. The algorithm does not store matrix A and only needs the action of A on q n . By induction, one can show that after n steps, the above algorithm produces a Hessenberg matrix H n of size n × n and an orthogonal matrix Q n of size N × n, satisfying

A Q n = Q n H n + h n+1,n q n+1 e t n = Q n+1 Hn ; Hn = H n 0 . . . 0 h n+1,n . (95)
At the d-th iteration, we have h d+1,d = 0, and we obtain

AQ d = Q d H d .

Figure 1 :

 1 Figure 1: Scattering of a planewave by various and non-overlapping obstacles.

[ 5 ]

 5 [Thm 5.22]. We have the equivalence between the BVP-s to the corresponding one with Ω ext replaced by B R , and the outgoing radiation condition replaced by[START_REF] Frayssé | A set of gmres routines for real and complex arithmetics on high performance computers[END_REF], c.f. [5][Lemma 5.24]. One next shows the existence and uniqueness of weak (variational) solutions to the equivalent problems, by showing the coercitivity of the resulting variational forms. For Dirichlet problem, see e.g. [5][p.105], for impedance problem [12][Lemma 2.1] under the same assumption as in Theorem 1, for transmission problem [12][Thm 3.1] under the assumption Re µ > 0 , Im (µ κ e ) ≥ 0 , Im (ρ κ e κ 2 I ) ≥ 0 . 3 Single Layer Potential formulation of the multi-scattering problem Notations : For 1 ≤ I ≤ N Obs , we have denoted by Ω I the region occupied by the I-th obstacle, and by Γ I its boundary, i.e Γ I = ∂ Ω I . We simplify the notations of the trace operator along Γ I (introduced in previous section) by writing γ 0,I,ext = γ 0,Γ I ,ext , 1 ≤ I ≤ N Obs ;

S

  IJ,κe ṽJ .

Proposition 2 .

 2 If (v I ) 1≤I≤N Obs satisfies (24), then {u, u 1 , . . . , u N Obs } defined by Equations (25) solves the transmission problem (5) with incident source u inc . Proof. The essential ideas of the proof follow those of [13][Theorem 5.1]. As defined, u satisfies the Helmholtz equation on exterior domain Ω ext and the outgoing radiation condition (1), while u I , 1 ≤ I ≤ N Obs , solves the Helmholtz equation on Ω I . It remains to verify that they satisfy transmission conditions.

2 e

 2 , u = 0 in Ω ext , for each case α = D, I, N and T . Dirichlet : That ψ is in the kernel of A D means that γ 0,I u = γ 0,I,int u = γ 0,I,ext u = 0 , I = 1, . . . , N Obs . This means u solves the exterior Helmholtz equation (-∆ -κ 2 e )u = 0 with zero homogeneous boundary. By the well-posedness of the exterior Dirichlet BVP for Helmholtz, this means u = 0 for Ω ext . Impedance (Neumann for λ = 0) : For α = I and N , that ψ is in the kernel of A α means that (γ 0,I,ext + i λ γ 1,I,ext ) u = 0 , I = 1, . . . , N Obs . This means u solves the exterior Helmholtz equation (-∆ -κ 2 e ) u = 0 with zero homogeneous boundary. By the well-posedness of the exterior Impedance (Neumann for λ = 0) BVP for Helmholtz, this means u = 0 for Ω ext .

  l interlace. Denote by N D and N N the set of obstacles whose interior wavenumbers are Dirichlet and Neumann EVs, correspondingly. The remaining obstacles form the set N T , i.e. {1, . . . , N Obs } = N D ∪ N N ∪ N T .

( 1 )

 1 l (κ e r I ) + κ e H

Inria

Figure 2 :

 2 Figure 2: Convergence curve with respect to order of approximation (10 randomly distributed obstacles). Obstacle Radius Wavelength

Figure 4 :

 4 Figure 4: Convergence curve with respect to order of approximation (10 randomly distributed obstacles). Obstacle Radius Wavelength

Figure 5 :Inria 5

 55 Figure 5: Convergence curve with respect to size of obstacles (10 randomly distributed obstacles). OrderApproximation : = 2 ; = 3 ; = 4 ; = 6

Parameters

  Angle of incidence of plane wave = 0 ; Domain of interest [-34, 34] × [-40, 60] ; Wavenumber κ = 1.0 ; Wavelength λ = 2π ; Number of obstacles = 6 ; out on the structured grid 200 × 200.

Figure 6 :

 6 Figure 6: Mesh used by Montjoie (created from geo2mesh).

Figure 7 :

 7 Figure 7: Numerical Convergence for 6 holes, showing consecutive relative error.

  L 2 Rel = 2.29 e-8

Figure 8 :

 8 Figure 8: Reference solutions with Precision tolerance at 1.e-12, for 6 obstacles.

Figure 9 :Figure 10 :

 910 Figure 9: Comparison at precision 1.e-6, for 6 obstacles, parameters in subsection 5.1.

Figure 11 :

 11 Figure 11: Numerical Convergence for 200 obstacles, showing consecutive error.

L 2

 2 Rel = 3.38 e-8

Figure 12 :

 12 Figure 12: Reference Solutions (at precision tolerance 1.e-8), for 200 obstacles, parameters in subsection 5.2.

Figure 13 :

 13 Figure 13: Comparison at precision 1.e-3 for 200 obstacles, parameters in subsection 5.2.

6. 1

 1 Test 200 holes Parameters Angle of incidence of plane wave = 90.0. Value for b and a (in the notation of the codes) : b = 0.15, a = 0.03. Domain of visualization [-1.05, 18.5] × [-1.05, 13.35]. Wavenumber κ = 10.0 ; Wavelength λ = 2π/10 ∼ 0.63. Number of obstacles = 200. Radius of obstacles = 0.03 ; Distance between two adjacent obstacles = 0.30 Ratio Obstacle Radius Wavelength = 0.03 2π/10 ∼ 0.048.

  [-1.05, 22.9] × [-1.05, 25.9]. Value for b and a (in the notation of the codes) : b = 0.15 ; a = 0.03. Wavenumber κ = 10.0 ; Wavelength λ = 2π/10 ∼ 0.63. Size matrix = 10 4 × 10 4 . Visualisation on structured grid of size 800 × 800. Number of obstacles = 2000. Radius of obstacles = 0.03 Distance between two adjacent obstacles =

  Relative difference in L 2 -norm compared between Exact and Inter (Parallel -n8) evaluation

Figure 14 :

 14 Figure 14: Reference solution (Mumps) for 2000 obstacles (with Parameters in subsection 6.3).

Inria

Figure 15 :Comparison 3 :

 153 Figure 15: Comparison 1: Relative difference in L 2 -norm compared with Mumps + Exact Eval for 2000 obstacles

Value for b 9 ]

 9 and a (in the notation of the codes) : b = 0.15 ; a = 0.03; Domain of visualization [-1.05, 22.9] × [-1.05, 25.Wavenumber κ = 10.0 ; Wavelength λ = 2π/10 ∼ 0.63. FS-SL method order =2. Size matrix = 10 4 × 10 4 . Number of obstacles = 2000. Radius of obstacles = 0.03 Distance between two adjacent obstacles =

Figure 17 :Figure 18 :

 1718 Figure 17: Lapack + Exact Eval on 800 × 800 structured grid for 2000 obstacles spaced apart, for Parameters listed in Subsection 6.5.

Figure 19 :

 19 Figure 19: Visualisation of solution for 10000 obstacles obtained by FS-SL + Scalapack on 1600 × 1600 structured grid. Relative L 2 difference between exact and interpolated version = 8.29 e -6.

Figure 20 :

 20 Figure 20: Visualisation of solution obtained by FS-SL + Scalapack on 800 × 800 structured grid. Relative L 2 difference between exact and interpolated version = 1.81 e -5.

f

  (r J (cos θ, sin θ)) e -imθ dt.The solution to(-∆ -κ 2 )u = F (x) , x ∈ R 2is given byG F (x) = i 4 R n H 0,(1) (κr x (y)) f (y) dy . For x ∈ R n , recall the polar coordinates centered at x by r x (y) := |y -x| , θ x (y) := y -x |y -x| .

1 ) 4 ∞

 14 |x(r, θ x ) -y|) f (y) dy dθ . (86) By Graf's addition theorem, c.f Prop 9, we haveH (1) 0 (κ r x (y)) = ∞ n=-∞ H (-n (κ|x -x|) e -inθx(x) J n (κ r x (y)) e inθx(y) , |y -x| < |x -x| ∞ n=-∞ J -n (κ|x -x|) e -inθx(x) H(1)n (κ r x (y)) e inθx(y) , |y -x| > |x -x| .In another word, we have(G F ) • ϕ Γx (θ x ) = i l=-∞ (-1) l J l (κ r) Nsrc K=1e -ilθx(s K ) H -l (κ|x -s K |) e ilθx .

Ax h -b 2 = 2 .

 22 Hn y -b 2 e 1 As a result, we are left with the minimization problem Hn y -b e 1 2

3 . 2 e 1 -y m 2 e 1 -

 32121 Compute the approximate solution x h = x 0 + Q m y m with y m minimizing |r 0 Hm Hm y 2 .

For j = 1 ,

 1 . . . , m, do (a) Compute z j = Aq j . (b) Compute w := P -1 z j .

3 .r 0 2 e 1 -y 2 . 4 .

 3124 Solving the least square problem for y m |r 0 2 e 1 -Hm y m 2 = min y ∈ Km(AP -1 , r0) Hm Compute the approximate solution x m = x 0 + P -1 Q m y m with residue Ax m -b 2 = r 0 2 e 1 -Hm y m 2 .

  V test at the Order of Approximation OrdApp test ; • and V ref at the Order of Approximation 16 ;

	Obstacle Radius Wavelength	∈ {0.0001, . . . , 0.1}
		{	1 10 p , p=1:0.1:4}
	we calculate
	•

  SL method, instead of evaluating the Hankel function exactly for 200 × 200 × 6 points in the interval [0.4, 68.96], we do a cubic Hermite interpolation on 10000 points with step size 6.86 e -3. Also see Figure9.

	Time comparison	FS-SL	Montjoie CG Q5
		BEM 3
	Size of linear system	42	435580
	Task	Duration of time (in secs)
	Construction of Coefficient matrix	1.08 e -4	1.61
	Construction of RHS	1.38 e -5	0.01
	Factorization of Coefficient Matrix	6.95 e -4	12.12
	Resolution of linear System (Direct Solver 1.03 e -4	0.35
	Total time	1.45 e -3	14.09

Inria

For FS-

  Post-processing : For FS-SL method, we have a choice of either evaluating the Hankel function exactly for 400 × 400 × 200 points in the interval [0.30, 156.93], or doing a cubic Hermite interpolation on 1000 points with step size 0.16, see alsoFigure 13 .

	Comparison at precision 1.e-3 :
	Time comparison for Pre-processing	FS-SL	Montjoie CG Q6
			BEM 2
	Size of linear system		1000		842677
	Task		Duration of time (in secs)
	Construction of Coefficient matrix	5.46 e -2	1.96
	Construction of RHS		4.34 e -5	0.0128
	Factorization of Coefficient Matrix	0.44		29.8
	Resolution of linear System (Direct Solver 2.91 e -3	0.35
	Total time		0.498	32.12
	Relative difference between Relative L 2 error
	FS-SL Ref	FS-SL 2	4.65 e -5
	FS-SL 2 Inter	FS-SL 2	1.76 e -5
	MJ Ref (Q17)	MJ Q6	6.52 e -4
	MJ Q6	FS-SL23	6.84 e -4
	FS-SL 2 Inter	MJ Q6	685 e -4
	Post-processing time	Using	Interpolate Montjoie
	Evaluation on 400 × 400 points Hankel Exact		Hankel
		26.2		4.30	0.72
	Total time (Pre + Post)	26.70		4.80	33.82

  For following comparison, the reference solution is produced by Mumps + Exact Eval. The parameters accompanying GMRES methods are (Error Tol, No of Iter Max, No Restart). See alsoFigure 18 and 17.

	Method Name		Method	Rel diff	Rel diff	No of	Pre	Post	Total
	Pre		Post		of density	in L 2		Iter	-proc	-proc	Time
	-processing		Proc		in H 1/2				time (secs) time (secs)	(secs)
	Mumps		Exact
	(-n1) Mumps		(-n16) Inter		0.0	0.0		n/a	2.51 e + 02 Number of obstacles = 2000. 9.60 e + 01 3.47 e + 02
	(-n1) Lapack		(-n16) Exact		0.0	1.34 e -05		n/a	2.51 e + 02 Radius of obstacles = 0.01 3.75 e + 01	2.89 e + 02
	(-n1) Lapack (-n1)		(-n16) Inter (-n16)	of visualization [27, 119] × 4.04 e -12 2.40 e -15 4.04 e -12 1.34 e -05	n/a n/a	7.99E + 01 Distance between two adjacent obstacles 9.60 e + 01 1.76E + 02 = 2.00 7.99E + 01 3.75 e + 01 1.18E + 02
	[27, 139]. Right LU-SGS (-n1) Exact (1.e -7, 5000, 500) (-n16)		2.89 e -04 1.20 e -07		57	Ratio Obstacle Radius Wavelength 3.75E + 01 9.60 e + 01 = 0.01 2π/10 ∼ 0.0159. 1.34E + 02
	Wavenumber κ = 10.0 ; wavelength λ = 2π/10 ∼ 0.63. Right LU-SGS (-n1) Inter (1.e -7, 5000, 500) (-n16) 2.89 e -04 1.34 e -05	57	Ratio Obstacle Radius Obs Distance = 0.01 2 = 0.005. 3.75E + 01 3.75 e + 01 7.53E + 01
	Right SGS (-n1) (1.e -7, 5000, 500) Right SGS (-n1) (1.e -7, 5000, 500)	Low definition visualization on 800 × 800 Exact (-n16) 3.50E -04 1.46 e -07 structured grid. Inter High definition visualization on 2400 × (-n16) 3.50E -04 1.34 e -05	56 56	Ratio Obstacle Distance Wavelength 3.70E + 01 9.60 e + 01 = 2 0.63 ∼ 3.18. 1.33E + 02 FS-SL method order =2. 3.70E + 01 3.75 e + 01 7.46E + 01
	Scalapack	2400 structured grid. Exact				Size matrix = 10000 × 10000.
	(-n16)		(-n16)		9.27E -12 3.87 e -15		n/a	3.49E + 01	9.60 e + 01	1.31E + 02
	Scalapack		Inter
	(-n16)		(-n16)		9.27E -12 1.34 e -05		n/a	3.49E + 01	3.75 e + 01	7.25E + 01
	Right LU-SGS (-n1)	Exact
	(1.e -9, 5000, 500)		(-n16)		3.24E -06 1.37 e -09		70	4.42E + 01	9.60 e + 01	1.42E + 02
	Right LU-SGS (-n1)	Inter
	(1.e -9, 5000, 500)		(-n16)		3.24E -06 1.34 e -05		70	4.42E + 01	3.75 e + 01	8.17E + 01
	Right SGS (-n1)		Exact
	(1.e -9, 5000, 500)		(-n16)		3.08E -06 1.31 e -09		70	4.41E + 01	9.60 e + 01	1.40E + 02
	Right SGS (-n1)		Inter
	(1.e -9, 5000, 500)		(-n16)		3.08E -06 1.34 e -05		70	4.41E + 01	3.75 e + 01	8.24E + 01

Inria

  As expected, the time cost reduces as the number of processors increases, see Figures19 and 20for a visualisation of the solution.

				High definition visualisation on 1600 ×
				1600 structured grid.
				Number of obstacles = 2000.
	[-1.05, 48.45] × [-1.05, 34.95].		Radius of obstacles = 0.03
	Wavenumber κ = 10.0 ;			Distance between two adjacent obstacles = 0.30
	Wavelength λ = 2π/10 ∼ 0.63.			Ratio Obstacle Radius Wavelength	= 0.03 2π/10 ∼ 0.048.
	Size matrix = 50000 × 50000. Low definition visualisation on 800 × 800	Ratio Obstacle Radius Obs Distance = 0.03 0.30 = 0.1.
	structured grid.			FS-SL method order =2.
	Method Name	Seq	Pre-proc	Post-proc	Total
		/ Par time (secs) time (secs) Time (secs)
	Scala + Exact (Low Def)	-n4	5.67 e + 03	1.48 e + 03	7.15 e + 03
	Scala + Inter (Low Def)	-n4	5.67 e + 03	5.62 e + 02	6.23 e + 03
	Scala + Exact (Low Def)	-n16	3.97 e + 03	4.84 e + 02	4.46 e + 03
	Scala + Inter (Low Def)	-n16	3.97 e + 03	1.82 e + 02	4.15 e + 03
	Scala + Exact (Low Def)	-n48	1.34E + 03	3.52E + 02	1.69 e + 03
	Scala + Inter (Low Def)	-n48	1.34E + 03	1.37E + 02	1.48E + 03
	Scala + Exact (High Def)	-n48	1.34E + 03	7.47E + 02	2.09E + 03
	Scala + Inter (High Def)	-n48	1.34E + 03	2.84E + 02	1.63E + 03

(among other choices of layer operators and modified ones, c.f.[START_REF] Kress | Boundary integral equations in time-harmonic acoustic scattering[END_REF])

However, this does not mean the restriction of single-layer to low-frequency problems; there have been studies of its potential in high-frequency problems, c.f.[3] and the references therein.RR n°8988

Our tests have been realized on the cluster Plafrim (www.plafrim.fr) Inria

However, we did not use a sparsified version of the preconditioner as is done in[3].RR n°8988

RR n°8988

The proof of the mapping property of Ams is given in Proposition 3, see alsoSubsection 3.4 Inria

Barucq & Chabassier & others

This corresponds to 18 points per wavelength on a domain of size 38 wavelengths × 43 wavelengths RR n°8988

However, this might change with the parallelization and optimization of the codes for the iterative solvers with SGS and LU-SGS solvers.Inria

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Case 1616 obstacles

The following experiments show that convergence is not obtained except for the family LUSGS and SGS. Even when convergence is obtained, the number of iterations increases drastically.

Parameters

Ratio Obstacle Radius

Obs Distance = 0.03 0.30 = 0.1. From Comparison 3 and 4, we see that lowering the GMRES precision improves the convergence for GMRES. Although the density error criteria indicate 10 0 error, the solution is still 'usable', expecting 10 -4 in L 2 error. However, the time cost for iterative solvers is so considerably high that for the current obstacle configuration, direct solvers are still better choices, even if one only requires low precision tolerance.

Ratio

Inria

Next we are interested in the trace of u m along ∂B(x, r) for B(x, r) that does not intersect B(x, r). We abbreviate notation by using (r, θ) and (r, θ) to denote the polar coordinates centered at x and x correspondingly,

Lemma 12.

1. Assume that x / ∈ B(x, r), for points x / ∈ B(x, r) i.e r x (x) > r, the single layer potential u m in terms of in the polar coordinates centered at x is given by

2. Assume the balls B(x, r) and B(x, r) do not intersect. The trace of u m along Γ = ∂B(x, r) is given by

3. The exterior first order trace of u m on Γ is then given by

l-m (κ |x -x|) e i(l-m)θx(x) J m (κ r) e im θ(x) .

(84)

Proof. Property 1 : For such points x / ∈ B(x, r), the value of u m in the polar coordinates centered at x is given by

The RHS of (82) is obtained by writing the factor e imθ H

m (κr) in the polar coordinates centered at x, using the addition formula Prop 9.

Property 2 : Since B(x, r) ∩ B(x, r) = ∅, for x ∈ Γ, we have r(x) = r < |x -x|. Hence, the value of u m along Γ is given by the first expression in (82). The trace of u m on Γ can now be obtained by letting r = r.

Property 3 : Since the convergence of the series is uniform on compact subsets, this allows us to take derivative in the radial variable.

B.5 Fourier Series of the trace of a general source

On [0, 2π], Fourier series inversion formula

.

Inria By (73), we have the identity e -inθx(x) = (-1) n e -inθx(x) .

As a result, we will split the integral with respect to y in (86) into two terms, one on B(x, r) and the other on the complement, The first double integral is equal to

-n (κ r)

By(x,r)

J n (κ r x (y)) e inθx(y) f (y) dy =(-1) l 2π H l,(1) (κr) r B(x,r)

J -l (κ r x (y)) e -ilθx(y) f (y) dy .

Under the assumption that Supp F ∩ B(x, r) = ∅, the above integral is zero.

The second double integral is equal to

H (1) n (κ r x (y)) e inθx(y) f (y) dy =(-1) l 2π J l (κ r)

-l (κ r x (y)) e -ilθx(y) f (y) dy .

Under the assumption that Supp F ∩ B(x, r) = ∅, the above integral is equal to

-l (κ r x (y)) e -ilθx(y) f (y) dy

-l (κ r x (y)) e -ilθx(y) f (y) dy

-l (κ r x (y)) e -ilθx(y) f (y) dy .

Point sources : Assumption on the position of sources :

Inria

C.4 Preconditioners

Intuitively, a preconditioner P is an operator such that its inverse P -1 approximates A -1 . In the algorithm of GMRES, we only need to describe the action of the inverse of the preconditioner i.e P -1 f . Denote by L, D, U the strictly lower, the diagonal and the strictly upper part of matrix A, we write

As a result, we have

Basic preconditioners :

1. P = M u is called the backward Gauss-Seidel (BGS) preconditioner.

2. P = M l is called the forward Gauss-Seidel (FGS) preconditioner .

3. P = D is the Jacobi preconditioner.

Second order Jacobi preconditioner We write

The action of P -1 can also be described as applying two Jacobi iterations; in particular, u = P -1 f is the solution to

We define

Remark 8. The above operator can be seen as a 'second-order' approximation of the Neumann series of A -1 via the splitting,

Formally, we write

Thus the second-order cut-off of the above Neumnann series gives back the operator defined in (90). The explicit definision for P is

Inria

Second order Forward Gauss-Seidel preconditioner Define P such that the action of its inverse is described as : u = P -1 f is the solution to Forward G-S :

We define

Remark 9. The above operator can be seen as a 'second-order' approximation of the Neumann series of A -1 via the splitting,

Formally, we write

Thus the second-order cut-off of the above Neumnann series gives back the operator defined in (91). The explicit definition for P is

Symmetric Gauss-Seidel (SGS) preconditioner Define P such that the action of its inverse is described as : u = P -1 f is the solution to Backward G-S :

As a result,

We define

Foward GS :

We define

C.5 Arnoldi Iteration

We define the n-th Krylov space associated with linear operator C and vector y as where Q is an orthonormal matrix. We refer readers to [START_REF] Liesen | Krylov subspace methods: principles and analysis[END_REF][Section 2.4.1] for detailed explanation of the Arnoldi process.

Intuition By writing AQ = QH we obtain the recursive formula,

Denote by Hn is the upper (n + 1) × n part of H. We also define Q n = (q 1 , . . . , q n ) ; Q n+1 = (q 1 , . . . , q n , q n+1 ) .

In terms Q n and Hn , we have AQ n = Q n+1 Hn .

Inria