Armin Lechleiter

Dinh Liem Nguyen

Spectral Volumetric Integral Equation Methods for Acoustic Medium Scattering in a Planar Homogeneous 3D-Waveguide

Scattering of acoustic waves from an inhomogeneous medium can be described by the Lippmann-Schwinger integral equation. For scattering problems in free space, Vainikko proposed a fast spectral solution method that exploits the convolution structure of this equation's integral operator by using the fast Fourier transform. In a planar 3-dimensional waveguide, the integral operator of the Lippmann-Schwinger integral equation fails to be a convolution. In this paper, we show that the separable structure of the kernel nevertheless allows to construct fast spectral collocation methods similar to Vainikko's technique. The numerical analysis of this method requires smooth material parameters; if the material parameters are, say, discontinuous, no theoretical statement on convergence is available. We show how to construct a Galerkin variant of Vainikko's method for which a rigorous convergence analysis is available even for discontinuous materials. For several distant scattering objects inside the 3-dimensional waveguide this discretization technique leads to a computational domain consisting of one large box containing all scatterers, and hence many unnecessary unknowns. However, the integral equation can be reformulated as a coupled system with unknowns defined on the different parts of the scatterer. Discretizing this coupled system by a combined spectral/multipole approach yields an efficient method for waveguide scattering from multiple objects.

Introduction

Propagation of acoustic signals inside the ocean is the basis for several modern marine technologies like SONAR (SOund Navigation And Ranging) and ocean-acoustic tomography. These techniques exploit that acoustic waves with a frequency less than a few hundred Hertz propagate inside the sea without (strong) attenuation. Consequently, low-frequency signals can propagate over huge distances along the ocean. In this paper we propose a class of spectral (i.e., Fourier-based) volumetric integral equation methods to compute such sound fields. We assume a couple of modeling assumptions on the ocean environment: The ocean has a constant height h > 0; thus, the domain of interest is a waveguide Ω = R 2 × (0, h). We further restrict ourselves to linear propagation of time-harmonic waves modeled by the Helmholtz equation ∆u

+ k 2 n 2 u = f in Ω, ( 1
)
where f is a source function with compact support, k > 0 is the constant wave number and n 2 is the refractive index. A further crucial assumption is that n 2 = 1 outside some bounded and open set D, meaning that n 2 models a local perturbation inside a homogeneous waveguide. Following [START_REF] Arens | Direct and inverse medium scattering in a 3D homogeneous planar waveguide[END_REF][START_REF] Santiago | Modified Green's functions for shallow water acoustic wave propagation[END_REF] we model the ocean-air and ocean-seabed interfaces by sound soft and sound hard boundaries, respectively, u = 0 on Γ -:= {x ∈ R 3 : x 3 = 0}, and ∂u ∂x 3 = 0 on Γ + := {x ∈ R 3 : x 3 = h}.

(

) 2
This model is reasonable for the description of underwater sound waves if the ocean depth and the water temperature are not too large. A better model would assume a layered background medium where the speed of sound depends on x 3 . The numerical methods that we develop are in principle able to incorporate such an x 3 -dependence, but this development is out of the scope of this paper.

The problem (1)-( 2) only determines the field u in a unique way if we impose additionally a radiation condition. To this end, we set

α m := π 2h (2m -1), k m := k 2 -α 2 m
for m ∈ N.

The square root in the latter definition is chosen such that k m = i α 2 m -k 2 for α m > k > 0. We assume that k m = 0 for all m ∈ N, which a kind of non-resonance condition. Let us additionally introduce the following notation for a point x in the stratified waveguide Ω,

x = (x 1 , x 2 , x 3 ) ⊤ = x x 3 with x = x 1 x 2 .
Using this notation, one can expand u by separation of variables as u(x) = ∞ m=1 sin(α m x 3 )u m (x) for |x| > R, for some R > 0 large enough such that the supports of the source f and the contrast q = n 2 -1 are contained in {|x| < R}. In view of [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] the modes u m need to satisfy ∆u m + k 2 m u m = 0 for |x| > R. Hence, as a radiation condition we impose Sommerfeld's radiation condition, lim r→∞ r 1/2 ∂u m ∂r -ik m u m = 0, where r := |x|.

There is an important connection between the source problem (1)-( 3) and scattering problems. Denote the medium's contrast by q := n 2 -1. Consider an incident wave field u i that solves the homogeneous Helmholtz equation in Ω subject to the boundary conditions [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF]. When u i hits the inhomogeneity in D there arises a scattered field u s such that the total field u = u i + u s solves ∆u + k 2 n 2 u = 0 in Ω, and both u and u s satisfy the boundary conditions [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF]. Additionally, u s satisfies the radiation conditions [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF]. The source problem (1)-( 3) hence describes the scattered field u s for the special choice f = -k 2 qu i .

Solving the waveguide scattering problem (1), ( 2) and ( 3) is equivalent to solve a volumetric integral equation of the second kind, the so-called Lippmann-Schwinger integral equation. For scattering in free space, one can exploit the convolution structure of the volume potential to construct fast spectral collocation methods [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF][START_REF]Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF][START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF][START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF] for the numerical solution of scattering problems. In contrast to scattering in free space, the volume potential defined via the waveguide Green's function G(x, y) is not a convolution; G can for instance be represented as a series of convolution operators in x weighted by certain trigonometric polynomials in x 3 . This separable structure allows to construct fast spectral integral equation methods for waveguide scattering problems via truncation and periodization of the waveguide Green's function in the lateral variables x. After periodization, a special class of trigonometric polynomials becomes the eigenfunctions of the (periodized) volume potential (see Theorem 3.5). The fast Fourier transform then allows to rapidly evaluate a spectral discretization of the potential and, using iterative solution methods, thereby yields fast methods for waveguide scattering.

If q is a smooth function, this spectral collocation method yields high-order approximations of the periodized Lippmann-Schwinger equation. However, if q is not smooth, no convergence analysis is available. To obtain rigorous convergence theory for non-smooth contrasts, we replace the spectral collocation method by a spectral Galerkin method. We show that this new discretization can still be written in an explicit discrete form, and that the method reaches optimal convergence rates in a range of Sobolev spaces W s for s > 0, whereas the collocation method can only be analyzed for s > 3/2 (see Theorem 5.4). The Galerkin method becomes more costly compared to the collocation method since it incorporates discrete Fourier transforms of larger size compared with the collocation method. Basically, the reason is that the product of two Fourier polynomials of degree l is a Fourier polynomial of degree 2l and hence discrete Fourier transforms of larger size are needed to implement this product exactly. Our numerical experiments indicate that the gain of accuracy of the Galerkin method at least corresponds to its larger memory and time demands.

We also consider a special multiple scattering problem where the penetrable scatterer is composed of several disconnected parts. Discretizing scattering problems for this geometry using the above approach yields one large computational domain containing all the scattering objects, and hence a lot of unnecessary unknowns. However, for this special setting the volumetric integral equation can be reformulated as a coupled system of integral equations where each component of the unknown is defined on one part of the scatterer. Discretization of this coupled system reduces the number of unknowns compared to the original approach. The key for the efficient discretization of the coupling terms are diagonal approximations of the waveguide Green's functions, relying on multipole expansions for complex wave numbers and uniform estimates of Bessel functions. (See Theorem 7.5 for the final spectral/multipole method.)

Alternative computational methods to solve the scattering problem (1)-( 3) include finite element methods with a non-reflecting boundary conditions, see [START_REF] Baronian | Transparent boundary conditions for the harmonic diffraction problem in an elastic waveguide[END_REF][START_REF] Zschiedrich | A new finite element realization of the perfectly matched layer method for Helmholtz scattering problems on polygonal domains in two dimensions[END_REF] for recent developments. Such methods have the advantage of local basis functions allowing for local and adaptive mesh refinement. On the other hand, 3D computations using finite element methods yield large (sparse) system matrices that require preconditioned iterative solution methods. The spectral integral approach presented here has the advantage that the (full) system matrix is never set up, preconditioners are easily constructed, and the system is solved by a simple fixed-point iteration. Since the integral equation is solved in the spectral domain the implementation of the method is easy; especially, no singular integrals have to be computed.

The structure of this paper is as follows. Section 2 provides background on the Lippmann-Schwinger integral equation. We define truncated Green's functions and corresponding integral operators in Section 3. After a reminder on trigonometric approximation in Section 4 we study spectral discretizations of volumetric integral equations in Section 5. In Section 6 we consider multiple scattering problems and diagonal approximations, yielding combined spectral/multipole methods in Section 7.

Notation: and| • | ∞ we denote the 1-the 2-and the ∞-norm on Euclidean vector spaces. For two real matrices A and B we write A ≤ B if A i,j ≤ M i,j . Following [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] we denote by J m , I m and K m the Bessel function of the first kind and the first and second modified Bessel function of order m, respectively. H

By | • | 1 , | • |,
(1) m is the Hankel function of the first kind of order m.

The Lippmann-Schwinger Integral Equation in a Waveguide

A function u that solves the Helmholtz equation ( 1), the waveguide boundary conditions [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF], and the radiation conditions (3) satisfies a volumetric integral equation of the second kind known as the Lippmann-Schwinger integral equation (see [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]). This integral equation uses the Green's function G(x, y) of the Helmholtz equation ∆u + k 2 u = 0 with constant coefficients, subject to the boundary conditions [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF] and the radiation conditions [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF]. Two explicit series representations of this Green's function are known, see [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF][START_REF] Santiago | Modified Green's functions for shallow water acoustic wave propagation[END_REF]. First, the modal representation

G(x, y) = i 2h ∞ m=1 sin(α m x 3 ) sin(α m y 3 )H (1) 0 (k m |x -ỹ|) , x = ỹ, (4)
has the advantage that it converges rapidly away from the line {x = z}, while on this line the expression is not defined. (Recall that we assumed that k m = 0 for all m ∈ N.) Second, the method of images yields an expression that is accurate near the singularity at x = z but slowly (and only conditionally) converging away from this point,

G(x, y) = 1 4π +∞ m=-∞ (-1) m    e ik|x-ym| |x -y m | - e ik ˛x-y ′ m |x -y ′ m |    , x = y, (5)
where the image source points are y m = (y 1 , y 2 , y 3 + 2mh) ⊤ and y ′ m = (y 1 , y 2 , -y 3 + 2mh) ⊤ . Formally, the volume potential V is defined by

Vf = D G(•, y)f (y) dy for f ∈ L 2 (D).
Using the series representation (5) one shows that the waveguide Green's function can be written as sum of the free-space fundamental solution Φ(x) = e ik|x| /(4π|x|) of the Helmholtz equation and an analytic function G(•, •) that solves the homogeneous Helmholtz equation in both variables, G(x, y) = Φ(x -y) + G(x, y) for x = y ∈ Ω. In consequence, the mapping properties of V are the same as those of the (free-space) volume potential with kernel Φ. From [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Chapter 6] it follows that V is a bounded operator from L 2 (D) into H 2 (B ∩ Ω) for any open ball B ⊂ R 3 . We denote this class of functions as H 2 loc (Ω). For the free-space Lippmann-Schwinger integral operator with kernel Φ and density f it is wellknown that the corresponding potential u solves ∆u + k 2 u = -f . Since G(x, y) solves the homogeneous Helmholtz equation we deduce that u = Vf ∈ H 2 loc (Ω) solves ∆u + k 2 u = -f in Ω, too. In the last equation, we understand f to be extended by zero outside of D, a convention that we will also use later on. By construction of the Green's function, the potential u also satisfies the boundary and radiation conditions (2)-( 3).

The Lippmann-Schwinger integral equation now arises in context of the scattering problem discussed in the introduction. For an incident wave u i , the scattered wave u s solves ∆u s + k 2 n 2 u s = -k 2 qu i and satisfies the boundary and radiation conditions (2) and (3). Hence ∆u s + k 2 u s = -k 2 q(u i + u s ), which means that u s -k 2 V(qu s ) = k 2 V(qu i ). Once a solution u s ∈ L 2 (D) to this equation is found, the scattered field u s in all of Ω is given by u s = k 2 V(q(u i + u s )). This is the Lippmann-Schwinger integral equation that we consider now in L 2 (D) for a right-hand side f ∈ L 2 (D),

u -k 2 V(qu) = f in L 2 (D). ( 6
)
Since the integral operator V maps L 2 (D) into H 2 (D) Riesz theory implies that existence of solution for (6) follows from uniqueness. However, in contrast to scattering in free-space, uniqueness of solution might fail for waveguide scattering problems due to resonance phenomena, see [START_REF] Morgenrother | On the instability of resonances in parallel-plane waveguides[END_REF]. In the sequel, we assume that uniqueness of solution holds, noting that the paper [START_REF] Arens | Direct and inverse medium scattering in a 3D homogeneous planar waveguide[END_REF] establishes a rather complete solution theory. In essence, uniqueness of solution always holds if the scatterer is absorbing, or under geometric non-trapping conditions on the contrast q. If none of these conditions is satisfied, analytic Fredholm theory implies that for fixed q the set of wave numbers k > 0 such that I -k 2 V(q•) has a non-trivial kernel is countable and possesses no finite accumulation point. Hence, non-uniqueness is a "rare" event, and our assumption that ( 6) is uniquely solvable in L 2 (D) for any right-hand side is reasonable.

Periodic Green's Functions and Integral Equations

Given f ∈ L 2 (D) and q ∈ L ∞ (D) we consider in this section the transformation of the Lippmann-Schwinger equation to a periodic integral equation and study properties of the transformed equation. Analysis (and, later on, computations) will be reduced to the domain

Λ ρ := {x ∈ Ω : |x| ∞ < ρ}, ρ > 0.
In analogy to our notation for points, we set Λρ = {x ∈ R 2 : |x| ∞ < ρ}. To this end, let us introduce trigonometric basis functions in L 2 (Λ ρ ). For n ∈ Z 3 + and ñ := (n 1 , n 2 ) ⊤ we set

ϕ n := 1 √ 2hρ sin(α n3 x 3 ) exp iπ ρ ñ • x , x ∈ Λ ρ , n ∈ Z 3 + := {n ∈ Z 3 : n 3 > 0}. (7)
Note that the closure of {ϕ n } n∈Z 3 + in L 2 (Λ ρ ) equals L 2 (Λ ρ ) and that ϕ n satisfies the waveguide boundary conditions [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF].

For v ∈ L 2 (Λ ρ ), v(n) := Λρ vϕ n dx , n ∈ Z 3 + ,
denotes the nth Fourier coefficient and since the functions ϕ n are orthonormal and their closure is dense,

each v ∈ L 2 (Λ ρ ) has a representation v = n∈Z 3 + v(n)v n . The basis function ϕ n is the product of v ñ(x) = 1 2ρ exp iπ ρ ñ • x and h n3 (x 3 ) = 2 h sin(α n3 x 3 ), n ∈ Z 3 + . (8)
Fractional Sobolev spaces in Λ ρ play an important role in the sequel. For s ∈ R we define W s to be the closure of the {ϕ n } n∈Z 3 + in the norm • H s (Λρ) ,

W s = {ϕ n } • H s (Λρ ) , u 2 H s (Λρ) = n∈Z 3 + 1 + |n| 2 s |û(n)| 2 .
Since differentiation becomes multiplication under the Fourier transform, one can show that for s ∈ N these spaces are subspaces of (standard) Sobolev spaces H s per of (laterally) 2ρ-periodic functions.

Let us from now on assume that the support D of the contrast q is included in

B ρ/2 = {x ∈ Ω : |x| < ρ/2}.
(We are going to see below why ρ/2 is the largest possible radius.) Then we can rewrite the Lippmann-Schwinger equation ( 6) as u -k 2 B ρ/2 G(•, y)q(y)u(y) dy = f in Λ ρ . The value of the integral in the last equation does not change if we redefine G(x, y) for |x -ỹ| > ρ, since for x ∈ B ρ/2 and y ∈ B ρ/2 the difference |x -ỹ| is always less than ρ. Motivated by this observation, let us define a function H ρ (x, k) for |x| ∞ < ρ and wave number k such that |k| > 0 and arg(k) ∈ [0, π/2],

H ρ (x, k) = H 1 0 (k|x|) if 0 < |x| < ρ, 0 else. ( 9
)
We extend this function 2ρ-periodically to R 2 by

H ρ (x + 2ρñ, k) = H ρ (x, k), for ñ ∈ Z 2 .
Using H ρ and the modal representation (4) we define a periodic Green's function G ρ (x, y) by

G ρ (x, y) = i 2h ∞ m=1 sin(α m x 3 ) sin(α m y 3 )H ρ (x -ỹ, k m ), x = ỹ ∈ Ω. ( 10
)
Note that G ρ (x, y) = G(x, y) for |x -ỹ| < ρ. In the next lemmas we provide estimates for the Fourier coefficients of t m , leading to a convergence result for the series in [START_REF] Frangi | On the application of the Fast Multipole Method to Helmholtz-like problems with complex wavenumber[END_REF]. For simplicity, we denote the series terms of G ρ by

t m (x, y) := i 2h sin(α m x 3 ) sin(α m y 3 )H ρ (x -ỹ, k m ). (11)
Lemma 3.1. Each term t m belongs to L 2 (Λ ρ × Λ ρ ) and the nth Fourier coefficient tm (n) of t m (x, •) is given by tm

(n) = iρ/ √ 2h δ m,n3 Ĥρ (n) sin(α m x 3 )v -ñ (x), n ∈ Z 3 + , with Ĥρ (n) =                      2iρ k 2 n3 ρ 2 -π 2 |ñ| 2 1 + iπ 2 2 |ñ|J 1 (π|ñ|)H (1) 0 (k n3 ρ) - iπ 2 ρk n3 J 0 (π|ñ|)H (1)
1 (k n3 ρ)

for ñ = 0, k n3 ρ = π|ñ|, 2i k 2 n3 ρ + π k n3 H (1) 1 (k n3 ρ) for ñ = 0, πρ 2 J 0 (π|ñ|)H (1)
0 (π|ñ|) + J 1 (π|ñ|)H

(1)

1 (π|ñ|) for k n3 ρ = π|ñ|. (12)
Proof. The function t m (x, y) belongs to

L 2 (Λ ρ × Λ ρ ) since H ρ (x -ỹ, k m ) has a (square-integrable) weak singularity at x = ỹ, Λρ |H ρ (x -ỹ, k m )| 2 dx = Λρ-ỹ |H ρ (x, k m )| 2 dx = Λρ |H ρ (x, k m )| 2 dx < ∞ (13)
since x → H ρ (x, k m ) is by construction 2ρ periodic in each argument. The second integral in (13) is independent of ỹ, H ρ (x -ỹ, k m ) is square-integrable in x and ỹ, and t m (x, y) belongs to L 2 (Λ ρ × Λ ρ ). We compute the Fourier coefficients tm (n) of t m (x, •),

tm (n) = i 2h Λρ sin(α m x 3 ) sin(α m y 3 )H ρ (x -ỹ, k m )ϕ n (y) dy = i √ 2h 3/2 sin(α m x 3 ) Λρ sin(α m y 3 )H ρ (x -ỹ, k m )v -ñ (ỹ) sin(α n3 y 3 ) dy = δ m,n3 i 2 √ 2h sin(α n3 x 3 ) Λρ H ρ (x -ỹ, k n3 )v -ñ (ỹ) dỹ
due to the orthogonality of the sine terms, h 0 sin(α m y 3 ) sin(α n3 y 3 ) dy = h/2 for m = n 3 and 0 else. Exploiting the periodicity of H ρ in its first argument we find that

Λρ H ρ (x -ỹ, k n3 )v -ñ (ỹ) dỹ = 1 2ρ Λρ H ρ (ỹ -x, k n3 ) exp - iπ ρ ñ • (ỹ -x) dỹ exp - iπ ρ ñ • x = 2ρ Λρ H ρ (z, k n3 )v -ñ (z) dz v -ñ (x) =: 2ρ Ĥρ (n) v -ñ (x). (14)
Hence, tm (n) = iρ/ √ 2h δ m,n3 Ĥρ (n) sin(α n3 x 3 ) v -ñ (x). In analogy to scattering problems in free space, see, e.g., [START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF]Section 3.8], the Fourier coefficients

Ĥρ (n) = 1 2ρ Λρ H ρ (z, k n3 ) exp - iπ ρ ñ • z dz = 1 2ρ |z|<ρ H (1) 0 (k n3 |z|) exp - iπ ρ ñ • z dz , n ∈ Z 3 + , (15)
of the truncated Hankel function H ρ (•, k n3 ) can be computed explicitly. Assume for a moment that

k 2 n3 ρ 2 = π 2 |ñ| 2 , such that λ n = ρ 2 /(k 2 n3 ρ 2 -π 2 |ñ| 2 ) is well-defined. Then ∆v ñ + k 2 n3 v ñ = λ -1 n v ñ and Green's second identity yields i 4 |z|<ρ H (1) 0 (k n3 |z|) v -ñ dz = iλ n 4 lim δ→0 δ<|z|<ρ H (1) 0 (k n3 |z|) (∆ + k 2 n3 )v -ñ dz = iλ n 4 |z|=ρ H (1) 0 (k n3 |z|) ∂v -ñ ∂ν -v -ñ ∂ ∂ν H (1) 0 (k n3 |z|) ds(z) -λ n i 4 lim δ→0 |z|=δ H (1) 0 (k n3 |z|) ∂v -ñ ∂ν -v -ñ ∂ ∂ν H (1)
0 (k n3 |z|) ds(z) .

For ñ = 0 we know from [28, Section 10. and for ñ = 0 we have |z|=ρ v 0 ds = π and |z|=ρ ∂v 0 /∂ν ds = 0. Moreover, it is well known (see, e.g., [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]) that i 4 lim

δ→0 |z|=δ H (1) 0 (k n3 |z|) ∂v -ñ ∂ν -v -ñ ∂ ∂ν H (1) 0 (k n3 |z|) ds(z) = v -ñ (0) = 1 2ρ . If k 2 n3 ρ 2 = π 2 |ñ| 2
we obtain the values for Ĥρ (n) given in [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF]. Finally, one uses L'Hospital's rule to compute the limiting value of the expression in the first line of [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF] for ρk n3 = π|ñ|. Lemma 3.2. For s < 1/2 the norms t m (x, •) W s are bounded by C(s)m s-3/2 for C(s) > 0 independent of m ∈ N and x ∈ Λ ρ .

Proof. Using the asymptotic expansions of Bessel and Hankel functions for large arguments in [1, Eqs. 9.2.1 & 9.2.3],

J ν (r) ∼ 2/(πr) cos(r -νπ/2 -π/4), r → ∞, H (1) ν (z) ∼ 2/(πz) exp i(z -νπ/2 -π/4) , |z| → ∞, arg(z) ∈ [0, π/2],
the growth of the coefficients Ĥρ (n) in [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF] of the periodic kernel G ρ can be bounded in terms of n.

Note that the third case of (12) only holds for a finite number of n. Assuming that the first condition in (12) holds, we find that

Ĥρ (n) ≤ C 1 + |ñ||H (1) 0 (k n3 ρ)| |J 1 (π|ñ|)| + |k n3 ||H (1) 1 (k n3 ρ)| |J 0 (π|ñ|)| α 2 n3 ρ 2 + π 2 |ñ| 2 -k 2 ρ 2 ≤ C 1 + |k n3 | -1/2 |ñ| 1/2 2 + |k n3 | 1/2 |ñ| -1/2 α 2 n3 ρ 2 + π 2 |ñ| 2 -k 2 ρ 2 ≤ C|ñ| -3/2 as |n| → ∞. (16)
Assuming that the second condition (ñ = 0) in [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF] holds, one verifies in the same way

| Ĥρ (n)| ≤ Cn -2 3 . Thus, | tm (n)| 2 ≤ ρ 2 /(2h) δ m,n3 | Ĥρ (n)| 2 ≤ Cδ m,n3 (1 + |n| 2 ) -3/2 as n → ∞,
where C is independent of m, n and x. For s < 1/2 each of the functions t m (x, •) is uniformly bounded in W s , because

t m (x, •) 2 W s = n∈Z 3 + , n3=m (1 + |n| 2 ) s | tm (n)| 2 ≤ C n∈Z 3 , n3=m (1 + |n| 2 ) s-3/2 ≤ C(s)m 2s-3 .
Proposition 3.3. The truncated Green's function G ρ belongs to W s for s < 1/2 in both variables x and y and the series in [START_REF] Frangi | On the application of the Fast Multipole Method to Helmholtz-like problems with complex wavenumber[END_REF] converges absolutely in W s for s < 1/2 as a function of x or y.

Proof. In view of the two last lemmas and the symmetry of

G ρ (x, y) = G ρ (y, x) it just remains to bound the W s -norm of G ρ (x, •): y → G ρ (x, y) W s ≤ ∞ m=1 t m W s ≤ C ∞ m=1 m s-3/2 < ∞ for s < 1/2. Remark 3.4.
The decay in x of the order 3/2 can be improved to 2 if one uses a smooth lateral cut-off of the Green's function, as it is described in [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]Page 324]. However, this procedure increases the cut-off parameter ρ and one looses the explicit knowledge of the Fourier coefficients of the kernel. Consequently, these coefficients need to be computed numerically. An efficient method to do so is described in [START_REF] Bruno | Higher-order Fourier approximation in scattering by twodimensional, inhomogeneous media[END_REF]. The scheme proposed in the latter paper can be seen as a variant of Vainikko's scheme.

For free-space scattering problems in three dimensions, a spherical cut-off of the Green's function by zero does not destroy second-order decay of the Fourier coefficients, see [START_REF]Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF][START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF]. However, the spherical cut-off used in these papers seems to be a bad choice here, in view of the explicit computations in Lemma 3.1.

The integral operator corresponding to

G ρ is V ρ f = Λρ G ρ (•, y)f dy , f ∈ L 2 (Λ ρ ). (17)
Fast methods for the solution of the Lippmann-Schwinger equation in free space [START_REF]Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF][START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF][START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF] exploit that a convolution operator becomes a multiplication operator under the Fourier transform. For our problem, V ρ lacks a convolution structure in the third coordinate, but the ϕ n still diagonalize V ρ .

Theorem 3.5. The integral operator V ρ is bounded from W s into W s+3/2 for all s ∈ R. The trigonometric basis functions ϕ n are the eigenfunctions of V ρ with corresponding eigenvalues iρ/2 Ĥρ (n

). If f ∈ L 2 (D), then V ρ f equals Vf in B ρ/2 .
Proof. We first show that V ρ diagonalizes on the basis functions ϕ n of L 2 (Λ ρ ) = W 0 and deduce boundedness of V ρ from the magnitude of the eigenvalues resulting from this computation. For fixed x ∈ Λ ρ ,

V ρ ϕ n (x) = i √ 2hh Λρ ∞ m=1 sin(α m x 3 ) sin(α m y 3 )H ρ (k n3 |x -ỹ|) sin(α n3 y 3 )v ñ(y) dy = i 2 √ 2h sin(α n3 x 3 ) Λρ H ρ (x -ỹ, k n3 )v ñ(ỹ) dỹ . (18)
In Lemma 3.3 we showed that the series defining G ρ (x, •) converges absolutely in L 2 (Λ ρ ). This validates permutation of integration and summation in [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]. In combination with ( 14) the last computation moreover implies

V ρ ϕ n = iρ √ 2h Ĥρ -ñ n3 sin(α n3 x 3 )v ñ(x) = iρ 2 Ĥρ (n)ϕ n ,
because the coefficients Ĥρ (n) depend only on the length |ñ| and n 3 , that is, Ĥρ

-ñ n3 = Ĥρ (n).
We have hence shown that V ρ diagonalizes on its eigenbasis ϕ n . The growth bound [START_REF] Lebedev | Special functions and their applications[END_REF]

shows that | Ĥρ (n)| 2 ≤ C(1 + |n| 2 ) -3/2 ,
and therefore

V ρ ϕ 2 W s = j∈Z 3 + (1 + |n| 2 ) s | (V ρ ϕ)(n)| 2 = ρ 2 4 j∈Z 3 + (1 + |n| 2 ) s | Ĥρ (n)| 2 | φ(n)| 2 ≤ C j∈Z 3 + (1 + |n| 2 ) s-3/2 | φ(n)| 2 = C ϕ 2 W s-3/2 .
This shows boundedness of V ρ from W s into W s+3/2 for s ∈ R. It remains to prove the last claim of the theorem. To this end, it is sufficient to show that for a smooth function

f ∈ C ∞ 0 (D) with compact support in D ⊂ B ρ/2 we have that V ρ (f ) equals V(f ) in B ρ/2 . Choose x ∈ B ρ/2 . Since |x -ỹ| < ρ, the definition of H ρ in (9) shows that (V ρ f )(x) = Λρ G ρ (x, y)f (y) dy = i 2h B ρ/2 ∞ m=1 sin(α m x 3 ) sin(α m y 3 )H ρ (x -ỹ, k m3 ) f (y) dy = i 2h B ρ/2 ∞ m=1 sin(α m x 3 ) sin(α m y 3 )H (1) 0 (k m3 |x -ỹ|) f (y) dy = (Vf )(x). By density of C ∞ 0 (D) in L 2 (D), V ρ equals V on B ρ/2 for any density f ∈ L 2 (D).
Due to the equality of the integral operators V and V ρ for densities supported in D stated in the last theorem, we can reformulate the Lippmann-Schwinger equation

u -k 2 V(qu)| D = f in L 2 (D) (19) as u -k 2 V ρ (qu)| D = f . Extending f by zero and u by k 2 V ρ (qu) to Λ ρ yields a function v ∈ L 2 (Λ ρ ) solving v -k 2 V ρ (qv) = f in L 2 (Λ ρ ). ( 20
) Proposition 3.6. Let f ∈ L 2 (D) and extend f by zero to Λ ρ . If u ∈ L 2 (D) solves (19) then v = V(qv) + f ∈ L 2 (Λ ρ ) solves (20). Furthermore, the restriction of a solution v ∈ L 2 (Λ ρ ) of (20) to D yields a function u in L 2 (D) that solves (19).
In the beginning of this paper we were interested to find the scattered field u s for a waveguide scattering problem with incident field u i . Setting f = k 2 V ρ (qu i ) yields a solution v to (20) that equals u s on D. To evaluate u s in the complement Ω \ D one uses the Lippmann-Schwinger integral equation (with smooth kernel!) u s (x) = k 2 D G(x, y)q(y)(u s (y) + u i (y)) dy for x ∈ D.

Trigonometric Projection and Interpolation

Spectral discretization of the periodic Lippmann-Schwinger equation ( 20) relies on the Fourier transform. Set

Z 3 l := {n ∈ Z 3 : -l j < n j ≤ l j , j = 1, 2, 1 ≤ n 3 ≤ l 3 } for l = (l 1 , l 2 , l 3 ) ⊤ ∈ N 3 .

The discrete subspace of trigonometric polynomials

T l := span ϕ n : n ∈ Z 3 l ⊂ L 2 (Λ ρ ) of dimension L := 4 l 1 l 2 l 3
will serve as approximation space for a solution to [START_REF] Nåsell | Inequalities for modified Bessel functions[END_REF]. We denote the orthonormal projection from L 2 (Λ ρ ) onto T l by

P l : L 2 (Λ ρ ) → T l ⊂ L 2 (Λ ρ ), P l v = n∈Z 3 l v(n)v n , where v(n) = Λρ vϕ -n dx
denotes the nth Fourier coefficient of v. Note here that Lemma 3.5 implies that the integral operator V ρ and the projection P l commute: V ρ (P l u) = P l V ρ (u) for u ∈ L 2 (Λ ρ ). Later on, we rely on the following classical approximation properties.

Lemma 4.1. For u ∈ W s and r ≤ s it holds that u -P l u W r ≤ min(l) r-s u W s , where min(l) = min{l 1 , l 2 , l 3 }. Proof. For u ∈ W s , u -P l u 2 W r = n ∈Z 3 l (1 + |n| 2 ) s (1 + |n| 2 ) r-s |û(n)| 2 ≤ 1 + min(l) 2 r-s u 2 H s (Λρ) .
The proof of the following lemma is technical and therefore presented in Appendix A. It uses the grid points

x (l) n = ρ n 1 l 1 , ρ n 2 l 2 , h n 3 l 3 ⊤ , n = (n 1 , n 2 , n 3 ) ⊤ ∈ Z 3 l . ( 21
) Lemma 4.2. There is a basis {ϕ * n } n∈Z 3 l of T l such that ϕ * n (x (l) j ) = δn, j for n, j ∈ Z 3 l .
In consequence, any function in u l ∈ T l can be uniquely represented by its grid values at the grid points x

(l) n , u l = n∈Z 3 l u l (x (l) n )ϕ *
n . Additionally, we can define an interpolation operator

Q l : C 0 (Λ ρ ) → T l , Q l (f ) = n∈Z 3 l f (x (l) n )ϕ * n ,
mapping a continuous function in Λ ρ to the unique trigonometric interpolation polynomial satisfying interpolation conditions at the grid points x (l)

n . We indicate the following approximation property without proof (see Theorem 8.3.1 in [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF] for an analogous result).

Lemma 4.3. For u ∈ W s with s > 3/2 it holds that u -Q l u W r ≤ C min(l) r-s u W s for 0 ≤ r ≤ s.
The transformation mapping the grid values v l (x (l) j ) of a function v l in T l to the Fourier coefficients vl of v l is called the discrete Fourier transform. Due to the structure of the basis functions ϕ l , this transform is the product of a two-dimensional discrete Fourier transform in the lateral variables, see [25, Section 10.5.4], and a discrete sine transform in the vertical variable (a type-3 discrete sine transform, see [START_REF] Martucci | Symmetric convolution and the discrete sine and cosine transforms[END_REF]). Combining formulas for the 2D discrete Fourier transform and the type-3 discrete sine transform, one arrives in the following formula connecting grid values v l (x

(l) j ) of v l ∈ T l with the Fourier coefficients vl (n) of v l , v(n) = ρ l 1 l 2 l 3 l1 j1=-l1+1 l2 j2=-l2+1 l3-1 j3=1 v l (x (l) j ) sin πj 3 (2n 3 -1) 2l 3 exp - iπ ρ j • xn + (-1) n3+1 ρ 2l 1 l 2 l 3 l1 j1=-l1+1 l2 j2=-l2+1 v l (x (l) j ) exp - iπ ρ j • xn . (22)
Define C l := {(c(j)) j∈Z 3 l , c(j) ∈ C}. The operator that maps grid values (v(j))

j∈Z 3 l ∈ C l to the Fourier coefficients vl = (v l (j)) j∈Z 3 l ∈ C l of the unique trigonometric interpolation polynomial v l ∈ T l that satisfies v l (x (l) j ) = v(j) is in the following denoted by S l : C l → C l , (v(j)) j∈Z 3 l → (v l (j)) j∈Z 3 l . (23)
It is well known that the transform S l can be computed in O L log(L) operations (recall: L = 4l 1 l 2 l 3 ), due to the fast Fourier transform, see [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF], and its variant, the fast sine transform.

Later on, we will once require a second Fourier transform where the sine transform in the vertical variable is replaced by a cosine transform. To distinguish this new transform from the above Fourier(sine) transform we call it the Fourier-cosine transform. The Fourier-cosine series of a function

v ∈ L 2 (Λ ρ ) is v(x) = 1 √ 2hρ n∈Z 3 ,n3≥0 v(n)ϕ ñ(x) cos π h n 3 x 3 , x ∈ Λ ρ , (24)
with Fourier-cosine coefficients v

(n) = ( √ 2hρ) -1 Λρ v(x)ϕ ñ (x) cos π h n 3 x 3 dx for n 3 = 0 and v(n) = ( √ 8hρ) -1 Λρ v(x)ϕ ñ(x)
dx for n 3 = 0. Again, it is well known that any function v ∈ L 2 (Λ ρ ) can be represented by its Fourier-cosine series which is norm-convergent in L 2 (Λ ρ ). The corresponding trigonometric polynomials are

v l (x) = 1 √ 2hρ l1 n1=-l1+1 l2 n2=-l2+1 l3-1 n3=0 v(n)ϕ ñ(x) cos π h n 3 x 3 , x ∈ Λ ρ . (25)
Since the summation over n 3 starts at 0, it is natural to introduce

Z 3 l,0 := {n ∈ Z 3 : -l j < n j ≤ l j , j = 1, 2, 0 ≤ n 3 ≤ l 3 -1}. Define vl = (v(j)) j∈Z 3 l,0 ∈ C l 0 := {(c(j)) j∈Z 3 l,0 , c(j) ∈ C}. (26)
For a polynomial v l of the form [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF] with Fourier-cosine coefficients vl , and m ∈ Z 3 + , we denote by C l,m : C l 0 → C m the operator that maps the Fourier-cosine coefficients vl to the values of the polynomial at the grid points (x (m) j

) j∈Z 3 m . The only case interesting for us will be when m > l. Roughly speaking, this operation can then be expressed as extension by zero of the Fourier-cosine coefficients, followed by a discrete inverse Fourier(-cosine) transform, see [START_REF] Strang | The discrete cosine transform[END_REF]. For m > l and M := 4m 1 m 2 m 3 , one can therefore evaluate C l,m in O(M log(M )) operations. However, since we never use the discrete inverse Fourier(-cosine) transform explicitly, we do neither give explicit formulas for this transform nor for C l,m .

Spectral Discretization and Error Estimates

The original spectral discretization approach for the Lippmann-Schwinger equation ( 20) in [START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF] is to multiply (20) by the contrast q and to solve for the new unknown u = qv in spaces of trigonometric polynomials. Whenever q is not smooth, solving for qv yields worse convergence rates for the discretized problem than solving for v, since multiplication by q destroys smoothness. Consequently, let us try to keep the contrast q inside the integral operator V ρ and solve for v instead of qv. A corresponding variant has been chosen, for different reasons, in [START_REF]Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF].

Assumption 5.1. Throughout this section we assume that the contrast q is compactly supported in B ρ/2 and belongs at least to P C(Λ ρ ), the class of piecewise continuous and everywhere defined functions in Λ ρ . By definition, for f ∈ P C(Λ ρ ) there are pairwise disjoint Lipschitz domains U j , j = 1, . . . , N with

Λ ρ = ∪ N j=1 U j and f | Uj is continuous.
The crucial point of any spectral discretization of ( 20) is to map the product qu l back into the space of trigonometric polynomials T l . Let us first use the interpolation operator Q l for this task and seek

u l ∈ T l such that u l -k 2 V ρ Q l (qu l ) = Q l f. ( 27
)
The discrete Fourier transform S l from Section 4 allows to state this problem matrix-vector form. To this end, we denote element-wise multiplication of two elements a and b in

C l by a • b, (a • b)(n) = a(n)b(n) for n ∈ Z 3 l .
Further, we abbreviate the Fourier coefficients of the integral kernel (multiplied by k 2 ) by

kl (j) = k 2 iρ 2 Ĥρ (j), j ∈ Z 3 l .
By q l ∈ C l and f l ∈ C l we denote the point values of q and f at the points x [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF]. (Obviously, one crucial assumption for the collocation method [START_REF] Strang | The discrete cosine transform[END_REF] is that these point values are welldefined.) Then the scheme ( 27) can be written in matrix-vector form in terms of the Fourier coefficients ûl ∈ C l of u l using the discrete Fourier transform S l from [START_REF] Pinsky | Fourier series of radial functions in several variables[END_REF],

(l) n , n ∈ Z 3 l , from
ûl -kl • S l q l • S -1 l (û l ) = S l f l . (28)
Error estimates for this scheme will only be available for q, f ∈ W s for s > 3/2. We can also discretize [START_REF] Nåsell | Inequalities for modified Bessel functions[END_REF] by applying the projection P l to all terms of the equation: Find

u l ∈ T l such that u l -k 2 V ρ P l (qu l ) = P l f. ( 29
)
Writing down the fully discrete form of this discrete problem requires some preparation. Let us denote the restriction of ûl

∈ C l to C m , m < l for j = 1, 2, 3, by R l,m , R l,m (û l ) = vm , vm (n) = ûl (n) for n ∈ Z 3 m . The extension operator E m,l extends ûm ∈ C m by zero to C l , m < l, E m,l (û m ) = vl , vl (n) = ûm (n) for n ∈ Z 3 m , 0 else. If q ∈ P C(Λ) ⊂ L 2 (Λ ρ
), then we can develop q into a Fourier-cosine series as in [START_REF] Santiago | Modified Green's functions for shallow water acoustic wave propagation[END_REF] that converges in L 2 (Λ ρ ). (Developping q into a Fourier(-sine) series as in [START_REF] Bruno | Higher-order Fourier approximation in scattering by twodimensional, inhomogeneous media[END_REF] does not allow to prove an analogue to the crucial Lemma (5.2) below.) For l ∈ Z 3 + we define a polynomial q l by

q l (x) = 1 2hρ 2 l1 m1=-l1+1 l2 m2=-l2+1 l3-1 m3=0 q(m)ϕ m(x) cos π h n 3 x 3 , l ∈ Z 3 + . ( 30
)
We extend the Fourier-cosine coefficients q(m) to all m ∈ Z 3 by setting q( m m3 ) = 0 for m 3 < 0. As in [START_REF] Sauter | Randelementmethoden[END_REF], we set ql = (q(m)) m∈Z 3 l,0 ∈ C l 0 . Recall the mapping C l,m : C l 0 → C m introduced in the end of Section 4.

Lemma 5.2. An equivalent fully discrete form of the projection method (29) is given by

ûl -kl • R 3l,l S 3l (C 2l,3l q2l ) • S -1 3l E l,3l ûl = fl . (31) Proof. For n ∈ Z 3 + , qu l (n) = Λρ qu l ϕ n dx = j∈Z 3 l ûl (j) Λρ qϕ n ϕ j dx = 1 2hρ 2 j∈Z 3 l ûl (j) m∈Z 3 ,m3≥0 q(m) h 0 cos π h m 3 x 3 sin(α n3 x 3 ) sin(α j3 x 3 ) dx 3 δ j-ñ, m.
Since 2 sin(α n3 x 3 ) sin(α j3 x 3 ) = cos((α n3 -α j3 )x 3 ) -cos((α n3 + α j3 )x 3 ) the latter integral reduces to

h 0 cos π h m 3 x 3 sin(α n3 x 3 ) sin(α j3 x 3 ) dx 3 = 1 2 h 0 cos π h m 3 x 3 cos π h (n 3 -j 3 )x 3 -cos π h (n 3 + j 3 -1)x 3 dx 3 = h 4 (δ m3,n3-j3 -δ m3,n3+j3-1 ) , m 3 ≥ 0, h 2 (δ m3,n3-j3 -δ m3,n3+j3-1 ) , m 3 = 0, and qu l (n) = 1 8ρ 2 j∈Z 3 l ûl (j)(1 + δ m3,0 ) q j-ñ n3-j3 -q j-ñ n3+j3-1 . (32)
Consequently, the Fourier-cosine coefficients q(m) for m ∈ Z 3 \ Z 3 2l,0 do not influence the computation of P l (q u l ). Using the definition of q 2l in (30) we find that P l (q u l ) = P l (q 2l u l ). Moreover, (32) with q replaced by q 2l also shows that q 2l u l (n) vanishes for n ∈ Z 3 + \ Z 3 3l . Therefore the Fourier coefficients of q 2l u l are given by S 3l applied to the grid values of this trigonometric function at the grid points (x

(3l) n ) n∈Z 3 3l .
For ûl , the grid values are given by S -1 3l E l,3l ûl . By construction, the grid values of q 2l at the grid points x (3l) n are C 2l,3l (q 2l ). Hence,

( q 2l u l )(j) = S 3l C 2l,3l (q 2l ) • S -1 3l E l,3l ûl , j ∈ Z 3 3l . Since V ρ (P l (qu l )
) is a Fourier polynomial, we can compute its Fourier coefficients S l V ρ (P l (qu l )) from the Fourier coefficients ûl of u l ∈ T l by

S l V ρ (P l (qu l )) = kl • P l (qu l ) = kl • P l (q 2l u l ) = kl • R 3l,l S 3l C 2l,3l (q 2l ) • (S -1
3l E l,3l ûl ) . Applying S l to (29), we obtain the fully discrete matrix-vector formulation that was announced in (31).

The advantage of the Galerkin (or projection) method (31) over the collocation method [START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF] is that convergence can even be shown for discontinuous q. The drawback is that the Galerkin variant requires larger memory space since it relies on Fourier transforms of higher dimension; however, our numerical experiments later on indicate that the method is also more accurate than its collocation counterpart. Further, it might be difficult to compute the Fourier coefficients q(n) of q l for complicated material configurations; see Remarks 5.6 and 5.7 concerning this issue.

Proposition 5.3. Assume that q ∈ P C(Λ ρ ) is such that the multiplication u → qu is continuous on W t for some t ≥ 0, qu W t ≤ C u W t for all u ∈ W t . Then there is a constant C independent of u and l such that V ρ (qu) -V ρ P l (qu) W s ≤ C min(l) s-t-3/2 u W t for s ≤ t + 3/2, and, if

t > 3/2, V ρ (qu) -V ρ Q l (qu) W s ≤ C min(l) s-t-3/2 u W t for 3/2 ≤ s ≤ t + 3/2. Proof. We estimate V ρ (qu) -V ρ P l (qu) W s ≤ C qu -P l (qu) W s-3/2 ≤ C min(l) s-t-3/2 qu W t ≤ C min(l) s-t-3/2 u W t for s -3/2 ≤ t. For t ≥ 3/2, s ≥ 3/2, and s -3/2 ≤ t, V ρ (qu) -V ρ Q l (qu) W s ≤ C qu -Q l (qu) W s-3/2 ≤ C min(l) s-t-3/2 qu W t ≤ C min(l) s-t-3/2 u W t .
Theorem 5.4. Assume that the Lippmann-Schwinger equation ( 6) is uniquely solvable in L 2 (Λ ρ ). Then (20) has a unique solution u ∈ L 2 (Λ ρ ) for any right-hand side f ∈ L 2 (Λ ρ ).

(a) Assume that f ∈ W t , t ≥ 0, and that q is such that the multiplication u → qu is continuous on W r for 0 ≤ r ≤ max(t -3/2, 0). Then for l ∈ Z 3 with min(l) large enough there is a unique solution u l ∈ T l of (29) and

u l -u W s ≤ C min(l) s-t f W t for 0 ≤ s ≤ t.
(b) Assume that f ∈ W t for t > 3/2 and that q is such that the multiplication u → qu is continuous on W r for 3/2 < r ≤ t. Then for l ∈ Z 3 with min(l) large enough there is a unique solution u l ∈ T l of [START_REF] Strang | The discrete cosine transform[END_REF] and

u l -u W s ≤ C min(l) s-t f W t for 0 ≤ s ≤ t.
Proof. Uniqueness of solution for the Lippmann-Schwinger equation ( 6) implies uniqueness of solution for the periodized version (6) due to Proposition (3.6). The assumption that u → qu is continuous on W r implies by the (periodized) Lippmann-Schwinger equation that for a right-hand side in W t the solution u ∈ L 2 (Λ ρ ) belongs to W t , too. E.g., under the assumptions of part (a) and for t ≤ 3/2,

u W t ≤ C qu L 2 (Λρ) + f W t ≤ C u L 2 (Λρ) + f W t ,
and for larger t one uses a bootstrap argument.

(a) For the projection method (29) we can apply a standard convergence results for such methods. For instance, under our assumptions [26, Satz 4.2.9] yields u l -u W s ≤ C min v l ∈T l v l -u W s and consequently

u l -u W s ≤ C min v l ∈T l v l -u W s ≤ C min(l) s-t u W t ≤ C min(l) s-t f W t .
(b) Compactness of V ρ on W s , the Fredholm alternative, boundedness of the multiplication operator u → qu on W r for the given range of r and unique solvability of (20

) in L 2 (Λ ρ ) imply that I -k 2 V ρ (q •) is an invertible operator on W s for 0 ≤ s ≤ t. Choose 3/2 < s ≤ t. Due to Proposition 5.3 I -k 2 V ρ (q •) -[I -k 2 V ρ Q l (q •)] W s →W s ≤ C min(l) -3/2 ,
and hence a Neumann series argument shows that I -k 2 V ρ Q l (q •) is invertible on W s and that (I -k 2 V ρ Q l (q •)) -1 is uniformly bounded for min(l) large enough. We apply I -k 2 V ρ Q l (q •) to the difference u l -u, where u l ∈ T l is the solution of the discrete problem and u is the solution of the continuous problem [START_REF] Nåsell | Inequalities for modified Bessel functions[END_REF], and find that

(I -k 2 V ρ Q l (q •))(u l -u) = (Q l -I)f + k 2 V ρ (Q l -I)(qu). ( 33
)
The uniform boundedness of (

I -k 2 V ρ Q l (q •)) -1 on W s yields u l -u W s ≤ C (Q l -I)f + V ρ (Q l -I)(qu) W s ≤ C min(l) s-t f W t + V ρ W s →W s (Q l -I)(qu) W s ≤ C min(l) s-t f W t + C V ρ W s →W s min(l) t-s qu W t , 3/2 < s ≤ t. (34)
Reformulating the left-hand side of (33), we obtain that

[I -k 2 V ρ (q•)](u l -u) -k 2 V ρ (Q l -I)(qu l -qu) = (Q l -I)f + k 2 V ρ (Q l -I)(qu)
and hence the following estimate for all s ∈ [0, t],

u l -u W s ≤ [I -k 2 V ρ (q•)] -1 W s →W s k 2 V ρ W s →W s (Q l -I)(qu l -qu) W s + (Q l -I)f W s + V ρ W s →W s (Q l -I)qu W s ≤ C min(l) t-s qu l -qu W t + f W t + qu W t ≤ C min(l) t-s f W t , 0 ≤ s ≤ t. (35)
The last inequality follows from the boundedness of the multiplication by q on W t and the bound

u l -u W t ≤ C u W t ≤ C f W t which is a consequence of (34).
Next we illustrate the above result by discussing smoothness assumptions for q. Let C r,1 (Λ ρ ), r ∈ N, be the space of all Lipschitz continuous functions on Λ ρ such that all partial derivatives up to order r are still Lipschitz continuous.

Example 5.5. We suppose that the general assumptions of Theorem 5.3 hold.

(a) If f ∈ W t , 0 ≤ t ≤ 3/2, and q ∈ P C(Λ ρ ), then Theorem 5.3(a) applies: Indeed, for q ∈ P C(Λ ρ ) the multiplication u → qu is bounded on L 2 (Λ ρ ) due to the Cauchy-Schwarz inequality.

(b) If f ∈ W t , t ≥ 0, and q ∈ C r,1 (Λ ρ ), r ∈ N such that r ≥ max(t -3/2, 0), then Theorem 5.3(a) applies: Multiplication by q is shown to be bounded on L 2 (Λ ρ ) as in part (a); thus, we assume now that t -3/2 > 0. Theorem 3.20 in [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF] shows that qu

W s-3/2 ≤ C q C r,1 u W s-3/2 for 0 ≤ s ≤ max(t -3/2, 0).
For f ∈ W t , t ≥ 3/2, and q ∈ C r,1 (Λ ρ ), r ∈ N such that r ≥ t, the same argument shows that Theorem 5.3(b) applies.

(c) If f, q ∈ W t , t > 3/2, then Theorem 5.3(b) applies: The multiplication by q is bounded on W s since there is

C > 0 such that uv W s ≤ C u W s v W s for u, v ∈ W s , s > 3/2.
This can be shown as in [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]Lemma 5.13.1], where the corresponding result in one dimension is proven.

If f ∈ W t for t > 3 then it is even sufficient to assume that q ∈ W t-3/2 to obtain the estimate u l -u W s ≤ C min(l) s-t f W t for 0 ≤ s ≤ t as in Theorem 5.3(b). To this end, we employ the inequality qu W t-3/2 ≤ C q W t-3/2 u W t from [13, Proof of Lemma 1]. Indeed, under these assumptions, the inverse of I -k 2 V ρ Q l (q•) is still uniformly bounded on W s for 3/2 < s ≤ t and (35) is replaced by

u l -u W s ≤ [I -k 2 V ρ (q•)] -1 W s →W s k 2 V ρ W s-3/2 →W s (Q l -I)(qu l -qu) W s-3/2 + (Q l -I)f W s + V ρ W s-3/2 →W s (Q l -I)qu W s-3/2 ≤ C min(l) t-s qu l -qu W t-3/2 + f W t + qu W t-3/2 ≤ C min(l) t-s f W t , 0 ≤ s ≤ t.
Consider again the scattering problem of finding the scattered field u s for an incident wave u i . The Lippmann-Schwinger equation describing the scattered field is u -k 2 V ρ (qu) = k 2 V ρ (qu i ). Since u i is smooth, the smoothness of the right-hand side depends only the smoothness of q. For q ∈ P C(Λ ρ ), V ρ (qu i ) belongs, e.g., to W 3/2 , and we obtain that u l -u W s ≤ C min(l) s-3/2 u i L 2 (Λρ) for the Galerkin approximation where s ≤ 3/2. If q ∈ C r,1 (Λ ρ ), then V ρ (qu i ) belongs to W r+3/2 and u l -u W s ≤ C min(l) s-r-3/2 u i L 2 (Λρ) , again for the Galerkin approximation u l . For r ≥ 1 the same holds for u l solving the collocation approximation.

The discrete systems ( 27) and ( 29) can be preconditioned by the inverse of the discrete system for a coarser discretization, which leads to a two-grid method. This idea is well-known for integral equations of the second kind and has been worked out in detail for Vainikko's method in [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF][START_REF]Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF][START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF]. For the above collocation and projection methods, preconditioners can be constructed in the very same way. Let us write the collocation discretization [START_REF] Strang | The discrete cosine transform[END_REF] or the Galerkin discretization (29) as

u l -K l u l = f l in T l , l ∈ Z 3 l , (36)
where K l is either k 2 V ρ Q l (q l u l ) or k 2 V ρ P l (q 2l u l ). For m < l, we multiply this equation by (I l -K m ) -1 . (I l is the identity on T l , K m is extended by zero from T m onto T l , and we assume that the inverse exists.)

Note that (I

l -K m ) -1 (K l -K m )u l = u l -(I l -K m ) -1 (I l -K l )u l .
This implies that u l , solution to (36), solves the fixed point equation

u l = (I l -K m ) -1 [(K l -K m )u l + f l ] .
If m < l is large enough, this is a fixed-point equation for a contracting operator and the successive approximations

u (0) l = (I l -K m ) -1 f l , u (j) l = (I l -K m ) -1 (K l -K m )u (j-1) l + f l , j = 1, 2, . . . (37)
converge to the solution u l to (36). Since, in Fourier space, K m only acts on Fourier coefficients with index in Z 3 m , (I l -K m ) -1 is best computed spectrally (see [START_REF] Hohage | On the numerical solution of a three-dimensional inverse medium scattering problem[END_REF] for details). For instance, the fixed-point iteration (37) formulated for the collocation method (28) reads û(j)

l = I l -E m,l R l,m + E m,l I m -km • S m [q l • S -1 m R l,m ] -1 kl • S l [q l • S -1 l û(j-1) l ] -E m,l km • S m [q l • (S -1 m R l,m û(j-1) l )] + S l f l . ( 38
)
The advantage of the fixed point iteration over solving (36) directly by an iterative solver is of course that (37) requires, at each iteration, the inversion a smaller linear system than (36). Therefore the memory requirements are considerably smaller than solving (36) directly by an iterative solver. For instance, on the fine grid the scheme (38) requires merely to evaluate matrix-vector products of the form kl

• S l [q l • S -1 l û(j-1)
l ] which requires O(L log(L)) (recall: L = 4l 1 l 2 l 3 ) operations using the FFT. The iteration requires to store the four vectors kl , q l , û(j-1) l and û(j) l and thus 16l 1 l 2 l 3 complex variables. We finish this section with a numerical experiment that confirms the convergence rates of Example 5.5(c). For these experiments, h = 1/2, ρ = 1/2, and the wave number equals k = 12.5 which corresponds to two propagating modes in the waveguide. We introduce the family of contrasts q α : Λ ρ → R for α ≥ 0 by

q α (x) = 1/16 -|x| 2 -(x 3 -1/2) 2 α if |x| 2 + (x 3 -1/2) 2 < 1/16, 0 else. ( 39
)
One can show that the function q α belongs to W s for all 0 ≤ s < α + 1/2. For better comparability of the results for different values of α we normalize the contrast q α to have L 2 (Ω)-norm equal to one in all computations.

Remark 5.6. The Fourier-cosine coefficients qα (n) can be computed (semi-)explicitly since q α is a radial function (with respect to (0, 0, 1/2) ⊤ ). Indeed, by rewriting the basis function cos(π/h n 3 x 3 ) as (exp(iπ/h n 3 x 3 ) + exp(-iπ/h n 3 x 3 ))/2 one finds a representation of qα (n) as an integral of a radial function times an exponential, and the formulas in [START_REF] Pinsky | Fourier series of radial functions in several variables[END_REF] (see the proof of the first Lemma of Section 6) allows to write qα (n) in terms of a one-dimensional integral that can either be computed explicitly (e.g., for α an integer) or approximated numerically.

The test problem that we consider is to compute the scattered field for an incident point source G(•, p) with source point p = (-2, 0, 1/4) ⊤ . We approximate the solution in T l where l = (2 n , 2 n , 2 n ) ⊤ for n = 2, . . . , 6. (The corresponding grid is uniform in all three directions with step width 2 -n .) For all computations, we use the two-grid scheme explained above, and we stop the fixed-point iteration of the scheme when the relative residual is less than 10 -8 . The preconditioner on the coarse grid is computed on T m where m = (2 ⌊n/2⌋ , 2 ⌊n/2⌋ , 2 ⌊n/2⌋ ) using GMRES and we also stop the GMRES iteration when the relative residual is less than 10 -8 . The choice of the discretization parameters as powers of two is not crucial, but it is sufficient for our purpose of checking convergence rates. The reference solution is computed using the collocation method for n = 8 and we iterate until the relative residual is less than 10 -10 . This computations was done on a workstation with 4 processors and 12 GB RAM using MATLAB. The computation of the reference solution just fitted into this RAM and the fixed-point iteration converged after about 360 seconds (independent of α). Note also that the Fourier transforms in MATLAB are computed using the FFTW package, and MATLAB executes these transforms in parallel on the four processors.

Figure 1 shows that the discrete L 2 -error of the projection method fits quite well to the theoretical convergence estimates of Example 5.5(c). (The theoretical convergence rates are indicated by a dasheddotted line.) For α = 3/2, q α belongs to W s for s < 2 and Example 5.5(c) indicates a convergence order 7/2 for the collocation method. For this α the numerical error indeed fits well to this predicted convergence rate. For α = 0, 1/2, 1 there is no convergence theory for the collocation method, and Figure 1 shows that the collocation scheme does not reach the convergence rate min(l) α+2 of the projection method. Table 1 indicates the run time for the experiments of Figure 1 in seconds. The run times for the projection method are about six to eight times larger than for the collocation method (for n = 6 quotients are between 5.8 and 6.2). On the other hand, the quotient of the relative L 2 -error of the collocation method and the projection method is, at least for m = 6, always larger than 8. This shows that the projection method has its own interest not only for the purpose of mathematical analysis: the more expensive numerical scheme leads to a corresponding gain in numerical precision, at least if the approximation space has high dimension. α=0 α=0.5 α=1 α=1.5 n=2 0.04 0.04 0.17 0.17 n=3 0.17 0.17 Remark 5.7. Since it is rather involved to obtain precise approximations for the Fourier coefficients q for general q, we also tested a simplification of the projection method (31) where C 2l,3l (q 2l ) is replaced by q evaluated at the grid points x (3l) j . For the above class of contrasts q α this simplification did not perturb the convergence rate of the method for α = 0, 1/2, 1, 3/2. However, we do not know whether it is possible to prove this observation rigorously.

Distributed Scatterers, a System of Integral Operators, and Diagonal Approximation of the Green's Function

In the next two sections we construct a variant of Vainikko's method designed for a scatterer that splits up into several disconnected parts. Our aim is to avoid to compute solutions for this geometry in a large box containing all parts of the obstacle. To this end, we develop a mixed spectral/multipole method generalizing the above technique in that the computational domains is the union of several boxes containing each one of the disconnected parts of the scatterer. In this section we first derive a reformulation of the Lippmann-Schwinger equation as a coupled system. Afterwards, we construct a diagonal approximation of the waveguide Green's function that relies on multipole expansions (see [START_REF] Greengard | A fast algorithm for particle simulations[END_REF]) for Hankel functions with complex arguments. We briefly derive these expansions and prove error estimates, since we did not find suitable multipole expansions for complex wave numbers in the literature (note, however, that the recent work [START_REF] Frangi | On the application of the Fast Multipole Method to Helmholtz-like problems with complex wavenumber[END_REF] provides numerical experiments and suggests parameter choices for multipole expansions with complex wave numbers in 3D). In the next section we prove convergence results and error estimates for the described extension of Vainikko's method. We assume that the contrast q = n 2 -1 splits up into J ∈ N parts,

q = q 1 + q 2 + • • • + q J ,
where the functions q j have disjoint support D j := supp(q j ) in Ω. Strengthening this assumption, we even suppose that there is o j ∈ R 3 and ρ j > 0 such that (recall:

B ρ = {|x| < ρ}) D j is compactly contained in the cylinder B (j) := o j + B ρj /2 , D j ⊂ B (j) = o j + B ρj /2 , j = 1, . . . , J,
and that the closures B (j) are mutually disjoint. Set

δ min = inf xi∈B (i) , xj ∈B (j) , 1≤i =j≤J |x i -x j | and δ max = sup xi∈B (i) , xj ∈B (j) , 1≤i =j≤J |x i -x j |.
We further introduce computational domains Λ (j) := o j + Λ ρj that compactly contain the domains D j . The Lippmann-Schwinger integral equation u -k 2 V(qu) = f posed in L 2 (D), D := ∪ J j=1 D j , is now decomposed into a system of equations. We split u = u 1 + • • •+ u J where supp(u j ) ⊂ D j for j = 1, . . . , J. Since the supports D j are disjoint, this splitting is unique. Set u = (u 1 , . . . , u J ) ⊤ , f = (f 1 , . . . , f J ) ⊤ ∈ ⊕ J j=1 L 2 (D j ) and q = (q 1 , . . . , q J ) ⊤ ∈ ⊕ J j=1 L ∞ (D j ). We denote the element wise multiplication of such vectors again by q • u, more precisely, q • u = (q 1 u 1 , . . . , q J u J ) ⊤ . Further, we define integral operators

V ij : L 2 (Λ (j) ) → L 2 (Λ (i) ) by V ij u = Λ (j)
G(•, y)u(y) dy

Λ (i) , i, j = 1, . . . , J. (40)
If u ∈ L 2 (D) solves ( 6), then u = (u 1 , . . . , u J ) ∈ ⊕ J j=1 L 2 (D j ) solves the integral equation

u -k 2 V 11 ••• V 1J . . . . . . . . . V J1 ••• V JJ (q • u) = f in ⊕ J j=1 L 2 (D j ). ( 41
)
Moreover, if u solves the latter vector-valued problem, then u ∈ L 2 (D) defined by u| Dj = u j solves (6).

We will now reformulate this equation in spaces of periodic functions and afterwards consider numerical schemes to solve it. Denote the fractional Sobolev spaces W s from now on as W s (Λ ρj ); further, W s = ⊕ J j=1 W s (Λ ρj ) for s ∈ R. The norm on W s defined as u W s = J j=1 u j W s (Λρ j ) . We also set

W = W 0 = ⊕ J j=1 L 2 (Λ ρj ). Since D j ⊂ B (j) there exist cut-off functions χ j ∈ C ∞ 0 (Λ (j)
) such that χ j = 1 on the support D j of q j , χ j = 0 in Λ (j) \ B (j) and 0 ≤ χ j ≤ 1. We use these cut-off functions to define truncation operators on Λ (j) that additionally shift the function into Λ ρj , see Figure 2,

T + j :L 2 (Λ ρj ) → L 2 (Λ (j) ), (T + j u)(x) = (χ j u)(x -o j ), x ∈ Λ (j) , j = 1, . . . , J, T - j :L 2 (Λ (j) ) → L 2 (Λ ρj ), (T - j v)(x) = (χ j v)(x + o j ), x ∈ Λ ρj , j = 1, . . . , J.
The component-wise action of T - j on u ∈ ⊕ J j=1 L 2 (Λ (j) ) is denoted as T ± u. As in Section (3), we consider the periodic kernels G ρj and the associated integral operators V ρj , and we abbreviate

K ij = T - i V ij T + j for i = j, K =   Vρ 1 ••• K 1J . . . . . . . . . K J1 ••• Vρ J   , (42)
that is, K has diagonal entries V ρj and off-diagonal entries = o 2 + B ρ2/2 are contained in the boxes Λ (1) and Λ (2) , respectively. The computational domain Λ ρ1 is the reference domain for Λ (1) and T ± 1 transport functions between the two domains.

K ij . 0 Λ ρ1 T - 1 T + 1 D 1 D 2 B (1) B (2) Λ (1) Λ (2)
Assumption 6.1. Throughout the rest of the paper we assume that the contrast q = (q 1 , . . . , q J ) belongs at least to P C := ⊕ J j=1 P C(Λ (j) ) and that each q j is compactly supported in

B (j) = o j + B ρj /2 . Theorem 6.2. Let f ∈ W . Then any solution v ∈ W to v -k 2 K(T -(q) • v) = f in W (43)
gives rise to a solution u ∈ ⊕ J j=1 L 2 (D j ) to (41) with right-hand side (T + f j ) J j=1 by setting u j = T + j (v j ) Dj , j = 1, . . . , J. Any solution u to (41) with right-hand side (T + f j ) J j=1 yields a solution v to (43

) by setting v = k 2 K(T -q • T -u) + f .
We omit the proof of Theorem 6.2, since it follows directly from Proposition 3.6. Note that Theorem 6.2 implies that the regularity theory for the Lippmann-Schwinger equation carries over to (43).

Our aim is to discretize (43) on product spaces of trigonometric polynomials. Discretization of the diagonal terms is done precisely as in Section 5. Discretization of the off-diagonal terms requires multipole expansions for the waveguide Green's function G in (4). Therefore we truncate the series (4) and investigate approximate diagonalization of Hankel functions, to finally arrive in an exponentially convergent diagonal approximation of G consisting only of a finite number of terms.

By m * we denote subsequently the first integer such that

k m = (k 2 -α 2 m ) 1/2 is purely imaginary. Lemma 6.3. There is a constant C = C(δ min , m * ) > 0 independent of 1 ≤ i = j ≤ J such that for all x ∈ B (i) and y ∈ B (j) G(x, y) - i 2h M m=1 sin(α m x 3 ) sin(α m y 3 )H (1) 0 (k m |x -ỹ|) ≤ Ce -πδminM/h , M ≥ m * . (44)
Proof. First, we use [8, Lemma 2.2], stating that |H

ν (z)| ≤ exp(-Im (z)(1 -Θ 2 /|z| 2 ) 1/2 )|H (1)
ν (Θ)| for ν ∈ N 0 , Re (z) ≥ 0, Im (z) ≥ 0 and 0 < Θ ≤ |z|, to obtain that

H (1) ν (k m |x -ỹ|) ≤ e -Im (km) |x-ỹ| " 1-˛km * km ˛2« 1/2 H (1) ν (|k m * ||x -ỹ|) , m ≥ m * , ν ∈ N 0 . (45)
Using the integral formula [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF]Pg. 441]

|H (1) ν (r)| 2 = J 2 ν (r) + Y 2 ν (r) = 8 π 2 ∞ 0 K 0 (2r sinh(t)) cosh(2νt) dt , ν ≥ 0, r > 0, (46)
and the integral representation [1, Eq. 9.6.24] for the modified Bessel function K 0 , it follows that |H

ν (r)| is monotonically decreasing in r > 0 for ν ∈ N 0 . For 1 ≤ i = j ≤ J, x ∈ B (i) and y ∈ B (j) , the difference x -ỹ is bounded from below by δ min > 0. Since

exp -Im (k m ) |x -ỹ|(1 -|k m * /k m | 2 ) 1/2 ≤ exp -Im (k m ) δ min (1 -|k m * /k m |) ≤ exp(-Im (k m ) δ min ) exp(δ min |k m * |)
we obtain H

(1)

0 (k m |x -ỹ|) ≤ C exp(-Im (k m )δ min ) H (1) 0 (|k m * |δ min ) for m ≥ m * where C = e δmin |k m * | . Moreover, Im (k m ) = (α 2 m -k 2 ) 1/2 ≥ α m -k ≥ πm/h -π/(2h) -k for m ≥ m * . Now we are ready to prove the claim, ∞ m=M+1 sin(α m x 3 ) sin(α m y 3 )H (1) 0 (k m |x -ỹ|) ≤ ∞ m=M+1 H (1) 0 (k m |x -ỹ|) ≤ C ∞ m=M+1 e -δmin(αm-k) ≤ C ∞ m=M+1 e -δmin( π(2m-1) 2h -k) ≤ C exp δ min k + π 2h 1 -exp(-πδ min /h) e -πδminM/h . (47)
For fixed m ∈ N, the Hankel functions H

0 (k m |x -ỹ|) can be approximately diagonalized. For x we write its cylindrical coordinates as x = (r x cos(ϕ x ), r x sin(ϕ x ))

⊤ .

Choose N ∈ N, x ∈ B (i) and y ∈ B (j) , i = j. Further, let a ∈ [-N, N ] such that a ≡ m (mod 2N + 1),

f ± n (x, k) = e ±ikrx cos( 2πn 2N +1 -ϕx) , and s n (x, k) = 1 2N + 1 N j=-N (-i) j H (1) j (kr x )e ij(ϕx-2πn 2N +1 ) . (48)
We introduce short-hand notation for cylindrical coordinates of some special vectors. For 1

≤ 1 = j ≤ J, o i -o j = r ij cos(ϕ ij ), r ij sin(ϕ ij ), 0 ⊤ ,
and for x ∈ B (i) and y ∈ B (i) we define functions r(x, y) and ϕ(x, y) by

(ỹ -o j ) -(x -o i ) = r(x, y) cos(ϕ(x, y)), r(x, y) sin(ϕ(x, y)), 0 ⊤ . (49)
Strictly speaking, r(x, y) = r ij (x, y) but since we always consider x ∈ Λ (i) and y ∈ Λ (j) in the following it will not cause confusion to suppress this dependence (also for ϕ(x, y)).

A multipole expansion for the two-dimensional fundamental solution to the Helmholtz equation with real wave number has been presented and analyzed in [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF][START_REF] Banjai | Hierarchical matrix techniques for low and high frequency Helmholtz equation[END_REF]. However, we did not find such expansions for complex wave numbers in the literature. Proposition 6.4. For x ∈ B (i) , y ∈ B (j) , 1 ≤ i = j ≤ J, and m ∈ N it holds that

H (1) 0 (k m |x -ỹ|) = 2N +1 n=1 f - n (x -o i , k m )s n (o i -o j , k m )f + n (ỹ -o j , k m ) + |n|>N J n (k m r(x, y))e -inϕ(x,y) H (1) n (k m r ij )e inϕij + i n-a H (1) a (k m r ij )e iaϕij . ( 50
)
Remark 6.5. In the next lemmas we show that the second line of (50) is small, and hence the first line yields an approximation of H

0 (k m |x -ỹ|) where x and ỹ are decoupled.

Proof. The wave number k in [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF] is restricted to be positive, however, inspecting the proof given in the latter paper one can verify that the expansion (50) extends to wave numbers k ∈ C, and thus especially to the modal wave numbers k m . Essentially, this follows from checking that the addition theorem and the integral identities used to derive (50) extends to complex wave numbers. To this end, one employs the following version of the addition theorem from [16, Ch. 5.12] which holds for all k ∈ C, H

0 (k|z|) = ∞ m=-∞ H (1) m (k|x|)J m (k|ỹ|)e im(θx-θỹ) , z = x + ỹ, |x| > |ỹ|. (1)
Consequently, splitting x -y = o i -o j + r(x, y)(cos(ϕ(x, y)), sin(ϕ(x, y))) ⊤ and noting that o i -

o j = |o i -o j |(cos(ϕ ij ), sin(ϕ ij )) ⊤ , we find that H (1) 0 (k m |x -ỹ|) = ∞ m=-∞ H (1) m (k m r ij )J m (k m r(x, y))e im(ϕij -ϕ(x,y)) . (51)
This identity replaces equation (3.2) in the proof of [3, Theorem 3.1] and allows to extend (50) to complex wave numbers. The rest of the proof is quite similar to the proof given in [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF]. One notes that

J m (k m r(x, y)) = 1 2πi m 2π 0 e i(kmr(x,y) cos(t)-mt) dt = 1 2πi m 2π 0
e ikmr(x,y) cos(t-ϕ(x,y)) e -im(t-ϕ(x,y)) dt by a change of variables t → t -ϕ(x, y). Then one exploits that

r(x, y) cos(t -ϕ(x, y)) = r(x, y) cos(t) cos(ϕ(x, y)) + sin(t) sin(ϕ(x, y)) = (x -o i -y + o j ) • cos(t) sin(t)
to get that

J m (k m r(x, y))e -imϕ(x,y) = 1 2πi m 2π 0 e ikm(x-oi-y+oj )• cos(t)
sin(t) e -imt dt .

Now one discretizes the integral appearing in the latter equation using the trapezoidal rule with 2N + 1 quadrature points and replaces the term J m (k m r(x, y)) exp(-imϕ(x, y)) in (51) by the resulting expression. The computation of the resulting error can be done precisely as in [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF] and yields the remainder term appearing in the second line of (50).

Lemma 6.6. Let m ∈ N and assume that there exists a constant η ∈ (0, 1/2) such that r(x, y) < ηr ij for all x ∈ B (i) and y ∈ B (j) , 1 ≤ i = j ≤ J. Then there is C = C(δ max , η) and

N 0 = N 0 (k m , δ max ) such that for 1 ≤ i = j ≤ J, x ∈ B (i) , and y ∈ B (j) , |n|>N J n (k m r(x, y))e -inϕ(x,y) H (1) n (k m r ij )e inϕij + i n-a H (1) a (k m r ij )e iaϕij ≤ Cη N , N ≥ N 0 .
Proof. For real k m the estimate for the remainder term that needs to be shown follows from [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF][START_REF] Banjai | Hierarchical matrix techniques for low and high frequency Helmholtz equation[END_REF]. However, we need to indicate how to extend this estimate to complex wave numbers k m for m > m * Due to πH

n (k m |x-ỹ|) = -2ie -nπi/2 K n (|k m ||x-ỹ|) and since the modified Bessel function K n is monotonic in n for real arguments due to [1, Eq. 9.6.24], we can estimate |n|>N J n (k m r(x, y))e -inϕ(x,y) H (1) n

(k m r ij )e inϕij + i n-a H (1) a (k m r ij )e iaϕij ≤ 2 |n|>N |J n (k m r(x, y))| |H (1) n (k m r ij )|, m > m * . (52)
Moreover, since J n (k m r(x, y)) = e nπi/2 I n (|k m |r(x, y)), and since I n and K n are positive and real for positive orders and real arguments we also have

|n|>N J n (k m r(x, y))e -inϕ(x,y) H (1) n (k m r ij )e inϕij + i n-a H (1) a (k m r ij )e iaϕij ≤ 4 π |n|>N I n (|k m |r(x, y)) I n (|k m |r ij ), m ≥ m * .
Due to [20, Proposition 1] we know that

0 < I n+1 (r) I n (r) ≤ r r + n , n ≥ 0, r > 0,
and [15, Theorem 2.1] states that

0 < K n+1 (r) K n (r) ≤ (n + 1) + (n + 1) 2 + r 2 r , n ≥ 0, r > 0.
Hence, because r(x, y) ≤ ηr ij ,

I n+1 (|k m |r(x, y))K n+1 (|k m |r ij ) I n (|k m |r(x, y))K n (|k m |r ij ) ≤ |k m |r(x, y) |k m |r(x, y) + n n + 1 + (n + 1) 2 + |k m | 2 r 2 ij |k m |r ij ≤ η n + 1 + (n + 1) 2 + |k m | 2 r 2 ij |k m |r(x, y) + n ≤ η   1 + 1 n + 1 + 1 n 2 + |k m | 2 δ max n   ≤ η   1 + 1 N + 1 + 1 + 1 n 2 + |k m | 2 δ max N + 1   .
Since we assumed that 0 < η < 1/2, we can choose N 0 = N 0 (k m , δ max ) so large that

η N := η   1 + 1 N + 1 + 1 + 1 N 2 + |k m | 2 δ max N + 1   < 1 for N ≥ N 0 .
With this choice of N 0 , for any N ≥ N 0 the series

n>N I n (|k m |r(x, y))K n (|k m |r ij ) ≤ I N +1 (|k m |r(x, y))K N +1 (|k m |r ij ) n>N η n-N -1 N ≤ 1 1 -η N I N +1 (|k m |r(x, y))K N +1 (|k m |r ij ) (53)
converges. From the asymptotic expansion of the Bessel functions I N and K N for large orders [1, Eq. 9.3.1] we infer that

I N +1 (|k m |r(x, y))K N +1 (|k m |r ij ) ≤ |J N +1 (|k m |r(x, y))| |H (1) N +1 (|k m |r ij )| ≤ C N + 1 e|k m |r(x, y) 2N + 2 N +1 e|k m |r ij 2N + 2 -N -1 = C N + 1 r(x, y) r ij N +1 ≤ C π η N +1 N + 1 . (54)
Plugging in (54) into (53) and combining it with (52) yields the result.

Combining Lemmas 6.3 and 6.6 yields the following error estimate for the approximation

G M,N (x, y) = i 2h M m=1 sin(α m x 3 ) sin(α m y 3 ) 2N +1 n=1 f - n (x -o i , k m )s n (o i -o j , k m )f + n (ỹ -o j , k m ) (55)
for x ∈ Λ (i) and y ∈ Λ (j) to the waveguide Green's function G.

Proposition 6.7. Assume that there exists a constant η ∈ (0, 1/2) such that r(x, y) < ηr ij for all x ∈ B (i) and y ∈ B (j) , 1 ≤ i = j ≤ J. Then there is C = C(δ min , δ max , m * , η) such that for 1 ≤ i = j ≤ J, x ∈ B (i) , and y ∈ B (j) , it holds

|G(x, y) -G M,N (x, y)| ≤ C(e -πδminM/h + M η N ) M ≥ m * , N ≥ N 0 (k M , δ max ).
Proof. By Lemma 6.3 we arrive in an exponentially convergent approximation to the waveguide's Green's function G that consists of M terms. The M Hankel functions in this approximation are then diagonalized using Lemma 6.6, introducing an error of the order of M η M .

We also need a corresponding error estimate for partial derivatives of G(x, y)-G M,N (x, y). To this end we use multiindex notation and denote the length of a multiindex β = (β 1 , β 2 , β 3 ) ⊤ by |β| 1 = β 1 +β 2 +β 3 . Proposition 6.8. Assume that there exists a constant η ∈ (0, 1/2) such that r(x, y) < ηr ij for all x ∈ B (i) and y ∈ B (j) ,

1 ≤ i = j ≤ J. Denote by β = (β 1 , β 2 , β 3 ) ∈ N 2 0 a multiindex with β 1 + β 2 ≤ 2, β 3 ∈ N 0 .
Then there is a constant C = C(δ min , δ max , η, m * , |β| 1 ) such that for 1 ≤ i = j ≤ J, x ∈ B (i) , and y ∈ B (j) , it holds

∂ |β|1 ∂x β G(x, y) -G M,N (x, y) ≤ C(e -πδminM/h + M η N ) M ≥ m * , N ≥ N 0 (k M , δ max ).
Proof. We proceed as for the estimate in Proposition 6.7 and first truncate the series defining the partial derivative of the Green's function G. Afterwards we use a multipole expansion for the Hankel functions appearing in the truncated series.

For 1 ≤ i = j ≤ J, x ∈ B (i) , and y ∈ B (j) we have that |x -ỹ| > δ min and therefore the series in (4) converges absolutely and uniformly,

∂ |β|1 ∂x β G(x, y) = i 2h ∞ m=1 sin(α m y 3 ) ∂ |β|1 ∂x β sin(α m x 3 )H (1) 0 (k m |x -ỹ|) , |x -ỹ| > δ min , (56)
and we need to obtain bounds for the remainder of the truncated series. Obviously, partial derivatives with respect to x 3 have no effect other than creating factors |α m | ≤ Cm. Again, we use (45) and monotonicity of |H Finally, we also need to estimate partial derivatives of the remainder term appearing in the second line of (50), for all m = 1, . . . , M . Indeed, after truncating the representation of the waveguide Green's function G in (56), we rely on the multipole expansion (50) to diagonalize the first M terms of the Green's function. We need to estimate |n|>N ∂ | β|1 ∂x β J n (k m r(x, y))e -inϕ(x,y) H (1) n (k m r ij )e inϕij + i n-a H (1) a (k m r ij )e iaϕij (57) for m = 1, . . . , M . We mimic the proof of [3, Lemma 3.2] and proceed by bounding partial derivatives of J n (kr(x, y))e -inϕ(x,y) in terms of higher-order Bessel functions. The definitions of r(x, y) and ϕ(x, y) in (49) imply that |∂r(x, y)/∂x 1,2 | ≤ 1, |∂ϕ(x, y)/∂x 1,2 | ≤ 1/r(x, y), and

∂ | β|1 ∂x β r(x, y) ≤ 2 r(x, y) , ∂ | β|1 ∂x β ϕ(x, y) ≤ 3 r(x, y) 2 for | β| ≤ 2, r(x, y) > 0.
Following the proof of [START_REF] Amini | Analysis of a diagonal form of the fast multipole algorithm for scattering theory[END_REF]Lemma 3.2], we arrive at

∂ | β|1 ∂x β J n (kr(x, y))e -inϕ(x,y) ≤ 3|k| 2 J n-1 (|k|r(x, y)), | β| ≤ 1, and
∂ | β|1 ∂x β J n (kr(x, y))e -inϕ(x,y) ≤ 4|k| 2 J n-2 (|k|r(x, y)), | β| ≤ 2, n ≥ 3.
These bounds imply the given error estimates for the remainder term in (57) by the same techniques employed in the proof of Lemma 6.6.

A Combined Spectral/Multipole Method for Multiple Scattering

In this section, we consider a fully discrete approximation to the system of integral equations (43) and prove convergence of the discrete solution to the solution of (43). The system of integral equations ( 43) is equivalent the (periodized) Lippmann-Schwinger integral equation [START_REF] Nåsell | Inequalities for modified Bessel functions[END_REF] considered in Section 5. However, if the scatterer contains multiple components, discretizing the system (43) results in a much smaller linear system then discretizing directly [START_REF] Nåsell | Inequalities for modified Bessel functions[END_REF].

Using the same notation as in the last section, we define V ij M,N to be the volume integral operator with kernel G M,N (see (55)),

V ij M,N u = Λ (j)
G M,N (•, y)u(y) dy

Λ (i) , i = j = 1, . . . , J,
and by

K ij M,N : L 2 (Λ ρj ) → L 2 (Λ ρi ) an associated operator defined by K ij M,N = T - i V ij M,N T + j .
We first prove the following operator approximation result. Proposition 7.1. Let 1 ≤ i = j ≤ J and assume that there exists a constant η ∈ (0, 1/2) such that r(x, y) < ηr ij for all x ∈ B (i) and y ∈ B (j) , 1 ≤ i = j ≤ J. Then

(K ij -K ij M,N )ϕ W s+2 (Λρ i ) ≤ C(e -πδminM/h + M η N ) ϕ W s (Λρ j ) (58)
for all ϕ ∈ W s (Λ ρj ), s ≥ 0, M ≥ m * and N ≥ N 0 (k M , δ max ) large enough.

Proof. Due to real interpolation theory for fractional Sobolev spaces (see, e.g., [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]Lemma 5.12.2]) it is sufficient to prove (58) for s ∈ N 0 . For s = 0, this estimate follows from Proposition 6.8 and the Cauchy-Schwartz inequality. For s a positive integer we make use of an equivalent norm on W s (Λ ρj ).

Namely, there are constants c 1 > 0 and c 2 > 0 such that

c 1 ϕ 2 W s (Λρ j ) ≤ |β|1≤s ∂ |β|1 u ∂x β 2 L 2 (Λρ j ) ≤ c 2 ϕ 2 W s (Λρ j ) , s ∈ N.
Moreover, we make use of the following property of G and G M,N : For β ∈ N 2 0 ,

∂ | β|1 ∂ x β G(x, y) = (-1) β ∂ | β|1 ∂ ỹ β G(x, y), and
∂ | β|1 ∂ x β G M,N (x, y) = (-1) β ∂ | β|1 ∂ ỹ β G M,N (x, y). (59)
Both relations hold since G(x, y) and G M,N (x, y) only depend on the difference x -ỹ. Let us denote

K(x, y) = G(x, y) -G M,N (x, y)
for the rest of this proof and assume that β ∈ N 3 0 is a multiindex of length less than or equal to s + 2. By (59),

∂ |β|1 ∂x β (T - i (V ij -V ij M,N )T + j ϕ) = ∂ |β|1 ∂x β χ i (x + o i )((V ij -V ij M,N )T + j ϕ)(x + o j ) = γ≤β β γ ∂ |β-γ|1 ∂x β-γ χ i (x + o i ) Λ (j) ∂ |γ|1 ∂x γ K(x, y) (χ j ϕ)(y -o j ) dy = γ≤β (-1) |γ| β γ ∂ |β-γ|1 ∂x β-γ χ i (x + o i ) Λ (j) ∂ γ3 ∂x γ3 3 ∂ |γ|1 ∂ ỹγ K(x, y) (χ j ϕ)(y -o j ) dy .
Now we integrate by parts up to s -2 times to transport derivatives with respect to ỹ onto ϕ. Hence, split γ = γ 1 + γ 2 into two multiindices where

|γ 1 | 1 ≤ 2 and |γ 2 | 1 ≤ s -2. Since χ j ∈ C ∞ 0 (Λ ρj )
vanishes on the boundary of Λ ρj , there arise no boundary terms from the partial integration,

Λ (j) ∂ γ3 ∂x γ3 3 ∂ |γ|1 ∂ ỹγ K(x, y) (χ j ϕ)(y -o j ) dy = (-1) |γ| Λ (j) ∂ γ3 ∂x γ3 3 ∂ |γ1|1 ∂ xγ1 K(x, y) ∂ |γ2|1 ∂ ỹγ2 (χ j ϕ)(y -o j ) dy ,
and using the Cauchy-Schwarz inequality and Proposition 6.8 we estimate

Λ (j) ∂ |γ|1 ∂x γ K(•, y) (χ j ϕ)(y -o j ) dy 2 L 2 (Λ (i) ) ≤ C Λ (i) Λ (j) ∂ |γ1|1 ∂ xγ1 K(x, y) 2 dy dx ϕ 2 W s-2 (Λρ j ) ≤ C(e -πδminM/h + M η N ) 2 ϕ 2 W s-2 (Λρ j ) for M ≥ m * and N ≥ N 0 (k M , δ max ). Consequently, ∂ |β|1 ∂ xβ (T - i (V ij -V ij M,N )T + j ϕ) L 2 (Λρ i ) ≤ C(e -πδminM/h + M η N ) ϕ W s-2 (Λρ j )
where C is independent of M ≥ m * , N ≥ N 0 (k M , δ max ), and ϕ ∈ W s-2 (Λ ρj ).

The last lemma yields discrete schemes for the spectral approximation of the off-diagonal terms T - i V ij T + j occurring in the operator K in (42). We discretize these operators by replacing the kernel G by G M,N and apply projection operators to the separable parts of G M,N depending on x and y. (Since there is no danger of confusion we do not yet denote the dependence of the projection operators on the domain Λ ρj explicitly.) This procedure yields finite-dimensional operators that we denote by

K ij M,N,l1,l2 : L 2 (Λ ρj ) → L 2 (Λ ρi ), (K ij M,N,l1,l2 ϕ)(x) = i 2h M m=1 2N +1 n=1 Q l1 χ i (x) sin(α m x 3 )f - n (x, k m ) (x) s n (o i -o j , k m ) Λ (j) Q l2 χ j (y -o j ) sin(α m y 3 )f + n (ỹ, k m ) (y) ϕ(y -o j ) dy (60) for x ∈ Λ ρi , ϕ ∈ L 2 (Λ ρj ), l 1,2 ∈ N 3 and 1 ≤ i = j ≤ J. Recall that we defined K ij = T - i V ij T + j in (42).
Proposition 7.2. Under the assumptions of Proposition 7.1,

K ij ϕ -K ij M,N,l1,l2 ϕ W s+2 (Λρ i ) ≤ C(e -πδminM/h + M η N + C 1 (N )C 2 (M )(min(l 1 , l 2 )) s+2-r ) ϕ W s (Λρ j ) (61) for all ϕ ∈ W s (Λ ρj ), s ≥ 0, M > m * , N ≥ N 0 (K m , δ min ), r > s + 2, and l 1,2 ∈ N 3 .
Remark 7.3. The proof shows that the constant C 1 (N ) grows exponentially in N , see (62). The behavior of C 2 (M ) for large M depends on the sign of max 1≤i =j≤J (ρ i + ρ j -δ min ). If this quantity is negative, then C 2 (M ) ∼ M 2r , if it is positive, then C 2 (M ) grows exponentially in M .

Proof. The estimate to be shown is based on Lemma (4.3), stating that u -Q l u W s+2 (Λρ j ) ≤ C min(l) s+2-r u W r (Λρ j ) for u ∈ W r (Λ ρj ) with r > 3/2 and 0 ≤ s + 2 ≤ r. Moreover, the error estimates for K ij M,N,l1,l2 depend on the smoothness of the kernel G M,N (x, y) in both variables x and y. The expansion (50) shows that the smoothness of G M,N (x, y) directly depends on the smoothness of the functions f ± n defined in (48),

f ± n (x, k m ) = e ±ikmrx cos( 2πn 2N +1 -ϕx) = e ±ikmrx[sin( 2πn 2N +1 ) cos(ϕx)-cos( 2πn 2N +1 ) sin(ϕx)] = e ±ikm[sin( 2πn 2N +1 )x1-cos( 2πn 2N +1 )x2] .
Hence, f ± n (x, k) is a plane wave with wave number k m , and hence especially a smooth function of x = (x 1 , x 2 ) ⊤ . Moreover, the partial derivatives of f ± n (•, k m ) are pointwise bounded by

∂ |β|1 ∂ xβ f ± n (x, k m ) ≤ |k m | β e Im (km)ρj , x ∈ Λ ρj , β ∈ N 2 0 .
Because |k m | grows like m as m → ∞ we have that χ j sin(α m x 3 )f ± n (x, k m ) W s (Λρ j ) ≤ Cm s e Im (km)ρj for s ∈ N 0 , where C is independent of j = 1, . . . , J, m ∈ N, and n ∈ Z. As in (45) one deduces that |H The norm of χ i (x)G M,N (x, y)χ j (x) in the space W r (Λ ρj ) × W r (Λ ρj ) is hence bounded by

χ i (x)G M,N (x, y)χ j (y) W r (Λρ j )×W r (Λρ j ) ≤ M m=1 2N +1 n=1 χ i sin(α m x 3 )f - n (x, k m ) W r (Λρ i ) |s n (o i -o j , k m )| χ j sin(α m y 3 )f + n (ỹ, k m ) W r (Λρ j ) ≤ C M m=1 m 2r e (ρi+ρj -δmin) Im (km) 2N +1 n=1 H (1) n (k min δ min ) (62) ≤ C(2N + 1) H (1) 2N +1 (k min δ min ) M m=1 m 2r e (ρi+ρj -δmin) Im (km) ≤ C (2N + 1) H (1) 2N +1 (k min δ min ) =:C1(N ) M 2r m * -1 + M m=m * e (ρi+ρj -δmin)(α 2 m -k 2 ) 1/2 =:C2(M)
.

The constant C 1 (N ) grows exponentially in N . The behavior of C 2 (M ) depends on the sign of

ρ i + ρ j - δ min . Recall that α m -k ≤ (α 2 m -k 2 ) 1/2 ≤ α m = π/(2h)(2m -1). If ρ i + ρ j -δ min ≥ 0 then C 2 (
M ) grows exponentially in M , too. More interesting is the case, ρ i + ρ j -δ min < 0, where

M m=m * e (ρi+ρj -δmin)(α 2 m -k 2 ) 1/2 ≤ e -k(ρi+ρj -δmin) M m=m * e π/h(ρi+ρj -δmin)m ≤ e (πm * /h-k)(ρi+ρj -δmin)
1 -e π/h(ρi+ρj -δmin) . for 0 ≤ s ≤ t, and an analogous estimate holds for P l replaced by Q l . Hence, by choosing l = l(M, N ) large enough, one can again use a Neumann series argument as in the proof of Theorem 5.4 to show that the discrete system (63) has a solution u l in T l for M , N , and min(l ♭ ) and min(l) large enough. The claimed error estimate for this approximate solution follows as in the proof of Theorem (5.3)(b).

Remark 7.6. (a) The proof shows that l ♭ can be chosen merely in dependence of M and N , but independent of l; especially, l ♭ can be chosen significantly smaller than l, yielding efficient discretizations in (63).

(b) Examples for smoothness assumptions on q can be constructed as in Example 5.5.

(c) Of course, one can also let the multipole parameters M and N depend on the domains Λ (j) .

In a first numerical experiment we test the accuracy of the multipole expansions by extending a solution to the discrete problem ( 27) from a uniform grid with step width 1/32 of the computational domain

[-1/4, 1/4] × [-1/4, n + 1/4] × [0, 1/2]
to the corresponding grid of the domain [n -1/4, n + 1/4] × [n -1/4, n + 1/4] × [0, 1/2], n = 1, 2, 4 or 6. The corresponding values of the multipole parameter η are 1/2, 1/4, 1/8 and 1/12, and the discretization parameter l hence equals l = (2 5 , 2 5 , 2 5 ) ⊤ . The wave number equals 17.5 and hence the first five modal wave numbers (rounded to four digits) are 17.22, 14.75, 7.714, 13.32i, and 22.21i. Table 2(a)-(b) shows the relative discrete L 2 error between the multipole expansion given in (60) (see (66) for the discrete form) and the exact extension of the solution using the waveguide Green's function. For η = 1/4 choosing M larger than three (corresponding to the number of propagating modes) does not improve the accuracy for all values of N . For η = 1/2 the domains between which we extend are closer. Consequently, incorporating the evanescent fourth mode into the multipole expansion improves the accuracy of the method. For M = 3, the minimal error is about 7e-09 while for M = 4 it is about 6e-10. Depending on the magnitude of the first complex wave number k m * this effect is more or less pronounced. If we set k = 20.8, the first complex modal wave numbers is k 4 = 7.14i. Since the magnitude of the first complex wave number is smaller than for k = 17.5, the corresponding mode is significant for the expansion error. If N is large enough, Table 2 n (k m r x ) in the multipole expansion (see (48)) causes cancellation and thus to a reduction of numerical accuracy due to numerical instability. This is reflected in the exponential growth of the constant C 1 (N ) appearing in Proposition (7.2) and already observable in Table 2. The accuracy of the multipole expansion for the Green's function of the Helmholtz equation in two dimensions has been investigated in [START_REF] Ohnuki | Truncation error analysis of multipole expansions[END_REF]. Define For a Dirichlet scattering problem with a scatterer of diameter d, the authors of [START_REF] Ohnuki | Truncation error analysis of multipole expansions[END_REF] claim that 10 -u(N ) is an estimate for the truncation error of the expansion (50), that is, u(N ) indicates the number of correct digits of the expansion in (50). Due to numerical instability caused by the exponential growth of the Bessel function H [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] n (k m r x ) for large n, an implementation of the multipole method can however not obtain arbitrary accuracy. In [START_REF] Ohnuki | Truncation error analysis of multipole expansions[END_REF], the authors claim that l(N ) is an upper bound for the number of correct digits that an implementation of the multipole method can reach in double precision.

In our experiments on multipole expansions for medium scattering we found that the parameter d should be taken to be q L 2 (Ω) , the L 2 norm of the contrast q. With this choice, the relative discrete L 2 error of the expansion is pretty well approximated by 10 -min(u(N ),l(N )) . Figure 3 In a second experiment we check the convergence rates given in Theorem 7.5 for scattering from two inhomogeneities. Again, we use the contrasts q α from (39). For o 1 = (0, 0, 0) ⊤ and o 2 = (1, 1, 0) ⊤ , the contrast q = (q 1 , q 2 ) is defined by q j (• -o j ) : Λ (j) = o j + Λ 1/2 → R, q j (• -o j ) = q α (•), j = 1, 2. The wave number k = 12.5, the height h = 1/2 and the discretization parameters are the same as in the numerical experiment for a single scatterer in the end of Section 5. The parameters M and N equal 4 and 22. To solve the linear system arising from the integral equation we employ a two-grid schemes as described in Section 5 with the same discretization in each reference domain as in Section 5, and also with the stopping criterion. The reference solution is computed for n = 9, the stopping criterion is a relative residual of 10 -8 and we choose N = 25. To compute the reference solution we use l ♭ = l/4. Still, the size of the multipole data structures for these parameters forces us to do compute the reference solution on a workstation with 48 GB RAM; these computations required roughly half of this RAM.

Figure 4 shows that the error curves for the multiple scattering problem behave precisely as to those for the single scatterer from Section 5. Especially, the error of the projection method behaves as predicted by theory. For α = 3/2, the error curve of the collocation method also fits well to the predicted rate, and for α = 0, 1/2, 1 the scheme does not reach the rate of the projection method.

As it is indicated by Theorem 7.5, one can discretize the coupling terms significantly sparser than the diagonal terms without spoiling the convergence of the scheme by choosing l ♭ < l. By exploiting this property, the diagonal terms of the operator matrices in (63) dominate the evaluation of the system matrix and hence also the complexity of the scheme. We demonstrate this feature by a numerical experiment that continues the above one. For α = 1, M = 4, N = 22 and l = (l, l) with l = (2 6 , 2 6 , 2 6 ) ⊤ we set l ♭ = (l ♭ , l ♭ ) with l ♭ = (2 m , 2 m , 2 m ) for m = 3, 4, 5. Figure 3 shows the relative error between the

|ñ|J 1

 1 (π|ñ|) and |z|=ρ v -ñ ds = πJ 0 (π|ñ|)

Figure 1 :

 1 Figure1: Relative L 2 -error of the collocation method and the projection method for waveguide medium scattering from inclusions with contrast q α ∈ W s for s < α + 1/2. Pluses, circles, dots, and crosses correspond to α = 0, 1/2, α = 1, and α = 3/2, respectively. The dotted lines indicate the convergence order 3/2 + α. (a) Collocation method: relative error versus n, l = (2 n , 2 n , 2 n ) ⊤ . (b) Projection method: relative error versus n, l = (2 n , 2 n , 2 n ) ⊤ .

Figure 2 :

 2 Figure 2: 2D-Sketch of a scatterer consisting of two parts D 1 and D 2 . The cylindersB (1) = o 1 + B ρ1/2 and B (2) = o 2 + B ρ2/2are contained in the boxes Λ(1) and Λ(2) , respectively. The computational domain Λ ρ1 is the reference domain for Λ(1) and T ± 1 transport functions between the two domains.

  )| in r to conclude that ∞ m=M sin (β3) (α m x 3 ) sin(α m y 3 )H (1) 0 (k m |x -ỹ|) ≤ C ∞ m=M |α m | β3 e -| Im (km)|δmin ≤ Ce -πδminM/hfor β 3 ∈ N. Taking partial derivatives with respect to x 2 or x 3 yields derivatives of H (1) 0 (k m |x -ỹ|). The resulting higher-order Bessel functions H (1) ν are then estimated again by (45) against an exponentially decaying sequence. The factor k m grows linearly in m and hence does not spoil this decay. Derivatives of x → a(x, y) = |x -ỹ| can be bounded by |∂a/∂x 1,2 | ≤ 1 and |∂ | β|1 a/∂x β | ≤ 2/δ min for any multiindex β ∈ N 2 0 of length | β| 1 = 2. Combining the latter estimates yields ∞ m=M sin(α m x 3 ) sin(α m y 3 ) ∂ | β|1 ∂x β H (1) 0 (k m |x -ỹ|) ≤ C exp(-δ min M ).

( 1 )

 1 n (k m r ij )| ≤ exp(δ min |k m * |) exp(-Im (k m )δ min )|H

( 1 )

 1 ν (|k m * |δ min )| for m ≥ m * and the monotonicity of H (1) ν , see (46) implies that |H (1) ν (k m r ij )| ≤ C exp(-Im (k m )δ min )|H (1) N (|k m * |δ min )| for m ≥ m * and |ν| ≤ N . Since 0 < k m * -1 < k, |H (1) ν (k m r ij )| ≤ |H (1) N (k m * -1 δ min )| for m < m * and |ν| ≤ N . Set k min = min(|k m * |, k m * -1 ).

4 N 4 N

 44 (c)-(d) indicates that for M = 4 the error is several orders of magnitude smaller than for M = 3. = 15 3.54e-05 3.54e-05 N = 25 6.94e-09 5.82e-10 N = 30 6.12e-07 6.12e-07 M = 3 M = 4 N = 10 4.36e-03 4.36e-03 N = 20 2.48e-08 2.48e-08 N = 30 4.63e-14 4.= 15 1.57e-04 1.54e-04 N = 25 2.85e-05 7.74e-10 N = 30 2.85e-05 9.02e-07 M = 3 M = 4 N = 10 9.78e-03 9.78e-03 N = 20 4.47e-07 4.47e-07 N = 30 1.13e-09 4.68e-13 (d) Table 2: Relative L 2 -error of multipole expansions for different truncation parameters M and N . (a) η = 1/2, k = 17.5 (b) η = 1/4, k = 17.5 (c) η = 1/2, k = 20.8 (d) η = 1/4, k = 20.8. Apart from η, k, M , and N , all other parameters are the same in (a)-(d). For M = 2, the relative errors are always ≥ 1.The exponential growth of the Hankel functions H

1 .N 3 2

 13 4kd < N ≤ 1.4kd + 12.3(kd) 1/3 ,15,N > 1.4kd + 12.3(kd) 1/3 , -kd/η 1.8(kd/η) 1/3 , N > kd/η.

  shows this error for the above-described experiment with wave number k = 17.5 for different values of η (1/2, 1/4, 1/8 and 1/12), corresponding to the extension from[-1/4, 1/4] × [-1/4, n + 1/4] × [0, 1/2] to into [n -1/4, n + 1/4] × [n -1/4, n + 1/4] × [0, 1/2],for n = 1, 2, 4 and 6. Setting d = q L 2 (Ω) always yields a good estimator N → 10 -min(u(N ),l(N )) for the numerical accuracy of the multipole expansion.

Figure 3 :

 3 Figure 3: Relative L 2 -error of the multipole expansion in dependence of N for four different values of η (all other parameters are the same in the three experiments; especially, M = 4 and l = (2 5 , 2 5 , 2 5 ) ⊤ throughout). Circles, crosses, triangles, and stars correspond to η = 1/2, 1/4, 1/8, and 1/12, respectively. The dotted curve N → 10 -u(N ) indicates an estimate for truncation error and the four dashed curves N → 10 -l(N ) for the above four different values of η indicate an estimate for the maximal numerical precision.

Table 1 :

 1 Computation times for the errors presented in Figure 1 in seconds. (a) Collocation method. (b) Projection method.

						α=0 α=0.5 α=1 α=1.5
					n=2	0.25	0.24	0.24	0.24
			0.93	0.93	n=3	0.93	0.93	0.93	0.93
	n=4 0.44	0.44	0.44	0.44	n=4	2.76	2.76	2.75	2.75
	n=5 1.62	1.63	1.64	1.64	n=5 13.36 12.93 13.17 13.00
	n=6 7.79	7.61	7.67	7.66	n=6 45.87 45.05 45.14 46.74
		(a)					(b)

Hence, C 2 (M ) grows only polynomially in M if ρ i + ρ j -δ min < 0, C(M ) ≤ M 2r m * -1 + exp((πm * /h -k)(ρ i + ρ j -δ min )) 1 -exp(π/h(ρ i + ρ j -δ min )) .

The kernel of K ij M,N,l1,l2 is obtained from χ i (x)G M,N (x + o i , y + o j )χ j (y) by application of Q l1 and Q l2 acting on the x-and on the y variables, respectively. Hence,

Λρ j (I -Q l2,y ) [G M,N (x + o i , y + o j )χ j (y + o j )] ϕ(y + o j ) dy

Thus,

where C is independent of M , N and l 1,2 . From Lemma 7.2 we know that (K ij -K ij M,N )ϕ W s+2 (Λρ i ) ≤ C(e -πδminM/h + M η N ) ϕ W s (Λρ j ) . Combining the last two estimates yields (61).

Next we derive error estimates for discretizations of the integral equation (43). As in Section 5, an important point is to project the product T -q • v l back into spaces of trigonometric polynomials. The two available options are the projection and interpolation operators P lj and Q lj into T l (Λ ρj ), the space of trigonometric polynomials T l defined on Λ ρj . For l = (l 1 , . . . , l J ) ∈ N 3×J we define a finite-dimensional product space T l = ⊕ J j=1 T lj (Λ ρj ). Let us set

For a discretization parameter l ∈ N 3×J we define a J × J matrix of operators K M,N,l with diagonal entries V ρi and off-diagonal entries K ij M,N,li,lj for 1 ≤ i, j ≤ J. The projection and collocation versions of the spectral/multipole discretization of (43) for the unknown v l ∈ T l are then

respectively. Here, l ♭ ∈ N 3×J can be chosen smaller than l to reduce the complexity of the scheme (see also Remark 7.4).

Let us now formulate the discrete scheme (63) in matrix-vector form. For the diagonal operators this has already been done in Section 5 and thus we concentrate now on the off-diagonal terms: Given the Fourier coefficients φ of ϕ ∈ T lj (Λ ρj ), we seek to compute the vector containing the Fourier coefficients of K ij M,N,li,lj ϕ. The interpolation projection of

for i, j = 1, . . . , J, m = 1, . . . , M , n = 1, . . . , 2N + 1. By b± m,n,i we denote the corresponding Fourier coefficients, considered as a column vector of size L i (the dimension of T li (Λ ρi )). The M (2N + 1) × L i matrix B ± i contains the vectors b± m,n,i as rows,

The column vector S ij , 1

The vector containing the Fourier coefficients of K ij M,N,li,lj ϕ can then be computed as

Thus, to apply the coupling operator K ij M,N,li,lj to ϕ we need to store three matrices B - i , S ij and B + j . Since S ij has M (2N + 1) complex entries and B ± i has M (2N + 1) × L i entries, which in total makes (2L i + 1)M (2N + 1) complex numbers. The matrix-vector products in (66) can be computed in O(M N (L i + L j )) operations.

Remark 7.4. The total number of operations to apply K M,N,l to ϕ ∈ T l of dimension L = 4Π J i=1 l i,j belongs to O(M N L log(L)). This shows the interest of discretizing the coupling terms of (63) sparser than the diagonal terms, i.e., to choose l ♭ < l. Setting L ♭ = 4Π J i=1 l ♭ i,j then the complexity of evaluating Let min(l) := min 1≤i≤3, 1≤j≤J l i,j for l ∈ R 3×J . Theorem 7.5. Assume that the Lippmann-Schwinger equation ( 6) is uniquely solvable in L 2 (Λ ρ ) for any right-hand side and that there exists a constant η ∈ (0, 1/2) such that r(x, y) < ηr ij for all x ∈ B (i) and y ∈ B (j) , 1 ≤ i = j ≤ J. Then the system (43) has a unique solution u ∈ W .

(a) Consider the Galerkin discretization in (63)(a). Assume that f ∈ W t , t ≥ 0, and that q ∈ P C is such that the multiplication u → T -q • u is continuous on W r for 0 ≤ r ≤ max(t -3/2, 0). Then for l ♭ ≤ l ∈ N 3×J with min(l), min(l ♭ ), M , and N large enough there is a unique solution u l ∈ T l of (63)(a) and

for 0 ≤ s ≤ t and arbitrary r > 0.

(b) Consider the collocation discretization in (63)(b). Assume that f ∈ W t for t > 3/2 and that q ∈ P C is such that the multiplication u → T -q • u is continuous on W r for 3/2 < r ≤ t. Then for l ♭ ≤ l ∈ N 3×J with min(l), min(l ♭ ), M , and N large enough there is a unique solution u l ∈ T l of (63)(b) and

for 0 ≤ s ≤ t and arbitrary r > 0.

Proof. If ( 6) is uniquely solvable in L 2 (Λ ρ ), then the same holds for the (equivalent) periodized vectorvalued system of integral equations (43). Under the stated smoothness assumption on q, the solution u to (43) belongs to

The remainder of the proof, for both cases, follows the proof of Theorem 5.4(b). One first exploits Proposition 7.2 to obtain an estimate for the approximation error in the off-diagonal terms of K M,N,l ♭ , and combines these estimates with the results of Theorem 5.3 on the approximation error in the diagonal terms. First,

for arbitrary r > 0. Hence, the operator norm of K -K M,N,l ♭ on W s can be made arbitrarily small by first choosing M and N large enough, and then choosing l ♭ = l ♭ (M, N ) large enough. Second, coupled multipole/spectral projection method for these parameters and the corresponding solution for m = 6. It is obvious that the error introduced by sparsifying the discretization of the coupling terms is way below the error of the solution of about 4e-06 even for m = 5. However, the memory requirements to store the multipole data structures are reduced from about 37 million for m = 6 to about 700000 complex numbers for m = 3. The scheme also speeds up as m is reduced.

Error 1.80e-12 1.77e-12 1.55e-12 Time 68.6 74.5 79.9

Table 3: Comparison of the coupled multipole/spectral projection method for l = (l, l) with discretization parameters l = (2 6 , 2 6 , 2 6 ) ⊤ and l ♭ = (l ♭ , l ♭ ), where l ♭ = (2 m , 2 m , 2 m ) ⊤ for m = 3, 4, 5, 6. The table shows the relative L 2 -error between the solutions for m = 3, 4, 5 and the solution for m = 6 (which corresponds to the same discretization for diagonal and off-diagonal terms), and the corresponding computation times in seconds.

A Proof of Lemma 4.2

We need to construct a basis {ϕ * n } n∈Z 3 l of T l such that ϕ * n (x (l) j ) = δ n,j . We construct ϕ * n as product of functions v * ñ depending on x and functions h * n3 depending on x 3 . For this construction we make used of the polynomials v ñ and h m from [START_REF] Chen | An adaptive perfectly matched layer technique for time-harmonic scattering problems[END_REF]. Denote V l = {v ñ : ñ ∈ Z 2 l } and define v * ñ by

It is shown in [START_REF] Saranen | Periodic integral and pseudodifferential equations with numerical approximation[END_REF]Sections 8.1 and 10.5] that the v * ñ form a basis of V l, and that v * ñ(x m) = δ ñ, m. To construct the special basis functions h * m of U l = {h n , 1 ≤ n ≤ l}, l ∈ N, we first consider trigonometric interpolation of 4h-periodic functions at grid points t m = h m/l, m = 0, 1, . . . , 4l -1. Note that the points t 1 , . . . , t l are the third components of the 3D grid points x (l) n . It is well known [START_REF] Kress | Linear Integral Equations[END_REF]Theorem 8.25

(68)

When given interpolation data p 1 , . . . , p l at t 1 , . . . , t l , we extend this data to all the 4l points t 0 , . . . , t 4l-1 by setting p 0 = 0 and then doing even and odd reflection at t l and t 2l , respectively,

Due to the two reflection symmetries, all β m and all γ m for even m in (68) vanish, that is, the interpolation polynomial satisfying p 2l (t m ) = p m , m = 0, 1, . . . , 4l -1, is j ) = δ n,j for n, j ∈ Z 3 l . The latter equation implies especially that the w * n are linear independent and since their number equals the dimension of T l they form a basis.