Topological optimality condition for the identification of the center of an inhomogeneity

The inverse scattering problem for inhomogeneous media is considered within the topology optimization framework. Varying the complex-valued refractive index we derive a zero-order necessary optimality condition in minimizing the L2 misfit cost functional of the far-field measurement. The topology asymptotic expansion of the optimality condition leads to an imaging operator, which is used to identify the center of the unknown inhomogeneity using few far-field measurements. Numerical tests show high precision and stability in the reconstruction using our optimality condition based imaging both in two and three dimensions.


Introduction
We consider the inverse scattering problem for the inhomogeneous media modeled by the Helmholtz equation in two and three dimensions. Such problems arise e.g. in non-destructive testing, medical imaging and many other areas of science and engineering. A popular class of methods for solving these problems, called qualitative methods, are able to quickly and accurately determine the shape and location of hidden objects and some information on the refractive index (at best), and require little a priori information about the objects. They are non-iterative in nature and do not require large scale wave simulations as oppose to iterative methods [15,20,26,27,35], however the latter recover everything (when possible) about the inhomogeneity. We refer the reader to [9][10][11]30] for a comprehensive account on the linear sampling methods and factorization method for the inhomogeneous media. Many other noniterative techniques are developed, including point source methods [39], enclosure method [24], MUSIC [4,38], and one wave methods [25,36,40] (our references are not exclusive since the literature on these methods is quite large).
A different point of view that has also led to non-iterative type techniques for solving inverse problem stems from the shape optimization. Our study falls in this category. Identification of small inhomogeneities in the context of shape optimization is developed in [8,18]. Based on topological derivatives [14,17,42], the optimization approach was adapted to inverse scattering problems in [2,5,6,13,19,41] (see also references therein). The first-order topological derivative specified for small inhomogeneities with real-valued refractive index can be found in [1,3], and the high-order topological expansions are derived in [7,22,31]. The asymptotic analysis in combination with perturbation theory is used in the framework of topological sensitivity for many inverse problems (see e.g [21,28,34,37] and references therein).
More specifically, in this study we consider an inhomogeneity of the arbitrary shape described by a bounded region D ⊂ R d , d = 2, 3 centered at x 0 and constant refractive index n. With respect to a parameter ε > 0, this domain is rescaled to fit in a ball of radius ε centered at x 0 ; this epsilon is referred to as the size of this scaled inhomogeneity. We setup the inverse problem as a minimization of the L 2 -least square misfit cost functional for a given measurement of the far-field pattern. The condition of feasible measurement expressed in a weak form (see (17)) provides us with the zero-order necessary optimality condition which is derivativefree. Then applying asymptotic arguments as the size ε 0 + to the optimality condition we derive an imaging function determined only by the incident field and far-field measurements. This enables us to accurately and stably reconstruct the center x 0 of the inhomogeneity with the minimum number of measurements needed being equal to the dimension d = 2, 3. This work extends to the case of inhomogeneous media the method proposed in [32] where similar type of imaging function was introduced to identify the center of an impenetrable impedance obstacle from boundary measurements. Our result could be used in combination with other perturbation techniques to augment the recovered information about small inhomogeneity (e.g. known the center of the inhomogeneity we could use the formulas in [12] to recover the refractive index). Our zero-order optimality condition could possibly motivate the construction of other imaging functions that could work for broader class of inhomogeneous media. We present here a variety of numerical example showing the accuracy of our reconstruction method as well as indicating what could happen in cases when our theory does not apply, e.g. in the case of multiple inhomogeneities or inhomogeneities with absorption.

Formulation of the scattering problem
Consider an infinite homogeneous background acoustic medium occupying all of R d with d = 2 or 3 and characterized by the constant wave velocity c 0 and the frequency ω > 0. Let k = ω/c 0 denote the background wave number in the background medium. Let D be an open, connected and bounded domain with Lipschitz boundary ∂D. We start by assuming that 0 ∈ D ⊂ B 1 (0), where B 1 (0) is the unit ball around origin 0. This assumption is to separate the far-field R d \ B 1 (0) from the near-field B 1 (0) \ D of the object D. Assuming that the minimal radius r = 1 among all bounding balls B r (0) ⊃ D, we denote by G D the set of such shapes D. After rescaling D ∈ G D by a size parameter ε > 0 this allows us to produce uniquely the geometric object located at a center which is inscribed in the ball B ε (x 0 ) of radius ε and center x 0 . Now in our problem D ε (x 0 ) denotes the support of a scattering inhomogeneity characterized by a real constant (we will comment later on the case when n can be a function or complex) index of refraction n = c/c 0 (with c denoting the wave velocity in D ε ). In general, in the presence of absorption n is complex such that Re(n) > 0 and Im(n) 0. We extended n to R d by setting n = 1 in R d \ D ε (x 0 ). The topology variables (D, ε, x 0 , n) ∈ G D × R + × R d × C + will be used for variation of inhomogeneities in our inverse scattering problem later on.
In the following formulation of the forward scattering problem D ε (x 0 ) and n are fixed. Considering an incident field u i that is a known solution of the unperturbed Helmholtz equation ∆u i + k 2 u i = 0 in R d called metaharmonic function [44] (typically such incident field is set to be a plane wave u i := e ıkx· z with incident direction z ∈ S and S denoting the unit circle if d = 2 or the unit sphere if d = 3). Then the forward acoustic scattering problem under consideration is to find the total field u ∈ H 1 where the scattered field u s satisfies the Sommerfeld radiation condition (1c) uniformly in x = x/|x| ∈ S. The latter condition implies the existence of a far-field pattern u ∞ ∈ L 2 (S) [15] such that The radiating fundamental solution Φ of the Helmholtz equation in R d is given by where H (1) 0 is the Hankel function of the first kind and order zero. The far-field pattern Φ ∞ ( x, y) of Φ(x, y), such that the asymptotic expansion holds, is given by Φ ∞ ( x, y) = γ d e −ık x·y for y ∈ R d , where the parameter γ d is such that It is well-known that the solution u of (1a)-(1c) satisfies the Lippmann-Schwinger equation (see e.g. [29, theorem 6.8]) where the volume potential defines the bounded linear operator (see [29, If |n − 1|k 2 T (D,ε,x0) < 1 implying a contraction mapping, then from (5) it yields the convergent Neumann series for the solution The far-field pattern u (D,ε,x0,n) ∞ of the scattered field corresponding to (1a)-(1c) admits the expression For the following consideration we introduce the Herglotz operator H : Note that Hg is an entire solution to the Helmholtz equation. For a given geometry D ε (x 0 ) we denote by H (D,ε,x0) : Then from (9) and (11) it follows the following expression for the far-field pattern u (D,ε,x0,n) ∞ of the scattered field corresponding to (1a)-(1c) which will be used later on in the proof of theorem 1.
The inverse scattering problem under consideration is to find the location x 0 of the inhomogeneity D ε (x 0 ) from a knowledge of measured far-field pattern corresponding to a few (to become precise later) incident waves without knowing its refractive index n. Next we derive a stable criteria for solving this inverse problem based on an optimality condition associated with the topology sensitivity functional.

Optimality condition for the topological sensitivity functional
Let the true (unknown) inhomogeneity be D ε (x 0 ) with refractive index n and let us denote the corresponding far-field measurement due to a given incident field by u ∞ ∈ L 2 (S). The latter is what we measure. For variation of the topology, a trial inhomogeneity D ε (x 0 ) with a trial refractive index n is put in the background medium. This constitutes a family of solutions u (D,ε,x0,n) of the forward problem (1a)-(1c) for the same incident field. We denote by u (D,ε,x0,n) ∞ the corresponding far field pattern. To identify the unknown inhomogeneity, we set the least square L 2 cost functional J : which should be minimized over the trial variables. This minimization problem reads: find We say that the measurement u ∞ is feasible with respect to the quadruple (D , ε , x 0 , n ), if the solution u (D ,ε ,x 0 ,n ) of the corresponding forward problem obeys exactly this far-field pattern In the feasible case (15), the trivial minimum in (14) is attained. Note that for uniqueness results on determining the support D ε (x 0 ) of the inhomogeneity (of polygonal shape or a ball) with one incident wave we refer to [16] and references therein. Now we formulate the zero-order necessary optimality condition for (14), which is derivative-free. Theorem 1. Assume that Im(n ) = 0 and the measurement u ∞ ∈ L 2 (S) is feasible. Then the following necessary optimality condition for (14) holds where the Herglotz operator H is defined in (10), is the solution of the forward scattering problem (1a)-(1c) for the true inhomogeneous medium, and the constant γ d is given in (4).
Proof. For the feasible measurement u ∞ , after multiplication with the complex conjugate u of a test function u and integrating over the unit ball boundary S, the identity (15) takes the following weak form for all test-functions u ∈ L 2 (S). Next we insert in (17) (11) and use the integral representation (9) for the far-field pattern u Using the identity (see [5, lemma 1]): and interchanging the order of integration we arrive at Since Im(n ) = 0, then D = D in the first equation in (18), which implies that the imaginary part of D is zero, and the second equation in (18) after complex conjugation follows (16), thus proving the assertion of the theorem. □ Note that the factor γ 3 is real-valued in 3d, and hence it can be omitted from formula (16) in three dimensions. (17), compared to its strong form (15), could potentially be used, by playing with test-functions, in order to seek other zero-order necessary optimality conditions in minimizing the least square L 2 cost misfit functional for given measurements.

Remark 1. The weak variational form of the optimality condition
Next we apply to the necessary optimality condition (16) the asymptotic argument as ε 0 + , i.e. the volume of the support of the inhomogeneity becomes arbitrarily small.

Asymptotic optimality condition
In order to arrive at an optimality condition that can be applied to solve the inverse problem we employ the topology asymptotic analysis using near-field expansions. To this end, in R d we assign to the center of the inhomogeneity x 0 a local d-dimensional coordinate system with the radial coordinate |x − x 0 | and d − 1 angular coordinates, i.e the polar coordinates in 2d, and the spherical coordinates in 3d. Furthermore x − x 0 := Since u i is an entire solution of the Helmholtz equation in an arbitrary ball B R (x 0 ) of radius R > 0 and center x 0 , it has the following expansion in the two-dimensional case (see [15, section 3.4]): with the Fourier coefficients K ∈ C and the Bessel functions J of order ∈ Z, whereas in the three-dimensional case (see [15, section 2.4]): where K m ∈ C, Y m are the spherical harmonics, and j are the spherical Bessel functions of order ∈ N 0 . From (19) we find that u i (x 0 ) = K 0 J 0 (0), whereas from (20) that u i (x 0 ) = K 0 0 j 0 (0)Y 0 0 , and in the both 2d and 3d cases, we infer the following near-field asymptotic representation of u i The Herglotz wave Hu ∞ defined in (10) is an entire solution of the Helmholtz equation, thus it satisfies Hence, for a small size parameter ε > 0, (21) and (22) in the ball of radius R = ε imply that The Born approximation stated in (7) and (8) now takes the form Plugging the asymptotic relations (23) and (24) into the necessary optimality condition (16) proves straightforwardly the following theorem.

Theorem 2.
Assume that Im(n ) = 0 and the measurement u ∞ ∈ L 2 (S) is feasible. Then the necessary optimality condition (16) has the asymptotic form as ε 0 + : Remark 2. The principal asymptotic term in (25) is determined by the incident field u i which is an entire solution to the Helmholtz equation and by the Herglotz wave function Hu ∞ which density is the far-field measurement corresponding to this incident field.
We give an interpretation of formula (25) from the point of view of a topological derivative following [5]. To this end, we first decompose the cost functional J in (13) as follows Interchanging the order of integration and using (9) implies that u (D,ε,x0,n) Hu ∞ dy.

F Cakoni and V A Kovtunenko Inverse Problems 34 (2018) 035009
Applying the asymptotic formula (23) and (24) we now get the two-parameter (ε, n)-asymptotic representation of the cost functional J where we have assumed that |n − 1|k 2 T (D,ε,x0) < 1 according to (8), and |D ε (x 0 )| = ε d |D| for the Hausdorff measure of the sets in R d . Hence from (26) it follows the following expression of the topological derivative of J : Theorem 3. For fixed n ∈ C + and ε 0 + we have that be the solution of the topology optim ization problem (14). This means that the optimal value of the cost functional J (D , ε , x 0 , n ) = 0 in (26). Varying the imaginary part of the optimal refractive index Im(n ), we observe that the second term in the expression of the topological derivative from (27): should be asymptotically zero according to theorem 2.
Based on theorem 2, we suggest to consider the real-valued imaging operator I : L 2 (S) → C(R d ) which, for given incident field u i is defined over the space of far-field patterns u ∞ as follows Due to (25), the continuous function Iu ∞ possesses the property that (Iu ∞ )(x 0 ) = 0 asymptotically at the center x 0 of the unknown inhomogeneity D ε (x 0 ). This implies that the zerolevel set L =0 of Iu ∞ contains (asymptotically) the inhomogeneity center x 0 : Consequently, this suggest the test center x 0 is looked for (asymptotically) as the intersection point of d zero-level sets: from (u ∞,1 , . . . , u ∞,d ) different far-field measurements in d dimensions.
An asymptotic model of similar type as (30) was justified numerically in [32] for the imaging of impedance obstacles from boundary measurements in the near-field. It was shown to produce high-accuracy and stable identification results. In the next section, we show that this is the case for the inhomogeneous media too, namely (30) produces accurate and stable location of the inhomogeneity center from far-field measurements.

Direct scattering problem
For computation of the solution of the forward scattering problem (1a)-(1c) we use its equivalent formulation in the form of Lippmann-Schwinger integral equation (5) with the weakly singular operator T (D,ε,x0) given in (6).
Without loss of generality, let Ω = (0, 1) d be the finite computational test domain, that is the unit square in 2d and the unite cube in 3d, containing a homogeneity D ε (x 0 ) with a refractive index n. For a fixed mesh size h ∈ R +, we denote the nodal points of the uniform grid (quadrilateral in 2d and polyhedral in 3d) in the closure Ω by Based on the fast solver from [43], we apply to (5) and (6) the following formula of numerical integration over x j ∈ N h : (31) The system matrix of (31) is sparse and non-symmetric. In [43, theorem 1], for fixed k = 1 and sufficiently small h, the unique solvability of the linear system (31) is established and in addition the O(h 2 | ln(h)|)-error estimate in the max-norm: is proved. We examine numerically the error defined in (32) with respect to the product kh following [23]. The error is presented here for a typical numerical test in 2d. In this example, the inhomogeneity D is circle-shaped of size ε = 2 −4 centered at x 0 = (0.25,0.75), with the complex-valued refractive index n = 2 + ı, which is illuminated in the direction of θ = π 3 by the plane wave u i (x) = e ık(x· z) with z = (cos θ, sin θ). The wave numbers are taken in the range of We observe that the discretization error in plots (a) and (b) of figure 1 is super-linear, and the pollution error (see the explanation in [23]) in plot (c) is moderate. From our numerical tests we confirm these features of the error Error(kh) for wave numbers k and mesh sizes h chosen such that kh 1.
In the far-field, one needs to discretize the unit circle in 2d or the unit sphere in 3d. Here M denotes the finite set of nodal points on the unit ball boundary x ∈ S in R d . Then the discretization of the far-field (see [43, section 3.5]) for x ∈ M is given by Finally, for far-field measurement u ∞ given at nodes x ∈ M , the Herglotz operator (10) can be discretized for x j ∈ N h as Proper weights w M ( x) in (34) should be determined depending on meshing of S. For example, on the unit circle in 2d we use equidistant nodes x l = (cos θ l , sin θ l ) with θ l = 2πl N θ for l = 1, . . . , N θ and fixed N θ 1, then all weights w M ( x l ) = 2π N θ . For the uniform latitude-longitude grid on the unit sphere in 3d such that x (l,m) = (cos θ l sin φ m , sin θ l sin φ m , cos φ m ) with θ l as above and φ m = π(m−1) N φ −1 for m = 1, . . . , N φ and fixed N φ 2, we use the cubature formula based on the trapezoidal rule with the positive weights

Optimality-based-imaging of the center
In our numerical tests, in order to compute the imaging function (28) for identifying the center of a given inhomogeneity with support D ε (x 0 ) and refractive index n , we use the numerical solution u (D ,ε ,x 0 ,n ) h given by (31), in order to synthesize the discrete far-field measurement u ∞ = (u (D ,ε ,x 0 ,n ) h ) ∞ according to formula (33) and the corresponding Herglotz wave H M u ∞ as per (34). In this setting from the discrete far-field measurement u ∞ ( x) given at nodes x ∈ M and the incident field u i (x j ) known at mesh points x j ∈ N h , depending on these two grids we calculate the discrete imaging function according to (28) for x j ∈ N h as follows In the following we show some examples as proof of concept showing the identification of the inhomogeneity center based on the imaging function from (35).
In figure 2 we present the result for the case when the inhomogeneity D is L-shaped of size ε = 2 −3 centered at x 0 = (0.125, 0.375), with the real-valued refractive index n = 0.5. The computation is done over the unit square Ω = (0, 1) 2 with the mesh sizes h = h = 2 −6 , and N θ = 10 equidistant points on the unit circle (relatively few measurement points).
Two far-field patterns u ∞,1 and u ∞,2 correspond to the case when the inhomogeneity is illuminated by a plane wave u i (x) = e ık(x· z) with z = (cos θ, sin θ) in the direction of θ 1 = π 6 and θ 2 = π 4 , respectively. Plots is unique and close (to become quantitative measure later) to the test center x 0 . Moreover, x (h,M) 0 = x 0 exactly, if the size is sufficiently small such that ε 2 −6 . Here zero-level sets L =0 (I (h,M) u ∞ ) are utilized numerically using the narrow strip technique, see e.g. [33]. More specifically, in a narrow strip of nodes where I (h,M) u ∞ changes its sign, we interpolate the discrete imaging function by multivariate polynomials which are bilinear in 2d (respectively, trilinear in 3d). This technique determines zero-level sets on a local grid which maybe finer than N h in the whole computational domain. Typically, we set it to N h in the narrow strip. It is important to stress that the wave number is chosen small: k < k 1 below the threshold (k 1 ≈ 2.25 in this example) when the reconstruction starts to break down. We observe that starting from approximately k = 1.35, there appear multiple lines of zero sets (see the line in the upper right corner in figure 2) and they grow when increasing k.
In figure 3 we show the same experiment for k = 2.25. In this case the level curves are not necessary lines any longer and there appear multiple intersection points. A possible remedy of this situation is to use an additional measurement, which confirms the true intersection point as illustrated in figure 3 for three measurements with θ 1 = π 6 , θ 2 = π 4 , and θ 3 = π. Since the asymptotic expansions in section 4 are carried out in the near-field of the inhomogeneity center x 0 , the imaging function cannot distinguish between multiple inhomogeneities put in the test domain. In figure 4 we give an illustrative example of triple inhomogeneities: the ball-shaped of size ε • = 2 −6 centered at x 0,• = (0.25, 0.75) with the refractive index n • = 2, the triangle-shaped of size ε = 2 −5 centered at x 0, = (2/3, 1/3) with n = 10 , and the L-shaped of size ε L = 2 −4 centered at x 0,L = (0.125, 0.375) with n L = 0.5. Here the wave number k = 1 and three plane waves propagate along the directions of θ ∈ { π 6 , π 3 , π}. From figure 4 we conclude that the zero-level sets intersection is closer to those inhomogeneity mostly contrasting with the background refractive index n = 1 (here the triangle-shaped inhomogeneity). While the remaining inhomogeneities can be considered as a perturbation of the homogeneous background medium. If the refractive indexes all are equal n • = n = n L = 2 here, then the intersection is close to a mass-center of the geometric set of inhomogeneities in the test domain, compare figures 4 and 5.
We carried out a series of numerical tests for inhomogeneities D ε (x 0 ) of various shapes, namely L-shape, ball, triangle, point (in the last case ε = 0), varying the location in the test domain Ω = (0, 1) 2 , illumination directions, the wave number, and the mesh sizes in the range   of h, h ∈ [2 −6 , . . . , 2 −3 ]. Although our theory covers only real refractive index n we test our imaging criteria (36) for refractive index with small imaginary part. In the following we report our numerical findings. Figure 6 presents how the error |x (h,M) 0 − x 0 | of center-identification depends on various parameters explained below. This is a 2d study.
(i) The intersection point of zero-level sets of the imaging function, calculated for a number of different measurements not less than the dimension, is uniquely determined for the wave numbers k 1 + δ with δ 0 depending on the problem parameters. (ii) The identification error depends super-linearly on the mesh size h of the grid used to synthesize the far-field measurement u ∞ . It is not sensitive to the mesh size h 2 −3 used for the imaging function, and to the grid M if chosen uniformly on the unit circle with N θ 5 nodes. (iii) If the test center x 0 coincides with a node of the grid M h , this leads to the better identification result than x 0 ∈ M h . For comparison see the two curves of the error depicted in the plot (a) in figure 6 in terms of the size ε . (iv) The identification error |x (h,M) 0 − x 0 | drops super-linearly with respect to decreasing the size ε of the unknown inhomogeneity as shown in figure 6(a). For comparison, the errorcurves for the ball of size ε = 2 −5 and for the equilateral triangle of size ε = 2 −6 are presented through the other plots (b)-(d).
Remark 3. We list some conclusions obtained from our two-dimensional numerical study: (1.) Numerical examples demonstrate high-precision of the identification of the center of the inhomogeneity as well as stability with respect to discretization and noise. (2.) In spite of the asymptotic arguments of theorem 2, we emphasize that the imaging technique (35) and (36) is applicable numerically also for medium with small absorption and even large inhomogeneities (up to 25% of the size of the test domain) depending on the specific parameters of the problem.