Detection of thin high contrast dielectrics from boundary measurements

We develop an efficient inversion method for thin high contrast scatterers when the contrast is on the order of the reciprocal of the thickness of the scatterer. We extend prior theory for the Helmholtz equation to arbitrary bounded domains and multiple scatterers in two and three dimensions to obtain a fast forward solver with complexity of one dimension lower. The lower-dimensional approximation is then paired with optimization to form the basis for parameter inversion. We show numerical results for the forward and inverse problems in two dimensions and describe extensions to Maxwell.


Introduction
It is well known in inverse scattering theory that if one has apriori information about the structure of the medium, aysymptotic approximation can be a tool to invert for relevant unknown parameters and may provide a natural regularization. Here we extend the results of one of the authors and her collaborators from [1] in this sense. There, Moskow, Santosa, and Zhang studied the Helmholtz equation in free space in the case of a thin, high contrast dielectric scatterer. In the limit as the height of the scatterer vanishes, simultaneously as the contrast goes to infinity, a limiting system of integral equations of one lower dimension was found to be satisfied. In the present contribution, we extend this theory to general bounded domains and multiple scatterers in two and three space dimensions and apply the lower-dimensional asymptotic approximation as a method for inversion. We find that the asymptotic results are similar in this case since the Green's functions, although different, have the same singularity.
Two of the authors, Ambrose and Moskow, previously worked with thin high contrast scatterers in the timeharmonic Maxwell equations (first analyzed in [2]), proving under some hypotheses that the thin high contrast problem has a somewhat different limiting system than for the Helmoltz case [3]. Ambrose, Gopalakrishnan, Moskow, and Rome subsequently used this limiting system as the basis of a method for inversion [4]. The inverse method in [4] was a direct method done under the assumption that the plane of the scatterer is known. In the present work, we also develop and implement (for Helmholtz) a method for inversion to locate the scatterer, but we do not make such an assumption and instead use an optimization algorithm.
In the case of low contrast, one approximates the scattered field by a known background or incident field which is then used in an integral formulation to yield the well known Born approximation, see for example [5]. A restatement of the Born approximation could be that a limiting problem, in which there is zero contrast, is used for inversion. We thus provide the high contrast analogue, loosely stated, by analyzing the case of infinite contrast, and using this for inversion. We find that unlike for the Born approximation, this high contrast analogue is nonlinear in the scatterer and captures interaction between multiple scatterers. This is also the case for other high contrast asymptotic limits such as in [6,7].
The paper is organized as follows. In section 2, we present the problem formulation for Helmholtz and describe the regime in which the analysis here will hold. Next, in section 3 we show the asymptotic results. In section 4 we describe the limiting solution for multiple scatterers and restate the main theorem for that case. In section 5 we describe the change of variables for an arbitrary slit and show numerical results for the forward approximations, while section 6 contains our inversion experiments. Finally in section 7 we remark on how this approach can potentially be extended to Maxwell's equations using the Lippmann-Schwinger-type formulation shown rigorously to hold in [8].

Problem formulation
Here we consider a scalar wave problem on a smooth bounded domain Ω in  d , for d=2,3, modeled by the Helmholtz equation in the presence of a thin high contrast scatterer. The Helmholtz equation is given by with Neumann boundary conditions where the coefficient ( ) n x h 2 models the squared index of refraction of a thin high contrast scatterer. We assume thatn 1 h 2 has support on a thin object S h ⊂Ω . For simplicity we assume that Ω contains the origin, and that =´- We separate out the variable in the transverse direction and write We assume that the scatterer is thin and high contrast in the sense that That is, the analysis here applies to the regime where the contrast of the thin scatterer is on the order of its thickness, as was studied in [1] for the free space case. We assume for simplicity that ( ) n x 0 is continuous, although piecewise continuous and bounded in ¥ ( ) L S 0 would suffice. Related to the problem (1) is the corresponding background problem Note that this solution corresponds to the field in the case when no scatterer exists in Ω, and is our analog to the incident wave in the free space setting. We assume also here that for given ( ) n x h 2 and fixed h, k 2 is not an eigenvalue for (1) nor for (3), that is, for both problems we have existence of a unique solution.
If the scatterer S h were to have bounded contrast, it is well known that for small thickness h the field u h would be close to the background field u i . However, in our high contrast regime, the field will instead be close to the solution of a crack problem. To see this, we use the background Green's function and integration by parts to first write a Lippmann-Schwinger form of (1), where the Green's function G satisfies Using the high contrast model (2), we also set ¢ = ¢ Solving for u 0 is a d−1 dimensional problem, and therefore yields a computationally efficient approximation to u h in all of Ω,

Estimates and convergence
We begin by formally deriving the first terms in an asymptotic expansion with respect to the parameter h, and then show convergence results. Taking z=x d /h which is as in [1] but with different notation, we define h h d and the scaled version of (6), We write the ansatz and expand and plug these into (10). The O(1) terms yield the equation (7) for u 0 . To obtain the O(h) term we also need to expand the term Let S=S 0 ×(−1/2, 1/2) be the scaled scatterer and define the operator 1 2 0 0 where we also abuse notation and use T to denote the same operator from ( ) L S 2 to L 2 (S). Indeed, T is bounded and compact in both cases, since the singularity in G is the same as that of the free space fundamental solution, and so even when restricted to the lower-dimensional S 0 the kernel is only weakly singular. A detailed proof of this for free space G can be found in [1]. We note that the range of T contains only functions which are constant in z, and we can simplify (13) to The following lemma gives us an explicit expression for v 1 .
Lemma 3.1. Assume k 2 is not an eigenvalue for the problems (7), (3) so that the solution u 0 exists and is unique and that n 0 is continuous on S 0 . If v h is defined by (11) and v 0 is defined by (12), then Proof. Let  be the free space fundamental solution for the Helmholtz equation. Note that when d=3 we have the formula and when d=2 we have instead for ( ) H 0 1 the Hankel function of the first kind. Then the difference  -G satisfies a homogeneous Helmholtz equation in Ω with Neumann data equal to the normal derivative of , which will be smooth for ¢ x bounded away from the boundary of Ω. For any smooth integrand f, we have by Taylor expansion and the fact that the linear term integrates to zero. Hence Since functions in the range of the operator T do not depend on the z variable, we see that if we write this will indeed solve (16) ifû 1 is the solution to the lower-dimensional integral equation given by . Assume also thatk 1 2 is not an eigenvalue for T given by (14), and let . Then there exists > C 0, independent of h, such that the following estimate holds: Furthermore, if u 0 is extended to all of Ω by the right hand side of (9) we have that  - Proof. The proof of the first result follows as in [1], but we give an abbreviated version for completeness. One sees that for T defined by (14),  x hz dx  Ch  sup  , 0 ,  , 0  ,  ,  , , 20 x z S S , 0 since  and G differ by a smooth function for . Using this and noting that the right-hand side of (19) is continuous, and we can bound the residual by where we have combined the last two terms of (19). Using the triangle inequality yields By our assumption that −1/k 2 is not an eigenvalue of T, the Fredholm theory implies that + -( ) I k T 2 1 is a bounded operator on C 0 (S). This, combined with moving the second term on the right-hand side of (22) over to the left-hand side yields the first result. For the second bound, first we abuse notation somewhat extend u 0 to be where inside the integral u 0 refers to thed 1 dimensional solution to (7). We invoke (6) to write the difference as . Assume also thatk 1 2 is not an eigenvalue for T given by (14).
, and let u 1 be defined by (17) whereû 1 is the unique solution to (18). Then we have that - Proof. We again describe briefly the proof which is as in [1]. One again can calculate the residual ò ò ò + - There are seven terms on the right-hand side of (25), and of these, the first, second, third, fourth, and seventh are O(h 2 ) in L 2 (S) because of the smoothness of u i , (20), and proposition 3.2. We therefore focus on the fifth and sixth terms, which we will now group together as one quantity; this is pointwise o(h) from lemma 3.1. If we divide by h, it is pointwise almost everywhere going to zero, and by (20) uniformly bounded. The Lebesgue dominated convergence theorem implies therefore that it is going to zero in L 2 (S). Hence the right hand side of (25) is o(h) in L 2 (S). Again by our assumption that −1/k 2 is not an eigenvalue of T and the Fredholm theory we have that + -( ) I k T 2 1 is a bounded operator on L 2 (S), and hence the result follows. ,

Multiple scatterers
For the Born approximation and other fixed contrast asymptotic regimes, in the first order approximations the scatterers are decoupled from one another, and so multiple scattering effects can not be captured without higher-order corrections. Furthermore, in many cases the scatterers need to be well separated from one another for the approximations to be valid. In the high contrast regime studied here, the limiting problem both captures multiple scattering and is valid for close and even crossing thin scatterers. We consider d=2 for simplicity in what follows, but similar results can be obtained when d=3. Assume that n h 2 takes the more general form x h x S h h and 0 h j j j j 0 j is also contained in Ω, where ν j is the unit normal to the line segment S j 0 . We assume that the line segments S j 0 are not overlapping for segments of positive measure in one dimension, but we allow them to possibly intersect at points, assuming that n 2 h is well-defined. We will let be the set of strips and the finite set of limiting line segments respectively. We assume also that n (j) 0 varies only in the direction tangential to S 0 j . Then as before we consider the solution to with n h defined by (26) and with Neumann boundary conditions Corresponding to this is the limiting problem to find Î ( ) By linear transformations one has the analogue of (20) on each ( ) S j 0 , and from the proof of (3.2) one has the following theorem: We remark that when the line segments S j 0 intersect at points, for finite h the strips S j h overlap for areas of O(h 2 ), which will not affect the limit of the Lippmann-Schwinger formulation. If the segments were to overlap for nontrivial portions of the line, the Lippmann-Schwinger formulation and its limiting problem (28) would double count these portions, which would correspond to a different coefficient n 2 h . For such a medium, the line segments should just be redefined so as not to overlap.

Numerical simulations for the forward problem
In this section we investigate the accuracy of the limiting solution u 0 given in general by (28). For the solutions of these limiting integral equations, we use a one-dimensional collocation method with piecewise linear basis functions and the FEniCS assemble function to compute the integrals which form the system matrix. These onedimensional solutions are inserted into the right hand side of (30) to obtain the approximate solution on all of Ω. For the brute force solution to (27) we use FEniCS to compute a piecewise linear FEM solution directly.
We first treat the case that the slit is at the origin, positioned horizontally, and of unit length. In this base case, the limiting one dimensional problem is to solve for u 0 (x 1 ) on the interval [−1/2, 1/2], where M is the number of nodes and f i are piecewise linear basis functions, and we set (31) to hold for x 1 at all node points. For a general configuration of the slit, we define the affine coordinate transformation from original coordinates y 1 and y 2 to yieldỹ 1 andỹ : As the parameter h approaches zero, the integral can likewise be reduced to the one-dimensional equation which we also solve by collocation. Once u 0 is found, to compute u(x) for any xäΩ we can use the formula

The high contrast inverse problem
The asymptotic approximation shown above provides a lower-dimensional and therefore computationally inexpensive forward solver for the field in the presence of a thin high contrast inhomogeneity. If one knows a priori that the medium contains a certain number (or even a maximum number) of inhomogeneities of this form, such a solver can be used in conjunction with optimization to recover them.
We assume here that  W = -´-Î ( ) ( ) 1, 1 1, 1 2 , n 0 ≡1, k=1. To start we assume we have only one slit with which is transformed by some unknown scaling, rotation, and translation. We do not assume we know h. For synthetic data, we will solve the full equation (27) with FEniCS FEM solutions as in the last section and read u on ∂Ωfor 4 piecewise constant Neumann sources. We then use the BFGS optimization routine in SciPy to image the slit, that is we optimize over the angle, length, and center point of the slit in the formula (34). We minimize where V(x; a, b) are the solutions computed on the boundaries using the dimensionally reduced (34) for the four sources, and f is our synthetic data. Note that we do not use nor do we assume we know h in the inversion procedure, but our plots of the resulting slits are thickened by h for illustration purposes. We do not recover h, either; to do that one would need a correction term, which will be the subject of future work. In figures 6, 7, 8, 9, and 10 we demonstrate results of the optimization procedure. In all cases we successfully find the location, angle and length of the slit.
In the final example we attempt to recover two slits, while still using only four sources. In figure 11 we initialize with two horizontal slits at the center and on top of each other. The algorithm separates them and finds their locations and angles successfully. We note that we expect the recovery of more complicated geometries and several slits will require more sources. This is the subject of future work. for Î x S 0 , and  the Helmholtz free space fundamental solution. In [4], a direct inversion method was proposed for this problem which made use of the next-order asymptotics. However, a more robust approach should result from inserting the two-dimensional approximation (38) into the formulation for Maxwell's equations of Lippmann-Schwinger type which was analyzed in [8]. Since the measurements are on the boundary (which is assumed to be bounded away from the scatterer), there will be no singularity in the integral even in limiting case. This will yield an approximation to the scattered field on the boundary which involves only twodimensional integrals an hence an efficient forward solver to do optimization as above. Indeed since for x in the