Tilo Arens

Drossos Gintides

AND Armin Lechleiter

DIRECT AND INVERSE MEDIUM SCATTERING IN A 3D HOMOGENEOUS PLANAR WAVEGUIDE

Keywords: scattering, inverse scattering, waveguide AMS subject classifications. 35P25, 45Q05, 65N21, 74J25, 76Q05, 86A05, 86A22

Time-harmonic acoustic waves in an ocean of finite height are modeled by the Helmholtz equation inside a layer with suitable boundary conditions. Scattering in this geometry features phenomena unknown in free space: resonances might occur at special frequencies and wave fields consist of partly evanescent modes. Inverse scattering in waveguides hence needs to cope with energy loss and limited aperture data due to the planar geometry. In this paper, we analyze direct wave scattering in a 3D planar waveguide and show that resonance frequencies do not exist for a certain class of bounded penetrable scatterers. More important, we propose the Factorization method for solving inverse scattering problems in the 3D waveguide. This fast inversion method requires near-field data for special incident fields and we rigorously show how to generate this data from standard point sources. Finally, we discuss our theoretical results in the light of numerical examples.

1. Introduction. Linear acoustics in a three-dimensional planar waveguide provides a simple model for the propagation of time-harmonic acoustic waves in the ocean [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF][START_REF] Buchanan | Marine Acoustics: Direct and Inverse Problems[END_REF][START_REF] Santiago | Modified Green's functions for shallow water acoustic wave propagation[END_REF][START_REF] Xu | The propagation solutions and far-field patterns for acoustic harmonic waves in a finite depth ocean[END_REF]. In this model, waves traveling inside a slab Ω = {x = (x 1 , x 2 , x 3 ) ⊤ ∈ R 3 : 0 < x 3 < h} of height h > 0 are governed by the Helmholtz equation ∆u + k 2 n 2 u = 0 in Ω, (

where k > 0 denotes the wave number and n the index of refraction. On the planar upper and lower boundaries of the waveguide, When choosing these boundary conditions for coherence with the literature, see, e.g., [START_REF] Bourgeois | Locating an obstacle in a 3d finite depth ocean using the convex scattering support[END_REF][START_REF] Xu | The propagation solutions and far-field patterns for acoustic harmonic waves in a finite depth ocean[END_REF], we note that the Dirichlet boundary condition models the sea surface, whereas the Neumann boundary condition models the ocean ground.

In our simple model, n is assumed to be the constant 1 throughout most of the waveguide, i.e. we are dealing with a homogeneous background medium. However, n may differ from this background value and be variable in space inside a bounded domain D ⊆ Ω. This model describes sound waves in the sea reasonably well if, e.g., the ocean depth and the water temperature are not too large.

Our primary interest lies in the inverse problem of determining the domain D from knowledge of the scattered field generated by given choices of incident fields. Further, we carefully analyze the direct scattering problems in the waveguide setting and show that uniqueness of solution for all wave numbers is available under special geometric conditions on the penetrable scattering object.

Inverse scattering in waveguides provides several special features that are not present in free-space inverse scattering problems. Part of the information is only present in evanescent fields that are undetectable away from the obstacle for all practical purposes. Due to the planar geometry the available scattering data is always of limited aperture. A large amount of measurement data is hence required for accurate reconstructions. Due to these theoretical challenges and the practical importance of the problem, inverse scattering problems in waveguides have received increasing attention in recent years [5, 7, 9-11, 13, 21, 25, 30]. The focus in most of these works is the scattering problem for a Dirichlet obstacle in a two-dimensional setting. The inhomogeneous medium scattering problem is considered in a two-dimensional setting in [START_REF] Buchanan | Marine Acoustics: Direct and Inverse Problems[END_REF] and references given therein, similarly in [START_REF] Dediu | Recovering inhomogeneities in a waveguide using eigensystem decomposition[END_REF] using a low frequency Born approximation.

Bearing in mind the above-mentioned difficulties, we will propose in this paper a Factorization method [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF] for carrying out the reconstruction. The factorization method is a fast reconstruction method that is able to provide estimates for the position and the shape of the scatterer, the quality of these estimates depending of course on the noise level of the data. The method's characterization principle is that a point z in the waveguide belongs to the support of the inhomogeneity in n if and only if the field generated by a point source at z is in the range of a certain linear integral operator defined by the data, termed the near field operator. The range of this operator can in turn be characterized by spectral properties. Our numerical examples show feasibility of this inversion method in three dimensions. A related method, the linear sampling method, has been investigated in [START_REF] Bourgeois | The linear sampling method in a waveguide: a modal formulation[END_REF][START_REF] Xu | Generalized dual space indicator method for underwater imaging[END_REF] for a two-dimensional waveguide and impenetrable obstacles. Though the derivation of the linear sampling method may appear simpler in many cases, we note that the Factorization method gives a mathematically exact characterization of the support of the inhomogeneity. Such a result is not available through linear sampling.

Solving an inverse problem requires a good understanding of the underlying direct problem. Hence we will start by describing an appropriate weak formulation of the direct problem and an equivalent formulation as a Lippmann-Schwinger type integral equation in Section 2. Analysis of the direct problem shows that in general uniqueness of solution cannot be expected for all positive wave numbers: there may be a sequence of wave numbers for which uniqueness of solution does not hold. However, under certain geometric conditions imposed on the index of refraction, uniqueness of solution to the direct problem does always hold, as we show in Section 3.

Analysis of the inverse problem under consideration starts in Section 4 where we present the factorization of the near field operator and present the characterization of the inhomogeneity through this factorization. As with characterization of obstacles through near field data in related problems [START_REF] Arens | The Factorization Method in inverse scattering from periodic structures[END_REF], the method requires data for incident fields that are not physical. For a practical implementation, the non-physical sources need to be generated by superpositions of physical fields. Although this has been proposed in various papers [START_REF] Arens | The Factorization Method in inverse scattering from periodic structures[END_REF][START_REF]An integral equation for Maxwell's equations in a layered medium with an application to the Factorization Method[END_REF], we present a first analysis of such an approximation using regularization techniques proving uniform approximation of the non-physical by the physical fields. The paper comes to a close with the presentation of numerical examples in Section 6.

2. Variational Solution of the Direct Scattering Problem. We will start this section with some notation used throughout the paper. Recall that the entire waveguide is denoted by Ω. The third coordinate axis is defined to be the one orthogonal to the waveguide, and we often combine the first two coordinates, writing

x = (x 1 , x 2 , x 3 ) ⊤ = x + (0, 0, x 3 ) ⊤ , x ∈ R 3 .
For discussions of the variational formulation of our scattering problem, we will mainly work in the bounded domain Ω R := {x ∈ Ω : |x| 2 < R 2 }, where the radius R is assumed to be large enough such that 1 + |x| 2 < R 2 for all x inside the scattering object D. This implies that D is contained in the interior of Ω R . The cylinder

C R := {x ∈ Ω : |x| 2 = R 2 }
denotes the part of the boundary of Ω R that is contained in Ω.

In the classical setting of an acoustic medium scattering problem the total field u is assumed to satisfy the Helmholtz equation (1.1). We will assume throughout the paper that the space dependent refractive index n 2 has positive real part, Re(n 2 (x)) > 0, x ∈ Ω, and non-negative imaginary part, Im(n 2 (x)) ≥ 0, x ∈ Ω. Some results presented are for the specific case Im(n 2 (x)) = 0. Moreover, we require the field u and its normal derivative to be continuous over interfaces where n 2 jumps. The closure of the scatterer D ⊂ Ω is defined as the support of the contrast q = n 2 -1. The set D is assumed to be a bounded Lipschitz domain with connected complement. In general, we do not require any regularity of n ∈ L ∞ (Ω) other than boundedness and measurability.

The incident field u i is assumed to satisfy the Helmholtz equation for constant refractive index, ∆u i + k 2 u i = 0 in Ω, subject to the boundary conditions (1.2). The direct scattering problem is then to find the total field u subject to (1.1), (1.2), such that u s := u -u i satisfies additionally a radiation condition. In the case of a waveguide problem, the radiation condition is usually obtained [START_REF] Gilbert | Dense sets and the projection theorem for acoustic waves in a homogeneous finite depth ocean[END_REF] by carrying out a separation of variables in x and x 3 , leading to a series expression

u s (x) = ∞ m=1 sin(α m x 3 )u m (x) for x ∈ Ω \ Ω R , (2.1)
with α m := (2m-1)π/(2h). The modes u m are required to satisfy the two-dimensional Helmholtz equation

∆u m + k 2 m u m = 0 |x| > R, with k m := k 2 -α 2 m , (2.2)
and the Sommerfeld radiation condition

lim r→∞ r 1/2 ∂u m ∂r -ik m u m = 0, where r := |x|, (2.3)
uniformly for all directions x/r. The modal wave numbers k m are only real for a finite number of values of m, say m ≤ M (k). These values correspond to modal frequencies of the waveguide. For m > M (k), the wave numbers k m are purely imaginary and correspond to exponentially decaying modes. The case k m = 0 for some m ∈ N corresponds to an exceptional frequency. In the sequel we always assume that the wave number k is such that k m = 0 for all m ∈ N.

(2.4)

Solutions to the waveguide scattering problem (1.1)-(2.1) are in the sequel understood in a variational sense and found using a variational formulation in W := {v ∈ H 1 (Ω R ) : v| Γ-= 0}. To incorporate the radiation condition in the variational formulation we rely on a so-called Dirichlet-to-Neumann operator Λ constructed and analyzed in [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF]. The trace space of W on C R , a closed subspace of H 1/2 (C R ) is denoted by V , and we set V ′ to be the dual of V for the inner product of L 2 (C R ). The boundary and radiation conditions (1.2) and (2.3) lead to a series representation of any solution u ∈ H 1 loc (Ω \ Ω R ) that satisfies the radiation conditions (2.3),

u(x) = ∞ m=1 l∈Z U m l H (1)
l (k m r) exp(ilϕ) sin(α m x 3 ), (2.5)

for cylindrical coordinates x = (r cos ϕ, r sin ϕ, x 3 ) ⊤ ∈ Ω R , where H

(1) l denotes the Hankel function of the first kind and of order l and U m l are the Fourier coefficients of u| CR , defined as U m l = π/(hR) CR u exp(-ilϕ) sin(α m x 3 ) ds for m ∈ N, l ∈ Z. From the representation (2.5), the normal derivative of u on C R can be computed,

∂u ∂r (x) = ∞ m=1 l∈Z k m U m l H (1) l ′ (k m R) exp(ilϕ) sin(α m x 3 ), x ∈ C R .
Therefore we set, for u ∈ V with Fourier coefficients U m l ,

Λu(x) := ∞ m=1 αm =k l∈Z k m U m l H (1) l ′ (k m R) H (1)
l (k m R) exp(ilϕ) sin(α m x 3 ). ( 2.6)
The following lemma is proved in Section 2 of [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF]. Lemma 2.1. Suppose k > 0 such that (2.4) is satisfied. The Dirichlet-to-Neumann operator Λ defined in (2.6) is a bounded operator from V to V ′ for all k > 0. Moreover, there is C = C(k) > 0 such that for all u ∈ V there holds

-Re CR u Λu ds ≥ -C u 2 L 2 (CR) , u ∈ V.
Finally, there is a neighborhood U ⊂ C of k such that Λ is an analytic operator-valued function of the wave number in U .

Using the Dirichlet-to-Neumann operator Λ, the variational formulation for the scattered field u s = u -u i is as follows: Find u s ∈ W such that, for all v ∈ W ,

B(u s , v) := ΩR ∇u s • ∇v -k 2 n 2 u s v dx - CR v Λu s ds = k 2 ΩR qu i v dx. (2.7) Problem 2.2. Given an incident field u i satisfying ∆u i +k 2 u i = 0 in Ω and (1.2), find u s ∈ W such that (2.7) holds for all v ∈ W .
Due to Lemma 2.1, the form B satisfies a Gårding inequality, implying by the Fredholm alternative that the variational problem is solvable whenever there is at most one solution.

Theorem 2.3. Suppose that Im(n 2 ) > 0 in some non-empty open subset of D. Then the scattering problem 2.2 is uniquely solvable for any incident field u i .

The theorem follows by setting v = u s in (2.7) and taking the imaginary part. In the proof of Lemmas 4.3 in [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF] the estimate Im CR u s Λu s ds ≥ 0 is shown. Hence u s = 0 allmost everywhere in D and by unique continuation in Ω R follows.

The case of real index of refraction n is much more difficult to treat and nonuniqueness phenomena might occur. Due to the analyticity properties of the Dirichletto-Neumann map mentioned in Lemma 2.1, analytic Fredholm theory [12, Theorem I.5.1] may be used to treat this case. It turns out that the variational problem 2.2 is uniquely solvable for all but possibly a discrete set of wave numbers that accumulate at most at infinity. We announce this result in the following theorem, that can be shown by arguments analogous to those provided in [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF].

Theorem 2.4. The medium scattering problem 2.2 is uniquely solvable for any incident field u i except possibly for a sequence of real wave numbers (k (j) ) such that k (j) → ∞ (j → ∞).

We continue by developing an integral formulation for the direct problem based on a volume potential defined in D, complementing the variational approach discussed above. In our later discussion of the Factorization method for the inverse problem, we require this integral equation formulation since it provides an analytic way to describe solution operators. The volume potential we consider here relies on the Green's function for the waveguide Ω with respect to the Dirichlet boundary condition on the bottom Γ -and the Neumann condition on Γ + . We denote this Green's function by G(•, •), defined on Ω × Ω \ {(x, x) : x ∈ Ω}. The following representation of G is given in [START_REF] Ahluwalia | Exact and asymptotic representations of the sound field in a stratified ocean[END_REF],

G(x, y) = i 2h ∞ m=1 sin (α m x 3 ) sin (α m y 3 ) H (1) 0 (k m |x -ỹ|) , x, y ∈ Ω, x = ỹ. (2.8)
An equivalent representation of the Green's function in D can be derived via the method of images,

G(x, y) = 1 4π +∞ m=-∞ (-1) m e ik|x-ym| |x -y m | - e ik|x-y ′ m | |x -y ′ m | (2.9)
where the source image points are given by y m = (y 1 , y 2 , y 3 ) ⊤ + (0, 0, 2mh) and y ′ m = (y 1 , y 2 , -y 3 ) + (0, 0, 2mh) , m ∈ Z.

From (2.9) it is obvious that G(•, •) can be written as the superposition of the fundamental solution Φ(x -y) = e ik|x-y| /(4π |x -y|) of the Helmholtz equation in free space and an analytic function G(•, •) that is,

G (x, y) = 1 4π e ik|x-y| |x -y| + G(x, y), x = y ∈ Ω. (2.10)
Consequently, the mapping properties of the volume potential

Vf = D G(•, y)f (y) dy
are the same as for the volume potential with kernel Φ. From the first two sections of [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Chapter 6] we infer that V is a bounded operator from L 2 (D) into H 2 loc (Ω). Moreover, the fact that G(•, •) is a Green's function for the boundary conditions (1.2) implies that u = Vf belongs to H 2 loc (Ω) for f ∈ L 2 (D), it solves the inhomogeneous Helmholtz equation ∆u + k 2 u = -f in Ω subject to the boundary conditions (1.2), and it satisfies the radiation conditions (2.1). In this last equation, we understand f to be extended by zero outside of D, a convention which we will also use in some of the later equations.

Due to these properties of the volume potential V, the Lippmann-Schwinger integral equation provides an alternative tool to solve the medium scattering problem 2.2. Even more generally, let us seek a solution u ∈ H 2 loc of the following source problem in Ω,

∆u + k 2 n 2 u = -k 2 qf, f ∈ L 2 (D), (2.11)
together with the waveguide boundary conditions (1.2) and the radiation condition (2.1). Here we used again the medium's contrast q = n 2 -1 that vanishes outside D. For f = u i D a solution to this problem provides a scattered field u s to the waveguide scattering problem 2.2.

Using Green's third identity [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Theorem 6.10], one finds that a solution u to the source problem (2.11) satisfies the Lippmann-Schwinger integral equation

u| D = k 2 D G(•, y)q(y)u(y) dy + k 2 D G(•, y)q(y)f (y) dy in D.
(2.12)

If u ∈ L 2 (D) solves (2.
12) then the right hand side of the equation defines an extension of u to the entire waveguide Ω. The mapping properties of the volume potential V imply that this extension is twice weakly differentiable. Since G is the waveguide Green's function the extension moreover solves the source problem (2.11). Thus, the variational approach to the scattering problem (or to the source problem) is equivalent to the integral equation approach. Furthermore, any solution to the scattering problem (2.7) belongs to H 2 loc (Ω). Writing the Lippmann-Schwinger equation (2.12) in the form

u -k 2 V(qu) = k 2 V(qf ) in L 2 (D), (2.13)
shows that we can apply the Riesz theory [START_REF] Kress | Linear integral equations[END_REF] to prove solvability of (2.13). In particular, existence and stability of a solution follow from uniqueness. However, for f = 0, the solution also solves the scattering problem 2.2 for u i = 0. Hence, under any of the conditions given in Theorems 2.3 and 2.4, (2.13) is uniquely solvable in H 2 (D) for f ∈ L 2 (D). Additionally, Corollary 3.3 below provides geometric conditions on q for uniqueness of (2.13).

3. Conditions on n for Uniqueness of the Direct Scattering Problem. The existence theory of the previous section does not exclude non-uniqueness or resonance phenomena that might possibly occur at a discrete set of wave numbers for purely real index of refraction n 2 . In general this result is optimal due to nonuniqueness examples for waveguide scattering problems [START_REF] Morgenrother | On the instability of resonances in parallel-plane waveguides[END_REF]. However, special classes of penetrable scattering objects do not feature resonances at any frequency. This is due to a Rellich identity and the asymptotic behavior of solutions of the homogeneous problem. The integral identities we use go back to Rellich [START_REF] Rellich | Darstellung der Eigenwerte von ∆u + λu = 0 durch ein Randintegral[END_REF]. They have been applied, e.g., to obtain frequency-dependent bounds for scattering problems in [START_REF] Cummings | Sharp regularity coefficient estimates for complex-valued acoutic and elastic Helmholtz equations[END_REF], and for existence and uniqueness theory for periodic surface scattering [START_REF] Bonnet-Bendhia | Guided waves by electromagnetic gratings and non-uniqueness examples for the diffraction problem[END_REF] and impenetrable obstacle scattering in waveguides [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF].

The next lemma proves that any solution and its first derivatives of the homogeneous scattering problem are exponentially decaying as the distance from the x 3 axes increases.

Lemma 3.1. Let Im(n 2 ) = 0. Assume that u ∈ W is such that B(u, v) = 0 for all v ∈ W .
We extend u uniquely by (2.5) to a radiating solution of the Helmholtz equation in all of Ω. Then there exist constants C, c > 0 such that

|u(x)| + |∇u(x)| ≤ Ce -c|x| , |x| > R.
The proof is similar to the proof of [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF]Lemma 4.3]. Next we present several Rellich type identities that are the key step for proving uniqueness of solution to the waveguide scattering problem under additional conditions. These identities involve integrals containing derivatives of the contrast q = n 2 -1. Since we suppose q ∈ L ∞ (Ω R ) these derivatives have to be understood in a distributional sense, i.e.,

ΩR ∂q ∂x j v dx := - ΩR q ∂v ∂x j dx for v ∈ C ∞ (Ω R ), j = 1, 2.
There arise no boundary terms when applying the divergence theorem since the support of q does not touch the cylinder C R and the components ν j of the normal vector satisfy ν j = 0, j = 1, 2, on Γ ± . The right-hand side in the latter definition can be continuously extended to v ∈ W 1,1 (Ω R ), that is, the derivatives of q are well-defined as elements in the dual of W 1,1 (Ω R ). We moreover introduce the following short hand notation: For x ∈ Ω and for a differentiable function u, we write

x = (x 1 , x 2 , 0) ⊤ and ∇u = (∂u/∂x 1 , ∂u/∂x 2 , 0) ⊤ . Lemma 3.2 (Rellich identity). Assume that n 2 ∈ L ∞ (Ω R ) with Im(n 2 ) = 0. Let u ∈ W ∩ H 2 (Ω R ) be a solution of ∆u + k 2 n 2 u = 0 and ∂u/∂ν = 0 on Γ + R .
Then, for q = n 2 -1, and j = 1, 2, the following identity holds,

k 2 ΩR x j ∂q ∂x j |u| 2 dx + 2 ΩR ∂u ∂x j 2 dx = CR x j ν j k 2 |u| 2 -|∇u| 2 + 2x j Re ∂u ∂x j ∂u ∂ν + u ∂u ∂ν ds. (3.1)
Proof. An application of Green's first identity implies

ΩR x 1 ∂u ∂x 1 ∆u dx = - ΩR ∇ x 1 ∂u ∂x 1 • ∇u dx + CR x 1 ∂u ∂x 1 ∂u ∂ν ds = - ΩR ∂u ∂x 1 2 + x 1 ∇ ∂u ∂x 1 • ∇u dx + CR x 1 ∂u ∂x 1 ∂u ∂ν ds.
Note that the integrals over the boundaries Γ ± appearing after Green's identity drop out due to the waveguide boundary conditions: ∂u/∂ν = 0 on Γ + and ∂u/∂x 1 = 0 on Γ -. By simple differentiation we get

∂ ∂x 1 |∇u| 2 = 2 Re ∇u • ∇ ∂u ∂x 1 = 2 Re ∇ ∂u ∂x 1 • ∇u . Therefore 2 Re ΩR x 1 ∂u ∂x 1 ∆u dx = - ΩR 2 ∂u ∂x 1 2 + x 1 ∂|∇u| 2 ∂x 1 dx + 2 Re CR x 1 ∂u ∂x 1 ∂u ∂ν ds = ΩR -2 ∂u ∂x 1 2 + |∇u| 2 dx + CR -x 1 |∇u| 2 ν 1 + 2x 1 Re ∂u ∂x 1 ∂u ∂ν ds.
On the other hand, from the equation ∆u + k 2 n 2 u = 0, which holds in

L 2 (Ω R ) since u ∈ H 2 (Ω R ), we derive the identity 2 Re ΩR x 1 ∂u ∂x 1 ∆u dx = -2k 2 ΩR x 1 Re ∂u ∂x 1 (q + 1) u dx = -k 2 ΩR x 1 (q + 1) ∂ ∂x 1 |u| 2 dx = -k 2 ΩR x 1 q ∂|u| 2 ∂x 1 dx -k 2 ΩR x 1 ∂|u| 2 ∂x 1 dx = k 2 ΩR ∂ ∂x 1 (x 1 q) |u| 2 dx + k 2 ΩR |u| 2 dx -k 2 CR x 1 |u| 2 ν 1 ds = k 2 ΩR (q + 1)|u| 2 + x 1 ∂q ∂x 1 |u| 2 dx -k 2 CR x 1 |u| 2 ν 1 ds.
As we explained above, the second term in the first integral has to be interpreted in a distributional sense. Again, the boundary terms on Γ ± arising from this partial integration cancel, since ν 1 = 0 on Γ ± . From the last two equations, and using again

n 2 = q + 1, we obtain ΩR -2 ∂u ∂x 1 2 + |∇u| 2 dx + CR -x 1 |∇u| 2 ν 1 + 2x 1 Re ∂u ∂x 1 ∂u ∂ν ds = k 2 ΩR n 2 |u| 2 + x 1 ∂q ∂x 1 |u| 2 dx -k 2 CR x 1 |u| 2 ν 1 ds.
Rearranging terms and using equation (2.7) in the special case v = u, we conclude

k 2 ΩR x 1 ∂q ∂x 1 |u| 2 dx + 2 ΩR ∂u ∂x 1 2 dx = CR x 1 ν 1 k 2 |u| 2 -|∇u| 2 + 2x 1 Re ∂u ∂x 1 ∂u ∂ν + u ∂u ∂ν ds.
The very same arguments show that the last equality remains correct if we replace x 1 and ν 1 by x 2 and ν 2 , respectively. The approach used in this proof cannot be used to establish (3.1) for j = 3, since the boundary terms on Γ + do not all vanish. Now we formulate the announced uniqueness statement.

Corollary 3.3. Assume that the contrast q satisfies either

Ω ∂q ∂x j x j v dx ≥ 0 for all v ∈ C ∞ (Ω), v ≥ 0 and j = 1 or 2, (3.2)
or else

Ω x • ∇q v dx ≥ 0 for all v ∈ C ∞ (Ω) with v ≥ 0. (3.3)
Then the variational problem 2.2 is uniquely solvable. Condition (3.3) is probably more intuitive when it is expressed using cylindrical coordinates (a 1 r cos(φ), a 2 r sin(φ), x 3 ) ⊤ with a 1 , a 2 > 0 . Since ∂q/∂r = r -1 ∇q •

x we can reformulate (3.3) as Ω ∂q/∂r r 2 v dr dφ ≥ 0 for all v ∈ W 1,1 (Ω) with v ≥ 0. Roughly speaking, this condition implies that the refractive index increases in radial direction.

Proof. Assume on the contrary that there is, under either of the above assumptions, u ∈ W, u = 0 such that the variational problem 2.2 is not uniquely solvable, that is, B(u, v) = 0 for all v ∈ W . We extend u to the unique radiating solution of the Helmholtz equation outside of Ω R . This extension is independent of the radius R of the domain Ω R and exponentially decaying at infinity due to Lemma 3.1. By well known elliptic regularity results, u ∈ H 2 (Ω). Therefore the term

CR x j ν j k 2 |u| 2 -|∇u| 2 + 2x j Re ∂u ∂x j ∂u ∂ν + u ∂u ∂ν ds
tends to zero as R → ∞. Using now (3.1) we find that

k 2 Ω x j ∂q ∂x j |u| 2 dx + 2 Ω ∂u ∂x j 2 dx = 0.
However, if assumption (3.2) holds, then ∂u/∂x j = 0 in W . Therefore u does not depend on the x j variable. Since u is exponentially decaying this implies that u vanishes. If assumption (3.3) holds, the same conclusion can be derived in an analogous manner, since by adding (3.1) for j = 1 and 2 one finds

k 2 ΩR (x • ∇q) |u| 2 dx + 2 ΩR ∇u 2 dx = CR R k 2 |u| 2 -|∇u| 2 + 2 Re (x • ∇u) ∂u ∂ν + 2u ∂u ∂ν ds.
4. The Inverse Problem and Factorizations of the Near Field Operators. The inverse problem that we consider is to find the support D of q from the scattered field u s on the boundary ∂B of some sufficiently smooth domain B containing the scatterer D for the incident field u i = G(•, y), y ∈ ∂B. We will use the notation u s (x, y), x, y ∈ ∂B to denote this field. In fact, we only require knowledge of this data for x, y in some open subset M of ∂B.

We start by reviewing a uniqueness result for this inverse problem shown in [START_REF] Ikehata | Inverse boundary value problem for ocean acoustics using point sources[END_REF]. In this paper, the authors consider the inverse problem of determining n from point source excitations and measurements on the surface of a cylinder containing the object. More precisely, the required data is the trace of the radiating Green's function for the Helmholtz equation ∆u + k 2 n 2 u on the surface of a vertical cylinder C R surrounding the inhomogeneity. The problem was solved using a similar approach as in Nachman's reconstruction method from boundary measurements [START_REF] Nachman | Reconstructions from boundary measurements[END_REF].

Theorem 4.1 (See [START_REF] Ikehata | Inverse boundary value problem for ocean acoustics using point sources[END_REF]). Assume that Problem 2.2 is uniquely solvable for any incident field u i and suppose further that the interior problem

∆u + k 2 n 2 u = 0 in Ω R , u = 0 on ∂Ω R \ Γ + , ∂u ∂x 3 = 0 on ∂Ω R ∩ Γ + ,
possesses only the trivial solution. If n 1 and n 2 are refractive indices such that the scattered field u s (x, y) is the same for x, y ∈ M , then n 1 = n 2 . Proof. Our assumptions mean that all assumptions of the uniqueness proof in [START_REF] Ikehata | Inverse boundary value problem for ocean acoustics using point sources[END_REF] are satisfied. Moreover, the Green's function of the problem needs to be known for all x, y ∈ C R . However, the Green's function in our notation is u s (x, y) -G(x, y). So, if we know u s (x, y) for x, y ∈ M , then by analytic continuation we know the Green's function of the problem for x, y ∈ C R and the result follows directly.

Even though this theorem indicates that we can uniquely reconstruct D from the data u s , we are not able to do so directly using the Factorization Method. The reason will become clear below. Hence, we modify the available data and, abusing notation, from now on denote by u s (x, y) the scattered field at y ∈ ∂B for the incident conjugate point source u i = G(•, y), y ∈ ∂B. As pointed out in the introduction, these scattered fields cannot be produced from a physical scattering process. However, we will address this issue in Section 5.

We now define the so-called near field operator N : L 2 (M ) → L 2 (M ) by setting

N ϕ = M u s (•, y) ϕ(y) ds(y), ϕ ∈ L 2 (M ) . (4.1)
In the Factorization method, a factorization of N is derived that serves in defining a function that characterizes the scatterer. To clarify the role of the individual operators in this factorization, we will derive such a factorization first in the simpler model of the Born approximation.

From Section 2, we know that the scattered field u s (•, y) corresponding to the incident field G(•, y) satisfies the Lippmann-Schwinger equation

u s (•, y) -k 2 D q(z) u s (z, y)G(•, z)dz = k 2 D q(z)G(z, y)G(•, z)dz (4.2)
for y belonging to the measurement surface M . This formula holds throughout Ω and is written in Born approximation for small k > 0 approximately as

u s Born (•, y) = k 2 D q(z)G(z, y)G(•, z)dz, y ∈ M.
The near field operator N Born : L 2 (M ) → L 2 (M ) in Born approximation is

N Born ϕ = M u s Born (•, z) ϕ(z) ds(z),
and a straight-forward calculation gives the factorization N Born = H * T Born H. Here, the Herglotz operator H :

L 2 (M ) → L 2 (D) is defined by (Hϕ)(x) = |q(x)| M G(x, y) ϕ(y) ds(y), x ∈ D, (4.3)
and represents propagation from the sources to the scatterer. Its adjoint

H * : L 2 (D) → L 2 (M ) is given by (H * g)(x) = D |q(y)|G(x, y) g(y) dy, x ∈ M.
and represents propagation from the scatterer to the receivers. Note that H * characterizes the scatterer, see Theorem 4.3 below. Finally, the operator T Born is defined as T Born f = k 2 sign(q) f for f ∈ L 2 (D) where the signum of q is defined as sign(q) = q/|q|. It represents the reflectivity of the scatterer.

For the functional analytic framework of the Factorization method it is essential that the two outer operators be adjoint of each other and that H * characterizes the scatterer. Hence we are given no choice for H which represents unphysical "incident" fields.

Let us now analyze the situation in the full model taking multiple scattering into account. Consider the solution operator G :

L 2 (D) → L 2 (M ) with f → Gf = v| M where v solves the source problem ∆v + k 2 (1 + q)v = -k 2 q |q| f in Ω, (4.4)
subject to the radiation conditions (2.3). From the linear superposition principle, it follows that N = GH where H is the Herglotz-like operator defined above. We rewrite (4.4) in the form

∆v + k 2 v = -k 2 q |q| f + |q| v =: -|q| T f
and conclude that v can be written as v = D |q(y)| G(•, y) T f (y) dy. Hence v| M = H * T f and thus Gf = H * T f . This yields

N = H * T H. (4.5) Note that T : L 2 (D) → L 2 (D) is given explicitly by T f = k 2 sign(q) f + |q|v| D
where v solves (4.4). Hence, we see exactly the modification over the Born approximation required to the central operator to account for multiple scattering, while the outer operators in the factorization do not change. We now use the factorization (4.5) to characterize the support of the contrast q. Although the factorization (4.5) holds for general contrasts, the Factorization method itself relies up to now on a positivity assumption, which will be guaranteed in our context by the following condition on the contrast q = n 2 -1. Recall that in the first section on the direct problem we already assumed that Re(n 2 ) = 1 + Re(q) > 0 and Im(n 2 ) = Im(q) ≥ 0. We need to strengthen these assumptions from now on for the rest of this work and suppose that (A1) Either there is c > 0 such that Re q ≥ c|q|, or (A2) there is c > 0 such that -Re q ≥ c|q|. Using the latter assumptions we will prove in this section that the real part of the middle operator T of the Factorization N = H * T H is the sum of a coercive and a compact operator and that the imaginary part of T is non-negative. We will use this property in the sequel to characterize the range of H * and provide thereby a characterization of points in the support D of the contrast q.

We call an operator A on some Hilbert space

X coercive if Re Af, f X ≥ c f 2 X
for all f ∈ X and some c > 0. Moreover, real and imaginary parts of A are defined as Re A = (A + A * )/2 and Im A := (A -A * ) /(2i).

Proposition 4.2. Under Assumption (A1), the operator Re T is a compact perturbation of a coercive operator and under Assumption (A2), -Re T is a compact perturbation of a coercive operator. In both cases, Im T is non-negative and T is an isomorphism.

Proof. Let v denote the unique radiating solution to (4.4) for f ∈ L 2 (D). We split T = T 0 + T 1 into two operators T 0 f = k 2 sign(q)f and T 1 f = k 2 sign(q) |q| v| D , and note that coercivity of T 0 and compactness of T 1 follow as in [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF]. For the imaginary part of q we compute

D f T f dx = k 2 D sign(q)w f dx = k 2 D sign(q)w w -|q|v dx = k 2 D sign(q) |w| 2 dx -k 2 D w q |q| v dx
where we set

w = f + |q| v| D . Note that ∆v + k 2 v = -k 2 qw/ |q|. Consequently, Green's first identity in Ω R \D implies D f T f dx = k 2 D sign(q) |w| 2 dx + |x|<R k 2 |v| 2 -|∇v| 2 dx + CR v ∂v ∂ν ds.
Taking the imaginary part of the latter equation results in

Im D f T f dx = k 2 D Im q |q| |w| 2 dx + Im CR v ∂v ∂ν ds.
Since Im q ≥ 0 the first integral is non-negative. For the second one, Proof. Since H * T = G and T is according to the last proposition an isomorphism, it is sufficient to show that φ z belongs to the range Rg(G) if and only if z ∈ D.

If z ∈ D assume that there exists

ψ ∈ L 2 (D) with Gψ = φ z | M = G(•, z)| M .
Then by the unique continuation principle D |q| G(x, y) T ψ(y) dy = G(x, z) for all x = z exterior to D. This is a contradiction because the integral is an analytic function in Ω \ D but G(•, z) has a singularity as at z.

In the remainder, we assume that (B2) holds and that z ∈ D and show φ z does belong to the range of H * , or, equivalently, to the range of G. For assumption (B1) the proof given in [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF] can easily be adapted to our setting.

Assuming that (B2) holds, choose z ∈ D. We further define D ε as in (B2) and choose ε small enough that q

-1 ∈ L 1 (D ε ) and that z ∈ D \ D ε . Let χ ∈ C ∞ (Ω) with values in [0, 1] be such that χ = 1 in (Ω \ D) ∪ D ǫ/2 and χ = 0 in D \ D ǫ . Then we define v(x) = χ(x) G(x, z) for x ∈ D. It holds v = G(x, z) outside D \ D ǫ/2 and in particular on M . Defining f by f = - √ q k 2 q ∆v + k 2 (1 + q)v
almost everywhere in D, we obtain from Assumption (B2) that f is square integrable in D. It follows that ∆v

+ k 2 (1 + q)v = -k 2 qf / √ q.
By the unique solvability of the source problem (4.4) for source terms in L 2 (D), which we assumed in the beginning of this section, we conclude that v = G(•, z) outside D and in particular that φ z = Gf ∈ Rg(G).

We have now shown all prerequisites necessary to apply [19, Theorem 2.1]. Theorem 4.4. Let X ⊂ Y ⊂ X * be a Gelfand triple with Hilbert space Y and reflexive Banach space X such that the embedding is dense. Furthermore, let Z be a second Hilbert space and F : Z → Z, H : Z → X and T : X → X * be linear and bounded operators with F = H * T H. We make the following assumptions:

(a) H is compact and injective. (b) Re T has the form Re T = T 0 + T 1 with some coercive operator T 0 and some compact operator T 1 : X → X * . (c) Im T is non-negative on X, i.e., Im T φ, φ ≥ 0 for all φ ∈ X. From this result on range identities in factorizations applied to the factorization N = H * T H we derive our final result of this section, which is an explicit reconstruction of D in terms of N . . Denote by (λ j , ψ j ) j∈N an eigensystem of the selfadjoint and positive operator N ♯ . Then

z ∈ D ⇔ ∞ j=1 φ z , ψ j L 2 (M) 2 λ j < ∞. (4.6)
Proof. Combination of Theorem 4.4 with Proposition 4.2 and the factorization of N yields that Rg N 1/2 ♯ = Rg(H * ). Theorem 4.3 then shows the first claim. The second characterization of D given in (4.6) follows from Picard's criterion [START_REF] Kirsch | An Introduction to the Mathematical Theory of Inverse Problems[END_REF].

5

. Approximating Unphysical Sources. We will in this section address the problem of obtaining the data for the construction of the operator N using physically relevant incident fields. Let us denote by N phys the linear integral operator defined as N in (4.1), but with u s (•, y) the scattered field due to an incident point source G(•, y) with y ∈ M . We present a regularization method to approximate the auxiliary operator N from N phys using a family of operators {P δ } δ>0 such that N phys P δ → N in the operator norm of L 2 (M ). This idea has been presented in [START_REF]An integral equation for Maxwell's equations in a layered medium with an application to the Factorization Method[END_REF], however, we improve the results of this reference substantially, proving uniform convergence of the approximation instead of pointwise convergence. The main ingredient in the construction is the a-priori knowledge of a bounded connected test domain B 0 which includes D but excludes the measurement surface M ⊂ Ω. Let us denote the boundary ∂B 0 by Σ and furthermore assume that one of the two following conditions holds:

(C1) M is an open non-empty subset of a vertical plane, or (C2) M is an open non-empty subset of a star-shaped cylinder enclosing D, and the cylinder's cross section is given by an analytic curve.

See Figure 5 for a sketch of this geometry. We assume furthermore that k 2 is not an eigenvalue for the Laplace operator in Construction of P δ requires a few definitions. First, for a smooth surface γ contained in Ω, we introduce the single layer potential

(SL γ ϕ)(x) := γ G(x, y) ϕ(y) ds(y), x ∈ Ω \ γ.
By SL γ ϕ we denote the corresponding potential with conjugate kernel G(x, y). We will use the choices γ = M and γ = Σ in the following arguments. Furthermore, we will make use of the evaluations on Σ of the potentials defined on M ,

V ϕ = SL M ϕ| Σ and Ṽ ϕ = SL M ϕ| Σ .
Lastly, we will use the corresponding single layer operator on Σ,

(S Σ ϕ)(x) := Σ G(x, y) ϕ(y) ds(y), x ∈ Σ.
The representation (2.9) implies that G(x, y) differs from the free-space fundamental solution Φ(x, y) by a smooth function, thus mapping properties of S Σ follow from mapping properties of the single layer operator for the free-space fundamental solution. Especially, S Σ :

H -1/2 (Σ) → H 1/2 (Σ) is bounded. The same argument shows that both V and Ṽ are bounded from L 2 (M ) into H 1/2 (Σ). Lemma 5.1. Under assumption (C1) or (C2), V , Ṽ : L 2 (M ) → H 1/2 (Σ) are injective with dense range.
Proof. Injectivity of V and Ṽ follows at once from the assumption that k 2 is not a Dirichlet eigenvalue of -∆ in B 0 . Furthermore, if the range of V is not dense in

H 1/2 (Σ), we find ψ = 0 such that V ϕ, ψ L 2 (Σ) = 0 for all ϕ ∈ L 2 (M ). By duality, ϕ, SL Σ ψ L 2 (M) = 0, that is, SL Σ ψ = 0 on M .
In the case that (C1) holds, that is, M is part of a vertical hyperplane Γ, we conclude by analytic continuation that SL Σ ψ vanishes on Γ ∩ Ω. The complex conjugate u := SL Σ ψ = SL Σ ψ moreover satisfies the radiation conditions (2.3).

The expansion (2.1) of u into u(x) = u m (x) sin(α m x 3 ) implies that the modes

u m (r, φ) = l∈Z U m l H (1)
l (k m r) exp(ilφ), which satisfy the Sommerfeld radiation condition, vanish on an entire straight line Γ = {x : x ∈ Γ}. Note that the scattering object lies on one side of Γ. Without loss of generality, we can assume that Γ is the x 1 -axis and that the scattering object is contained in {x 2 < 0}. Since u m satisfies the Sommerfeld radiation condition in the upper half-space, the odd reflection u r m defined by u r m

(x) = u m (x) in {x 2 > 0}, u r m (x) = -u m (x 1 , -x 2 ) in {x 2 <
0} is an entire radiating solution to the Helmholtz equation. Therefore u r m vanishes and due to the unique continuation property both u m and u = SL Σ ψ vanish in the exterior of B 0 . Together with our assumption on k 2 not being an eigenvalue of -∆, we conclude by the jump relations for the single layer potential that ψ vanishes -a contradiction.

In the case that (C2) holds, that is, M is part of a star-shaped cylinder Γ with analytic cross-section, we conclude by analytic continuation that SL Σ ψ vanishes on Γ ∩ Ω. The function u := SL Σ ψ = SL Σ ψ thus satisfies a homogeneous Dirichlet scattering problem in the exterior domain to Γ. The Rellich identity for Dirichlet scatterers shown in [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF]Lemma 4.1] implies that such a solution must vanish (see [START_REF] Arens | Variational formulations for scattering in a 3-dimensional acoustic waveguide[END_REF]Theorem 4.4] for a proof) and we conclude as above that ψ vanishes.

The very same technique shows that the range of Ṽ is dense in L 2 (Σ), too, and we omit the proof.

To approximate the unphysical incident field Hϕ by a radiating field on the obstacle D we approximate Ṽ ϕ by V ψ on Σ. A continuous dependence results will then imply approximation of the irradiating field by the radiating one in all of B 0 which contains D. However, since the integral equation of the first kind to find ψ in L 2 (M ) such that V ψ = Ṽ ϕ for given ϕ ∈ L 2 (M ) is ill-posed, we apply Tikhonov regularization to this equation and define

P δ : L 2 (M ) → L 2 (M ), P δ ϕ = (δ + V * V ) -1 V * Ṽ ϕ, δ > 0.
Proposition 5.2. N phys P δ converges to N pointwise in L 2 (M ).

Proof. From standard regularization theory [START_REF] Kirsch | An Introduction to the Mathematical Theory of Inverse Problems[END_REF] for ill-posed problems we obtain that V P δ ϕ → Ṽ ϕ pointwise for all ϕ ∈ L 2 (M ) as δ → 0. As the boundary values V P δ ϕ of the single layer potential SL M P δ ϕ converge in H 1/2 (Σ) to the boundary values Ṽ ϕ of SL M ϕ, our assumption that k 2 is no Dirichlet eigenvalue of -∆ in B 0 implies that SL M P δ ϕ → SL M ϕ pointwise in H 1 (B 0 ). Consequently, the incident fields used to define N phys converge pointwise to the incident fields used to define N . Linearity of problem (2.11) yields that N phys P δ ϕ → N ϕ for all ϕ ∈ L 2 (M ).

In the next theorem we improve this convergence result substantially, proving N phys P δ → N uniformly. The physical incident fields can be represented through a Herglotz operator H phys : L 2 (M ) → L 2 (D) with H phys ϕ defined as the restriction of SL M ϕ to D. Then we see, as in Section 4 for N , that N phys satisfies the factorization

N phys = H * T H phys .
(5.1)

The next result shows that N phys P δ → N in norm. Hence, the eigenvalues and eigenspaces of N phys P δ converge to those of N . Therefore it makes sense to apply (a regularized version of) the series criterion (4.6) to N phys P δ for δ > 0 small enough. We skip details on this regularization procedure and refer to [START_REF] Lechleiter | A regularization technique for the factorization method[END_REF]. Theorem 5.3. N phys P δ converges to N in the operator norm of L 2 (M ). Proof. From Lemma 5.2 we note that N phys P δ ϕ → N ϕ for all ϕ ∈ L 2 (M ) and hence P * δ N * phys ⇀ N * converges weakly in L 2 (M ) as δ → 0. Since (δ + V * V ) -1 is selfadjoint, the representation P * δ = Ṽ * V (δ + V * V ) -1 holds. Using a singular system of the compact operator V , one computes that V (δ + V * V ) -1 = (δ + V V * ) -1 V . Consequently, from (5.1) we obtain

P * δ N * phys = Ṽ * (δ + V V * ) -1 V H * phys T * H . (5.2)
Define the operator H Σ : H -1/2 (Σ) → L 2 (D) by H Σ ϕ = SL Σ ϕ| D . Also, note that our assumption that k 2 is not an eigenvalue of -∆ in B 0 implies that S Σ is boundedly invertible.

The following observation is crucial: The operator H * phys : L 2 (D) → L 2 (M ) can be factorized as

H * phys = V * S -1 * Σ H * Σ , or equivalently, H phys = H Σ S -1 Σ V .
To show this equality we choose ϕ ∈ L 2 (M ) and set ψ = V ϕ. Since SL Σ S -1

Σ ψ Σ = ψ we find that the radiating function SL Σ S -1

Σ ψ equals V ϕ on Σ, and hence SL M ϕ = SL Σ S -1 Σ ψ on Σ. Exploiting again the fact that k 2 is not a Dirichlet eigenvalue in B 0 , both potentials necessarily equal each other in B 0 . Thus, their evaluations on D equal, too. This implies H phys = H Σ S -1 Σ V . From the factorization of H * phys together with (5.2) we have [START_REF] Kress | Linear integral equations[END_REF] implies that P * δ N * phys converges in norm. As P * δ N * phys ϕ ⇀ N * ϕ weakly and as weak and strong limit coincide, the limit of P * δ N * phys is N * . 6. Numerical Examples. In this section we present numerical reconstructions of penetrable inclusions in planar waveguides. Before going into details concerning the inverse problem, let us briefly indicate that we numerically solve the direct medium scattering problem (2.13) using a volumetric integral equation approach that is constructed along the lines of Vainikko's method [START_REF] Vainikko | Fast solvers of the Lippmann-Schwinger equation[END_REF]. After suitable periodization, the integral operator V diagonalizes on trigonometric polynomials and the corresponding eigenvalues can be computed analytically in terms of Bessel functions. Therefore one can implement the application of the discrete periodized Lippmann-Schwinger integral operator in O(N log(N )) operations, where N is the number of unknowns. Being able to compute matrix-vector multiplications with the underlying system matrix, one then applies an iterative solver to the linear system. For further details see [START_REF] Lechleiter | Spectral volumetric integral equation methods for acoustic medium scattering in a planar homogeneous 3D wave guide[END_REF].

P * δ N * phys = Ṽ * (δ + V V * ) -1 V V * S -1 * Σ H * Σ T * H . (δ + V V * ) -1 V V * converges pointwise to the identity. Since S -1 * Σ H * Σ S -1 * H is com- pact, Theorem 10.6 in
For our numerical reconstructions we use two penetrable inclusions that we call the submarine and the double cylinder. The submarine consists of four intersecting (oblate) ellipsoids of revolution that model the main body, tower, and two fins of a submarine. The refractive index is piecewise constant and partially absorbing in a layer next to the submarine's boundary. It equals 4 + 2i in the absorbing parts while in the non-absorbing part it equals 6. Figure 6(a) shows a plot of the submarine inside a guide of height 0.5. The double cylinder consists of two upright cylindrical inhomogeneities with constant refractive index equal to 3 and 6, see Figure 6(c). Both cylinders are 0.15 high with radius 0.125 and are centered at (0.125, 0.125, 0.075) and (-0.125, -0.125, 0.075).

Our first numerical experiment provides a proof of concept for the factorization method (Theorem 4.5) and the theory on the projection operator P δ in combination with data from physical sources (Theorem 5.3). We use the submarine scatterer, the wave length equals π/10 ≈ 0.31, and we use 470 incident point sources on a cylinder of radius 1.5 surrounding the object. The measurement cylinder is about 4 (a) (b) wave lengths away from the submarine, thus, information contained in the evanescent modes is not measured reliably (even without artificial noise added to the synthetic data). Figure 6.2 shows reconstructions by the Factorization method from conjugate incident point sources (N ), from physical incident point sources (N phys ), and from physical incident point sources combined with the projection operator (N phys P δ ). For the plots in the first row no artificial noise has been added to the data. For the plots in the second row uniformly distributed random noise has been added elementwise to the measurements, and the relative noise level equals one percent in the spectral norm. For the plots in the third row, the relative noise level equals 5 percent. The reconstructions provided in Figure 6.2 show that without noise several details of the submarine's shape can be well reconstructed from incident conjugate point sources. The fins of the submarine are too small to be reconstructed, but at least the reconstruction allows to distinguish front and back of the submarine. For most of the reconstructions, the scatterer throws a shadow towards the Dirichlet waveguide boundary. For physical incident point sources the quality of the reconstructions is substantially worse while the reconstructions using N phys P δ are of the same quality as those using N . This observation can be explained by considering the discrepancy N phys P δ -N / N that is about 10 -4 for small δ > 0, see Figure 6.3(a). The convergence result from Theorem 6.3 is clearly validated numerically by the error curve corresponding to noise-free data in Figure 6.3(a). When adding five percent of relative noise to N and Ñ , N phys P δ -N / N tends roughly to 10 -1 as δ → 0.

In our second experiment we restrict the horizontal aperture of the measurements to 45 degrees, see Figure 6.3(b). Apart from the restriction of aperture, the setting is the same as in the above-described first experiment. For this experiment the application of the projector P δ did not improve the reconstructions for physical incident point sources, possibly due to the small number of source points, thus we do not show reconstructions involving P δ . We note that the reconstructions using physical sources do merely allow a rough estimate of position and shape of the submarine. As in the first experiment, the reconstructions for conjugate sources are better, especially, size and position of the submarine seem to be more stably reconstructed.

In our third experiment we consider the double cylinder and we are interested in separating the two objects. We consider two wave lengths λ 1 = 2π/10 ≈ 0.63 and λ 2 = 2π/15 ≈ 0.42 and use 45 conjugate point sources on a cylinder of radius two. Figure 6.5(a) shows that the Factorization method separates the two cylinders using the data taken at small wave length λ 1 when no artificial noise has been added to that data. However, the separation fails when we add one percent of relative noise to the data, see Figure 6.5(b) and (c). For the data taken at large wavelength λ 2 one can add substantially more noise to the data while still preserving separability. In all reconstructions of this experiment the cylinder where the contrast q equals 5 appears to be larger compared to the one where q = 2.

Γ

  + := {x ∈ R 3 : x 3 = h} and Γ -:= {x ∈ R 3 : x 3 = 0}, u is supposed to satisfy a Neumann and a Dirichlet boundary condition, respectively, ∂u ∂x 3 = 0 on Γ + , and u = 0 on Γ -. (1.2)

(

  B1) Either for each closed ball B ⊂ D there is C B > 0 with |q| ≥ C B a.e. in B, or (B2) there is ǫ > 0 such that Dǫ dx |q| < ∞, where D ǫ = {x ∈ D : dist(x, ∂D) < ǫ}.

2 ,

 2 Proof ofLemma 4.3]. Hence the first claim of the proposition follows and we conclude that T is a Fredholm operator of index 0. Showing injectivity of T in the next step hence also implies that T is an isomorphism and finishes the proof. Assume that T f = 0 for some f ∈ L 2 (D), that is, f + |q|v = 0 in D. We plug this relation in the equation ∆v + k 2 (1 + q)v = -k 2 qf / |q| which defines v and find that ∆v + k 2 v = 0 in Ω. However, since v is radiating we conclude by uniqueness of solution for the direct problem that v vanishes in Ω. In turn f vanishes in D and hence T is injective.The next theorem characterizes points of D by the range of H * . Theorem 4.3. For z ∈ D, we defineφ z = G(•, z)| D ∈ L 2 (D).Under assumption (B1) or (B2), it holds that z ∈ D if and only if φ z ∈ Rg(H * ).

  (d) T is injective. Then the operator F ♯ := | Re F | + Im F is positive, self-adjoint and compact, and the ranges of H * : X * → Z and F 1/2 ♯ : Z → Z coincide. Here | Re F | is defined in the standard way via the spectral representation for compact self-adjoint operators.

Theorem 4 . 5 .

 45 For z ∈ D define φ z := G(•, z)| M ∈ L 2 (M ). Then z ∈ D if, and only if, φ z ∈ Rg N 1/2 ♯

Figure 5 . 1 .

 51 Figure 5.1. Sketch of the geometry of the auxiliary domain B 0 with boundary Σ in a horizontal slice. B 0 includes the scatterer D and excludes the measurement surface M .

Figure 6 . 1 .

 61 Figure 6.1. Two scattering objects. (a) The submarine. (b) The double cylinder.

Figure 6 . 2 .

 62 Figure 6.2. Experiment 1: Reconstructions of the submarine for different types of incident fields. Blue lines in the plots indicates the wave length. (a) Conjugate incident point sources, no noise. (b) Physical incident point sources, no noise. (c) Physical incident point sources, the projection operator P 10 -12 has been used, no noise. (d) Conjugate incident point sources, 1% relative noise. (e) Physical incident point sources, 1% relative noise. (f ) Physical incident point sources, the projection operator P 10 -12 has been used, 1% relative noise. (g) Conjugate incident point sources, 5% relative noise. (h) Physical incident point sources, 5% relative noise. (i) Physical incident point sources, the projection operator P 10 -12 has been used, 5% relative noise.

Figure 6 . 3 .

 63 Figure 6.3. (a) The error N phys P δ -N / N versus δ, in the setting used for the reconstructions shown in Figure 6.2. Norms of matrices are measured in the spectral norm. Circles (o) correspond to errors computed without adding artificial noise to N and N phys . Pluses (+) correspond to errors computed from noisy versions of N and N phys at a noise level of 5%. (b) Horizontal positions of the sources (red circles) for experiment 2.

Figure 6 . 4 .

 64 Figure 6.4. Experiment 2: Reconstructions of the submarine from limited aperture data (1/8 horizontal aperture). Blue lines in the plots indicates the wave length. (a) Conjugate incident point sources, no noise. (b) Conjugate incident point sources, 1% relative noise. (c) Conjugate incident point sources, 5% relative noise. (d) Physical incident point sources, no noise. (e) Physical incident point sources, 1% relative noise. (f ) Physical incident point sources, 5% relative noise.

Figure 6 .

 6 Figure 6.5(e) and (f) show reconstructions for a relative noise level of 5 percent that clearly separate the two objects.In all reconstructions of this experiment the cylinder where the contrast q equals 5 appears to be larger compared to the one where q = 2.

Figure 6 . 5 .

 65 Figure 6.5. Experiment 3: Reconstructions of the double cylinder for two wave lengths λ 1 = 2π/10 and λ 2 = 2π/15 (indicated by blue lines in the plots), using conjugate incident point sources. (a) Large wave length λ 1 , no noise. (b) Wave length λ 1 , 1% relative noise. (c) Wave length λ 1 , slice plot of reconstruction in (b) along {x 1 = x 2 }. (d) Short wave length λ 2 , no noise. (e) Wave length λ 2 , 5% relative noise. (f ) Wave length λ 2 , slice plot of reconstruction in (e) along {x 1 = x 2 }.