A globally convergent numerical method for a 3D coefficient inverse problem with a single measurement of multi-frequency data

The goal of this paper is to reconstruct spatially distributed dielectric constants from complex-valued scattered wave field by solving a 3D coefficient inverse problem for the Helmholtz equation at multi-frequencies. The data are generated by only a single direction of the incident plane wave. To solve this inverse problem, a globally convergent algorithm is analytically developed. We prove that this algorithm provides a good approximation for the exact coefficient without any \textit{a priori} knowledge of any point in a small neighborhood of that coefficient. This is the main advantage of our method, compared with classical approaches using optimization schemes. Numerical results are presented for both computationally simulated data and experimental data. Potential applications of this problem are in detection and identification of explosive-like targets.


Introduction
We are interested in this paper in the inverse problem of recovering the spatially distributed dielectric constant of the Helmholtz equation from a single boundary measurement of its solution. This inverse problem is also called a coefficient inverse problem with a single measurement data (CIP). CIPs are both ill-posed and highly nonlinear. Therefore, an important and challenging question to address in a numerical treatment of a CIP is: How to obtain at least one point in a sufficiently small neighborhood of the exact coefficient without any advanced knowledge of this neighborhood? As soon as this point is obtained, one can refine it using one of well established techniques, such as, e.g., a Newton-like method or a gradient-like method, see, e.g., [4]. It is well known that a rigorous guarantee of convergence of such a technique to the exact solution can be obtained only if the starting point of iterations belongs to a sufficiently small neighborhood of that solution.
We call a numerical method for a CIP approximately globally convergent (globally convergent in short, or GCM) if: (1) a theorem is proved, which claims that, under a reasonable mathematical assumption, this method delivers at least one point in a sufficiently small neighborhood of the exact solution without any advanced knowledge of this neighborhood, (2) the error of the approximation of the true solution depends only on the level of noise in the boundary data and on some discretization errors. Such an assumption is necessary due to the well known tremendous challenge of the goal to develop those numerical methods for CIPs, which positively address the question posed in the previous paragraph. This is especially true for CIPs with single measurement data, which are known to be one of the most challenging ones. In our particular case that assumption amounts dropping a small term of an asymptotic expansion, and this is used only on the zero iteration of our method, see Section 4.3. Due to our target application, the CIP considered here is a CIP with single measurement data. We refer to [5,6] for detailed discussions of the notion of the approximate global convergence.
Our CIP has many applications in sonar imaging, geographical exploration, medical imaging, nearfield optical microscopy, nano-optics, see, e.g., [13]. However, the target application of this paper is in imaging of dielectric constants of antipersonnel mines and improvised explosive devices (IEDs). We model the so-called "stand off" detection, which means that we use only backscatter data. Currently the radar community relies only on the intensity of radar images [37] to see the geometrical information such as, e.g., sizes and shapes of targets. Thus, we hope that the additional information about values of dielectric constants of targets, which our method delivers, might lead in the future to a better differentiation between explosives and clutter.
CIPs of wave propagation are a part of a bigger subfield, Inverse Scattering Problems (ISPs). ISPs attract a significant attention of the scientific community. There is a large body of the literature on imaging methods for reconstructing geometrical information about targets, such as their shapes, sizes and locations. We refer to, e.g., [2,3,9,12,19,21,30,31] and references therein for well-known imaging techniques in inverse scattering such as level set methods, sampling methods, expansion methods, and shape optimization methods.
The most popular approach to solve CIPs is nonlinear optimization methods minimizing some cost functionals, see, e.g., [4,10,15,17,40]. However, in general, such functionals are non-convex and might have multiple local minima and ravines, see, e.g., Figure 1 in [36] for a convincing example of the latter. Therefore, a good initial guess for the true solution is required. Such a first guess might be found by solving the linearized problem using the Born approximation. As a results, the relative target/background contrast is assumed to be small, while it is actually high in many real world applications (see, e.g. the table of dielectric constants on the web site of Honeywell with shortened link by Google https://goo.gl/kAxtzB). This is one of the reasons for us to establish a GCM to solve this inverse problem without knowing a good initial approximation of the true solution.
In [5,6] an approximately globally convergent method was established for a CIP with single measurement data for a wave-like PDE using Laplace transform. That method was verified on experimental data, see [23,27,38,39], and references cited therein. The GCM of [5,6] was extended in [11] to a CIP for a different wave equation.
The goal of this paper is to develop a new GCM for a CIP for the Helmholtz equation in R 3 with a single measurement of multi-frequency data. In other words, we now work in the frequency domain. Formally, the Helmholtz equation can be obtained via the Fourier transform with respect to time of the solution of the hyperbolic equation considered in [5,6]. One of the reasons for us not to do so is because the kernel of the Laplace transform, used in the previous GCM [5,6], decays exponentially fast and therefore some significant information might be lost. The main difficulty in developing this new GCM is to work with complex-valued functions where the maximum principle, which plays an important role in the previous GCM, is no longer applicable. We first develop the theory for the method. Next, we conduct a numerical study for computationally simulated data. At the end, we present one numerical example for experimental data. In fact, this is one of results of the work [32] of our group.
In [32] the method of the current publication was successfully tested on experimental data collected in the frequency domain. In addition, in [26] this method was tested on experimental data collected in the time domain: after a certain preprocessing procedure, these data were "moved" in the frequency domain via the Fourier transform. In both works [26,32], results of the reconstructions of locations and dielectric constants of inclusions were quite accurate. The 1D version of the method under consideration along with its performance for experimental data can be found in [24].
As to the CIPs with multiple measurement, i.e., the Dirichlet-to-Neumann map data, we mention the recent works [1,20,33] and references cited therein, where reconstruction procedures are developed, which do not require a priori knowledge of a small neighborhood of the exact coefficient.
The paper is organized as follows. In Section 2, we state the forward and inverse problems. In Section 3, we discuss some facts about the Lippmann-Schwinger equation. We need these facts for the analysis of our numerical method. In Section 4, we establish an integral differential equation for the CIP and derive an initial approximation for the target coefficient. In Section 5, we present our globally convergent numerical method and formulate the corresponding algorithm. In Section 6, we prove the main theorem of this paper: the global convergence theorem. In addition, we discuss this theorem in Section 6 and present the main discrepancies between the theory and computations. In Section 7, we present our numerical results for both computationally simulated and experimental data. Section 8 is devoted to a short summary.
be the ball of the radius R > 0 with the center at {x = 0} . Let Ω 1 ⋐ Ω ⋐ B(R) be convex bounded domains and ∂Ω ∈ C 2+α . Here and below C m+α are Hölder spaces where m ≥ 0 is an integer and α ∈ (0, 1) . Let c(x) be a function defined in R 3 . Throughout the paper we assume that The Riemannian metric generated by the function c(x) is This metric generates geodesic lines. Let be a > R an arbitrary number. Consider the plane P a = x ∈ R 3 : x 3 = −a . Note that since c (x) = 1 outside of the domain Ω 1 , then geodesic lines are straight lines outside of Ω 1 . In particular, since our incident plane wave propagates along the positive direction of the x 3 −axis (see this section below), we consider for x ∈ {x 3 < −R} only those parts of geodesic lines which are parallel to the x 3 −axis. Then these lines intersect the plane P a orthogonally. Throughout the paper we use the following assumption about the regularity of c(x) and geodesic lines generated by the function c(x).
. Furthermore, for any point x ∈ R 3 there exists unique geodesic line L(x) connecting x and the plane P a and such that L (x) intersects the plane P a orthogonally.
The length of the geodesic line L (x) is A sufficient condition for the validity of Assumption was found in [35], Remark 2.1.
1. The conditions in assumption 2.1 are technical ones. We use them only to derive the asymptotic expansion of the function u (x, k) with respect to k in Section 4.1. This expansion was proved in [25], where these conditions were used quite essentially. However, they are not imposed in computations, see also Section 6.2.
2. In the case of propagation of electromagnetic waves, as in the experimental data of [26,32], c (x) is the spatially distributed dielectric constant. We assume its independence on the wavenumber k since we end up working on a small k−interval. In this paper, we assume that the total field satisfies the Helmholtz equation rather than the full Maxwell's system. This is reasonable when the dielectric constant varies so slowly that it is almost constant over the distances of the order of the wavelength, see [7,Chapter 13] for more details. The good accuracy of our results for experimental data of [26,32] as well as in Section 7.8 of this paper points towards the validity of this assumption.

The Coefficient Inverse Problem
Below k > 0 is the wavenumber. Consider the function u 0 = exp (ikx 3 ), which represents the incident plane wave propagating along the positive direction of the x 3 −axis. Let u(x, k) be the solution of the following problem: (2.5) The functions u, u 0 and u sc are called the total field, incident field and scattered field respectively. It is well-known (see [13,Theorem 8.7]) that problem (2.5) has unique solution u(x, k) ∈ C 2 (R 3 ). In fact, u(x, k) ∈ C 2+α (R 3 ) for any α ∈ (0, 1), see [16,Theorem 6.17]. In this paper, we consider the following problem: Coefficient Inverse Problem (CIP). Let k and k be two positive numbers and k < k. Assume that the function g(x, k), is known where u(x, k) is the solution of the problem (2.5). Determine the function c(x) for x ∈ Ω.
In the current paper, we develop the theory of our numerical method for this CIP. We also numerically test our method for a CIP with partial data. Since we have only a single direction of the propagation of the incident plane wave, our CIP is an inverse problem with the data from a single measurement event. In this paper, we focus only on the numerical method for our CIP. As to the question of the uniqueness of the solution, it is worth mentioning that currently uniqueness for the single measurement case can be proved only if the homogeneous equation in (2.5) would be replaced with a non-homogeneous one, in which the right hand side would be nonzero everywhere in Ω. The proof of this fact was first introduced in [8] by using the Carleman estimate. There were many works of many authors since then exploring the method of [8]. We refer the reader to review papers [22,43]. In the case when the right hand side of equation in (2.5) has some zeros, an attempt to address the uniqueness question faces fundamental challenges which are yet unclear how to tackle. Hence, we assume below uniqueness for our CIP.
It is well known in the theory of Ill-Posed Problems that one should assume the existence of the exact solution, which corresponds to the ideal, noiseless data [4,5,40]. Hence, throughout the paper, c * (x) is that exact coefficient and we use the superscript * to indicate functions related to c * (x). Naturally, we assume that assumption 2.1 is valid for the function c * (x).

Using the Lippmann-Schwinger equation
In this section we use the Lippmann-Schwinger equation to derive some important facts, which we need both for our algorithm and for the convergence analysis. We use here results of Chapter 8 of the book [13]. Denote In this section the function β ∈ C α R 3 and it also satisfies conditions (2.1), (2.2). The Lippmann-Schwinger equation for the function u (x) := u(x, k) is If the function u (x) satisfies equation (3.1) for x ∈ Ω, then we can extend it for x ∈ R 3 \ Ω via substitution these points x in the right hand side of (3.1). Hence, to solve (3.1), it is sufficient to find the function u (x) only for points x ∈ Ω. Consider the linear operator K β defined as It follows from Theorem 8.1 of [13] that Here and below B = B(β, k, Ω 1 , Ω) > 0 denotes different constants which depend only on listed parameters. Therefore, the operator K β maps C α (Ω) to C α (Ω) as a compact operator. Hence, the Fredholm theory is applicable to equation (  Lemma 3.2. There exists unique solution u ∈ C 2+α (R 3 ) of the problem (2.5). Consequently (Lemma 3.1), there exists unique solution u ∈ C 2+α (R 3 ) of the problem (3.1) and these two solutions coincide. Furthermore, by the Fredholm theory the following estimates hold Proof. We need to prove only estimate (3.5). By (3.3) and (3.4) On the other hand, by (3.1) Thus, estimate (3.5) follows from (3.6) and (3.7). Lemma 3.3 follows from Lemmata 3.1 and 3.2 as well as from results of Chapter 9 of the book [41].
For all x ∈ Ω, k > 0 the solution u (x, k) of the problem (2.5) is infinitely many times differentiable with respect to k. Furthermore, ∂ n k u ∈ C 2+α (Ω) and Let the function χ ∈ C 2 (R 3 ) be such that (3.8) The existence of such functions χ is well known from the Real Analysis course. While Lemmata 3.1 and 3.2 are formulated for the case when the function β (x) is real valued, we now consider the case of a complex valued analog of this function. We will need this for our algorithm. Consider a complex valued function ρ(x) ∈ C α (Ω). Let ρ (x) = χ (x) ρ (x). Then We have the lemma.

Integral differential equation formulation
In this section, we eliminate the target coefficient c(x) from the governing equation (2.5). This process leads to an integral differential equation. Solving this equation is equivalent to solving our CIP. Note that this approach is different from locally convergent methods with iterative optimization processes which typically rely on least-squares formulation. Furthermore, under some mathematical assumptions, we rigorously derive a good initial approximation for c(x). The analysis relies on the high frequency asymptotic behavior of the total wave in the next section.

The asymptotics with respect to k → ∞
Let τ (x) be the function defined in (2.4). It was proven in [25, Theorem 1 and its consequence (4.26)] that for x ∈ B(R) the asymptotic expansion of the function u(x, k) with respect to k is where the function A(x) > 0 is defined in B(R) and the number B 1 = B 1 (B(R), c) > 0 depends only on listed parameters. By (4.1) and (4.2), given a coefficient c(x), one can choose a number We note that by (4.1), (4.2) we need only an estimate from the below for the number κ 0 rather than its exact value. The development of the theory of this paper would become much more complicated if we would work with the numbers κ 0 , k, k, which would depend on the function c. Hence, from now on, we assume that these numbers are the same for all functions c which we consider here. This assumption is supported from the Physics standpoint by our study of experimental data [26,32]. Thus, below numbers k, k, κ 0 , m 0 are the same for all functions c (x) considered here and also

The tail function and some auxiliary functions
We obtain from (4.3), (4.4) and a straightforward calculation that the vector ∇u( where C is a constant. The smoothness of the function u(x, k) and (4.5) imply that V (x) ∈ C 2+α (B(R)). Choosing C = 1, we have found a function V ∈ C 2+α (B(R)) such that , is determined uniquely up to the addition of a constant 2πni, n ∈ Z. We can choose n = 0. This choice determines the function V (x) uniquely. We name it the tail function.
2. The choice of n does not effect the results of our reconstruction method since we only use the for all x ∈ B(R).
We next eliminate the function c(x) from equation (4.9). To this end, define the function q(x, k), Differentiating equation (4.9) with respect to k, we obtain Plugging the function v(x, k) of (4.7) into (4.11), we obtain Lemma 4.2. Keeping in mind our numerical method, it is convenient to consider in this lemma that x ∈ Ω rather than the above x ∈ B (R) .
Lemma 4.2 (integral differential equation). The function q (x, k) satisfies the following nonlinear integral differential equation Now, solving our CIP is equivalent to solving problem (4.11)-(4.13). To do so, we iteratively approximate both functions q(x, k) and V (x), using the predictor-corrector scheme. Here approximations for V (x) are used as predictors and approximations for q(x, k) are used as correctors. We first start from an approximation V 0 (x) for the tail function V (x).

The initial approximation of the tail function
In this section we introduce a reasonable mathematical assumption mentioned in Section 1. The idea is to find an initial approximation for V using the asymptotic behavior of u(x, k) as k → ∞. Let c = c * and let condition (4.3) be fulfilled. We set in (4.1) u (x, k) = u * (x, k) and then ignore the term O (1/k) for k ≥ k. It follows from (4.2) that we can do this uniformly for all x ∈ B(R). Hence, we obtain the following approximate formula (4.14) Define the function v * (x, k) for k ≥ k as in (4.7). Then, the result in Lemma 4.1 is still valid for k ≥ k.
1. Formulas (4.14)-(4.18) form the first part of our reasonable mathematical assumption mentioned in Section 1. The second part is that we assume that numbers k 0 , k, k are the same for all functions c we consider here, see (4.3), (4.4) and the paragraph above them. We point out that formulas (4.14)-(4.18) are used only on the first iteration of our method. However, we do not use them on follow up iterations.

It follows from estimate (4.24) that the first approximation
belongs to a small neighborhood of c * (x), provided that the error in the data is sufficiently small. Therefore, the global convergence property is achieved just at the start of our iterative process.
3. We will prove in Theorem 6.1 that the global convergence property is still valid after certain further iterations. The latter property is important because it is well-known that the more data are used, up to a certain point, the better reconstruction is obtained, also see section 6.2.

The algorithm
We first consider the discretization of problem (4.12)-(4.13) with respect to the wavenumber k. Then, we prove the well-posedness of the resulting problem. In the last part of this section, we propose the globally convergent numerical algorithm.

The discretization of the interval of wavenumbers
To find q and V from problem (4.12)-(4.13) using an iterative process, we consider a discretization of this problem with respect to the wavenumber k. We divide the interval [k, k] into N subintervals with the uniform step size as follows We approximate the function q (x, k) as a piecewise constant function with respect to k, Equation (4.12) becomes We use here V n−1 instead of V since we will update the tail function in our algorithm. By Lemma 3.3 So, assuming that the discretization step size h is sufficiently small, we replace in (5.4) ∇V n−1 ∇q n with ∇V n−1 ∇q n−1 . Using again the fact that h is sufficiently small, we also drop in (5.4) the nonlinear term −2hk n (∇q n ) 2 as well as the terms −2∆(hq n ), −2∇(hq n ). However, we do not drop terms ∇(hQ n−1 ), ∆(hQ n−1 ) since Q n−1 is a sum of other terms, see (5.3). These approximations enable us to simplify the analysis as well as the computations of our method. Using (5.4), we obtain . Then (4.13) and (5.2) lead to the following Dirichlet boundary condition for the function q n (x) q n (x) = ψ n (x), x ∈ ∂Ω. In this section we study the existence and uniqueness of the solution q n ∈ C 2+α (Ω) of the Dirichlet boundary value problem (5.5)-(5.6). In the case of real-valued functions, existence and uniqueness of the solution follow from [16, Theorem 6.14]. The proof is based on the maximum principle and Schauder Theorem [16,28]. However, complex-valued functions cause some additional difficulties, since the maximum principle does not work in this case. Still, we can get our desired results for the case of complex-valued functions using the assumption that the norm h∇Q n−1 C 1+α is sufficiently small. Keeping in mind the convergence analysis in the next section, it is convenient to consider here problem (5.5)-(5.6) in a more general form Lemma 5.1. Assume that in (5.7) all functions are complex valued ones and also that p ∈ C 1+α (Ω), f ∈ C α (Ω), µ ∈ C 2+α (∂Ω). Consider a number γ ∈ (0, 1) . Then there exists a constant C 1 = C 1 (Ω) > 0 and a sufficiently small number σ = σ(C 1 , γ) ∈ (0, 1) such that if C 1 σ ≤ γ and ∇p C α (Ω) ≤ σ, then problem (5.7) has a unique solution w ∈ C 2+α (Ω). Moreover, Proof. Below in this paper C 1 = C 1 (Ω) > 0 denotes different positive constants depending only on the domain Ω. Let the complex-valued function v ∈ C 2+α (Ω). Consider the following Dirichlet boundary value problem with respect to the function U : The Schauder Theorem implies that there exists unique solution U ∈ C 2+α (Ω) of the problem (5.9) and Hence, for each fixed pair f ∈ C α (Ω), µ ∈ C 2+α (∂Ω), we can define a map which sends the function v ∈ C 2+α (Ω) in the solution U ∈ C 2+α (Ω) of the problem (5.9), say U = S f,µ (v) . Hence, S f,µ : Since the operator S f,µ is affine and C 1 σ < 1, (5.10) implies that S f,µ is a contraction mapping. Let the function w = S f,µ (w) be its unique fixed point. Then the function w ∈ C 2+α (Ω) is the unique solution of the problem (5.7). In addition, by (5.10) which implies (5.8).
Suppose that q s C 2+α (Ω) ≤ Y where Y > 0 is a constant. Let γ ∈ (0, 1) and C 1 = C 1 (Ω) be the numbers of Lemma 5.1. Denote a = k − k. Suppose that the number a is so small that Then there exists a unique solution q n,i ∈ C 2+α (Ω) of the problem (5.5)-(5.6). Moreover, Now, we are ready to establish the globally convergent algorithm. iii. Find u n,i (x, k) by solving the Lippmann-Schwinger equation (4.9) in which the function ρ(y) is replaced with χ(y)(c n,i − 1)(y). iv. Update ∇V n,i (x) = ∇u n,i (x, k)/u n,i (x, k). (c) Set c n = c n,m , ∇V n = ∇V n,m , q n = q n,m and Q n = Q n−1 + q n .

The algorithm
3. Let N ∈ {1, . . . , N } be the optimal number for the stopping criterion. Set the function c N as the computed solution of CIP.
Remark 5.1. The stopping criterion in step 3 and the choice of m in step 2b in the algorithm above are addressed computationally. We refer to Section 7.5 for a corresponding discussion. We also specify the computation of (5.12) in Section 7.3.

Global convergence
In this section we formulate and prove the global convergence theorem, which is the main theorem of this paper. Next, we present a discussion of this result.

The global convergence theorem
Let δ > 0 be the level of the error in the boundary data ψ n (x) in (5.6) in the following sense Following (4.10), let q * (x, k) = ∂ k u * (x, k) /u * (x, k) . All approximations for the function q * (x, k) and associated functions with the accuracy O (h) as h → 0, which are used in this section below, can be justified by Lemma 3.3. For n = 1, . . . , N , denote q * n (x) = q * (x, k n ) . For k ∈ [k n , k n−1 ) we obtain Since ∆v * n = div (∇v * n ) , by (4.9) where the function F * n (x) is clarified in this section below. While equation (5.4) is precise, we have obtained equation (5.5) using some approximations whose error is O (h) as h → 0. This justifies the presence of the term G * n (x) in the following analog of the Dirichlet boundary value problem (5.5)-(5.6): k n ∆q * n − 2k n ∇hQ * n−1 ∇q * n = −2k n ∇V * ∇q * n−1 + 2∆(−hQ * n−1 + V * ) + 2(∇(−hQ * n−1 + V * )) 2 + G * n (6.4) where F * n (x) and G * n (x) are error functions. Let M > 0 be a generic constant which we will specify later. We assume that q * n C 2+α (Ω) ≤ M, n = 1, . . . , N. (6.6) Now, we prove the main theorem of the paper. For brevity and also to emphasize the main idea of the proof, we formulate and prove Theorem 6.1 only for the case when inner iterations are absent, i.e. for the case m = 1. The case m ≥ 2 is a little bit more technical and is, therefore, more space consuming, while the idea is still the same. Thus, any function f nj in above formulae should be f n below in this section. Before stating the theorem, we recall that V 0 (x) is the first approximation of the tail as in Section 4.2 and that c n , n ≥ 1, is computed as in Section 5.3. where N ∈ [1, N ] is the optimal number for the stopping criterion as in Algorithm 5.1.
Remark 6.1. We refer to Section 6.2 for a discussion of Theorem 6.1.

Discussion of Theorem 6.1
It seems that estimate (6.7) indicates that the function c 1 should provide the best approximation for the exact coefficient c * . However, our computational experience tells us that more iterations should be performed. In fact, it is well-known that the use of more frequencies (i.e. more data) provides more stable reconstruction results. In addition, the convergence estimate (6.7) ensures that functions c n are located in a sufficiently small neighborhood of the exact coefficient c * , as long as h + δ is sufficiently small and the number of iterations n does not exceed N . At the same time, these estimates do not guarantee that the distances between functions c n and the function c * decrease when the iteration number n increases from n = 1 to n = N . Due to this phenomenon, the stopping rule for our algorithm should be chosen computationally. A similar phenomenon for Landweber iterations for a linear ill-posed problem is observed in Lemma 6.2 on page 156 of the book [15]. Note that a similar statement can be found on page 157 of [15].
As it was mentioned in the previous paragraph, we have obtained functions c n , n = 1, . . . , N , in a small neighborhood of the true coefficient c * . This is achieved without any advanced knowledge of a small neighborhood of c * . By (6.8) the approximation accuracy depends only on the level of the error δ in the boundary data and the discretization step size h with respect to the wavenumber k. Hence, Theorem 6.1 implies the global convergence of our algorithm. Furthermore, the global convergence is established for one of the most challenging types of CIPs: a CIP with the data resulting from a single measurement event. Two latter features form the key advantage of Theorem 6.1. On the other hand, the global convergence property is achieved within the framework of a reasonable mathematical approximation indicated in Section 4.3. Hence, to be more precise, this is the approximate global convergence property, as defined in Introduction and in [5,6].
It can be seen both from our algorithm and from the proof of Theorem 6.1 that the approximations (4.14)-(4.18) are used only to obtain the first approximation V 0 for the tail function and are not used on follow up iterations with n = 1, . . . , N . Recall that the number of iterations N can be considered sometimes as a regularization parameter in the theory of ill-posed problems [5,15,40].

Main discrepancies between the theory and numerical implementation
In our opinion, it is hard to anticipate that some discrepancies between the theory and computations would not occur for such a challenging problem as the one considered in this paper. So, as it often happens with other numerical methods, some conditions are relaxed in computations, compared with ones in the convergence theorem. Indeed, it is well known that the theory is usually more pessimistic than the computational results are. This is because the theory cannot grasp all features of computations. The main discrepancies are: 1. The Assumption 2.1 is used only to obtain the asymptotic formula (4.1), which was derived in [25]. However, since our target applications are in imaging antipersonnel mines and IEDs, the function c(x), in our numerical testing in Section 7 as well as in [26,32], is discontinuous on the inclusion/background interface.
2. Instead of computing the initial approximation of the tail as in Section 4.3, we find its derivatives directly by solving problem (7.7)-(7.8). This enables us to avoid the error arising from numerical differentiation. We recall that our algorithm requires the derivatives of the tail function rather than this function itself. In addition, the right hand side of formula (5.12) might be nonpositive. Hence, we amend this formula in computations via using formula (7.6) and a smoothing procedure (Section 7.3).
3. Our theory works only for the case when the data are given on the entire boundary of the domain Ω. However, due to our target applications in detecting and identifying explosives, we are also interested in the case of backscatter data. Thus, we complement the backscatter data as in (7.16).
4. We use the so-called "data propagation procedure" in Section 7.2. This is a heuristic procedure, which, nevertheless, is very effective and it is widely used in the optics community, see, e.g., [34,Chapter 2]. In particular, this procedure helps us to reduce the search domain as in (7.5). We also refer to [14] for a time-reversal technique that was exploited for reducing the computational domain of a coefficient inverse problem.

Numerical study
We describe in this section some details of the numerical implementation of the above globally convergent algorithm. Also, we present some reconstruction results for both computationally simulated and experimental data. It is convenient to denote in this section points x = (x 1 , x 2 , x 3 ) as x = (x, y, z) .

Dimensionless variables
We want to demonstrate in our computations that our algorithm works with ranges of parameters, which are realistic in our desired application in imaging of antipersonnel mines and IEDs. Hence, we start from variables, which have dimensions and make them dimensionless then. We are guided by our experimental data in [32].
Let 1 cm = 1 centimeter. Consider the dimensionless variable x = x/ (10 cm) . We do not change notations for brevity. For example, the dimensionless 0.6 means 6 cm. The experimental data of [32] are stable only in a small neighborhood of 3.1 GHz. Motivated by this fact, we will work on a corresponding interval of wavenumbers [6.2, 6.7] in our numerical study. Following the measurement arrangement of [32], we assume in our numerical study that the backscatter data are measured on a 100cm × 100cm rectangle, which we call "measurement plane". We assume that this plane is 76 cm away from the front face of the target to be imaged. Hence, in dimensionless variables we have measurements on a 10 × 10 rectangle.

Data propagation
The data propagation process is applied to the function u sc (x, k) on the measurement plane and it aims to approximate the data for u sc (x, k) on a plane, which is closer to the target. In other words, we "move" the data closer to the target. In doing so, we assume that we know in advance some rough estimates of locations of targets. However, we do not assume that we know exact locations. Such an advanced knowledge of those estimates goes along well with the main principles of the regularization theory [4,5,15,40].
The data propagation process firstly reduces the computational domain for our algorithm and secondly makes the scattered field data look more focused. The latter feature also provides us a reasonable estimate for the xy−location of the target, which is important for the numerical implementation. The data propagation is done using the so called angular spectrum representation, which is a well known technique in optics, see, e.g., [34,Chapter 2]. We also refer to [26,32,39], where this process has been exploited for the preprocessing of the experimental data. We display in Figure 1a and Figure 1b the absolute value of the noisy scattered field u sc (x, k) before and after propagation, respectively. The scatterer in this case consists of two different inclusions as in Figure 4. From the propagated data in Figure 1b we can see two clear peaks at the locations of the inclusions. Comparison of Figure 1a with Figure 1b demonstrates the effectiveness of the data propagation procedure.

Some computational details
In our computations, [k, k] = [6.2, 6.7] and N = 9. This means that we have 10 wavenumbers in this interval and from now on we refer to them as k n for n = 0, . . . , 9. All computationally simulated data considered below are perturbed by 15% of random artificial noise. More precisely, for each k n and for each grid point x on the measurement plane, the total field data were perturbed as g noise (x, k n ) = g(x, k n ) + 0.15 g(x, k n ) L 2 σ(x, k n ) σ(x, k n ) L 2 , (7.1)  where σ(x, k n ) = σ 1 (x, k n ) + iσ 2 (x, k n ), and σ 1 (x, k n ) and σ 2 (x, k n ) are random numbers, uniformly distributed on (−1, 1). The L 2 norm is taken over the measurement plane. Now assume that f (x, k) is the data obtained after the propagation. More precisely, let f (x, k) be the sum of the propagated scattered field and the incident field. In our algorithm we need the derivative of f (x, k) with respect k. We compute it via the finite difference as Even though the differentiation of the noisy data is an ill-posed problem, this approximation turned out to work well in our numerical implementation. This is partly because the propagated data ( Figure 1b) have less noise than the given data ( Figure 1a). Each of our targets to be imaged has its front face at the plane {z = 0}. The backscatter data u sc (x, k) are measured at the plane {z = −7.6} on the above mentioned rectangle We have observed that the data propagation process also provides a very good assistance in determining the locations of the targets of interest in the xy−plane. More precisely, the peaks that can be seen in magnitude of the propagated data (Figure 1b), are located exactly at the location of those targets in the xy−plane. Therefore, we can consider our search domain in the xy−plane as where J is the number of peaks of |f (x, k)| on the set Γ and O j is a neighborhood of the point where |f (x, k)| attains its j−th peak. These neighborhoods should not intersect with each other. The truncation value 0.7 was chosen by trial-and-error tests. When running our algorithm, we do not assume the knowledge of the fact that the front face of any target is located on the plane {z = 0} . Still, as we have stated above, we assume the knowledge of some rough estimates of linear sizes of targets of interest. In the z−direction we search for the target in the range (−0.75, 1). This search range is motivated by our desired applications for detection of explosives. Indeed, z ∈ (−0.75, 1) means that the linear size of the target in the z−direction does not exceed 17.5 cm while linear size of antipersonnel mines and IEDs are typically between 5 cm and 15 cm. Note that similar ranges were used in [26,32,39,38]. Thus, our search domain for the target is Ω T × (−0.75, 1). In addition, the right hand side of the formula (5.12) for the function c n,i might be complex valued in real computations. Hence, at each iteration we apply the following truncation to the coefficient c n,i (x) found by (5.12) c n,i (x) := max (|c n,i (x)|, 1) , x ∈ Ω T × (−0.75, 1), 1, elsewhere. (7.6) Next, the function c n,i (x) is smoothed by smooth3 in MATLAB with the Gaussian option. As a result, a new function c n,i (x) is obtained. We do not change notation here for brevity. Then, the so obtained function c n,i (x) is used for solving the Lippmann-Schwinger equation in the algorithm of Section 5.3. We solve the Lippmann-Schwinger integral equations (3.10) using a spectral method developed in [29,42]. This method relies on the periodization technique introduced by Vainikko, which enables the use of the fast Fourier transform in the approximation scheme. For solving the integral equation in our numerical implementation, the domain Ω defined in (7.4) is discretized with the uniform mesh size 0.067 in three dimensions. The code for this numerical solver was done in MATLAB. The boundary value problem (5.5)-(5.6) is solved by the finite element method. We use tetrahedral mesh with a local adaptive refinement for the domain Ω. We coded the numerical solver using FreeFem++ [18], which is a standard software designed with a focus on solving partial differential equations using finite element methods. We refer to www.freefem.org for more information about FreeFem++. All of the figures that are displayed below in this paper were done by visualization tools in MATLAB.

The first tail function V 0 (x)
It is clear from the algorithm (at n = 1) that we use only the gradient ∇V 0 , instead of the first tail function V 0 . Since we set ∆V 0 = 0 in Section 4.3, we find the derivatives of V 0 by solving the following problem instead of the problem (4.21) , on ∂Ω, j = 1, 2, 3. (7.8) Here, boundary condition (7.8) is from (4.5). Solving (7.7)-(7.8) also helps us avoid the error associated with the numerical differentiation of the first tail function V 0 (x) For the case where only the backscatter data are given on Γ, we simply complete the data on the other parts of the boundary by the incident plane wave e ikx 3 , see (7.16). In doing so, we approximately assume that u (x, k) = e ikx 3 for points x located in the domain Ω near ∂Ω \ Γ. Therefore, this completion leads to the boundary condition for ∂ x j V 0 on ∂Ω \ Γ as ∂ x j V 0 = 0, on ∂Ω \ Γ, j = 1, 2, (7.9) It is seen from condition (7.8) that we need the boundary data for functions ∂ x j u(x, k), x ∈ Γ, j = 1, 2, 3. The derivative with respect to x 3 for x ∈ Γ is calculated as whereg(x, k) is the propagated backscatter field. The function ∂ x 3g (x, k) was found by propagating the backscatter field to the two nearby planes P p = {x 3 = −0.75} and P pε = {x 3 = −0.75+ε}, subtracting the results from each other and dividing by ε, where we took ε = 0.1. Derivatives ∂ x j u(x, k), j = 1, 2, for x ∈ Γ were calculated using FreeFem++. The propagated data looked very smooth. Again, we have not observed any instabilities.

The stopping criteria and the choice of the final result
We address in this section the rules for stopping the iterations and the choice of the final result for our algorithm. These rules are guided by the convergence estimate (6.7) of Theorem 6.1 and also by the trial-and-error testing. We recall that estimate (6.7) only claims that functions c n,i are located in a sufficiently small neighborhood of the true solution if the number of iterations is not large. However, this theorem does not claim that these functions tend to the exact solution, also see Theorem 2.9.4 in [5] and Theorem 5.1 in [6] for similar results. Therefore, our stopping criteria are derived computationally. We note that the same stopping rules have been used in [32] for the experimental data. For the convenience of the reader we present these rules, which are taken from [32].
We have the stopping criterion for the inner iterations and another one for the outer iterations. Denote by e n,i the relative error between the two computed coefficients corresponding by two consecutive inner iterations of the n−th outer iteration e n,i = c n,i − c n,i−1 L 2 (Ω) c n,i−1 L 2 (Ω) , i = 2, 3, . . . (7.11) We consider the n−th and (n + 1)−th outer iterations which contains I 1 and I 2 inner iterations respectively. The sequence of relative errors associated to these two outer iterations is defined by e n,2 , . . . , e n,I 1 , e n+1,1 , e n+1,2 , . . . , e n+1,I 2 , (7.12) where e n+1,1 = c n+1,1 − c n,I 1 L 2 (Ω) c n,I 1 L 2 (Ω) .
The inner iterations with respect to i of the n−th outer iteration in the above algorithm are stopped when either e n,2 < 10 −6 or i = 3. We have observed in our numerical experiments that the reconstruction results are essentially the same when we use either 3 or 5 for the maximal number of inner iterations.
Concerning the outer iterations with respect to n, it can be seen from the stopping rule for the inner iterations that each outer iteration consists of at least two and at most three inner iterations. Equivalently, the error sequence (7.12) has at least three and at most five elements. We stop the outer iterations if there are two consecutive outer iterations for which their error sequence (7.12) has three consecutive elements less than or equal to 5 × 10 −4 . Again, we have no rigorous justification for these stopping rules; they rely on the content of the convergence theorem Theorem 6.1 and on trial-and-error testing. We point out, however, that we do trial-and-error only for one target, which we call "reference target". Then we use the same rules for all other targets.
We choose the final result for c comp (x) by taking the average of its approximations c n,i (x) corresponding to the relative errors in (7.12) that meet the stopping criterion for outer iterations. The computed dielectric constant is determined as the maximal value of the computed c comp (x), c max = max c comp (x). (7.13) We have observed in our numerical studies that we need no more than five outer iterations to obtain the final result.

Numerical examples with complete data
The target we consider for this example is a cube of the size 0.6 defined by We define the coefficient c(x) as c(x) := 5, x ∈ D, 1, elsewhere. (7.14) We assume that the backscatter data u sc (x, k) are given on the rectangle (7.2), and on this rectangle the function g(x, k) in (7.1) is replaced with the function u sc (x, k). We also assume that the function u(x, k) := g(x, k) in (7.1) is given on the part ∂Ω \ Γ of the domain Ω, where Γ and Ω are defined in (7.3) and (7.4) respectively. First, we use the data propagation and obtain the function u sc (x, k) on the rectangle Γ this way. Next, we consider Ω as the computational domain with the measurement data given on its entire boundary ∂Ω. Figure 2 displays the reconstruction result using our algorithm. It can be seen that both the maximal value and the location of the coefficient c(x) are well reconstructed. More precisely, for the computed coefficient c comp (x) we have c max = 4.9, (7.15) where c max is defined in (7.13). Hence, it follows from (7.14) and (7.15) that the relative error between the computed maximal value and the exact maximal value of the coefficient c (x) is only 2%. We visualize the shape and location of the exact target and the reconstructed one using isosurface in MATLAB. The isovalue was chosen as 50% of the maximal value of c comp (x), and this isovalue is applied to all the other examples in this paper when using isosurface visualization.

Numerical examples with backscatter data
We examine in this section the performance of our algorithm for the case of backscatter data. This case is of our main interest, which is motivated by our desired application in detecting and identifying of mines and IEDs, where one typically uses only the backscatter measurement in the stand off detection. More precisely, by backscatter data we mean that we are only given the function u sc (x, k) on the rectangle (7.2). Using the data propagation, we can consider the computational domain Ω in (7.4), where we have the backscatter data on the rectangle Γ defined in (7.3).
Since our method theoretically needs the data on the entire boundary of Ω, we complete the missing data on the other parts of Ω by the corresponding solution of the forward problem in the homogeneous medium, where c(x) = 1. In other words, assuming that g(x, k) = e ikz + u sc (x, k) is the data given on Γ, we extend it on the entire boundary ∂Ω as g(x, k) := g(x, k), x ∈ Γ, e ikz , x ∈ ∂Ω \ Γ.  We remark that data completion methods are widely used for inverse problems with incomplete data. The data completion (7.16) is a heuristic technique relying on the successful experiences of [26,32,39,38] when working with globally convergent numerical methods for experimental data.
As the first example for the case of backscatter data, we present in Figure 3 the reconstruction results of the coefficient c(x) defined in (7.14) for this case. It can be easily seen again that the location is well reconstructed. For the maximal value we have c max = 4.65, which implies that the relative error between the exact value and the computed value is 7%. Thus, the error has increased from 2% for the case of complete data to 7% for the case of backscatter data only.
In the second example we consider the target which consists of two similar but separate cubes. The length of the side of each cube is 0.4. The coefficient c(x) is defined as c(x) := 5, x ∈ D 1 ∪ D 2 , 1, elsewhere, (7.17)

Numerical examples with experimental data
In the last numerical example we address the performance of our algorithm for experimental multifrequency data. The data were collected by a microwave scattering facility at the University of North Carolina at Charlotte. These measured data are the backscatter signals corresponding to a single incident plane wave e ikz . The setting for the measurement is similar to that of the above considered backscatter data case.
It is worth mentioning that the raw data are contaminated by a significant amount of noise, see [32] for a detailed discussion. Furthermore, standard denoising techniques are inapplicable in this case since the raw data have a richer content than the computationally simulated data. Therefore, a heuristic data preprocessing procedure was applied to the raw data in [32]. The preprocessed data were used then as the input for our globally convergent algorithm. The aim of the data preprocessing of [32] is to make the resulting data look somewhat similar to the computationally simulated data. We refer to the paper [32] for all the details about the experimental setup as well as the data preprocessing.
The target here is a tennis ball of the radius 0.31 (see Figure 5(a)). The directly measured average dielectric constant of this target at 3 GHz (or equivalently k = 6.28) is 3.80, with 13% measurement error given by the standard deviation. This dielectric constant was independently measured by the physicists Professor Michael A. Fiddy and Mr. Steven Kitchin from the Department of Physics and Optical Science at the University of North Carolina Charlotte. Our reconstruction result displayed in   Figure 5 indicates that our method accurately reconstructs the location of the target. The maximal value for the computed coefficient is c max = 4.00, which is 5.3% error, compared with the direct measurement.
In [32] the above algorithm was tested for five sets of experimental data. Table 1 presents reconstruction results of the value c max for all five cases; "std.dev" means standard deviation of the error in the a posteriori measured dielectric constants at the frequency of 3 GHz. One can observe that the relative error of the reconstruction is between 2% and 10.3% in cases 1-4, and it is unclear in the case number 5.

Summary
We have analytically developed a new numerical method for solving a Coefficient Inverse Problem with single measurement data for the 3D Helmholtz equation. The approximate global convergence of this method is proved. Our targets of interest mimic antipersonnel mines and improvised explosive devices. Numerical results for both computationally simulated and experimental data demonstrate a   good accuracy in reconstructed locations and dielectric constants of targets. It is worth mentioning that, in the case of experimental data, this good accuracy is achieved regardless on the above mentioned high level of noise in the data.