Reciprocity gap functional methods for potentials/sources with small volume
support for two elliptic equations

Govanni Granados and Isaac Harris
Department of Mathematics, Purdue University, West Lafayette, IN 47907
Email<EMAIL_ADDRESS>and<EMAIL_ADDRESS>

###### Abstract

In this paper, we consider inverse shape problems coming from diffuse optical
tomography and inverse scattering. In both problems, our goal is to
reconstruct small volume interior regions from measured data on the exterior
surface of an object. In order to achieve this, we will derive an asymptotic
expansion of the reciprocity gap functional associated with each problem. The
reciprocity gap functional takes in the measured Cauchy data on the exterior
surface of the object. In diffuse optical tomography, we prove that a MUSIC-
type algorithm can be used to recover the unknown subregions. This gives an
analytically rigorous and computationally simple method for recovering the
small volume regions. For the problem coming from inverse scattering, we
recover the subregions of interest via a direct sampling method. The direct
sampling method presented here allows use to accurately recover the small
volume region from one pair of Cauchy data. We also prove that the direct
sampling method is stable with respect to noisy data. Numerical examples will
be presented for both cases in two dimensions where the measurement surface is
the unit circle.

Keywords: Diffuse Optical Tomography $\cdot$ Inverse Scattering $\cdot$ MUSIC
Algorithm $\cdot$ Direct Sampling

MSC: 35J05, 35J25

##  1 Introduction

The two problems we consider in this paper are motivated by diffuse optical
tomography (DOT) and inverse scattering theory. In both problems, the goal is
to reconstruct interior subregions of small volume from known Cauchy data on
the boundary of the given bounded open set $\Omega$ in $\mathbb{R}^{2}$ or
$\mathbb{R}^{3}$. These are inverse shape problems where the knowledge of the
solution to a partial differential equation on the boundary is used to recover
unknown interior regions. Here we are interested in reconstructing a subregion
$D\subset\Omega$ such that dist$(D,\partial\Omega)>0$. In our models, a
Dirichlet condition is imposed on the exterior boundary $\partial\Omega$ and
the corresponding Neumann condition is measured. For the entirety of this
paper, we assume that $D$ is a collection of small volume subregions such that
$|D|=\mathcal{O}(\epsilon^{d})$, where $d$ = 2 or 3 is the dimension. To fix
the notation, we let

$D=\bigcup_{j=1}^{J}D_{j}\quad\text{with}\quad D_{j}=(x_{j}+\epsilon
B_{j})\quad\text{such that}\quad\text{dist}(x_{i},x_{j})\geq c_{0}>0$ (1)

for $i\neq j$ where the parameter $0<\epsilon\ll 1$ and $B_{j}$ is a domain
with Lipschitz boundary centered at the origin such that
$|B_{j}|=\mathcal{O}(1)$. We also assume that the individual regions $D_{j}$
are disjoint. See Figure 1 for a visual representation of the described set
up.

Figure 1: Here is an example of the type of small volume regions $D$ described
in (1) contained in a circle $\Omega$ that will be considered throughout the
paper.

In DOT, the propagation of light through a medium is modeled by the steady-
state diffusion equation. Inside the medium, we consider the case where the
absorption coefficient is zero except in the small volume subregions. In this
case, the Cauchy data corresponds to inward and outward light fluxes across
the medium’s surface. For a comprehensive description of DOT see topical
reviews [4, 14]. In our inverse scattering problem, a forcing term will be
applied to the direct scattering problem where the source term is zero except
in the small volume subregions. Here, the Dirichlet condition represents the
scattered field on the surface of the exterior boundary $\partial\Omega$. See
[26, 33, 35, 36] for more discussion on the theory and applications of this
inverse scattering problem.

In order to solve both inverse shape problems, we will develop reconstruction
algorithms that fall under the category of qualitative methods. In many
applications, qualitative methods are optimal since one of their main
advantages is that they generally require little a priori knowledge of the
unknown region $D$. Whereas iterative methods usually require a priori
information to construct a “good” initial estimate for the unknown region
and/or parameters to ensure that the iterative process will converge to the
unique solution of the inverse problem. Iterative methods can also be
computationally expensive as well as highly ill-conditioned. To avoid
requiring a priori knowledge of the small volume regions, we will analyze two
qualitative methods. From the given Dirichlet data, we will assume that we
have the corresponding normal derivative on the surface $\partial\Omega$ and
analyze its asymptotic expansion with respect to the small parameter
$\epsilon>0$.

In our DOT problem, we reconstruct $D$ by developing a MUSIC-type algorithm.
This method has been used in many imaging modalities such as acoustic [3, 6,
25, 32], electromagnetic [7, 8, 31], and elastic [15, 34] inverse scattering.
Recently, the factorization method was applied to this problem for recovering
extended regions in [19]. For our second problem associated with inverse
scattering, we derive a direct sampling method which is similar to the
orthogonality sampling method and reverse time migration, to reconstruct $D$.
This method has been widely studied for far-field measurements for several
inverse scattering problems, see for e.g. [20, 28, 30]. These methods have
also been applied to problems in DOT [10] and Electrical Impedance Tomography
[11]. The aforementioned approaches will be in combination with the so called
reciprocity gap functional defined as the surface integral

$R[v]\mathrel{\mathop{\ordinarycolon}}=\int_{\partial\Omega}v\partial_{\nu}u-u\partial_{\nu}v\,\text{d}s.$
(2)

In the two problems we are considering, the solution $u\in H^{1}(\Omega)$ with
$L^{2}$ Laplacian represents the respective fields and $v\in H^{1}(\Omega)$
with $L^{2}$ Laplacian represents a solution to the problem without the small
volume regions. This functional has been studied in [12] for another inverse
scattering problem. We utilize this functional in our asymptotic analysis in
order to reconstruct the unknown region with little a priori knowledge.

The rest of the paper is organized as follows. In Section 2, we consider an
inverse shape problem in DOT and develop the analytic framework for the MUSIC-
type algorithm, which requires multiple measurements. To do this, we apply the
reciprocity gap functional to a harmonic lifting of the Dirichlet data and the
measured Cauchy data in order to derive an imaging functional. We proceed in
Section 3 by considering the problem from inverse scattering where we derive
and analyze a direct sampling imaging functional. This method only requires a
single Cauchy pair and we show its stability with respect to error. Here the
reciprocity gap functional is applied to a plane wave and the measured Cauchy
data. In Sections 2 and 3 numerical examples are presented in $\mathbb{R}^{2}$
to validate the analysis of the constructed imaging functionals which is based
on the asymptotic expansion of the Neumann data. Lastly, in Section 4 we
provide a summary of the results of this paper and briefly discuss potential
directions for future research.

##  2 An Application to Diffuse Optical Tomography

We begin by considering the direct problem associated with DOT. This problem
stems from semiconductor theory where boundary measurements are used to
determine the existence of an interior structure. Recall, that we are
concerned with the case where these interior structures are of small volume.
We assume that the domain $\Omega\subset\mathbb{R}^{d}$ (for $d=2,3$) is a
bounded simply connected open set with Lipschitz boundary $\partial\Omega$
with unit outward normal $\nu$. We let $D\subset\Omega$ with Lipschitz
boundary $\partial D$ satisfying (1).

Now, we let $u\in H^{1}(\Omega)$ be the unique solution to

$-\Delta u+\rho\chi_{D}u=0\quad\text{in}\quad\Omega\quad\text{and}\quad
u\big{\rvert}_{\partial\Omega}=f$ (3)

for any given $f\in H^{1/2}(\partial\Omega)$ where $\chi_{(\cdot)}$ denotes
the indicator function. We assume the absorption coefficient $\rho\in
L^{\infty}(D)$. For analytical purposes of well-posedness for the direct
problem and the upcoming analysis of the inverse problem, we assume that there
are constants $\rho_{\text{min}}$ and $\rho_{\text{max}}$ such that

$0<\rho_{\text{min}}\leq\rho\leq\rho_{\text{max}}$ for a.e. $x\in D$.

One may easily verify that (3) is well-posed by considering its variational
formulation (see for e.g. [13]). Thus, one can show that for some $C>0$ that
is independent of $0<\epsilon\ll 1$ we have that

$\|{u}\|_{H^{1}(\Omega)}\leq C\|{f}\|_{H^{1/2}(\partial\Omega)}.$

By equation (3) we have that the Cauchy data is such that
$(f,\partial_{\nu}u)\in H^{1/2}(\partial\Omega)\times
H^{-1/2}(\partial\Omega)$.

In this section, we will develop the MUSIC Algorithm for solving the inverse
problem under consideration. The goal is to first derive an asymptotic
expansion for the Neumann data. Being motivated by analysis in [25, 27], we
will derive an analog of the multi-static response matrix derived from the
reciprocity gap functional (2) for this problem.

###  2.1 MUSIC Algorithm

We begin, by proving an asymptotic expansion of the Neumann data
$\partial_{\nu}u$ on $\partial\Omega$ in terms of the parameter $0<\epsilon\ll
1$. To this end, we let $u_{0}\in H^{1}(\Omega)$ be the harmonic lifting of
the Dirichlet data such that

$-\Delta u_{0}=0\quad\text{in}\quad\Omega\quad\text{and}\quad
u_{0}\big{\rvert}_{\partial\Omega}=f.$ (4)

In other words, $u_{0}$ satisfies the background problem associated with (3)
without the absorption coefficient with the same Dirichlet data $f\in
H^{1/2}(\partial\Omega)$. We continue by defining the Dirichlet Green’s
function for the negative Laplacian on the known domain $\Omega$ as
$\mathbb{G}(\cdot,z)\in H^{1}_{loc}(\Omega\setminus\left\\{z\right\\})$, which
is the unique solution to the boundary value problem

$-\Delta\mathbb{G}(\cdot,z)=\delta(\cdot,z)\enspace\text{in}\enspace\Omega\quad\text{and}\quad\mathbb{G}(\cdot,z)\big{\rvert}_{\partial\Omega}=0.$

For any fixed $z\in\Omega$, we appeal to Green’s 2nd Theorem to obtain the
representation

$\displaystyle-(u-u_{0})(z)=\int_{\Omega}(u-u_{0})(x)\Delta\mathbb{G}(x,z)\,\text{d}x$
$\displaystyle=\int_{\Omega}\mathbb{G}(x,z)\rho(x)\chi_{D}(x)u(x)\,\text{d}x$
$\displaystyle=\int_{D}\mathbb{G}(x,z)\rho(x)u(x)\,\text{d}x$

where we have used the fact that the absorption coefficient is zero outside of
the region $D$. By taking the normal derivative, we have that for all
$z\in\partial\Omega$

$\displaystyle\partial_{\nu}(u-u_{0})(z)$
$\displaystyle=-\int_{D}\rho(x)u(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x$
$\displaystyle=-\int_{D}\rho(x)u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x-\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x$
(5)

where the integrands are well defined due to the fact that
$z\in\partial\Omega$ and $\partial_{\nu(z)}$ denotes the normal derivative on
$\partial\Omega$ with respect to $z$. Given that the region $D$ satisfies (1),
we claim that (2.1) is dominated by the first integral. In other words, the
Neumann data can be approximated by the harmonic lifting $u_{0}$ restricted to
the small volume subregions, instead of the unknown photon density $u$.

The following estimate derived in Theorem 3.1 of [5] will help us in our
asymptotic analysis of (2.1). It states that for all $\varphi\in
H^{1}(\Omega)$ with $D\subset\Omega$ such that
$|D|=\mathcal{O}(\epsilon^{d})$, we have that

$\|{\varphi}\|_{L^{2}(D)}\leq
C\epsilon^{\frac{d}{2}\big{(}1-\frac{2}{p}\big{)}}\|{\varphi}\|_{H^{1}(\Omega)}$
(6)

where $p\geq 2$ in $d=2$ and $2\leq p\leq 6$ in $d=3$. This estimate is proven
using the Sobolev embedding of $H^{1}(\Omega)\hookrightarrow L^{p}(\Omega)$
(see for e.g. Chapter 5 of [1]). Using (6), we prove that $u_{0}$ approximates
$u$ when $D$ has small volume.

######  Lemma 2.1.

For all $f\in H^{1/2}(\partial\Omega)$, let $u$ and $u_{0}$ be the solutions
to (3) and (4), respectively. Then, we have that

$\|{u-u_{0}}\|_{H^{1}(\Omega)}\leq
C\epsilon^{d\big{(}1-\frac{2}{p}\big{)}}\|{f}\|_{H^{1/2}(\partial\Omega)}$

provided that $|D|=\mathcal{O}(\epsilon^{d})$ where $p\geq 2$ in $d=2$ and
$2\leq p\leq 6$ in $d=3$.

###### Proof.

Since, $u-u_{0}\in H_{0}^{1}(\Omega)$ we have that
$\|u-u_{0}\|_{H^{1}(\Omega)}\leq C\|\nabla(u-u_{0})\|_{L^{2}(\Omega)}$ by the
Poincaré inequality. By appealing to equations (3) and (4) along with Green’s
1st Theorem, we have that

$\displaystyle\int_{\Omega}\big{|}\nabla(u-u_{0})\big{|}^{2}\,\text{d}x=-\int_{D}\rho
u\overline{(u-u_{0})}\,\text{d}x$
$\displaystyle\leq\rho_{\text{max}}\|{u}\|_{L^{2}(D)}\|{u-u_{0}}\|_{L^{2}(D)}.$

Now, by using the estimate in (6), we obtain that

$\displaystyle\|{u}\|_{L^{2}(D)}\|{u-u_{0}}\|_{L^{2}(D)}$ $\displaystyle\leq
C\epsilon^{d\left(1-\frac{2}{p}\right)}\|{u}\|_{H^{1}(\Omega)}\|{u-u_{0}}\|_{H^{1}(\Omega)}.$

This proves the claim by appealing to the well-posedness of (3). ∎

From the above lemma, we have shown that $u$ can be approximated by $u_{0}$ in
norm when $|D|=\mathcal{O}(\epsilon^{d})$ is small. Under the same assumption,
we will use the previous lemma along with (6) to compare the magnitudes of the
two integrals in equation (2.1). We start by analyzing the second integral and
provide the following results.

######  Lemma 2.2.

For $z\in\partial\Omega$ and $|D|=\mathcal{O}(\epsilon^{d})$, we have that

$\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x=\mathcal{O}(\epsilon^{d+1})\quad\text{as}\quad\epsilon\rightarrow
0.$

###### Proof.

In order to prove the claim, we must estimate

$\displaystyle\left|\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x\right|$
$\displaystyle\leq
C\epsilon^{\frac{d}{2}\left(1-\frac{2}{p}\right)}\|{u-u_{0}}\|_{H^{1}(\Omega)}\|{\partial_{\nu(z)}\mathbb{G}(\cdot,z)}\|_{L^{2}(D)}$
$\displaystyle\leq
C\epsilon^{\frac{3d}{2}\left(1-\frac{2}{p}\right)}\|{f}\|_{H^{1/2}(\partial\Omega)}\|{\partial_{\nu(z)}\mathbb{G}(\cdot,z)}\|_{L^{2}(D)}$

where we have used (6), and Lemma 2.1 in order. We also have that

$\displaystyle\|{\partial_{\nu(z)}\mathbb{G}(\cdot,z)}\|_{L^{2}(D)}\leq\|{\nabla\mathbb{G}(\cdot,z)}\|_{L^{2}(D)}\leq
C\epsilon^{d/2}\|{\mathbb{G}(\cdot,z)}\|_{C^{1}(\Omega^{*})}$

by the fact that $\nu$ is a unit vector and the symmetry of the Green’s
function. The region $\Omega^{*}$ satisfies that
$D\subset\Omega^{*}\subset\Omega$ for all $0<\epsilon\ll 1$ with
dist$(\partial\Omega,\Omega^{*})>0$. Thus, we have that for all
$z\in\partial\Omega$

$\displaystyle\left|\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x\right|\leq
C\epsilon^{d\left(2-\frac{3}{p}\right)}\|{f}\|_{H^{1/2}(\partial\Omega)}.$ (7)

For $d=2$, we recall that $p\geq 2$. In order to prove the claim, we impose
the condition that

$3=2\left(2-\frac{3}{p}\right)\quad\text{making the exponent of $\epsilon$
equal to $d+1$ in \eqref{powereps}},$

which yields that $p=6$. From the above inequality we get that

$\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x\leq
C\epsilon^{3}\|{f}\|_{H^{1/2}(\partial\Omega)}.$

Similarly, for $d=3$, we recall that $2\leq p\leq 6$. Again, to prove the
claim we impose that

$4=3\left(2-\frac{3}{p}\right)\quad\text{again making the exponent of
$\epsilon$ equal to $d+1$ in \eqref{powereps}},$

which yields that $p=4.5$. Thus, we have that

$\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x\leq
C\epsilon^{4}\|{f}\|_{H^{1/2}(\partial\Omega)}.$

Therefore, for both $d=2$ and $d=3$ taking $p=6$ and $p=4.5$, respectively, we
have that

$\int_{D}\rho(x)(u-u_{0})(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x=\mathcal{O}\big{(}\epsilon^{d+1}\big{)}\quad\text{as}\quad\epsilon\to
0$

which proves the claim. ∎

Next, we show that the first integral in (2.1) is of order $\epsilon^{d}$.
From this, equation (2.1) will imply that the first integral is the leading
term, rendering the second integral as negligible. This is proven in the
following result.

######  Lemma 2.3.

For all $z\in\partial\Omega$ where $D$ is given by (1), we have that as
$\epsilon\rightarrow 0$

$\int_{D}\rho(x)u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})u_{0}(x_{j})\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})$

where Avg$(\rho_{j})$ is the average value of $\rho$ in $D_{j}$.

###### Proof.

By (1), we have that $x\in D_{j}$ if and only if $x=x_{j}+\epsilon y$ for some
$y\in B_{j}$. Now, recall that both $u_{0}$ and
$\partial_{\nu(z)}\mathbb{G}(\cdot,z)$ are smooth in the interior of of
$\Omega$ since $z\in\partial\Omega$ by standard elliptic regularity.
Therefore, we have that for all $x\in D_{j}$

$u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,z)=u_{0}(x_{j}+\epsilon
y)\partial_{\nu(z)}\mathbb{G}(x_{j}+\epsilon
y,z)=u_{0}(x_{j})\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon)$

as $\epsilon\rightarrow 0$ by appealing to Taylor’s Theorem. From this, we
obtain that

$\displaystyle\int_{D}\rho(x)u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x$
$\displaystyle=\sum_{j=1}^{J}\int_{D_{j}}\rho(x)u_{0}(x_{j}+\epsilon
y)\partial_{\nu(z)}\mathbb{G}(x_{j}+\epsilon y,z)\,\text{d}x$
$\displaystyle=\sum_{j=1}^{J}\Big{(}u_{0}(x_{j})\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon)\Big{)}\int_{D_{j}}\rho(x)\,\text{d}x.$

This implies that

$\int_{D}\rho(x)u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,z)\,\text{d}x=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})u_{0}(x_{j})\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})$

as $\epsilon\rightarrow 0$ where $\text{Avg}(\rho_{j})$ denotes the average
value of $\rho$ in $D_{j}$ as well as using the fact that
$|D_{j}|=\epsilon^{d}|B_{j}|$. ∎

Using Lemmas 2.2 and 2.3, it is clear that for a specified
$z\in\partial\Omega$, the normal derivative of the difference of $u$ and
$u_{0}$ is dominated by the first integral from equation (2.1). Therefore, we
have proven an asymptotic expansion for the boundary data $\partial_{\nu}u$ in
therms of the known harmonic lifting and Green’s function. Similar results
have been proven in [2, 18] using boundary integral operators.

######  Theorem 2.1.

For any $z\in\partial\Omega$ we have that

$\partial_{\nu(z)}u(z)=\partial_{\nu(z)}u_{0}(z)-\epsilon^{d}\sum_{j=1}^{J}u(x_{j})\text{Avg}(\rho_{j})|B_{j}|\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})\quad\text{as}\quad\epsilon\rightarrow
0$

provided that $D$ satisfies (1).

We use this asymptotic expansion to develop an algorithm that detects the
centers of the defective regions. To achieve this, we study the MUSIC
algorithm which can be considered as a discrete analogue of the factorization
method (see for e.g. [9, 24, 25]). To this end, we let $u_{0}(\cdot,g)$ and
$u_{0}(\cdot,f)$ denote the harmonic liftings with Dirichlet data $g$ and
$f\in H^{1/2}(\partial\Omega)$, respectively (see for e.g. Chapter 2 of [13]).
Thus, using Theorem 2.1 we can approximate the reciprocity gap functional (2)
with Cauchy data $(u=f,\partial_{\nu}u)$ on $\partial\Omega$ and input
$v=u_{0}(\cdot,g)$. Therefore, we have that

$\displaystyle R_{f}\big{[}u_{0}(\cdot,g)\big{]}$
$\displaystyle=\int_{\partial\Omega}u_{0}(\cdot,g)\partial_{\nu}u(\cdot,f)-u(\cdot,f)\partial_{\nu}u_{0}(\cdot,g)\,\text{d}s$
$\displaystyle=\int_{\partial\Omega}u_{0}(\cdot,g)\Big{[}\partial_{\nu}u_{0}(\cdot,f)-\epsilon^{d}\sum_{j=1}^{J}u(x_{j})\text{Avg}(\rho_{j})|B_{j}|\partial_{\nu(z)}\mathbb{G}(x_{j},z)\Big{]}\,\text{d}s$
$\displaystyle\mathrel{\phantom{=}}{}-\int_{\partial\Omega}u_{0}(\cdot,f)\partial_{\nu}u_{0}(\cdot,g)\,\text{d}s+\mathcal{O}(\epsilon^{d+1})$
$\displaystyle=-\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})u_{0}(x_{j},f)\int_{\partial\Omega}u_{0}(\cdot,g)\partial_{\nu(z)}\mathbb{G}(x_{j},z)\,\text{d}s+\mathcal{O}(\epsilon^{d+1})$

where we used the fact that
$u\big{\rvert}_{\partial\Omega}=u_{0}(\cdot,f)\big{\rvert}_{\partial\Omega}=f$,
as well as $u_{0}(\cdot,f)$ and $u_{0}(\cdot,g)$ being harmonic in $\Omega$.
Furthermore, since $z\in\partial\Omega$, we have that

$\int_{\partial\Omega}u_{0}(\cdot,g)\partial_{\nu(z)}\mathbb{G}(x_{j},z)\,\text{d}s=-u_{0}(x_{j},g).$

With this, we obtain the expansion

$R_{f}\big{[}u_{0}(\cdot,g)\big{]}=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})u_{0}(x_{j},f)u_{0}(x_{j},g)+\mathcal{O}(\epsilon^{d+1})$
(8)

as $\epsilon\rightarrow 0$ where we have made the dependance on $f$ explicit.

In order to derive the MUSIC algorithm, we assume that $\Omega$ is the unit
circle for $d=2$ where we let $g=\text{e}^{\text{i}m\theta}$ and
$f=\text{e}^{\text{i}n\theta}$ for $m,n=0,\cdots,N$ for some fixed
$N\in\mathbb{N}$. Here $\theta$ denotes the angle formed by points on
$\partial\Omega$ when converted to polar coordinates. Using only the leading
order term of (8), we define the matrix

$\textbf{F}_{n,m}=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})u_{0}(x_{j},f_{n})u_{0}(x_{j},g_{n}).$

From the definition of F, we see that it can be factorized by the matrices
$\textbf{U}\in\mathbb{C}^{(N+1)\times J}$ and
$\textbf{T}\in\mathbb{C}^{J\times J}$ that are given by

$\text{{U}}_{m,j}=u_{0}(x_{j},f_{m})\quad\text{ and
}\quad\text{{T}}=\text{diag}\big{(}\epsilon^{d}|{B_{j}}|\text{Avg}(\rho_{j})\big{)}.$

Therefore, it is easy to see that F=UTU⊤. Notice, that by our assumptions on
$\rho$, all the diagonal entries of the matrix T are non-zero. We now define
the vector $\boldsymbol{\phi}_{x}\in\mathbb{C}^{N+1}$ for any
$x\in\mathbb{R}^{d}$ by

$\boldsymbol{\phi}_{x}=\left(u_{0}(x,f_{0}),\cdots,u_{0}(x,f_{N})\right)^{\top}.$
(9)

The goal is to prove that the vector $\boldsymbol{\phi}_{x}$ is in the range
of $\text{{FF}}^{*}$ if and only if $x$ is contained in the set
$\left\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\ldots,J\right\\}$ as
similarly done in [16]. This is a discretized version of the factorization
method initially studied for this problem in [19]. See for e.g. [9, 25] for
the connection of the factorization method and MUSIC algorithm.

We now construct an imaging functional derived from the leading order term in
the asymptotic expansion of the reciprocity gap functional. To this end, we
need to show that for each sampling point $x\in\Omega$ we have that
$\boldsymbol{\phi}_{x}$ is in the range of FF∗ if and only if
$x\in\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\ldots,J\\}$. This result
has been proven in Theorem 3.2 of [16]. To avoid repetition, we will state the
following result and reference the proofs in Section 3 of [16] for details.

######  Theorem 2.2.

Assume that $N+1>J$. Then for all $x\in\Omega$ being given by the unit circle

$\boldsymbol{\phi}_{x}\in Range(\text{{FF}}^{*})\quad\text{if and only
if}\quad
x\in\left\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\ldots,J\right\\}$

where $\boldsymbol{\phi}_{x}$ is defined as in (9). Moreover, the rank of the
matrix $\text{{FF}}^{*}$ is given by $J$.

Notice, that the matrix F can be approximated by the known reciprocity gap
functional. This implies that Theorem 2.1 can be used to recover the centers
of the subregions $\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\ldots,J\\}$.
To this end, we must verify whether $\boldsymbol{\phi}_{x}\in$
Range$(\text{{FF}}^{*})$. This is equivalent to
$\text{{P}}\boldsymbol{\phi}_{x}$ = 0 where P is the orthogonal projection
onto the Null$(\text{{FF}}^{*})$.

###  2.2 Numerical Validation for the MUSIC Algorithm

We now provide some numerical examples of recovering locations the
$\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\ldots,J\\}$ using Theorem 2.2.
All of our numerical experiments are done with the software MATLAB 2020a. We
will let $\Omega$ be given by the unit circle in $\mathbb{R}^{2}$ and we need
to compute the Neumann data $\partial_{\nu}u$. It is clear that the Neumann
data can be approximated by the harmonic lifting $u_{0}$ with the same
Dirichlet data. Indeed, Lemmas 2.2 and 2.3 imply that for all
$z\in\partial\Omega$

$\partial_{\nu}u\approx\partial_{\nu}u_{0}-\int_{D}\rho(x)u_{0}(x)\partial_{\nu(z)}\mathbb{G}(x,\cdot)\,\text{d}x.$

This can be seen as an analog to the Born approximation used in scattering
theory (see for e.g. [25]). It is clear that for Dirichlet data
$\text{e}^{\text{i}n\theta}$ the harmonic lifting is given by

$u_{0}(x,\text{e}^{\text{i}n\theta})=|x|^{n}\text{e}^{\text{i}n\theta}\quad\text{
for all }\quad n\in\mathbb{N}\cup\\{0\\}.$

It is also well known that the normal derivative of $\mathbb{G}(x,z)$ is given
by

$\partial_{\nu(z)}\mathbb{G}\big{(}x,z\big{)}\big{|}_{\partial\Omega}=\frac{1}{2\pi}\left[\frac{1-|x|^{2}}{|x|^{2}+1-2|x|\text{cos}(\theta-\theta_{z})}\right]$

for $z\in\partial\Omega$. Therefore, we can compute $\partial_{\nu}u$ using
the ‘integral2’ command in MATLAB.

Given the Dirichlet data

$f_{m}=\text{e}^{\text{i}m\theta}\quad\text{ and its corresponding Neumann
data}\quad\partial_{\nu}u(\cdot,\text{e}^{\text{i}m\theta}),$

we can easily approximate the reciprocity gap functional as given by (2) using
64 equally spaced points on the unit circle for
$v=u_{0}(\cdot,\text{e}^{\text{i}n\theta})$. We have that

$R_{f_{m}}\big{[}u_{0}(\cdot,\text{e}^{\text{i}n\theta})\big{]}=\int_{\partial\Omega}\text{e}^{\text{i}n\theta}\partial_{\nu}u(\cdot,\text{e}^{\text{i}m\theta})-u(\cdot,\text{e}^{\text{i}m\theta})n\text{e}^{\text{i}n\theta}\,\text{d}s$

which is approximated via a Riemann sum using the ‘dot’ command in MATLAB for
$m,n=0,\cdots,20$. By appealing to the asymptotic result in Theorem 2.1, we
have that

$\textbf{F}_{n,m}\approx\int_{\partial\Omega}\text{e}^{\text{i}n\theta}\partial_{\nu}u(\cdot,\text{e}^{\text{i}m\theta})-u(\cdot,\text{e}^{\text{i}m\theta})n\text{e}^{\text{i}n\theta}\,\text{d}s\quad\text{
such that }\quad\textbf{F}\in\mathbb{C}^{21\times 21}.$

Once F has been approximated we can use Theorem 2.2 to recover the locations
of the components of $D$. We only need to check if the vector
$\boldsymbol{\phi}_{x}$ is in the range of $\text{{FF}}^{*}$. Therefore, we
compute the norm

$\|\text{{P}}\boldsymbol{\phi}_{x}\|^{2}_{2}=\sum\limits_{\ell=r+1}^{21}\left|\big{(}\boldsymbol{\phi}_{x},{\bf
u}_{\ell}\big{)}\right|^{2}$

where the vectors ${\bf u}_{\ell}$ are the orthonormal eigenvectors for ${\bf
F}{\bf F}^{*}$ and $r=$Rank$\big{(}{\bf F}{\bf F}^{*}\big{)}$. Recall, that
the vector

$\boldsymbol{\phi}_{x}=\bigg{(}1,|x|\text{e}^{\text{i}\theta}\,,\,\ldots\,,\,|x|^{20}\text{e}^{20\text{i}\theta}\bigg{)}^{\top}$

by equation (9) where $\theta$ is the polar angle for the sampling point
$x\in\Omega$. Here P denotes the orthogonal projection onto the
Null$(\text{{FF}}^{*})$. Therefore, the imaging functional for recovering the
centers is given by

$W_{\text{MUSIC}}(x)=\left[\sum\limits_{\ell=r+1}^{21}\left|\big{(}\boldsymbol{\phi}_{x},{\bf
u}_{\ell}\big{)}\right|^{2}\right]^{-1}\quad\text{ for any }\quad x\in\Omega$

which has the property that $W_{\text{MUSIC}}(x)\gg 1$ for $x=x_{j}$ and
$W_{\text{MUSIC}}(x)=\mathcal{O}(1)$ for $x\neq x_{j}$. We will plot the
imaging functional to provided a numerical approximation of the centers
$x_{j}$ for $j=1,\cdots,J$.

In Figures 2, 4, and 5 we use the imaging functional $W_{\text{MUSIC}}(x)$
given above to recover the locations of the two components of the region $D$.
In theses experiments, the region

$D=\big{(}x_{1}+\epsilon B(0,1)\big{)}\bigcup\big{(}x_{2}+\epsilon
B(0,1)\big{)}$

with $B(0,1)$ being the unit ball centered at the origin. The points $x_{1}$
and $x_{2}$ are points contained in the region $\Omega$. We will take the
forcing term to be given by $\rho=1$ on both components of $D$. Here, we take
$\epsilon=0.01$ as well as adding $\delta=5\%$ random noise to the computed
normal derivative of the difference of $u$ and its harmonic lifting to
simulate error in measured data.

Example 1:
In our first example presented here, we let

$x_{1}=(-0.25,0.25)\quad\text{and}\quad x_{2}=(0.25,-0.25)$

for the reconstruction in Figure 2. Here we let $\delta=5\%$ and $\rho=1$ in
both subregions. Presented is a contour and surface plot of the imaging
functional $W_{\text{MUSIC}}(x)$. As we can see, the imaging functional is
elevated in the general region around the centers.

Figure 2: Initial reconstruction of the locations $x_{1}=(-0.25,0.25)$ and
$x_{2}=(0.25,-0.25)$ via the MUSIC algorithm with the faulty rank of
$\textbf{FF}^{*}$. Contour plot on the left and Surface plot on the right of
the imaging function $W_{\text{MUSIC}}(x)$.

Recall, that the imaging functional $W_{\text{MUSIC}}(x)$ depends on the rank
of $\textbf{FF}^{*}$, which was calculated using the rank function in MATLAB.
However, $\|{\textbf{FF}^{*}}\|\ll 1$ and the default tolerance of the rank
function produces an overestimation of the true rank since the singular values
of $\textbf{FF}^{*}$ are very small. For the rest of our numerical
experiments, we improve the rank calculation by computing the singular values
of $\textbf{FF}^{*}$ and ad hoc checking when a singular value decreases by at
least 3 orders of magnitude from the previous one.

Figure 3: Semi-log plot of the svd($\textbf{FF}^{*}$). Labeled are its 2nd and
3rd singular values.

Figure 3, suggests that the actual rank of $\textbf{FF}^{*}$ is 2, as expected
by Theorem 2.2 since there are 2 components of $D$. Throughout the remaining
examples of this section, we will continue to heuristically calculate the rank
of $\textbf{FF}^{*}$ with this method. With this new method of calculating the
rank, Figure 4 demonstrates a clearer reconstruction of the 2 subregions
centered at the locations $x_{1}=(-0.25,0.25)$ and $x_{2}=(0.25,-0.25)$. As we
can see from the data tips, the improved imaging functional has spikes at the
points

$\widetilde{x}_{1}=(-0.2462,0.2462)\quad\text{and}\quad\widetilde{x}_{2}=(0.2462,-0.2462).$

Here we see that the locations $\widetilde{x}_{1}$ and $\widetilde{x}_{2}$
given by the MUSIC algorithm provide a good approximation for the actual
locations of the components of the region D.

Figure 4: Reconstruction of the locations $x_{1}=(-0.25,0.25)$ and
$x_{2}=(0.25,-0.25)$ via the MUSIC algorithm. Contour plot on the left and
Surface plot on the right of the imaging functional $W_{\text{MUSIC}}(x)$.

Example 2:
In our second example presented here, we let

$x_{1}=(-0.25,0.25)\quad\text{and}\quad x_{2}=(-0.25,-0.25)$

for the reconstruction in Figure 5. Presented is a contour and surface plot of
the imaging functional $W_{\text{MUSIC}}(x)$. We again let $\delta=5\%$ and
$\rho=1$ in both subregions. As we can see from the data tips, the imaging
functional has spikes at the points

$\widetilde{x}_{1}=(-0.2462,0.2462)\quad\text{and}\quad\widetilde{x}_{2}=(-0.2462,-0.2462).$

Again, in this example we see that the locations of $\widetilde{x}_{1}$ and
$\widetilde{x}_{2}$ provide an approximation for the locations of the
components of the region $D$.

Figure 5: Reconstruction of the locations $x_{1}=(-0.25,0.25)$ and
$x_{2}=(-0.25,-0.25)$ via the MUSIC algorithm. Contour plot on the left and
Surface plot on the right of the imaging functional $W_{\text{MUSIC}}(x)$.

In our final two examples of this section, we similarly let the region

$D=\bigcup_{j=1}^{J}\big{(}x_{j}+\epsilon B(0,1)\big{)}$

where $J=3$ and $J=4$, respectively, with $B(0,1)$ being the unit ball
centered at the origin. The points $x_{j}$ for $j=1,\cdots,J$ are contained in
the region $\Omega$ and we once again let $\epsilon=0.01$. However, we vary
the value of the forcing term $\rho$ on each of the components of $D$.
Furthermore, we add $\delta=10\%$ random noise to the approximated normal
derivative of the difference of $u$ and its harmonic lifting to simulate error
in measured data.

Example 3:
In our third example presented here, we let

$x_{1}=(-0.75,0),\quad x_{2}=(0.25,0.5),\quad\text{and}\quad
x_{3}=(-0.3,-0.4)$

for the reconstruction in Figure 6. Presented is a contour and surface plot of
the imaging functional $W_{\text{MUSIC}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(-0.7487,-0.0151),\quad\widetilde{x}_{2}=(0.2563,0.4975),\quad\text{and}\quad\widetilde{x}_{3}=(-0.2965,-0.407).$

In this example we see that the locations of $\widetilde{x}_{1}$,
$\widetilde{x}_{2}$, and $\widetilde{x}_{3}$ provide an approximation for the
locations of the components of the region $D$. For this example we let
$\delta=10\%$ where $\rho=1/4$ in the region centered at $x_{1}$, $\rho=1$ in
the region centered at $x_{2}$, and $\rho=2$ in the region centered at
$x_{3}$. Notice that this example suggests that the MUSIC algorithm gives
sharper reconstructions when the regions are well separated.

Figure 6: Reconstruction of the locations $x_{1}=(-0.75,0)$,
$x_{2}=(0.25,0.5)$, and $x_{3}=(-0.3,-0.4)$ via the MUSIC algorithm. Contour
plot on the left and Surface plot on the right of the imaging functional
$W_{\text{MUSIC}}(x)$.

Example 4:
In our final example presented here, we let

$x_{1}=(0,0.75),\quad x_{2}=(-0.25,0.25),\quad
x_{3}=(0.25,-0.25)\quad\text{and}\quad x_{4}=(0.2,-0.6)$

for the reconstruction in Figure 7. Presented is a contour and surface plot of
the imaging functional $W_{\text{MUSIC}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(-0.005,0.7487),\quad\widetilde{x}_{2}=(-0.2362,0.2362),\quad\widetilde{x}_{3}=(0.2362,-0.2562),$

and

$\widetilde{x}_{4}=(0.206,-0.5879).$

In this example, we see that the reconstructed locations provide an
approximation for the locations of the components of the region $D$. For this
example, we let $\delta=10\%$ where $\rho=3/4$ in the region centered at
$x_{1}$, $\rho=1$ in the region centered at $x_{2}$, $\rho=3/2$ in the region
centered at $x_{3}$, and $\rho=1/2$ in the region centered at $x_{4}$.

Figure 7: Reconstruction of the locations $x_{1}=(0,0.75)$,
$x_{2}=(-0.25,0.25),$ $x_{3}=(0.25,-0.25)$ and $x_{4}=(0.2,-0.6)$ via the
MUSIC algorithm. Contour plot on the left and Surface plot on the right of the
imaging functional $W_{\text{MUSIC}}(x)$.

##  3 An Application to Inverse Scattering

We now consider the direct problem in inverse scattering where the governing
physical equation is the Helmholtz equation. Inverse scattering has many
scientific applications in medical imagining, non-destructive testing, as well
as geophysics. We are particularly concerned with detecting small volume
hidden objects within a complex media in the case where one can only make
measurements on an exterior surface. Just as in the previous section, we
assume that the domain $\Omega\subset\mathbb{R}^{d}$ (for $d=2,3$) is a
bounded, simply connected open set with Lipschitz boundary $\partial\Omega$
with unit outward normal $\nu$. We let $D\subset\Omega$ with Lipschitz
boundary $\partial D$ satisfying (1). Now, let the scattered field $u^{s}\in
H^{1}(\Omega)$ satisfy

$\Delta
u^{s}+k^{2}u^{s}=\rho\chi_{D}\quad\text{in}\quad\Omega\quad\text{and}\quad
u^{s}\big{\rvert}_{\partial\Omega}=f$ (10)

for any given $f\in H^{1/2}(\partial\Omega)$ where once again $\chi_{(\cdot)}$
denotes the indicator function. We let $k$ denote the wavenumber where we
assume $k^{2}$ is not a Dirichlet eigenvalue of $-\Delta$ in $\Omega$. With
this assumption on the wave number, we have that (10) is well-posed provided
that the source $\rho\in L^{\infty}(D)$. By equation (10) we have that the
Cauchy data is such that $(f,\partial_{\nu}u^{s})\in
H^{1/2}(\partial\Omega)\times H^{-1/2}(\partial\Omega)$.

In this section, we will develop a direct sampling method for solving the
inverse shape problem. This method has been employed for other imaging
modalities such as DOT [10] and Electrical Impedance Tomography [11]. See also
[21, 22, 29] for applications with near field measurements. MUSIC-type
algorithms has also been extensively used for similar shape reconstruction
problems in [8, 16, 31, 32]. However, our method only requires pair of Cauchy
data to recover the support of the source and also avoids matrix operations.
Lastly, our method is also highly tolerant to noise.

###  3.1 Direct Sampling Method

We denote $u_{0}^{s}\in H^{1}(\Omega)$ as the lifting which solves the
Helmholtz equation such that

$\Delta
u_{0}^{s}+k^{2}u_{0}^{s}=0\enspace\text{in}\enspace\Omega\quad\text{and}\quad
u_{0}^{s}\big{\rvert}_{\partial\Omega}=f.$ (11)

Therefore, $u_{0}^{s}$ satisfies the background problem (10) (i.e. without the
forcing term) with Dirichlet data $f\in H^{1/2}(\partial\Omega)$ and
wavenumber $k$. By our assumption on the wave number we have that (11) is also
well-posed. We proceed by defining the Dirichlet Green’s function for the
Helmholtz equation for the known domain $\Omega$ as
$\mathbb{G}_{k}(\cdot,z)\in H^{1}_{loc}(\Omega\setminus\left\\{z\right\\})$,
which is the unique solution to the boundary value problem

$\Delta\mathbb{G}_{k}(\cdot,z)+k^{2}\mathbb{G}_{k}(\cdot,z)=-\delta(\cdot,z)\enspace\text{in}\enspace\Omega\quad\text{and}\quad\mathbb{G}_{k}(\cdot,z)\big{\rvert}_{\partial\Omega}=0.$

Here, we again assume that the wavenumber $k$ is as in (10) and (11). For any
fixed $z\in\Omega$, we appeal to Green’s 2nd Theorem to obtain the
representation

$\displaystyle-(u^{s}-u_{0}^{s})(z)$
$\displaystyle=\int_{\Omega}(u^{s}-u_{0}^{s})(x)\big{[}\Delta\mathbb{G}_{k}(x,z)+k^{2}\mathbb{G}_{k}(x,z)\big{]}\,\text{d}x$
$\displaystyle=\int_{D}\mathbb{G}_{k}(x,z)\rho(x)\,\text{d}x$

where we used the indicator function from our source term. By taking the
normal derivative, we have that for all $z\in\partial\Omega$

$\partial_{\nu}(u^{s}-u_{0}^{s})(z)=-\int_{D}\rho(x)\partial_{\nu(z)}\mathbb{G}_{k}(x,z)\,\text{d}x$
(12)

where the integrand is well defined since $z\in\partial\Omega$. Again, we let
$\partial_{\nu(z)}$ denote the normal derivative on $\partial\Omega$ with
respect to $z$. We now begin our asymptotic analysis of the normal derivative
where $D$ is the finite union of small volume regions as given by (1). The
following lemma is key in deriving the asymptotic expansion.

######  Lemma 3.1.

For all $z\in\partial\Omega$ where $D$ is given by (1), we have that as
$\epsilon\rightarrow 0$

$\int_{D}\rho(x)\mathbb{G}_{k}(x,z)\,dx=-\epsilon^{d}\sum_{j=1}^{J}|B_{j}|Avg(\rho_{j})\partial_{\nu(z)}\mathbb{G}_{k}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})$

where $Avg(\rho_{j})$ is the average value of $\rho$ in $D_{j}$.

###### Proof.

By (1), we have that $x\in D_{j}$ if and only if $x=x_{j}+\epsilon y$ for some
$y\in B_{j}$. Since $z\in\partial\Omega$, then
$\partial_{\nu(z)}\mathbb{G}(\cdot,z)$ is smooth in the interior of $\Omega$
by elliptic regularity. Therefore, we have that for all $x\in D_{j}$

$\partial_{\nu(z)}\mathbb{G}_{k}(x,z)=\partial_{\nu(z)}\mathbb{G}_{k}(x_{j}+\epsilon
y,z)=\partial_{\nu(z)}\mathbb{G}_{k}(x_{j},z)+\mathcal{O}(\epsilon)$

as $\epsilon\rightarrow 0$ by appealing to Taylor’s Theorem. From this, we
obtain that

$\displaystyle\int_{D}\rho(x)\partial_{\nu(z)}\mathbb{G}_{k}(x,z)\,\text{d}x$
$\displaystyle=\sum_{j=1}^{J}\int_{D_{j}}\rho(x)\partial_{\nu(z)}\mathbb{G}_{k}(x_{j}+\epsilon
y,z)\,\text{d}x$
$\displaystyle=\sum_{j=1}^{J}\Big{(}\partial_{\nu(z)}\mathbb{G}_{k}(x_{j},z)+\mathcal{O}(\epsilon)\Big{)}\int_{D_{j}}\rho(x)\,\text{d}x$

Therefore, we have that

$\int_{D}\rho(x)\partial_{\nu(z)}\mathbb{G}_{k}(x,z)\,\text{d}x=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})\partial_{\nu(z)}\mathbb{G}_{k}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})$

as $\epsilon\rightarrow 0$ where we used the fact that
$|D_{j}|=\epsilon^{d}|B_{j}|$ and $\text{Avg}(\rho_{j})$ denotes the average
value of $\rho$ in $D_{j}$. ∎

From the above lemma, it is clear that for a specified $z\in\partial\Omega$,
the normal derivative of the difference of $u^{s}$ and the lifting $u_{0}^{s}$
is approximated by the centers of the inclusions.

######  Theorem 3.1.

For any $z\in\partial\Omega$ we have that

$\partial_{\nu(z)}u^{s}(z)=\partial_{\nu(z)}u_{0}^{s}(z)-\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})\partial_{\nu(z)}\mathbb{G}_{k}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})\quad\text{as}\quad\epsilon\rightarrow
0$

provided that $D$ satisfies (1).

With this approximation to the Neumann data, we develop an algorithm that
detects the centers of small volume regions within our domain. We now study a
direct sampling method. This is done by using Theorem 3.1 and evaluating the
reciprocity gap functional $R[v]$ given by (2), where the Cauchy data
$(u^{s}=f,\partial_{\nu}u^{s})$ on $\partial\Omega$ is fixed. Recall, that we
assume that $v\in H^{1}(\Omega)$ solves the Helmholtz equation in $\Omega$
which gives that

$\displaystyle R[v]$
$\displaystyle=\int_{\partial\Omega}v\partial_{\nu}u^{s}-u^{s}\partial_{\nu}v\,\text{d}s$
$\displaystyle=\int_{\partial\Omega}v\Big{[}\partial_{\nu}u_{0}^{s}-\epsilon^{d}\sum_{j=1}^{J}\text{Avg}(\rho_{j})|B_{j}|\partial_{\nu(z)}\mathbb{G}(x_{j},z)+\mathcal{O}(\epsilon^{d+1})\Big{]}-u_{0}^{s}\partial_{\nu}v\,\text{d}s$
$\displaystyle=-\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})\int_{\partial\Omega}v\partial_{\nu(z)}\mathbb{G}(x_{j},z)\,\text{d}s+\mathcal{O}(\epsilon^{d+1})$
$\displaystyle=\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})v(x_{j})+\mathcal{O}(\epsilon^{d+1})$

where we used (11) as well as the fact that
$u^{s}\big{\rvert}_{\partial\Omega}=u_{0}^{s}\big{\rvert}_{\partial\Omega}=f$.

Notice, we can take $v=\text{e}^{\text{i}kz\cdot\hat{y}}$, which is clearly a
solution to the Helmholtz equation for all $z\in\mathbb{R}^{d}$, when
$\hat{y}\in\mathbb{S}^{d-1}$(i.e. unit circle/sphere). We proceed by defining
the imaging functional
$W(z)\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq
0}$ as

$\displaystyle
W(z)=\left|\big{(}R[\text{e}^{\text{i}kz\cdot\hat{y}}],\text{e}^{\text{i}kz\cdot\hat{y}}\big{)}_{L^{2}(\mathbb{S}^{d-1})}\right|.$
(13)

This functional can be used to recover the region $D$ by plotting it’s values
in $\Omega$. To prove this fact, we will study the resolution analysis for
this imaging functional. This will involve using the asymptotic expansion
derived in Theorem 3.1 to write the functional in terms of Bessel functions.
To this end, notice that

$\displaystyle W(z)$
$\displaystyle=\left|\left(\epsilon^{d}\sum_{j=1}^{J}|B_{j}|\text{Avg}(\rho_{j})\text{e}^{\text{i}kx_{j}\cdot\hat{y}}+\mathcal{O}(\epsilon^{d+1}),\text{e}^{\text{i}kz\cdot\hat{y}}\right)_{L^{2}(\mathbb{S}^{d-1})}\right|$
$\displaystyle=\left|\epsilon^{d}\sum_{j=1}^{J}\text{Avg}(\rho_{j})|B_{j}|\int_{\mathbb{S}^{d-1}}\text{e}^{\text{i}k(x_{j}-z)\cdot\hat{y}}\,\text{d}s(\hat{y})\right|+\mathcal{O}(\epsilon^{d+1})\
$

where we used straightforward calculations and the asymptotic expansion of the
reciprocity gap functional. Now, we will recall the Funk–Hecke integral
identity

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

see for e.g. [17, 28]. Therefore, it is now clear that for all $z\in\Omega$,

$W(z)=\begin{cases}\pm\epsilon^{2}\,2\pi\displaystyle\sum_{j=1}^{J}\text{Avg}(\rho_{j})|B_{j}|J_{0}(k|x_{j}-z|)+\mathcal{O}(\epsilon^{3}),&d=2\\\
\\\
\pm\epsilon^{3}\,4\pi\displaystyle\sum_{j=1}^{J}\text{Avg}(\rho_{j})|B_{j}|j_{0}(k|x_{j}-z|)+\mathcal{O}(\epsilon^{4}),&d=3\end{cases}$
(14)

where $J_{0}$ represents the zeroth order Bessel function of the first kind
and $j_{0}$ represents the zeroth order spherical Bessel function of the first
kind. This allows us to provide our main result of this section.

######  Theorem 3.2.

Up to leading order, if $\text{Avg}(\rho_{j})\neq 0$ we have that for all
$z\in\mathbb{R}^{d}\setminus D$,

$W(z)=\mathcal{O}\Big{(}\text{dist}(z,\mathcal{X})^{\frac{1-d}{2}}\Big{)}\quad\text{as}\quad\text{dist}(z,\mathcal{X})\rightarrow\infty$

provided that the region $D$ satisfies (1), where the set
$\mathcal{X}=\left\\{x_{j}\mathrel{\mathop{\ordinarycolon}}1,\ldots,J\right\\}$.

###### Proof.

In order to prove the result, we use the fact that

$J_{0}(|z-x|)=\mathcal{O}\left(|z-x|^{-1/2}\right)\quad\text{ and }\quad
j_{0}(|z-x|)=\mathcal{O}\left(|z-x|^{-1}\right)$

as $|z-x|\to\infty$ for the case when $d=2$ or $3$ along with the expansion in
(14). ∎

Thus, Theorem 3.2 can be used to recover the centers of the subregions since
the imaging functional $W(z)$ attains a local maximum at each of the centers.
We now introduce a lemma regarding the stability for the reciprocity gap
functional. This will help us obtain a stability estimate for $W(z)$.

######  Lemma 3.2.

For added random noise $0<\delta<1$, we have that for any solution $v\in
H^{1}(\Omega)$ to the Helmholtz equation,

$\big{|}R[v]-R^{\delta}[v]\big{|}\leq C\delta\|{v}\|_{H^{1}(\Omega)}$

where $R[v]$ is given by (2) and the perturbed reciprocity gap functional is
given by

$R^{\delta}[v]=\int_{\partial\Omega}v\partial_{\nu}u_{\delta}^{s}-u_{\delta}^{s}\partial_{\nu}\,ds$
(15)

provided that there are positive constants $C_{1}$ and $C_{2}$ such that

$\|\partial_{\nu}(u_{\delta}^{s}-u^{s})\|_{H^{-1/2}(\partial\Omega)}\leq
C_{1}\delta\quad\text{ and
}\quad\|u_{\delta}^{s}-u^{s}\|_{H^{1/2}(\partial\Omega)}\leq C_{2}\delta.$

###### Proof.

By simply subtracting the expressions, we have that

$\displaystyle R[v]-R^{\delta}[v]$
$\displaystyle=\int_{\partial\Omega}v(\partial_{\nu}u^{s}-\partial_{\nu}u_{\delta}^{s})-(u^{s}-u_{\delta}^{s})\partial_{\nu}v\,\text{d}s$

Thus, we can estimate the above quantity such that

$\displaystyle|R[v]-R^{\delta}[v]|$
$\displaystyle\leq\Big{(}\|{u_{\delta}^{s}-u^{s}}\|_{H^{1/2}(\partial\Omega)}\|{\partial_{\nu}v}\|_{H^{-1/2}(\partial\Omega)}$
$\displaystyle\hskip
28.90755pt+\|{v}\|_{H^{1/2}(\partial\Omega)}\|{\partial_{\nu}(u_{\delta}^{s}-u^{s})}\|_{H^{-1/2}(\partial\Omega)}\Big{)}\leq
C\delta\|v\|_{H^{1}(\Omega)}$

by the dual-pairing of $H^{1/2}(\partial\Omega)$ and
$H^{-1/2}(\partial\Omega)$. We have also used Trace Theorems and the fact that
$v$ solves Helmholtz equation in $\Omega$. This proves the claim. ∎

We are now able to present the following theorem on the error estimate for the
imaging functional $W(z)$.

######  Theorem 3.3.

For added random noise $0<\delta<1$, we have that for any
$z\in\mathbb{R}^{d}$,

$|W(z)-W^{\delta}(z)|=\mathcal{O}(\delta)\quad\text{ as }\quad\delta\to 0$
(16)

such that the perturbed imaging functional is defined as

$W^{\delta}(z)=\Big{|}\big{(}R^{\delta}[e^{ikz\cdot\hat{y}}],e^{ikz\cdot\hat{y}}\big{)}_{L^{2}(\mathbb{S}^{d-1})}\Big{|}$

where $R^{\delta}[\cdot]$ is defined as in (15).

###### Proof.

By the Triangle and Cauchy-Schwarz inequalities, we have that

$|W(z)-W^{\delta}(z)|\leq\|{R[\text{e}^{\text{i}kz\cdot\hat{y}}]-R^{\delta}[\text{e}^{\text{i}kz\cdot\hat{y}}]}\|_{L^{2}(\mathbb{S}^{d-1})}\|{\text{e}^{\text{i}kz\cdot\hat{y}}}\|_{L^{2}(\mathbb{S}^{d-1})}.$

Note, that by the previous result in Lemma 15, we have that

$\|{R[\text{e}^{\text{i}kz\cdot\hat{y}}]-R^{\delta}[\text{e}^{\text{i}kz\cdot\hat{y}}]}\|_{L^{2}(\partial\Omega)}\leq
C\delta\|\text{e}^{\text{i}kz\cdot\hat{y}}\|_{H^{1}(\Omega)}.$

Furthermore, we have that both
$\|\text{e}^{\text{i}kz\cdot\hat{y}}\|_{H^{1}(\Omega)}$ and
$\|{\text{e}^{\text{i}kz\cdot\hat{y}}}\|_{L^{2}(\mathbb{S}^{d-1})}$ are
bounded and independent of the parameter $\delta$. Thus, we have that

$\displaystyle|W(z)-W^{\delta}(z)|\leq C\delta\quad\text{ as }\quad\delta\to
0$

which proves the claim. ∎

This result demonstrates that the imaging functional $W(z)$ is stable with
respect to error in the measured Cauchy data. This implies that plotting the
imaging functional is an analytically rigorous as well as computationally
simple and stable.

###  3.2 Numerical Validation for the Direct Sampling Algorithm

In this section, we provide some numerical examples for recovering the
locations of the unknown source given by
$\left\\{x_{j}\mathrel{\mathop{\ordinarycolon}}j=1,\cdots,J\right\\}$ using
Theorem 3.2. Just as in the previous section, all of our numerical experiments
are once again done with the software MATLAB 2020a. For simplicity, we let
$\Omega$ be given by the unit ball in $\mathbb{R}^{2}$. In order to do so, we
first need a way to calculate the corresponding scattered field $u^{s}$
solving (10). To this end, we can take the radiation scattered field for all
$x\in\mathbb{R}^{2}$ given by

$x\longmapsto-\int_{\mathbb{R}^{2}}\rho(y)\chi_{D}(y)\Phi_{k}(x,y)\,\text{d}y.$

This scattered field solves the associated source problem in all of
$\mathbb{R}^{2}$ where $\Phi_{k}$ denotes the radiating fundamental solution
to the Helmholtz equation. Since $\chi_{(\cdot)}$ denotes the indicator
function, we have that for all $x\in\Omega$

$u^{s}(x)=-\int_{D}\rho(y)\frac{\text{i}}{4}H_{0}^{(1)}(k|x-y|)\,\text{d}y$
(17)

solves (10) with the corresponding Dirichlet data. It is a well known fact
that the fundamental solution is given by

$\Phi_{k}(x,y)=\frac{\text{i}}{4}H_{0}^{(1)}(k|x-y|)$

where $H_{0}^{(1)}$ represents the first kind Hankel function of order zero.

Next, we compute the normal derivative of the scattered field. It is
straightforward to conclude that the normal derivative on $\partial\Omega$ of
the solution $u^{s}$ is given by

$\partial_{\nu}u^{s}(x)=\int_{D}\rho(y)\frac{\text{i}k}{4}H_{1}^{(1)}(k|x-y|)\bigg{[}\frac{1-x\cdot
y}{|x-y|}\bigg{]}\,\text{d}y$ (18)

where $H_{1}^{(1)}$ represents the first kind Hankel function of order one. We
calculate the scattered field and its normal derivative as given by (17) and
(18), respectively, using the ‘integral2’ command in MATLAB. Here, we evaluate
the reciprocity gap functional $R[\text{e}^{\text{i}kx\cdot\hat{y}}]$ for 64
equally spaced points $\hat{y}$ on the unit circle. By appealing to our
asymptotic result in (14), the imaging functional is given by

$W_{\text{DIRECT}}(x)=\left|\big{(}R[\text{e}^{\text{i}kx\cdot\hat{y}}],\text{e}^{\text{i}kx\cdot\hat{y}}\big{)}_{L^{2}(\mathbb{S}^{1})}\right|^{p}\quad\text{for
any}\quad x\in\Omega$

which is approximated via a Riemann sum using the ‘dot’ command in MATLAB. In
our calculations $p>0$ is a fixed chosen parameter to sharpen the resolution
of the imaging functional. We also normalize the values of the imaging
functional and pick $p=4$ in our calculations such that
$W_{\text{DIRECT}}(x)=\mathcal{O}(1)$ for $x=x_{j}$ and
$W_{\text{DIRECT}}(x)\approx 0$ for $x\neq x_{j}$.

In all our examples, we use the imaging functional $W_{\text{DIRECT}}(x)$
given above to recover the location of the components of the region $D$. In
these experiments, the region

$D=\bigcup_{j=1}^{J}\big{(}x_{j}+\epsilon B(0,1)\big{)}$

with $B(0,1)$ being the unit circle centered at the origin. Here, we take
$\epsilon=0.01$ as well as adding random noise level $\delta$ to the simulate
data $u^{s}$ and $\partial_{\nu}u^{s}$ on $\partial\Omega$. We let the wave
number $k=25$ and the points $x_{j}$ are points contained in the region
$\Omega$. In Examples 1 and 2, $D$ is composed of two regions. In Example 3,
$D$ is composed of three regions, and in Example 4, $D$ is composed of four
regions.

Example 1:
In our first example presented here, we let

$x_{1}=(0,0.75)\quad\text{and}\quad x_{2}=(0.5,0)$

for the reconstruction in Figure 8. Presented is a contour and surface plot of
the imaging functional $W_{\text{DIRECT}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(0.0101,0.7374)\quad\text{and}\quad\widetilde{x}_{2}=(0.4949,0.0101).$

We can see that the locations of $\widetilde{x}_{1}$ and $\widetilde{x}_{2}$
given by the Direct Sampling Algorithm provide an approximation for the
locations of the components of the region $D$. Here we let noise level
$\delta=1\%$ and $\rho=1$ in both subregions.

Figure 8: Reconstruction of the locations $x_{1}=(0,0.75)$ and $x_{2}=(0.5,0)$
via the imaging function $W_{\text{DIRECT}}(x)$. Contour plot on the left and
Surface plot on the right.

For the rest of the examples of this section, we vary the value of the forcing
term $\rho$ on each of the components of $D$. Furthermore, we also increment
the random noise level $\delta$ to demonstrate the stability of the method.

Example 2:
For our second example presented here, we let

$x_{1}=(0.15,0.5)\quad\text{and}\quad x_{2}=(0.35,0.2)$

for the reconstruction in Figure 9. Presented is a contour and surface plot of
the imaging functional $W_{\text{DIRECT}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(0.1717,0.4747)\quad\text{and}\quad\widetilde{x}_{2}=(0.3333,0.2121).$

Again, in this example we see that the locations of $\widetilde{x}_{1}$ and
$\widetilde{x}_{2}$ given by the direct sampling method provide an
approximation for the locations of the components of the region $D$. Here we
let noise level $\delta=10\%$ where $\rho=0.9$ in the region centered at
$x_{1}$ and $\rho=1$ in the region centered at $x_{2}$. In this example,
notice that we have reduced the distance between $x_{1}$ and $x_{2}$ and
incremented noise level $\delta$ from Example 1. Thus, the sharp
reconstruction of $D$ as shown in Figure 9 illustrates the stability and
robustness of this method.

Figure 9: Reconstruction of the locations $x_{1}=(0.15,0.5)$ and
$x_{2}=(0.35,0.2)$ via the imaging functional $W_{\text{DIRECT}}(x)$. Contour
plot on the left and Surface plot on the right.

Example 3:
In our third example presented here, we let

$x_{1}=(-0.5,-0.5),\quad x_{2}=(0,0)\quad\text{and}\quad x_{3}=(0.5,0.25)$

for the reconstruction in Figure 10. Presented is a contour and surface plot
of the imaging functional $W_{\text{DIRECT}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(-0.5152,-0.5152),\quad\widetilde{x}_{2}=(0.0101,0.0101)\quad\text{and}\quad\widetilde{x}_{3}=(0.4949,0.2525).$

In this example we see that the locations of $\widetilde{x}_{1}$,
$\widetilde{x}_{2}$ and $\widetilde{x}_{3}$ provide an approximation for the
locations of the components of the region $D$. For this example we let
$\delta=20\%$ where $\rho=0.8$ in the region centered at $x_{1}$, $\rho=1.1$
in the region centered at $x_{2}$, and $\rho=0.9$ in the region centered at
$x_{3}$.

Figure 10: Reconstruction of the locations $x_{1}=(-0.5,-0.5)$, $x_{2}=(0,0)$
and $x_{3}=(0.5,0.25)$ via the imaging functional $W_{\text{DIRECT}}(x)$.
Contour plot on the left and Surface plot on the right.

Example 4:
In our final example presented here, we let

$x_{1}=(0,0.5),\quad x_{2}=(0.25,0.25),\quad
x_{3}=(-0.25,-0.25)\quad\text{and}\quad x_{4}=(0,-0.75)$

for the reconstruction in Figure 7. Presented is a contour and surface plot of
the imaging functional $W_{\text{DIRECT}}(x)$. As we can see from the data
tips, the imaging functional has spikes at the points

$\widetilde{x}_{1}=(0.0101,0.4949),\quad\widetilde{x}_{2}=(0.2525,0.2727),\quad\widetilde{x}_{3}=(-0.2525,-0.2727)$

and

$\widetilde{x}_{4}=(-0.0101,-0.7576).$

In this example we see that the reconstructed locations provide an
approximation for the locations of the region $D_{j}$. For this example we let
$\delta=25\%$ where $\rho=0.95$ in the region centered at $x_{1}$, $\rho=1$ in
the region centered at $x_{2}$, $\rho=0.9$ in the region centered at $x_{3}$,
and $\rho=1.1$ in the region centered at $x_{4}$.

Figure 11: Reconstruction of the locations $x_{1}=(0,0.5)$,
$x_{2}=(0.25,0.25)$, $x_{3}=(-0.25,-0.25)$ and $x_{4}=(0,-0.75)$ via the
imaging functional $W_{\text{DIRECT}}(x)$. Contour plot on the left and
Surface plot on the right.

##  4 Conclusions

In this paper, we studied the use of qualitative methods for small volume
inverse shape problems in DOT and inverse scattering. In both cases, we
analyzed the asymptotic expansion of the reciprocity gap functional (2) in
order to construct an imaging functional to recover the region of interest
$D$. For the DOT problem, we have studied the MUSIC algorithm. Whereas in the
inverse scattering problem, we derived a direct sampling method. We note that
the analysis provided here can be used to study the inverse scattering
problems in $\mathbb{R}^{d}$ for $d=2,3$, where one can use (17) and the
asymptotic analysis presented here. Both algorithms allow for fast and
accurate reconstruction with little a priori knowledge of $D$. A future
direction for this project, in the area of inverse scattering can be to study
the problem in Section 3 for the case of electromagnetic and elastic
scattering. Another interesting project would be to develop a direct sampling
method as in [10] for the DOT problem presented in Section 2.

Acknowledgments: The research of G. Granados and I. Harris is partially
supported by the NSF DMS Grant 2107891.

## References

  * [1] R. Adams,  _“Sobolev Spaces”_ ,  1st edition, Academic Press London, 1975.
  * [2] H. Ammari, H. Kang, E. Kim, K. Louati, and M. Vogelius,  A MUSIC-type algorithm for detecting internal corrosion from electrostatic boundary measurements.  Numer. Math., 108, (2008), 501–528.
  * [3] H. Ammari, E. Iakovleva, and D. Lesselier,  A MUSIC Algorithm for Locating Small Inclusions Buried in a Half-Space from the Scattering Amplitude at a Fixed Frequency.  Multiscale Model. Simul., 3:3, (2005), 597–628
  * [4] S. R. Arridge,  Optical tomography in medical imaging.  Inverse Problems, 15:R41-R93, (1999)
  * [5] F. Cakoni, I. Harris, and S. Moskow,  The Imaging of Small Perturbations in an Anisotropic Media.  Comp. Math. App., 74:11, (2017), 2769–2783
  * [6] F. Cakoni and J. Rezac,  Direct imaging of small scatterers using reduced time dependent data.  J. Comp. Physics, 338, (2017), 371–387
  * [7] D. Challa, G. Hu and M. Sini,  Multiple scattering of electromagnetic waves by finitely many point-like obstacles.  Math. Models Methods in Appl. Sci., 24:5, (2014), 863–899.
  * [8] X. Chen and Y. Zhong,  MUSIC electromagnetic imaging with enhanced resolution for small inclusions.  Inverse Problems, 25, (2009), 015008
  * [9] M. Cheney,  The linear sampling method and the MUSIC algorithm.  Inverse Problems, 17, (2001), 591595
  * [10] Y.T. Chow, K. Ito, K. Liu and J. Zou,  Direct Sampling Method for Diffusive Optical Tomography.  SIAM J. Sci. Comput., 37:4, (2015), A1658–A1684.
  * [11] Y.T. Chow, K. Ito, K. Liu and J. Zou,  Direct Sampling Method for Electrical Impedance Tomography.  Inverse Problems, 30, (2014), 095003.
  * [12] D. Colton and H. Haddar,  An application of the reciprocity gap functional to inverse scattering theory.  Inverse Problems, 21, (2005)
  * [13] L. Evans,  _“Partial Differential Equation”_ ,  2nd edition, AMS Providence RI, 2010.
  * [14] A.P. Gibson, J.C. Hebden, and S.R. Arridge,  Recent advances in diffuse optical imaging.  Phys. Med. Biol., 50:R1-R43, (2005)
  * [15] D. Gintides, M. Sini and N. Thanh,  Detection of point-like scatterers using one type of scattered elastic waves.  J. Comp. App. Math., 236, (2012), 2137–2145.
  * [16] G. Granados and I. Harris,  Reconstruction of small and extended regions in EIT with a Robin transmission condition.  Inverse Problems, 38, (2022), 105009.
  * [17] R. Griesmaier,  Multi-frequency orthogonality sampling for inverse obstacle scattering problems.  Inverse Problems, 27, (2011), 085005.
  * [18] M. Hanke,  A note on the MUSIC algorithm for impedance tomography.  Inverse Problems, 33, (2017), 025001.
  * [19] I. Harris,  Regularization of the Factorization Method applied to diffuse optical tomography.  Inverse Problems, 37, (2021), 125010.
  * [20] I. Harris and D-L. Nguyen,  Orthogonality Sampling Method for the Electromagnetic Inverse Scattering Problem.  SIAM J. Sci. Comp., 42:3, (2020), B722-B737.
  * [21] I. Harris, D.-L. Nguyen and T.-P. Nguyen,  Direct sampling methods for isotropic and anisotropic scatterers with point source measurements.  Inverse Problems and Imaging, 16(5), (2022), 1137–1162.
  * [22] K. Ito, B. Jin, and J. Zou,  A direct sampling method to an inverse medium scattering problem.  Inverse Problems, 28, (2012), 025003.
  * [23] X. Ji and X. Liu,  Identification of Point-Like Objects with Multifrequency Sparse Data.  SIAM J. Sci. Comput., 42:4, (2020), A2325–A2343.
  * [24] A. Kirsch A and N. Grinberg,  “The Factorization Method for Inverse Problems”.  1st edition Oxford University Press, Oxford 2008.
  * [25] A. Kirsch,  The MUSIC-algorithm and the factorization method in inverse scattering theory for inhomogeneous media.  Inverse Problems, 18, (2002), 1025–1040.
  * [26] R. Kress and W. Rundell,  Reconstruction of extended sources for the Helmholtz equation.  Inverse Problems, 29, (2013), 035005.
  * [27] A. Lechleiter,  The MUSIC algorithm for impedance tomography of small inclusions from discrete data.  Inverse Problems, 31, (2015), 095004.
  * [28] X. Liu,  A novel sampling method for multiple multiscale targets from scattering amplitudes at a fixed frequency.  Inverse Problems, 33, (2017), 085011.
  * [29] X. Liu, S. Meng and B. Zhang,  Modified sampling method with near field measurements.  SIAM J. App. Math., 82:1 (2022) 244–266.
  * [30] D-L. Nguyen,  Direct and inverse electromagnetic scattering problems for bi-anisotropic media.  Inverse Problems, 35, (2019), 124001.
  * [31] W. Park,  Asymptotic properties of MUSIC-Type Imaging in Two-Dimensional Inverse Scattering from Thin Electromagnetic Inclusions.  SIAM J. Appl. Math., 75:1, (2015), 209–228.
  * [32] W. Park and D. Lesselier,  MUSIC-type imaging of a thin penetrable inclusion from its multi-static response matrix.  Inverse Problems, 25, (2009), 075002.
  * [33] K. Ren and Y. Zhong,  Imaging point sources in heterogeneousenvironments.  Inverse Problems, 35, (2019), 125003.
  * [34] T. Yin, G. Hu and L. Xu,  Near-field Imaging Point-like Scatterers and Extended Elastic Solid in a Fluid.  Commun. Comput. Phys., 19:5, (2016), 1317–1342.
  * [35] D. Zhang and Y. Guo,  Fourier method for solving the multi- frequency inverse source problem for the Helmholtz equation.  Inverse Problems, 31, (2015), 035007.
  * [36] D. Zhang, Y. Guo, J. Li and H. Liu  Locating Multiple Multipolar Acoustic Sources Using the Direct Sampling Method.  Commun. Comput. Phys., 25:5, (2019), 1328–1356.


# Signs behaviour of sums of weighted numbers of partitions

Filip Gawron and Maciej Ulas

###### Abstract.

Let $A$ be a subset of positive integers. By $A$-partition of $n$ we
understand the representation of $n$ as a sum of elements from the set $A$.
For given $i,n\in\mathbb{N}$, by $c_{A}(i,n)$ we denote the number of
$A$-partitions of $n$ with exactly $i$ parts. In the paper we obtain several
result concerning sign behaviour of the sequence
$S_{A,k}(n)=\sum_{i=0}^{n}(-1)^{i}i^{k}c_{A}(i,n)$, where $k\in\mathbb{N}$ is
fixed. In particular, we prove that for a broad class $\mathcal{A}$ of subsets
of $\mathbb{N}_{+}$ we have that for each $A\in\mathcal{A}$ we have
$(-1)^{n}S_{A,k}(n)\geq 0$ for each $n,k\in\mathbb{N}$.

###### Key words and phrases:

sums involving partitions; sign

Research of the authors was supported by a grant of the National Science
Centre (NCN), Poland, no. UMO-2019/34/E/ST1/00094

##  1\. Introduction

Let $\mathbb{N}$ and $\mathbb{N}_{+}$ denotes the set of non-negative integers
and the set of positive integers, respectively. Let $A\subset\mathbb{N}_{+}$
be given and take $n\in\mathbb{N}$. By an $A$-partition
$\lambda=(\lambda_{1},\ldots,\lambda_{k})$ with parts in $A$, of a non-
negative integer $n$, we mean a representation of $n$ in the form

$n=\lambda_{1}+\ldots+\lambda_{k},$

where $\lambda_{i}\in A$. The representations of $n$ differing only in order
of the terms are counted as one. By $p_{A}(n)$ we denote the number of all
$A$-partitions of $n$. We also use the standard convention $p_{A}(0)=1$. The
famous partition function $p(n):=p_{\mathbb{N}_{+}}(n)$ was introduced by
Euler. However, a true explosion of works dedicated to the sequence
$(p(n))_{n\in\mathbb{N}}$ was initiated by Ramanujan and his collaborations
with Hardy. The literature concerning various arithmetic properties of $p(n)$
and others partition functions associated with specific subsets of
$\mathbb{N}_{+}$ is enormous. The standard reference covering various aspects
of the theory partitions is the monograph of Andrews [1].

The common starting point in investigations of the sequence
$(p_{A}(n))_{n\in\mathbb{N}}$ is the ordinary generating function

$\sum_{n=0}^{\infty}p_{A}(n)x^{n}=\prod_{a\in A}\frac{1}{1-x^{a}}.$

The above identity can be seen twofold: as a identity in the formal power
series ring and also as a equality of two analytic functions (in the disc
$\\{z\in\mathbb{C}:\;|z|<1\\}$) represented by corresponding sides.

For given $A\subset\mathbb{N}_{+}$ let us write

$F_{A}(t,x)=\prod_{a\in
A}\frac{1}{1-tx^{a}}=\sum_{n=0}^{\infty}f_{A,n}(t)x^{n}$

and observe that the polynomial $f_{A,n}(t)$ is a natural generalization of
the number $p_{A}(n)=f_{A,n}(1)$ which counts the number of $A$-partitions of
the number $n$.

Let us write

$f_{A,n}(t)=\sum_{i=0}^{n}c_{A}(i,n)t^{i},$

and note that $c_{A}(i,n)$ is the number of $A$-partitions of $n$ with exactly
$i$ parts. In particular $f_{A,n}(1)=\sum_{i=0}^{n}c_{A}(i,n)$ is the number
of $A$-partitions of $A$. In the sequel, we will call the polynomial
$f_{A,n}(t)$ as a $n$-th $A$-partition polynomial.

In a recent paper [3] the case of $A$-partitions polynomials for
$A=\\{2^{i}:\;i\in\mathbb{N}\\}$ were investigated. Among many properties of
the $A$-partition polynomials, the authors proved that for $k=0,1$, the
expression $(-1)^{n}\sum_{i=0}^{n}(-1)^{i}i^{k}c_{A}(i,n)$ is non-negative.
Moreover, the conjecture was stated that for each $k\in\mathbb{N}$ this
expression is non-negative. This suggest the following general question.

######  Question 1.1.

Let $A\subset\mathbb{N}_{+}$ and $k\in\mathbb{N}$ be fixed. What can be said
about the sign behaviour of the sequence

$S_{A,k}(n)=\sum_{i=0}^{n}(-1)^{i}i^{k}c_{A}(i,n)?$

In particular, under which conditions on $A$ we have that
$(-1)^{n}S_{A,k}(n)\geq 0$ for given $k$ and all but finitely many $n$?

Let us also note that the value $S_{A,0}(n)=\sum_{i=0}^{n}(-1)^{i}c_{A}(i,n)$
has a nice combinatorial interpretation as the difference between the number
of partitions of $n$ with even number of parts and odd number of parts. One
can think that in general the sign behaviour of $S_{A,0}(n)$ should be
alternating, that means that $(-1)^{n}S_{A,0}(n)\geq 0$. The aim of this paper
is to investigate this expectation and prove several results which in many
cases allow the proof of the non-negativity of $(-1)^{n}S_{A,k}(n)$ for a
broad class of sets $A$ and all $k,n\in\mathbb{N}$. However, as we will see
this is not true in general. Our investigations fits in various studies
concerning behaviour of various counting objects involving parts of partitions
(see [2]).

Let us describe the content of the paper in some details. In Section 2 we
collect some basic observations concerning the sequence
$(S_{A,k}(n))_{n\in\mathbb{N}}$. In particular, by investigating so called
$\delta$ operator we introduce two family of polynomials
$p_{A,n}(t),q_{A,n}(t)$ closely related to the logarithmic derivative of the
generating function of $A$-partition polynomials and obtain recurrence
relations satisfied by terms of the double sequence
$(S_{A,k}(n))_{k,n\in\mathbb{N}}$. In Section 3 we obtain several results
which allow the proof of nonnegativity of $S_{A,k}(n)$ for each $k$ and $n$
and several sets $A$. In particular, if $A$ contains only odd numbers then
$(-1)^{n}S_{A,k}(n)\geq 0$. The same statement is true for the set
$A=\mathbb{N}_{+}$ and many others. Finally, in the last section we state
several problems and conjectures which appeared during our investigations and
we were unable to solve. We hope that this collection of problems will
stimulate further activity in this area.

##  2\. Basic observations

In order to investigate the sign behaviour of the sequence
$(S_{A,k}(n))_{n\in\mathbb{N}}$, where $k\in\mathbb{N}$ is fixed, we consider
the differential operator $\delta:\;\mathbb{R}[t]\rightarrow\mathbb{R}[t]$
defined in the following way: for $f\in\mathbb{R}[t]$ we put

$\delta(f(t))=t\frac{df(t)}{dt}.$

For $k=0$ we put $\delta^{(0)}(f)=f$ and for $k\in\mathbb{N}_{+}$ we define by
induction a $k$-th power of the operator $\delta$ as $\delta^{(k)}(f)$, where

$\delta^{(k)}(f)=\delta(\delta^{(k-1)}(f)).$

The basic properties of the operator $\delta$ can be summarized in the
following well known

######  Lemma 2.1.

  1. (1)

$\delta$ is linear operator, i.e., for all
$u,v\in\mathbb{R},f,g\in\mathbb{R}[t]$ we have
$\delta(uf+vg)=u\delta(f)+v\delta(g)$;

  2. (2)

if $f\in\mathbb{R}[t]\setminus\mathbb{R}$ then
$\operatorname{deg}\delta(f)=\operatorname{deg}f$;

  3. (3)

for $k\in\mathbb{N}$ we have

$\delta^{(k)}(fg)=\sum_{i=0}^{k}\binom{k}{i}\delta^{(i)}(f)\delta^{(k-i)}(g).$

  4. (4)

if $f(t)=\sum_{i=0}^{n}a_{i}t^{i}\in\mathbb{R}[t]$ and $k\in\mathbb{N}$ is
given then

$\delta^{(n)}(f(t))=\sum_{i=0}^{n}a_{i}i^{n}t^{i}$

###### Proof.

The first two properties are easy consequences of the definition of the
operator $\delta$. The proof of the third property is just the application of
the Leibnitz rule together with induction on $k$. The last property follows
from the linearity of the operator $\delta$ and the equality
$\delta^{(k)}(t^{i})=i^{k}t^{i}$. We omit the simple details. ∎

As an immediate consequence of the fourth property of the operator $\delta$ we
get the equality

$S_{A,k}(n)=\sum_{i=0}^{n}(-1)^{i}i^{k}c_{A}(i,n)=\delta^{(k)}(f_{A,n}(t))\mid_{t=-1}.$

In order to simplify the notation for a given $A\subset\mathbb{N}_{+}$ and
$n\in\mathbb{N}_{+}$ we write

$A(n):=\\{d\in A:\;d|n\\},$

i.e., the set $A(n)$ contains those divisors of $n$ which are lying in $A$.

We prove certain identities which will be useful in our investigations and
allow to express action of $\delta$ operator on $A$-partition polynomial in
terms of $A$-partitions polynomials and a related sequence of polynomials.
More precisely, we prove the following.

######  Lemma 2.2.

Let $n\in\mathbb{N}$. Then

$\displaystyle\sum_{i=1}^{n}p_{A,i}(t)f_{A,n-i}(t)$
$\displaystyle=tf_{A,n}^{\prime}(t)=\delta(f_{A,n}(t)),$
$\displaystyle\sum_{i=1}^{n}q_{A,i}(t)f_{A,n-i}(t)$
$\displaystyle=nf_{A,n}(t),$ $\displaystyle tq^{\prime}_{A,n}(t)$
$\displaystyle=np_{A,n}(t),$

where

$p_{A,n}(t)=\sum_{a\in A(n)}t^{n/a},$

and

$q_{A,n}(t)=\sum_{a\in A(n)}at^{n/a}.$

###### Proof.

Applying logarithmic differentiation with respect to the variable $t$, we get
the following expansion

$\frac{t}{F_{A}(t,x)}\frac{\partial F_{A}(t,x)}{\partial
t}=t\frac{\partial}{\partial t}(\log F_{A}(t,x))=\sum_{a\in
A}\frac{tx^{a}}{1-tx^{a}}=\sum_{n=1}^{\infty}p_{A,n}(t)x^{n}.$

Thus, the first identity is a consequence of the equality

$t\frac{\partial F_{A}(t,x)}{\partial
t}=\left(\sum_{n=1}^{\infty}p_{A,n}(t)x^{n}\right)F_{A}(t,x)=\sum_{n=0}^{\infty}\left(\sum_{i=1}^{n}p_{A,i}(n)f_{A,n-i}(t)\right)x^{n}$

by a comparison of coefficients near $x^{n}$ on far ends.

Similarly, applying logarithmic differentiation with respect to the variable
$x$, we get

$\frac{x}{F_{A}(t,x)}\frac{\partial F_{A}(t,x)}{\partial
x}=x\frac{\partial}{\partial x}(\log F_{A}(t,x))=\sum_{a\in
A}\frac{atx^{a}}{1-tx^{a}}=\sum_{n=1}^{\infty}q_{A,n}(t)x^{n},$

and thus, form the equality

$x\frac{\partial F_{A}(t,x)}{\partial
x}=\left(\sum_{n=1}^{\infty}q_{A,n}(t)x^{n}\right)F_{A}(t,x)=\sum_{n=0}^{\infty}\left(\sum_{i=1}^{n}q_{A,i}(n)f_{A,n-i}(t)\right)x^{n}$

and comparison of coefficients, we get the second identity form the statement.

The third identity is a direct consequence of the from of the polynomials
$p_{A,n},q_{A,n}$.

∎

######  Remark 2.3.

Let us note that the polynomial $p_{A,n}(t)$ can be seen as a polynomial
generalization of the number of $A$-divisors function $\tau_{A}(n)=\sum_{d\in
A(n)}1=p_{A,n}(1)$. Similarly, the polynomial $q_{A,n}(t)$ can be seen as a
polynomial analogue of the sum of $A$-divisors function
$\sigma_{A}(n)=\sum_{d\in A(n)}d=q_{A,n}(1)$.

As an application of the above result we prove a double recurrence relation
satisfied by the sequence $(S_{A,k}(n))_{k,n}$. More precisely, we have the
following.

######  Corollary 2.4.

Let $k\in\mathbb{N}$. We have the following identities

$S_{A,k+1}(n)=\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\delta^{(j)}(p_{A,i}(t))\mid_{t=-1}S_{A,k-j}(n-i).$

and

$\delta^{(k)}(p_{A,n}(t))=\sum_{a\in
A(n)}\left(\frac{n}{a}\right)^{k}t^{\frac{n}{a}}.$

###### Proof.

As a consequence of the first identity from Lemma 2.2 we get the following
chain of identities

$\displaystyle S_{A,k+1}(n)$
$\displaystyle=\delta^{(k)}(\delta(f_{A,n}(t)))=\delta^{(k)}(tf_{A,n}^{\prime}(t))\mid_{t=-1}$
$\displaystyle=\delta^{(k)}\left(\sum_{i=1}^{n}p_{A,i}(t)f_{A,n-i}(t)\right)\mid_{t=-1}$
$\displaystyle=\sum_{i=1}^{n}\delta^{(k)}(p_{A,i}(t)f_{A,n-i}(t))\mid_{t=-1}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\delta^{(j)}(p_{A,i}(t))\mid_{t=-1}\delta^{(k-j)}(f_{A,n-i}(t))\mid_{t=-1}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\delta^{(j)}(p_{A,i}(t))\mid_{t=-1}S_{A,k-j}(n-i)$

and hence the result.

The second identity from the statement is a simple consequence of the last
property from Lemma 2.1.

∎

##  3\. Results

In this section we prove several results which allow to prove non-negativity
of $(-1)^{n}S_{A,k}(n)$ for a broad class of sets $A$. Moreover, we show that
the sequence $(\operatorname{sign}(S_{A,k}(n)))_{n\in\mathbb{N}}$, where
$A=\\{a,b\\},a<b,\gcd(a,b)=1$, is periodic of period $2a(b-a)$.

We start with the following result.

######  Lemma 3.1.

Let $A_{1},A_{2}\subset\mathbb{N}_{+}$ and suppose that $A_{1}\cap
A_{2}=\emptyset$. If for each $k\in\mathbb{N}$ we have
$(-1)^{n}S_{A_{i},k}(n)\geq 0$ for $i=1,2$, then $(-1)^{n}S_{A_{1}\cup
A_{2},k}(n)\geq 0$.

###### Proof.

We are interesting in the value of

$S_{A_{1}\cup A_{2},k}(n)=\delta^{(k)}(f_{A_{1}\cup A_{2},n}(t))\mid_{t=-1}.$

Because $A_{1}\cap A_{2}=\emptyset$ we have the identity $F_{A_{1}\cup
A_{2}}(t,x)=F_{A_{1}}(t,x)F_{A_{2}}(t,x)$. As a consequence we deduce that

$f_{A_{1}\cup A_{2},n}(t)=\sum_{i=0}^{n}f_{A_{1},i}(t)f_{A_{2},n-i}(t).$

Thus, from Lemma 2.1 we obtain

$\delta^{(k)}(f_{A_{1}\cup
A_{2},n}(t))=\sum_{i=0}^{n}\sum_{j=0}^{k}\binom{k}{j}\delta^{(j)}(f_{A_{1},i}(t))\delta^{(k-j)}(f_{A_{2},n-i}(t)).$

Finally, we get

$\displaystyle(-1)^{n}S_{A_{1}\cup A_{2},k}(n)$
$\displaystyle=(-1)^{n}\delta^{(k)}(f_{A_{1}\cup A_{2},n}(t))\mid_{t=-1}$
$\displaystyle=\sum_{i=0}^{n}\sum_{j=0}^{k}\binom{k}{j}\\{(-1)^{i}S_{A_{1},j}(i)\\}\\{(-1)^{n-i}S_{A_{2},k-j}(n-i)\\}\geq
0$

and from our assumptions we get the result. ∎

######  Corollary 3.2.

Let $\mathcal{A}=\\{A_{1},\ldots,A_{m}\\}$ be a possible infinite family,
i.e., we allow $m=+\infty$, of pairwise disjoint subsets of $\mathbb{N}_{+}$
and suppose that for each $k,n\in\mathbb{N}$ we have
$(-1)^{n}S_{A_{i},k}(n)\geq 0$ for $i=1,\ldots,m$. If $A=\cup_{i=1}^{m}A_{i}$,
then for each $k,n\in\mathbb{N}$ we have $(-1)^{n}S_{A,k}(n)\geq 0$.

###### Proof.

If $m<+\infty$ then the statement follows from the simple induction on $m$
using Lemma 3.1 with $A=A_{1}\cup\ldots\cup A_{m-1},B=A_{m}$. We omit the
simple details.

If $m=+\infty$ we proceed as follows. Let $n\in\mathbb{N}_{+}$. Because the
family $\mathcal{A}$ is infinite and contains pairwise-disjoint sets there is
an integer $s$ such that $n<\operatorname{min}A_{i}$ for each $i>s$. Because
for $(-1)^{n}S_{A_{i},k}(n)\geq 0$ for each $i\leq s$ the same is true for the
set $\bigcup_{i=1}^{s}A_{i}$. Next, due to the fact that $n$ can not be
represented as a sum of elements from the set $\bigcup_{i=m+1}^{\infty}A_{i}$
we get that $(-1)^{n}S_{B,k}(n)\geq 0$ \- and hence the result.

∎

######  Lemma 3.3.

Let $A\subset 2\mathbb{N}+1$. Then, for each $k\in\mathbb{N}$ we have
$(-1)^{n}S_{A,k}(n)\geq 0$.

###### Proof.

We proceed by induction on $k\in\mathbb{N}$. The generating function for the
sequence $(S_{A,0}(n))_{n\in\mathbb{N}}$ is

$F_{A}(-1,x)=\sum_{n=0}^{\infty}S_{A,0}(n)x^{n}=\prod_{a\in
A}\frac{1}{1+x^{a}}=\prod_{a\in
A}\frac{1}{1-(-x)^{a}}=\sum_{n=0}^{\infty}(-1)^{n}p_{A}(n)x^{n},$

i.e., $S_{A,0}(n)=(-1)^{n}p_{A}(n)$ and our statement holds for $k=0$. Let us
suppose that our statement is true up to fixed $k\in\mathbb{N}$. We prove it
for $k+1$. For $j\in\\{0,\ldots,k\\}$ we have that

$\delta^{(j)}(p_{A,n}(t))\mid_{t=-1}=\sum_{a\in
A(n)}\left(\frac{n}{a}\right)^{j}(-1)^{\frac{n}{a}}=\begin{cases}\begin{array}[]{ll}-\sum_{a\in
A(n)}\left(\frac{n}{a}\right)^{j},&n\equiv 1\mkern 4.0mu({\operator@font
mod}\mkern 6.0mu2)\\\ \sum_{a\in A(n)}\left(\frac{n}{a}\right)^{j},&n\equiv
0\mkern 4.0mu({\operator@font mod}\mkern 6.0mu2).\end{array}\end{cases}$

Thus, for $n\in\mathbb{N}$, we have that
$(-1)^{n}\delta^{(j)}(p_{A,n}(t))\mid_{t=-1}\geq 0$. To get the statement for
$k+1$ we use Corollary 2.4 and note that

$\displaystyle(-1)^{n}S_{A,k+1}$
$\displaystyle(n)=(-1)^{n}\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\delta^{(j)}(p_{A,i}(t))\mid_{t=-1}S_{A,k-j}(n-i)$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\\{(-1)^{i}\delta^{(j)}(p_{A,i}(t))\mid_{t=-1}\\}\\{(-1)^{n-i}S_{A,k-j}(n-i)\\}\geq
0.$

The first expression in curly parentheses is non-negative due to our
observation above and the second expression in curly parentheses is non-
negative from our induction hypothesis. Hence, the statement is true for $k+1$
and our theorem is proved.

∎

We generalize the above result and prove the following:

######  Theorem 3.4.

Let $A\subset 2\mathbb{N}+1,m\in\mathbb{N}_{+}$ and ${\bf p}=(p_{i})_{i\leq
m}$ be an increasing possibly infinite sequence of non-negative integers with
$p_{0}=0$, i.e., we allow $m=+\infty$. Let us put

$B=\bigcup_{j=0}^{m}2^{p_{j}}A.$

Then for each $n,k\in\mathbb{N}$ we have $(-1)^{n}S_{B,k}(n)\geq 0$.

###### Proof.

First let us note that

$B(n)=\\{b\in
B:\;b|n\\}=\left\\{b\in\bigcup_{j=0}^{m}2^{p_{j}}A:\;b|n\right\\}=\bigcup_{j=0}^{m}A(n/2^{p_{j}}).$

Although we allow $m=+\infty$, for any given $n\in\mathbb{N}$, the sum on the
right side is finite. Indeed, if $p_{j}>\nu_{2}(n)$ then
$A(n/2^{p_{j}})=\emptyset$. To get the result we consider two cases:
$m<+\infty$ and $m=+\infty$.

Suppose that $m<+\infty$. Let us define $q_{i}=p_{i+1}-p_{i}$. From our
assumption on the sequence ${\bf p}$ we know that $q_{i}\geq 1$ for
$i=0,\ldots,m-1$.

We note that the sign of $S_{B,0}(n)=f_{B,n}(-1)$ is the sign of the $n$-th
coefficient of

$\displaystyle\prod_{a\in B}\frac{1}{1+x^{a}}$ $\displaystyle=\prod_{a\in
A}\prod_{j=0}^{m}\frac{1}{1+x^{2^{p_{j}}a}}=\prod_{a\in
A}\prod_{j=0}^{m}\frac{1-x^{2^{p_{j}}a}}{1-x^{2^{p_{j}+1}a}}=\prod_{a\in
A}\left(\frac{1}{1-x^{2^{p_{m}+1}a}}\right)F(x),$

where

$\displaystyle F(x)=$ $\displaystyle\prod_{a\in
A}\left(\left(1-x^{2^{p_{0}}a}\right)\left(\prod_{i=0}^{m-1}\frac{1-x^{2^{p_{i+1}}a}}{1-x^{2^{p_{i}+1}a}}\right)\right)$
$\displaystyle\prod_{a\in
A}\left(\left(1-x^{2^{p_{0}}a}\right)\left(\prod_{i=0}^{m-1}\prod_{j=1}^{q_{i}-1}\left(1+x^{2^{p_{i}+j}a}\right)\right)\right)$
$\displaystyle\prod_{a\in
A}\left(\left(1+(-x)^{2^{p_{0}}a}\right)\left(\prod_{i=0}^{m-1}\prod_{j=1}^{q_{i}-1}\left(1+(-x)^{2^{p_{i}+j}a}\right)\right)\right).$

Now let us note that $F(-x)$ is a generating function for the partition
function, say $pd_{A}({\bf p},n)$, which counts partitions of $n$ into
distinct parts from the set

$\bigcup_{i=0}^{k}\bigcup_{j=1}^{q_{i}-1}2^{p_{i}+j}A.$

Consequently we get that

$\displaystyle\prod_{a\in B}\frac{1}{1+x^{a}}$
$\displaystyle=\left(\sum_{n=0}^{\infty}p_{A}(n)x^{2^{p_{k}+1}a}\right)\left(\sum_{n=0}^{\infty}(-1)^{n}pd_{A}({\bf
p},n)x^{n}\right)$ $\displaystyle=\sum_{n=0}^{\infty}\left(\sum_{i=0}^{\lfloor
n/2^{p_{k}+1}\rfloor}(-1)^{n-2^{p_{k}+1}i}p_{A}(i)pd_{A}({\bf
p},n-2^{p_{k}+1}i)\right)x^{n}$
$\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}\left(\sum_{i=0}^{\lfloor
n/2^{p_{k}+1}\rfloor}p_{A}(i)pd_{A}({\bf p},n-2^{p_{k}+1}i)\right)x^{n}$

and it is clear that $(-1)^{n}S_{B,0}(n)\geq 0$.

Next observe that for $n$ odd we have

$\displaystyle\delta^{(i)}(p_{B,n}(-1))=\sum_{a\in
B(n)}\left(\frac{n}{a}\right)^{i}(-1)^{\frac{n}{a}}.$

and this sum is smaller than zero because $(-1)^{n/a}=-1$. If $n$ is even,
i.e.,. $n=2^{\nu_{2}(n)}m$, where $\nu_{2}(n)$ is the highest power of 2 which
divides $n$, then we can group elements of $B_{2^{\nu_{2}(n)}m}$ into
$M+1$-tuples $(2^{p_{0}},a,2^{p_{1}}a,\ldots,2^{p_{M}}a)$ with $a$ odd and

$M=\operatorname{max}\\{i\in\\{1,\ldots,k\\}:\;p_{i}\leq\nu_{2}(n)\\}.$

We thus have

$\displaystyle\delta^{(i)}(p_{B,n}(-1))$ $\displaystyle=\sum_{a\in
A(m)}\left(\sum_{s=0}^{M}\left(\frac{n}{2^{p_{s}}a}\right)^{i}(-1)^{\frac{n}{2^{p_{s}}a}}\right).$

It is clear that for any given $n$, the only possible value of
$s\in\\{0,\ldots,M\\}$ for which $\frac{n}{2^{p_{s}}a}$ may be an odd integer
is $s=M$ and additionally the equality $p_{M}=\nu_{2}(n)$ need to be
satisfied. We thus get that for each $i,n\in\mathbb{N}$ we have the inequality
$(-1)^{n}\delta^{(i)}(p_{B,n}(t))\mid_{t=-1}\geq 0$.

We are ready to finish the proof that $(-1)^{n}S_{B,k}(n)\geq 0$. We proceed
by induction on $k$. We already proved that the statement is true for $k=0$.
Let us assume that it holds for $j=0,1,\ldots,k$. Therefore,

$\displaystyle(-1)^{n}S_{B,k+1}(n)$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=0}^{k}\binom{k}{j}\\{(-1)^{i}\delta^{(j)}(p_{B,i}(t))\mid_{t=-1}\\}\\{(-1)^{n-i}S_{B,k-j}(n-i)\\}.$

As we observed, the first expression in the curly bracket is non-negative. The
second is non-negative because of our induction assumption, and hence we get
the result.

To get the result for $m=+\infty$ we need the following observation. For a
given $A\subset\mathbb{N}_{+}$ and given $n\in\mathbb{N}$ we put

$U(n)=\\{a\in A:\;a\leq n\\}=A\cap[1,n].$

Let us observe that for a given $N\in\mathbb{N}_{+}$ and $n\leq N$ we have the
identity

$f_{A,n}(t)=f_{U(N),n}(t).$

In consequence, for a given $k\in\mathbb{N},n\leq N$, we have the equality
$S_{A,k}(n)=S_{U(N),k}(n)$.

As a consequence of the above property we see that if $m=+\infty$ and $n$ is
given, there is a $s\in\mathbb{N}_{+}$ such that $n<2^{p_{i}}$ for $i\geq s$.
Indeed, the sequence ${\bf p}$ is increasing. Hence, if we put ${\bf
p}^{\prime}=(p_{i})_{i\leq s}$ then we have the equality

$f_{A,n}(t)=f_{U(2^{p_{s}}),n}(t)$

and hance we can apply our statement for the finite set ${\bf p}^{\prime}$ and
get the required result.

∎

The combination of Lemma 3.3 and Theorem 3.4 is a very useful tool which can
be used to prove sign alternating property of the sequence
$(S_{A,k}(n))_{n\in\mathbb{N}}$ for several sets $A$. We start with the
following.

######  Theorem 3.5.

Let $m\in\mathbb{N}_{+}$ and put $A=\\{1,2,\ldots,m\\}$. Then, for each
$k,n\in\mathbb{N}$ we have $(-1)^{n}S_{A,k}(n)\geq 0$.

###### Proof.

To get the result it is enough to observe the equality of sets

$A=\\{1,\ldots,m\\}=\bigcup_{i=1}^{\lceil
m/2\rceil}(2i-1)\\{1,2,\ldots,2^{\lfloor\log_{2}(m/(2i-1))\rfloor}\\}$

and apply Theorem 3.4 together with Lemma 3.2. ∎

######  Corollary 3.6.

Let $A=\mathbb{N}_{+}$. Then for each $k\in\mathbb{N}$ we have
$(-1)^{n}S_{A,k}(n)\geq 0$.

###### Proof.

To get the result it is enough to note the equality

$\mathbb{N}_{+}=\bigcup_{i=0}^{\infty}(2\mathbb{N}+1)2^{i}$

and apply Theorem 3.4 with $p_{i}=i$ and $A=2\mathbb{N}+1$-the set of odd
numbers. ∎

######  Corollary 3.7.

Let $m\in\mathbb{N}_{\geq 2}$ be odd or be a power of two. Let us put
$B=\\{m^{i}:\;i\in\mathbb{N}\\}$. Then for each $k\in\mathbb{N}$ we have
$(-1)^{n}S_{B,k}(n)\geq 0$.

###### Proof.

If $m$ is odd then we directly apply Lemma 3.3. If $m=2^{s}$ for
$s\in\mathbb{N}_{+}$ then we apply the decomposition into pairwise disjoint
sets

$B=\bigcup_{j=0}^{\infty}\\{1\\}2^{js}.$

and apply Theorem 3.4 with $p_{j}=js$ for $j=1,2\ldots$ and $A=\\{1\\}$. ∎

######  Remark 3.8.

Let us note that if $m=2$ then we get that the [3, Conjecture 5.1] is true.

Our results shows that the alternating behaviour of the signs of the sums
$S_{A,k}(n)$ is quite typical. This may suggest that for all sets $A$, each
$k\in\mathbb{N}$ and sufficiently large $n\in\mathbb{N}$ we have
$(-1)^{n}S_{A,k}(n)\geq 0$. This is clearly not the case. We prove the
following.

######  Theorem 3.9.

Let $k\in\mathbb{N},a,b\in\mathbb{N}_{+}$ and take $A=\\{a,b\\}$ with $a<b$
and $(a,b)=1$. If $k=0$ then the sequence
$(\operatorname{sign}(S_{A,k}(n)))_{n\in\mathbb{N}}$ is eventually periodic of
period $2ab$.

If $k>0$ then the sequence
$(\operatorname{sign}(S_{A,k}(n)))_{n\in\mathbb{N}}$ is eventually periodic of
period $2a(b-a)$.

###### Proof.

Let us recall that $S_{A,k}(n)=\sum_{i=0}^{n}i^{k}c_{A}(i,n)$, where
$c_{A}(i,n)$ is the number of $A$-partitions of $n$ with exactly $i$-parts.
Because we have only two elements so we can calculate $c_{A}(i,n)$ directly as
the number of non-negative solutions $(x,y)$ of the linear system of
equations:

(1) $\begin{cases}\begin{array}[]{ll}xa+yb&=n\\\ x+y&=i\end{array}\end{cases}$

The necessary condition for solvability of (1) is the condition $ia\leq n\leq
ib$. Equivalently we have $\frac{n}{b}\leq i\leq\frac{n}{a}.$ Moreover, we
have

$\displaystyle n$ $\displaystyle=xa+yb=(i-y)a+yb=y(b-a)+ia,$ $\displaystyle y$
$\displaystyle=\frac{n-ia}{b-a}.$

Consequently, the solution exists if and only if $i\equiv n/a\equiv
r(n)\mod(b-a)$ and $\frac{n}{b}\leq i\leq\frac{n}{a}$. Note that, the
congruence $ai\equiv n\mkern 4.0mu({\operator@font mod}\mkern 6.0mub-a)$ has
always a solutions under due to the fact that $\gcd(a,b-a)=1|n$. Thus, we get
the formula

$c_{A}(i,n)=\begin{cases}\begin{array}[]{ll}1&\textrm{if $i\in[n/b,n/a]$ and
$i\equiv r(n)\mod(b-a)$}\\\ 0&\textrm{in the remaining
cases}\end{array}.\end{cases}$

From our reasoning we obtain the expression for $f_{A,n}(t)$, i.e.,

$f_{A,n}(t)=\sum_{\begin{subarray}{c}j\in[n/b,n/a]\\\ j\equiv
r(n)\mod(b-a)\end{subarray}}t^{j}$

and

$\delta^{(k)}(f_{A,n})(t)=\sum_{\begin{subarray}{c}j\in[n/b,n/a]\\\ j\equiv
r(n)\mod(b-a)\end{subarray}}j^{k}t^{j}$

First let’s consider the case $k>0.$ Assume that $n$ is big enough such that
the segment $[n/b,n/a]$ contains at least $3(b-a)$ elements. Now substitute
$-1$ for $t.$ If $b-a$ is even then the parity of $j$ is the same as the
parity of $r(n)$ so every term $j^{k}(-1)^{j}$ has the same sign. In
particular, the sign of the whole sum is the same as the sign of the last
element in the sum. If $b-a$ is odd then we obtain an alternating sum where
every next term is bigger then the previous one, so the sign of this sum will
be also the sign of the last element in the sum. Hence

$\operatorname{sign}(\delta^{(k)}(f_{A,n})(-1))=\operatorname{sign}(j^{k}(-1)^{j})=\operatorname{sign}((-1)^{j}),$

where $j$ is the biggest number in the set $[n/b,n/a]$ satisfying the
condition $j\equiv r(n)\mkern 4.0mu({\operator@font mod}\mkern 6.0mu(b-a))$.
Now if we add $a(b-a)$ to $n$ then $n/a$ will increase by $b-a$ and $r(n)$
will not change. So our maximal $j$ will be now our previous maximal $j$ plus
$b-a$ and from our reasoning it is clear that the period of our sum is a
divisor of the number $2a(b-a).$

For $k=0$ looking at the last term is not sufficient, because if two divides
$\\#(\mathbb{Z}\cap[n/a,n/b])$ then our sum can be zero. But that is the only
difference. As in $k>0$ case, if $b-a$ is even then the parity of our summand
$j$ is the same and the sign of the whole sum is the same as the sign of the
last element in the sum. If $b-a$ is odd and $2|\\#(\mathbb{Z}\cap[n/a,n/b])$
then our sum is zero. Otherwise it has the same sign as the last element. Now,
if we add $2ab$ to $n$ then the number $r(n)$, the parity of
$\\#(\mathbb{Z}\cap[n/a,n/b])$, and the parity of its last integer element
stay the same. Thus we get that the period of our sign sequence is a divisor
of the number $2ab$. ∎

##  4\. Some problems and conjectures

In this section we collect some problems and conjectures which appeared during
our investigations.

######  Conjecture 4.1.

For any given $k\in\mathbb{N}$ and finite $A\subset\mathbb{N}_{+}$ with
$\gcd(A)=1$, the sequence $(S_{A,k}(n))_{n\in\mathbb{N}}$ is eventually
periodic.

Let us consider the set $A=\\{1,2,6\\}$. One can show that the sequence of
signs of $S_{A,k}(n)$ is eventually periodic with period 4 for $k=0,1$. More
precisely,

$\displaystyle(-1)^{\left\lfloor\frac{n+1}{2}\right\rfloor}S_{0,A}(n)>0\quad\mbox{for}\quad
n\geq 3,$
$\displaystyle(-1)^{\left\lfloor\frac{n+1}{2}\right\rfloor}S_{1,A}(n)>0\quad\mbox{for}\quad
n\geq 8.$

This can be easily proved from the explicit expressions for the sequences
$(S_{A,i}(n))_{n\in\mathbb{N}}$ for $i=0,1$. Indeed, if $i=0$, then
$S_{A,0}(n)=f_{A,n}(-1)$ is just the $n$-th coefficient in the power series
expansion of the rational function $F_{A}(-1,x)=1/(1+x)(1+x^{2})(1+x^{6})$.
Similarly, if $i=1$ we need to investigate coefficients of the power series
expansion of the rational function

$\sum_{n=0}^{\infty}S_{A,1}(n)x^{n}=\partial_{t}F(t,x)\mid_{t=-1}=-\frac{x\left(3x^{6}+2x^{5}-x^{4}-x^{3}+x^{2}+x+1\right)}{(1+x)^{2}(1+x^{2})(1+x^{6})^{2}}.$

However, one can observe that $(-1)^{n}S_{A,k}(n)>0$ for $k\geq 2$ and $n\gg
0$. This and related numerical observations suggest the following

######  Problem 4.2.

Let $A\subset\mathbb{N}_{+}$ and suppose that for some $k\in\mathbb{N}$ we
have $(-1)^{n}S_{A,k}(n)\geq 0$ for $n\gg 0$. Then for each $i\in\mathbb{N}$
we have $(-1)^{n}S_{A,k+i}(n)\geq 0$ for $n\gg 0$.

During our investigations we obtained several results concerning periodic
behaviour of the sequence of signs of the sequence
$(S_{A,k}(n))_{n\in\mathbb{N}}$. It is an interesting question whether we can
find a set $A\subset\mathbb{N}_{+}$ such that the sign behaviour of the
corresponding sequence (for some $k\in\mathbb{N}$) is not periodic.

######  Conjecture 4.3.

Let $A=\\{3n+1:\;n\in\mathbb{N}\\}$. The sequence
$(\operatorname{sign}(f_{A,n}(-1)))_{n\in\mathbb{N}}$ is not periodic.

We expect that the following is true

######  Conjecture 4.4.

Let $m\in\mathbb{N}_{\geq 2}$ and take $U_{m}=\\{i!:\;i=1,\ldots,m\\}$. Then,
for each $k\in\\{0,\ldots,m-2\\}$

$(-1)^{\left\lfloor\frac{n+1}{2}\right\rfloor}S_{U_{m},k}(n)\geq
0\quad\mbox{for}\quad n\gg 0,$

i.e., the sequence of signs of the sequence
$(S_{U_{m},k}(n))_{n\in\mathbb{N}}$ is eventually periodic of period 4.
Moreover, if $k\geq m-1$ then $(-1)^{n}S_{U_{m},k}(n)\geq 0$ for $n\gg 0$.

If the above conjecture is true then clearly implies that if
$A=\\{i!:\;i\in\mathbb{N}_{+}\\}$, then for each $k\in\mathbb{N}$, the
sequence $(\operatorname{sign}(S_{A,0}(n)))_{n\in\mathbb{N}}$ is eventually
periodic of period 4.

######  Problem 4.5.

Explain the behaviour of $(\operatorname{sign}(S_{A,k}(n))_{n\in\mathbb{N}}$
for $A=\\{m+1,\ldots,m+k\\}$, where $k$ is fixed.

Numerical calculations suggest that for $k\geq 3$ the behaviour is
alternating.

One can also go one step further and note that $S_{A,k}(n)$ is just the value
of the polynomial $\delta^{(k)}(f_{A,n}(t))$ at the simplest root of unity
$t=-1$. This suggest to investigate the following general

######  Problem 4.6.

Let $d\in\mathbb{N}_{\geq 3}$ and $\zeta_{d}$ be the $d$-th primitive root of
unity. Let us write

$\delta^{(k)}(f_{A,n}(t))\mid_{t=\zeta_{d}}=\sum_{i=0}^{\varphi(d)}u^{(d)}_{A,k}(i,n)\zeta_{d}^{i}.$

What can be said about the sign behaviour of the family of sequences
$(u^{(d)}_{A,k}(i,n))_{n\in\mathbb{N}}$, where $k\in\mathbb{N}$ and
$i=0,1,\ldots,n-1$ are fixed? Especially interesting is the case of $k=0$.

Numerical calculations suggest that we can not expect periodic behaviour even
in the basic case when $A=\mathbb{N}_{+}$. In fact, based on our calculations
we expect that there are arbitrary long sequences of consecutive positive and
negative values of $u^{(4)}_{A,0}(i,n)$ for $t=\zeta_{4}$, i.e.,
$\zeta_{4}^{2}+1=0$, $i=0,1$. On the other hand, if
$A=\\{2^{i}:\;i\in\mathbb{N}\\}$ then in [3, Corollary 4.2] the authors
proved, among other things, that

$\sum_{i=0}^{4n+2}(-1)^{i}u^{(4)}_{A,0}(2i,8n+4)=\sum_{i=0}^{4n+1}(-1)^{i}u^{(4)}_{A,0}(2i+1,8n+4)=0,$

which is a consequence of the vanishing of the value at $t=\zeta_{4}$ of
$8n+4$-th $A$-partition polynomial. Thus, a deeper study of Problem 4.6 may
reveals interesting combinatorial identities.

## References

  * [1] G. E. Andrews, The Theory of Partitions, Cambridge University Press, Cambridge, 1998.
  * [2] M. Merca, On the sum of parts in the partitions of $n$ into distinct parts, Bull. Aust. Math. Soc. 104, No. 2 (2021), 228–237.
  * [3] M. Ulas, B. Żmija, On arithmetic properties of binary partition polynomials, Adv. Appl. Math. 110 (2019), 153–179.

Maciej Ulas, Jagiellonian University, Faculty of Mathematics and Computer
Science, Institute of Mathematics, Łojasiewicza 6, 30-348 Kraków, Poland;
email<EMAIL_ADDRESS>
Filip Gawron, Jagiellonian University, Faculty of Mathematics and Computer
Science, Institute of Mathematics, Łojasiewicza 6, 30-348 Kraków, Poland;
email<EMAIL_ADDRESS>