Analysis of the monotonicity method for an anisotropic scatterer with a
conductive boundary

Victor Hughes, Isaac Harris, and Heejin Lee
Department of Mathematics, Purdue University, West Lafayette, IN 47907
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>

###### Abstract

In this paper, we consider the inverse scattering problem associated with an
anisotropic medium with a conductive boundary. We will assume that the
corresponding far–field pattern is known/measured and we consider two inverse
problems. First, we show that the far–field data uniquely determines the
boundary coefficient. Next, since it is known that anisotropic coefficients
are not uniquely determined by this data we will develop a qualitative method
to recover the scatterer. To this end, we study the so–called monotonicity
method applied to this inverse shape problem. This method has recently been
applied to some inverse scattering problems but this is the first time it has
been applied to an anisotropic scatterer. This method allows one to recover
the scatterer by considering the eigenvalues of an operator associated with
the far–field operator. We present some simple numerical reconstructions to
illustrate our theory in two dimensions. For our reconstructions, we need to
compute the adjoint of the Herglotz wave function as an operator mapping into
$H^{1}$ of a small ball.

Keywords: Inverse Scattering $\cdot$ Monotonicity Method $\cdot$ Far–Field
Operator $\cdot$ Anisotropic Scatterer

MSC: 35J05, 35J25

##  1 Introduction

Here, we will study the inverse shape problem of recovering the scatterer as
well as the inverse boundary parameter problem from given far–field data. For
the inverse shape problem, we analyze the monotonicity method for
reconstructing an anisotropic medium with a conductive boundary. The
monotonicity method is a qualitative reconstruction scheme for the unknown
region that is comparable to the so–called linear sampling [11, 16] and
factorization methods [29, 30]. This is do to the fact that each method
derives an imaging function that requires the eigenvalue(or singular) value
decomposition of an operator associated with the known/measured far–field
operator. The monotonicity method does not require us to have any restriction
on the wave number i.e. valid even at transmission eigenvalues(see for e.g.
[6, 12, 28, 34]). This is not true for the linear sampling and factorization
methods but here we need to assume we know the sign of the anisotropic
contrast. Our main contribution will be to study the application of the
monotonicity method for an anisotropic material with a conductive boundary on
an unbounded domain. This problem has not been studied before and will require
original analysis to complete. Indeed, one of the main calculations needed to
apply the monotonicity method to an anisotropic scatterer is to compute the
adjoint of the Herglotz wave function as an operator mapping into
$L^{2}(\mathbb{S}^{d-1})\to H^{1}\big{(}B(z,\epsilon)\big{)}$. Here
$B(z,\epsilon)$ denotes the ‘sampling ball’ at the ‘sampling point’
$z\in\mathbb{R}^{d}$. The main result of this paper states that for a known
operator associated with the measurements only has finitely many negative
eigenvalues provided that the sampling ball is contained in the scatterer.
This gives us a way to recover the scatterer via computing the eigenvalues of
a known operator with little a priori information on the geometry. Since the
imaging function is given by computing the eigenvalues of known operators,
this would imply stability by the continuity of the spectrum.

Here the conductive boundary condition is modeled by a Robin condition that
states the normal/co-normal derivative of the total field has a jump across
the boundary of the scatterer that is proportional to the total field. We
remark that this is different from the model presented in [5] but was recently
studied in [28, 34]. The aforementioned papers consider the associated
transmission eigenvalue problem. We wish to fill the current gap in knowledge
i.e. to study the associated inverse problems for this model. For this model,
we will need to derive a factorization of the far–field operator. See [31] for
a a factorization of the far–field operator with no conductivity parameter. We
note that our analysis is still valid for the case when there is no conductive
boundary parameter. This is then used in conjunction with the main result in
[18] to derive a monotonicity method. In [18] the authors give a theoretical
blueprint for applying the monotonicity method to inverse scattering problems.
This method has been studied in [20] for other inverse media scattering
problems as well as in [1] for some inverse obstacle scattering problems. The
monotonicity method has been used in shape reconstruction for the $p$–Laplace
[3], electrical impedance tomography [23] as well as other
models/configurations see for e.g. [2, 8, 17, 22, 33]. This paper adds to the
growing literature on the monotonicity method applied to many imaging
modalities.

The rest of the paper is organized as follows: In Section 2, we formulate the
direct and inverse scattering problem for an anisotropic material with a
conductive boundary. Next, in Section 3 we study the uniqueness of the inverse
parameter problem for determining the boundary parameter. We are able to prove
that the far–field pattern uniquely determines the boundary parameter. In
section 4, we decompose the far–field operator into a symmetric factorization,
this is necessary for the use of the monotonicity method. Once an an
appropriate factorization is derived, in Section 5 we provide an explanation
on how to apply the monotonicity method to this problem given our
factorization of the far field operator using the theory developed in [18].
Finally, in Section 6, we provide numerical examples to validate our
theoretical results.

##  2 Formulation of the Problem

We start by discussing the direct scattering problem for an anisotropic
material with a conductive boundary. Let $D\subset\mathbb{R}^{d}$ for $d=$ 2
or 3 be a simply connected open set with Lipschitz boundary $\partial D$,
where $\nu$ represents the unit outward normal vector. The region $D$ denotes
the scatterer that we illuminate with an incident plane wave
$u^{i}\mathrel{\mathop{\mathchar 58\relax}}=\text{e}^{\text{i}kx\cdot\hat{y}}$
such that $k>0$ is the wave number and $\hat{y}\in\mathbb{S}^{d-1}$ denotes
the incident direction. We assume that we have a symmetric matrix-valued
function $A(x)\in\mathscr{C}^{1}(\overline{D},\mathbb{C}^{d\times d})$ that is
uniformly positive definite in $D$ satisfying

$\overline{\xi}\cdot\text{Re}(A(x))\xi\geq A_{\text{min}}|\xi|^{2}$ and
$\overline{\xi}\cdot\text{Im}(A(x))\xi\leq 0$ for a.e. $x\in D$ and
$\xi\in\mathbb{C}^{d}$

where $A_{\text{min}}$ is a positive constant. We also assume that the
refractive index $n(x)\in L^{\infty}(\mathbb{R}^{d})$ satisfies

$\text{Im}(n(x))\geq 0$ for a.e. $x\in D$.

We assume that the parameters $A(x)$ and $n(x)$ for an anisotropic material
satisfy

$A(x)=I\quad\text{and}\quad n(x)=1\quad\text{for
}x\in\mathbb{R}^{d}\setminus\overline{D}.$

We also consider the conductivity parameter $\eta\in L^{\infty}(\partial D)$,
with the condition

$\text{Im}(\eta(x))\geq 0$ for a.e. $x\in\partial D$.

We have that the direct scattering problem for an anisotropic material with a
conductive boundary condition is formulated as follows: find $u^{s}\in
H^{1}_{loc}(\mathbb{R}^{d})$ such that

$\displaystyle\Delta u^{s}+k^{2}u^{s}=0\quad\text{in
$\mathbb{R}^{d}\setminus\overline{D}$}\quad\text{ and }\quad\nabla\cdot
A\nabla u+k^{2}nu=0\quad$ in $D$ (1)
$\displaystyle(u^{s}+u^{i})^{+}=u^{-}\quad\text{ and
}\quad\partial_{\nu}(u^{s}+u^{i})^{+}=\nu\cdot A\nabla u^{-}-\eta u\quad$ on
$\partial D$ (2)

where $u=u^{s}+u^{i}$ denotes the total field. Here the superscripts $+$ and
$-$ demonstrate approaching the boundary from the outside and inside of $D$,
respectively. To close the system, we assume that the scattered field
satisfies the radiation condition

$\displaystyle\partial_{r}u^{s}-\text{i}ku^{s}=\mathcal{O}\Big{(}\frac{1}{r^{(d+1)/2}}\Big{)}\quad\text{
as }\quad r\to\infty$ (3)

uniformly with respect to $\hat{x}=\frac{x}{r}$ where $r=|x|$. With the above
assumptions on the coefficients we can prove well-posedness of (1)–(2) along
with the radiation condition. This can be done in a similar manner to the
analysis in [4] where a variational formulation of (1)–(2) for the scattered
field $u^{s}$ is used to study the system.

We are interested in recovering the scatterer $D$ and/or the boundary
parameter $\eta$ from measured far–field data. To this end, recall that the
asymptotic expansion of the scattering field $u^{s}(\cdot\,,\hat{y})$ has the
form

$\displaystyle
u^{s}(x,\hat{y})=\gamma_{d}\frac{\text{e}^{\text{i}k|x|}}{|x|^{\frac{d-1}{2}}}\Bigg{\\{}u^{\infty}(\hat{x},\hat{y})+\mathcal{O}\Bigg{(}\frac{1}{|x|}\Bigg{)}\Bigg{\\}}$

where $\gamma_{d}$ is given by

$\gamma_{2}=\frac{\text{e}^{\text{i}\pi/4}}{\sqrt{8\pi
k}}\quad\text{and}\quad\gamma_{3}=\frac{1}{4\pi}.$

Here we let $\hat{x}=x/|x|$ and $u^{\infty}(\hat{x},\hat{y})$ is the
corresponding far–field pattern associated with (1)–(2), which is dependent on
the incident direction $\hat{y}$ and the observation direction $\hat{x}$. Due
to the fact that $u^{s}$ is a radiating solution to the Helmholtz equation in
$\mathbb{R}^{d}\setminus\overline{D}$ we have that it can be written using
Green’s representation formula [10]. Therefore, we have that the far–field
pattern has the following integral representation

$\displaystyle
u^{\infty}(\hat{x},\hat{y})=-\int_{\partial\Omega}\text{e}^{-\text{i}k\hat{x}\cdot
z}\partial_{\nu}u^{s}-u^{s}\partial_{\nu}\text{e}^{-\text{i}k\hat{x}\cdot
z}\,\mathrm{d}s(z)$

for any region $\Omega$ such that $\overline{D}\subseteq\Omega$. With the
far–field pattern, we now define the far–field operator
$F\mathrel{\mathop{\mathchar 58\relax}}L^{2}(\mathbb{S}^{d-1})\longrightarrow
L^{2}(\mathbb{S}^{d-1})$ as

$\displaystyle(Fg)(\hat{x})\mathrel{\mathop{\mathchar
58\relax}}=\int_{\mathbb{S}^{d-1}}u^{\infty}(\hat{x},\hat{y})g(\hat{y})\,\text{d}s(\hat{y})\quad\text{
where }g\in L^{2}(\mathbb{S}^{d-1}).$

It is well known that the far–field operator $F$ is compact.

We are interested in the inverse problem of recovering $D$ from the given the
far–field data. The main goal is to apply the monotonicity method to perform
this reconstruction of $D$, given we have the far–field data. Another inverse
problem that we will consider is the uniqueness for recovering the boundary
parameter $\eta$ from the known far–field data.

##  3 Uniqueness of the Inverse Parameter Problem

In this section, we show that the far–field pattern uniquely determines the
boundary parameter $\eta$ provided that $D$ is known. In general, the
anisotropic material parameter $A$ is not uniquely recovered by the far–field
pattern(see for e.g. [21]). Therefore, we will consider a qualitative method
to recover $D$ in the preceding section. Here we will prove a uniqueness
result for the inverse boundary parameter problem. We are considering the
injectivity of the mapping

$\eta\in L^{\infty}(\partial D)\longmapsto u^{\infty}(\hat{x},\hat{y})\in
L^{2}\big{(}\mathbb{S}^{d-1}\times\mathbb{S}^{d-1}\big{)}.$

To this end, in addition to the assumption on $\eta$ in Section 2, we assume
that

$\displaystyle\varphi\in H^{1}_{0}(D)\quad\text{satisfies}\quad\nabla\cdot
A\nabla\varphi+k^{2}n\varphi=0\,\,\text{in $D$}\quad\text{implies
$\varphi=0\,\,$ a.e. in $D$ }$ (4)

and we need a density result as in [13, 15]. Notice, that (4) is true for all
$k\in\mathbb{R}$ provided that $A$ and $n$ are complex–valued.

######  Lemma 3.1.

Under the assumption (4), the set

$\displaystyle\mathcal{U}\coloneqq\left\\{u(\cdot,\hat{y})|_{\partial
D}\mathrel{\mathop{\mathchar 58\relax}}u\textrm{ is the total field of
\eqref{directacbc1}--\eqref{directacbc2} for
}\hat{y}\in\mathbb{S}^{d-1}\right\\}$

is dense in $L^{2}(\partial D)$.

###### Proof.

To prove the claim, we will show that $\mathcal{U}^{\perp}$ is trivial. For
any $\phi\in\mathcal{U}^{\perp}$ and let $p\in H^{1}_{loc}(\mathbb{R}^{d})$ be
the unique radiating solution of

$\displaystyle\Delta p+k^{2}p=0\quad\text{in
$\mathbb{R}^{d}\setminus\overline{D}$}\quad\text{ and }\quad\nabla\cdot
A\nabla p+k^{2}np=0\quad$ in $D$ $\displaystyle p^{+}=p^{-}\quad\text{ and
}\quad\partial_{\nu}p^{+}-\nu\cdot A\nabla p^{-}+\eta p=\phi\quad$
$\displaystyle\text{on $\partial D$}.$

If $u$ is the total field of (1)–(2), from the boundary conditions, we have
that

$\displaystyle\int_{\partial
D}p\partial_{\nu}u^{+}-u\partial_{\nu}p^{+}\,\mathrm{d}s=\int_{\partial
D}p(\nu\cdot A\nabla u^{-})-u(\nu\cdot A\nabla p^{-})-u\phi\,\mathrm{d}s.$

Notice that from Green’s 2nd identity and the symmetry of $A$,

$\displaystyle\int_{\partial D}p(\nu\cdot A\nabla u^{-})-u(\nu\cdot A\nabla
p^{-})\mathrm{d}s=\int_{D}\nabla p\cdot A\nabla u-\nabla u\cdot A\nabla
p\,\mathrm{d}x=0.$

Therefore, combining the above two equations, we have that

$\displaystyle\int_{\partial
D}p\partial_{\nu}u^{+}-u\partial_{\nu}p^{+}\,\mathrm{d}s=-\int_{\partial
D}u\phi\,\mathrm{d}s=0,$

that is,

$\displaystyle\int_{\partial
D}p\partial_{\nu}u^{i}-u^{i}\partial_{\nu}p^{+}\,\mathrm{d}s=\int_{\partial
D}u^{s}\partial_{\nu}p^{+}-p\partial_{\nu}u^{s}\,\mathrm{d}s.$

Applying Green’s 2nd identity to $u^{s}$ and $p$ in
$B_{R}\setminus\overline{D}$ where $B_{R}$ is the ball of radius $R>0$
centered at the origin and containing $D$ in its interior, we have

$\displaystyle\int_{\partial
D}u^{s}\partial_{\nu}p^{+}-p\partial_{\nu}u^{s}\,\mathrm{d}s=-\int_{B_{R}}u^{s}\partial_{\nu}p-p\partial_{\nu}u^{s}\,\mathrm{d}s\rightarrow
0\,\,\text{ as }\,\,R\rightarrow\infty,$

since $u^{s}$ and $p$ are satisfy the radiation condition (3). Thus,

$\displaystyle\int_{\partial
D}p\partial_{\nu}u^{i}-u^{i}\partial_{\nu}p^{+}\,\mathrm{d}s=0.$ (5)

Recall, that the fundamental solution $\Phi(x,y)$ of the Helmholtz equation is
given by

$\displaystyle\Phi(x,y)=\frac{\text{i}}{4}H^{(1)}_{0}(k|x-y|)\text{ for
}d=2\quad\text{and}\quad\Phi(x,y)=\frac{\text{e}^{\text{i}k|x-y|}}{4\pi|x-y|}\text{
for }d=3,$ (6)

where $H_{0}^{(1)}$ is the Hankel function of the first kind and of order $0$.
Now, consider the incident field
$u^{i}(x,-\hat{y})=\text{e}^{-\text{i}kx\cdot\hat{y}}$ which is the far field
pattern of $\Phi(y,x)$. From the Green’s representation formula(see for e.g.
[10]),

$\displaystyle p(y)=\int_{\partial
D}p\partial_{\nu}\Phi(\cdot,y)-\partial_{\nu}p^{+}\Phi(\cdot,y)\,\mathrm{d}s\quad\text{
for }\,\,\,y\in\mathbb{R}^{d}\setminus\overline{D}.$

This implies that left-hand side of (5) represents the far–field pattern
$p^{\infty}(\hat{y})$ of the function $p(y)$. Since $p^{\infty}(\hat{y})=0$
from (5), we obtain that $p(y)$ vanishes in
$\mathbb{R}^{d}\setminus\overline{D}$ by Rellich’s lemma(see for e.g. [10]).
This implies that $p^{-}=0$ on $\partial D$. By assumption (4), we conclude
that $p\equiv 0$ in $D$ which completes the proof. ∎

We showed the density of the set $\mathcal{U}$ in $L^{2}(\partial D)$. With
this, we will show that the far–field pattern uniquely determines the boundary
parameter under the same assumption as in Lemma 3.1.

######  Theorem 3.1.

Under the assumption (4), the knowledge of the far–field pattern
$u^{\infty}(\hat{x},\hat{y})$ for all $\hat{x}$ and
$\hat{y}\in\mathbb{S}^{d-1}$ uniquely determines the boundary parameter
$\eta\in L^{\infty}(\partial D)$.

###### Proof.

Let us denote $u_{j}(\cdot\,,\hat{y})$ as the total field to (1)–(2) with the
boundary parameter $\eta$ replaced by $\eta_{j}$ for each $j=1,2$, and let
$u_{j}^{s}$ be the corresponding scattered field. Assuming that the far–field
patterns of $u_{1}^{\infty}(\hat{x},\hat{y})=u_{2}^{\infty}(\hat{x},\hat{y})$
for all $\hat{x}$ and $\hat{y}\in\mathbb{S}^{d-1}$. Then by Rellich’s lemma we
have that $u_{1}^{s}({x},\hat{y})=u_{2}^{s}({x},\hat{y})$ for all
$x\in\mathbb{R}^{d}\setminus\overline{D}$ and $\hat{y}\in\mathbb{S}^{d-1}$,
which implies that $u_{1}(\cdot\,,\hat{y})=u_{2}(\cdot\,,\hat{y})$ on
$\partial D$. With assumption (4), $u_{1}(\cdot\,,\hat{y})$ and
$u_{2}(\cdot\,,\hat{y})$ coincide in $\mathbb{R}^{d}$.

Now, let us denote $u_{1}$ and $u_{2}$ by $u$. From the boundary condition,

$\displaystyle 0=\partial_{\nu}u^{+}-\nu\cdot A\nabla
u^{-}+\eta_{1}u=\partial_{\nu}u^{+}-\nu\cdot A\nabla u^{-}+\eta_{2}u\text{ on
}\partial D.$

Then, $(\eta_{1}-\eta_{2})u(\cdot\,,\hat{y})=0$ on $\partial D$. Since
$(\eta_{1}-\eta_{2})\in L^{\infty}(\partial D)\subset L^{2}(\partial D)$, it
follows from Lemma 3.1 that $\eta_{1}=\eta_{2}$ a.e. on $\partial D$. ∎

With this result, we know that the inverse boundary parameter problem is
unique with respect to the data. In the following section we will solve the
inverse shape problem i.e. recover the unknown scatterer $D$ from the given
far–field pattern/operator.

##  4 Factorization of Far Field Operator

In order to utilize the monotonicity method for reconstructing the region $D$,
we will first need to obtain a symmetric factorization of the far–field
operator. To this end, motivated by the original direct scattering problem
(1)–(2), we consider the problem of finding $u\in H^{1}_{loc}(\mathbb{R}^{d})$
for any 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$ $\displaystyle\text{ in }\mathbb{R}^{d}\setminus\partial
D$ (7) $\displaystyle u^{+}-u^{-}=0\quad\text{ and
}\quad\partial_{\nu}u^{+}-\nu\cdot A\nabla u^{-}=\nu\cdot(A-I)\nabla
v-\eta(u+v)\quad$ $\displaystyle\text{ on }\partial D$ (8)
$\displaystyle\partial_{r}u-\text{i}ku=\mathcal{O}\Big{(}\frac{1}{r^{(d+1)/2}}\Big{)}\quad$
$\displaystyle\text{ as }r\to\infty.$

Note that $I-A$ and $1-n$ are supported in $D$. This system is equivalent to
(1)–(2) when you use the fact the incident field solves the Helmholtz equation
in all of $\mathbb{R}^{d}$. The variational formulation for the above system
is given by

$\displaystyle-\int_{B_{R}}A\nabla
u\cdot\nabla\overline{\phi}-k^{2}nu\overline{\phi}\,\text{d}x$
$\displaystyle+\int_{\partial B_{R}}\overline{\phi}\Lambda_{k}u\,\text{d}s$
$\displaystyle=\int_{D}\nabla\overline{\phi}\cdot(A-I)\nabla
v-k^{2}(n-1)v\overline{\phi}\,\text{d}x-\int_{\partial
D}\overline{\phi}\eta(u+v)\,\text{d}s$ (9)

for all $\phi\in H^{1}(B_{R})$, where $B_{R}$ is the open ball of radius $R$
such that $\overline{D}\subset B_{R}$. Here we let

$\Lambda_{k}\mathrel{\mathop{\mathchar 58\relax}}H^{1/2}(B_{R})\longrightarrow
H^{-1/2}(B_{R})$

denote the Dirichlet to Neumann map on $\partial B_{R}$ for the radiating
solution to the Helmholtz equation on the exterior of $B_{R}$, defined as

$\Lambda_{k}f=\partial_{\nu}\varphi\big{|}_{\partial B_{R}}\quad\text{ where
}\quad\Delta\varphi+k^{2}\varphi=0\text{ in
$\mathbb{R}^{d}\setminus\overline{B_{R}}$}\quad\text{ and
}\quad\varphi\big{|}_{\partial B_{R}}=f,$

along with the radiation condition (3). From Theorem 5.22 of [9] we know that
$\Lambda_{k}$ is a bounded linear operator. This is a direct consequence of
the well-posedness of the above Dirichlet problem along with the (Neumann)
Trace Theorem.

To get an initial factorization of $F$ we first define the source to far–field
pattern operator

$\displaystyle G\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\longrightarrow
L^{2}(\mathbb{S}^{d-1})\quad\text{given by}\quad Gv=u^{\infty}$ (10)

where $u\in H^{1}_{loc}(\mathbb{R}^{d})$ is the unique solution to (7)–(8) for
a given $v\in H^{1}(D)$. Next we define the bounded linear operator

$\displaystyle H\mathrel{\mathop{\mathchar
58\relax}}L^{2}(\mathbb{S}^{d-1})\longrightarrow H^{1}(D)\quad\text{given
by}\quad(Hg)(x)=\int_{\mathbb{S}^{d-1}}\text{e}^{\text{i}kx\cdot\hat{y}}g(\hat{y})\,\text{d}s(\hat{y})\Big{|}_{D}.$
(11)

By the superposition principle, we have that the far–field operator associated
with (1)–(2) is given by $F=GH$. In order to use the monotonicity method, we
need to have a symmetric factorization of far–field operator. To achieve this,
we compute the adjoint of the operator $H$.

######  Theorem 4.1.

The operator $H^{*}\mathrel{\mathop{\mathchar
58\relax}}H^{1}(D)\longrightarrow L^{2}(\mathbb{S}^{d-1})$ is given by

$H^{*}f=-w^{\infty}\quad\text{for all }\quad f\in H^{1}(D)$

where $w^{\infty}$ is the far–field pattern for the unique radiating solution
$w$ that satisfies

$\displaystyle-\int_{B_{R}}\nabla
w\cdot\nabla\overline{\phi}-k^{2}w\overline{\phi}\,\mathrm{d}x+\int_{\partial
B_{R}}\overline{\phi}\Lambda_{k}w\,\mathrm{d}s=(f,\phi)_{H^{1}(D)}$ (12)

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

###### Proof.

To prove the claim, we let $\psi=Hg$ be extended for all $x\in\mathbb{R}^{d}$.
Then we have that

$\displaystyle(H^{*}f,g)_{L^{2}(\mathbb{S}^{d-1})}=(f,Hg)_{H^{1}(D)}=-\int_{B_{R}}\nabla
w\cdot\nabla\overline{\psi}-k^{2}w\overline{\psi}\,\mathrm{d}x+\int_{\partial
B_{R}}\overline{\psi}\partial_{\nu}w\,\mathrm{d}s.$

Using Green’s 1st Identity and along with the fact that $\psi$ is a solution
to the Helmholtz equation in $\mathbb{R}^{d}$ we now obtain

$\displaystyle(H^{*}f,g)_{L^{2}(\mathbb{S}^{d-1})}$
$\displaystyle=\int_{B_{R}}w\overline{(\Delta\psi+k^{2}\psi)}\,\mathrm{d}x+\int_{\partial
B_{R}}\overline{\psi}\partial_{\nu}w-w\partial_{\nu}\overline{\psi}\,\mathrm{d}s$
$\displaystyle=\int_{\partial
B_{R}}\overline{\psi}\partial_{\nu}w-w\partial_{\nu}\overline{\psi}\,\mathrm{d}s$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\overline{g(\hat{y})}\Bigg{[}\int_{\partial
B_{R}}\text{e}^{-\text{i}kx\cdot\hat{y}}\partial_{\nu}w-w\partial_{\nu}\text{e}^{-\text{i}kx\cdot\hat{y}}\,\mathrm{d}s(x)\Bigg{]}\,\mathrm{d}s(\hat{y})$
$\displaystyle=-\int_{\mathbb{S}^{d-1}}\overline{g(\hat{y})}w^{\infty}(\hat{y})\,\mathrm{d}s(\hat{y})=(-w^{\infty},g)_{L^{2}(\mathbb{S}^{d-1})}.$

This implies that $H^{*}f=-w^{\infty}$, proving the claim. ∎

Now that we have the adjoint of our operator $H$, we proceed as is done in
[14, 27] to show that $G=-H^{*}T$, for some operator $T$. To this end, from
(7) we have that

$\Delta u+k^{2}u=\nabla\cdot(I-A)\nabla(u+v)+k^{2}(1-n)(u+v)\quad\text{ in
}\mathbb{R}^{d}\setminus\partial D$

where $u\in H^{1}_{loc}(\mathbb{R}^{d})$ is the solution of (7)–(8) for any
given $v\in H^{1}(D)$. The corresponding variational form is given by

$\displaystyle-\int_{B_{R}}\nabla u\cdot\nabla\overline{\phi}$
$\displaystyle-k^{2}u\overline{\phi}\,\mathrm{d}x+\int_{\partial
B_{R}}\overline{\phi}\Lambda_{k}u\,\mathrm{d}s$
$\displaystyle=\int_{D}\nabla\overline{\phi}\cdot(A-I)\nabla(u+v)-k^{2}(n-1)(u+v)\overline{\phi}\,\mathrm{d}x-\int_{\partial
D}\overline{\phi}\eta(u+v)\,\mathrm{d}s$ (13)

for any $\phi$. By wellposedness of the direct scattering problem, we have
that $v\longmapsto u$ is a bounded linear mapping. With this fact we can use
the Riesz Representation Theorem on the right hand side of the variational
form (13) which implies that there exists a bounded linear operator
$T\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\longrightarrow H^{1}(D)$ such
that for all $v\in H^{1}(D)$

$\displaystyle(Tv,\phi)_{H^{1}(D)}=\int_{D}\nabla\overline{\phi}\cdot(A-I)\nabla(u+v)-k^{2}(n-1)(u+v)\overline{\phi}\,\mathrm{d}x-\int_{\partial
D}\overline{\phi}\eta(u+v)\,\mathrm{d}s.$ (14)

It is clear by (13), that we have

$\displaystyle-\int_{B_{R}}\nabla
u\cdot\nabla\overline{\phi}-k^{2}u\overline{\phi}\,\mathrm{d}x+\int_{\partial
B_{R}}\overline{\phi}\Lambda_{k}u\,\mathrm{d}s=(Tv,\phi)_{H^{1}(D)}.$

By appealing to the definition of $G$ and $H^{*}$ as defined in (10) and (12)
we have that

$u^{\infty}=Gv\quad\text{ and }\quad u^{\infty}=-H^{*}Tv\quad\text{ for any
$v\in H^{1}(D)$ }$

which implying $G=-H^{*}T$. Now, recalling that we have the initial
factorization $F=GH$, which implies the desired symmetric factorization
$F=-H^{*}TH$.

######  Lemma 4.1.

The far-field operator $F\mathrel{\mathop{\mathchar
58\relax}}L^{2}(\mathbb{S}^{d-1})\longrightarrow L^{2}(\mathbb{S}^{d-1})$ for
the scattering problem (1)–(2) has the factorization $F=-H^{*}TH$.

With this result we have the factorization required to apply the theory of the
monotonicity method to our inverse problem. In the next section we will
develop the required theoretical results. This will require us to study the
analytical properties of the operators $H$ and $T$ that are used to factorize
the far-field operator.

##  5 The Monotonicity Method

In this section, we connect the support of the scattering object $D$ to the
spectrum of an operator associated with the far–field operator $F$. We recall
the symmetric factorization of the far–field operator $F=-H^{*}TH$, this
factorization is important for applying the theory of the monotonicity method.
We now provide the following definition that will give some background for our
analysis going forward.

######  Definition 5.1.

Let $A$, $B$ be self-adjoint compact operators on a Hilbert space. We write

$A\leq_{\mathrm{fin}}B$

If $B-A$ had finitely many negative eigenvalues.

The main result of the monotonicity method that we will employ here is given
in Theorem 4.2 of [18]. This result outlines how the number of negative
eigenvalues of an operator associated with the far-field operator will
increase provided that the sampling point is in the interior of the scatterer
$D$. For completeness we state a version of Theorem 4.2 in [18] for our
setting.

######  Theorem 5.1.

(Theorem 4.2 of [18]) Assume we have compact operators
$\mathcal{F}_{j}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{Y}\rightarrow\mathbb{Y}$ acting on the Hilbert space
$\mathbb{Y}$ with factorizations

$\mathcal{F}_{j}=H^{*}_{j}T_{j}H_{j}\quad\text{ such that }\quad
H_{j}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{Y}\to\mathbb{X}_{j}\quad\text{ and }\quad
T_{j}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{X}_{j}\to\mathbb{X}_{j}$

with $T_{j}$ and $H_{j}$ being bounded linear operators with Hilbert spaces
$\mathbb{X}_{j}$ for $j=1$ and 2.

  1. 1.

Assume that

    1. (a)

$\mathrm{Re}(T_{1})$ is the sum of a positive coercive operator and self-
adjoint compact operator.

    2. (b)

There exists a compact operator $R\mathrel{\mathop{\mathchar
58\relax}}\mathbb{X}_{1}\to\mathbb{X}_{2}$ such that $H_{2}=RH_{1}$.

Then we have that
$\mathrm{Re}(\mathcal{F}_{2})\leq_{\mathrm{fin}}\mathrm{Re}(\mathcal{F}_{1}).$

  2. 2.

Assume that

    1. (a)

$\mathrm{Re}(T_{2})$ is the sum of a positive coercive operator and self-
adjoint compact operator.

    2. (b)

$\mathrm{Ran}(H_{1}^{*})\cap\mathrm{Ran}(H_{2}^{*})=\\{0\\}$ where
Dim$\big{(}\mathrm{Ran}(H_{2}^{*})\big{)}=\infty$.

Then we have that
$\mathrm{Re}(\mathcal{F}_{2})\nleq_{\mathrm{fin}}\mathrm{Re}(\mathcal{F}_{1}).$

We note that the real part of a compact operator
$\mathcal{F}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{Y}\rightarrow\mathbb{Y}$ is given by

$\mathrm{Re}(\mathcal{F})=\frac{1}{2}\left(\mathcal{F}+\mathcal{F}^{*}\right).$

Notice, that this is always a self-adjoint compact operator. To apply the
above result to our problem we first notice we already have the needed
factorization $F=-H^{*}TH$. The next thing to accomplish is to show that
$\pm\mathrm{Re}(T)$ is the sum of a positive coercive operator and self-
adjoint compact operator where $T$ is defined by (14).

Going forward, we want to prove that we can apply the monotonicity method
result under the analytical assumption that

$\displaystyle\text{Re}(A)-I>0\quad\text{uniformly in }D$

or for some constant $\alpha>0$

$\displaystyle I-\text{Re}(A)-\alpha\text{Im}(A)>0\quad\text{uniformly in
}D\,\,\textrm{ and}\quad\text{Re}(A)-\frac{1}{\alpha}\text{Im}(A)\geq
0\quad\text{in }D.$

The notation $M>0$ ($M\geq 0)$ denotes that the matrix $M$ is a positive
definite (semi-positive definite) matrix. With this in mind we will prove that
the operator $T$ has the desired properties to apply the above result.

Recall, the definition of the operator $T\mathrel{\mathop{\mathchar
58\relax}}H^{1}(D)\longrightarrow H^{1}(D)$ for any $v\in H^{1}(D)$

$\displaystyle(Tv,\phi)_{H^{1}(D)}=\int_{D}\nabla\overline{\phi}\cdot(A-I)\nabla(u+v)-k^{2}(n-1)(u+v)\overline{\phi}\,\text{d}x-\int_{\partial
D}\overline{\phi}\eta(u+v)\,\text{d}s$

where $u\in H^{1}_{loc}(\mathbb{R}^{d})$ satisfies the auxiliary problem
(7)–(8) with $v\in H^{1}(D)$ given. Therefore, we have that for any $v_{j}$
there exists a unique $u_{j}$ that satisfies the variational formulation (13)
for $j$=1 and 2. By definition we have that

$\displaystyle(Tv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{v_{2}}\,\text{d}x-\int_{\partial
D}\overline{v_{2}}\eta(u_{1}+v_{1})\,\text{d}s$
$\displaystyle=\int_{D}\nabla\overline{(u_{2}+v_{2})}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}x$
$\displaystyle\hskip
72.26999pt-\int_{D}\nabla\overline{u_{2}}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{u_{2}}\,\text{d}x$
$\displaystyle\hskip 72.26999pt-\int_{\partial
D}\overline{(u_{2}+v_{2})}\eta(u_{1}+v_{1})\,\text{d}s+\int_{\partial
D}\overline{u_{2}}\eta(u_{1}+v_{1})\,\text{d}s$

Now using the variational form (13) corresponding to $(u_{1},v_{1})$ with
$\phi=u_{2}$ we have

$\displaystyle-\int_{B_{R}}\nabla
u_{1}\cdot\nabla\overline{u_{2}}-k^{2}u_{1}\overline{u_{2}}\,\text{d}x$
$\displaystyle+\int_{\partial
B_{R}}\overline{u_{2}}\Lambda_{k}u_{1}\,\text{d}s$
$\displaystyle=\int_{D}\nabla\overline{u_{2}}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{u_{2}}\,\text{d}x-\int_{\partial
D}\overline{u_{2}}\eta(u_{1}+v_{1})\,\text{d}s$

Plugging this into the expression for $T$ we get

$\displaystyle(Tv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{(u_{2}+v_{2})}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}x$
$\displaystyle-\int_{\partial
D}\eta(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}s+\int_{B_{R}}\nabla
u_{1}\cdot\nabla\overline{u_{2}}-k^{2}u_{1}\overline{u_{2}}\,\text{d}x-\int_{\partial
B_{R}}\overline{u_{2}}\Lambda_{k}u_{1}\,\text{d}s$

With these calculations, we show that our operator $T$ has the desired
decomposition for us to use the monotonicity method given by Theorem 5.1.

######  Theorem 5.2.

Let $T\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\longrightarrow H^{1}(D)$
be as defined in (14).

  1. 1.

If $\mathrm{Re}(A)-I>0$ uniformly in $D$, then $\mathrm{Re}(T)$ is the sum of
a positive coercive operator and self-adjoint compact operator.

  2. 2.

If $I-\mathrm{Re}(A)-\alpha\mathrm{Im}(A)>0$ uniformly in $D$ and
$\mathrm{Re}(A)-\frac{1}{\alpha}\mathrm{Im}(A)\geq 0$ in $D$, then
$-\mathrm{Re}(T)$ is the sum of a positive coercive operator and self-adjoint
compact operator.

###### Proof.

(1) Using Riesz Representation Theorem on the above calculations for the
operator $T$, we can define the following bounded linear operators $S$,
$K\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\longrightarrow H^{1}(D)$ as
follows

$\displaystyle(Sv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{(u_{2}+v_{2})}\cdot(A-I)\nabla(u_{1}+v_{1})+(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}x$
$\displaystyle\hskip 36.135pt+\int_{B_{R}}\nabla
u_{1}\cdot\nabla\overline{u_{2}}+u_{1}\overline{u_{2}}\,\text{d}x-\int_{\partial
B_{R}}\overline{u_{2}}\Lambda_{k}u_{1}\,\text{d}s$ (15)

and

$\displaystyle(Kv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=-\int_{D}k^{2}(n-1)(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}x-\int_{D}(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}x$
$\displaystyle\hskip 36.135pt-\int_{\partial
D}\eta(u_{1}+v_{1})\overline{(u_{2}+v_{2})}\,\text{d}s+\int_{B_{R}}(k^{2}+1)u_{1}\overline{u_{2}}\,\text{d}x.$
(16)

With this it is clear that $T=S+K$ which implies that
$\text{Re}(T)=\text{Re}(S)+\text{Re}(K)$. Notice, by virtue of the compact
embeddings $H^{1}(D)$ into $L^{2}(D)$, $H^{1}(B_{R})$ into $L^{2}(B_{R})$, and
$H^{1/2}(\partial D)$ into $L^{2}(\partial D)$(see for e.g. [32]), we have
that $K$, and thus $\text{Re}(K)$ is also compact operator. It is also clear
to see that $\text{Re}(K)$ is a self-adjoint operator.
To complete the proof, we will show that $\text{Re}(S)$ is a positive coercive
operator provided that $\text{Re}(A)-I>0$ uniformly in $D$. To this end,
notice that

$\displaystyle(\text{Re}(S)v_{1},v_{1})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{(u_{1}+v_{1})}\cdot(\text{Re}(A)-I)\nabla(u_{1}+v_{1})+|u_{1}+v_{1}|^{2}\,\text{d}x$
$\displaystyle\hskip 36.135pt+\int_{B_{R}}|\nabla
u_{1}|^{2}+|u_{1}|^{2}\,\text{d}x-\int_{\partial
B_{R}}\overline{u_{1}}\text{Re}(\Lambda_{k})u_{1}\,\text{d}s$ (17)

where $\text{Re}(\Lambda_{k})$ is a non-positive operator (see for e.g. [10]).
We want to show that $\text{Re}(S)$ is positive coercive i.e. there is a
$\beta>0$ independent of $v_{1}\in H^{1}(D)$ such that

$(\text{Re}(S)v_{1},v_{1})_{H^{1}(D)}\geq\beta||v_{1}||^{2}_{H^{1}(D)}.$

To this end, we proceed by considering a contradiction argument. We assume
$\text{Re}(S)$ has no such constant, then there exists a sequence $v^{n}\in
H^{1}(D)$ with corresponding $u^{n}\in H^{1}_{loc}(\mathbb{R}^{d})$ satisfying
(7)–(8) such that

$||v^{n}||_{H^{1}(D)}=1\quad\text{and}\quad(\text{Re}(S)v^{n},v^{n})_{H^{1}(D)}\leq\frac{1}{n}.$

Since $\text{Re}(A)-I$ is positive definite and $\text{Re}(\Lambda_{k})$ is
non-positive, as $n\rightarrow\infty$, we have

$(\text{Re}(S)v^{n},v^{n})_{H^{1}(D)}\rightarrow 0\quad\text{implying}\quad
u^{n}\rightarrow 0\quad\text{and}\quad v^{n}\rightarrow 0\quad\text{in
}H^{1}(D)$

This contradicts the normalization of $v^{n}$ in $H^{1}(D)$. Therefore,
$\text{Re}(S)$ is a positive coercive operator. This proves the claim for
$\text{Re}(A)-I>0$ uniformly in $D$.

(2) Since we have new assumptions on the coefficients, we will have to derive
a new variational expression of the operator $T$. To this end, recall that by
(14) for a given $v_{j}\in H^{1}(D)$ with corresponding $u_{j}\in
H^{1}_{loc}(\mathbb{R}^{d})$ satisfying (7)–(8) we have that

$\displaystyle(Tv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(A-I)\nabla(u_{1}+v_{1})-k^{2}(n-1)(u_{1}+v_{1})\overline{v_{2}}\,\text{d}x-\int_{\partial
D}\eta(u_{1}+v_{1})\overline{v_{2}}\,\text{d}s$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(A-I)\nabla
v_{1}-k^{2}(n-1)v_{1}\overline{v_{2}}\,\text{d}x-\int_{\partial
D}\eta(u_{1}+v_{1})\overline{v_{2}}\,\text{d}s$ $\displaystyle\hskip
108.405pt+\int_{D}\nabla\overline{v_{2}}\cdot(A-I)\nabla
u_{1}-k^{2}(n-1)u_{1}\overline{v_{2}}\,\text{d}x.$

In an effort to rewrite the above expression for $T$, we use the variational
form (9) with respect to $(v_{2},u_{2})$ and let $\phi=u_{1}$, obtaining

$\displaystyle-\int_{B_{R}}A\nabla u_{2}\cdot\nabla\overline{u_{1}}$
$\displaystyle-k^{2}nu_{2}\overline{u_{1}}\,\text{d}x+\int_{\partial
B_{R}}\overline{u_{1}}\Lambda_{k}u_{2}\,\text{d}s$
$\displaystyle=\int_{D}\nabla\overline{u_{1}}\cdot(A-I)\nabla
v_{2}-k^{2}(n-1)v_{2}\overline{u_{1}}\,\text{d}x-\int_{\partial
D}\eta(u_{2}+v_{2})\overline{u_{1}}\,\text{d}s$

By taking the conjugate of the above equation, we can rewrite the definition
of $-T$ we obtain

$\displaystyle(-Tv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(I-A)\nabla
v_{1}+k^{2}(n-1)v_{1}\overline{v_{2}}\,\text{d}x+\int_{\partial D}\eta
v_{1}\overline{v_{2}}\,\text{d}s$ $\displaystyle\hskip
36.135pt-\int_{D}\nabla\overline{v_{2}}\cdot(A-\overline{A})\nabla
u_{1}-k^{2}(n-\overline{n})u_{1}\overline{v_{2}}\,\text{d}x+\int_{\partial
D}(\eta-\overline{\eta})u_{1}\overline{v_{2}}\,\text{d}s$ $\displaystyle\hskip
36.135pt+\int_{B_{R}}\nabla\overline{u_{2}}\cdot\overline{A}\nabla
u_{1}-k^{2}\overline{n}u_{1}\overline{u_{2}}\,\text{d}x-\int_{\partial
D}\overline{\eta}u_{1}\overline{u_{2}}\,\text{d}s-\int_{\partial
B_{R}}u_{1}\overline{\Lambda_{k}u_{2}}\,\text{d}s.$

Now, we compute
$-(\text{Re}(T)v_{1},v_{2})_{H^{1}(D)}=-\frac{1}{2}\big{(}(T+T^{*})v_{1},v_{2}\big{)}_{H^{1}(D)}$
we can determine

$(T^{*}v_{1},v_{2})_{H^{1}(D)}=\overline{(Tv_{2},v_{1})}_{H^{1}(D)}$

by the above calculations to obtain

$\displaystyle(-\text{Re}(T)v_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(I-\text{Re}(A))\nabla
v_{1}+k^{2}(\text{Re}(n)-1)v_{1}\overline{v_{2}}\,\text{d}x+\int_{\partial
D}\text{Re}(\eta)v_{1}\overline{v_{2}}\,\text{d}s$
$\displaystyle+\int_{B_{R}}\nabla\overline{u_{2}}\cdot\text{Re}(A)\nabla
u_{1}-k^{2}\text{Re}(n)u_{1}\overline{u_{2}}\,\text{d}x-\int_{\partial
D}\text{Re}(\eta)u_{1}\overline{u_{2}}\,\text{d}s$
$\displaystyle-\frac{1}{2}\int_{\partial
B_{R}}\overline{u_{2}}\Lambda_{k}u_{1}+u_{1}\overline{\Lambda_{k}u_{2}}\,\text{d}s$
$\displaystyle-\text{i}\int_{D}\nabla\overline{v_{2}}\cdot\text{Im}(A)\nabla
u_{1}-k^{2}\text{Im}(n)u_{1}\overline{v_{2}}\,\text{d}x+\text{i}\int_{\partial
D}\text{Im}(\eta)u_{1}\overline{v_{2}}\,\text{d}s$
$\displaystyle+\text{i}\int_{D}\nabla\overline{u_{2}}\cdot\text{Im}(A)\nabla
v_{1}-k^{2}\text{Im}(n)v_{1}\overline{u_{2}}\,\text{d}x-\text{i}\int_{\partial
D}\text{Im}(\eta)v_{1}\overline{u_{2}}\,\text{d}s.$

Using Riesz Representation Theorem on the definition of the operator
$-\text{Re}(T)$, we can define the following bounded linear operators
$S,K\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\longrightarrow H^{1}(D)$ as

$\displaystyle(Sv_{1},v_{2})_{H^{1}(D)}$
$\displaystyle=\int_{D}\nabla\overline{v_{2}}\cdot(I-\text{Re}(A))\nabla
v_{1}+v_{1}\overline{v_{2}}\,\text{d}x$
$\displaystyle+\int_{B_{R}}\nabla\overline{u_{2}}\cdot\text{Re}(A)\nabla
u_{1}\,\text{d}x-\frac{1}{2}\int_{\partial
B_{R}}\overline{u_{2}}\Lambda_{k}u_{1}+u_{1}\overline{\Lambda_{k}u_{2}}\,\text{d}s$
$\displaystyle-\text{i}\int_{D}\nabla\overline{v_{2}}\cdot\text{Im}(A)\nabla
u_{1}\,\text{d}x+\text{i}\int_{D}\nabla\overline{u_{2}}\cdot\text{Im}(A)\nabla
v_{1}\,\text{d}x$ (18)

and

$(Kv_{1},v_{2})_{H^{1}(D)}=(-\text{Re}(T)v_{1},v_{2})_{H^{1}(D)}-(Sv_{1},v_{2})_{H^{1}(D)}.$

Notice that $(Kv_{1},v_{2})_{H^{1}(D)}$ contains only $L^{2}$ terms in $D$,
$B_{R}$, and on $\partial D$, therefore the compactness of $K$ follows
similarly to the other case. It is also clear by the variational formulation
that $K$ is a self-adjoint operator.

Now we wish to show that $S$ is a positive coercive operator. To this end,
first we use Young’s Inequality to get the estimate

$\Big{|}\int_{D}\nabla\overline{v}\cdot\text{Im}(A)\nabla
u\,\text{d}x\Big{|}\leq\frac{\alpha}{2}(\text{Im}(A)\nabla v,\nabla
v)_{H^{1}(D)}+\frac{1}{2\alpha}(\text{Im}(A)\nabla u,\nabla u)_{H^{1}(D)}$

for any $\alpha>0$. Therefore we have

$\displaystyle(Sv_{1},v_{1})_{H^{1}(D)}$
$\displaystyle\geq\Big{(}\big{[}I-\text{Re}(A)-\alpha\text{Im}(A)\big{]}\nabla
v_{1},\nabla v_{1}\Big{)}_{L^{2}(D)}+(v_{1},v_{1})_{L^{2}(D)}$
$\displaystyle+\Big{(}\big{[}\text{Re}(A)-\frac{1}{\alpha}\text{Im}(A)\big{]}\nabla
u_{1},\nabla u_{1}\Big{)}_{H^{1}(D)}-\int_{\partial
B_{R}}\overline{u_{1}}\text{Re}(\Lambda_{k})u_{1}\,\text{d}s$
$\displaystyle\geq C||v_{1}||^{2}_{H^{1}(D)}$

meaning that we have completed this part of the proof provided $\alpha>0$ is a
constant such that $I-\text{Re}(A)-\alpha\text{Im}(A)>0$ uniformly in $D$ and
$\text{Re}(A)-\frac{1}{\alpha}\text{Im}(A)\geq 0$ in $D$. ∎

Now, we turn attention to satisfying Theorem 5.1. To this end, we need to find
another operator with a similar factorization to the far–field operator $F$.
With this in mind, being motivated by the work in [20] we now define
$B\mathrel{\mathop{\mathchar 58\relax}}=B(z,\epsilon)$ to be the ball of
radius $\epsilon>0$ centered at the point $z\in\mathbb{R}^{d}$, and the
operator

$H_{B}\mathrel{\mathop{\mathchar
58\relax}}L^{2}(\mathbb{S}^{d-1})\longrightarrow L^{2}(B)\quad\text{ given by
}\quad(H_{B}g)(x)=\int_{\mathbb{S}^{d-1}}\text{e}^{\text{i}kx\cdot\hat{y}}g(\hat{y})\,\text{d}s(\hat{y})\Big{|}_{B}.$

Notice, that when $B(z,\epsilon)\subset D$ the operator $H_{B}$ is the
restriction of the operator $H$ defined in (11) on the ball
$B(z,\epsilon)\subset D$ in the $L^{2}(B)$. This implies that $H_{B}=R_{B}H$
where

$\displaystyle R_{B}\mathrel{\mathop{\mathchar 58\relax}}H^{1}(D)\rightarrow
L^{2}(B)\quad R_{B}f=f\big{|}_{B}\quad\text{provided that $B\subset D$}$ (19)

is the restriction operator onto $B(z,\epsilon)$. Again, by appealing to the
compact embedding of $H^{1}$ into $L^{2}$ we have that $R_{B}$ is a compact
operator.

With this in mind, we now have just about all we need to state the main
result. The last piece of the puzzle to apply Theorem 5.1 to our scattering
problem is to consider the intersection of the range of $H^{*}$ and
$H^{*}_{B}$. In order to apply the result we need to show that the
corresponding ranges are disjoint provided that the ‘sampling ball’
$B=B(z,\epsilon)$ and $D$ are disjoint. This will ultimately allow us to test
whether or not the ‘sampling point’ $z\in\mathbb{R}^{d}$ is in the support of
the scatterer $D$.

######  Theorem 5.3.

If for a given sampling point $z\in\mathbb{R}^{d}$ and radius $\epsilon>0$ we
have that $B(z,\epsilon)\cap D=\emptyset$ then
$\mathrm{Ran}(H^{*}_{B})\cap\mathrm{Ran}(H^{*})=\\{0\\}$.

###### Proof.

To begin, we will assume $B(z,\epsilon)\cap D=\emptyset$, and let
$g\in\mathrm{Ran}(H^{*}_{B})\cap\mathrm{Ran}(H^{*})$. By Theorem 4.1 we have
that for some $f_{1}\in H^{1}(D)$ by (12) there exists $u_{D}\in
H^{1}_{loc}(\mathbb{R}^{d})$ such that $g=H^{*}f_{1}=u_{D}^{\infty}$ that
satisfy

$\Delta u_{D}+k^{2}u_{D}=0\quad\text{in }\mathbb{R}^{d}\setminus\overline{D}.$

Similarly, we have that for some $f_{2}\in L^{2}(B)$ just as in [20] there
exists $u_{B}\in H^{1}_{loc}(\mathbb{R}^{d})$ such that
$g=H_{B}^{*}f_{2}=u_{B}^{\infty}$ that satisfy

$\Delta u_{B}+k^{2}u_{B}=0\quad\text{in }\mathbb{R}^{d}\setminus\overline{B}.$

Here both $u_{D}$ and $u_{B}$ satisfy the radiation condition (3). By
Rellich’s Lemma, we have that $u_{B}=u_{D}$ in
$\mathbb{R}^{d}\setminus\big{(}\overline{D}\cup\overline{B}\big{)}$. With
this, we can now define $w\in H^{1}_{loc}(\mathbb{R}^{d})$ as

$w=\begin{cases}u_{D}&\quad\text{in }D,\\\ u_{B}&\quad\text{in
}\mathbb{R}^{d}\setminus\overline{D}.\end{cases}$

Notice, that $w$ is a radiating solution to Helmholtz equation in all of
$\mathbb{R}^{d}$, and thus $w=0$ in $\mathbb{R}^{d}$. This implies that
$0=w^{\infty}=u_{D}^{\infty}$, proving the claim since $g=u_{D}^{\infty}=0$. ∎

Notice, that by the above results we can apply Theorem 5.1 to our inverse
problem. Indeed, let us consider the case when $\text{Re}(A)-I>0$. Then we
would have that in Theorem 5.1 we consider

$\mathcal{F}_{1}=-F=H^{*}TH\quad\text{and}\quad\mathcal{F}_{2}=H^{*}_{B}H_{B}.$

By Theorem 5.2 gives that $\text{Re}(T)$ is the sum of a positive coercive
operator and a self-adjoint compact operator. Recall, that $H_{B}=R_{B}H$ and
by the definition given in (19) we have that $R_{B}$ is compact. Therefore, by
appealing to (1) of Theorems 5.1 and 5.2 we have that provided the sampling
ball $B(z,\epsilon)\subset D$ then

$H^{*}_{B}H_{B}\leq_{\mathrm{fin}}-\text{Re}(F).$

Now, by appealing to (2) of Theorem 5.1 as well as Theorem 5.3 we have that
provided the sampling ball $B(z,\epsilon)\cap D=\emptyset$ then

$H^{*}_{B}H_{B}\nleq_{\mathrm{fin}}-\text{Re}(F).$

Here we use that fact that $H_{B}$ is injective which implies that $H^{*}_{B}$
has a dense range i.e. the range is infinite dimensional. Similarly for the
case when $I-\text{Re}(A)-\alpha\text{Im}(A)>0$ uniformly in $D$ and
$\text{Re}(A)-\frac{1}{\alpha}\text{Im}(A)\geq 0$ in $D$. For this case we
note that in Theorem 5.1 we let

$\mathcal{F}_{1}=F=H^{*}(-T)H\quad\text{and}\quad\mathcal{F}_{2}=H^{*}_{B}H_{B}.$

This gives the main result of the paper that derives a monotonicity based
reconstruction method to recover the unknown scatterer.

######  Theorem 5.4.

Let $D\subset\mathbb{R}^{d}$ be a simply connected open set with a Lipschitz
boundary and $B=B(z,\epsilon)$ be the ball of radius $\epsilon$ centered at
$z\in\mathbb{R}^{d}$, we have the following:

  1. 1.

If $\mathrm{Re}(A)-I>0$, then

$B(z,\epsilon)\subset D\quad\iff\quad
H_{B}^{*}H_{B}\leq_{\mathrm{fin}}-\mathrm{Re}(F).$

  2. 2.

If $I-\mathrm{Re}(A)-\alpha\mathrm{Im}(A)>0$ and
$\mathrm{Re}(A)-\frac{1}{\alpha}\mathrm{Im}(A)\geq 0$ for some $\alpha>0$,
then

$B(z,\epsilon)\subset D\quad\iff\quad
H_{B}^{*}H_{B}\leq_{\mathrm{fin}}\mathrm{Re}(F).$

Theorem 5.4 implies that the eigenvalues of the self-adjoint operator
$\pm\mathrm{Re}(F)-H_{B}^{*}H_{B}$ can be used to recover $D$. In the next
section, we will provide some numerical examples of this to recover circular
regions in two dimensions with respect to noisy data.

##  6 Numerical Results

In previous sections, we have proved that we can use the monotonicity method
to perform reconstructions of the scatterer $D$. Now, we will provide
numerical examples to validate our theoretical results. It is worth noting
that the monotonicity method given in Theorem 5.4 is a qualitative
reconstruction method that is valid for all wave number $k$. Unlike the linear
sampling and factorization methods that have to avoid wave numbers that are
so–called transmission eigenvalues. Notice, that even though the monotonicity
method is theoretically valid for all wave numbers, one does need to have the
a priori information on the matrix–valued coefficient to apply this method.

In order to apply the monotonicity method to recover $D$, we need to compute
$H^{*}_{B}H_{B}$ where the sampling ball $B=B(z,\epsilon)$. To this end, in
the following subsection we derive an integral representation for the operator
$H^{*}_{B}H_{B}$. We will do the calculation in both 2 and 3 dimensions.
Recall, that the operator

$H^{*}_{B}H_{B}\mathrel{\mathop{\mathchar
58\relax}}L^{2}(\mathbb{S}^{d-1})\longrightarrow L^{2}(\mathbb{S}^{d-1})$

so we will write it as an integral operator with an explicit kernel function.

###  6.1 Computing the Integral Representation of $H^{*}_{B}H_{B}$

In this section, we provide the calculations for writing $H^{*}_{B}H_{B}$ as
an integral operator mapping $L^{2}(\mathbb{S}^{d-1})$ into itself. Note, that
here we consider $H_{B}$ as an operator mapping into $H^{1}(B)$ which implies
that $H^{*}_{B}$ is given in Theorem 4.1 and is associated with the
variational formulation given in (12). This is due to the fact that in our
investigation the reconstructions are more accurate for this case. Recall,
that $H^{*}_{B}f=-w^{\infty}$ where $w$ solves (12) for any given $f\in
H^{1}(D)$. To this end, we will need to determine the strong form of the PDE
associated with the equivalent variational formulation given in (12).

By appealing to Green’s Identities we have that $w\in
H^{1}_{loc}(\mathbb{R}^{d})$ satisfying the variational formulation (12) will
be the solution to

$\displaystyle\Delta w+k^{2}w=f-\Delta f\quad\text{ in
}B=B(z_{j},\epsilon)\quad$ $\displaystyle\text{ and }\quad\Delta
w+k^{2}w=0\quad\text{ in }\mathbb{R}^{d}\setminus\overline{B}$ $\displaystyle
w^{+}-w^{-}=0\quad$ $\displaystyle\text{ and
}\quad\partial_{\nu}w^{+}-\partial_{\nu}w^{-}=\partial_{\nu}f\quad\text{ on
}\partial B$

for any $f\in H^{1}(B)$ such that $\Delta f\in L^{2}(B)$, closed with the
Sommerfeld radiation condition. From this, we get the expression for
$H^{*}_{B}$, given by

$\displaystyle(H^{*}_{B})f(\hat{x})=-w^{\infty}(\hat{x})=\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(f-\Delta f)\,\text{d}z+\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\partial_{\nu}f\,\text{d}s(z)$

just as in [25]. Now, using the fact that $H_{B}g$ is a smooth solution to the
Helmholtz equation in $B$ we have that

$\displaystyle(H^{*}_{B}H_{B})g(\hat{x})$
$\displaystyle=\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot z}(H_{B}g-\Delta
H_{B}g)\,\text{d}z+\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\partial_{\nu}H_{B}g\,\text{d}s(z)$
$\displaystyle=\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(k^{2}+1)(H_{B}g)\,\text{d}z+\int_{\partial
B}\text{e}^{-\text{i}k\hat{x}\cdot z}\partial_{\nu}(H_{B}g)\,\text{d}s(z).$

In order to proceed, we will need to evaluate the volume and surface integral
above.

With this, we will start by considering the volume above to obtain that

$\displaystyle\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(k^{2}+1)(H_{B}g)\,\text{d}z$
$\displaystyle=\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(k^{2}+1)\Bigg{[}\int_{\mathbb{S}^{d-1}}\text{e}^{\text{i}k\hat{y}\cdot
z}g(\hat{y})\,\text{d}s(\hat{y})\Bigg{]}\,\text{d}z$
$\displaystyle=(k^{2}+1)\int_{\mathbb{S}^{d-1}}\Bigg{[}\int_{B}\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\,\text{d}z\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$

Using polar/spherical coordinates to represent the inner integral over
$B=B(z_{j},\epsilon)$, we obtain that

$\displaystyle\int_{B}\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\,\text{d}z$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\int_{B(0,\epsilon)}\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\,\text{d}z$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\int_{\mathbb{S}^{d-1}}\int_{0}^{\epsilon}\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}r^{d-1}\,\text{d}r\,\text{d}s(\hat{z})$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\int_{0}^{\epsilon}\Bigg{[}\int_{\mathbb{S}^{d-1}}\text{e}^{-\text{i}k\hat{z}\cdot
r(\hat{x}-\hat{y})}\,\text{d}s(\hat{z})\Bigg{]}r^{d-1}\,\text{d}r$

where $z=r\hat{z}$ and $r=|z|$. Recall, the Funk–Hecke integral identities(see
for e.g. [24]) that are given as follows

$\displaystyle\int_{\mathbb{S}^{d-1}}\text{e}^{-\text{i}k\hat{z}\cdot
r(x-y)}\,\text{d}s(\hat{z})=\begin{cases}2\pi J_{0}(kr|x-y|)&\quad\text{for
}d=2,\\\ \\\ 4\pi j_{0}(kr|x-y|)&\quad\text{for }d=3\end{cases}$

where $J_{m}$ are Bessel functions of the first kind of order $m$ and $j_{m}$
are Spherical Bessel functions of the first kind of order $m$. Using well
known recurrence relationships for Bessel functions we can calculate the
following anti–derivates

$\int tJ_{0}(t)\,\text{d}t=tJ_{1}(t)\quad\text{and}\quad\int
t^{2}j_{0}(t)\,\text{d}t=t^{2}j_{1}(t).$

Therefore, for $d=2$ we have

$\displaystyle\int_{B}\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\,\text{d}z$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\int_{0}^{\epsilon}J_{0}(kr|\hat{x}-\hat{y}|)r\,\text{d}r$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\Bigg{(}\frac{2\pi\epsilon}{k|\hat{x}-\hat{y}|}J_{1}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi\epsilon^{2}\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}$

where we have used the recurrence relationship
$(2/t)J_{1}(t)=J_{0}(t)+J_{2}(t)$ along with a substitution in the above
calculations. We obtain that the volume integral is given by

$\displaystyle\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(k^{2}+1)(H_{B}g)\,\text{d}z$ $\displaystyle\hskip
72.26999pt=(k^{2}+1)\int_{\mathbb{S}^{1}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi\epsilon^{2}\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y}).$

Similar calculation for $d=3$, we obtain that

$\displaystyle\int_{B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(k^{2}+1)(H_{B}g)\,\text{d}z$ $\displaystyle\hskip
72.26999pt=(k^{2}+1)\int_{\mathbb{S}^{2}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\frac{4}{3}\pi\epsilon^{3}\Bigg{(}j_{0}(k\epsilon|\hat{x}-\hat{y}|)+j_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$

by appealing to the recurrence relationship $(3/t)j_{1}(t)=j_{0}(t)+j_{2}(t)$
along with a substitution in the above calculations.

Now, we turn our attention to the surface integral to continue our derivation
of an integral representation of $H^{*}_{B}H_{B}$. Notice that

$\displaystyle\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(\partial_{\nu}H_{B}g)\,\text{d}s(z)$ $\displaystyle=\int_{\partial
B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\Bigg{[}\int_{\mathbb{S}^{d-1}}\partial_{\nu}\text{e}^{\text{i}kz\cdot\hat{y}}g(\hat{y})\,\text{d}s(\hat{y})\Bigg{]}\,\text{d}s(z)$
$\displaystyle=\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\Bigg{[}\int_{\mathbb{S}^{d-1}}\partial_{r}\text{e}^{\text{i}kr\hat{z}\cdot\hat{y}}g(\hat{y})\,\text{d}s(\hat{y})\Bigg{]}\,\text{d}s(z)$
$\displaystyle=\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\Bigg{[}\int_{\mathbb{S}^{d-1}}(\text{i}k\hat{z}\cdot\hat{y})\text{e}^{\text{i}kr\hat{z}\cdot\hat{y}}g(\hat{y})\,\text{d}s(\hat{y})\Bigg{]}\,\text{d}s(z)$

Here we have used the fact that we are integrating over $\partial B$, which
would imply that $\partial_{\nu}=\partial_{r}$ when written in polar/spherical
coordinates. Continuing with the calculations we see that

$\displaystyle\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(\partial_{\nu}H_{B}g)\,\text{d}s(z)$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\Bigg{[}\int_{\partial
B}(\text{i}k\hat{z}\cdot\hat{y})\text{e}^{-\text{i}k{z}\cdot(\hat{x}-\hat{y})}\text{d}s(z)\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\int_{\partial
B(0,\epsilon)}(\text{i}k\hat{z}\cdot\hat{y})\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\text{d}s(z)\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\text{i}k\hat{y}\epsilon^{d-1}\cdot\int_{\mathbb{S}^{d-1}}\hat{z}\text{e}^{-\text{i}k\hat{z}\cdot\epsilon(\hat{x}-\hat{y})}\,\text{d}s(\hat{z})\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$

We now use the following Funk–Hecke integral identities(see for e.g. [24])

$\displaystyle\int_{\mathbb{S}^{d-1}}\hat{z}\text{e}^{-\text{i}k\hat{z}\cdot
r(x-y)}\,\text{d}s(\hat{z})=\begin{cases}2\pi\frac{(x-y)}{\text{i}|x-y|}J_{1}(kr|x-y|)&\quad\text{for
}d=2,\\\ \\\ 4\pi\frac{(x-y)}{\text{i}|x-y|}j_{1}(kr|x-y|)&\quad\text{for
}d=3\end{cases}$

to get the following expression

$\displaystyle\int_{\partial
B}(\text{i}k\hat{z}\cdot\hat{y})\text{e}^{-\text{i}kz\cdot(\hat{x}-\hat{y})}\,\text{d}s(z)$
$\displaystyle=2\pi\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\text{i}k\epsilon\hat{y}\cdot\frac{(\hat{x}-\hat{y})}{\text{i}|\hat{x}-\hat{y}|}J_{1}(k\epsilon|\hat{x}-\hat{y}|)$
$\displaystyle=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi
k^{2}\epsilon^{2}(\hat{x}\cdot\hat{y}-1)\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}$

for $d=2$. Notice, that we have used that
$\hat{y}\cdot(\hat{x}-\hat{y})=\hat{x}\cdot\hat{y}-1$ to now get the
expression

$\displaystyle\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}(\partial_{\nu}H_{{}_{B}}g)\,\text{d}s(z)$ $\displaystyle\hskip
36.135pt=\int_{\mathbb{S}^{1}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi
k^{2}\epsilon^{2}(\hat{x}\cdot\hat{y}-1)\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y}).$

Again, similar calculations for $d=3$, we have that

$\displaystyle\int_{\partial B}\text{e}^{-\text{i}k\hat{x}\cdot
z}\partial_{\nu}H_{B}g\,\text{d}s(z)$ $\displaystyle\hskip
36.135pt=\int_{\mathbb{S}^{2}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\frac{4}{3}\pi
k^{2}\epsilon^{3}(\hat{x}\cdot\hat{y}-1)\Bigg{(}j_{0}(k\epsilon|\hat{x}-\hat{y}|)+j_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y}).$

With this, we can get the integral representation for $H^{*}_{B}H_{B}$, given
by

$\displaystyle(H^{*}_{B}H_{B})g(\hat{x})=\int_{\mathbb{S}^{1}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi\epsilon^{2}(1+k^{2}\hat{x}\cdot\hat{y})\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$

for $d=2$. Similarly, by combining our calculations we have

$\displaystyle(H^{*}_{B}H_{B})g(\hat{x})=\int_{\mathbb{S}^{2}}\Bigg{[}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\frac{4}{3}\pi\epsilon^{3}(1+k^{2}\hat{x}\cdot\hat{y})\Bigg{(}j_{0}(k\epsilon|\hat{x}-\hat{y}|)+j_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}\Bigg{]}g(\hat{y})\,\text{d}s(\hat{y})$

for $d=3$.

Notice, that by our calculations we obtained that

$(H^{*}_{B}H_{B})g(\hat{x})=\int_{\mathbb{S}^{d-1}}h(\hat{x},\hat{y})g(\hat{y})\,\text{d}s(\hat{y})$

where the kernel function is given by

$h(\hat{x},\hat{y})=\begin{cases}\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi\epsilon^{2}(1+k^{2}\hat{x}\cdot\hat{y})\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}&\quad\text{for
}d=2,\\\ \\\
\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\frac{4}{3}\pi\epsilon^{3}(1+k^{2}\hat{x}\cdot\hat{y})\Bigg{(}j_{0}(k\epsilon|\hat{x}-\hat{y}|)+j_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}&\quad\text{for
}d=3\end{cases}$

where we see that $h(\hat{x},\hat{y})=\overline{h(\hat{y},\hat{x})}$. This is
to be expected due to the fact that $H^{*}_{B}H_{B}$ is a self-adjoint
operator. With this we have all we need to provide some numerical
reconstructions using the monotonicity result in Theorem 5.4.

###  6.2 Numerical Reconstruction

Now, we will provide some numerical examples for reconstruction of the
scatterer $D$ via the monotonicity method. To this end, we will use separation
of variables to obtain an integral representation for the far–field operator
$F$. For simplicity, we consider the direct scattering problem with constant
coefficients on $B(0,R)$, where assume $A=aI$, with constant $a\in\mathbb{C}$
along with constants $n$, and $\eta$. Therefore, we have that (1)–(2) is given
by

$\displaystyle\Delta u^{s}+k^{2}u^{s}=0\quad\text{in
}\mathbb{R}^{d}\setminus\overline{B(0,R)}\quad\text{ and }\quad\Delta
u+\frac{n}{a}k^{2}u=0\quad$ $\displaystyle\text{in }B(0,R)$
$\displaystyle(u^{s}+u^{i})^{+}=u^{-}\quad\text{and}\quad
a\partial_{r}u^{-}=\partial_{r}(u^{s}+u^{i})^{+}+\eta(u^{s}+u^{i})\quad$
$\displaystyle\text{on }\partial B(0,R)$

along with the radiation condition as $|x|\to\infty$ for the scattered field
$u^{s}$. Recall that we can express the incident field
$u^{i}=\text{e}^{\text{i}kx\cdot\hat{y}}$, using the Jacobi–Anger expansion is
given by

$u^{i}(r,\theta)=\sum_{p\in\mathbb{Z}}\text{i}^{p}J_{p}(kr)\text{e}^{\text{i}p(\theta-\phi)}$

where $x=r\big{(}\cos(\theta),\sin(\theta)\big{)}$ and
$\hat{y}=\big{(}\cos(\phi),\sin(\phi)\big{)}$.

To solve for the scattered field $u^{s}$, and the total field $u$ for this
problem we make the assumption that they can be written as

$u^{s}(r,\theta)=\sum_{p\in\mathbb{Z}}\text{i}^{p}\text{u}^{s}_{p}H_{p}^{(1)}(kr)\text{e}^{\text{i}p(\theta-\phi)}\quad\text{and}\quad
u(r,\theta)=\sum_{p\in\mathbb{Z}}\text{i}^{p}\text{u}_{p}J_{p}\Big{(}k\sqrt{\frac{n}{a}}r\Big{)}\text{e}^{\text{i}p(\theta-\phi)}$

where $H^{(1)}_{p}$ are Hankel functions of the first kind of order $p$. Here
$\text{u}^{s}_{p}$ and $\text{u}_{p}$ can be seen as the coefficients that are
to be determined by using the boundary conditions at $r=R$. Therefore, we have
that $\text{u}^{s}_{p}$ and $\text{u}_{p}$ satisfy the $2\times 2$ linear
system

$\displaystyle\begin{pmatrix}H_{p}^{(1)}(kR)&-J_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\\\
kH_{p}^{(1)^{\prime}}(kR)+\eta
H_{p}^{(1)}(kR)&-k\sqrt{na}J^{\prime}_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\end{pmatrix}\begin{pmatrix}\text{u}^{s}_{p}\\\
\text{u}_{p}\end{pmatrix}=\begin{pmatrix}-J_{p}(kR)\\\
-kJ^{\prime}_{p}(kR)-\eta J_{p}(kR)\end{pmatrix}$

for each $p\in\mathbb{Z}$. We can solve for the coefficients
$\text{u}^{s}_{p}$ by using Cramer’s Rule, we have

$\displaystyle\text{u}^{s}_{p}=\frac{\text{Det}(M_{x})}{\text{Det}(M)}\quad\text{
with matrices }\quad$ $\displaystyle
M=\begin{pmatrix}H_{p}^{(1)}(kR)&-J_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\\\
kH_{p}^{(1)^{\prime}}(kR)+\eta
H_{p}^{(1)}(kR)&-k\sqrt{na}J^{\prime}_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\end{pmatrix}$
and $\displaystyle
M_{x}=\begin{pmatrix}-J_{p}(kR)&-J_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\\\
-kJ^{\prime}_{p}(kR)-\eta
J_{p}(kR)&-k\sqrt{na}J^{\prime}_{p}\Big{(}k\sqrt{\frac{n}{a}}R\Big{)}\end{pmatrix}.$

With this computation of $\text{u}^{s}_{p}$, we get an expression for the
far–field pattern is given by

$\displaystyle
u^{\infty}(\hat{x},\hat{y})\approx\frac{4}{\text{i}}\sum_{|p|=0}^{15}\text{u}^{s}_{p}\text{e}^{\text{i}p(\theta-\phi)}\quad\text{
where }\quad\text{u}^{s}_{p}=\frac{\text{Det}(M_{x})}{\text{Det}(M)}.$

In all our numerical experiments we use the above approximation of the far-
field pattern.

In order to numerically compute the eigenvalues of
$\pm\text{Re}(F)-H^{*}_{B}H_{B}$, we approximate the integral representations
of $F$ and $H^{*}_{B}H_{B}$ with a 64–point Riemann sum. The $64\times 64$
discretization of the operators $F$ and $H^{*}_{B}H_{B}$ are computed as

$\displaystyle\textbf{F}=\left[\frac{2\pi}{64}u^{\infty}(\hat{x}_{i},\hat{y}_{j})\right]_{i,j=1}^{64}\quad\text{
and
}\quad\textbf{H}^{*}_{B}\textbf{H}_{B}=\left[\frac{2\pi}{64}h(\hat{x}_{i},\hat{y}_{j})\right]_{i,j=1}^{64}\quad\text{
respectively}.$

We recall, that the kernel function to define
$\textbf{H}^{*}_{B}\textbf{H}_{B}$ is given by

$h(\hat{x},\hat{y})=\text{e}^{-\text{i}kz_{j}\cdot(\hat{x}-\hat{y})}\pi\epsilon^{2}(1+k^{2}\hat{x}\cdot\hat{y})\Bigg{(}J_{0}(k\epsilon|\hat{x}-\hat{y}|)+J_{2}(k\epsilon|\hat{x}-\hat{y}|)\Bigg{)}.$

Here we discretize the $\mathbb{S}^{1}$ such that

$\hat{x}_{i}=\hat{y}_{i}=(\cos(\theta_{i}),\sin(\theta_{i}))\quad\text{ where
}\quad\theta_{i}=2\pi(i-1)/64\,\,\text{ for}\,\,i=1,\dots,64.$

In order to test the stability we will add random noise to the discretized
far–field operator $\mathbf{F}$ such that

$\mathbf{F}_{\delta}=\Big{[}\mathbf{F}_{i,j}(1+\delta\mathbf{E}_{i,j})\Big{]}^{64}_{i,j=1}\hskip
14.22636pt\text{where}\hskip 14.22636pt\|\mathbf{E}\|_{2}=1.$

Here, the matrix $\mathbf{E}\in\mathbb{C}^{64\times 64}$ is taken to have
random entries and $0<\delta<1$ is the relative noise level added to the data.

In our numerical computations, we consider recovering $D=B(0,R)$ for different
values of $R\in(0,1]$. Therefore, begin by discretizing a rectangular region
$\Omega=[-2,2]\times[-2,2]$ that completely encompasses our region $D=B(0,R)$,
and we fix a small radius $\epsilon$ of our ball $B$. In all our numerical
examples we pick the radius $\epsilon=10^{-4}$. According to the theory, as we
move our ball $B=B(z,\epsilon)$ along the grid $\Omega$, we should have that
there are more positive eigenvalues of the self-adjoint matrix

$\pm\text{Re}(\textbf{F}_{\delta})-\textbf{H}^{*}_{B}\textbf{H}_{B}\quad\text{
where
}\quad\text{Re}(\textbf{F}_{\delta})=\frac{1}{2}\big{(}\textbf{F}_{\delta}+\textbf{F}_{\delta}^{*}\big{)}$

provided that $z\in D$ and less positive eigenvalues if $z\not\in D$, roughly
speaking. Therefore, at each point in the region $\Omega$, we compute the
imaging function

$W(z)=\\#\Big{\\{}\lambda_{j}\in\text{eig}\big{(}\pm\mathrm{Re}(\textbf{F}_{\delta})-\textbf{H}^{*}_{B}\textbf{H}_{B}\big{)}\,\,\mathrel{\mathop{\mathchar
58\relax}}\,\,\lambda_{j}>0\Big{\\}}.$

This is just the cardinality of the set of positive eigenvalues. To evaluate
the imaging function we use the built in eig command in MATLAB. This gives us
a contour plot that will serve as our reconstruction for $D$. We remark that
in our experiments we do not attempt to regularize the imaging function(see
for e.g. [19]). This is due to the simplicity of the geometry we are
considering.

Example 1: For this reconstruction, we will provide two examples with minimal
amounts of noise in the data. We provide reconstructions with added $1\%$
noise i.e. $\delta=0.01$ in our calculations. The reconstructions are shown in
Figure 1.

Figure 1:  Reconstruction of domains with $1\%$ noise. Here $D=B(0,3/4)$ with
$k=10$, $a=2-\text{i}$, $n=1/2$, $\eta=2$ (left) and $D=B(0,1)$ with $k=2\pi$,
$a=1/4$, $n=2$, $\eta=2$ (right).

Example 2: We provide reconstruction examples for the same scatterers in the
previous example with more noise added added to the data i.e. $5\%$ noise. The
reconstructions are shown in Figure 2.

Figure 2:  Reconstruction of domains with $5\%$ noise. Here $D=B(0,3/4)$ with
$k=10$, $a=2-\text{i}$, $n=1/2$, $\eta=2$ (left) and $D=B(0,1)$ with $k=2\pi$,
$a=1/4$, $n=2$, $\eta=2$ (right).

Example 3: Now, we provide an example with added $10\%$ noise with slightly
different refractive indices $n$ in the scatterer. The reconstructions are
shown in Figure 3.

Figure 3:  Reconstruction of domains with $10\%$ noise. Here $D=B(0,3/4)$ with
$k=10$, $a=2-\text{i}$, $n=1/2$, $\eta=2$ (left) and $D=B(0,1)$ with $k=2\pi$,
$a=1/4$, $n=2$, $\eta=2$ (right).

Example 4: Lastly, we provide an example where there is no conductivity
parameter $\eta$. Also, we are inspired by [1] and consider finding the radius
of the region $D$ by plotting the number of positive eigenvalues where the
sampling point is fixed but the sampling ball $B(0,\epsilon)$ has radius
$\epsilon\in(0,2]$. As we see in Figure 4, the number of positive eigenvalues
of our operator decrease as the sampling radius increases. Here, the radius
$R=1/2$ for the scatterer $D$ and the number of positive eigenvalues decreases
almost monotonically when the radius of the sampling ball $\epsilon>1/2$. This
is consistent with the theoretical results.

Figure 4:  Reconstruction of domains with $2\%$ noise. Here $D=B(0,1/2)$ with
$k=10$, $a=4$, $n=1/2$, $\eta=0$ (left) and the plot of the number of positive
eigenvalues as a function of $\epsilon$ where the sampling ball is
$B(0,\epsilon)$ for $\epsilon\in(0,2]$ (right).

From our examples, we have faithful reconstructions for the circular regions
under either set of assumptions on the coefficients. We also see that the
method is resilient to added noise in these examples. With this we see that
the monotonicity method is applicable to reconstructing anisotropic scatterers
with or without a conductive boundary condition.

##  7 Conclusions

In this paper, we studied the monotonicity method for an anisotropic material
with a conductive boundary on an unbounded domain. We obtained the necessary
symmetric factorization of the far–field to apply the monotonicity method. Our
main contribution was adopting the theory of the monotonicity method for the
anisotropic material under two different assumptions on its physical
parameters. We then provided some numerical examples of using the method on
different circular domains in two dimensions with various levels of noise. We
see that this method creates faithful reconstructions that are resilient to
noise. With that, there is still room for more extensive numerical tests of
the monotonicity method studied in this paper. Another area to be investigated
is applying this analysis for near–field data [27] as well as scattering in
waveguides [2, 7] and by periodic structures [26] since the factorization of
the data operators has been studied.

Acknowledgments: The research of V. Hughes, I. Harris and H. Lee is partially
supported by the NSF DMS Grant 2107891.

## References

  * [1] A. Albicker and R. Griesmaier,  Monotonicity in inverse obstacle scattering on unbounded domains,  Inverse Problems, 36 (2020), 085014.
  * [2] T. Arens, R. Griesmaier and R. Zhang, Monotonicity-based shape reconstruction for an inverse scattering problem in a waveguide  Inverse Problems, 36 (2023), 075009.
  * [3] T. Brander, B. Harrach, M. Kar, and M. Salo,  Monotonicity and enclosure methods for the $p$-Laplace equation,  SIAM J. Appl. Math., 78 (2018), 742–758.
  * [4] O. Bondarenko and X. Liu,  The factorization method for inverse obstacle scattering with conductive boundary condition,  Inverse Problems, 29 (2013), 095021.
  * [5] O. Bondarenko and A. Kirsch,  The factorization method for inverse scattering by a penetrable anisotropic obstacle with conductive transmission conditions,  Inverse Problems, 32 (2016), 105011.
  * [6] O. Bondarenko, I. Harris, and A. Kleefeld,  The interior transmission eigenvalue problem for an inhomogeneous media with a conductive boundary,  Applicable Analysis, 96(1) (2017), 2–22.
  * [7] L. Borcea and S. Meng,  Factorization method versus migration imaging in a waveguide,  Inverse Problems, 35 (2019), 124006.
  * [8] V. Candiani, J. Darde, H. Garde, and N. Hyvonen,  Monotonicity-based reconstruction of extreme inclusions in electrical impedance tomography,  SIAM J. Math. Anal., 52 (2020), 6234–6259.
  * [9] F. Cakoni, D. Colton,  “A Qualitative Approach to Inverse Scattering Theory”,  Springer, New York, (2013).
  * [10] F. Cakoni, D. Colton, and H. Haddar,  “Inverse Scattering Theory and Transmission Eigenvalues”,  CBMS Series, SIAM 88, Philadelphia, (2016).
  * [11] F. Cakoni, D. Colton, and P. Monk,  “The linear Sampling Method in Inverse Electromagnetic Scattering”,  CBMS Series, SIAM Publications 80, (2011).
  * [12] F. Cakoni and A. Kirsch,  On the interior transmission eigenvalue problem,  Int. Jour. Comp. Sci. Math., 3 (2010), 142–167.
  * [13] F. Cakoni, H. Lee, P. Monk and Y. Zhang,  A spectral target signature for thin surfaces with higher order jump conditions,  Inverse Problems and Imaging, 16(6) (2022), 1473–1500.
  * [14] F. Cakoni, S. Meng and H. Haddar,  The factorization method for a cavity in an inhomogeneous medium,  Inverse Problems, 30 (2014), 045008.
  * [15] S. Chaabane, B. Charfi and H. Haddar,  Reconstruction of discontinuous parameters in a second order impedance boundary operator.  Inverse Problems 32(10) (2016), 105004.
  * [16] M. Cheney,  The linear sampling method and the MUSIC algorithm,  Inverse Problems, 17 (2001), 591–595.
  * [17] T. Daimon, T. Furuya and R. Saiin,  The monotonicity method for the inverse crack scattering problem,  Inverse Problems in Sci. and Eng., 28 (2020), 1570–1581.
  * [18] T. Furuya,  Remarks on the factorization and monotonicity method for inverse acoustic scatterings,  Inverse Problems, 37 (2021), 065006.
  * [19] H. Garde and S. Staboulis,  Convergence and regularization for monotonicity-based shape reconstruction in electrical impedance tomography,  Numer.Math., 135 (2017), 1221–1251.
  * [20] R. Griesmaier and B. Harrach,  Monotonicity in Inverse Medium Scattering on Unbounded Domains,  SIAM J. Appl. Math., 78 (2018), 2533–2557.
  * [21] F. Gylys-Colwell,  An inverse problem for the Helmholtz equation,  Inverse Problems 12 (1996), 139-156.
  * [22] B. Harrach and Y-H. Lin,  Monotonicity-based inversion of the fractional Schrodinger equation I. Positive potentials,  SIAM J. Math. Anal., 51 (2019), 3092–3111.
  * [23] B. Harrach and M. Ullrich,  Monotonicity-based shape reconstruction in electrical impedance tomography,  SIAM J. Math. Anal., 45 (2013), 3382–3403.
  * [24] I. Harris,  Direct methods for recovering sound soft scatterers from point source measurements,  Computation, 9(11) (2021), 120.
  * [25] I. Harris,  Regularized factorization method for a perturbed positive compact operator applied to inverse scattering,  Inverse Problems, 39 (2023), 115007.
  * [26] I. Harris, D.-L. Nguyen, J. Sands, and T. Truong,  On the inverse scattering from anisotropic periodic layers and transmission eigenvalues.  Applicable Analysis 101(8) (2022), 3065–3081.
  * [27] I. Harris and S. Rome,  Near field imaging of small isotropic and extended anisotropic scatterers,  Applicable Analysis, 96(10) (2017) 1713–1736.
  * [28] V. Hughes, I. Harris, and J. Sun,  The anisotropic transmission eigenvalue problem with a conductive boundary. (arXiv:2311.00526)
  * [29] A. Kirsch, The MUSIC algorithm and the factorization method in inverse scattering theory for inhomogeneous media,  Inverse Problems 18  (2002), 1025–1040.
  * [30] A. Kirsch and N. Grinberg,  “The Factorization Method for Inverse Problems”. 1st edition Oxford University Press, Oxford 2008.
  * [31] A. Kirsch and X. Liu,  The Factorization Method for Inverse Acoustic Scattering by a Penetrable Anisotropic Obstacle.  Math. Methods in the App.Sci. 37(8) (2014), 1159–1170.
  * [32] S. Salsa,  “Partial Differential Equations in Action From Modelling to Theory”,  Springer Italia, Milano, (2008).
  * [33] A. Tamburrino,  Monotonicity based imaging methods for elliptic and parabolic inverse problems  J. Inverse Ill-Posed Probl., 14 (2006), 633–642.
  * [34] J. Xiang and G. Yan,  The interior transmission eigenvalue problem for an anisotropic medium by a partially coated boundary,  Acta Mathematica Scientia 44 (2024), 339–354.