Near field imaging of small isotropic and extended anisotropic scatterers
Isaac Harris111Department of Mathematics, Texas A $\&$ M University, College
Station, Texas 77843-3368, (Corresponding Author) E-mail:
<EMAIL_ADDRESS>and Scott Rome222Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider two time-harmonic inverse scattering problems of
reconstructing penetrable inhomogeneous obstacles from near field
measurements. First we appeal to the Born approximation for reconstructing
small isotropic scatterers via the MUSIC algorithm. Some numerical
reconstructions using the MUSIC algorithm are provided for reconstructing the
scatterer and piecewise constant refractive index using a Bayesian method. We
then consider the reconstruction of an anisotropic extended scatterer by a
modified linear sampling method and the factorization method applied to the
near field operator. This provides a rigorous characterization of the support
of the anisotropic obstacle in terms of a range test derived from the measured
data. Under appropriate assumptions on the material parameters, the derived
factorization can be used to determine the support of the medium without
a-priori knowledge of the material properties.

Keywords: inverse scattering, factorization method, MUSIC algorithm.

AMS subject classifications: 35J05, 45Q05, 78A46

##  1 Introduction

Qualitative/Sampling methods have been used to solve multiple inverse problems
such as parameter identification and shape reconstruction. These methods have
been used to solve inverse boundary value problems for elliptic, parabolic and
hyperbolic partial differential equation. See [21] and [20] for examples of
the linear sampling methods applied to parabolic and hyperbolic equations,
respectively. The vast available literature is a representative of the many
directions that this research has taken (see [8], [15], [16], [25] and the
references therein).

The factorization method is a qualitative method which develops a range test
to determine the support of the unknown scatterer. In general, we have that

$\phi_{z}\in\text{Range}\big{(}N_{\sharp}^{1/2}\big{)}\iff z\,\,\text{is
inside the obstacle}$

where $\phi_{z}$ is known and $N_{\sharp}$ is a positive selfadjoint compact
operator defined by the measurements. By appealing to the range test, we have
that any two obstacles admitting the same data must be equal, proving
uniqueness of the inverse problem. Since the support of the scatterer is
connected to the range of a compact operator which is known, evaluating the
indicator function derived from the factorization method amounts to applying
Picard s criteria, which only requires the singular values and functions of a
known operator. This makes the factorization method computationally cheap to
implement and analytically rigorous. Using the characterization of the
obstacle given by the factorization method we propose a modified linear
sampling method as an alternative sampling algorithm for reconstructing the
obstacle. The indicator function proposed for the modified linear sampling
method can be shown to satisfy a similar analytically rigors characterization
of the obstacle.

In this paper, we consider two problems associated with inverse scattering.
The first problem we consider is the reconstruction of small isotropic
obstacles and the refractive index from near field data. To do so, we derive
the Multiple-Signal-Classification (MUSIC) algorithm which gives an
analytically rigorous method of reconstructing small obstacles. Once we have
reconstructed the obstacles, we will give a few example of how one can
reconstruct constant refractive indices using a Bayesian method. Next we
consider the reconstructing an anisotropic obstacle from near field data via
the factorization method and modified linear sampling method. This problem has
been considered in [9] and [26] for far field measurements using the linear
sampling and factorization methods, respectively. The factorization method for
an inhomogeneous isotropic media with near field measurements was studied in
[22]. We analyze the factorization method applied to the near field operator
for an anisotropic obstacle to derive a simple numerical algorithm to
reconstruct the scatterer. For the case of an anisotropic media, the unique
determination of the support is optimal because a matrix-valued coefficient is
in general not uniquely determined by the scattering data (see for e.g. [18]).
The analysis done in this paper is of a similar flavor to the work done in
[12] and [29].

##  2 Recovering small volume isotropic scatterers and the refractive index

###  2.1 Scattering by small isotropic scatterers

In this section, we formulate the direct time-harmonic scattering problem in
$\mathbb{R}^{d}$ for small isotropic scatterers. The scattering obstacle may
be made up of multiple simply components, but is not necessarily simply
connected. We are interested in using near field measurements for ‘small’
obstacles embedded in a homogenous media. Let $u^{i}(\cdot\,,y)$ be the
incident field that originates from the source point $y$ located on the curves
$C$. We will assume that the curve $C$ is class $\mathcal{C}^{2}$ and that the
incident field satisfies the Helmholtz equation in
$\mathbb{R}^{d}\setminus\\{y\\}$ for $d=2$ or $3$. Now let

$D=\bigcup\limits_{m=1}^{M}D_{m}\quad\textrm{ with }\quad
D_{m}=\big{(}z_{m}+\varepsilon B_{m}\big{)},$

where $\varepsilon>0$ is a small parameter and $B_{m}$ is a domain with a
piece-wise smooth boundary that is centered at the origin. Outside of the
scattering obstacle(s) $D_{m}$, the material parameters are homogeneous
isotropic with refractive index scaled to one. Here we denote $n(x)$ the
refractive index of the homogenous background with the (possibly
inhomogeneous) obstacle $D$ in ${\mathbb{R}}^{d}$ by

$n(x)=1\big{(}1-\chi_{D}\big{)}+\sum\limits_{m=1}^{M}n_{m}(x)\chi_{D_{m}}$

where $\chi$ is the characteristic function. The refractive index of the
obstacle(s) $D_{m}$ is represented by a continuous scalar function $n_{m}$
such that the real and imaginary parts satisfy

$\inf\limits_{m=1\cdots M}\Re\big{(}n_{m}(x)\big{)}\geq
n_{min}>0\quad\text{and}\quad\Im\left(n_{m}\right)\geq 0$

for almost every $x\in D_{m}$ where we assume that $n_{min}>1$ or $n_{min}<1$
. Now the radiating scattered field $u^{s}(\cdot\,,y)\in
H^{1}_{loc}(\mathbb{R}^{d})$ given by the point source incident field
$u^{i}(\cdot\,,y)=\Phi(\cdot\,,y)$ is the unique solution to

$\displaystyle\Delta_{x}u^{s}+k^{2}nu^{s}=k^{2}(1-n)u^{i}\quad\textrm{ in
}\quad\mathbb{R}^{d}$ (1)
$\displaystyle{\partial_{r}u^{s}}-iku^{s}=\mathcal{O}\left(\frac{1}{r^{(d+1)/2}}\right)\quad\text{
as }\quad r\rightarrow\infty$ (2)

where $r=|x|$ and the radiation condition (2) is satisfied uniformly in all
directions $\hat{x}=x/|x|$. Recall the fundamental solution to Helmholtz
equation

$\Phi(x,y)=\left\\{\begin{array}[]{lr}\frac{\text{i}}{4}H^{(1)}_{0}(k|x-y|)&\,\text{in}\,\,\mathbb{R}^{2}\\\
&\\\
\frac{1}{4\pi}\frac{\text{exp}({\text{i}k|x-y|})}{|x-y|}&\,\text{in}\,\,\mathbb{R}^{3},\end{array}\right.$

where $H^{(1)}_{0}$ is the first kind Hankel function of order zero. The
corresponding total field for the scattering problem is given by
$u(\cdot\,,y)=u^{s}(\cdot\,,y)+u^{i}(\cdot\,,y)$ satisfies

$u^{(-)}=u^{(+)}\quad\textrm{ and }\quad\frac{\partial
u^{(-)}}{\partial\nu}=\frac{\partial u^{(+)}}{\partial\nu}\quad\textrm{ on
}\partial D_{m},$

where the superscripts $+$ and $-$ for a generic function indicates the trace
on the boundary taken from the exterior or interior of its surrounding domain,
respectively. We have that the scattered field is given by the solution to the
Lippmann-Schwinger Integral Equation (see for e.g. [15])

$u^{s}(x,y)=k^{2}\int_{D}\big{(}n(z)-1\big{)}\Phi(x,z)u(z,y)\,dz.$

Assume that $\varepsilon$ is sufficiently small such that

$k^{2}\left|\int_{D}\big{(}n(z)-1\big{)}\Phi(x,z)\,dz\right|\ll 1.$

Therefore, we can conclude that the Born approximation, which is the first
term in the Neumann series solution to the Lippmann-Schwinger Integral
Equation, given by

$\displaystyle
u_{B}^{s}(x,y)=k^{2}\int_{D}\big{(}n(z)-1\big{)}\Phi(x,z)u^{i}(z,y)\,dz$ (3)

is a ‘good’ approximation for the scattered field $u^{s}(x,y)$. Therefore, we
assume that we have the measured near field data that is given by the Born
approximation $u_{B}^{s}(x,y)$. Let $\Omega$ be the bounded simply connected
open set such that satisfies $\overline{D}\subset\Omega$ and
$C=\partial\Omega$. Now assume that we have the data set of near field
measurements $u_{B}^{s}(x,y)$ for all $x,y\in C$. The inverse problem we
consider is to reconstruct the scatterer $D$ and the refractive index $n$ from
a knowledge of the Born approximation of the near field measurements.

###  2.2 The MUSIC Algorithm

The MUSIC Algorithm can be considered as the discrete analogue of the
factorization method (see for e.g. [14]). In particular, we connect the
locations of the obstacles $D_{m}$ given by $\\{z_{m}\,:\,m=1,\dots,M\\}$ to
the range of the so-called multi-static response matrix denoted by ${\bf N}$
that is defined by the near field measurements $u_{B}^{s}(x,y)$ for all
$x,y\in C$. To arrive at the MUSIC algorithm, we will need to define our
measurements operator. To this end, we let the incident field be give by
$u^{i}(x,y)=\Phi(x,y)$. Notice that by (3) we obtain that

$\displaystyle
u_{B}^{s}(x,y)=k^{2}\int_{D}\big{(}n(z)-1\big{)}\Phi(x,z)\Phi(z,y)\,dz$
$\displaystyle\qquad\qquad=\sum\limits_{m=1}^{M}k^{2}\big{(}n(z_{m})-1\big{)}\varepsilon^{d}|B_{m}|\Phi(x,z_{m})\Phi(z_{m},y)+o(\varepsilon^{d})$

since $x,y\in C$ the integrand is continuous with respect to the variable
$z\in D$. Assume we have a finite number of incident and observation
directions where $N\geq M$, with $M$ being the number of small obstacles and
$N$ being the number of incident and observation directions given by
$x_{i}\,,y_{j}\in C$. Now define the multi-static response matrix

${\bf
N}_{i,j}=\sum\limits_{m=1}^{M}k^{2}\big{(}n(z_{m})-1\big{)}\varepsilon^{d}|B_{m}|\Phi(x_{i},z_{m})\Phi(z_{m},y_{j}).$

Notice that for $\varepsilon$ sufficiently small we have that the multi-static
response matrix is approximated by $u^{s}_{B}(x_{i}\,,y_{j})$ therefore we
assume that ${\bf N}$ is known. For convenience, we will assume that the
sources and receivers are placed at the same points (i.e. $x_{i}=y_{i}$). We
now factorize the multi-static response matrix, therefore we define the
matrices ${\bf U}\in\mathbb{C}^{N\times M}$ and
${\bf\Sigma}\in\mathbb{C}^{M\times M}$, where the matrix ${\bf U}$ is given by

${\bf U}_{i,m}=\Phi(x_{i},z_{m})$

and the matrix ${\bf\Sigma}=\text{diag}(\sigma_{m})$ where
$\sigma_{m}=k^{2}\big{(}n(z_{m})-1\big{)}\varepsilon^{d}|B_{m}|\neq 0$.
Therefore, we have that ${\bf N}={\bf U}{\bf\Sigma}{\bf U}^{\top}$ and it
follows that

${\bf N}{\bf N}^{*}={\bf U}{\widetilde{\bf\Sigma}}{\bf U}^{*}\quad\text{ with
}\quad{\widetilde{\bf\Sigma}}={\bf\Sigma}{\bf U}^{\top}\,\overline{{\bf
U}}{\bf\Sigma}^{*}.$ (4)

Now define the vector ${\bf\phi}_{z}\in\mathbb{C}^{N}$ for any point
$z\in\mathbb{R}^{d}$ by

$\displaystyle{\bf\phi}_{z}=\big{(}\,\Phi(x_{1},z),\,\dots,\Phi(x_{N},z)\,\big{)}^{\top}.$
(5)

The goal of this section is to show that ${\bf\phi}_{z}$ is in the range of
the matrix ${\bf N}{\bf N}^{*}$ if and only if
$z\in\\{z_{m}\,:\,m=1,\dots,M\\}$, which is a discrete reformulation of the
result of the factorization method stated in the Introduction. Since we are
interested in finding obstacles $D_{m}$, it is sufficient to prove the result
only for values of $z\in\Omega$. Recall that without loss of generality, we
assume that the sources and receivers are placed at the same locations. We
will follow the analytic framework in Section 4.1 of [25]. We now give a
result that can be proven by using standard arguments from linear algebra (see
proof of Theorem 3.1 in [17] for details).

######  Lemma 2.1.

Let the matrix ${\bf N}$ have the following factorization ${\bf N}={\bf
U}{\bf\Sigma}{\bf U}^{\top}$ where ${\bf U}\in\mathbb{C}^{N\times M}$ and
${\bf\Sigma}\in\mathbb{C}^{M\times M}$ with $N>M$. Assume that

  1. 1.

the matrix ${\bf U}$ has full rank $M$

  2. 2.

the matrix ${\bf\Sigma}$ is invertible

then Range$\big{(}{\bf U}\big{)}=$Range$\big{(}{\bf N}{\bf N}^{*}\big{)}$.

We now construct an indicator function derived from Lemma 2.1 to reconstruct
the locations of the scattering obstacles. To this end, for each sampling
point $z\in\Omega$ we will show that the vector ${\bf\phi}_{z}$ is in the
range of ${\bf N}{\bf N}^{*}$ if and only if
$z\in\\{z_{m}\,:\,m=1,\dots,M\\}$. We now prove an auxiliary result that
connects the location of the scattering obstacles to the range of the matrix
${\bf U}$.

######  Theorem 2.1.

Assume that the set $S=\\{{x}_{i}\in C\,:\,i\in\mathbb{N}\\}$ is dense in $C$.
Let $z\in\Omega$ then there is a number $N_{0}\in\mathbb{N}$ such that for all
$N\geq N_{0}$, we have that

  1. 1.

the matrix ${\bf U}$ has full rank,

  2. 2.

${\bf\phi}_{z}\in$ Range$\big{(}{\bf U}\big{)}$ if and only if
$z\in\\{z_{m}\,:\,m=1,\dots,M\\}$.

###### Proof.

It is clear that ${\bf\phi}_{z_{m}}$ is in the range of ${\bf U}$ since
${\bf\phi}_{z_{m}}$ is the $m$-th column of ${\bf U}.$ To prove that for
$z\in\Omega\setminus\\{z_{m}\,:\,m=1,\dots,M\\}$ that ${\bf\phi}_{z}$ is not
in the range of ${\bf U}$ for some $N$ sufficiently large, we proceed by way
of contradiction. Let $z\in\Omega\setminus\\{z_{m}\,:\,m=1,\dots,M\\}$ and
assume that there does not exist such a $N_{0}$, then we have that there
exists $\alpha_{m}^{N}$, $\alpha^{N}\in\mathbb{C}$ such that

$\big{|}\alpha^{N}\big{|}+\sum\limits_{m=1}^{M}\left|\alpha_{m}^{N}\right|=1$
(6)

and

$\alpha^{N}\Phi({x}_{i},z)+\sum\limits_{m=1}^{M}\alpha_{m}^{N}\Phi({x}_{i},z_{m})=0\quad\text{for
all}\quad 1\leq i\leq N.$

We then conclude that (up to a subsequence) $\alpha_{m}^{N}$,
$\alpha^{N}\rightarrow\alpha_{m}$, $\alpha$ as $N\rightarrow\infty$. This
gives that due to the density of $S$ and since $\Phi({x},z)$ is analytic with
respect to $x$ for $z\neq x$, we obtain that

$\alpha\Phi({x},z)+\sum\limits_{m=1}^{M}\alpha_{m}\Phi({x},z_{m})=0\quad\text{for
all}\quad{x}\in C.$

Since the mapping

$x\,\longmapsto\,\alpha\,\Phi({x},z)+\sum\limits_{m=1}^{M}\alpha_{m}\Phi({x},z_{m})$

is a radiation exterior solution to the Helmholtz equation in
$\mathbb{R}^{d}\setminus\overline{\Omega}$ we conclude that

$\alpha\,\Phi({x},z)+\sum\limits_{m=1}^{M}\alpha_{m}\Phi({x},z_{m})=0\quad\text{for
all}\quad{x}\in\mathbb{R}^{d}\setminus\\{z\\}\cup\\{z_{m}\,:\,m=1,\dots,M\\}$

by uniqueness to the exterior Dirichlet problem and unique continuation. Now
letting $x\rightarrow z,z_{m}$ gives that $\alpha_{m}=\alpha=0$ for $1\leq
m\leq M$, but by (6) and the convergence of the coefficients $\alpha_{m}^{N}$,
$\alpha^{N}$ as $N\rightarrow\infty$ we have that

$|\alpha|+\sum\limits_{m=1}^{M}|\alpha_{m}|=1$

which gives a contradiction. Moreover, the fact that ${\bf U}$ has full rank
is a consequence of the given arguments. ∎

Now by combining this with Lemma 2.1 and Theorem 2.1 we have the following
range test which can be used to reconstruct the set
$z\in\\{z_{m}\,:\,m=1,\dots,M\\}.$

######  Theorem 2.2.

Assume that the set $S=\\{{x}_{i}\in C\,:\,i\in\mathbb{N}\\}$ is dense in $C$.
If $z\in\Omega$ then there is a number $N_{0}\in\mathbb{N}$ such that for all
$N\geq N_{0}$

${\bf\phi}_{z}\in\text{Range}\big{(}{\bf N}{\bf N}^{*}\big{)}\quad\text{ if
and only if }\quad z\in\\{z_{m}\,:\,m=1,\dots,M\\}.$

Now we are ready to construct an indicator function which characterizes the
locations of the small obstacles. To this end, let ${\bf
P}:\mathbb{C}^{N}\mapsto\text{Null}\big{(}{\bf N}{\bf N}^{*}\big{)}$ be the
orthogonal projection onto $\text{Null}\big{(}{\bf N}{\bf N}^{*}\big{)}$.
Therefore, by Theorem 2.2 we have that

${\bf P}{\bf\phi}_{z}=0\quad\text{ if and only if }\quad
z\in\\{z_{m}\,:\,m=1,\dots,M\\}.$

Notice that since ${\bf N}{\bf N}^{*}$ is a self-adjoint matrix we have that
it is orthogonally diagonalizable. So we let ${\bf w}_{j}$ be the $j$-th
orthonormal eigenvector of ${\bf N}{\bf N}^{*}$, we also let
$r$=Rank$\big{(}{\bf N}{\bf N}^{*}\big{)}$. This gives that ${\bf w}_{j}$ for
$r+1\leq j\leq N$ is an orthonormal basis for $\text{Null}\big{(}{\bf N}{\bf
N}^{*}\big{)}$, and therefore we have that ${\bf P}$ is given by

${\bf P}{\bf\phi}_{z}=\sum\limits_{j=r+1}^{N}\big{(}{\bf\phi}_{z},{\bf
w}_{j}\big{)}_{\ell^{2}}{\bf w}_{j}.$

This leads to the following result:

######  Lemma 2.2.

Let ${\bf w}_{j}$ be the $j$-th orthonormal eigenvector of ${\bf N}{\bf
N}^{*}$ and let $r$=Rank$\big{(}{\bf N}{\bf N}^{*}\big{)}$. Assume that the
set $S=\\{{x}_{i}\in C\,:\,i\in\mathbb{N}\\}$ is dense in $C$. If $z\in\Omega$
then there is a number $N_{0}\in\mathbb{N}$ such that for all $N\geq N_{0}$

$\mathcal{I}(z)=\left[\sum\limits_{j=r+1}^{N}\left|\big{(}{\bf\phi}_{z},{\bf
w}_{j}\big{)}\right|_{\ell^{2}}^{2}\right]^{-1}<\infty\quad\text{ if and only
if }\quad z\in\\{z_{m}\,:\,m=1,\dots,M\\}.$ (7)

###  2.3 Numerical validation of the MUSIC algorithm

We now present a few examples of reconstructing small obstacles using the
MUSIC algorithm developed in this section in ${\mathbb{R}}^{2}$. To simulate
the data, we use the Born approximation of the scattered field and a $2d$
Gaussian quadrature to approximate the volume integral. We now let the
collection curve $C$ be the boundary of a disk of radius one. In the following
calculations we use $32$ different sources and receivers at the points

$x_{i}=y_{i}=\big{(}\cos(2\pi(i-1)/32),\sin(2\pi(i-1)/32)\big{)}^{\top}\quad\text{for}\quad
i=1,\dots,32.$

This leads to an approximated multi-static response matrix

${\bf
N}\approx\big{[}u^{s}_{B}({x}_{i},y_{j})\big{]}_{i,j=1}^{32}\quad\textrm{ such
that }\quad
u^{s}_{B}({x},y)=k^{2}\sum_{p}\omega_{p}\big{(}n(z_{p})-1\big{)}\Phi(x,z_{p}){\Phi(z_{p},y)},$

where $\omega_{p}$ and $z_{p}$ are the weights and quadrature points for the
numerical approximation of the volume integral by a 2$d$ Gaussian quadrature
method. We give examples with random noise added to the simulated data for
$u^{s}_{B}({x}_{i},y_{j})$. The random noise level is given by $\delta$ where
the noise is added to the multi-static response matrix such that

$\Big{[}u^{s}_{B}({x}_{i},y_{j})\big{(}1+\delta
E_{i,j}\big{)}\Big{]}_{i,j=1}^{32}\quad\text{with the random matrix}\quad
E\quad\text{satisfying }\quad\|E\|_{2}=1.$

The set $\\{z_{m}\,:\,m=1,\dots,M\\}$ can be visualized by plotting the
corresponding indicator function given in (7)

$\mathcal{I}(z)=\left[\sum\limits_{j=r+1}^{32}\left|\big{(}{\bf\phi}_{z},{\bf
w}_{j}\big{)}\right|_{\ell^{2}}^{2}\right]^{-1}$

where ${\bf w}_{j}\in\mathbb{C}^{32}$ are the eigenvectors of the approximated
matrix ${\bf N}{\bf N}^{*}$ and $r$=Rank$\big{(}{\bf N}{\bf N}^{*}\big{)}$ is
computed by the Matlab command rank. In the following examples, we fix the
wave number $k=1$.

Figure 1:  Reconstruction of the location of the two scatterers, the circle of
radius 0.2 centered at $(-0.5,0.5)$ and the ellipse with $a=0.2$ and $b=0.1$
centered at $(0.5,-0.5)$. The refractive index $n=5$ in both scatterers.
Contour plot of the indicator function $z\mapsto\mathcal{I}(z)$ where the red
lines are the boundary of the scatterer.  Figure 2:  Reconstruction of the
location of one scatterer, the ellipse with $a=0.2$ and $b=0.1$ centered at
$(0.5,-0.5)$. The refractive index $n=2+\text{i}$ in the scatterer. Contour
plot of the indicator function $z\mapsto\mathcal{I}(z)$ where the red lines
are the boundary of the scatterer.  Figure 3:  Reconstruction of the location
of the a square scatterer, $[-0.2,0.2]^{2}$. The refractive index
$n(x_{1},x_{2})=x_{1}^{2}+2$ in the scatterer. Contour plot of the indicator
function $z\mapsto\mathcal{I}(z)$ where the red lines are the boundary of the
scatterer. Remark: Even though both figures look identical, different data we
used in the reconstructions.

###  2.4 Bayesian approach to approximating the Refractive Index

In this section, we consider the the second part of our inverse problem which
is to reconstruct (i.e. approximate) the refractive index of the scatterer(s).
In particular, the idea is to use a hierarchical Bayesian model for estimating
the refractive index by a constant in the scatterer(s) which has been
reconstructed, possibly by the MUSIC algorithm proposed in the previous
section. Since we assume that $n(x)$ is continuous in the scatterer the small
size of the obstacle implies that there will be only small deviations from the
value at the center of the obstacle.

Here we assume that we that an approximation of the scatterer is known using
the reconstruction from the previous section or by another qualitative method.
We use a Bayesian collocation method to approximate the value of the
refractive index at the center of the obstacle. For this discussion, assume
that the obstacle is centered at the point $x_{0}$. Rather than solving for
the value of the unknown $n(x_{0})$ exactly, we solve for a probability
distribution of the possible values of $n(x_{0})$. To do this, we must assume
a prior distribution for possible values of $n(x_{0})$ and a data model for
the scattered wave, which we will discuss below. Then, Markov chain Monte
Carlo methods allow us to sample from the posterior distribution of $n(x_{0})$
given the scattered wave and the prior distribution. The mode of the resulting
posterior distribution (called the Maximum a posteriori estimator or MAP). The
MAP in certain settings can be shown to be equivalent to the standard linear
regression estimator or the frequentist statistics maximum likelihood
estimator. In this context, our approach provides more information on possible
values of $n(x_{0})$, and because the scatterer is small with a continuous
refractive index.

To proceed, we define the necessary notation. Assume we have a reconstruction
of $D$ via the MUSIC algorithm, denoted $\widehat{D}$. Furthermore, we assume
that the measured scattered field denoted $U^{s}_{B}$ is given by the Born
approximation on $D$, that is

$U^{s}_{B}(x,y)=k^{2}\int_{\widehat{D}}\big{(}n(z)-1\big{)}\Phi(x,z){\Phi(z,y)}\,dz+\mathcal{O}(\delta)$

where $\delta$ is the amount of error added for the reconstruction of $D$ and
standard measurement error. We assume the error term is modeled by a normal
distribution with zero mean and standard deviation $\delta$, then we may write
the data model, as a statistical distribution where values of $U^{s}_{B}$ are
drawn from

$U^{s}_{B}(x,y)\sim
k^{2}\int_{\widehat{D}}\big{(}n(z)-1\big{)}\Phi(x,z){\Phi(z,y)}\,dz+\mathcal{N}(0,\delta).$
(8)

Now, recall that we have readings of the scattered field $U^{s}_{B}$ at points
$x_{i}\in C$ with the location of the point source given by $y_{i}\in C$ with
$C$ being the collection curve discussed in the previous section. This forms
the data set
$R=\big{\\{}(x_{i},y_{i},U^{s}_{B}(x_{i},y_{i}))\big{\\}}_{i=1}^{N}$ where $N$
indicates the number of sources and receivers.

We now motivate the subsequent Bayesian model. Since the goal to reconstruct
the value of the refractive index at the center of the scatterer, using
Taylor’s Theorem we have that $n(x)=n(x_{0})+\mathcal{O}(h)$ where
$h=|\widehat{D}|$ is small. Therefore, the constant term $n(x_{0})$ of the
Taylor polynomial centered at the center of $\widehat{D}$ is approximately
close to all values of $n(x)$ on $\widehat{D}$. Based on these assumptions, we
define the following auxiliary function

$\eta(x)=n(x)-1=n(x_{0})-1+\mathcal{O}(h).$ (9)

Therefore, we describe the following Bayesian model to solve for the posterior
distribution of $n(x_{0})$ given readings of the scatterer $R$, i.e. we solve
for the conditional pdf function $P(\gamma|R,w_{p},z_{p})$

$U^{s}_{B}(x,y)\sim\mathcal{N}(\mu,\delta^{2})\quad\text{with}\quad\mu(x,y)=k^{2}\sum_{p}\omega_{p}\eta(z_{p})\Phi(x,z_{p})\Phi(z_{p},y)$

where $\omega_{p}$ and $z_{p}$ are the weights and quadrature points for the
numerical approximation of the volume integral over $\widehat{D}$. Here we
assume that the error from the numerical integration is negligible since one
can take a quadrature rule to approximate the integral to some preferred
tolerance. By the definition of $\eta$ we have that

$\eta(z)\sim\mathcal{N}(\gamma,h^{2})\quad\text{ where
}\quad\gamma\sim\mathcal{U}(-\infty,\infty).$

We may solve the model for posterior distribution for $\gamma$, which is given
according to the proportion:

$P\big{(}\gamma,\eta|R,\omega_{p},z_{p}\big{)}\propto\ell\big{(}\mu|R,\gamma,\eta,\omega_{p},z_{p}\big{)}P(\eta|\gamma,\omega_{p},z_{p})P(\gamma),$
(10)

where $\ell$ is the likelihood function arising from the data model for
$y_{i}$ and $P(\nu|\gamma),P(\nu)$ are the assumed priors given above. One may
marginalize over $\eta$ obtain $P\big{(}\gamma|R,\omega_{p},z_{p}\big{)}$ if
desired. By using the Metropolis-Hastings algorithm, we can derive this
posterior distribution, and thus values for $\gamma$. We note that this model
is hierarchical with care taken to model the sources of error. This should
lead to a more realistic solution.

For our calculation, we use the implementation in the Python package PyMC3
[32]. Therefore we will have $\gamma$ which is the average value of $\eta$
which will be the approximation on $n(x_{0})-1$. For our numerical example, it
is a common practice to use a normal distribution with a large standard
deviation ($10^{5}$) rather than a flat prior, $\mathcal{U}(-\infty,\infty)$.
For our example we consider the square scatterer $[-0.2,0.2]^{2}$ with
refractive index $n(x_{1},x_{2})=x^{2}_{1}+2$ where we add $15\%$ noise to the
measurements. For the reconstruction we use the domain actual scatterer
$D=[-0.2,0.2]^{2}$ as well as $\widehat{D}=[-0.265,0.265]^{2}$ which is an
approximation taken from the MUSIC Algorithm’s reconstruction of the obstacle.
Using the methods described in this section we are able to determine the
distribution of possible values of $\gamma$ which should approximately be
$n(x_{0})-1$.

Figure 4: Recovery from measurements in the case when
$n(x_{1},x_{2})=x^{2}_{1}+2$. The distribution is the possible values of
$\gamma\approx n(x_{0})-1=1$ in the reconstructed scatterer
$D=[-0.2,0.2]^{2}$. Figure 5: Recovery from measurements in the case when
$n(x_{1},x_{2})=x^{2}_{1}+2$. The distribution is the possible values of
$\gamma\approx n(x_{0})-1=1$ in the reconstructed scatterer
$\widehat{D}=[-0.265,0.265]^{2}$.

Similar approaches utilizing Bayesian statistics and the Born approximation
have been used when applying compressive sensing techniques to inverse
scattering [33],[34],[35],[36]. In those context, there is an a-priori
assumption on the sparsity of the scatterer. The approach has also been
applied to continuous random media [37]. However, in this case, the scatterer
is not sparse and it is not a continuous random medium. Moreover, we extend
the approach by using the MUSIC algorithm to find the support of the
scatterer, and use that information in our recovery of the refractive index.

##  3  Sampling methods for anisotropic scatterers

###  3.1 Scattering by extended anisotropic scatterers

In this section we formulate the direct and inverse time-harmonic scattering
problems under consideration. To this end, we let the point source located at
the point $y\in C$ be given by the fundamental solution to Helmholtz equation
$\Phi(x,y)$ where $C$ is assumed to be a closed curve (in $\mathbb{R}^{2}$) or
surface (in $\mathbb{R}^{3}$) that is class $\mathcal{C}^{2}$ or sufficiently
smooth, such that $H^{1/2}(C)$ is compactly embedded in $L^{2}(C)$. We
consider the scattering of the (non-standard) incident field
$u^{i}(\cdot\,,y)=\overline{\Phi(\cdot\,,y)}$ for $y\in C$ by a penetrable
anisotropic obstacle. Assume that the obstacle $D\subset\mathbb{R}^{d}$ (for
$d=2,3$) is a bounded simply connected open set having piece-wise smooth
boundary $\partial D$ with $\nu$ being the unit outward normal to the
boundary.

We assume that the constitutive parameters of the scatterer $D$ are
represented by a symmetric matrix
$\tilde{A}\in\mathcal{C}^{2}\left(D,\mathbb{C}^{d\times d}\right)$ and a
scalar function $\tilde{n}\in\mathcal{C}(D)$ such that

$\overline{\xi}\cdot\Re\left(\tilde{A}(x)\right)\xi\geq
a_{min}|\xi|^{2}>0\quad\text{and}\quad\overline{\xi}\cdot\Im\left(\tilde{A}(x)\right)\xi\leq
0,$

where as

$\Re\big{(}\tilde{n}(x)\big{)}\geq
n_{min}>0\quad\text{and}\quad\Im\big{(}\tilde{n}(x)\big{)}\geq 0$

for almost every $x\in D$ and all $\xi\in{\mathbb{C}}^{d}$. Outside the
obstacle $D$, the material parameters are homogeneous isotropic with
refractive index scaled to one. We denote by $A$ and $n$ the constitutive
parameters of the homogenous background with the anisotropic obstacle $D$ in
${\mathbb{R}}^{d}$ by

$A(x)=\tilde{A}(x)\chi_{D}+I(1-\chi_{D})\quad\text{and}\quad
n(x)=\tilde{n}(x)\chi_{D}+1(1-\chi_{D})$

where $I$ denotes the $d\times d$ identity matrix and $\chi_{D}$ is the
characteristic function on $D$. Note that the support of the contrasts $A-I$
and $n-1$ is $\overline{D}$. The radiating scattered field $u^{s}(\cdot\,,y)$
given by the (non-standard) point source $u^{i}(\cdot\,,y)$ satisfies

$\displaystyle\nabla_{x}\cdot
A\nabla_{x}u^{s}+k^{2}nu^{s}=\nabla_{x}\cdot(I-A)\nabla_{x}u^{i}+k^{2}(1-n)u^{i}\quad\textrm{
in }\,\mathbb{R}^{d}$ (11)
$\displaystyle{\partial_{r}u^{s}}-iku^{s}=\mathcal{O}\left(\frac{1}{r^{(d+1)/2}}\right)\quad\text{
as }\quad r\rightarrow\infty,$ (12)

where $r=|x|$ and (12) is satisfied uniformly in all directions
$\hat{x}=x/|x|$. It can be shown that the scattering problem (11)-(12) has a
unique solution $u^{s}\in H^{1}_{loc}(\mathbb{R}^{d})$ by using a variational
approach (see for e.g. [8]). The total field corresponding to the scattering
problem (11) is defined as $u(\cdot\,,y)=u^{s}(\cdot\,,y)+u^{i}(\cdot\,,y)$
and satisfies

$u^{(-)}=u^{(+)}\quad\textrm{ and }\quad\frac{\partial
u^{(-)}}{\partial\nu_{A}}=\frac{\partial u^{(+)}}{\partial\nu}\quad\textrm{ on
}\partial D,$

where the superscripts $+$ and $-$ for a generic function indicates the trace
on the boundary taken from the exterior or interior of the domain,
respectively. Now assume that we have the data set of near field measurements
$u^{s}(x,y)$ for all $x,y\in C$ and define the Near Field operator

${N}:L^{2}(C)\longmapsto L^{2}(C)\quad\text{given
by}\quad({N}g)(x)=\int_{C}u^{s}(x,y)g(y)\,ds_{y}.$

The inverse problem we consider is to reconstruct the scattering obstacle $D$
from a knowledge of $u^{s}(x,y)$ for all $x,\,y\in C$. We will let $\Omega$ be
the bounded simply connected open set such that $\overline{D}\subset\Omega$
and $C=\partial\Omega$.

In the coming sections, we will see that to obtain a symmetric factorization
of the near field operator we require that the incident wave be the complex
conjugate of a point source. However, it is known that these non-physical
sources can be approximated by a linear combination of physical point sources
(see [23]). In [22], the authors are able to avoid using non-physical sources
but must construct the so-called outgoing-to-incoming operator. If $C$ is
given by the boundary of a disk/sphere centered at the origin, then the
outgoing-to-incoming operator is given by a boundary integral operator that
does not depend on the type of scatterer.

###  3.2 Factorization of the near field operator

This section is dedicated to constructing a suitable factorization of the near
field operator $N$ so that we are able to use the theoretical results in [24]
and/or [27], in order to derive indicator functions for the support of $D$ in
terms of the near field measurements. Here we consider a modified linear
sampling method as well as the factorization method and just as in [4] we will
us the range characterization from the factorization method to show that these
two methods give equivalent indicator functions. To this end, motivated by the
source problem given in equation (11)-(12) for the scattered field due to the
(non-standard) point source, we consider the problem of finding $u\in
H^{1}_{loc}(\mathbb{R}^{d})$ for a given $v\in H^{1}(D)$ such that

$\displaystyle\nabla\cdot A\nabla u+k^{2}nu=\nabla\cdot(I-A)\nabla
v+k^{2}(1-n)v\quad\textrm{ in }\,\mathbb{R}^{d}$ (13) $\displaystyle\hskip
28.45274pt{\partial_{r}u^{s}}-iku^{s}=\mathcal{O}\left(\frac{1}{r^{(d+1)/2}}\right)\quad\text{
as }\quad r\rightarrow\infty.$

It can be shown that for any given $v$, the solution operator that maps
$v\mapsto u$ is a bounded linear mapping from $H^{1}(D)$ to
$H^{1}_{loc}(\mathbb{R}^{d})$. At this point, let us recall the exterior
Dirichlet-to-Neumann mapping for Helmholtz equation

$\mathbb{T}_{k}:H^{1/2}(\partial B_{R})\longmapsto H^{-1/2}(\partial
B_{R})\quad\text{given
by}\quad\mathbb{T}_{k}f=\frac{\partial\varphi}{\partial\nu}\quad\text{on}\,\,\,{\partial
B_{R}},$

where $\varphi\in H^{1}_{loc}(\mathbb{R}^{d}\setminus\overline{B_{R}})$
satisfies

$\displaystyle\Delta\varphi+k^{2}\varphi=0\quad\quad\text{in
}\mathbb{R}^{d}\setminus\overline{B}_{R}\quad\text{and}\quad\varphi=f\quad\text{on
}\partial B_{R}$

along with the radiation condition (12) where
$B_{R}=\\{x\in\mathbb{R}^{d}:|x|<R\\}$. With the help of the Dirichlet-to-
Neumann operator we can write (13) in its equivalent variational form: find
$u\in H^{1}(B_{R})$ such that

$\displaystyle-\int\limits_{B_{R}}A\nabla
u\cdot\nabla\overline{\varphi}-k^{2}nu\overline{\varphi}\,dx\,+\int\limits_{\partial
B_{R}}\overline{\varphi}\,\mathbb{T}_{k}u\,ds$ $\displaystyle\hskip
85.35826pt=-\int\limits_{D}(I-A)\nabla
v\cdot\nabla\overline{\varphi}-k^{2}(1-n)v\overline{\varphi}\,dx,\quad\forall\varphi\in
H^{1}(B_{R}),$ (14)

where the boundary integral over $\partial B_{R}$ is understood as the dual
pairing between $H^{1/2}(\partial B_{R})$ and $H^{-1/2}(\partial B_{R})$. It
is clear by well-posedness that if we take
$v=\overline{\Phi(\cdot\,,y)}\big{|}_{D}$ then we have that the scattered
field $u^{s}(\cdot\,y)=u$ for all $x\in\mathbb{R}^{d}$. Now define the source
to trace operator

$G:H^{1}(D)\longmapsto L^{2}(C)\quad\text{given by}\quad Gv=u\big{|}_{C}$

as well as the bounded linear operator

$H:L^{2}(C)\longmapsto H^{1}(D)\quad\text{given
by}\quad(Hg)(x)=\int_{C}\overline{\Phi(\cdot\,,y)}g(y)\,ds_{y}\big{|}_{D}.$

By superposition, we have that near field operator is given by $N=GH.$ In
order to use the theory developed in [25] and/or [27], we must construct a
symmetric factorization of the near field operator. Therefore, we compute the
adjoint of the operator $H$ which is given in the following result.

######  Theorem 3.1.

The operator $H^{*}:H^{1}(D)\longmapsto L^{2}(C)$ is given by

$H^{*}v=\tilde{v}\big{|}_{C}\quad\text{for all}\quad v\in H^{1}(D),$

where $\tilde{v}\in H^{1}_{loc}(\mathbb{R}^{d})$ is the unique radiating
solution to

$\displaystyle-\int\limits_{B_{R}}\nabla\tilde{v}\cdot\nabla\overline{\varphi}-k^{2}\tilde{v}\overline{\varphi}\,dx+\int\limits_{\partial
B_{R}}\overline{\varphi}\mathbb{T}_{k}\tilde{v}\,ds=(v,\varphi)_{H^{1}(D)}$
(15)

for all $\varphi\in H^{1}(B_{R})$.

###### Proof.

Let the radius of $B_{R}$ be sufficiently large so that
$\overline{\Omega}\subset B_{R}$. Given any $v\in H^{1}(D)$ we can construct a
unique radiating field $\tilde{v}\in H^{1}_{loc}(\mathbb{R}^{d})$ that
satisfies (15). We now let $w\in H^{1}_{loc}(\mathbb{R}^{d})$ be defined as

$w=\int_{C}\overline{\Phi(x\,,y)}g(y)\,ds_{y}\quad\text{for}\,\,\,\,x\in\mathbb{R}^{d}.$

Therefore, we have that

$\displaystyle(H^{*}v,g)_{L^{2}(C)}=(v,Hg)_{H^{1}(D)}$ $\displaystyle\hskip
42.67912pt=-\int\limits_{B_{R}}\nabla\tilde{v}\cdot\nabla\overline{w}-k^{2}\tilde{v}\overline{w}\,dx+\int\limits_{\partial
B_{R}}\overline{w}\mathbb{T}_{k}\tilde{v}\,ds$ $\displaystyle\hskip
42.67912pt=-\int\limits_{B_{R}\setminus\overline{\Omega}}\nabla\tilde{v}\cdot\nabla\overline{w}-k^{2}\tilde{v}\overline{w}\,dx-\int\limits_{{\Omega}}\nabla\tilde{v}\cdot\nabla\overline{w}-k^{2}\tilde{v}\overline{w}\,dx+\int\limits_{\partial
B_{R}}\overline{w}\mathbb{T}_{k}\tilde{v}\,ds.$

Using integration by parts on the volume integrals and the fact that $w$
solves Helmholtz equation in $\mathbb{R}^{d}\setminus C$ we have that

$\displaystyle(H^{*}v,g)_{L^{2}(C)}=\int\limits_{C}\tilde{v}\left(\frac{\partial}{\partial\nu}\overline{w}^{(-)}-\frac{\partial}{\partial\nu}\overline{w}^{(+)}\right)ds+\int\limits_{\partial
B_{R}}\overline{w}\frac{\partial}{\partial\nu}\tilde{v}-\tilde{v}\frac{\partial}{\partial\nu}\overline{w}\,ds.$

Now using the jump relations for the normal derivative of the single layer
potential for Helmholtz equation gives that

$\frac{\partial}{\partial\nu}\overline{w}^{(-)}-\frac{\partial}{\partial\nu}\overline{w}^{(+)}=\overline{g}\quad\text{on}\quad
C.$

Notice that both $\tilde{v}$ and $\overline{w}$ satisfy the radiation
condition (12) therefore letting $R\rightarrow\infty$ gives that

$(H^{*}v,g)_{L^{2}(C)}=(\tilde{v},g)_{L^{2}(C)},$

proving the claim. ∎

Just as in [12] notice that (13), implies that

$\displaystyle\Delta
u+k^{2}u=\nabla\cdot(I-A)\nabla(v+u)+k^{2}(1-n)(v+u)\,\textrm{ in
}\,\mathbb{R}^{d},$ (16)

where $u\in H^{1}_{loc}(\mathbb{R}^{d})$ is the solution to (13) for a given
$v\in H^{1}(D)$. The variational form of (16) is given by

$\displaystyle-\int\limits_{B_{R}}\nabla
u\cdot\nabla\overline{\varphi}-k^{2}u\overline{\varphi}\,dx+\int\limits_{\partial
B_{R}}\overline{\varphi}\mathbb{T}_{k}u\,ds=-\int\limits_{D}(I-A)\nabla(v+u)\cdot\nabla\overline{\varphi}\,dx$
$\displaystyle\hskip
56.9055pt+k^{2}\int\limits_{D}(1-n)(v+u)\overline{\varphi}\,dx\qquad\forall\varphi\in
H^{1}(B_{R}).$

Now appealing to the Riesz representation theorem, we can define the bounded
linear operator $T:H^{1}(D)\longmapsto H^{1}(D)$ such that for all $v\in
H^{1}(D)$

$(Tv,\varphi)_{H^{1}(D)}=-\int\limits_{D}(I-A)\nabla(v+u)\cdot\nabla\overline{\varphi}-k^{2}(1-n)(v+u)\overline{\varphi}\,dx.\,\,\,$
(17)

Notice for $u\in H^{1}_{loc}(\mathbb{R}^{d})$ the unique solution to (13) that
by the definition of $G$ we have that $u\big{|}_{C}=Gv$. Now by the definition
of $H^{*}$ and $T$ we conclude that $u\big{|}_{C}=H^{*}Tv$. We conclude that
$Gv=H^{*}Tv$ for all $v\in H^{1}(D)$. Now recalling that $N=GH$ gives the
desired factorization.

######  Lemma 3.1.

The near field operator $N:L^{2}(C)\longmapsto L^{2}(C)$ for the scattering
problem (11)-(12) has the factorization $N=H^{*}TH$.

###  3.3 The solution to the inverse problem

The main goal of this section is to connect the support of $D$ to the range of
an operator defined by the measured near field operator. We make this
connection by analyzing the factorization of the near field operator developed
in the previous section. We recall from the previous section that we have the
following factorization $N=H^{*}TH$. Under some appropriate assumptions on $H$
and $T$, the factorization method gives that
Range$(H^{*})$=$\text{Range}(N_{\sharp}^{1/2})$. Here $N_{\sharp}$ will be
given by $N_{\sharp}=|\Re({N})|+|\Im({N})|$ for a non-absorbing material or
$N_{\sharp}=\sigma\Im({N})$ with $|\sigma|=1$ for an absorbing material where

$\Re(N)=\frac{N+N^{*}}{2}\quad\text{and}\quad\Im(N)=\frac{N-N^{*}}{2\text{i}}.$

Furthermore for a generic self-adjoint compact operator $B$ on a Hilbert space
$U$, by appealing to the Hilbert-Schmidt theorem we can define $|B|$ in terms
of the spectral decomposition such that

$|B|x=\sum|\lambda_{j}|(x,\psi_{j})\psi_{j}$

for all $x\in U$ where $(\lambda_{j},\psi_{j})\in\mathbb{R}\times U$ is the
orthonormal eigensystem of $B$.

####  3.3.1 Analysis of the operator $H$

Now we turn our attention to showing that the operator $H$ has the properties
needed to apply Theorems A.1. To this end, let us define the interior
transmission problem as finding a pair $(w,v)\in H^{1}(D)\times H^{1}(D)$ such
that for a given pair $(f,h)\in H^{1/2}(\partial D)\times H^{-1/2}(\partial
D)$ satisfies

$\displaystyle\nabla\cdot A\nabla w+k^{2}nw=0\quad\text{and}\quad\Delta
v+k^{2}v=0\quad\text{ in }\,\,D$ (18) $\displaystyle
w-v=f\quad\text{and}\quad\frac{\partial w}{\partial\nu_{A}}-\frac{\partial
v}{\partial\nu}=h\quad\text{ on }\,\,\partial D$ (19)

######  Definition 3.1.

The values of $k\in\mathbb{C}$ for which the homogeneous interior transmission
problem$\big{(}$i.e. (18)-(19) with $(f,h)=(0,0)\big{)}$, has nontrivial
solutions are called transmission eigenvalues for $A,\,n$ and $D$.

It is known that under appropriate assumptions on the coefficients that the
interior transmission problem (18)-(19) is Fredholm of index zero, i.e if $k$
is not a transmission eigenvalue then (18)-(19) is well-posed (see for e.g.
[7] and [11]). The following results are know for the transmission eigenvalue
problem (see [7], [8] and [13])

  1. 1.

If the scatterer is absorbing such that $\Im(A)\leq 0$ and $\Im(n)>0$ in $D$
then there are no real transmission eigenvalues.

  2. 2.

If the scatterer is non-absorbing (i.e. $\Im(A)=0$ and $\Im(n)=0$) then the
set of real transmission eigenvalues is at most discrete provided that $A-I$
is positive (or negative) definite uniformly in $D$.

  3. 3.

If the scatterer is non-absorbing then the set of real transmission
eigenvalues is at most discrete provided that $A-I$ is positive (or negative)
definite uniformly and $n-1$ is positive (or negative) in a neighborhood of
the boundary $\partial D$.

######  Assumption 3.1.

The wave number $k\in\mathbb{R}$ is not a transmission eigenvalue for
(18)-(19).

We are now ready to connect the support of $D$ to the range of $H^{*}$. Since
we assume that the obstacle $D$ is contained in the interior of $\Omega$ we
only need to consider sampling points $z\in\Omega$.

######  Theorem 3.2.

Under Assumption 3.1, we have that for $z\in\Omega$

$\Phi(\cdot\,,z)\in\text{Range}(H^{*})\quad\text{ if and only if}\quad z\in
D.$

###### Proof.

We start with the case where $z\in\Omega\setminus\overline{D}$ and assume that
$v\in H^{1}(D)$ is such that $H^{*}v=\Phi(\cdot\,,z)$. By the definition of
$H^{*}$ there exists a $\tilde{v}\in H^{1}_{loc}(\mathbb{R}^{d})$ satisfying
(15) which implies that

$\Delta\tilde{v}+k^{2}\tilde{v}=0\quad\text{in}\quad\mathbb{R}^{d}\setminus\overline{D}\quad\text{and}\quad\tilde{v}=\Phi(\cdot,z)\quad\text{on}\,\,C$

along with the Sommerfeld radiation condition (12). Now since
$z\in\Omega\setminus\overline{D}$, we have that $\Phi(\cdot\,,z)$ is in
$H^{1}_{loc}(\mathbb{R}^{d}\setminus\overline{\Omega})$, and by the uniqueness
of solutions to the exterior Dirichlet problem for Helmholtz equation, we have
that $\tilde{v}=\Phi(\cdot,z)$ in $\mathbb{R}^{d}\setminus\overline{\Omega}$.
Therefore by Holmgren’s theorem and unique continuation, we have that
$\tilde{v}=\Phi(\cdot,z)$ in $\Omega\setminus(\overline{D}\cup\\{z\\})$. We
then conclude that $\Phi(\cdot\,,z)\notin\text{Range}(H^{*})$ since

$\Phi(\cdot,z)\notin
H^{1}\big{(}B_{r}(z)\big{)}\quad\text{and}\quad\tilde{v}\in
H^{1}\big{(}B_{r}(z)\big{)},$

for any disk $B_{r}(z)$ centered at $z$ of radius $r>0$.

Now assume that $z\in D$ then we have that $\Phi(\cdot\,,z)\in
H^{1}_{loc}(\mathbb{R}^{d}\setminus\overline{D})$. Since $k$ is not a
transmission eigenvalue in $D$ we have that there is a unique solution to the
interior transmission problem (18)-(19) with

$(f,h)=\left(\Phi(\cdot\,,z)\,,\,\displaystyle{\frac{\partial}{\partial\nu}}\Phi(\cdot\,,z)\right)\quad\text{denoted}\quad(w_{z},v_{z})\in
H^{1}(D)\times H^{1}(D).$

Now define the function

$u_{z}:=\left\\{\begin{array}[]{rcl}w_{z}-v_{z}&\mbox{in}&D\\\ \\\
\Phi(\cdot\,,z)&\mbox{in}&\mathbb{R}^{d}\setminus\overline{D}.\end{array}\right.$

Therefore we have that $u_{z}\in H^{1}_{loc}(\mathbb{R}^{d})$ with
$u_{z}\big{|}_{C}=\Phi(\cdot\,,z)\big{|}_{C}$ such that

$\displaystyle\Delta u_{z}+k^{2}u_{z}=\nabla\cdot(I-A)\nabla
w_{z}+k^{2}(1-n)w_{z}\,\textrm{ in }\,\mathbb{R}^{d}.$

The variational form of the above equation gives that for all $\varphi\in
H^{1}(B_{R})$ with $\mathbb{T}_{k}$ being the Dirichlet to Neumann mapping

$\displaystyle-\int\limits_{B_{R}}\nabla
u_{z}\cdot\nabla\overline{\varphi}-k^{2}u_{z}\overline{\varphi}\,dx+\int\limits_{\partial
B_{R}}\overline{\varphi}\mathbb{T}_{k}u_{z}\,ds=$ $\displaystyle\hskip
56.9055pt-\int\limits_{D}(I-A)\nabla
w_{z}\cdot\nabla\overline{\varphi}-k^{2}(1-n)w_{z}\overline{\varphi}\,dx.$
(20)

Now let $\phi_{z}\in H^{1}(D)$ be defined from the right hand side of (20) by
appealing to the Riesz representation theorem, hence we have

$-\int\limits_{B_{R}}\nabla
u_{z}\cdot\nabla\overline{\varphi}-k^{2}u_{z}\overline{\varphi}\,dx+\int\limits_{\partial
B_{R}}\overline{\varphi}\mathbb{T}_{k}u_{z}\,ds=(\phi_{z},\varphi)_{H^{1}(D)}.$

Thus we now conclude that $H^{*}\phi_{z}=\Phi(\cdot,z)$ by the definition of
$H^{*}$ giving the result. ∎

Next we show that $H$ satisfies that necessary properties to apply Theorems
A.1.

######  Theorem 3.3.

The operator $H:L^{2}(C)\longmapsto H^{1}(D)$ is compact and injective.

###### Proof.

Notice that the mapping $v\in H^{1}(D)\mapsto\tilde{v}\in
H^{1}_{loc}(\mathbb{R}^{d})$ is bounded and by the Trace theorem the mapping
$\tilde{v}\in H^{1}_{loc}(\mathbb{R}^{d})\mapsto\tilde{v}\big{|}_{C}\in
H^{1/2}(C)$ is also bounded. Now using that $H^{1/2}(C)$ is compactly embedded
in $L^{2}(C)$ gives that $H^{*}$ is compact and therefore so is $H$.

Now assume that $g\in$ Null$(H)$, therefore we define

$w=\int_{C}{\Phi(x\,,y)}\overline{g(y)}\,ds_{y}\quad\text{for}\,\,\,\,x\in\mathbb{R}^{d}.$

We have that $w$ is a radiating solution to Helmholtz equation in
$\mathbb{R}^{d}\setminus C$, and by the definition of $H$ we have that $w=0$
in $D$, and therefore by unique continuation we have that $w=0$ in $\Omega$.
The continuity of the trace of $w$ on $C$ as well as the uniqueness of
solutions to the exterior Dirichlet problem for Helmholtz equation implies
that $w=0$ in $B_{R}$ for any radius sufficiently large. Now assume that
$\overline{\Omega}\subset B_{R}$ using the jump relations for the normal
derivative for the single layer potential for Helmholtz equation on $C$ gives
that

$0=\left(\frac{\partial}{\partial\nu}{w}^{(-)}-\frac{\partial}{\partial\nu}{w}^{(+)}\right)=\overline{g}$

proving the claim. ∎

####  3.3.2 Analysis of the operator $T$

Next we analyze the properties of the middle operator $T$ defined
variationally by equation (17). The operator $T$ is defined similarly to the
operator used in the factorization of the far field operator discussed in
[12], this gives that the analysis to show that $T$ satisfies the assumptions
in Theorem A.1 provided that the material is non absorbing (see [12] for
details). The following result can be obtained from the analysis in [12].

######  Theorem 3.4.

Assume $\Im(A)=0$, $\Im(n)=0$ in $D$ and $A-I$ is positive (or negative)
definite uniformly in $D$ and $n-1$ is negative (or positive) uniformly in
$D$. Then the operator $T:H^{1}(D)\longmapsto H^{1}(D)$ satisfies that
following:

  1. 1.

$T$ is injective

  2. 2.

$\Im(T)$ is compact and non-positive

  3. 3.

If $A-I$ is uniformly positive definite in $D$ then $\Re(T)$ is the sum of a
compact operator and a self-adjoint coercive operator

  4. 4.

If $A-I$ uniformly negative definite in $D$ then $-\Re(T)$ is the sum of a
compact operator and a self-adjoint coercive operator

The above result gives that the operator $T$ for a non absorbing media
satisfies the assumptions of Theorem A.1. We now show that for an absorbing
media that the imaginary part of $T$ is coercive. Recall that if $u\in
H^{1}_{loc}(\mathbb{R}^{d})$ is a radiating solution to the exterior Helmholtz
equation then we have that $u$ satisfies that following expansion

$\displaystyle
u(x)=\frac{\text{e}^{\text{i}k|x|}}{|x|^{\frac{d-1}{2}}}\left\\{u^{\infty}(\hat{x})+\mathcal{O}\left(\frac{1}{|x|}\right)\right\\}\;\textrm{
as }\;|x|\to\infty,$ (21)

where $\hat{x}:=x/|x|$ and $u^{\infty}(\hat{x})$ is the corresponding far
field pattern (see [15] in ${\mathbb{R}}^{3}$ and [8] in ${\mathbb{R}}^{2}$).

######  Theorem 3.5.

Assume that $\Im(A)$ is uniformly negative definite and $\Im(n)$ is uniformly
positive in $D$. Then the operator $-\Im(T):H^{1}(D)\longmapsto H^{1}(D)$ is
coercive.

###### Proof.

By (13) we have that for any $v_{j}\in H^{1}(D)$ there is a unique solution
$u_{j}\in H^{1}_{loc}(\mathbb{R}^{d})$. Now letting $\phi_{j}=v_{j}+u_{j}$ and
using (17) we have that

$\displaystyle(Tv_{1},v_{2})_{H^{1}(D)}=-\int\limits_{D}(I-A)\nabla\phi_{1}\cdot\nabla\overline{(\phi_{2}-u_{2})}-k^{2}(1-n)\phi_{1}\overline{(\phi_{2}-u_{2})}\,dx$
$\displaystyle\hskip
28.45274pt=-\int\limits_{D}(I-A)\nabla\phi_{1}\cdot\nabla\overline{\phi_{2}}-k^{2}(1-n)\phi_{1}\overline{\phi_{2}}\,dx$
$\displaystyle\hskip
85.35826pt+\int\limits_{D}(I-A)\nabla\phi_{1}\cdot\nabla\overline{u_{2}}-k^{2}(1-n)\phi_{1}\overline{u_{2}}\,dx.$

Now by equation (13), we have that

$\displaystyle\Delta
u_{1}+k^{2}u_{1}=\nabla\cdot(I-A)\nabla\phi_{1}+k^{2}(1-n)\phi_{1}\quad\text{
in }\quad\mathbb{R}^{d},$

and therefore multiplying by $\overline{u_{2}}$ and integrating by parts over
$B_{R}$ such that $D\subset B_{R}$ gives that

$\displaystyle-\int\limits_{B_{R}}A\nabla
u_{1}\cdot\nabla\overline{u_{2}}-k^{2}nu_{1}\overline{u_{2}}\,dx+\int\limits_{\partial
B_{R}}\overline{u_{2}}\frac{\partial u_{1}}{\partial\nu}\,ds$
$\displaystyle\hskip
142.26378pt=-\int\limits_{D}(I-A)\nabla\phi_{1}\cdot\nabla\overline{u_{2}}-k^{2}\int\limits_{D}(1-n)\phi_{1}\overline{u_{2}}\,dx.$

This gives that

$\displaystyle(Tv_{1},v_{2})_{H^{1}(D)}=-\int\limits_{D}(I-A)\nabla\phi_{1}\cdot\nabla\overline{\phi_{2}}-k^{2}(1-n)\phi_{1}\overline{\phi_{2}}\,dx$
$\displaystyle\hskip 85.35826pt+\int\limits_{B_{R}}\nabla
u_{1}\cdot\nabla\overline{u_{2}}-k^{2}u_{1}\overline{u_{2}}\,dx-\int\limits_{\partial
B_{R}}\overline{u_{2}}\frac{\partial u_{1}}{\partial\nu}\,ds.$ (22)

Now take the imaginary part of equation (22) with $u_{2}=u_{1}$ and
$\phi_{2}=\phi_{1}$. Furthermore by letting $R\rightarrow\infty$, we obtain
that

$\big{(}\Im(T)v_{1},v_{1}\big{)}_{H^{1}(D)}=\int\limits_{D}\Im(A)\nabla\phi_{1}\cdot\nabla\overline{\phi_{1}}-k^{2}\Im(n)|\phi_{1}|^{2}\,dx-k\int\limits_{\mathbb{S}}|u^{\infty}_{1}|^{2}\,ds(\hat{x})$
(23)

where the far field pattern $u_{1}^{\infty}$ is defined by the asymptotic
expansion in (21) and $\mathbb{S}=\\{{x}\in\mathbb{R}^{d}\,:\,|{x}|=1\\}$ is
the unit circle or sphere.

We are going to prove coercivity by using a contradiction argument, therefore
assume that $\exists\,v_{\ell}\in H^{1}(D)$ such that

$\|v_{\ell}\|_{H^{1}(D)}=1\quad\text{where}\quad\big{|}(\Im(T)v_{\ell},v_{\ell})_{H^{1}(D)}\big{|}\longrightarrow
0\quad\text{as}\quad\ell\rightarrow\infty.$

This gives that $u_{\ell}$ defined by solving (13) with source $v=v_{\ell}$ is
bounded in $H^{1}_{loc}(\mathbb{R}^{m})$. Using equation (23) as well as the
fact that $\Im(A)$ is uniformly negative definite and $\Im(n)$ is uniformly
positive in $D$ we conclude that $\phi_{\ell}\rightarrow 0$ in $H^{1}(D)$. By
the well-posedness of (16) we have that

$\|u_{\ell}\|_{H^{1}(B_{R})}\leq C\|\phi_{\ell}\|_{H^{1}(D)}\longrightarrow
0\quad\text{as}\quad\ell\rightarrow\infty.$

Now since $v_{\ell}=\phi_{\ell}-u_{\ell}$ we have that $v_{\ell}\rightarrow 0$
in $H^{1}(D)$ which gives that result by contradiction since
$||v_{\ell}||_{H^{1}(D)}=1$. ∎

######  Corollary 3.1.

  1. 1.

Assume $\Im(A)=0$, $\Im(n)=0$ in $D$ and $A-I$ is positive (or negative)
definite uniformly in $D$ and $n-1$ is negative (or positive) uniformly in
$D$. Then $N_{\sharp}=\big{|}\Re({N})\big{|}+\big{|}\Im({N})\big{|}$ is
positive, compact, injective and has a dense range.

  2. 2.

Assume that $\Im(A)$ is uniformly negative definite and $\Im(n)$ is uniformly
positive in $D$. Then $N_{\sharp}=\sigma\Im({N})$ is positive, compact,
injective and has a dense range where $\sigma=-1$.

####  3.3.3 Characterization of the scattering object

We now have all we need to prove the main results of this paper. Using Theorem
A.1 along with Picard’s criteria we will construct an indicator function to
reconstruct the support of $D$ via the factorization method which depends on
the analysis in the previous section. We also connect a modified linear
sampling to the factorization method. The indicator function proposed for
modified linear sampling here is motivated by the work done in [2]. It is
shown in [2] that the ‘minimizer’ $g_{z}^{\varepsilon}$ of the functional

$\mathcal{J}_{\varepsilon}\big{(}\Phi(\cdot\,,z)\,;\,g_{z}\big{)}=\varepsilon\big{(}N_{\sharp}\,g_{z},g_{z}\big{)}_{L^{2}(C)}+\|Ng_{z}-\Phi(\cdot\,,z)\|^{2}_{L^{2}(C)}$

(recall that the operator $N_{\sharp}$ is positive and injective) can be used
to derive an indicator function given by

$\big{(}N_{\sharp}\,g_{z}^{\varepsilon},g_{z}^{\varepsilon}\big{)}_{L^{2}(C)}=\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}$

which tends to infinity as the regularization parameter $\varepsilon$ tends to
zero for $z\notin D$ and is bounded otherwise. The square root of the operator
is defined by the spectral decomposition where

$\displaystyle N_{\sharp}^{1/2}g=\sum\lambda_{j}^{1/2}(g,\psi_{j})\psi_{j}$
(24)

for all $g\in L^{2}(C)$ where $(\lambda_{j},\psi_{j})\in\mathbb{R}^{+}\times
L^{2}(C)$ is the orthonormal eigensystem of $N_{\sharp}$. Here we extend this
idea and show that any regularized solution to

$N_{\sharp}\,g_{z}=\Phi(\cdot\,,z)$

satisfies that
$\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}_{L^{2}(C)}$ is bounded as
the regularization parameter $\varepsilon$ tends to zero if and only if $z\in
D$. In this section we will consider two cases, i.e. for an absorbing and non-
absording obstacle. We will state the main results separately for each case
for simplicity.

######  Theorem 3.6 (Range test for a non-absorbing media).

Assume that $\Im(A)=0$ and $\Im(n)=0$ in $D$. If $A-I$ is positive (or
negative) definite uniformly and $n-1$ is negative (or positive) uniformly in
$D$ then for $z\in\Omega$

$\Phi(\cdot\,,z)\in\text{Range}\big{(}N_{\sharp}^{1/2}\big{)}\quad\text{ if
and only if}\quad z\in D$

where $N_{\sharp}=|\Re(N)|+|\Im(N)|$.

###### Proof.

The result follows from the fact that Theorem A.1 along with Theorems 3.2 and
3.4 gives that $\text{Range}(N_{\sharp}^{1/2})=\text{Range}(H^{*})$ and then
by appealing to Theorem 3.3 proves the claim. ∎

We now give a range test for an absorbing media.

######  Theorem 3.7 (Range test for an absorbing media).

Assume that $\Im(A)$ is uniformly negative definite and $\Im(n)$ is uniformly
positive in $D$ then we have that for $z\in\Omega$

$\Phi(\cdot\,,z)\in\text{Range}\big{(}N_{\sharp}^{1/2}\big{)}\quad\text{ if
and only if}\quad z\in D$

where $N_{\sharp}=\sigma\Im({N})$ with $\sigma=-1$.

###### Proof.

Similarly we have that by using Theorem 4.1 in [24] for the imaginary part of
the operator. Then the result follows from Theorems 3.2 and 3.5 gives that
$\text{Range}(N_{\sharp}^{1/2})=\text{Range}(H^{*})$ and then by appealing to
Theorem 3.3 proves the claim. ∎

In either case we now let $\lambda_{j}\in\mathbb{R}^{+}$ and $\psi_{j}\in
L^{2}(C)$ be an orthonormal eigensystem of the positive, self-adjoint compact
operator ${N}_{\sharp}$ then applying Picard’s criterion (see for e.g. Theorem
2.7 of [8]) to Theorems 3.6 or 3.7 gives the following.

######  Corollary 3.2.

Assume that the assumptions of Theorem 3.6 or 3.7 are satisfied then for all
$z\in\Omega$

$\sum\limits_{j=1}^{\infty}\frac{\big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\big{|}^{2}}{\lambda_{j}}<\infty\quad\text{if
and only if }\quad z\in D.$

Corollary 3.2 gives a way to characterize the scatterer from the spectral data
of $N_{\sharp}$ and gives that

$W(z)=\left[\sum\limits_{j=1}^{\infty}\frac{\big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\big{|}^{2}}{\lambda_{j}}\right]^{-1}>0\quad\text{if
and only if }\quad z\in D.$

So one can reconstruct $D$ by plotting the function $W(z)$ since the support
of the function coincided with the scatterer.

We are now ready show that the modified linear sampling gives an equivalent
indicator function for reconstructing $D$. To this end, we define the modified
near field equation

$\displaystyle N_{\sharp}\,g_{z}=\Phi(\cdot\,,z)$ (25)

and the regularized solution to (25) given by

$\displaystyle
g^{\varepsilon}_{z}=\sum\limits_{j=1}^{\infty}\lambda_{j}f_{\varepsilon}(\lambda_{j}^{2})\,\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\,\psi_{j}$
(26)

where $\lambda_{j}\in\mathbb{R}^{+}$ and $\psi_{j}\in L^{2}(C)$ be an
orthonormal eigensystem of the positive, self-adjoint compact operator
${N}_{\sharp}$. The filter
$f_{\varepsilon}(t):\big{(}0,\lambda_{1}^{2}\big{]}\longmapsto\mathbb{R}^{+}$
for the regularization scheme is a piece-wise continuous sequence of functions
that satisfies for $0<t\leq\lambda_{1}^{2}$

$\lim\limits_{\varepsilon\to 0}f_{\varepsilon}(t)=\frac{1}{t}\quad\text{ and
}\quad tf_{\varepsilon}(t)\leq C_{\text{reg}}\quad\text{ for all
}\varepsilon>0$

where $\lambda_{1}$ is the largest eigenvalue of the operator $N_{\sharp}$.
Some common filter-functions are given by

$\displaystyle
f_{\varepsilon}(t)=\frac{1}{t+\varepsilon},\quad\quad\displaystyle{f(t)=\left\\{\begin{array}[]{lr}\frac{1}{t}&\,\,t>\varepsilon\\\
&\\\ 0&\,\,t\leq\varepsilon\end{array}\right.}\quad\text{and}\quad
f_{\varepsilon}(t)=\frac{1-(1-at)^{1/\varepsilon}}{t}$ (30)

where the constant $a<\lambda_{1}^{2}$. The filter functions in (30) are for
Tiknohov’s regularization, Spectral cutoff and Landweber iteration,
respectively. Notice that by (24) and (26)

$N_{\sharp}^{1/2}g_{z}^{\varepsilon}=\sum\limits_{j=1}^{\infty}\lambda^{3/2}_{j}f_{\varepsilon}(\lambda_{j}^{2})\,\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\,\psi_{j}$

which implies that

$\displaystyle\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}$
$\displaystyle=\sum\limits_{j=1}^{\infty}\lambda^{3}_{j}\big{[}f_{\varepsilon}(\lambda_{j}^{2})\big{]}^{2}\,\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}=\sum\limits_{j=1}^{\infty}\lambda^{4}_{j}\big{[}f_{\varepsilon}(\lambda_{j}^{2})\big{]}^{2}\,\frac{\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}}{\lambda_{j}}$
$\displaystyle\leq
C^{2}_{\text{reg}}\sum\limits_{j=1}^{\infty}\,\frac{\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}}{\lambda_{j}}.$

By Corollary 3.2 we have that
$\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}_{L^{2}(C)}$ is bounded as
$\varepsilon$ tends to zero for all $z\in D$. We now consider the case when
$z\notin D$ and notice that

$\displaystyle\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}$
$\displaystyle=\sum\limits_{j=1}^{\infty}\lambda^{3}_{j}\big{[}f_{\varepsilon}(\lambda_{j}^{2})\big{]}^{2}\,\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}$
$\displaystyle\geq\sum\limits_{j=1}^{M}\lambda^{3}_{j}\big{[}f_{\varepsilon}(\lambda_{j}^{2})\big{]}^{2}\,\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}.$

where $M\in\mathbb{N}$ is finite. Now, letting $\varepsilon\to 0$ then we
conclude that

$\displaystyle\lim\inf\limits_{\varepsilon\to
0}\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}\geq\sum\limits_{j=1}^{M}\,\frac{\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}}{\lambda_{j}}$

which implies that $\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}$ is
unbounded as $\varepsilon$ tends to zero since the right hand side of the
inequality becomes unbounded as $M$ increases for $z\notin D$. It is clear
that if $C_{\text{reg}}=1$ which is the case for the filter functions given in
(30) (for $\varepsilon$ small) then

$\sum\limits_{j=1}^{M}\,\frac{\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}}{\lambda_{j}}\leq\,\lim\limits_{\varepsilon\to
0}\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}\,\leq\sum\limits_{j=1}^{\infty}\frac{\Big{|}\big{(}\Phi(\cdot\,,z)\,,\psi_{j}\big{)}_{L^{2}(C)}\Big{|}^{2}}{\lambda_{j}}.$

This gives that in the limit the modified linear sampling method’s indicator
given by

$\lim\inf\limits_{\varepsilon\to
0}\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}=\lim\inf\limits_{\varepsilon\to
0}\big{(}N_{\sharp}\,g_{z}^{\varepsilon},g_{z}^{\varepsilon}\big{)}_{L^{2}(C)}$

is equivalent to the factorization method. The analysis done here follows for
any scatterer where either the near field or far field operator where the
factorization method is valid. Here one does not need the corresponding near
field or far field operator to have an eigenvalue decomposition as in [4] and
[6]. The above analysis gives that following results.

######  Theorem 3.8.

Assume that the assumptions of Theorem 3.6 or 3.7 is satisfied then the
regularized solution $g_{z}^{\varepsilon}$ of (25) for all $z\in\Omega$
satisfies that

$\lim\inf\limits_{\varepsilon\to
0}\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{2}_{L^{2}(C)}=\lim\inf\limits_{\varepsilon\to
0}\big{(}N_{\sharp}g_{z}^{\varepsilon},g_{z}^{\varepsilon}\big{)}_{L^{2}(C)}<\infty\quad\text{if
and only if }\quad z\in D.$

######  Corollary 3.3.

Assume that the assumptions of Theorem 3.6 or 3.7 is satisfied then the
regularized solution $g_{z}^{\varepsilon}$ of (25) for all $z\in\Omega$
satisfies that

$\lim\inf\limits_{\varepsilon\to
0}\big{\|}g_{z}^{\varepsilon}\big{\|}_{L^{2}(C)}=\infty\quad\text{ for all
}\quad z\notin D.$

Corollary 3.3 gives that any regularized solution to (25) becomes unbounded as
the regularization parameter $\varepsilon$ tends to zero. In practice the
indicator function for the modified linear sampling method can be taken to be
either

$P(z)=\big{\|}N_{\sharp}^{1/2}g_{z}^{\varepsilon}\big{\|}^{-1}_{L^{2}(C)}\quad\text{
or }\quad I(z)=\big{\|}g_{z}^{\varepsilon}\big{\|}^{-1}_{L^{2}(C)}$

where the support of these functions can be used to reconstruct the scatterer
$D$.

Due to the rigorous range characterization given by the factorization method,
we have the following uniqueness result for the inverse problem.

######  Corollary 3.4 (Uniqueness for the Inverse Problem).

Assume that $D_{1}$ and $D_{2}$ are two penetrable anisotropic obstacles
having material properties $A_{1}$, $n_{1}$ and $A_{2}$, $n_{2}$ that satisfy
the assumptions for either Theorem 3.6 or 3.7 in the interior of $D_{1}$ and
$D_{2}$, respectively. If the corresponding scattered fields are such that
$u^{s}_{1}(x,y)=u^{s}_{2}(x,y)$ for all $x,\,y\in C$ for a fixed wave number
$k$ (that satisfies assumption 3.1) then $D_{1}$ = $D_{2}$.

A numerical examples for unit ball in $\mathbb{R}^{2}$: We will now give an
explicit example of Theorem 3.6 and 3.8 where we assume that the coefficients
are isotropic homogeneous in a disk of radius one. This is a simple example to
numerically see the validity of the equivalence between the factorization
method and the modified linear sampling method. This is not an extensive
numerical study but we give this numerical evidence to illustrate the analysis
in this section. Therefore, the coefficient matrix is given by $A=aI$ where
both $a$ and $n$ are constant. Now notice that the scattering problem (11) can
be written as find $u\in H^{1}(B_{1})$ and $u^{s}\in
H^{1}_{loc}(\mathbb{R}^{2}\setminus\overline{B}_{1})$ such that

$\displaystyle a\Delta u+k^{2}nu=0\,\,\text{in
}\,\,B_{1}\quad\text{and}\quad\Delta u^{s}+k^{2}u^{s}=0\,\,\,$ in
$\displaystyle\,\,\mathbb{R}^{2}\setminus\overline{B}_{1}$ $\displaystyle
u^{s}-u=-\overline{\Phi(x,y)}\quad\textrm{ and }\quad\frac{\partial
u^{s}}{\partial r}-a\frac{\partial u}{\partial r}=-\frac{\partial}{\partial
r}\overline{\Phi(x,y)}\quad$ on $\displaystyle\partial B_{1},$

where $u$ is the total field in $B_{1}$ and $u^{s}$ is the radiating scattered
field in $\mathbb{R}^{2}\setminus\overline{B}_{1}$. We make the ansatz that
the solutions can be written as the following series

$u^{s}(x,y)=\sum\limits_{|m|=0}^{\infty}\alpha_{m}H^{(1)}_{m}(k|x|)\text{e}^{\text{i}m\hat{x}}\quad\text{and}\quad
u(x,y)=\sum\limits_{|m|=0}^{\infty}\beta_{m}\text{J}_{m}\left(k\sqrt{\frac{n}{a}}|x|\right)\text{e}^{\text{i}m\hat{x}}.$

Here we adopt the notation that points in $\mathbb{R}^{2}$ are given in by the
polar representation

$x=|x|\big{(}\cos(\hat{x}),\sin(\hat{x})\,\big{)}^{\top},$

where $\hat{x}\in[0,2\pi)$. Recall that for $|y|>|x|$ the fundamental solution
has the expansion

$\Phi(x,y)=\frac{\text{i}}{4}\sum\limits_{|m|=0}^{\infty}H^{(1)}_{m}(k|y|)\text{J}_{m}\left(k|x|\right)\text{e}^{\text{i}m(\hat{y}-\hat{x})}$

(see for e.g. [15]) where $\text{J}_{m}$ is the Bessel function and
$H^{(1)}_{m}$ is the first kind Hankel function of order $m$. We will assume
that $C=\partial B_{2}$ and therefore applying the boundary conditions we
obtain that

$\begin{bmatrix}H^{(1)}_{m}(k)&-\text{J}_{m}\left(k\sqrt{\frac{n}{a}}\right)\\\
(H^{(1)}_{m}(k))^{\prime}&-\sqrt{na}\text{J}{{}^{\prime}}_{m}\left(k\sqrt{\frac{n}{a}}\right)\end{bmatrix}\begin{bmatrix}\alpha_{m}\\\
\beta_{m}\end{bmatrix}=\frac{\text{i}}{4}H^{(2)}_{m}(2k)\text{e}^{-\text{i}m\hat{y}}\begin{bmatrix}\text{J}_{m}\left(k\right)\\\
\text{J}{{}^{\prime}}_{m}(k)\end{bmatrix}.$

Solving the linear system gives that

$\alpha_{m}=\frac{\text{i}}{4}H^{(2)}_{m}(2k)\text{e}^{-\text{i}m\hat{y}}\frac{\text{J}_{m}\left(k\sqrt{\frac{n}{a}}\right)\text{J}{{}^{\prime}}_{m}(k)-\sqrt{na}\text{J}{{}^{\prime}}_{m}\left(k\sqrt{\frac{n}{a}}\right)\text{J}_{m}(k)}{\text{J}_{m}\left(k\sqrt{\frac{n}{a}}\right)\left(H^{(1)}_{m}(k)\right)^{\prime}-\sqrt{na}\text{J}{{}^{\prime}}_{m}\left(k\sqrt{\frac{n}{a}}\right)H^{(1)}_{m}(k)}$

with $H^{(2)}_{m}$ being the second kind Hankel function of order $m$. Simple
manipulations give that the scattered field for $x,y\in\partial B_{2}$ can be
written as

$u^{s}(x,y)=\frac{\text{i}}{4}\sum\limits_{|m|=0}^{\infty}\sigma_{m}\big{|}H^{(1)}_{m}(2k)\big{|}^{2}\text{e}^{\text{i}m(\hat{x}-\hat{y})},$

where the polar angles $\hat{x},\hat{y}\in[0,2\pi)$ and we define

$\sigma_{m}=\frac{\text{J}_{m}\left(k\sqrt{\frac{n}{a}}\right)\text{J}{{}^{\prime}}_{m}(k)-\sqrt{na}\text{J}{{}^{\prime}}_{m}\left(k\sqrt{\frac{n}{a}}\right)\text{J}_{m}(k)}{\text{J}_{m}\left(k\sqrt{\frac{n}{a}}\right)\left(H^{(1)}_{m}(k)\right)^{\prime}-\sqrt{na}\text{J}{{}^{\prime}}_{m}\left(k\sqrt{\frac{n}{a}}\right)H^{(1)}_{m}(k)}.$

Given the above expression for the scattered field we have that the near field
operator can be written as $N:L^{2}(0,2\pi)\longmapsto L^{2}(0,2\pi)$

$(Ng)(\hat{x})=\frac{\text{i}}{4}\int_{0}^{2\pi}\sum\limits_{|m|=0}^{\infty}\sigma_{m}\big{|}H^{(1)}_{m}(2k)\big{|}^{2}\text{e}^{\text{i}m(\hat{x}-\hat{y})}\,g(\hat{y})\,d\hat{y}.$

In our examples we approximate the near field operator by truncating the
series. Let $N_{M}$ be the truncated approximation of the near field operator
for some $M\in\mathbb{N}$ given by

$(N_{M}g)(\hat{x})=\frac{\text{i}}{4}\int_{0}^{2\pi}\sum\limits_{|m|=0}^{M}\sigma_{m}\big{|}H^{(1)}_{m}(2k)\big{|}^{2}\text{e}^{\text{i}m(\hat{x}-\hat{y})}\,g(\hat{y})\,d\hat{y}.$

We apply Theorem 3.6 and 3.8 to the discretized version of the truncated near
field operator ${\bf N}_{M}$ where a simple 64 point Riemann sum approximation
is used to discretize the integral. Therefore we will plot the function

$W_{M}(z)=\left[\sum\limits_{i=1}^{64}\frac{\big{|}(\phi_{z}\,,\psi^{M}_{i})_{\ell^{2}}\big{|}^{2}}{\lambda^{M}_{i}}\right]^{-1}\quad\text{for}\quad\phi_{z}=\Phi(\cdot\,,z)$

where $\lambda^{M}_{i}$ and $\psi^{M}_{i}$ are the eigenvalues and vectors of
the matrix

$\big{(}{\bf N}_{M}\big{)}_{\sharp}=\big{|}\Re\big{(}{\bf
N}_{M}\big{)}\big{|}+\big{|}\Im\big{(}{\bf N}_{M}\big{)}\big{|}.$

Now, let the vector $g_{z}^{\varepsilon}$ be the solution to the discretized
modified near field equation by a spectral cutoff where we use the numerical
rank of the matrix as the cutoff parameter with indicator function

$P_{M}(z)=\Big{|}\Big{(}\big{(}{\bf
N}_{M}\big{)}_{\sharp}\,g_{z}^{\varepsilon},g_{z}^{\varepsilon}\Big{)}_{\ell^{2}}\Big{|}^{-1}.$

In the following examples we plot the functions $W_{M}(z)$ and $P_{M}(z)$ with
$M=20$.

Figure 6:  Reconstruction via the factorization and modified linear sampling
methods where $A=(1/2)I$ and $n=5$. The dotted red line is the actual boundary
of the unit circle. Figure 7: Reconstruction via the factorization and
modified linear sampling methods where $A=(3-\text{i})I$ and
$n=1/4+2\text{i}$. The dotted line is the actual boundary of the unit circle.

##  Appendix A Abstract theory of the factorization method

We now recall the main theoretical results from [27] which we use to develop
the range test in the previous section. To this end, let $X\subset U\subset
X^{*}$ be a Gelfand triple with a Hilbert space $U$ and a reflexive Banach
space $X$ such that the embedding is dense. Furthermore, let $Y$ be a second
Hilbert space and let ${N}:Y\mapsto Y$, $H:Y\mapsto X$ and $T:X\mapsto X^{*}$
be linear bounded operators such that $N=H^{*}TH$.

######  Theorem A.1.

(Theorem 2.1 in [27]) Assume that

  1. 1.

$H$ is compact and injective.

  2. 2.

$\Re(T)$ is the sum of a compact operator and a self adjoint coercive
operator.

  3. 3.

$\Im(T)$ is non-negative(or non-positive) on $X$.
Moreover assume that either of the following is satisfied:

  4. 4.

$T$ is injective.

  5. 5.

$\Im(T)$ is strictly positive(or negative) on the null space of $\Re(T)$.

Then the operator $N_{\sharp}=|\Re(\text{e}^{\text{i}t}N)|+|\Im(N)|$ is
positive, and the range of the operators $H^{*}:X^{*}\mapsto Y$ and
$N_{\sharp}^{1/2}:Y\mapsto Y$ coincide.

## References

  * [1]
  * [2] L. Audibert and H. Haddar,  A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements  Inverse Problems 30 (2014) 035011
  * [3]
  * [4] T. Arens,  Why linear sampling method works, Inverse Problems 20 (2004), 163-173 (2004).
  * [5]
  * [6] T. Arens and A. Lechleiter, Indicator Functions for Shape Reconstruction Related to the Linear Sampling Method, SIAM Journal on Imaging Sciences 8(1) (2015) 513-535.
  * [7] AS. Bonnet-BenDhia, L. Chesnel, H. Haddar  On the use of t-coercivity to study the interior transmission eigenvalue problem.  C. R. Acad. Sci., Ser. I 340: 647-651 (2011).
  * [8] F. Cakoni and D. Colton,  A Qualitative Approach to Inverse Scattering Theory Springer, Berlin 2014.
  * [9] F. Cakoni, D. Colton and H. Haddar,The linear sampling method for anisotropic media, J.Comput. Appl. Math., 146, 285-299 (2002).
  * [10] F. Cakoni, D. Colton, and H. Haddar,  The interior transmission problem for regions with cavities.  SIAM J. Math. Anal., 42(1):145–162, (2010).
  * [11] F. Cakoni, H. Haddar and I. Harris,  Homogenization of the transmission eigenvalue problem for periodic media and application to the inverse problem,  Inverse Problems and Imaging Volume 9, No. 4, 2015, 1025-1049.
  * [12] F. Cakoni and I. Harris, The factorization method for a defective region in an anisotropic material, Inverse Problems, 31, 025002 (2015).
  * [13] F. Cakoni and A. Kirsch,  On the interior transmission eigenvalue problem.  Int. Jour. Comp. Sci. Math. 3:142–167, (2010).
  * [14] M. Cheney, The linear sampling method and the MUSIC algorithm, Inverse Problems, 17 (2001) 591 595
  * [15] D. Colton and R. Kress,  Inverse Acoustic and Electromagnetic Scattering Theory.  Springer, New York, 3nd edition, 2013.
  * [16] D. Colton, J. Coyle and P. Monk, Recent developments in inverse acoustic scattering theory, SIAM Rev. 42 no. 3, 369 414, (2000).
  * [17] D. Gintides, M. Sini and N. Thanh  Detection of point-like scatterers using one type of scattered elastic waves  Journal of Computational and Applied Mathematics  236, 2137-2145 (2012).
  * [18] F. Gylys-Colwell, An inverse problem for the Helmholtz equation, Inverse Problems 12 (1996) 139 156.
  * [19] L. Evans,  _“Partial Differential Equations.”_ 2nd edition, AMS 2010.
  * [20] H. Haddar, A. Lechleiter, S. Marmorat,  An improved time domain linear sampling method for Robin and Neumann obstacles,  Applicable Analysis , 93, issue 2, 2014.
  * [21] H. Heck, G. Nakamura and H. Wang,  Linear sampling method for identifying cavities in a heat conductor,  Inverse Problems, 28, (2012) 075014.
  * [22] G. Hu, J. Yang, B. Zhang and H. Zhang,  Near-field imaging of scattering obstacles with the factorization method  Inverse Problems 30 (2014) 095005.
  * [23] A. Kirsch, An integral equation for Maxwell s equations in a layered medium with an application to the factorization method, Journal of Integral Equations and Applications 19, (2007) 333-357.
  * [24] A. Kirsch, The MUSIC algorithm and the factorization method in inverse scattering theory for inhomogeneous media, Inverse Problems 18 , (2002) 1025-1040.
  * [25] A. Kirsch and N. Grinberg, _The Factorization Method for Inverse Problems_. Oxford University Press, Oxford 2008.
  * [26] A. Kirsch and X. Liu, The factorization method for inverse acoustic scattering by a penetrable anisotropic obstacle, Math. Meth. Appl. Sci., 37: 1159-1170.
  * [27] A. Lechleiter, The factorization method is independent of transmission eigenvalues, Inverse Problems and Imaging 3 (2009), 123 138.
  * [28] X. Liu and J. Sun, Reconstruction of Neumann eigenvalues and support of sound hard obstacles Inverse Problems 30  (2014) 065011
  * [29] S. Meng, F. Cakoni and H. Haddar, The factorization method for a cavity in an inhomogeneous medium, Inverse Problems, 30, 045008 (2014)
  * [30] C. Miranda, Partial Differential Equations of Elliptic Type, Springer, Berlin, 1970.
  * [31] J. C. Nedelec, Acoustic and Electromagnetic Equations Integral Representations for Harmonic Problems (Applied Mathematical Sciences vol 144). New York: Springer 2001
  * [32] Salvatier J, Wiecki TV, Fonnesbeck C. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:55 (2016) https://doi.org/10.7717/peerj-cs.55
  * [33] L. Poli, G. Oliveri, and A. Massa. Microwave imaging within the first-order Born approximation by means of the contrast-field Bayesian compressive sensing. IEEE Trans. Antennas Propagat. 60(6): 2865 2879, June 2012.
  * [34] E. A. Marengo, R. D. Hernandez, Y. R. Citron, F. K. Gruber, M. Zambrano, and H. Lev-Ari. Compressive sensing for inverse scattering. Proc. XXIX URSI Gen. Assem. Jan. 7-16 (2008).
  * [35] E. A. Marengo. Compressive sensing and signal subspace methods for inverse scattering including multiple scattering. Proc. of the IEEE Geoscience Remote Sensing Symp. Jul. 7-11 2008.
  * [36] E. A. Marengo. Subspace and Bayesian compressive sensing methods in imaging. Proc. Prog. Electromagn. Res. Symp. Jul. 2-6 (2008).
  * [37] A. E. Fouda and F. L. Teixeira, Bayesian compressive sensing for ultrawideband inverse scattering in random media. Inv. Probl. 30(11): 114017 (2014). http://arxiv.org/abs/1401.1092