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

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–11, 30] for a comprehensive account on the linear sampling methods and factorization method for the inhomogeneous media. Many other non-iterative 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\subset \mathbb{R}^d$, d  =  2, 3 centered at x₀ and constant refractive index n. With respect to a parameter $\varepsilon>0$, this domain is rescaled to fit in a ball of radius ε centered at x₀; this epsilon is referred to as the size of this scaled inhomogeneity. We setup the inverse problem as a minimization of the L²-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 derivative-free. Then applying asymptotic arguments as the size $\varepsilon\searrow0^+$ 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₀ 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.

Consider an infinite homogeneous background acoustic medium occupying all of $\mathbb{R}^d$ with d  =  2 or 3 and characterized by the constant wave velocity c₀ and the frequency $\omega>0$. Let $k=\omega/c_0$ denote the background wave number in the background medium. Let D be an open, connected and bounded domain with Lipschitz boundary $\partial D$. We start by assuming that $0\in D\subset B_1(0)$, where B₁(0) is the unit ball around origin 0. This assumption is to separate the far-field $\mathbb{R}^d\setminus B_1(0)$ from the near-field $B_1(0)\setminus D$ of the object D. Assuming that the minimal radius r  =  1 among all bounding balls $B_r(0)\supset D$, we denote by $\mathfrak{G}_D$ the set of such shapes D. After rescaling $D\in \mathfrak{G}_D$ by a size parameter $\varepsilon>0$ this allows us to produce uniquely the geometric object located at a center $x_0\in\mathbb{R}^d$

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle {D_\varepsilon(x_0) =\bigl\{x\in\mathbb{R}^d: {\frac{x-x_0}{\varepsilon}}\in D\bigr\}\subset B_\varepsilon(x_0)} \nonumber \end{align*}$$

which is inscribed in the ball $B_\varepsilon(x_0)$ of radius ε and center x₀. Now in our problem $D_\varepsilon(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₀ (with c denoting the wave velocity in $D_\varepsilon$). In general, in the presence of absorption n is complex such that ${\rm Re}(n)>0$ and ${\rm Im}(n)\geqslant 0$. We extended n to $\mathbb{R}^d$ by setting n  =  1 in $\mathbb{R}^d\setminus \overline{D_\varepsilon(x_0)}$. The topology variables $(D, \varepsilon, x_0, n) \in\mathfrak{G}_D\times \mathbb{R}_+\times \mathbb{R}^d \times\mathbb{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_\varepsilon(x_0)$ and n are fixed. Considering an incident field $u^{\rm i}$ that is a known solution of the unperturbed Helmholtz equation $\Delta u^{\rm i} +k^2 u^{\rm i} =0$ in $\mathbb{R}^d$ called metaharmonic function [44] (typically such incident field is set to be a plane wave $u^{\rm i} :={\rm e}^{\imath kx\cdot \widehat{z}}$ with incident direction $\widehat{z}\in \mathbb{S}$ and $\mathbb{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\in H^1_{{\rm loc}}(\mathbb{R}^d)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \qquad \Delta u +k^2n u =0 \qquad \qquad \qquad \qquad {\rm in }~~\mathbb{R}^d, \label{HelmTot} \nonumber \end{align} \tag{ 1a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u =u^{\rm i} +u^{\rm s}, \nonumber \end{align} \tag{ 1b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial u^{\rm s}}{\partial \vert x\vert } -\imath k u^{\rm s} = \mathcal{O} (\vert x\vert ^{-(d+1)/2}) \qquad \,{\rm for }~~~\vert x\vert \to\infty, \qquad \label{Somm} \nonumber \end{align} \tag{ 1c }$$

where the scattered field $u^{\rm s}$ satisfies the Sommerfeld radiation condition (1c) uniformly in $\widehat{x}=x/\vert x\vert \in \mathbb{S}$. The latter condition implies the existence of a far-field pattern $u^\infty\in L^2(\mathbb{S})$ [15] such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Ocal}{\mathcal{O}} \displaystyle \label{ffu} u(x) -u^{\rm i}(x)=\frac{{\rm e}^{\imath k\vert x\vert }}{\vert x\vert ^{(d-1)/2}}u_\infty(\widehat{x}) +\Ocal(\vert x\vert ^{-(d+1)/2}) \quad {\rm for }\quad \vert x\vert \to\infty. \nonumber \end{align} \tag{ 2 }$$

The radiating fundamental solution Φ of the Helmholtz equation in $\newcommand{\Rbb}{\mathbb{R}} \Rbb^d$ is given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Phi(x,y)={\frac{\imath}{4}} H_0^{(1)}(k\vert x-y\vert){\rm \ (\,for~ }d=2),\qquad \Phi(x,y)=\frac{{\rm e}^{\imath k\vert x-y\vert }}{4\pi\vert x-y\vert }{\rm \ (\,for}~ d=3), \nonumber \end{align*}$$

where $H_0^{(1)}$ is the Hankel function of the first kind and order zero. The far-field pattern $\Phi_\infty(\widehat{x}, y)$ of $\Phi(x, y)$, such that the asymptotic expansion

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Ocal}{\mathcal{O}} \displaystyle \label{ffPhi} \Phi(x,y)=\frac{{\rm e}^{\imath k\vert x\vert }}{\vert x\vert ^{(d-1)/2}} \,\Phi_\infty(\widehat{x}, y) +\Ocal(\vert x\vert ^{-(d+1)/2}) \quad {\rm for }\quad \vert x\vert \to\infty \nonumber \end{align} \tag{ 3 }$$

holds, is given by $\Phi_\infty(\widehat{x}, y) =\gamma_d{\rm e}^{-\imath k \widehat x\cdot y}$ for $y\in\mathbb{R}^d$, where the parameter $\gamma_d$ is such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{gammad} \gamma_2={\frac{1+\imath}{4\sqrt{\pi k}}} \qquad {\rm and} \qquad \gamma_3=1/(4\pi). \nonumber \end{align} \tag{ 4 }$$

It is well-known that the solution u of (1a)–(1c) satisfies the Lippmann–Schwinger equation (see e.g. [29, theorem 6.8])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{6} u =u^{\rm i} +(n-1) k^2 T_{(D, \varepsilon, x_0)} u, \nonumber \end{align} \tag{ 5 }$$

where the volume potential

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{7} (T_{(D, \varepsilon, x_0)} u) (x) :=\int_{D_\varepsilon(x_0)} \Phi(x,y) u(y) \,{\rm d}y \nonumber \end{align} \tag{ 6 }$$

defines the bounded linear operator (see [29, section 6.2]) $T_{(D, \varepsilon, x_0)}: L^\infty(D_\varepsilon(x_0)) \mapsto L^\infty(D_\varepsilon(x_0))$ with the operator norm

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{8} \Vert T_{(D, \varepsilon, x_0)}\Vert =\max_{x\in D_\varepsilon(x_0)} \Bigl\vert \int_{D_\varepsilon(x_0)} \Phi(x,y) \,{\rm d}y \Bigr\vert =\mathcal{O}(\varepsilon^2 \vert \ln\varepsilon\vert ^{3-d}). \nonumber \end{align} \tag{ 7 }$$

If $\vert n-1\vert k^2 \Vert T_{(D, \varepsilon, x_0)}\Vert <1$ implying a contraction mapping, then from (5) it yields the convergent Neumann series for the solution

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{9} u^{(D, \varepsilon, x_0,n)} :=u =\sum_{\ell=0}^\infty \bigl( (n-1) k^2 T_{(D, \varepsilon, x_0)} \bigr)^\ell u^{\rm i}. \nonumber \end{align} \tag{ 8 }$$

The far-field pattern $u^{(D, \varepsilon, x_0, n)}_\infty$ of the scattered field corresponding to (1a)–(1c) admits the expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{11} u^{(D, \varepsilon, x_0,n)}_\infty (\widehat{x}) :=\gamma_d (n-1) k^2 \int_{D_\varepsilon(x_0)} {\rm e}^{-\imath k (\widehat{x}\cdot y)} u^{(D, \varepsilon, x_0,n)}(y) \,{\rm d}y. \nonumber \end{align} \tag{ 9 }$$

For the following consideration we introduce the Herglotz operator $H:L^2(\mathbb{S})\mapsto H^1_{\rm loc}(\mathbb{R}^d)$ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{12} (Hg)(y) :=\int_{{\mathbb S}} {\rm e}^{\imath k (\widehat{x}\cdot y)} g(\widehat x) \,{\rm d}S_x. \nonumber \end{align} \tag{ 10 }$$

Note that Hg is an entire solution to the Helmholtz equation. For a given geometry $D_\varepsilon(x_0)$ we denote by $H_{(D, \varepsilon, x_0)}: L^2(\mathbb{S})\mapsto L^2(D_\varepsilon(x_0))$ its restriction

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (H_{(D, \varepsilon, x_0)})g :=Hg\vert _{D_\varepsilon(x_0)}. \nonumber \end{align*}$$

The L²-adjoint Herglotz operator $H^*_{(D, \varepsilon, x_0)}: L^2(D_\varepsilon(x_0))\mapsto L^2(\mathbb{S})$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{13} (H^*_{(D, \varepsilon, x_0)} \varphi) (\widehat{x}) :=\int_{D_\varepsilon(x_0)} {\rm e}^{-\imath k (\widehat{x}\cdot y)} \varphi(y) \,{\rm d}y. \nonumber \end{align} \tag{ 11 }$$

Then from (9) and (11) it follows the following expression for the far-field pattern $u^{(D, \varepsilon, x_0, n)}_\infty$ of the scattered field corresponding to (1a)–(1c)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{far-exp} u^{(D,\varepsilon,x_0,n)}_\infty (\widehat{x}) =\gamma_d (n-1) k^2 (H^*_{(D, \varepsilon, x_0)} u^{(D, \varepsilon, x_0, n)}) (\widehat{x})\qquad {\rm for~}\widehat x\in {{\mathbb S}} \nonumber \end{align} \tag{ 12 }$$

which will be used later on in the proof of theorem 1.

The inverse scattering problem under consideration is to find the location x₀ of the inhomogeneity $D_\varepsilon(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.

Let the true (unknown) inhomogeneity be $D^\star_{\varepsilon^\star}(x^\star_0)$ with refractive index $n^\star$ and let us denote the corresponding far-field measurement due to a given incident field by $u^\star_\infty\in L^2(\mathbb{S})$. The latter is what we measure. For variation of the topology, a trial inhomogeneity $D_\varepsilon(x_0)$ with a trial refractive index n is put in the background medium. This constitutes a family of solutions $u^{(D, \varepsilon, x_0, n)}$ of the forward problem (1a)–(1c) for the same incident field. We denote by $u^{(D, \varepsilon, x_0, n)}_\infty$ the corresponding far field pattern. To identify the unknown inhomogeneity, we set the least square L² cost functional $\mathcal{J}: \mathfrak{G}_D\times \mathbb{R}_+\times \mathbb{R}^d \times\mathbb{C}_+\mapsto \mathbb{R}_+$ of the misfit far-field patterns by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{14} \mathcal{J}(D,\varepsilon,x_0,n) :={\frac{1}{2}} \int_{{\mathbb S}} \vert u^{(D,\varepsilon,x_0,n)}_\infty -u^\star_\infty\vert ^2 \,{\rm d}S_x, \nonumber \end{align} \tag{ 13 }$$

which should be minimized over the trial variables. This minimization problem reads: find true variables $(D^\star, \varepsilon^\star, x^\star_0, n^\star)\in \mathfrak{G}_D\times \mathbb{R}_+\times \mathbb{R}^d \times\mathbb{C}_+$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{15} \mathcal{J}(D^\star, \varepsilon^\star, x^\star_0,n^\star) =\min_{(D,\varepsilon,x_0,n)\in \mathfrak{G}_D\times \mathbb{R}_+\times \mathbb{R}^d \times\mathbb{C}_+} \mathcal{J}(D,\varepsilon,x_0,n). \nonumber \end{align} \tag{ 14 }$$

We say that the measurement $u^\star_\infty$ is feasible with respect to the quadruple $(D^\star, \varepsilon^\star, x^\star_0, n^\star)$, if the solution $u^{(D^\star, \varepsilon^\star, x^\star_0, n^\star)}$ of the corresponding forward problem obeys exactly this far-field pattern

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{16} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}_\infty (\widehat{x}) =u^\star_\infty(\widehat{x})\quad {\rm for}~ \widehat{x}\in{{\mathbb S}}. \nonumber \end{align} \tag{ 15 }$$

In the feasible case (15), the trivial minimum in (14) is attained. Note that for uniqueness results on determining the support $D^\star_{\varepsilon^\star}(x^\star_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 ${\rm Im}(n^\star)=0$ and the measurement $u^\star_\infty\in L^2(\mathbb{S})$ is feasible. Then the following necessary optimality condition for (14) holds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{17} {\rm Im}\Bigl\{\gamma_d \int_{D^\star_{\varepsilon^\star}(x^\star_0)} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)} \overline{H u^\star_\infty} \,{\rm d}x\Bigr\} =0, \nonumber \end{align} \tag{ 16 }$$

where the Herglotz operator H is defined in (10), $u^{(D^\star, \varepsilon^\star, x^\star_0, n^\star)}\in H^1_{\rm loc}(\mathbb{R}^d)$ is the solution of the forward scattering problem (1a)–(1c) for the true inhomogeneous medium, and the constant $\gamma_d$ is given in (4).

Proof. For the feasible measurement $u^\star_\infty$, after multiplication with the complex conjugate $\overline{u}$ of a test function u and integrating over the unit ball boundary ${{\mathbb S}}$, the identity (15) takes the following weak form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{18} \int_{{\mathbb S}} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}_\infty (\widehat{x})\, \overline{u}(\widehat{x}) \,{\rm d}S_x =\int_{{\mathbb S}} u^\star_\infty(\widehat{x})\, \overline{u}(\widehat{x}) \,{\rm d}S_x \nonumber \end{align} \tag{ 17 }$$

for all test-functions $u\in L^2(\mathbb{S})$.

Next we insert in (17) $u =\gamma_d H^*_{(D, \varepsilon, x_0)} u^{(D^\star, \varepsilon^\star, x^\star_0, n^\star)}$ with the adjoint Herglotz operator $H^*_{(D, \varepsilon, x_0)}$ defined in (11) and use the integral representation (9) for the far-field pattern $u^{(D^\star, \varepsilon^\star, x^\star_0, n^\star)}_\infty$ to derive from (17) that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{D} :&=\int_{{\mathbb S}} \gamma_d (n^\star-1) k^2 \left( \int_{D^\star_{\varepsilon^\star}(x^\star_0)} {\rm e}^{-\imath k (\widehat{x}\cdot y)} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}(y) \,{\rm d}y\right. \nonumber \\ &\quad\left. \overline{\gamma_d} \int_{D^\star_{\varepsilon^\star}(x^\star_0)} {\rm e}^{\imath k (\widehat{x}\cdot z)} \overline{u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}(z)} \,{\rm d}z\right) \,{\rm d}S_x \nonumber \\ &=\int_{{\mathbb S}} u^\star_\infty\, \overline{\gamma_d H^*_{(D, \varepsilon, x_0)} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}} \,{\rm d}S_x. \nonumber \end{align*}$$

Using the identity (see [5, lemma 1]):

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \int_{{\mathbb S}} {\rm e}^{\imath k \widehat{x}\cdot(z-y)} \,{\rm d}S_x =\beta_d\, {\rm Im}\{\Phi(z,y)\}\quad {\rm with}~ \beta_2 =8\pi, \beta_3 =(4\pi)^2/k \nonumber \end{align*}$$

and interchanging the order of integration we arrive at

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{19} \mathcal{D} =&\vert \gamma_d\vert ^2 \beta_d (n^\star-1) k^2 \int_{D^\star_{\varepsilon^\star}(x^\star_0)} \int_{D^\star_{\varepsilon^\star}(x^\star_0)} {\rm Im}\{\Phi(z,y)\} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}(y)\nonumber \\ &\overline{u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}(z)} \,{\rm d}z \,{\rm d}y =\int_{D^\star_{\varepsilon^\star}(x^\star_0)} (H u^\star_\infty)(z)\, \overline{\gamma_d u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}}(z) \,{\rm d}z. \nonumber \end{align} \tag{ 18 }$$

Since ${\rm Im}(n^\star)=0$, then $\overline{\mathcal{D}} =\mathcal{D}$ in the first equation in (18), which implies that the imaginary part of $\mathcal{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 $\gamma_3$ is real-valued in 3d, and hence it can be omitted from formula (16) in three dimensions.

Remark 1. The weak variational form of the optimality condition (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² cost misfit functional for given measurements.

Next we apply to the necessary optimality condition (16) the asymptotic argument as $\varepsilon^\star\searrow0^+$, i.e. the volume of the support of the inhomogeneity becomes arbitrarily small.

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 $\mathbb{R}^d$ we assign to the center of the inhomogeneity $x^\star_0$ a local d-dimensional coordinate system with the radial coordinate $\vert x-x^\star_0\vert$ and d  −  1 angular coordinates, i.e the polar coordinates in 2d, and the spherical coordinates in 3d. Furthermore $\widehat{x-x^\star_0} :=\frac{x-x^\star_0}{\vert x-x^\star_0\vert }$ denotes d-coordinate on the unit ball boundary $\partial B_1(x^\star_0)$ centered at $x^\star_0$.

Since $u^{\rm i}$ is an entire solution of the Helmholtz equation in an arbitrary ball $B_R(x^\star_0)$ of radius R  >  0 and center $x^\star_0$, it has the following expansion in the two-dimensional case (see [15, section 3.4]):

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{20} u^{\rm i}(x) =\sum_{\ell=-\infty}^\infty K_\ell J_\ell(k\vert x-x^\star_0\vert)\, \bigl( (\widehat{x-x^\star_0})_1 +\imath (\widehat{x-x^\star_0})_2 \bigr)^\ell, \nonumber \end{align} \tag{ 19 }$$

with the Fourier coefficients $\newcommand{\e}{{\rm e}} K_\ell\in\mathbb{C}$ and the Bessel functions $\newcommand{\e}{{\rm e}} J_\ell$ of order $\newcommand{\e}{{\rm e}} \ell\in\mathbb{Z}$, whereas in the three-dimensional case (see [15, section 2.4]):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{21} u^{\rm i}(x) =\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell K^m_\ell j_\ell(k\vert x-x^\star_0\vert) Y^m_\ell (\widehat{x-x^\star_0}), \nonumber \end{align} \tag{ 20 }$$

where $\newcommand{\e}{{\rm e}} K^m_\ell\in\mathbb{C}$, $\newcommand{\e}{{\rm e}} Y^m_\ell$ are the spherical harmonics, and $\newcommand{\e}{{\rm e}} j_\ell$ are the spherical Bessel functions of order $\newcommand{\e}{{\rm e}} \ell\in\mathbb{N}_0$.

From (19) we find that $u^{\rm i}(x^\star_0) =K_0 J_0(0)$, whereas from (20) that $u^{\rm i}(x^\star_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^{\rm i}$ as $\vert x-x^\star_0\vert \searrow0^+$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{22} u^{\rm i}(x) =u^{\rm i}(x^\star_0) +\mathcal{O}(\vert x-x^\star_0\vert)\;\; \quad {\rm in~}B_R(x^\star_0)~ \; {\rm with~} R>0. \nonumber \end{align} \tag{ 21 }$$

The Herglotz wave $H u^\star_\infty$ defined in (10) is an entire solution of the Helmholtz equation, thus it satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{23} (H u^\star_\infty)(x) =(H u^\star_\infty)(x^\star_0) +\mathcal{O}(\vert x-x^\star_0\vert)\;\; \quad {\rm in~} B_R(x^\star_0) \; {\rm ~with }~R>0. \nonumber \end{align} \tag{ 22 }$$

Hence, for a small size parameter $\varepsilon^\star>0$, (21) and (22) in the ball of radius $R =\varepsilon^\star$ imply that in $D^\star_{\varepsilon^\star}(x^\star_0)\subset B_{\varepsilon^\star}(x^\star_0)$ it holds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{24} u^{\rm i}(x) =u^{\rm i}(x^\star_0) +\mathcal{O}(\varepsilon^\star),\quad (H u^\star_\infty)(x) =(H u^\star_\infty)(x^\star_0) +\mathcal{O}(\varepsilon^\star). \nonumber \end{align} \tag{ 23 }$$

The Born approximation stated in (7) and (8) now takes the form

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{25} u^{(D^\star, \varepsilon^\star, x^\star_0,n^\star)}(x) &=u^{\rm i}(x) +\mathcal{O}\bigl( \vert n^\star-1\vert k^2 \Vert T_{(D^\star, \varepsilon^\star, x^\star_0)}\Vert \bigr)\nonumber \\ &=u^{\rm i}(x) +\mathcal{O}\bigl( \vert n^\star-1\vert k^2 (\varepsilon^\star)^2 \vert \ln\varepsilon^\star\vert ^{3-d} \bigr).\nonumber \end{align} \tag{ 24 }$$

Plugging the asymptotic relations (23) and (24) into the necessary optimality condition (16) proves straightforwardly the following theorem.

Theorem 2. Assume that ${\rm Im}(n^\star)=0$ and the measurement $u^\star_\infty\in L^2(\mathbb{S})$ is feasible. Then the necessary optimality condition (16) has the asymptotic form as $\varepsilon^\star\searrow0^+$:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{26} 0 ={\rm Im}\bigl\{u^{\rm i}(x^\star_0) \gamma_d \overline{(H u^\star_\infty)(x^\star_0)} \bigr\}+\mathcal{O}\bigl( \varepsilon^\star +\vert n^\star -1\vert k^2 (\varepsilon^\star)^2 \vert \ln\varepsilon^\star\vert ^{3-d} \bigr). \nonumber \end{align} \tag{ 25 }$$

Remark 2. The principal asymptotic term in (25) is determined by the incident field $u^{\rm i}$ which is an entire solution to the Helmholtz equation and by the Herglotz wave function $H u^\star_\infty$ 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 $\mathcal{J}$ in (13) as follows

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{J}(D,\varepsilon,x_0,n) &=\mathcal{J}_0 -{\rm Re}\Bigl\{\int_{{\mathbb S}} u^{(D,\varepsilon,x_0,n)}_\infty \overline{u^\star_\infty} \,{\rm d}S_x \Bigr\} +{\frac{1}{2}} \int_{{\mathbb S}} \vert u^{(D,\varepsilon,x_0,n)}_\infty\vert ^2 \,{\rm d}S_x,\nonumber \\ &{\rm where}\; \mathcal{J}_0 := {\frac{1}{2}} \int_{{\mathbb S}} \vert u^\star_\infty\vert ^2 \,{\rm d}S_x. \nonumber \end{align*}$$

Interchanging the order of integration and using (9) implies that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \int_{{\mathbb S}} u^{(D,\varepsilon,x_0,n)}_\infty \overline{u^\star_\infty} \,{\rm d}S_x =\gamma_d (n-1) k^2 \int_{D_\varepsilon(x_0)} u^{(D,\varepsilon,x_0,n)} \overline{H u^\star_\infty} \,{\rm d}y. \nonumber \end{align*}$$

Applying the asymptotic formula (23) and (24) we now get the two-parameter $(\varepsilon, n)$-asymptotic representation of the cost functional ${{\mathcal J}}$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{27} \mathcal{J}(D,\varepsilon,x_0,n) &=\mathcal{J}_0 -k^2 \vert D_\varepsilon(x_0)\vert \, {\rm Re}\bigl\{(n-1) u^{\rm i}(x_0) \gamma_d \overline{(H u^\star_\infty)(x_0)}\bigr\}\nonumber \\ &+\mathcal{O}\bigl( \vert D_\varepsilon(x_0)\vert \vert n-1\vert k^2 (\varepsilon +\vert n-1\vert k^2 \varepsilon^2 \vert \ln\varepsilon\vert ^{3-d}) \nonumber \\ &+\vert D_\varepsilon(x_0)\vert ^2 \vert n-1\vert ^2 k^4 \bigr), \nonumber \end{align} \tag{ 26 }$$

where we have assumed that $\vert n-1\vert k^2 \Vert T_{(D, \varepsilon, x_0)}\Vert <1$ according to (8), and $\vert D_\varepsilon(x_0)\vert =\varepsilon^d \vert D\vert$ for the Hausdorff measure of the sets in $\mathbb{R}^d$.

Hence from (26) it follows the following expression of the topological derivative of $\mathcal{J}$:

Theorem 3. For fixed $n\in\mathbb{C}_+$ and $\varepsilon\searrow0^+$ we have that

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{28} \mathcal{J}^\prime(D,x_0,n) :=\lim_{\varepsilon\searrow 0^+} {\frac{\mathcal{J}(D,\varepsilon,x_0,n) -\mathcal{J}_0}{\varepsilon^d}} =-k^2 \vert D\vert \, {\rm Re}\bigl\{(n-1) u^{\rm i}(x_0) \gamma_d \overline{(H u^\star_\infty)(x_0)} \bigr\}. \nonumber \end{align} \tag{ 27 }$$

Now let $(D^\star, \varepsilon^\star, x^\star_0, n^\star)\in \mathfrak{G}_D\times \mathbb{R}_+\times \mathbb{R}^d \times\mathbb{C}_+$ be the solution of the topology optimization problem (14). This means that the optimal value of the cost functional $\mathcal{J} (D^\star, \varepsilon^\star, x^\star_0, n^\star) =0$ in (26). Varying the imaginary part of the optimal refractive index ${\rm Im}(n^\star)$, we observe that the second term in the expression of the topological derivative from (27):

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{J}^\prime(D^\star,x^\star_0,n^\star) =-k^2 \vert D^\star\vert \Bigl( &{\rm Re}(n^\star-1)\, {\rm Re}\bigl\{u^{\rm i}(x^\star_0) \gamma_d \overline{(H u^\star_\infty)(x^\star_0)} \bigr\}\nonumber \\ -&{\rm Im}(n^\star)\, {\rm Im}\bigl\{u^{\rm i}(x^\star_0) \gamma_d \overline{(H u^\star_\infty)(x^\star_0)} \bigr\} \Bigr) \nonumber \end{align*}$$

should be asymptotically zero according to theorem 2.

Based on theorem 2, we suggest to consider the real-valued imaging operator $\mathcal{I}:L^2(\mathbb{S}) \mapsto C(\mathbb{R}^d)$ which, for given incident field $u^{\rm i}$ is defined over the space of far-field patterns $u^\star_\infty$ as follows

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{29} (\mathcal{I} u^\star_\infty)(x) :={\rm Im}\bigl\{u^{\rm i}(x) \gamma_d \overline{(H u^\star_\infty)(x)} \bigr\}\quad {\rm for~}x\in \mathbb{R}^d. \nonumber \end{align} \tag{ 28 }$$

Due to (25), the continuous function $\mathcal{I} u^\star_\infty$ possesses the property that $(\mathcal{I} u^\star_\infty)(x^\star_0) =0$ asymptotically at the center $x^\star_0$ of the unknown inhomogeneity $D^\star_{\varepsilon^\star}(x^\star_0)$. This implies that the zero-level set $\mathcal{L}_{=0}$ of $\mathcal{I} u^\star_\infty$ contains (asymptotically) the inhomogeneity center $x^\star_0$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{30} x^\star_0\in \mathcal{L}_{=0}(\mathcal{I} u^\star_\infty) :=\{x\in\mathbb{R}^d:\; (\mathcal{I} u^\star_\infty)(x) =0\}. \nonumber \end{align} \tag{ 29 }$$

Consequently, this suggest the test center $x^\star_0$ is looked for (asymptotically) as the intersection point of d zero-level sets:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{31} x^\star_0 =\mathcal{L}_{=0} (\mathcal{I} u^\star_{\infty,1}) \cap\dots\cap \mathcal{L}_{=0} (\mathcal{I} u^\star_{\infty,d}) \nonumber \end{align} \tag{ 30 }$$

from $(u^\star_{\infty, 1}, \dots, u^\star_{\infty, 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.

5.1. 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, \varepsilon, x_0)}$ given in (6).

Without loss of generality, let $\Omega =(0, 1){\hspace{0pt}}^d$ be the finite computational test domain, that is the unit square in 2d and the unite cube in 3d, containing a homogeneity $D_\varepsilon(x_0)$ with a refractive index n. For a fixed mesh size $h\in\mathbb{R}_+$, we denote the nodal points of the uniform grid (quadrilateral in 2d and polyhedral in 3d) in the closure $\overline{\Omega}$ by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle N_h :=\{x^{\,j}=h z,\, z\in\mathbb{Z}^d: \quad x^{\,j}\in \overline{\Omega} \}. \nonumber \end{align*}$$

Based on the fast solver from [43], we apply to (5) and (6) the following formula of numerical integration over $x^{\,j}\in N_h$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{32} u^{(D, \varepsilon, x_0,n)}_h(x^{\,j}) =u^{\rm i}(x^{\,j}) +(n-1) k^2 h^d \sum_{\scriptsize\begin{array}{@{}rrclll@{}} x^l\in N_h\cap \overline{D_\varepsilon(x_0)} \nonumber \\ \,x^l\not =x^{\,j} \qquad \end{array}} \Phi(x^{\,j},x^l) u^{(D, \varepsilon, x_0,n)}_h(x^l). \nonumber \end{align} \tag{ 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 $\mathcal{O}(h^2\vert \ln(h)\vert)$-error estimate in the max-norm:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{33} {\rm Error} :=\max_{x^{\,j}\in N_h} \bigl\vert u^{(D, \varepsilon, x_0,n)}(x^{\,j}) -u^{(D, \varepsilon, x_0,n)}_h(x^{\,j})\bigr\vert \nonumber \end{align} \tag{ 32 }$$

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 $\varepsilon=2^{-4}$ centered at x₀  =  (0.25,0.75), with the complex-valued refractive index $n=2+\imath$, which is illuminated in the direction of $\theta=\frac{\pi}{3}$ by the plane wave $u^{\rm i}(x) ={\rm e}^{\imath k(x\cdot\widehat{z})}$ with $\widehat{z} =(\cos\theta, \sin\theta)$. The wave numbers are taken in the range of $k\in[2, \dots, 2^7]$, and the mesh size $h\in[2^{-7}, \dots, 2^{-3}]$. Since there is no any analytical solution available, we substitute the exact solution $u^{(D, \varepsilon, x_0, n)}$ with the discrete solution $u^{(D, \varepsilon, x_0, n)}_{\widetilde{h}}$ calculated on the finest mesh with $\widetilde{h}=2^{-7}$. In figure 1 the values of ${\rm error}$ are depicted in dependence on the parameter kh, respectively:

• for fixed waves numbers $k=\{2^1, 2^2, 2^3\}$ when varying h in plot (a),
• for fixed mesh sizes $h=\{2^{-6}, 2^{-5}, 2^{-4}, 2^{-3}\}$ when varying k in plot (b),
• for fixed $kh=\{2^{-2}, 2^{-1}, 2^0\}$ when increasing k in plot (c).

Zoom In Zoom Out Reset image size

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 ${\rm Error}(kh)$ for wave numbers k and mesh sizes h chosen such that $kh\leqslant1$.

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 $\widehat{x}\in{{\mathbb S}}$ in $\mathbb{R}^d$. Then the discretization of the far-field (see [43, section 3.5]) for $\widehat{x}\in M$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{34} (u^{(D, \varepsilon, x_0,n)}_h)_\infty (\widehat{x}) =\gamma_d (n -1) k^2 h^d \sum_{x^{\,j}\in N_h\cap \overline{D_\varepsilon(x_0)}} {\rm e}^{-\imath k (\widehat{x}\cdot x^{\,j})} u^{(D, \varepsilon, x_0,n)}_h(x^{\,j}). \nonumber \end{align} \tag{ 33 }$$

Finally, for far-field measurement $u_\infty$ given at nodes $\widehat{x}\in M$, the Herglotz operator (10) can be discretized for $x^{\,j}\in N_h$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{35} (H_M u_\infty) (x^{\,j}) =\sum_{\widehat{x}\in M} w_M(\widehat{x})\, {\rm e}^{\imath k (\widehat{x}\cdot x^{\,j})} u_\infty (\widehat{x}) \quad {\rm with~ }w_M(\widehat{x})\in\mathbb{R}. \nonumber \end{align} \tag{ 34 }$$

Proper weights $w_M(\widehat{x})$ in (34) should be determined depending on meshing of ${{\mathbb S}}$. For example, on the unit circle in 2d we use equidistant nodes $\widehat{x}^l =(\cos\theta_l, \sin\theta_l)$ with $\theta_l=\frac{2\pi l}{N_\theta}$ for $l=1, \dots, N_\theta$ and fixed $N_\theta\geqslant1$, then all weights $w_M(\widehat{x}^l)=\frac{2\pi}{N_\theta}$. For the uniform latitude-longitude grid on the unit sphere in 3d such that $\widehat{x}^{(l, m)} =(\cos\theta_l\sin\phi_m, \sin\theta_l\sin\phi_m, \cos\phi_m)$ with $\theta_l$ as above and $\phi_m=\frac{\pi(m-1)}{N_\phi-1}$ for $m=1, \dots, N_\phi$ and fixed $N_\phi\geqslant2$, we use the cubature formula based on the trapezoidal rule with the positive weights

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle &w_M(\widehat{x}^{(l,m)}) ={\frac{2\pi}{N_\theta}} \sin^2({\frac{\pi}{2(N_\phi-1)}}) \quad {\rm for~} m=\{1,N_\phi\},\nonumber \\ {\rm and~if~}N_\phi\geqslant3:\quad &w_M(\widehat{x}^{(l,m)}) ={\frac{2\pi}{N_\theta}} \sin({\frac{\pi}{N_\phi-1}}) \sin({\frac{\pi(m-1)}{N_\phi-1}})\quad {\rm for }~m=\{2,\dots,N_\phi-1\}. \nonumber \end{align*}$$

5.2. 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^\star_{\varepsilon^\star}(x^\star_0)$ and refractive index $n^\star$, we use the numerical solution $u^{(D^\star, \varepsilon^\star, x_0^\star, n^\star)}_{h^\star}$ given by (31), in order to synthesize the discrete far-field measurement $u^\star_\infty =(u^{(D^\star, \varepsilon^\star, x_0^\star, n^\star)}_{h^\star})_\infty$ according to formula (33) and the corresponding Herglotz wave $H_M u^\star_\infty$ as per (34). In this setting from the discrete far-field measurement $u^\star_\infty(\widehat{x})$ given at nodes $\widehat{x}\in M$ and the incident field $u^{\rm i}(x^{\,j})$ known at mesh points $x^{\,j}\in N_h$, depending on these two grids we calculate the discrete imaging function according to (28) for $x^{\,j}\in N_h$ as follows

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{36} (\mathcal{I}_{(h,M)} u^\star_\infty)(x^{\,j}) ={\rm Im}\bigl\{u^{\rm i}(x^{\,j}) \gamma_d \overline{(H_M u^\star_\infty)(x^{\,j})} \bigr\}. \nonumber \end{align} \tag{ 35 }$$

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^\star$ is L-shaped of size $\varepsilon^\star=2^{-3}$ centered at $x_0^\star =(0.125, 0.375)$, with the real-valued refractive index $n^\star=0.5$. The computation is done over the unit square $\Omega=(0, 1){\hspace{0pt}}^2$ with the mesh sizes $h=h^\star=2^{-6}$, and $N_\theta=10$ equidistant points on the unit circle (relatively few measurement points).

Zoom In Zoom Out Reset image size

Two far-field patterns $u^\star_{\infty, 1}$ and $u^\star_{\infty, 2}$ correspond to the case when the inhomogeneity is illuminated by a plane wave $u^{\rm i}(x) ={\rm e}^{\imath k(x\cdot\widehat{z})}$ with $\widehat{z} =(\cos\theta, \sin\theta)$ in the direction of $\theta_1=\frac{\pi}{6}$ and $\theta_2=\frac{\pi}{4}$, respectively. Plots (a) and (b) depict the corresponding imaging functions $\mathcal{I}_{(h, M)} u^\star_{\infty, 1}$ and $\mathcal{I}_{(h, M)} u^\star_{\infty, 2}$ together with the zero-level sets $\mathcal{L}_{=0}(\mathcal{I}_{(h, M)} u^\star_{\infty, 1})$ and $\mathcal{L}_{=0}(\mathcal{I}_{(h, M)} u^\star_{\infty, 2})$, which form straight lines crossing the inhomogeneity. The zero-level sets are depicted also in the computational test domain Ω in the plot (c), their intersection point:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{37} x^{(h,M)}_0 =\mathcal{L}_{=0}(\mathcal{I}_{(h,M)} u^\star_{\infty,1}) \cap\mathcal{L}_{=0}(\mathcal{I}_{(h,M)} u^\star_{\infty,2}) \nonumber \end{align} \tag{ 36 }$$

is unique and close (to become quantitative measure later) to the test center $x_0^\star$. Moreover, $x^{(h, M)}_0=x_0^\star$ exactly, if the size is sufficiently small such that $\varepsilon^\star\leqslant2^{-6}$. Here zero-level sets $\mathcal{L}_{=0}(\mathcal{I}_{(h, M)} u^\star_{\infty})$ are utilized numerically using the narrow strip technique, see e.g. [33]. More specifically, in a narrow strip of nodes where $\mathcal{I}_{(h, M)} u^\star_{\infty}$ 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^\star}$ in the narrow strip. It is important to stress that the wave number is chosen small: k  <  k₁ below the threshold ($k_1 \approx2.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 $\theta_1=\frac{\pi}{6}$, $\theta_2=\frac{\pi}{4}$, and $\theta_3=\pi$.

Zoom In Zoom Out Reset image size

Since the asymptotic expansions in section 4 are carried out in the near-field of the inhomogeneity center $x_0^\star$, 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 $\varepsilon^\star_{\circ}=2^{-6}$ centered at $x^\star_{0, {\circ}} =(0.25, 0.75)$ with the refractive index $n^\star_{\circ}=2$, the triangle-shaped of size $\newcommand{\tr}{{\rm tr}} \varepsilon^\star_\triangle=2^{-5}$ centered at $\newcommand{\tr}{{\rm tr}} x^\star_{0, \triangle} =(2/3, 1/3)$ with $\newcommand{\tr}{{\rm tr}} n^\star_\triangle=10$, and the L-shaped of size $\varepsilon^\star_{\rm L}=2^{-4}$ centered at $x^\star_{0, {\rm L}} =(0.125, 0.375)$ with $n^\star_{\rm L}=0.5$. Here the wave number k  =  1 and three plane waves propagate along the directions of $\theta\in\{\frac{\pi}{6}, \frac{\pi}{3}, \pi\}$. 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 $\newcommand{\tr}{{\rm tr}} n^\star_{\circ} =n^\star_\triangle= n^\star_{\rm 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We carried out a series of numerical tests for inhomogeneities $D^\star_{\varepsilon^\star}(x_0^\star)$ of various shapes, namely L-shape, ball, triangle, point (in the last case $\varepsilon^\star=0$), varying the location in the test domain $\Omega=(0, 1){\hspace{0pt}}^2$, illumination directions, the wave number, and the mesh sizes in the range of $h, h^\star\in[2^{-6}, \dots, 2^{-3}]$. Although our theory covers only real refractive index $n^\star$ 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 $\vert x^{(h, M)}_0 -x_0^\star\vert$ 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\leqslant1+\delta$ with $\delta\geqslant 0$ depending on the problem parameters.
• (ii)  
  The identification error depends super-linearly on the mesh size $h^\star$ of the grid used to synthesize the far-field measurement $u^\star_{\infty}$. It is not sensitive to the mesh size $h\leqslant2^{-3}$ used for the imaging function, and to the grid M if chosen uniformly on the unit circle with $N_\theta\geqslant5$ nodes.
• (iii)  
  If the test center $x_0^\star$ coincides with a node of the grid $M_{h^\star}$, this leads to the better identification result than $x_0^\star\not\in M_{h^\star}$. For comparison see the two curves of the error depicted in the plot (a) in figure 6 in terms of the size $\varepsilon^\star$.
• (iv)  
  The identification error $\vert x^{(h, M)}_0 -x_0^\star\vert$ drops super-linearly with respect to decreasing the size $\varepsilon^\star$ of the unknown inhomogeneity as shown in figure 6(a). For comparison, the error-curves for the ball of size $\varepsilon^\star=2^{-5}$ and for the equilateral triangle of size $\varepsilon^\star=2^{-6}$ are presented through the other plots (b)–(d).
• (v)  
  Figures 6(b) and (c) show the error with respect to the refractive index, namely, for ${\rm Re}(n^\star)\in[0.01, 100]$ and ${\rm Im}(n^\star)=0$, and ${\rm Im}(n^\star)\in[0, 0.1]$ when either ${\rm Re}(n^\star)=2$ for $\varepsilon^\star=2^{-5}$ or ${\rm Re}(n^\star)=10$ for $\varepsilon^\star=2^{-6}$. We observe here that the choice of $\varepsilon^\star=2^{-6}$ provides $\vert x^{(h, M)}_0 -x_0^\star\vert <h^\star$, where the value of $h^\star$ is drawn with dotted line through the plots.
• (vi)  
  In figure 6(d) the stability of the identification is tested with respect to a Gaussian noise with the standard deviation $\sigma=1-12\%$ computed over 200 realizations of the far-field measurement. The error here is less then $4\%$ for $\varepsilon^\star=2^{-6}$.

Zoom In Zoom Out Reset image size

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.

Finally, in figure 7 we show a three-dimensional example, namely the identification of the center of a ball-shaped inhomogeneity of size $\varepsilon^\star=2^{-5}$ centered at $x_0^\star =(0.75, 0.375, 0.125)$ with the refractive index $n^\star=2$. Computations are done for the mesh sizes $h=h^\star=2^{-4}$ over the uniform grid with $N_\theta=10$ longitude and $N_\phi=10$ latitude lines on the unit sphere. The inhomogeneity is illuminated by 3 plane waves $u^{\rm i}(x) ={\rm e}^{\imath k((x-(1/2, 1/2, 1/2))\cdot\widehat{z})}$ with $k=1/\sqrt{2}$ and $\widehat{z} =(\cos\theta \sin\phi, \sin\theta \sin\phi, \cos\phi)$ in three perpendicular directions of $(\theta_1, \phi_1)=(0, 0)$, $(\theta_2, \phi_2)=(0, \frac{\pi}{2})$, and $(\theta_3, \phi_3)=(\frac{\pi}{2}, \frac{\pi}{2})$. In the plots (a)–(c) we show the corresponding zero-level sets forming planes which cross the inhomogeneity, and their intersection point $x^{(h, M)}_0$ coincides with the center of the ball as shown in the plot (d).

Zoom In Zoom Out Reset image size

FC is supported by AFOSR grant FA9550-17-1-0147, NSF Grant DMS-1602802 and Simons Foundation Award 392261.

VK is supported by the Austrian Science Fund (FWF) project P26147-N26: 'Object identification problems: numerical analysis' (PION) and the Austrian Academy of Sciences (OeAW).