Precisely, we can reconstruct the obstacle by the eigenvalues and eigen- functions of the far-field operator

ANISOTROPIC INHOMOGENEOUS BACKGROUND

PU-ZHAO KOW AND JENN-NAN WANG

Abstract. In this paper, we are interested in the problem of determining the shape of sound-soft or impedance obstacles in the acoustic wave scattering with an anisotropic inhomogeneous medium. The main theme is to extend the factorization method to this setting. Precisely, we can reconstruct the obstacle by the eigenvalues and eigen- functions of the far-field operator. We also provide some numerical simulations in the paper.

1. Introduction

Let an obstacle be embedded inside an anisotropic inhomogeneous background. This system causes the acoustic wave to scatter. In this work, we are interested in re- constructing the shape of the obstacle by the far-field information assuming that the inhomogeneity is known. There are several existing methods for other types of inverse scattering problems. Here we will focus on the factorization method, which was first introduced by Kirsch in [Ki98, Ki99] to treat the inverse scattering problems for the scalar Helmholtz equation. For the detailed development of the factorization method for the Helmholtz and time-harmonic Maxwell equations, we refer the reader to the monograph [KG08]. The factorization method shares the same spirit as the linear sam- pling method, where both methods try to reconstruct the shape of the scatterer (an obstacle or the support of the inhomogeneity) by determining whether a point is inside or outside the scatterer.

We now briefly describe the problem considered in the paper. Let A(x) = (aij(x)) be a real-symmetric matrix with C^∞ entries, which satisfies the following uniform elliptic condition: there exists a constant 0 < c < 1 such that

c|ξ|² ≤X

ij

a_ij(x)ξ_iξ_j ≤ c⁻¹|ξ|²

for all x ∈ R^d, ξ ∈ R^d, d = 2 or 3. Moreover, we assume that supp(A − I) is compact in R^d. Suppose that the acoustic refraction index n ∈ L^∞(R^d) with n ≥ c, such that supp(n − 1) is compact in R^d. Let B_R be the ball centered at origin with radius R > 0 which contains supp(A − I) and supp(n − 1). Let D be an open bounded domain where ∂D is Lipschitz, D ⊂ BR, and R^d\ D is connected. Let u^to_D(x, ˆz) = u^inc_A,n(x, ˆz) + u^sc_D(x, ˆz) with an appropriate incident field u^inc_A,n(x, ˆz) and the corresponding scattered field u^sc_D(x, ˆz), ˆz = z/|z|, satisfy the following acoustic equation







∇ · (A(x)∇u^to_D) + k²n(x)u^to_D = 0 in R^d\ D,

Bu^to_D = 0 on ∂D,

u^sc_D satisfies the Sommerfeld radiation condition at |x| → ∞,

(1.1)

1

where Bu^to_D is either Dirichlet or impedance boundary conditions. Let u^∞_D(ˆx, ˆz) be the far-field pattern of u^sc_D. The inverse obstacle problem is to determine D from the knowledge of u^∞_D(ˆx, ˆz) for all ˆx, ˆz ∈ S^d−1.

This type of inverse obstacle problem including theoretical and numerical investiga- tions has been studied extensively. We will not try to exhaust all the related works here. Instead, we refer to recent monographs [CCH16, KG08] and references therein for the detailed development of the problems and reconstruction methods. To put our problem in perspective, we mention several closely related results.

• Penetrable obstacle in an inhomogeneous background. In [KP98], the authors considered the acoustic scattering from a penetrable obstacle inside an inhomo- geneous structure. In their work, they showed that, if the reference medium is known, then the penetrable obstacle can be uniquely determined by the far-field data at a fixed energy. The approach used in [KP98] was only for the uniqueness proof. Later, in [GMMR12], the authors designed a reconstruction algorithm in the spirit of the factorization method. For the case of anisotropic penetrable ob- stacle embedded inside of a homogeneous medium, a factorization method was developed in [KL13]. On the other hand, when there was an unknown cavity hidden inside of the penetrable obstacle, a reconstruction result based on the factorization method was investigated in [YZZ13].

• Impenetrable obstacle in a homogeneous background. A factorization method for reconstructing an impenetrable obstacle in a homogeneous medium (Helmholtz equation) using the spectral data of the far-field operator was developed in [Ki98]

(see also [KG08, Chapter 1]).

• Impenetrable obstacle in an inhomogeneous background. When the inhomoge- neous background was described the Schrödinger equation with an inhomoge- neous potential function, the theoretical study of the factorization method was carried out in [NPT07]. Moreover, the obstacle can be reconstructed without knowing the boundary condition on the boundary of the obstacle (either Dirich- let or Robin boundary conditions). Our work is closely related to this paper.

In our paper, we consider an impenetrable obstacle in an anisotropic inhomoge- neous background.

We would like to remark that, in [GMMR12, KG08, KP98], the far-field operator is defined in terms of the far-field pattern induced by the incident plane waves. In our case, we take a different incident field. Roughly speaking, the incident field used in our approach is the far-field pattern of the outgoing fundamental solution for the background equation. Indeed, in the case where the reference medium is homogeneous, the far-field pattern of the standard outgoing fundamental solution is simply the plane waves. Since we assume that the background medium is known, the far-field pattern of the outgoing fundamental solution for the background equation can be determined.

However, from the numerics perspective, the computation of the far-field of the outgoing fundamental solution directly is not efficient. Luckily, thanks to the mixed reciprocity relation (see Lemma 2.2), we can compute the far-field of the outgoing fundamental solution by solving the total field of the background equation with the incident plane waves.

Another focal point of this work is the numerical demonstration of increasing stability phenomenon with respect to the wave number k in the reconstruction of the obstacle.

It is known that most of inverse problems are ill-posed including the problem we discuss here. However, in the case of identifying the potential in the Schrödinger equation by the Dirichlet-to-Neumann map, the stability increases as we increase the wave number (the logarithmic part decreases as the wave number increases) [INUW14]. A similar property was also established for the determination of the refraction index in the acous- tic equation [NUW13]. We need to point out a major difference between [INUW14] and [NUW13]. In [INUW14], the constant appears in the Hölder part of the stability es- timate depends polynomially in k. Hence, the stability of determining the potential in the Schrödinger equation tends a Hölder type as k increases. While, in [NUW13], the constant in the Hölder part of the stability estimate depends exponentially in k.

Consequently, the stability estimate in the acoustic case remains a logarithmic type when k is too large.

Considering the discussion above, even though the inverse obstacle problems for the Schrödinger equation and for the acoustic equation at a fixed wave number are mathematically equivalent, it is necessary in practice to study the acoustic equation separately at least from the viewpoint of stability. We believe that it is an important and challenging problem to solve the inverse obstacle problem in the acoustic equation stably at high wave numbers.

The paper is organized as follows. In Section 2, we first review some properties of the outgoing fundamental solution for the background equation and its far-field pattern.

In Section 3, we prove our main reconstruction theorem for the sound-soft (Dirichlet) obstacle. The same ideas can be generalized to the obstacle with Neumann or impedance boundary conditions. Thus, in Section4, we state the reconstruction theorems without proofs corresponding to the cases of Neumann and impedance boundary conditions.

In Section 5, some numerical simulations are presented to show the efficiency of our method. Finally, in Section 6, we compare the reconstruction results for the acoustic and the Schrödinger equations at different wave numbers to demonstrate the increasing stability phenomena.

2. Fundamental solution and Herglotz wave functions

In this section, we discuss some properties of the outgoing fundamental solution to the background equation (that is, the reference medium).

Definition 2.1. We say that a function u ∈ H_loc¹ (R^d) satisfies the Sommerfeld radiation condition, if

∂_ru(x) − iku(x) = o(|x|^−d−12 ) as |x| → ∞, where ∂_r= ˆx · ∇.

For each x ∈ R^d, let ΦA,n(z, x) (for x 6= z) be the outgoing fundamental solution of the following anisotropic inhomogeneous acoustic equation

(∇_z· (A(z)∇_zΦ_A,n(z, x)) + k²n(z)Φ_A,n(z, x) = δ(z − x) ∀ z ∈ R^d, Φ_A,n(z, x) satisfies the Sommerfeld radiation condition at |z| → ∞.

Let u^inc_ref(x, ˆz) = e^ikx·ˆz be an incident plane wave. The presence of inhomogeneity (with- out the obstacle) in B_R gives rise to a unique scattered field u^sc_ref ∈ H_loc¹ (R^d) which

satisfies the Sommerfeld radiation condition. Then the total field u^to_ref = u^sc_ref + u^inc_ref satisfies

∇ · (A(x)∇u^to_ref) + k²n(x)u^to_ref = 0, ∀ x ∈ R^d. (2.1) We first prove the following mixed reciprocity principle.

Lemma 2.2. Let Φ^∞_A,n(ˆz, x) be the far-field pattern of Φ_A,n(z, x). Then we have Φ^∞_A,n(ˆz, x) = u^to_ref(x, −ˆz) ∀ x ∈ R^d, ˆz ∈ S^d−1,

where u^to_ref(·, −ˆz) is given in (2.1).

Proof. Here we adopt the proof from [KG08]. Let Φ(z, x) = Φ_I,1(z, x) be the outgoing fundamental solution of the standard Helmholtz equation. Define

Ψ(z, x) := Φ_A,n(z, x) − Φ(z, x).

Then Ψ(z, x) is a smooth at all z ∈ R^d and

∆_zΨ(z, x) + k²Ψ(z, x)

= −(∇_z· ((A(z) − I)∇_zΦ_A,n(z, x)) + k²(n(z) − 1)Φ_A,n(z, x)).

We can write Ψ(z, x)

= − Z

R^d

Φ(z, y)(∇_y · ((A(y) − I)∇_yΦ_A,n(y, x)) + k²(n(y) − 1)Φ_A,n(y, x)) dy.

Since both A − I and n − 1 have compact supports, the integral equation above is well- defined. Note that the far-field pattern of Φ(z, y) is Φ^∞(ˆz, y) = e^−ikˆz·y = u^inc_ref(y, −ˆz).

We then see that the far-field pattern of Ψ(z, x) is Ψ^∞(ˆz, x)

= − Z

R^d

e^−ikˆz·y(∇_y · ((A(y) − I)∇_yΦ_A,n(y, x)) + k²(n(y) − 1)Φ_A,n(y, x)) dy. (2.2) Since u^to_ref(x, −ˆz) satisfies

∇_x· (A(x)∇_xu^to_ref(x, −ˆz)) + k²n(x)u^to_ref(x, −ˆz) = 0 for all x ∈ R^d, we have

∇_x· (A(x)∇_xu^sc_ref(x, −ˆz)) + k²n(x)u^sc_ref(x, −ˆz)

= − (∇_x· (A(x)∇_xu^inc_ref(x, −ˆz)) + k²n(x)u^inc_ref(x, −ˆz))

= − (∇_x· ((A(x) − I)∇_xu^inc_ref(x, −ˆz)) + k²(n(x) − 1)u^inc_ref(x, −ˆz)).

In view of the symmetry property ΦA,n(x, y) = ΦA,n(y, x), using integration by parts, we can derive

u^sc_ref(x, −ˆz)

= − Z

R^d

Φ_A,n(y, x)(∇_y· ((A(y) − I)∇_yu^inc_ref(y, −ˆz)) + k²(n(y) − 1)u^inc_ref(y, −ˆz)) dy

= − Z

R^d

e^−ikˆz·y(∇_y · ((A(y) − I)∇_yΦ_A,n(y, x)) + k²(n(y) − 1)Φ_A,n(y, x)) dy. (2.3)

Combining (2.2) and (2.3) implies that

Ψ^∞(ˆz, x) = u^sc_ref(x, −ˆz), and hence

Φ^∞_A,n(ˆz, x) = Ψ^∞(ˆz, x) + Φ^∞(ˆz, x) = u^sc_ref(x, −ˆz) + u^inc_ref(x, −ˆz) = u^to_ref(x, −ˆz),

which is our desired lemma. 

Definition 2.3. For each g ∈ L²(S^d−1), we define the modified Herglotz wave function as

u_A,n,g(x) :=

Z

S^d−1

Φ^∞_A,n(ˆz, x)g(ˆz) ds(ˆz) for all x ∈ R^d. It is well-known that

g 7→

Z

S^d−1

e^ikx·ˆzg(ˆz) ds(ˆz) = Z

S^d−1

u^inc_ref(x, ˆz)g(ˆz) ds(ˆz)

is injective. Since each u^inc_ref(x, ˆz) leads a unique scattered field u^sc_ref(x, ˆz), and thus induces a unique total field u^to_ref(x, ˆz). Combining these discussions and by the super- position property, we can prove the following lemma.

Lemma 2.4. For coefficients A and refraction index n prescribed above, the mapping g 7→ u_A,n,g is injective.

3. Reconstruction of sound-soft obstacle

By Lemma 2.2, we can choose the incident field u^inc_A,n(x, ˆz) = Φ^∞_A,n(ˆz, x) instead of the usual plane wave u^inc_ref(x, ˆz) = e^ikx·ˆz. By choosing this incident field, no scattering occurs if there is no obstacle. Now let u^sc_Dir(x, ˆz) be the corresponding scattered field, then the total field u^to_Dir(x, ˆz) = u^inc_A,n(x, ˆz) + u^sc_Dir(x, ˆz) satisfies







∇ · (A(x)∇u^to_Dir) + k²n(x)u^to_Dir = 0 in R^d\ D,

u^to_Dir = 0 on ∂D,

u^sc_Dir satisfies Sommerfeld radiation condition at |x| → ∞,

(3.1)

that is,







∇ · (A(x)∇u^sc_Dir) + k²n(x)u^sc_Dir = 0 in R^d\ D,

u^sc_Dir = −u^inc_A,n on ∂D,

u^sc_Dir satisfies Sommerfeld radiation condition at |x| → ∞.

(3.2)

Denote u^∞_Dir be the far-field pattern of u^sc_Dir, which satisfies u^sc_Dir(x) = γ_d e^ik|x|

|x|^d−12 u^∞_Dir(ˆx) + O(|x|^−d+12 ) as |x| → ∞ uniformly for all directions ˆx = x/|x|, where

γ_d :=

(e^iπ/4/√

8πk if d = 2,

1/4π if d = 3, (3.3)

see e.g. [GMMR12].

The well-posedness of u^sc_D is guaranteed by the following well-known fact:

Lemma 3.1. Given any f ∈ H^1/2(∂D), there exists a unique v ∈ H_loc¹ (R^d\ D) such

that 





∇ · (A(x)∇v) + k²n(x)v = 0 in R^d\ D,

v = f on ∂D,

v satisfies Sommerfeld radiation condition at |x| → ∞.

(3.4)

3.1. The factorization of the far-field operator. For each g ∈ L²(S^d−1), recall the modified Herglotz wave function is defined by

uA,n,g(x) :=

Z

S^d−1

Φ^∞_A,n(ˆz, x)g(ˆz) ds(ˆz)

= Z

S^d−1

u^inc_A,n(x, ˆz)g(ˆz) ds(ˆz), (3.5) which is simply the superposition of the incident fields. So, the far-field of the scattered field is the superposition of u^∞_Dir. Thus, we can define the far-field operator as follows:

Definition 3.2. The far-field operator F_Dir: L²(S²) → L²(S²) is given by F_Dir(g)(ˆx) =

Z

S^d−1

u^∞_Dir(ˆx, ˆz)g(ˆz) ds(ˆz) for all ˆx ∈ S^d−1.

With Lemma 3.1, we also define an operator closely related to the far-field operator:

Definition 3.3. The data-to-pattern operator G_Dir : H^1/2(∂D) → L²(S^d−1) is defined by G_Dirf = v^∞, where v^∞ is the far-field of v ∈ H_loc¹ (R^d\ D) which satisfies (3.4).

To discuss the relation between the far-field operator and the data-to-pattern oper- ator, we introduce the single-layer potential operator S defined by

(Sφ)(x) :=

Z

∂D

ΦA,n(x, z)φ(z) ds(z) for all x ∈ ∂D.

We first establish the following mapping property of S.

Lemma 3.4. The potential S : H^−1/2(∂D) → H^1/2(∂D) is continuous.

Proof. Let ΦA,1(x, z) be the outgoing fundamental solution of

(∇x· (A(x)∇xΦA,1(x, z)) + k²ΦA,1(x, z) = δ(x − z) ∀ x ∈ R^d, Φ_A,1(x, z) satisfies the Sommerfeld radiation condition at |x| → ∞.

Since the coefficient A is smooth, Φ_A,1(x, z) is C^∞ for x 6= z. By the techniques of pseudodifferential operators (see [McL00, Theorem 6.11]), let

( ˜Sφ)(x) :=

Z

∂D

ΦA,1(x, z)φ(z) ds(z) for all x ∈ ∂D, then ˜S : H^−1/2(∂D) → H^1/2(∂D) is continuous. Note that

∇_x· (A(x)∇_x(Φ_A,n− Φ_A,1)) + k²n(x)(Φ_A,n− Φ_A,1) = k²(1 − n)Φ_A,1.

Observe that k²(1−n)ΦA,1 ∈ L²(R^d) (uniformly in z). By the elliptic estimate theorem, we have that

(Φ_A,n− Φ_A,1)(·, z) ∈ H^2,2(D).

From the trace theorem and the symmetry of the fundamental solutions, it follows that

|∇_x(Φ_A,n− Φ_A,1)(x, z)| = |∇_z(Φ_A,n− Φ_A,1)(x, z)| ∈ H^1/2(∂D) uniformly for x ∈ ∂D. Therefore, we can show that the operator

S⁰ : φ 7→

Z

∂D

(Φ_A,n− Φ_A,1)(x, z)φ(z)ds(z) (3.6) maps H^−1/2(∂D) to H¹(∂D), in particular, S⁰ : H^−1/2(∂D) → H^1/2(∂D) is continuous.

Since S = ˜S + S⁰, Lemma 3.4 follows.

 Following the ideas in [KG08, Theorem 1.15], we can obtain the following proposition, which contains the main idea of this work.

Proposition 3.5. Let G^∗_Dir : L²(S^d−1) → H^−1/2(∂D) be the adjoint of G_Dir and S^∗ : H^−1/2(∂D) → H^1/2(∂D) be the adjoint of S. Then

F_Dir = −G_DirS^∗G^∗_Dir. (3.7) Proof. We consider an auxiliary operator H_Dir : L²(S^d−1) → H^1/2(∂D) defined as the restriction of the modified Herglotz wave function on ∂D, that is,

(H_Dirg)(x) :=

Z

S^d−1

Φ^∞_A,n(ˆz, x)g(ˆz) ds(ˆz) for all x ∈ ∂D.

From (3.2), we know that

FDir = −GDirHDir. (3.8)

Also, it is not difficult to check that the adjoint H_Dir^∗ : H^−1/2(∂D) → L²(S^d−1) of H_Dir is given by

(H_Dir^∗ φ)(ˆz) = Z

∂D

Φ^∞_A,n(ˆz, x)φ(x) dx

for all φ ∈ H^−1/2(∂D). Observe that (H_Dir^∗ φ)(ˆz) is nothing but the far-field pattern of v(z) =

Z

∂D

ΦA,n(z, x)φ(x) dx, which shows that H_Dir^∗ = G_DirS, that is,

H_Dir = S^∗G^∗_Dir. (3.9)

Putting together (3.8) and (3.9) implies (3.7). 

Similar to [KG08, Theorem 1.12], we can show that

Proposition 3.6. Φ^∞_A,n(·, z) ∈ range(G_Dir) if and only if z ∈ D.

Proof. First of all, we assume that Φ^∞_A,n(·, z) ∈ range(G_Dir), i.e., there exists f ∈ H^1/2(∂D) such that G_Dirf = Φ^∞_A,n(·, z). Let w ∈ H_loc¹ (R^d \ D) be the unique solu- tion to (3.4). Since both

w and Φ_A,n(·, z) have the same far-field pattern Φ^∞_A,n(·, z), by Rellich’s lemma, there exists r0 > 0 such that

w(x) = Φ_A,n(x, z) for |x| > r₀.

Since R^d\ D is connected, and both w and Φ_A,n satisfy the anisotropic inhomogeneous acoustic equation (2.1) in R^d\ D, by the unique continuation property of this equation, we have that

w(x) = Φ_A,n(x, z) for all x ∈ R^d\ D.

Note that w ∈ H_loc¹ (R^d\ D). In other words, w has no singularity in R^d\ D and thus z ∈ D.

Conversely, let z ∈ D. Define w = Φ_A,n(·, z) in R^d\ D and f = w|_∂D, we can easily see that GDirf = Φ^∞_A,n(·, z), that is, Φ^∞_A,n(·, z) ∈ range(GDir).  With Proposition 3.5 and 3.6 at hand, it is helpful to study the operators S, G_Dir, and F_Dir.

3.2. The properties of the single-layer potential S. In this section, we shall study the properties of the far-field operator. With the help of [McL00, Theorem 7.5], we can further improve Lemma3.4.

Lemma 3.7. If k² is not an eigenvalue of

(∇ · (A(x)∇u) + k²n(x)u = 0 in D,

u = 0 on ∂D, (3.10)

then S : H^−1/2(∂D) → H^1/2(∂D) is an isomorphism.

Next, by mimicking the proofs of [KG08, Lemma 1.14(b),(c),(d)], we can obtain the following lemmas. Since we are considering a different equation, we include the proofs here for the sake of completeness.

Lemma 3.8. If k² is not an eigenvalue of (3.10), then =hφ, Sφi ≤ 0 for all φ ∈ H^1/2(∂D), where h·, ·i is the H^−1/2(∂D) × H^1/2(∂D) duality pair and = denotes the imaginary part of a complex number. Moreover,

=hφ, Sφi = 0 if and only if φ ≡ 0.

Proof. For each φ ∈ H^−1/2(∂D), define v(x) = (Sφ)(x) =

Z

∂D

Φ_A,n(x, y)φ(y) ds(y) for x ∈ R³\ ∂D.

Then, by [McL00, Theorem 6.11], we have v ∈ H¹(D) ∩ H_loc¹ (R³\ D) and the traces on

∂D with the outer normal normal ν(x), v_±(x) = lim

h→0+

v(x ± hν(x)), ν(x) · A(x)∇v_±(x) = lim

h→0+

ν(x) · A(x ± hν(x))∇v(x ± hν(x)), exist and satisfy the jump condition

v = v± = Sφ for all φ ∈ H^−1/2(∂D).

Choose φ = ν · A∇v−− ν · A∇v₊, then hφ, Sφi = hν · A∇v−− ν · A∇v₊, vi

= Z

∂D

(ν · A∇v−)v ds − Z

∂D

(ν · A∇v₊)v ds

= Z

∂D

(ν · A∇v−)v ds + Z

∂(BR\D)

(ν · A∇v₊)v ds − Z

|x|=R

(ν · ∇v)v ds

= Z

D∪(BR\D)

(∇ · (A∇v))v dx + Z

D∪(BR\D)

∇¯v · A∇v dx − Z

|x|=R

(∂rv)v ds.

Since ∇ · (A∇v) + k²n(x)v = 0 in R³ \ ∂D, together with the Sommerfeld radiation condition, then

hφ, Sφi = Z

D∪(BR\D)

(∇¯v · A∇v − k²n(x)|v|²) dx − Z

|x|=R

(∂_rv)v ds

= Z

D∪(BR\D)

(∇¯v · A∇v − k²n(x)|v|²) dx − ik Z

|x|=R

|v|²ds + o(R^−d−12 )

(3.11)

as R → ∞. Then

=hφ, Sφi = −k Z

|x|=R

|v|²ds + o(R^−d−12 ), and thus

=hφ, Sφi = −k lim

R→∞

Z

|x|=R

|v|²ds = −k lim

R→∞

Z

S^d−1

|R^(d−1)/2v|²ds(ˆx)

= −|γd|²k Z

S^d−1

|v^∞|²ds(ˆx) ≤ 0,

(3.12)

where γ_d is given in (3.3).

Finally, we want to show that =hφ, Sφi = 0 leads to φ ≡ 0. If =hφ, Sφi = 0, then v^∞ ≡ 0 from (3.12). Since R³ \ D is connected, from Rellich’s lemma and unique continuation property, we know that v = 0 in R³\ D. Therefore, by the trace theorem, we have Sφ = 0. By Lemma3.7, we know that S is an isomorphism and thus φ ≡ 0.  Lemma 3.9. Let S_i be the single-layer operator corresponding to the wave number k = i = √

−1. Then Si : H^−1/2(∂D) → H^1/2(∂D) is self-adjoint and coercive. Moreover, S − S_i : H^−1/2(∂D) → H^1/2(∂D) is compact.

Proof. Substituting k = i into (3.11) gives hφ, S_iφi =

Z

D∪(BR\D)

(∇¯v · A∇v + n(x)|v|²) dx + Z

|x|=R

|v|²ds + o(R^−d−12 ) as R → ∞. By taking R → ∞, we reach

hφ, S_iφi = Z

R³\∂D

(∇¯v · A∇v + n(x)|v|²) dx, which shows that Si is self-adjoint.

Since A is uniform elliptic and n(x) ≥ c, then

hφ, Siφi ≥ ckvk²_H1(R³\∂D).

By trace theorem, together with Lemma3.7, there exist constants c₁, c₂ with 0 < c₂ <

c₁ < c such that

hφ, S_iφi ≥ c₁kvk²_H1/2(∂D) = c₁kS_iφk²_H1/2(∂D) ≥ c₂kφk²_H−1/2(∂D),

that is, S_i is coercive. Note that the operator S − S_i has the same mapping property as S⁰ in (3.6). We thus immediately conclude that

S − S_i : H^−1/2(∂D) → H^1/2(∂D)

is compact since the embedding H¹(∂D) into H^1/2(∂D) is compact.  3.3. The properties of the data-to-pattern operator GDir. We can prove following two lemmas about G_Dir.

Lemma 3.10. G_Dir : H^1/2(∂D) → L²(S^d−1) is injective.

Proof. This is a direct consequence of Rellich’s lemma and unique continuation property.

 Lemma 3.11. G_Dir : H^1/2(∂D) → L²(S^d−1) is compact.

Proof. This lemma can be proved following the line of [KG08, Lemma 1.13]. Let R > 1 be a large number. For each f ∈ H^1/2(∂D), let v ∈ H_loc¹ (R^d\ D) satisfy (3.4). By the representation formula

v^∞(ˆx) = Z

∂D



(ν · A(z)∇Φ^∞_A,n(ˆx, z))v(z) − Φ^∞_A,n(ˆx, z)(ν · A(z)∇v(z))



ds(z) (3.13) for ˆx ∈ S^d−1, we can decompose GDir into GDir = G2◦ G1, where

G₁ : H^1/2(∂D) → C(∂B_R) × C(∂B_R) G₂ : C(∂B_R) × C(∂B_R) → L²(S^d−1) are defined by

G₁f =



v|_∂BR, ν · A∇v|_∂BR

 , G₂(g, h) =

Z

∂BR



(ν · A∇Φ^∞_A,n(ˆx, z))g(z) − Φ^∞_A,n(ˆx, z)h(z)

 ds(z).

By interior regularity results for elliptic equation, we know that G1 : H^1/2(∂D) → C(∂B_R) × C(∂B_R) is bounded. Moreover, from the analyticity of Φ^∞_A,n(ˆx, z) in ˆx, it immediately follows that G2 : C(∂BR) × C(∂BR) → L²(S^d−1) is compact and the proof

is completed. 

Modify the ideas in [AC06, Theorem 3.1], we can establish the denseness of range(GDir).

Lemma 3.12. ran(GDir) is dense in L²(S^d−1).

Proof. Let the orthogonal complement of range(GDir)

ran(G_Dir)^⊥ := g ∈ L²(S^d−1) (g, ψ)_L²_(S^d−1₎ = 0 for all ψ ∈ ran(G_Dir) . Our aim is to show that ran(G_Dir)^⊥ = 0. Recall that

ran(G_Dir)^⊥ = ker(G^∗_Dir).

Given any f ∈ H^1/2(∂D), let v ∈ H_loc¹ (R^d\ D) be the unique solution to (3.4) and GDirf = v^∞ be its far-field pattern. Again, using the representation formula (3.13), we have

(g, G_Dirf )_L2(S^d−1)= (g, v^∞)_L2(S^d−1)

= Z

S^d−1

g(ˆx)

Z

∂D

(ν · A(z)∇Φ^∞_A,n(ˆx, z))v(z) ds(z)

 ds(ˆx)

− Z

S^d−1

g(ˆx)

Z

∂D

Φ^∞_A,n(ˆx, z)(ν · A(z)∇v(z)) ds(z)

 ds(ˆx)

= Z

∂D

ν · A(z)∇

 Z

S^d−1

Φ^∞_A,n(ˆx, z)g(ˆx) ds(ˆx)



v(z) ds(z)

− Z

∂D

 Z

S^d−1

Φ^∞_A,n(ˆx, z)g(ˆx) ds(ˆx)



ν · A∇v(z) ds(z)

= Z

∂D



(ν · A(z)∇u_A,n,g)f − u_A,n,g(ν · A(z)∇v)



ds. (3.14) Let u be the unique solution to (3.4) with boundary data u = u_A,n,g on ∂D, that is,







∇ · (A(x)∇u) + k²n(x)u = 0 in R^d\ D,

u = u_A,n,g on ∂D,

u satisfies Sommerfeld radiation condition at |x| → ∞.

(3.15)

Since both v and u satisfies ∇·(A(x)∇u)+k²n(x)u = 0 in R^d\D and satisfy Sommerfeld radiation condition, then

0 = Z

∂D



(ν · (A∇u))f − u_A,n,g(ν · (A∇v))



ds (3.16)

Subtracting (3.14) by (3.16) yields (g, G_Dirf )_L²_(S^d−1₎=

Z

∂D

(ν · (A∇(u_A,n,g − u)))f ds, that is,

G^∗_Dirg = ν · A∇(u_A,n,g− u).

If g ∈ ker(G^∗_Dir), then ν · A∇u_A,n,g = ν · A∇u. Combining this with (3.15), we can extend u in D by defining u = u_A,n,g, and then u satisfies

∇ · (A∇u) + k²n(x)u = 0 in R³.

Since u satisfies Sommerfeld radiation condition at |x| → ∞, by uniqueness result, we see that u ≡ 0, which gives u_A,n,g = 0 in D. By the unique continuation property, we hence have uA,n,g ≡ 0. Finally, Lemma 2.4 implies g ≡ 0, that is, ran(GDir)^⊥ = 0.  3.4. The properties of the far-field operator F_Dir. In view of the argument in [KG08, Theorem 1.8], we can prove an important property of F_Dir.

Lemma 3.13. The far-field operator F_Dir : L²(S^d−1) → L²(S^d−1) satisfies

F_Dir− F_Dir^∗ = i2kγ²F_Dir^∗ F_Dir, (3.17)

where F_Dir^∗ : L²(S^d−1) → L²(S^d−1) is the adjoint operator of F_Dir, and the constant γ is given in (3.3).

Proof. Given any g, h ∈ L²(S^d−1), choose v^inc = u_A,n,g and w^inc = u_A,n,h, which are modified Herglotz wave functions given in (3.5). Let v and w be the solutions of the scattering problem (3.1), with corresponding scattered fields v^sc = v − vⁱⁿ and w^sc= w − wⁱⁿ, which are solution to (3.2), and the far-field patterns of v^sc and w^sc are v^∞ and w^∞, respectively. By definition of F_Dir, we have

v^∞= FDirg and w^∞= FDirh.

Note that 0 =

Z

BR\D

(v∇ · (A∇w) − w∇ · (A∇v)) dx = Z

|x|=R



v(∂_rw) − w(∂_rv)



ds. (3.18) Clearly, for each R > 0,

Z

|x|=R



v^inc(∂rw^inc) − w^inc(∂rv^inc)



ds = 0. (3.19)

By the Sommerfeld radiation condition and the asymptotic expansion of the scattered field, we can compute

v^sc(x)(∂_rw^sc(x)) − w^sc(x)(∂_rv^sc(x))

= −ikv^sc(x)w^sc(x) − ikw^sc(x)v^sc(x) + o(|x|^−d−12 )

= −ik

 γ e^ik|x|

|x|^d−12 v^∞(ˆx) + O(|x|^−d+12 )



γe^−ik|x|

|x|^d−12 w^∞(ˆx) + O(|x|^−d+12 )



− ik



γe^−ik|x|

|x|^d−12 w^∞(ˆx) + O(|x|^−d+12 )



γ e^ik|x|

|x|^d−12 v^∞(ˆx) + O(|x|^−d+12 )



+ o(|x|^−d−12 )

= −i2kγ²

|x|^d−1v^∞(ˆx)w^∞(ˆx) + O(|x|^−d) uniformly for all ˆx = x/|x|. Hence we have

R→∞lim Z

|x|=R



v^sc(∂_rw^sc) − w^sc(∂_rv^sc)

 ds

= −i lim

R→∞

Z

|x|=R

2kγ²

|x|^d−1v^∞(ˆx)w^∞(ˆx) ds

= −i lim

R→∞

2kγ² R^d−1

Z

|x|=R

v^∞(ˆx)w^∞(ˆx) ds

= −i lim

R→∞

2kγ² R^d−1

Z

S^d−1

v^∞(ˆx)w^∞(ˆx)R^d−1ds

= −i2kγ² Z

S^d−1

v^∞(ˆx)w^∞(ˆx) ds

= −i2kγ²(F_Dirg, F_Dirh)_L2(S^d−1). (3.20)

Next, for each R > 1, by Fubini’s theorem and the representation formula of the far-field pattern, we can obtain that

Z

|x|=R



v^inc(∂rw^sc) − w^sc(∂rv^inc)

 ds(x)

= Z

S^d−1

g(ˆz)

 Z

|x|=R

Φ^∞_A,n(ˆz, x)(∂rw^sc) − w^sc(∂rΦ^∞_A,n(ˆz, x)) ds(x)

 ds(ˆz)

= − Z

S^d−1

g(ˆz)w^∞(ˆz) ds(ˆz)

= −(g, F_Dirh)_L²_(S^d−1₎, (3.21)

and similarly,

Z

|x|=R



v^sc(∂_rw^inc) − w^inc(∂_rv^sc)

 ds(x)

= (F_Dirg, h)_L2(S^d−1). (3.22)

Combining (3.18), (3.19), (3.20), (3.21), and (3.22) gives

−i2kγ²(F_Dirg, F_Dirh)_L2(S^d−1)− (g, F_Dirh)_L2(S^d−1)+ (F_Dirg, h)_L2(S^d−1) = 0,

which is our desired result. 

With the help of (3.17), we can investigate the range property of the far-field operator F_Dir.

Lemma 3.14. Let SDir : L²(S^d−1) → L²(S^d−1) be the scattering operator defined by S_Dir = I + i2kγ²F_Dir,

then S_Dir is unitary, i.e., S_Dir^∗ S_Dir = S_DirS_Dir^∗ = I. Moreover, F_Dir is normal.

Proof. It is straightforward to see that S_Dir^∗ S_Dir =



I − i2kγ²F_Dir^∗



I + i2kγ²F_Dir



= I + i2kγ²(F_Dir− F_Dir^∗ ) + 4k²γ⁴F_Dir^∗ F_Dir.

(3.23)

Combining this and (3.17), we obtain S_Dir^∗ S_Dir = I. Since F_Dir : L²(S^d−1) → L²(S^d−1) is a compact operator, then S_Dir is a compact perturbation of the identity, and so S_Dir is bijective. Hence, S_Dir⁻¹ = S_Dir^∗ , which shows that S_Dir is unitary. Note that

S_DirS_Dir^∗ =



I + i2kγ²F_Dir



I − i2kγ²F_Dir^∗



= I + i2kγ²(F_Dir− F_Dir^∗ ) + 4k²γ⁴F_DirF_Dir^∗ .

(3.24)

By comparing (3.23) and (3.24), we conclude that F_Dir is normal.  Lemma 3.15. If k² is not an eigenvalue of (3.10), then F_Dir : L²(S^d−1) → L²(S^d−1) is injective and ran(F_Dir) is dense in L²(S^d−1).

Proof. We first establish the injectivity. Given g ∈ ker(FDir). Recall that F_Dirg is the far-field pattern corresponding to the incident field

v^inc(x) = u_A,n,g(x).

Let v^scbe the associated scattered field. As above, applying Rellich’s lemma and unique continuation property, we obtain v^sc = 0 in R^d. This shows that v^inc satisfies (3.10).

Since k² is not an eigenvalue for (3.10), then v^inc = 0 in D. By the unique continuation property for the acoustic equation, we then have v^inc = uA,n,g ≡ 0 in R^d. The injectivity of g → uA,n,g in Lemma 2.4 implies g ≡ 0.

Next, to prove that ran(F_Dir) is dense in L²(S^d−1), it suffices to show that F_Dir^∗ : L²(S^d−1) → L²(S^d−1) is injective. For each g, h ∈ L²(S^d−1), the reciprocity relation for the far-field pattern gives

(g, F_Dirh)_L²_(S^d−1₎= Z

S^d−1

g(ˆx)

Z

S^d−1

v_Dir^∞(ˆx, ˆz)h(ˆz) ds(ˆz)

 ds(ˆx)

= Z

S^d−1

g(ˆx)

Z

S^d−1

v_Dir^∞(−ˆz, −ˆx)h(ˆz) ds(ˆz)

 ds(ˆx)

= Z

S^d−1

 Z

S^d−1

v_Dir^∞(−ˆz, −ˆx)g(ˆx) ds(ˆx)



h(ˆz) ds(ˆz)

=

 Z

S^d−1

v^∞_Dir(−ˆz, −ˆx)g(ˆx) ds(ˆx), h



L²(S^d−1)

, that is,

(F_Dir^∗ g)(ˆz) = Z

S^d−1

v_Dir^∞(−ˆz, −ˆx)g(ˆx) ds(ˆx).

Therefore, if g ∈ ker(F_Dir^∗ ), then by the injectivity of F , we conclude that g ≡ 0.  3.5. The (F^∗F )^1/4-method. Now we are going to derive our main reconstruction method. The result is based on the following theorem in [KG08].

Theorem 3.16. [KG08, Theorem 1.23] Let H be a Hilbert space and X be a reflexive Banach space. Assume that the compact operator F : H → H have a factorization of the form

F = B ˜AB^∗ with operators B : X → H and ˜A : X^∗ → X such that

(1) =hϕ, ˜Aϕi 6= 0 for all 0 6= ϕ ∈ ran(B^∗).

(2) ˜A = ˜A₀ + C for some compact operator C and some self-adjoint operator ˜A₀ which is coercive on ran(B^∗).

If F is injective and I + irF is unitary for some r > 0, then ran(B) = ran((F^∗F )^1/4).

Putting together Proposition3.5, Lemma3.8–3.11,3.14,3.15, and applying the above theorem, we can immediately obtain the following lemma.

Lemma 3.17. If k² is not an eigenvalue of (3.10), then ran(GDir) = ran((F_Dir^∗ FDir)^1/4).

Finally, we want to characterize ran((F_Dir^∗ F_Dir)^1/4) as in [KG08] for the purpose of numerical simulations. Since F_Dir : L²(S^d−1) → L²(S^d−1) is compact and normal, there exists a set of (complex) eigenvalues {λ^j_Dir}_j∈N with corresponding eigenfunctions {φ^j_Dir}_j∈N in L²(S^d−1). Furthermore, the set of eigenfunctions forms an orthonormal

basis of L²(S^d−1), see e.g. [Zim90, Corollary 3.2.9]. For g ∈ L²(S^d−1), we can write (F_Dir^∗ F_Dir)^1/4g =X

j∈N

q

|λ^j_Dir|(g, φ^j_Dir)_L2(S^d−1)φ^j_Dir, Φ^∞_A,n(·, z) =X

j∈N

(Φ^∞_A,n(·, z), φ^j_Dir)_L²_(S^d−1₎φ^j_Dir.

Therefore, if Φ^∞_A,n(·, z) ∈ ran((F_Dir^∗ F_Dir)^1/4), then there exists a g ∈ L²(S^d−1) such that X

j∈N

q

|λ^j_Dir|(g, φ^j_Dir)_L²_(S^d−1₎φ^j_Dir=X

j∈N

(Φ^∞_A,n(·, z), φ^j_Dir)_L²_(S^d−1₎φ^j_Dir. Equivalently,

q

|λ^j_Dir|(g, φ^j_Dir)_L²_(S^d−1₎ = (Φ^∞_A,n(·, z), φ^j_Dir)_L²_(S^d−1₎ for all j ∈ N, or

(g, φ^j_Dir)_L²_(S^d−1₎ = (Φ^∞_A,n(·, z), φ^j_Dir)_L²_(S^d−1₎ q

|λ^j_Dir|

for all j ∈ N.

Hence, Φ^∞_A,n(·, z) ∈ ran((F_Dir^∗ F_Dir)^1/4) if and only if X

j∈N

|(Φ^∞_A,n(·, z), φ^j_Dir)_L2(S^d−1)|²

|λ^j_Dir| < ∞.

Combining this with Lemma3.17and Proposition3.6, we get the following key theorem.

Theorem 3.18. If k² is not an eigenvalue of (3.10), then the following are equivalent:

(1) z ∈ D;

(2) Φ^∞_A,n(·, z) ∈ ran((F_Dir^∗ F_Dir)^1/4);

(3) WDir(z) :=

 X

j∈N

|(Φ^∞_A,n(·, z), φ^j_Dir)_L²_(S^d−1₎|²

|λ^j_Dir|

−1

> 0.

In other words, the characteristic function of the sound-soft obstacle D is given by χ_D(z) = sign(W_Dir(z)).

4. The determination of the obstacle with other boundary conditions We can extend the factorization developed above to the obstacle with other types of boundary conditions on ∂D, i.e., Neumann or impedance boundary conditions. In the case of Helmholtz equation (homogeneous medium) with the obstacle having these two types of boundary conditions, the factorization method has been discussed in detail in [KG08, Chapter 1,2]. For the acoustic equation considered here, once we use Φ^∞_A,n(ˆz, x) as the incident field, the factorization method for obstacles with Neumann or impedance boundary conditions can be established following exactly the arguments in [KG08, Chapter 1,2]. Since we have presented the detailed proof for the Dirichlet case in Section3, we will not repeat the arguments here. Instead, we state the main theorems of the reconstruction method without proofs.

We first discuss the Neumann case (sound-hard). The (F^∗F )^1/4-method can be easily modified to treat this case as in [KG08, Chapter 1]. Let F_Neu be the far-field operator corresponding to the sound-hard obstacles, which is defined similarly as F_Dir.

Theorem 4.1. If k² is not an eigenvalue of

(∇ · (A(x)∇u) + k²n(x)u = 0 in D,

ν · (A(x)∇u) = 0 on ∂D, (4.1)

then the following are equivalent:

(1) z ∈ D;

(2) Φ^∞_A,n(·, z) ∈ ran((F_Neu^∗ F_Neu)^1/4);

(3) W_Neu(z) :=

 X

j∈N

|(Φ^∞_A,n(·, z), φ^j_Neu)_L²_(S^d−1₎|²

|λ^j_Neu|

−1

> 0, where {λ^j_Neu, φ^j_Neu} is an eigen-system of the operator F_Neu.

In other words, the characteristic function of the sound-hard obstacle D is given by χ_D(z) = sign(W_Neu(z)).

Now we consider the obstacle with impedance condition. Let FImp be the corre- sponding far-field operator. In this case, F_Imp fails to be normal in general, unless λ is real-valued (i.e. Robin boundary condition). However, the problem can be overcomed by the F_]-method, which can be found in [KG08, Chapter 2]. The self-adjoint operator F_] is defined as follows:

F_]= |<F_Imp| + |=F_Imp|, where

<FImp := 1

2(FImp+ F_Imp^∗ ) and =FImp = 1

2i(FImp− F_Imp^∗ ).

In this case, the result reads:

Theorem 4.2. Let λ ∈ L^∞(∂D) be a complex-valued function with non-negative the imaginary part on ∂D. If k² is not an eigenvalue of

(∇ · (A(x)∇u) + k²n(x)u = 0 in D,

ν · (A(x)∇u) + λ(x)u = 0 on ∂D, (4.2)

then the following are equivalent:

(1) z ∈ D;

(2) Φ^∞_A,n(·, z) ∈ ran(F_]^1/2);

(3) WImp(z) :=

 X

j∈N

|(Φ^∞_A,n(·, z), φ^j_Imp)_L2(S^d−1)|

|λ^j_Imp|

−1

> 0, where {λ^j_Imp, φ^j_Imp} is an eigen-system of the operator F].

5. Numerical results

In this section, we present some numerical simulation results to show the efficiency of the method developed in Section 3. For simplicity, we treat the acoustic equation in two dimensions. The simulations are obtained using MATLAB 2020a (with PDE Toolbox). We consider four shapes of sound-soft obstacles whose parametric equations are listed in Table1.

Obstacle shape Parametrization (anti-clockwise oriented) circle x(t) = 0.5 cos t

0 ≤ t ≤ 2π y(t) = 0.5 sin t

peanut [YZZ13] x(t) = cos t ·√

cos²t + 0.25 sin²t

0 ≤ t ≤ 2π y(t) = sin t ·√

cos²t + 0.25 sin²t kite [YZZ13] x(t) = 0.5 cos t + 0.325 cos 2t − 0.325

0 ≤ t ≤ 2π y(t) = 0.75 sin t

heart x(t) = −0.5 sin³t

0 ≤ t ≤ 2π (Wolfram Mathworld) y(t) = ¹³₃₂cos t −₃₂⁵ cos 2t −₁₆¹ cos 3t − ₃₂¹ cos 4t

Table 1. Parametrization of the obstacles

To carry out the simulations, we first need to compute the incident field u^inc_A,n(x, ˆz) = Φ^∞_A,n(ˆz, x) = u^to_ref(x, −ˆz).

By the mixed reciprocity relation described above, it suffices to compute the total field for the acoustic equation (without obstacle) with plane incident fields. After generating the required incident field, we then simulate the scattered field and compute the far- field pattern for the acoustic equation with a sound-soft obstacle. The simulation is explained in the following algorithm.

Algorithm 1 Simulation of the scattered field u^sc_Dir and the far-field pattern u^∞_Dir

1: Choose an incident direction ˆz ∈ S¹;

2: Compute the scattered field u^sc_ref(x, −ˆz) for the acoustic equation (1.1) without ob- stacle by the incident field u^inc_ref(x, −ˆz) = e^{ikx·(−ˆz)};

3: Compute the total field u^to_ref(x, −ˆz) = u^sc_ref(x, −ˆz) + u^inc_ref(x, −ˆz);

4: Compute the scattered field u^sc_Dir(x, ˆz) for the acoustic equation (3.2) with incident field u^inc_A,n(x, ˆz) = u^to_ref(x, −ˆz);

5: Compute the far-field of u^∞_Dir and generate the far-field operator FDir.

We now discuss the algorithm in more detail. In our simulation, we take A ≡ I and n(x, y) =

( 1 + e⁻

1

1+x2+y2 for x²+ y² < 1

1 otherwise

which has jump discontinuities at x²+ y² = 1. The original problem was formulated in R², together with the Sommerfeld radiation condition at infinity. Here, we restrict the computational domain in the ball x² + y² < 4 and approximate the Sommerfeld radiation condition by the impedance condition:

(∆˜u^sc_ref + k²n˜u^sc_ref = k²(1 − n)u^inc_ref in x²+ y² < 4,

∂_ru˜^sc_ref − ik ˜u^sc_ref = 0 on x²+ y² = 4.

We solve this boundary valued problem FEM with mesh size ≤ 0.1. Then the approxi- mated total field ˜u^to_ref = ˜u^sc_ref + u^inc_ref and the needed incident field

˜

u^inc_I,n(x, ˆz) = ˜u^to_ref(x, −ˆz).