\section{Introduction}
The direct and inverse medium scattering problems have been investigated extensively in the past few decades due to their important practical applications in geophysics, nondestructive testing, biological studies, evaluation and medicine, see, for instance, \cite{MSR, medical}.
A large variety of numerical reconstruction methods have been developed and applied in practice for different kinds of models. For instance, some effective methods are developed for the homogeneous background model: the linear sampling or probing methods (LSM) \cite{CK}, the time-reversal multiple signal classification method (MUSIC) \cite{iakov1,iakov2,music}, the direct sampling method (DSM) \cite{Jin12}, the multilevel sampling method (MSM) \cite{LZ}, the contrast source inversion method (CSIM) \cite{BerBro}, etc. Moreover, several efficient algorithms are applied in various practical situations, e.g., in the stratified ocean waveguide model: the simple method (SM) \cite{L1}, extended direct sampling methods (EDSMs) \cite{CILZ, DSM11, DSM12, DSM13, DSM14, DSM15, L2, LXZ2, LXZ3},  the extended MUSIC method (EMUSIC) \cite{AIK, extended1},  the topological based imaging  functional \cite{topological}, etc. Furthermore, several efficient approaches are introduced in the two-layered medium model: the extended multilevel sampling method (EMSM) \cite{YL}, the extended linear sampling method (ELSM) \cite{DLLY}, etc.

Most reconstruction methods involve an indicator function (called also imaging functional) to be evaluated at each sampling point inside a sampling region for the recovery of the shape (and, if possible, the parameters therein) of the unknown inclusions.
In practice, the sampling region is usually chosen to be large enough to contain all the unknown inclusions; then a considerably fine mesh is used in the sampling region to reconstruct the details of the objects. However, an extremely fine mesh would drastically increase the computational complexity in practice, especially in the three-dimensional case.  In fact, the computation complexity of most reconstruction methods (see, for instance, DSM, MUSIC, LSM, CSIM, etc.) is of $\mathcal{O}(N^2)$ in  two dimensions and $\mathcal{O}(N^3)$ in three dimensions, where $N$ is the number of sampling points along one direction.  Reconstruction methods such as MSM and the EMSM  has complexity of order of $N \log N$ in two dimensions and $N^2\log N$ in  three dimensions. Hence, they are still quite demanding if $N$ is large.
Therefore, an optimal mesh size for the sampling and the reconstruction is crucial to reduce the computation complexity of these algorithms and at the same time to reconstruct details of the inclusions in practical applications.
Moreover, several iterative methods (i.e., EMSM, CSIM, etc.) apply certain cut-off values as a stopping criterion of the method.  Such a choice of stopping criterion may heavily depend on the reconstruction quality under a coarser mesh and a subjective choice of the cutoff value; which sometimes causes the reconstruction algorithm to entirely miss some small inclusions.  Choosing an optimal mesh size can improve the quality of the reconstruction  and practically avoid the aforementioned issue.

To the best of the authors' knowledge, although there are quite a lot of discussions of resolution analysis (e.g., Abbe-Rayleigh resolution limit, Sparrow's resolution limit, etc) in inverse scattering problems, e.g., \cite{resol1, resol2, FWA, GP}, a systematic choice of the optimal mesh size for numerical reconstruction algorithms in inverse scattering problems has not yet been addressed in the literature.
We would like to remark that the resolution analysis aims at understanding the quality of an image given a fixed grid size; whereas here we aim at choosing an optimal mesh size to ensure certain quality of the reconstruction.
This paper is hence devoted to an optimal selection of the mesh size in reconstruction methods. It is worth emphasizing that, in this paper, instead of focusing on one particular reconstruction method, we are interested in an optimal choice of the mesh size in a numerical reconstruction scheme assuming a maximal resolving power.  We shall provide an optimal selection of mesh size for the identifications of unknown inclusions from both the far-field and the near-field data in two and three dimensions.

This paper is organized as follows. In section \ref{sec:1}, we state the forward problems that we focus on.  We then provide in section \ref{sec:2} an  optimal mesh size selection given far-field data; whereas in section \ref{sec:3}, we consider the mesh size problem when near-field data are available.
In section \ref{sec:refine}, we describe how our optimal mesh size function may help to refine our mesh in order to resolve clustered inclusions.
In section \ref{sec:numerics}, we show that our choice of mesh size is indeed optimal by comparing several reconstructions from a typical numerical method with different mesh sizes.
Finally, some concluding remarks are presented in section \ref{con}. The results of this paper can help minimizing the computation complexity of methods designed for solving inverse scattering problems.  A similar analysis can be performed for electromagnetic and elastic inverse wave scattering problems.

\section{Problem description} \label{sec:1}
In this section, we shall describe the direct and inverse scattering problems in the presence of inhomogeneous inclusions (inside a sampling domain $\Omega$; see Figure \ref{fig:demo0}) embedded in a homogeneous background medium $\mathbb{R}^M (M \geq 2)$.   Consider an inhomogeneous inclusion represented by a (variable) contrast $q \in L^{\infty} ( \mathbb{R}^M ) $ such that $\text{supp}(q) \subset \Omega$, i.e., it vanishes outside $\Omega$.

\subsection{Forward problem}

We are now ready to state the considered forward problems.  For numerical simulations of solutions of these forward problems, we refer the reader, for instance, to \cite{LZ, LXZ3}.

\subsubsection{Plane wave illumination} \noindent Suppose that $u^i$ is an incident plane given by
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
u^i(x,d_y) := e^ {  i k  \, d_y \cdot x } \,,
\eqnx
where $d_y \in \mathbb{S}^{M}$ is the direction of incidence with $k$ being the wave number and $\mathbb{S}^{M}$ being the unit sphere in $\mathbb{R}^{M}$. Then, $u^i$ satisfies the Helmholtz equation
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
\Delta u^i+k^2 u^i = 0\quad \text{for}\quad x\in\mathbb{R}^M \,.
\eqnx
We now consider the total field $u$ as the solution to the following equation,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\Delta u+k^2\Big(1+q(x)\Big)u = 0\quad \text{for}\quad x\in\mathbb{R}^M \,,
\label{scattering1}
\eqn
subject to the outgoing Sommerfeld radiation condition:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
   \f{\partial u^s}{\partial r} - i k u^s = o\Big(\frac{1}{r^{(M-1)/2}}\Big) \quad \text{as}\quad r\to\infty,
    \label{sommerfield}
\eqn
which holds uniformly in all directions $x/|x|$, where $r=|x|$ and  $u^s := u - u^i$ is the scattered field.

From \eqref{scattering1}, due to the presence of inclusions in $\Omega$, the Lippmann-Schwinger representation formula of the total field $u$ yields
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
u^s(x, d_y)= u(x,d_y)- u^i(x,d_y) = k^2\int_\Omega q(z)G(x,z)u(z,d_y) \, dz,\quad d_y\in\mathbb{S}^{M}\,,
\label{compute1}
\eqn
where $G(x,z) , x,z \in \,\mathbb{R}^M \text{ and } x \neq z,$ is the Green function associated with the homogeneous background,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
G(x,z)= C_{M} k^{\frac{2-M}{2}} |x-z|^{\frac{2-M}{2}} H^{(1)}_{\frac{M-2}{2}}(k |x-z|) \,.
\eqn
Here,  $H_\nu^{(1)}$ is the Hankel function  of the  first kind of order $\nu$ and $C_M = -i (2 \pi)^{\frac{2-M}{2}} /4$. In particular, we have
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
G(x,z)= \left\{
\begin{aligned}
&\frac{i}{4}H^{(1)}_0(k|x-z|) &\text{for}\quad M=2,\\[2mm]
&\frac{e^{ik|x-z|}}{4\pi |x-z|} &\text{for}\quad M=3.
\end{aligned}
\right.
\eqn

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.55\textwidth]{./demo0.eps}\hfill{}
 \vskip -1truecm
 \caption{\label{fig:demo0}  \small{\emph{Plane wave illumination.}}}
 \end{figurehere}

\subsubsection{Point source illumination} \noindent On the other hand, suppose that $u^i$ is generated from a point source $ y \in \Gamma_s$, where $\Gamma_s$ is a surface. Then,
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
u^i(x,y) := G(x,y)= C_{M} k^{\frac{2-M}{2}} |x-z|^{\frac{2-M}{2}} H^{(1)}_{\frac{M-2}{2}}(|x-y|) \,,
\eqnx
which now satisfies
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
\Delta u^i+k^2 u^i = \delta_y  \quad \text{for}\quad x\in\mathbb{R}^M \,.
\eqnx
Consider the total field $u$ as the solution to the following equation:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\Delta u+k^2\Big(1+q(x)\Big)u = \delta_y  \quad \text{for}\quad x\in\mathbb{R}^M \,.
\label{scattering2}
\eqn
with $u^s := u - u^i$ satisfying the outgoing Sommerfeld radiation
condition \eqref{sommerfield}.

From \eqref{scattering2}, due to the presence of inclusions in $\Omega$, the total field $u$ satisfies the following equation:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
u(x,y)=u^i(x,y)+k^2\int_\Omega q(z)G(x,z)u(z,y) \, d z,\quad x\in\mathbb{R}^M\;\;\text{and}\;\; y\in\Gamma_s\,.
\label{compute2}
\eqn
Because the corresponding scattered field $u^s$ is measured by the receivers on the surface $\Gamma_r$, it has the following representation:
 \begin{eqnarray}}  \def\eqn{\end{eqnarray}
u^s(x,y)=k^2\int_\Omega q(z)G(x,z)u(z,y) \, d z,\quad x\in\mathbb{R}^M\;\;\text{and}\;\; y\in\Gamma_s\,.
\label{us2}
\eqn
We refer to Figure \ref{fig:demo} for a schematic configuration of the inverse medium scattering problem.
We may again consider measurements $u(x,y)$ along a surface $x \in \Gamma_r$ where $y \in \Gamma_s$.

\begin{figurehere}
 \vskip 0.5truecm
 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./demo.eps}\hfill{}
 \vskip -0.2truecm
 \caption{\label{fig:demo}  \small{\emph{Point source illumination.}}}
 \end{figurehere}

\subsection{Inverse problem} \label{example_haha}

The inverse problem of our interest is to recover the physical features of the inclusions inside $\Omega$ from knowledge of the scattered field $u^s$.

For {plane wave illuminations}, we consider the following problem:  given the measurements of $u^s(x,d_y)$ along a surface $x \in \Gamma_r$, where the distance $\text{dist} ( \Gamma_r, \Omega)  \gg 1$, and $d_y\in\mathbb{S}^{M}$, where $u$ satisfies \eqref{compute1}, we aim at recovering $q(x), x \in \Omega$.

In the case of {point source illuminations}, we consider the following problem:  given the measurements of $u^s(x,y)$ along a surface $x \in \Gamma_r$ and $y\in \Gamma_s $, where $u$ satisfies \eqref{compute2}, we aim at recovering $q(x), x \in \Omega$.

Note that there are many numerical methods addressing the number of required directions of incidence, see, for instance, \cite{BerBro, music, CK,  LXZ3, LXZ2}.
Since our work does not focus on any particular type of reconstruction methods, but aims at discussing how to choose the mesh size in general, we do not advocate any particular choice of numerical reconstruction methods. Nevertheless, in this subsection, for the sake of completeness, we state two examples of indicator functions in order to illustrate our main ideas.  In what follows, we provide a brief description of two types of indicator functions: index functions provided in the multilevel sampling method (MSM) \cite{L1, L2, LZ}, and the direct sampling method (DSM), see, for instance, \cite{Jin12, L1, LXZ3}.  Both methods work either with plane wave or point source illuminations; nonetheless, for the sake of simplicity, we only consider the case of plane wave illuminations. \\

\noindent \textbf{Example 1:}
We would like to briefly describe the index function shown in the MSM.
Given a choice of sampling points $x \in \Omega$.
Consider an incident direction $d_y \in \mathbb{S}^{M-1}$. We first introduce the contrast source function $\phi$, which is defined as follows,
\[\phi(x,d_y):= k^2 q(x) \, u(x,d_y) \quad \text{for}\;\;x\in \Omega.\]
The scattered field $u^s$ can thus be expressed into the following compact form
\[u^s(x_r,d_y)=\mathcal{P}_r\phi \,(x_r,d_y)  \quad \text{for}\;\;x_r\in \Gamma_r,\]
where the integral operator $\mathcal{P}_r$ is given, for $x_r \in \Gamma_r$, by
\[\mathcal{P}_r\varphi\,(x_r,d_y):=\int_\Omega G(x,x_r)\phi(x,d_y) dx.\]
The MSM aims at approximating the contrast source function at sampling points $x \in \Omega$ given by a mesh with multiple choices of $d_y$ via a backpropagation method; and then refine the choices of mesh sizes in order to obtain iteratively another set of $x \in \Omega$.
A backpropagation $\phi_b(\cdot, d_y)$, which is defined as follows, is employed to provide an (initial) approximation of $\phi(x,d_y)$:
\[
\phi_b : =\frac{\big\lVert \mathcal{P}_r^* u^s  \big\rVert_{L^2(\Omega)}^2}{\big\lVert  \mathcal{P}_r\,\mathcal{P}_r^* u^s  \big\rVert_{L^2(\Gamma_r)}^2}\mathcal{P}_r^*u^s.
\]
Here, $\mathcal{P}_r^*$ is known as the backpropagation operator which is represented as
\[\mathcal{P}_r^*\,\phi(x; d_{y}):=\int_{\Gamma_r}\overline{G(\xi;x)}\phi(\xi;d_y) \, d \sigma_\xi \quad \text{for}\;\; x\in \Omega.\]
In practice, the expression of contrast source function is numerically obtained by
\begin{equation}\label{eq:varphi}
 \tilde{\phi}_{b,n} (x) :=\frac{\big\lVert \mathcal{P}_r^* u^s(\cdot,  d_{y,n})  \big\rVert_{L^2(\Omega)}^2}{\big\lVert  \mathcal{P}_r\,\mathcal{P}_r^* u^s(\cdot,  d_{y,n})  \big\rVert_{L^2(\Gamma_r)}^2}\mathcal{P}_r^*u^s(x,  d_{y,n})  \quad \text{for}\;\; x\in \Omega,
 \end{equation}
for a finite number of incidences $d_{y,n} \in \mathbb{S}^{M-1}$, where $n=1,2\cdots,N$. With this in hand, we now define an index function or indicator function $\eta(x)$, which aims at numerically approximating $k^2 q(x)$ from the measured values of $\{ u^s(\cdot ,d_{y,n}) \}_{n=1}^N $, as the least-square solution to the system
\begin{equation}
 \tilde{\phi}_{b,n} = \eta(x) \, u(x,d_{y_n}),
 \end{equation}
for $n=1,2\cdots,N$. Here, the points $ x\in \Omega$ are the chosen sampling points in $\Omega$, and $\eta(x)$ can be explicitly given by
\begin{equation}\label{eq:eta}
\eta(x) :=\Re\left\{\frac{\sum\limits_{n=1}^N  \tilde{\phi}_{b,n} (x) \overline{\big(u^i(x, d_{y,n})+\mathcal{P}_\Omega  \tilde{\phi}_{b,n} \big)}}{\sum\limits_{n=1}^N\Big|u^i(x, d_{y,n})+\mathcal{P}_\Omega  \tilde{\phi}_{b,n}\Big|^2}\right\}\quad \text{for}\;\; x\in\Omega,
\end{equation}
where $\Re\{\cdot\}$ denotes the real part of a complex number, $ \tilde{\phi}_{b,n}$ are defined in \eqref{eq:varphi}, and the integral operator $\mathcal{P}_\Omega$ is defined, for all $x \in \Omega$, by
\[\mathcal{P}_\Omega\phi\,(x, d_y):=\int_\Omega G(\xi,x)\phi(\xi,d_y) d\xi.\]
Now the algorithm chooses a certain threshold of cutoff to numerically identify the support of the index function or indicator function $\eta$, and then the choices of $x \in \Omega$ are re-sampled, and the whole process is preformed again and iterated until a fine reconstruction profile is achieved. \\

\noindent \textbf{Example 2:}
Another example of indicator functions is provided by the DSM.  For a given incidence $d_y \in \mathbb{S}^M$, the index function is defined by
\begin{equation}\label{eqn:indexfcn}
I(x, d_y):=\frac{\Big| \langle u^s(\cdot, d_y),G(\cdot,x) \rangle_{L^2(\Gamma_r)}\Big|}{||u^s(\cdot, d_y)||_{L^2(\Gamma_r)}||G(\cdot,x)||_{L^2(\Gamma_r)}}, \,\quad \forall\,x\in \Omega\,,
\end{equation}
where the inner product $ \langle \cdot, \cdot \rangle_{L^2(\Gamma_r)}$ is given by
\[ \langle u^s(x_r, d_y),G(x_r,x) \rangle_{L^2(\Gamma_r)} : =\int_{\Gamma_r}u^s(x_r, d_y) \overline{G(x_r,x)}ds\,\]
and
\[||u^s(\cdot, d_y)||_{L^2(\Gamma_r)}=\bigg(\int_{\Gamma_r}\Big|u^s(x_r,d_y)\Big|^2 \, d \sigma_{x_r} \bigg)^{\frac{1}{2}},\quad ||G(\cdot,x)||_{L^2(\Gamma_r)}=\bigg(\int_{\Gamma_r}\Big|G(x_r,x)\Big|^2 \, d \sigma_{x_r} \bigg)^{\frac{1}{2}}.\]
In order to fully utilize measurements from multiple measurement events, i.e., with incidences $d_{y,n} \in \mathbb{S}^{M}, n = 1,...,N$, we define an index function $I(x)$ to combine different values of $I(x, d_{y,n})$ in the following way:
\[I(x)=\max_{n=1,..., N}\Big\{I(x, d_{y,n})\Big\}, \quad \forall x\in \Omega.\]
Here, $I_{m}$ means the index function for the $m$th data set.
 It has been discussed theorectically and illustrated numerically in \cite{LXZ3} that the function in \eqref{eqn:indexfcn} has
a relative large value when a sampling point $x$ is close to an inhomogeneous inclusion in $\Omega$ and decays rapidly as $x$ moves away from any of the inclusions. \\

We notice that in the two aforementioned examples, the choices of sampling points $x$ places a crucial row in constructing the indicator functions (either $\eta(x)$ or $I(x)$ in the previous two examples), and appropriate choices of the mesh size (or the density of sampling points) are of  importance for insuring the best possible quality of the reconstruction procedure.

\section{Optimal mesh size with far-field data under plane wave illumination} \label{sec:2}

In this section, we shall consider the optimal choice of the mesh size for the reconstruction of inclusions inside $\Omega$ from the far-field data.
We consider the case when $\Gamma_r = R \,  \mathbb{S}^{M-1}$ where $R \gg 1$. In this case, instead of having the scattered field $u^s$ as measurement, we may assume that we have obtained our measurement as the far-field pattern $u_\infty$ given by
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
u_\infty (d_x,d_y) = \lim_{|x| \rightarrow \infty} |x |^{\frac{M-1}{2}} \exp^{-i k |x|}  u_s( |x| d_x, d_y ),
\eqn
where $d_x$ and $d_y$ belong to the unit sphere $\mathbb{S}^{M-1}$.

Hence, we assume that the measurements $u_\infty (d_x,d_y)$, where $d_x, d_y \in \mathbb{S}^{M-1}$,  are performed.   For small wave number $k$, the Born approximation shows that the far-field pattern can be expressed as
\begin{eqnarray}}  \def\eqn{\end{eqnarray}\label{eq:far-field}
u_\infty (d_x,d_y) = \tilde{C}_{k,M} k^2 \int_\Omega q(z)e^{ik(d_x-d_y)\cdot z}dz + O (k^{4} |C_{k,M} |^2),
\eqn
where there error term is in $L^2$, and $$\tilde{C}_{k,M} = \frac{-i}{\sqrt{8 \pi}} \left( \frac{k}{2 \pi} \right)^{\frac{M-2}{2}} \exp \left( - \frac{(M-1)}{4} \pi i \right)$$ is a constant depending on the space dimension $M$ and  $k$.
We have the following result from the Nyquist-Shannon sampling theorem.
\begin{Theorem}
Let $k$ be fixed and small.
Denote by $\mathfrak{F}(q)(\xi)$ the Fourier transform of $q(y)$:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\mathfrak{F}(q)(\xi):=\int_{\Omega}q(z)e^{i \xi\cdot z}dz\quad \text{for} \quad \xi \in \mathbb{R}^M\,.
\eqn
Suppose that $\mathfrak{F}(q)(\xi) = 0 $ for all  $|\xi| > 2 k$. Then, the reconstruction of $q$ from $u_\infty (d_x,d_y)$, where $d_x, d_y \in \mathbb{S}^{M-1}$,  is with an error of order $O(k^2 |\tilde{C}_{k,M} |)$.  The reconstruction can be performed with a choice of the mesh size $h <  \f{\lambda}{2}$, where $\lambda :=2\pi/k$ is the operating wavelength.
\end{Theorem}

\begin{proof}
It is direct to observe that
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\mathfrak{F}(q)(\xi) =  \frac{1}{\tilde{C}_{k,M} k^2}  u_\infty (d_x,d_y) + \mathcal{O} (k^{2} |\tilde{C}_{k,M}|) ,
\eqn
where $\xi:= k ( d_x - d_y)  \in\{\xi:|\xi|\leq 2 k\}$. Hence, from $u_\infty (d_x,d_y), \, d_x, d_y \in \mathbb{S}^{M-1}$, we obtain the values of $\mathfrak{F}(q)(\xi)$ for $|\xi|\leq 2 k$ up to an error of the order $\mathcal{O} (k^{2} |\tilde{C}_{k,M}|)$ in $L^2$. Since $\mathfrak{F}(q)(\xi) = 0 $ for all  $|\xi| > 2 k$, we have
$\mathfrak{F}(q)(\xi)$  up to an error of the order $\mathcal{O} (k^{2} |\tilde{C}_{k,M}|)$ in $L^2$.
Applying the inverse Fourier transform, we obtain $q$ up to an error of the order $\mathcal{O} (k^{2} |\tilde{C}_{k,M}|)$ in $L^2$.  An application of the Nyquist-Shannon sampling theorem yields the last statement in the theorem with a choice of $h < \f{\lambda}{2}$.
\end{proof}
Notice that the above result shows that we may maximize our choice of the mesh size by taking $h = \f{\lambda}{2} - \epsilon$ for a small $\epsilon$ in order to minimize the computation complexity.

In order to recover the contrast function $q(x)$, we should discretize \eqref{eq:far-field} into the following form:
\[u_\infty\approx \gamma k^2\sum_{j=1}^N s_j \,q(z_j)e^{ik(d_x-d_y)\cdot z_j}.\]
Here,  the optimal $N$ is $\big\lceil 4|\Omega|/\lambda^2 \big\rceil$ when $M = 2$, where $\lceil x\rceil$ means the largest integer that is less than or equal to $x$, $s_j$ denotes the $j$-th grid point in $\Omega$, and $|\Omega|$ represents the area (volume) of $\Omega$.

\section{Optimal mesh size with near-field data under point source illumination}  \label{sec:3}
In this section, optimal mesh size for contrast reconstruction from near-field data will be discussed.
Based on the first term of the Neumann series solution to the Lippmann-Schwinger integral equation which is known as the Born approximation, the scattered field $u^s$ is approximated by
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
u^s_B(x,y)=k^2\int_{\Omega} q(z)G(x,z)u^i(y,z) \, d z + O(k^4),\quad x\in\Gamma_r\;\;\text{and}\;\; y\in\Gamma_s\,.
\eqn
Furthermore, the incident field $u^i$ in this case is provided by the Green function. Thus, the estimated scattered field $u^s_B$ can be written into the following multi-static response (MSR):
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
[\text{MSR}(q)](x,y)  := \frac{1}{k^2} u^s_B(x,y)  =\int_{\Omega} q(z)G(x,z)G(z,y) \, d z + O(k^2) ,\quad x\in\Gamma_r\;\;\text{and}\;\; y\in\Gamma_s\,.
\eqn

In practice, we again discretize the sampling domain $D$ into $N$ sub-regions and denote the center of each sub-region  by $z_j$ for $j=1,2,\cdots, N$. Based on this discretization, we have the following MSR approximation operator $A$ up to a normalization constant $\hat{C}_{k,M}$ depending on $k$ and $M$, which is the dominant part of the near-field data,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}{\label{A}}
A :=
\sum \limits_{j=1}^N q(z_j)  |z_j-x|^{\frac{2-M}{2}}  |z_j-y|^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z_j-x|) H^{(1)}_{\frac{M-2}{2}}(k|z_j -y|)   \quad \text{for} \quad x\in\Gamma_r\;\;\text{and}\;\; y\in\Gamma_s\,,
\eqn
where $z_j \in \Omega$.

Let us consider for simplicity the case when $\Gamma_r = \Gamma_s.$
In order to investigate the optimal mesh size for the reconstruction process, we need to determine the number of observable singular values for the MSR operator $A$ when $z_j$ are close to each other.

For this purpose, we first define a conjugation map on $L^2(\Gamma)$ as follows,
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
C: L^2(\Gamma) \rightarrow L^2(\Gamma) \quad f \mapsto \bar{f} \,.
\eqnx
Now, from the fact that the operator $A$ is compact and $C$-complex symmetric, namely, $A^* = C A C$. We have the following Danciger's variational principle for the singular values $\{ \sigma_i \}$ of $A$ \cite{Garcia},
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\sigma_n = \min \limits_{\text{codim}_{\mathbb{R}} V = n} \max \limits_{\substack{p \in V \\ || p || = 1}} {\Re}[A p, p] \,,
\label{varvar}
\eqn
where $[Ap, p]=p^TAp$ and $\sigma_0 \geq \sigma_1 \geq ... \geq 0$.

Before we state the following main result, we denote the locations of receivers and sources by $\Gamma:=\Gamma_r\cup\Gamma_s$.
Firstly, we consider the simplest case when $N=1$. Because $q$ is a constant, we investigate the following operator $A_1$ which is the key part of the operator $A$:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
A_1 := a  |z-x|^{\frac{2-M}{2}}   |z-y|^{\frac{2-M}{2}}  H^{(0)}_{\frac{M-2}{2}}(k|z -x|) H^{(1)}_{\frac{M-2}{2}}(k|z -y|) = w \otimes w^T\quad \text{for}\quad x,y\in\Gamma,
\eqn
where $a= || |z-x|^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z - x |)  ||^2_{L^2(\Gamma)}$, $\otimes$ means the Kronecker product and $$w :=  |z-x|^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z -x|) \bigg/ \Big\| |z-x|^{\frac{2-M}{2}}   H^{(1)}_{\frac{M-2}{2}}(k|z - x |)  \Big\|_{L^2(\Gamma)}$$ is a unit vector.
It is easy to see that
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
{\Re}[A_1 p, p] = {\Re} (p^T A_1 p ) = {\Re} ( w^T p )^2
\begin{cases}
= 0 & \text{ if }\; p \perp_{\mathbb{R}} \overline{w} \,, \\[2mm]
\neq 0 & \text{ otherwise},
\end{cases}
\eqn
and that ${\Re}[A_1 p, p] \leq 1$. We infer from Danciger's principle \cite{Dan} that the singular values of $A_1$ are as follows: $\sigma_0 = 1$ and $\sigma_i=0\, (i\geq1)$.

Secondly, we investigate the case when $N=2$.
With the help of  \eqref{A}, the main part of the MSR matrix is of the form
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
A_2 := b_1w_1 \otimes w_1^T + b_2w_2 \otimes w_2^T , \label{A2}
\eqn
where $b_j= ||   |z_j-x|^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z_j - x |)  ||^2_{L^2(\Gamma)}$, and $$w_j:= |z_j-x|^{\frac{2-M}{2}}    H^{(1)}_{\frac{M-2}{2}}(k|z_j -x|) \bigg/ \Big\|  |z_j-x|^{\frac{2-M}{2}}    H^{(1)}_{\frac{M-2}{2}}(k|z_j - x |)  \Big\|_{L^2(\Gamma)},$$ for $j = 1,2$, are unit vectors.
Similarly, we are able to directly obtain that
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
{\Re}[A_2 p, p] = {\Re} (p^T A_2 p ) \begin{cases}
= 0 & \text{ if } p \perp_{\mathbb{R}} \text{span}_{\mathbb{R}}\{ \overline{w_1}, \overline{w_2} \}\,, \\[2mm]
\neq 0 & \text{  otherwise}\,.
\end{cases}
\label{var2}
\eqn
Now, let us define $\{ \overline{h_1} , \overline{h_2} \} := \Big\{ \frac{\overline{w_1  + w_2}}{\|w_1  + w_2\|_{L^2(\Gamma)}} , \frac{\overline{w_1 - w_2}}{\|w_1 - w_2\|_{L^2(\Gamma)}} \Big\}$. Then, we have
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\| h_1 \|_{L^2(\Gamma)} = \| h_2 \|_{L^2(\Gamma)} = 1 \,, \quad {\Re} \langle h_1, h_2 \rangle_\Gamma  = 0 \,,
\label{var22}
\eqn
and that $\text{span}_{\mathbb{R}}\{ \overline{h_1}, \overline{h_2} \} = \text{span}_{\mathbb{R}}\{ \overline{w_1}, \overline{w_2} \}$.  Hence, the operator $A_2$ can be written as,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
A_2 =\frac{1}{2}\|w_1  + w_2\|_{L^2(\Gamma)}^2 \,  h_1 \otimes h_1^T + \frac{1}{2} \|w_1  - w_2\|_{L^2(\Gamma)}^2 \, h_2 \otimes h_2^T.
\label{var222}
\eqn
Considering \eqref{var2}-\eqref{var222}, Danciger's principle yields that the singular values of $A_2$ satisfy
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\sigma_0 &=& \frac{1}{2}\max \Big\{ \|w_1  + w_2\|_{L^2(\Gamma)}^2 , \|w_1  - w_2\|_{L^2(\Gamma)}^2 \Big\} = 1 + |\Re\langle w_1,w_2\rangle|\, , \\[2mm]
\sigma_1 &=& \frac{1}{2}\min \Big\{ \|w_1  + w_2\|_{L^2(\Gamma)}^2 , \|w_1  - w_2\|_{L^2(\Gamma)}^2 \Big\} = 1 - |\Re\langle w_1,w_2\rangle|\, , \\[2mm]
\sigma_l &=& 0,  \quad \quad  \forall\, l >1,
\eqn
where $\langle w_1,w_2\rangle=\int_\Gamma \overline{w_1}w_2d\sigma$. The above results can be summarized in the following theorem.
\begin{Theorem}
Suppose that the operator $A$ can be represented as follows,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
A := w_1 \otimes w_1^T + w_2 \otimes w_2^T \,,
\eqn
where $w_j$ are unit vectors in $L^2(\Gamma)$ for $j= 1,2$. We have the following expression for the singular values of $A$:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\sigma_0 = 1 + |{\Re}\langle w_1,w_2\rangle|\, ,
\sigma_1 = 1 - |{\Re}\langle w_1,w_2\rangle|\, ,
\sigma_2 = \sigma_3 =\cdots= 0,
\eqn
and the quotient $\alpha$ between the two non-zero singular values is expressed as,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\alpha := \frac{\sigma_1}{\sigma_0} =  \frac{1 - |{\Re}\langle w_1,w_2\rangle|}{1 + |{\Re}\langle w_1,w_2\rangle|}\,.
\eqn
\end{Theorem}

We have noticed that the ratio $\alpha$ represents the relative observability of the second singular value with respect to the first one.
Therefore, it quantifies the ability to distinguish between two inhomogeneities from the measurements in the form of the MSR matrix.  Consequently, we introduce the following definition of the $\alpha$-distinguishable with respect to $\Gamma$, for $0 < \alpha < 1$.
\begin{Definition}
Two point scatterers at locations $z_1$ and $z_2$ with equal weights
are $\alpha$-distinguishable with respect to $\Gamma_r$ if the following inequality is satisfied:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
1 - \alpha \geq  \frac{1 - |{\Re}\langle {w}_1,{w}_2\rangle|}{1 + |{\Re}\langle {w}_1, {w}_2\rangle|},
\label{alpha1}
\eqn
where  $$w_j:= \left( k |z_j-x| \right)^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z_j -x|) \bigg/ \Big\|  \left( k |z_j-x| \right)^{\frac{2-M}{2}}  H^{(1)}_{\frac{M-2}{2}}(k|z_j - x |)  \Big\|_{L^2(\Gamma)}$$  with $j=1,2$ and $\forall{z}_j\in D$, $x\in\Gamma_r$.
\end{Definition}

Without loss of generality, we consider two point scatterers located at $z_1 = z$ and $z_2 = z + \varepsilon v$, where $v$ is a unit vector and $\varepsilon$ is a distance parameter which can be viewed as the mesh size. We acquire the interval of $\varepsilon$ when the two point scatterers are $\alpha$-distinguishable with respect to $\Gamma$.
From the series expansion of $(1 + a)^{- \beta}$ for all $\beta \neq 0 $ w.r.t. $a$, we derive that if $ \varepsilon < |z-x| $, then
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
|z + \varepsilon v - x |^{- \beta} = \sum_{r=0}^\infty \sum_{s=0}^r \frac{\Gamma(\frac{\beta}{2} + r)}{\Gamma(\frac{\beta}{2}) s! (r-s) ! }  (-1)^r  \varepsilon^{r + s} \cos^{r-s} (\theta_{z-x,v} )   |z - x |^{ - \beta - r - s},
\eqnx
where $\theta_{z-x,v}$ denotes the relative angle between $z-x$ and $v$ and that the series converges absolutely and uniformly when $ \varepsilon < (1-\gamma) \min \{ 1, |z-x| \} $,  for some $0 < \gamma < 1$.
Together with Graf's fromula \cite{Watson}, we have if $ \varepsilon <  (1-\gamma) \, \min\{ 1, \text{dist}(\Gamma, z) \} \, $, then the following expansion holds:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
& & {\Re} \left(\int_{\Gamma}  |z-x|^{\frac{2-M}{2}}    |z + \varepsilon v -x|^{\frac{2-M}{2}}  H^{(2)}_{\frac{M-2}{2}}(k|z -x|) H^{(1)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) \, d \sigma_x \right) \notag  \\[2mm]
&= &
\begin{cases}
\sum \limits_{n\in \mathbb{Z}} p_{n} J_n(k \varepsilon )  & \text{ if } M = 2, \\[3mm]
\sum \limits_{n\in \mathbb{Z}}  \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r   \, p_{nrs} \,  J_n(k \varepsilon ) \, \varepsilon^{r+s} & \text{ if } M > 2,
\end{cases}
\label{seriesA}
\eqn
where $p_n$ and $p_{n,r,s}$ are explicitly given by
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
p_n &:=& {\Re} \left(\int_{\Gamma}  H^{(2)}_{0}(k|z -x|) H^{(1)}_{n}(k|z -x|) e^{- i n \theta_{z-x,v}} \, d \sigma_x \right)\,
\label{pn} \\[2mm]
p_{nrs} &:= & (-1)^r   \frac{\Gamma(\frac{M-2}{4} + r)}{\Gamma(\frac{M-2}{4}) s! (r-s) ! } \,  {\Re} \bigg(\int_{\Gamma}  |z-x|^{2-M - r - s}   H^{(2)}_{\frac{M-2}{2}}(k|z -x|) \notag \\[2mm]
      &&\cdot H^{(1)}_{\frac{M-2}{2} + n}(k|z -x|) \cos^{r-s} (\theta_{z-x,v} )  e^{- i n \theta_{z-x,v}} \, d \sigma_x \bigg)\,\label{pnrs}
\eqn
}with $\theta_{z-x,v}$ denoting the relative angle between $z-x$ and $v$. Likewise, we have
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
& &\int_{\Gamma}  |z + \varepsilon v -x|^{2-M}  H^{(2)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) H^{(1)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) \, d \sigma_x   \notag\\[2mm]
&= &
\begin{cases}
\sum \limits_{n,m\in \mathbb{Z}} h_{nm} J_n(k \varepsilon )  J_m(k \varepsilon ) & \text{ if } M = 2 ,\\[3mm]
\sum \limits_{n,m\in \mathbb{Z}} \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r  h_{nmrs}  \, J_n(k \varepsilon )  J_m(k \varepsilon ) \, \varepsilon^{r+s} & \text{ if } M > 2,
\end{cases}
\label{seriesB}
\eqn
where $h_{nm}$ and $h_{nmrs}$ are given by
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
h_{nm} &:=& {\Re} \left(\int_{\Gamma} H^{(2)}_{m}(k|z -x|) H^{(1)}_{n}(k|z -x|) e^{- i (n-m)\theta_{z-x,v}} \, d \sigma_x \right)\,,
\label{hnm} \\[2mm]
h_{nmrs}  &:=&(-1)^r
\frac{\Gamma(\frac{M-2}{2} + r)}{\Gamma(\frac{M-2}{2}) s! (r-s) ! }
 {\Re} \bigg(\int_{\Gamma} |z - x |^{ 2 - M - r - s} H^{(2)}_{\frac{M-2}{2}+m}(k|z -x|)\notag\\[2mm]
 &&\cdot H^{(1)}_{\frac{M-2}{2} + n}(k|z -x|)  \cos^{r-s} (\theta_{z-x,v} )  e^{- i (n-m)\theta_{z-x,v}} \, d \sigma_x \bigg)\,.\label{hnmrs}
\eqn
}Assume that $0 <R_0 \leq \text{dist}(z,\Gamma) \leq R_1$ for some $R_0, R_1 \in \mathbb{R}$, and consider the asymptotic property of $J_m$ \cite{handbook}. We can infer that both the series \eqref{seriesA} and \eqref{seriesB} converge absolutely.
If we further impose the restriction that $\frac{k \, \varepsilon}{2} \leq 1$, then from the following well-known Taylor series of the Bessel functions \cite{Watson},
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
J_n(k \varepsilon ) = (\text{sgn}(n))^{n} \sum_{l \in \mathbb{N}} \frac{(-1)^{l}}{l! (l+|n|)!} \left(\frac{k \,  \varepsilon}{2}\right)^{2 l + |n|} \,,
\eqn
we immediately obtain the following expansion of the inner product:
{\small
\begin{eqnarray}
& & {\Re} \left(\int_{\Gamma}  |z-x|^{\frac{2-M}{2}}    |z + \varepsilon v -x|^{\frac{2-M}{2}}  H^{(2)}_{\frac{M-2}{2}}(k|z -x|) H^{(1)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) \, d \sigma_x \right)
\notag \\[2mm]
&=&
\begin{cases}
\sum \limits_{l \in \mathbb{N}} \sum \limits_{n\in \mathbb{Z}} (\text{sgn}(n))^{n} p_n \frac{(-1)^{l}}{l! (l+|n|)!} \left(\frac{k \,  \varepsilon}{2}\right)^{2 l + |n|}  \;\text{ for } M = 2 ,\\[3mm]
\,
\sum \limits_{l \in \mathbb{N}}  \sum \limits_{n\in \mathbb{Z}}  \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r   \, (- \text{sgn}(n))^{n} p_{nrs}  \frac{(-1)^{l}}{l! (l+|n|)!} \left(\frac{k \,  \varepsilon}{2}\right)^{2 l + |n|} \varepsilon^{r+s} \;  \text{ for } M > 2,
\end{cases}
\label{series1}
\\[2mm]
& &\int_{\Gamma}  |z + \varepsilon v -x|^{2-M}  H^{(2)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) H^{(1)}_{\frac{M-2}{2}}(k|z + \varepsilon v -x|) \, d \sigma_x  \notag \\[2mm]
&=&
\begin{cases}
 \sum \limits_{l, p \in \mathbb{N}} \sum \limits_{n,m \in \mathbb{Z}}  (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m} h_{nm} \frac{(-1)^{l+p}}{l!p! (l+|n|)! (p+|m|)!} \left(\frac{k \,  \varepsilon}{2}\right)^{2 (l+p) + |n| + |m| } \; \text{ for } M = 2 , \\[3mm]
 \sum \limits_{l, p \in \mathbb{N}} \sum \limits_{n,m \in \mathbb{Z}} \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r   (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m} h_{nmrs} \frac{(-1)^{l+p}}{l!p! (l+|n|)! (p+|m|)!} \left(\frac{k \,  \varepsilon}{2}\right)^{2 (l+p) + |n| + |m| } \varepsilon^{r+s} \; \text{ for } M > 2,
\end{cases}
\label{series2}
\end{eqnarray}
}where the two series at the right-hand side converge absolutely for $\frac{k \, \varepsilon}{2} \leq 1$ and  $ \varepsilon <   (1-\gamma)\, \min \{ 1, \text{dist}(\Gamma, z) \}$.

\subsection{$\alpha$-distinguishability when $M=2$}

For simplicity, let us first focus on the case when $M = 2$. We obtain the following estimate,
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\left| {\Re} \left(\int_{\Gamma} H^{(2)}_{0}(k|z -x|) H^{(1)}_{0}(k|z + \varepsilon v -x|) \, d \sigma_x \right)- p_0 + \frac{k \,  \varepsilon}{2} ( p_1 - p_{-1})  \right| &\leq C_1 \frac{k ^2\,  \varepsilon^2}{4}, \quad \quad  \label{approx1}\\[2mm]
\left| \int_{\Gamma} H^{(2)}_{0}(k|z + \varepsilon v -x|) H^{(1)}_{0}(k|z + \varepsilon v -x|) \, d \sigma_x -
h_{00} - \frac{k \,  \varepsilon}{2} \sum_{|m|+|n|=1} (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m}  h_{nm}  \right| &\leq C_2 \frac{k^2 \,  \varepsilon^2}{4},\quad \quad  \label{approx2}
\eqn}
where $p_0 = h_{00} =  || H^{(1)}_{0}(k|z - x |)  ||^2_{L^2(\Gamma)}$ and $C_1$ and $C_2$ are two finite constants having the following forms:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C_1 &=& \sum_{l>0} \sum \limits_{n\in \mathbb{Z}}  \frac{|p_n|}{l! (l+n)!} + \sum \limits_{n>1}  \frac{|p_n|}{ n !}
\label{constant1} \,,\\[2mm]
C_2 &=& \sum_{l+p >0} \sum \limits_{m,n\in \mathbb{Z}}  \frac{|h_{nm}|}{l!p! (l+n)! (p+m)!} +\sum \limits_{m+n>1} \frac{|h_{nm}|}{ n! m!} \,
\label{constant2}
\eqn
if  $\frac{k \, \varepsilon}{2} \leq 1$ and  $ \varepsilon < (1-\gamma)\, \min \{ 1, \text{dist}(\Gamma, z) \} $.
From \eqref{approx1} and \eqref{approx2}, we have the following approximation:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
1 - \left| {\Re}\langle w_1,w_2\rangle \right|^2 = 1 - \frac{\left(p_0 - ( p_1 - p_{-1}) \frac{k \,  \varepsilon}{2} + \mathcal{O}\left(\frac{k^2 \,  \varepsilon^2}{4} \right) \right)^2}{ p_0 \left(p_0 + \frac{k \,  \varepsilon}{2}\sum\limits_{|m|+|n|=1}  (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m}   h_{nm}  + \mathcal{O}\left(\frac{k^2 \,  \varepsilon^2}{4}\right)  \right) } \,.
\label{taylortaylor}
\eqn
We further assume that
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
    \frac{k \,  \varepsilon}{2} \leq  \sqrt{ \frac{p_0^2 k_1^2 }{4 C_2^2} + \frac{p_0 (1 - \delta)}{C_2} } -  \frac{p_0 |k_1| }{2 C_2}
    \label{bounddelta}
\eqn
for some $ 0 < \delta < 1$, where
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
k_1= \frac{\sum_{|m|+|n|=1} (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m}  h_{nm} }{p_0} \,.
\label{k1}
\eqn
Then, we have
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\max \left\{ \frac{k \,  \varepsilon}{2}   \left| k_1 - \frac{C_2}{p_0} \frac{k \,  \varepsilon}{2}  \right|,  \frac{k \,  \varepsilon}{2}   \left| k_1  +  \frac{C_2}{p_0} \frac{k \,  \varepsilon}{2}  \right| \right\} \leq 1 - \delta,
\eqn
and thus, we can derive from the Taylor expansion and \eqref{taylortaylor} that
\begin{eqnarray}}  \def\eqn{\end{eqnarray} \label{estima}
\left| 1 - |{\Re}\langle w_1,w_2\rangle|^2 - \left( \frac{ 2 (p_1 - p_{-1})}{p_0} + k_1 \right)\frac{k \,  \varepsilon}{2}
\right| \leq C \frac{k^2 \,  \varepsilon^2}{4}
\,
\eqn
with the constant $C$ given by
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C = \frac{C_2}{p_0}  \left(1 + 2 \left|k_1 \right| \delta^{-1} \right) +  k_1^2 \delta^{-1} +  \frac{C_2^2}{p_0^2} \delta^{-1}
+ \frac{ 2 ( | p_1 | + |p_{-1}| ) | k_1| }{p_0}
+ \frac{ ( | p_1 | + |p_{-1}|  )^2 + 2 (p_0 +  | p_1 | + |p_{-1}|  )C_1 + C_1^2}{p_0^2}, \quad
\label{constant3}
\eqn
}when $0<\delta<1$.
We can rewrite the condition in \eqref{alpha1} as follows,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\frac{4  (1- \alpha) }{(2 - \alpha)^2}  \geq 1 - |{\Re}\langle {w}_1,{w}_2\rangle|^2.
\eqn
With the help of \eqref{estima}, the two point inclusions with the same magnitude
are $\alpha$-distinguishable if the following inequality holds:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\frac{4 (1- \alpha)}{(2 + \alpha)^2}  \geq  \left( \frac{ 2 ( p_1 - p_{-1})}{p_0} + k_1 \right)\frac{k \,  \varepsilon}{2} + C \left( \frac{k \,  \varepsilon}{2}\right)^{2} \,.
\eqn
This directly infers the following theorem.
Note that we may always assume that $k \varepsilon < 1/2$.

\begin{Theorem} \label{important_theorem}
When the dimension $M = 2$, consider a point $z$ such that $0 <R_0 \leq \text{dist}(z,\Gamma) \leq R_1$ for some $R_0, R_1 \in \mathbb{R}$, where $\Gamma = \Gamma_s = \Gamma_r$, $ \frac{k \,  \varepsilon}{2} \leq  \min \left\{ \frac{1}{4} , \sqrt{ \frac{p_0^2 k_1^2 }{4 C_2^2} + \frac{p_0 (1 - \delta)}{C_2} } -  \frac{p_0 |k_1| }{2 C_2} \right\}$ for some $0< \delta <1 $,  $ \varepsilon < (1-\gamma) \, \min \left\{ 1,  \text{dist}(\Gamma, z) \right \} $ for some $ 0< \gamma <1 $, $p_n$ are defined as in \eqref{pn}, $h_{nm}$ as in \eqref{hnm}, and let $C_1,C_2,C, k_1$ be stated in \eqref{constant1}, \eqref{constant2} \eqref{constant3} and \eqref{k1}, respectively.
The two point scatterers with equal weights located at $z$ and $z+ \varepsilon v$ are $\alpha$-distinguishable with respect to $\Gamma$ if the following inequality holds:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
 \frac{k \,  \varepsilon}{2}
\leq
- \left( \frac{ 2 p_1 - 2 p_{-1} + k_1 p_0 }{ 2 C p_0}  \right) + \sqrt{ \left( \frac{ 2 p_1 - 2 p_{-1} + k_1 p_0}{ 2 C p_0} \right)^2 +  \frac{16 (1- \alpha)}{ C (2 + \alpha)^2}  } \,.
\label{conclusion}
\eqn
\end{Theorem}

Unfortunately, both the conditions and conclusions of the above theorem are too complicated to be practical for choosing an optimal mesh size. In order to obtain an optimal mesh size for practical use, we provide  approximations for the coefficients in the above inequalities: we estimate $C_1$ and $C_2$  by
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C_1 \approx p_0
\, , \quad
C_2 \approx 2 p_0 \, .
\label{constantapprox}
\eqn
By considering \eqref{series1} and  \eqref{series2}, and further imposing that $\delta = 1/2$, the following approximation of the constant $C$ holds:
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C \approx  13  + 8 \left|k_1 \right|  +  2 k_1^2
+ \frac{ ( | p_1 | + |p_{-1}|  )^2 + 2 ( | p_1 | + |p_{-1}|  ) (1 + |k_1| )p_0 }{p_0^2}.
\label{constantapprox2}
\eqn}
The condition \eqref{bounddelta} with $k \varepsilon /2 < 1/4$ can also be simplified to be the above approximate condition
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
    \frac{k \,  \varepsilon}{2} <   \min \left\{ \frac{1}{4}, \sqrt{ \frac{k_1^2 }{16} + \frac{1}{4} } -  \frac{|k_1| }{4} \right\}.
\eqn
Combining the above estimations of the constants $C_1,C_2,$ and $C$ with \eqref{conclusion}, we obtain the following approximate description of an \textbf{anisotropic} optimal density  at any point $z$ and any direction $v$ for contrast reconstruction from the near-field data in terms of the $\alpha$-distinguishably,
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
h_{z,v}(\alpha) :=  \frac{2}{k}\min \left\{ \frac{1}{4} , \sqrt{ \frac{k_1^2 }{16} + \frac{1}{4} } -  \frac{|k_1| }{4}  ,
    \left( \frac{ 2 p_1 - 2 p_{-1} + k_1 p_0 }{ 2 C p_0}  \right) - \sqrt{ \left( \frac{ 2 p_1 - 2 p_{-1} + k_1 p_0}{ 2 C p_0} \right)^2 + \frac{16 (1- \alpha)}{ C (2 + \alpha)^2}  }
    \right\}.
\label{finalfunction}
\eqn
Note that  $p_0, p_1, p_{-1}, k_1,$ and $C$ depend on $z$, $v,$ and $\Gamma$.

\begin{Remark}
A few remarks are in order:
\begin{enumerate}
\item[(i)]
In practice, the measurements are taken over a finite number of points on $\Gamma = \Gamma_r = \Gamma_s$ and hence, the integrals in \eqref{pn},  \eqref{pnrs}, \eqref{hnm}, and \eqref{hnmrs} can be evaluated numerically as a sum over the points.
\item[(ii)]
In practice, small regions with fine details are usually anticipated to be recovered with high resolution. Thus, a uniform discretization in the whole sampling domain may not be the best.
Hence, we may apply adaptivity, and use the optimal mesh size in \eqref{finalfunction} in the region with the requirements of high resolution and a coarse mesh to the other regions.
\item[(iii)]
To get the interpretation of the optimal mesh size correct, actually it means if you refine further than the mesh size given, it will not improve the resolution.  Hence, if the object is far away, then it is almost identical to using far-field reconstruction, and we do not need to refine better than $\lambda/2$.
\item[(iv)]
A similar analysis can be performed in the general case when $\Gamma_r$ and $\Gamma_s$ may not coincide, but we skip the analysis for the sake of simplicity.
\end{enumerate}
\end{Remark}

\subsection{$\alpha$-distinguishability when $M > 2$}

Similarly, when $M > 2$, we have,
{\footnotesize
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\left| {\Re} \left(\int_{\Gamma} H^{(2)}_{0}(k|z -x|) H^{(1)}_{0}(k|z + \varepsilon v -x|) \, d \sigma_x \right)- p_{000} + \frac{k \,  \varepsilon}{2}   \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  \,   \,  \right| &\leq C_{1,M} \frac{k ^2\,  \varepsilon^2}{4}, \quad \quad   \label{approx1M}\\[3mm]
\left| \int_{\Gamma} H^{(2)}_{0}(k|z + \varepsilon v -x|) H^{(1)}_{0}(k|z + \varepsilon v -x|) \, d \sigma_x -
h_{0000} + \frac{k \,  \varepsilon}{2}  \left(\sum_{|m|+|n|=1} (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m}  h_{nm00} +  \frac{2 h_{0010} }{k} \right)  \right| &\leq C_{2,M} \frac{k^2 \,  \varepsilon^2}{4}, \quad \quad  \label{approx2M}
\eqn
}where $p_{000} = h_{0000} =  \Big\| | z-x |^{\frac{M-2}{2}} H^{(1)}_{\frac{M-2}{2}}(k|z - x |)  \Big\|^2_{L^2(\Gamma)}$ and $C_{1,M}$ and  $C_{2,M}$ are two constants given if  $\frac{k \, \varepsilon}{2} \leq 1$ and  $ \varepsilon < (1-\gamma)\, \min \{ 1, \text{dist}(\Gamma, z) \} $ by
{\footnotesize
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C_{1,M} &=& \sum \limits_{l>0}   \sum \limits_{n\in \mathbb{Z}}  \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r    \frac{|p_{nrs}|}{l! (l+n)!} + \sum \limits_{n>1}   \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r    \frac{|p_{nrs}|}{ n !}
\label{constant1M} \,,\\[3mm]
C_{2,M} &=& \sum_{l+p >0} \sum \limits_{m,n\in \mathbb{Z}}  \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r     \frac{|h_{nmrs}|}{l!p! (l+n)! (p+m)!} +\sum \limits_{m+n>1}  \sum \limits_{r \in \mathbb{N}} \sum \limits_{s=0}^r    \frac{|h_{nmrs}|}{ n! m!} \,.
\label{constant2M}
\eqn
}Hence, a similar analysis can be performed as in the previous subsection provided that $ \varepsilon < (1-\gamma)\, \min \{ 1, \text{dist}(\Gamma, z) \} $ and
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
    \frac{k \,  \varepsilon}{2} \leq  \min\left\{1, \sqrt{ \frac{p_{000}^2 k_1^2 }{4 C_{2,M}^2} + \frac{p_{000} (1 - \delta)}{C_{2,M}} } -  \frac{p_{000} |k_{1,M}| }{2 C_{2,M}} \right\}
    \label{bounddeltaM}
\eqn
for some $ 0 < \delta < 1$, where
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
k_{1,M}=  \left(\sum_{|m|+|n|=1} (\text{sgn}(n))^{n}  (\text{sgn}(m))^{m}  h_{nm00} +  \frac{2 h_{0010} }{k} \right) \bigg/p_{000} \,.
\label{k1M}
\eqn
Then, by a Taylor expansion, we readily see that
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
\left| 1 - |{\Re}\langle w_1,w_2\rangle|^2 - \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right)\frac{k \,  \varepsilon}{2}
\right| \leq C \frac{k^2 \,  \varepsilon^2}{4}
\,
\eqnx
with the constant $C_M$ being given by
{\footnotesize
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
C_M &=& \frac{C_{2,M}}{p_{000}}  \left(1 + 2 \left|k_{1,M} \right| \delta^{-1} \right) +  k_{1,M}^2 \delta^{-1} +  \frac{C_{2,M}^2}{p_{000}^2} \delta^{-1}
+ \frac{ 2  \left(|p_{1\,0\,0}| + |p_{-1\,0\,0}|  - \frac{2 |p_{0\,1\,0}|}{k}  \right)  ) | k_{1,M}| }{p_{000}} \notag \\[2mm]
& &
+ \frac{  \left(|p_{1\,0\,0}| + |p_{-1\,0\,0}|  - \frac{2 |p_{0\,1\,0}|}{k}  \right) ^2 + 2 \left(p_{000} +  \left(|p_{1\,0\,0}| + |p_{-1\,0\,0}|  - \frac{2 |p_{0\,1\,0}|}{k}  \right)     \right) C_{1,M} + C_{1,M}^2}{p_{000}^2}, \quad
\label{constant3M}
\eqn
}when $0<\delta<1$.
Hence, condition \eqref{alpha1} is fulfilled if
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
\frac{4 (1- \alpha)}{(2 + \alpha)^2}  \geq \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right) \frac{k \,  \varepsilon}{2} + C_M \left( \frac{k \,  \varepsilon}{2}\right)^{2} \,.
\eqn
Likewise, adding the constraint $k \varepsilon < 1/4$, this directly infers the following theorem.
\begin{Theorem}
For $M > 2$, consider a point $z$ such that $0 <R_0 \leq \text{dist}(z,\Gamma) \leq R_1$ for some $R_0, R_1 \in \mathbb{R}$, where $\Gamma = \Gamma_s = \Gamma_r$,  and that $ \frac{k \,  \varepsilon}{2}$ satisfies \eqref{bounddeltaM} for some $0< \delta <1$ and $k \varepsilon < 1/4$,  $ \varepsilon < (1-\gamma) \, \min \left\{ 1,  \text{dist}(\Gamma, z) \right \} $ for some $0 < \gamma <1 $, where $p_{nrs}$ are defined as in \eqref{pnrs}, $h_{nmrs}$ as in \eqref{hnmrs}, and let $C_{1,M},C_{2,M},C_{M}, k_{1,M}$ be as  in \eqref{constant1M}, \eqref{constant2M} \eqref{constant3M}, and \eqref{k1M}, respectively.
The two point scatterers with equal weights located at $z$ and $z+ \varepsilon v$ are $\alpha$-distinguishable with respect to $\Gamma$ if the following inequality holds:
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
 \frac{k \,  \varepsilon}{2}
\leq
\ \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right)  - \sqrt{  \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right)^2 +  \frac{16 (1- \alpha)}{ C_{M} (2 + \alpha)^2}  }   \,.
\label{conclusionM}
\eqn
\end{Theorem}

In order to make the conditions stated in the above theorem practical, we assume that $\delta = 1/2$ and provide, by considering \eqref{series1} and  \eqref{series2},  approximations for the coefficients with
{\footnotesize
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
C_{1,M} &\approx& p_{000}
\, , \quad
C_{2,M} \approx 2 p_{000} \, ,
\quad \\[2mm]
C_M &\approx &
13  + 8\left|k_{1,M} \right| +  2 k_{1,M}^2
+ \frac{  \left(|p_{1\,0\,0}| + |p_{-1\,0\,0}|  - \frac{2 |p_{0\,1\,0}|}{k}  \right) ^2 + 2 \left(|p_{1\,0\,0}| + |p_{-1\,0\,0}|  - \frac{2 |p_{0\,1\,0}|}{k}  \right) ( p_{000} + |k_{1,M} | ) }{p_{000}^2}. \quad
\eqnx
}  With the above approximations of the constants $C_{1,M},C_{2,M},$ and $C_M$ and estimate \eqref{conclusionM}, we obtain an \textbf{anisotropic} optimal density at any point $z$ and at any direction $v$ for contrast reconstruction from the near-field data in terms of $\alpha$-distinguishably as follows:
{\small
\begin{eqnarray}}  \def\eqn{\end{eqnarray}
 h_{z,v,M}(\alpha) &:= & \frac{2}{k}\min \Bigg\{ \frac{1}{4}, \sqrt{ \frac{k_{1,M}^2 }{16} + \frac{1}{4} } -  \frac{|k_{1,M}| }{4} ,
  \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right)  \notag \\[3mm]
   &&    - \sqrt{  \left( \frac{ 2  \left(p_{1\,0\,0} - p_{-1\,0\,0}  - \frac{2 p_{0\,1\,0}}{k}  \right)  }{p_{000}} + k_{1,M} \right)^2 +  \frac{16 (1- \alpha)}{ C_{M} (2 + \alpha)^2}  }
    \Bigg\}. \label{finalfunctionM}
\eqn
}The remarks in the previous subsection are also in order in this subsection.

\subsubsection{Combining the optimal mesh size function along two opposite directions}

Practically, one may combine the optimal mesh-size for $v$ and $-v$ since they are refined in the same direction, and we may consider the following function instead,
$$\widetilde{h_{z,v}} (\alpha) := \min\big\{ h_{z,v}(\alpha) , h_{z,-v}(\alpha) \big\} \, ,$$
as well as
$$\widetilde{h_{z,v,M}} (\alpha) := \min \big\{ h_{z,v,M}(\alpha) , h_{z,-v,M}(\alpha) \big\} \,. $$
This is to forget the information about separating a point forward or backward along the direction $v$.

\section{Adaptive refinement using the optimal mesh function} \label{sec:refine}

We can apply the function $\widetilde{h_{z,v}}(\alpha)$ to refine a mesh at a point $z$ along a direction $v$ (which we may refer to as a polarization direction).
Suppose we have a finite number of inhomogeneous inclusions  inside $\Omega$, and we implement the following procedures to recover all  the inhomogeneities:
\begin{enumerate}
\item[(i)] Select a coarse mesh with mesh size $h = \lambda/2 $;
\item[(ii)] Locate where the reconstruction has an elongated shape at a sampling point $z$, say, a rod with direction $v$;
\item[(iii)] Refine the mesh $\widetilde{h_{z,v}}(\alpha)$ at the point $z$ and the polarization direction $v$ for further reconstruction;
\item[(iv)] Repeat step (iii)  until all inclusions are separated.
\end{enumerate}
It is worth emphasizing that the purpose of steps (ii) and (iii) is to determine whether the inhomogeneities are actually two objects or only one long rod.

\section{Numerical illustrations} \label{sec:numerics}

\subsection{Behavior of the mesh-size function $\widetilde{h_{z,v}}(\alpha)$}

In order to understand better the behavior of the optimal mesh-size function, let us numerically illustrate some of the behaviors of $\widetilde{h_{z,v}}(\alpha)$ defined as in \eqref{finalfunction} (see Theorem \ref{important_theorem}).
Assume that the sources $x_s$ and receivers $x_r$ are located at $\Gamma= \{ ( \cos( 4 \pi n/14) , 4 \sin( 2 \pi n/14)  )\} $ for $n=0,1,\cdots, 7$, as shown in Figure \ref{fig:gamma}(a), and let $\Omega = B_4(0)$.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.4\textwidth]{./gamma_location.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.4\textwidth]{./optimal_mesh_size_z_v_fixed_alpha_varies.eps}\hfill{}

  \hfill{}(a)\hfill{} \hfill{}(b)\hfill{}

 \caption{\label{fig:gamma}  \small{\emph{(a) $\Gamma$ in red. (b) The graph of $h_{z,v}(\alpha)$ with respect to $\alpha$ when $k = 2$, $z = (1,3.5)$ and the polarization direction $v = (1,0)$.  Measurement points in red and grid points in blue. }}}
 \end{figurehere}

We assume that $k = 2$, $z = (1,3.5)$, and the polarization direction $v = (1,0)$. The mesh-size function $h_{z,v}(\alpha)$ with respect to $\alpha$ is illustrated in Figure \ref{fig:gamma}(b).
We may also fix a large $\alpha$ to obtain a density function and uneven mesh with a mesh-size satisfying the density provided by $$ z \mapsto p_{v,\alpha} (z) :=  \frac{ \big(\widetilde{h_{z,v}}(\alpha)\big)^{-2} } { \int_{\Omega}   \big(\widetilde{h_{z,v}}(\alpha) \big)^{-2} dz}$$ for a fixed $v$.
The reason to take this function is that the optimal density is proportional to $ \big( \widetilde{h_{z,v}}(\alpha) \big)^{-M}$ where $M = 2$. Figure \ref{fig:meshsize2} shows the density function $ p_{v,\alpha} (z) $ at $v = (1,0)$ and $v=(0,1)$, respectively,  when $\alpha = 0.995$.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./density_v1.eps}\hfill{}
 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./density_v2.eps}\hfill{}
 \vskip -0.5truecm
 \caption{\label{fig:meshsize2}  \small{\emph{The function $ p_{v,\alpha} (z) $ at $v = (1,0)$ and $v=(0,1)$, respectively, when $\alpha = 0.99$. Measurement points in red and grid points in blue.  }}}
 \end{figurehere}

The above plots show that the refinement is unnecessary when the sampling points are at a distance half the radius from  the sources/receivers.
Now, we fix a particular $\alpha$ and generate a set of grid points with a mesh-size satisfying $h_{z,v}$ along the direction $v$.
Figure \ref{fig:meshsize3} illustrates the grid points with $v = (1,0)$ and $v=(0,1)$,  respectively,  when $\alpha = 0.55$.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_v1.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_v2.eps}\hfill{}
 \vskip -0.5truecm
 \caption{\label{fig:meshsize3}  \small{\emph{Grid points with mesh-size satisfying $h_{z,v}(0.1)$ along the direction $v$, where $v = (1,0)$ and $v=(0,1)$, respectively. }}}
 \end{figurehere}

One may observe that the refinements are done in the directions $v$ near the source/receiver points and are left coarser when the sampling points are further away.
In practice, we may combine the two mesh sizes and generate a tensor mesh with a mesh-size $\widetilde{h_{z,e_1}}(\alpha) $ along  $e_1 = (1,0)$ and $\widetilde{h_{z,e_2}}(\alpha) $ along $e_2 = (0,1)$.
Figure \ref{fig:meshsize4} respectively shows the grid points of the tensor mesh with  $\widetilde{h_{z,e_1}}(\alpha) $ along  $e_1 = (1,0)$ and $\widetilde{h_{z,e_2}}(\alpha) $ along $e_2 = (0,1)$ when $\alpha = 0.55$.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_v12.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_v3.eps}\hfill{}

  \hfill{}(a)\hfill{} \hfill{}(b)\hfill{}
 \caption{\label{fig:meshsize4}  \small{\emph{(a) Grid points with mesh-size satisfying $h_{z,v}(0.1)$ along the direction $v$, where $v = (1,0)$ and $v=(0,1)$, respectively. (b) Unrefined mesh.  Measurement points in red and grid points in blue. }}}
 \end{figurehere}

\subsection{Use of the mesh-size function $\widetilde{h_{z,v}}(\alpha)$ for improving the reconstruction methods}

This paper aims at getting an optimal mesh-size which maximizes the resolution of the reconstruction procedure with minimal numerical cost.  As discussed before, this choice of mesh is not restricted to one particular numerical method, and the numerical methods given in subsection \ref{example_haha} are just several examples to illustrate the concept.   In this subsection, we illustrate  the efficiency of the optimal meshing procedure using the DSM, which is introduced in subsection \ref{example_haha}.  For further details on the DSM, we refer to \cite{L1, L2, LZ,CILZ, L2, LXZ2, LXZ3, DSM11, DSM12, DSM13, DSM14, DSM15}.

In what follows, we consider an example of a near-field measurement corresponding to only one illumination.  Notice that since the case considered here is more pathological than the assumptions given in Theorem \ref{important_theorem}, the mesh-size function $\widetilde{h_{z,v}}(\alpha)$ defined as in \eqref{finalfunction} shall serve as a bottom line for the reconstruction.

We consider the case  $k = \pi^2 $ when $q(x)$ in \eqref{scattering1} is given by
\begin{eqnarray*}} \def\eqnx{\end{eqnarray*}
q(x) = 1 + 2 \mathcal{I}_{Q_1}(x) + 2 \mathcal{I}_{Q_2}(x),
\eqnx
where $Q_1 = \{ x \in \mathbb{R}^2 ; -0.275<x_1<-0.175, -0.275<x_2<-0.175  \} $ and
$Q_2 = \{ x \in \mathbb{R}^2 ; -0.15<x_1<-0.05, -0.15<x_2<-0.05  \} $
are two squares  inside $\mathbb{R}^2$, and $ \mathcal{I}_A$ is an indicator function of a set $A$ such that the value is $1$ if $x\in A$ and is $0$ otherwise.  We take a plane wave incidence with angle $\theta_y = \pi/4$, i.e., $d_y= (\cos(\theta_y), \sin(\theta_y))$.  As shown in the following figures,
$20$ measurement points along a circle of radius $0.5$ are employed in this example.
We perform the DSM with different mesh sizes to illustrate the efficiency of using $\widetilde{h_{z,v}}(\alpha)$ for meshing.
Note that we intentionally define $Q_1$ and $Q_2$ closer to the boundary measurement points, otherwise the density given by $\widetilde{h_{z,v}}(\alpha)$ will be similar to the situation when we have far-field measurements, and that will not serve the purpose of illustrating the efficiency of $\widetilde{h_{x,v}}(\alpha)$.

Moreover, we also intentionally put the distance between $Q_1$ and $Q_2$ smaller than half of the wavelength which is approximately $0.16$.  This helps to test if the meshing given by $\widetilde{h_{x,v}}(\alpha)$ can separate the two obstacles.
One remark is that, since we are taking near-field measurements, the DSM would have its kernal possessing ``heavy tail", (see, for instance, \cite{L} for more details), and hence the reconstruction will not be quite sharp and are expected to be more fussy.  We do not anticipate clear reconstructions since only one single measurement is applied in the recovery procedure.
This section, however, illustrates that, taking a mesh-size given by $\widetilde{h_{x,v}}(\alpha)$ helps to obtain as much details as those we can recover from our measurements.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_recon_1.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.45\textwidth]{./recon_1.eps}\hfill{}

 \hfill{}(a)\hfill{} \hfill{}(b)\hfill{}
 \vskip -0.2truecm
 \caption{\label{fig:reconstruction}  \small{\emph{(a) Grid points on a coarse mesh.  (b) reconstruction given by $ [I(x,dy) ]^2$, where $I_{\text{coarse}} := I(x,dy)$ in \eqref{eqn:indexfcn} with the grid point given by the coarse mesh. }}}
 \end{figurehere}

Figure \ref{fig:reconstruction} (a) illustrates the coarse mesh that we considered, with its mesh-size being half the wavelength.
 Figure \ref{fig:reconstruction} (b) shows the reconstruction given by the index function $I_{\text{coarse}} := I(x,dy)$ with this coarse mesh.  As we may expect, we cannot distinguish two inclusions in the image, and basically we only observe a large patch at the lower bottom part of the reconstruction.

Now, let us consider a refinement of the mesh as a tensor grid with mesh-sizes $\widetilde{h_{z,e_i}}(\alpha)$ along the two directions $e_1=(1,0)$ and $e_2=(0,1)$, where we take $\alpha = 0.9$.  We refine the mesh in the whole sampling domain $\Omega$, and obtain a fine mesh given as in Figure \ref{fig:reconstruction2} (a).  We then perform the DMS and Figure \ref{fig:reconstruction2} (b) is the approximation provided by the index function $I(x,dy)$ with this fine sampling.  We notice that this sampling helps to identify whether or not there are indeed two inclusions, with a clear dip separating the two inclusions.  However, we do not expect a very fine reconstruction since the problem is severely ill-posed.  Nonetheless, we are able to separate the two inclusions in the image.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_recon_2.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.45\textwidth]{./recon_2.eps}\hfill{}

 \hfill{}(a)\hfill{} \hfill{}(b)\hfill{}
 \vskip -0.2truecm
 \caption{\label{fig:reconstruction2}  \small{\emph{(a) Grid points on a refined mesh given by a tensor grid with mesh-sizes $\widetilde{h_{x,e_i}}(\alpha)$ along $e_1=(1,0), e_2=(0,1)$, $\alpha = 0.9$, and $x$ is taken in the whole sampling domain $\Omega$.  (b) A reconstruction given by $ [I(x,dy) ]^2$, where $I(x,dy)$ in \eqref{eqn:indexfcn} with a refined grid in the left. }}}
 \end{figurehere}

Next, we consider a refinement of the mesh only at the points where the value of $ I_{\text{coarse}} (x) > 0.995 $.
Again, we refine as a tensor grid with mesh-sizes $\widetilde{h_{z,e_i}}(\alpha)$ along the two directions $e_1=(1,0)$ and $e_2=(0,1)$, where we take $\alpha = 0.9$.  The resulting adaptive refinement is shown in Figure \ref{fig:reconstruction3} (a).  We perform the DSM. Figure \ref{fig:reconstruction3} (b) illustrates the recovery obtained by the index function $I(x,dy)$ with this adaptive fine sampling.  Comparing the reconstruction in Figures \ref{fig:reconstruction2} and \ref{fig:reconstruction3}, we can clearly observe  that the qualities of the reconstructions  are similar, i.e., both of them separate the two inclusions, but adaptive refinement requires a significantly smaller number of sampling points in the reconstruction.

\begin{figurehere}

 \hfill{}\includegraphics[clip,width=0.45\textwidth]{./mesh_recon_3.eps}\hfill{}
\hfill{}\includegraphics[clip,width=0.45\textwidth]{./recon_3.eps}\hfill{}

 \hfill{}(a)\hfill{} \hfill{}(b)\hfill{}
 \vskip -0.2truecm
 \caption{\label{fig:reconstruction3}  \small{\emph{(a) Grid points on adaptively-refined mesh given by a tensor grid with mesh-sizes $\widetilde{h_{x,e_i}}(\alpha)$ along $e_1=(1,0), e_2=(0,1)$, $\alpha = 0.9$, and $x$ is taken only when $I_{\text{coarse}}(x) > 0.995$.  (b) A reconstruction given by $ [I(x,dy) ]^2$, where $I(x,dy)$ in \eqref{eqn:indexfcn} with adaptively-refined grid on the left. }}}
 \end{figurehere}

\section{Conclusion} \label{con}

In this paper, we have discussed the choice of optimal mesh size for the reconstruction of inclusions from both the near-field and the far-field data in two and three dimensions.  The computation complexity would be sharply reduced with an optimal mesh size for state-of-the-art algorithms (i.e., LSM, DSM, CSIM, etc,) since the computation complexity is usually at least of the order of $\mathcal{O}(N^{M-1})$ in $\mathbb{R}^M$, where $N$ is the number of sampling points along one direction.  Moreover, an optimized mesh size can help improving reconstructions of some iterative reconstruction approaches (i.e., MSM, EMSM, etc,)
by avoiding choices of stopping criterion which may heavily depend on the quality of the reconstruction under a coarser mesh and a subjective choice of the cut-off value.
We can conclude that our results are expected to have important implications in solving inverse scattering problems.

\section{Acoknowledgement}
The work of Yat Tin Chow was supported by the ONR grant N000141712162 and the NSF grant DMS-1720237.
The work of Keji Liu was supported by the NNSF of China under grants Nos. 11601308 and 11771349.


\section{Introduction}
Several fields come together in private machine learning: cryptography, machine learning, distributed systems, and high-performance computing. As a result, researchers and practitioners developing scalable solutions may find themselves faced with a task requiring many diverse skill sets and expertise.

Adapting machine learning models in a way that allows for privacy-preserving prediction and training is complex and non-trivial on its own, and often requires deeper insight into cryptography or machine learning than experts from either field have about the other.
For instance, one may want to mix secure computation protocols with complementary properties and experiment with various cryptographic optimizations such as vectorization and other specializations. Simultaneously, one may want to find variants of machine learning components that give sufficient accuracy yet utilize the fact that some secure operations are significantly faster than others. Thus, modularity and extensibility are crucial for effective experimentation.

One may additionally be faced with a large implementation challenge, potentially having to stitch together several otherwise independent systems or re-implement complex methods from both fields. One must not only optimize local computations by taking advantage of diverse multi-core architectures, but also the distributed processes inherent in secure computation, including how best to orchestrate the joint execution on multiple machines and minimize the overhead of sending large amounts of data across the network. Building all of this at the right level of abstraction can be overwhelming, and often comes at the cost of extensibility and code readability. This makes experimentation harder, potentially harming accessibility and correctness.

Finally, for even small-scale solutions it is highly valuable to have access to tools for visual inspection, debugging, and profiling such as TensorBoard in order to identify issues and bottlenecks in both the protocol and the machine learning model. Lack of such tools represents an additional barrier to entry.

In this paper we illustrate the benefits of embedding protocols directly into TensorFlow in order to show that it can serve as a platform for easily experimenting with secure computation for private machine learning.

\setlength{\intextsep}{0pt}%
\setlength{\columnsep}{15pt}%
\begin{wrapfigure}{R}{0.38\textwidth}
\centering
\lstset{
    language=Python,
    morekeywords={with},
    basicstyle=\fontsize{6}{7}\selectfont\ttfamily,
    keywordstyle=\color{blue}\ttfamily,
    stringstyle=\color{red}\ttfamily,
    commentstyle=\color{darkgreen}\ttfamily,
    breaklines=false,
}
\vspace{-5px}
\lstinputlisting{pred.txt}
\caption{Example using \texttt{tf-encrypted} for private prediction,
with the prediction input known only in plaintext by the client and the model weights only by the owner.}
\label{fig:pred}
\vspace{-5px}
\end{wrapfigure}

\subsection{Contributions}

We present \texttt{tf-encrypted}, an open source library built on top of TensorFlow with the aim of making private machine learning more accessible to researchers and practitioners coming from either cryptography or machine learning, and without the need to be an expert in distributed systems or high-performance computing. A small example is given in Figure~\ref{fig:pred}, where a logistic regression model is used for making private predictions. Note that the methods processing input and output specify on which party they should run, and that usual TensorFlow operations may be used within to perform local plaintext operations.

To this end we adapt and implement a state-of-the-art secure computation protocol for tensor oriented applications (Section~\ref{sec:pond}). We report on benchmarks using common models from the literature (Section~\ref{sec:experiments}) and highlight additional properties of this approach that we find of value:
\textbf{usability:} by leveraging TensorFlow we obtain a familiar and comprehensive platform for building scalable solutions;
\textbf{integration:} by reducing all secure computations to TensorFlow graphs, it becomes straight-forward to mix these with ordinary computations for hybrid approaches;
\textbf{extensibility:}
using TensorFlow's high-level abstractions makes it easier to experiment with and develop new secure protocols on top of optimized primitives while maintaining code readability;
\textbf{performance:} we reach high runtime efficiency without sacrificing other properties via TensorFlow's distributed execution engine heavily optimized for networking, parallel execution, and scalability;
\textbf{benchmarking:} combining all of the above we obtain a common framework for comparable private machine learning.

\subsection{Related Work}

Several freely available implementations of secure computation protocols exist, including those of~\cite{SCALEMAMBA, FRESCO, DSZ15, OblivC}, yet all of these are standalone frameworks that do not provide integration with existing machine learning platforms, and arguably aim more at general purpose secure computation than our focus on private machine learning. As a result, users are faced with a lower level interface and may have to implement basic machine learning components. We are furthermore not aware of any development tools for these frameworks outside of general purpose debuggers.
To the best of our knowledge no machine learning support tools exist for these platforms.

The works of~\cite{SecureML, Gazelle, SecureNN, ABY3} focus on adapting secure computation protocols to private machine learning, in some cases using similar optimizations as the protocol presented in Section~\ref{sec:pond}. However, to the best of our knowledge none of these have openly available implementations and hence require a significant investment from anyone wanting to apply them. To a large extent our aim is to provide a common platform based on TensorFlow for implementing and experimenting with protocols such as these.
The concurrent work of~\cite{PySyft} takes a somewhat similar approach yet currently focuses on PyTorch as opposed to TensorFlow.

Another line of work has focused on using differential privacy for privacy-preserving machine learning~\cite{PAEGT17, PSMRTE18}, with reference implementations available in TensorFlow. These works remain orthogonal to our approach and do not employ any form of secure computation.

\section{Secure Computations in TensorFlow}
\label{sec:pond}

TensorFlow as described by~\cite{TensorFlow16} is among the leading frameworks for constructing and deploying machine learning models, offering an optimized engine for executing local and distributed computations as well as a high-level interface for expressing these. Importantly, the latter abstracts away lower-level operations such as networking while remaining powerful enough to succinctly express computations with operations and tensors pinned to specified machines (see for example Figure~\ref{fig:mul}). This combines for a powerful platform for giving efficient implementations of distributed computations, including complex secure computations where it is crucial for privacy that some data remain known only on select machines.

Moreover, the link between TensorFlow's engine and high-level interface takes the form of stateful dataflow graphs, with nodes for performing tensor operations in the distributed setting. This means that the engine can not only take advantage of optimizations such as lazy evaluation and multi-core processing tailored for the specific runtime environment, but also optimize the graph itself and chose a node execution order based on both static and runtime information (see e.g.~\cite{TensorFlow15}). For instance, we observe that large network transfers from Beaver triple generation are often automatically batched and moved to the beginning of the execution.

As an example of a secure computation protocol implemented in TensorFlow, we here outline our variant of the well-known SPDZ protocol by~\cite{SPDZ12} with two servers $\server{0}$ and $\server{1}$ that is used for benchmarking in Section~\ref{sec:experiments}. The protocol is vectorized to improve performance of applications relying heavily on tensor operations and we use generalized triples produced by an independent third server\footnote{This server is essentially the crypto producer as used in e.g.~\cite{SecureML} that can run entirely offline as long as the function to compute is known, which can be fully determined at compile time due to TensorFlow's use of static dataflow graphs.} to avoid sending redundant data when possible. Any number of input providers and output receivers are supported, holding e.g. training data or prediction inputs.  We currently ensure passive (honest-but-curious) security under a single corruption and rely on two cryptographic primitives, namely
additive secret sharing and secure channels between all players.

\setlength{\intextsep}{0pt}%
\setlength{\columnsep}{15pt}%
\begin{wrapfigure}{R}{0.42\textwidth}
\centering
\lstset{
    language=Python,
    morekeywords={with},
    basicstyle=\fontsize{6}{7}\selectfont\ttfamily,
    keywordstyle=\color{blue}\ttfamily,
    stringstyle=\color{red}\ttfamily,
    commentstyle=\color{darkgreen}\ttfamily,
    breaklines=false,
}
\vspace{-5px}
\lstinputlisting{mul.txt}
\caption{Secure multiplication implemented in TensorFlow: \texttt{tf.device} pin data and operations to specific machines to multiply two masked tensors, and takes care of implicitly adding nodes for transmitting \texttt{ab0} and \texttt{ab1}.}
\label{fig:mul}
\vspace{-5px}
\end{wrapfigure}

Following typical practice we use a fixed point encoding for the floating point numbers commonly used in machine learning, i.e. we scale by a fixed factor and treat the result as an integer. To represent these we support both a fixed \texttt{int64} and a CRT-based \texttt{int100} tensor type, the former offering higher performance and the latter higher precision. To maintain precision after multiplications we implement both the conservative truncation protocol of~\cite{CS10} requiring one round of communication, as well as the optimistic non-interactive protocol of~\cite{SecureML} that may fail with a small probability.

As in other SPDZ variants, \textit{private tensors} $\private{x}$ are secret shared into two tensors $\share{0}{x}$ and $\share{1}{x}$, held by server $\server{0}$ and $\server{1}$ respectively, such that $x = \share{0}{x} + \share{1}{x}$ yet either share on its own reveals nothing about $x$. However, unlike other variants, we also rely on a \textit{masked tensor} $\masked{x}$ which in addition to $\share{0}{x}$ and $\share{1}{x}$ also includes a random tensor $a$ held by $\server{2}$, shares $\share{0}{a}$ and $\share{1}{a}$ of it held by $\server{0}$ and $\server{1}$, and an $\alpha$ held by both $\server{0}$ and $\server{1}$ such that $\alpha = x - a$. While this is simply an explicit representation of the intermediate state of a SPDZ multiplication, having it in this form allows us to easily extend and generalize triples in order to reduce computation and networking.
Converting a tensor from private to masked takes one round of \emph{offline} communication where $a$ is sampled by $\server{2}$ and shares $\share{0}{a}$ and $\share{1}{a}$ of it are sent to $\server{0}$ and $\server{1}$, and one round of \emph{online} communication where $\share{i}{x} - \share{i}{a}$ is sent by $\server{i}$ to $\server{1-i}$ for $i \in \{0, 1\}$ and $(\share{0}{x} - \share{0}{a}) + (\share{1}{x} - \share{1}{a}) = \alpha$ is computed by both.

Secure computation on tensors proceed as in other SPDZ variants, albeit with some operations first converting private tensors into masked tensors. For instance,
multiplication is only implemented for masked tensors as $\private{z} = \mathsf{mul}(\masked{x}, \masked{y})$ with $\share{i}{z} = i \alpha^x \alpha^y + \alpha^x \share{i}{a}^y + \share{i}{a}^x \alpha^y + \share{i}{a^x a^y}$ and where $a^x a^y$ is computed and shared by $\server{2}$ (see Figure~\ref{fig:mul}).
Note that the result of a multiplication is a private tensor, meaning $z$ will have to be masked before it can be used as input to another multiplication; however, both $\masked{x}$ and $\masked{y}$ can readily be used again as is, meaning every tensor only needs to be masked once.
While special triples for squaring are no longer needed with explicit masking, we can still optimize e.g. matrix multiplications and convolutions via specialized triples and avoid reducing everything to (scalar) multiplications, in turn reducing networking further. Finally, operations such as transposing and stacking are done locally without interaction by letting all three servers operate on the values associated with private and masked tensors.

\section{Experiments}
\label{sec:experiments}

We benchmark the protocol of Section~\ref{sec:pond} for private inference on the typical MNIST handwritten digit classification task. Using TensorFlow we train each of the neural networks\footnote{These are variants of models studied in \cite{SecureML, Gazelle, SecureNN}. We use polynomials interpolated to fit ReLU on the interval $[-3,3]$ for activation functions and perform the final argmax and softmax on logits in plaintext to avoid computing these securely.} in Figure~\ref{fig:networks} on the plaintext training set and then run private inferences on the remaining test set, keeping both prediction input and model weights private. We perform all experiments on the Google Cloud Platform using instances in the same region (\texttt{us-east1}) and with 36 vCPUs/60 GB memory each. We note that TensorBoard was invaluable in this process for picking the right approximation intervals for activation functions and inspecting overall correctness, and that performing all operations in TensorFlow simplified the process of model handling significantly.

\begin{wrapfigure}{L}{0.55\textwidth}
\small
\centering
\input{experiments-networks}
\caption{Neural network architectures. The convolutional layers are denoted by $(\textit{field size}, \textit{channels}, \textit{stride}, \textit{padding})$ and average pooling layers by $\textit{window size}$}
\label{fig:networks}
\end{wrapfigure}

The left part of Table~\ref{tab:perf} summarizes the combined offline and online runtime averaged over 100 inferences. As seen we get reasonable performance that may already be adequate for concrete applications. The fact that \texttt{int100} also achieves good performance suggests that we are not limited by the lower precision of \texttt{int64} when looking at larger models. Additional experiments furthermore indicate sub-linearly scaling with respect to typical batch sizes, leading to an interesting trade-off between latency and through-put; for instance, Network C with batch sizes 1, 10, and 100 take respectively 124, 182, and 541ms.

We also compare our accuracy over the entire testing set against plaintext TensorFlow. As seen in the right part of Table~\ref{tab:perf} we obtain almost identical accuracy using both \texttt{int64} and \texttt{int100}, with indications that the latter may give a slightly better output distribution according to the mean KL divergence, in turn potentially justifying its higher runtime cost. We note that these models achieve good performance despite using approximations, further underlining the importance of being able to adapt models to the encrypted setting.

Based on the above we conjecture that our approach scales to larger models, and is at least within an order of magnitude of related work in terms of runtime performance (at what we believe to be a lower implementation cost). We defer proper exploration of both topics to the full version of this paper.

\begin{table}
\small
\centering
\input{runtime-benchmark}
\caption{Runtime benchmarks with \texttt{int64} and \texttt{int100}.}
\label{tab:perf}
\end{table}

\section{Conclusion}

We have proposed an open source framework for experimenting with secure computation in TensorFlow, and illustrated how implementation of such protocols can be easily expressed using high level abstractions (Figure~\ref{fig:mul}). This additionally allows private machine learning to be expressed in an interface similar to ordinary TensorFlow (Figure~\ref{fig:pred}) while maintaining good performance.

In the full version of this paper we elaborate on these results, and present a modular extension of the concrete protocol presented here that adds features from~\cite{SecureNN} in order to compute exact ReLU and MaxPooling.

\bibliographystyle{plainnat}


\section{Introduction}
In thermoelectric materials, the concentration of free charge carriers needs to be highly tunable to allow optimization of the thermoelectric figure of merit \textit{zT}.\cite{toberer2008,gorainrm2017} Therefore, the ability to dope a material to the desired charge carrier type ($n$ or $p$) and carrier concentration plays a vital role. Beyond thermoelectrics, dopability of semiconductors is critical for optoelectronic\cite{nakamura1994,hosono2017} and power electronic applications.\cite{hudgins2003}

Mg$_3$Sb$_2$-based thermoelectrics have recently risen to prominence almost exclusively due to the success in $n$-type doping that was key to achieving high \textit{zT}.\cite{tamaki2016,ohno2017,zhang2017exp} Additionally, the earth-abundance and non-toxicity of the constituent elements make Mg$_3$Sb$_2$ a very attractive material for thermoelectric applications. The $n$-type doping of Mg$_3$Sb$_2$, which has been achieved only recently\cite{tamaki2016,ohno2017,zhang2017exp}, requires Mg-rich conditions during growth to suppress the formation of compensating, acceptor-behaving Mg vacancies. In addition, one needs to introduce extrinsic dopants, most commonly chalcogens such as Se and  Te, which behave as donor defects when substituting Sb. Combined with Bi alloying to reduce the lattice thermal conductivity, a high \textit{zT}\textgreater1.5 has been demonstrated in Te-doped Mg$_3$Sb$_2$-based materials.\cite{ohno2017,imasato2018}

Using first-principles defect calculations, we have previously examined the role of native defects and growth conditions in the dopability of Mg$_3$Sb$_2$.\cite{ohno2017,Gorai2018} We found that under Mg-poor conditions, the predominant defects are acceptor Mg vacancies that inhibit $n$-type doping, consistent with the observations in previous works\cite{tamaki2016}. In contrast, under Mg-rich conditions, donor-behaving Mg interstitials are the predominant defects. As a result of the native defect energetics, we found that Mg$_3$Sb$_2$ is natively self-doped $p$-type under Mg-poor and $n$-type under Mg-rich growth conditions. More importantly, the low concentrations of acceptor Mg vacancies under Mg-rich conditions allows extrinsic $n$-type doping with Te\cite{ohno2017,imasato2018}, with predicted free electron concentrations exceeding 10$^{19}$ cm$^{-3}$ at synthesis temperatures of $\sim$900 K, consistent with the experimental findings\cite{ohno2017}.

In addition to Sb substitution with Se and Te, we have previously investigated\cite{Gorai2018} other $n$-type doping strategies for Mg$_3$Sb$_2$, including: (a) Sb substitution with halides Br and I, (b) Mg substitution with trivalent cations La, Al, Nb, and Ga, and (c) insertion of cation interstitials Li, Zn, Cu, and Be. We found La to be the most effective $n$-type dopant of all the studied elements, and Mg substitution the most effective doping strategy. According to our calculations, La preferentially substitutes the Mg atoms (labelled Mg(1)) located between the [Mg$_2$Sb$_2$]$^{2-}$ layers, as shown in Figure \ref{structure}. The predicted electron concentration in La-doped Mg$_3$Sb$_2$ is $\sim$5$\times$10$^{20}$ cm$^{-3}$ (at 900K), which is almost an order of magnitude higher than for Te doping. This theoretical prediction has been experimentally validated;\cite{imasato2018} La is found to be an effective $n$-type dopant in Mg$_3$Sb$_{1.5}$Bi$_{0.5}$ with free electron concentrations nearly twice of those achieved through Te doping. Addition of La is shown to also significantly improve the thermal stability of $n$-type Mg$_3$Sb$_{1.5}$Bi$_{0.5}$. \cite{imasato2018}

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{mg3sb2structure.pdf}
\caption{\label{structure}
Mg$_3$Sb$_2$ contains two unique Mg Wyckoff positions denoted by Mg(1) and Mg(2), and one unique Sb Wyckoff position. The crystal structure of Mg$_3$Sb$_2$ can be visualized as slabs of [Mg$_2$Sb$_2$]$^{2-}$ intercalated with Mg$^{2+}$ cations. An effective $n$-type doping strategy involve Mg substitution with tri-valent group-3 elements (La, Sc, Y).}
\end{figure}

In this work, we extend our investigation of extrinsic $n$-type doping of Mg$_3$Sb$_2$ to other group-3 elements such as Sc and Y. We predict that, in addition to La, both Sc and Y are effective $n$-type dopants. We find that both Sc and Y are more effective than Te in terms of the free electron concentrations, which are predicted to exceed 10$^{20}$ cm$^{-3}$ (900 K); Y is as effective as La while Sc only slightly lesser. Overall, we find that out of the 15 different extrinsic dopants we have studied, Mg substitution with group-3 elements (Sc, Y, La) is by far the most effective $n$-type doping strategy.

\section{Computational Methods}
First-principles point defect calculations provide formation energies of native defects and extrinsic dopants as functions of the Fermi energy. The achievable dopant and charge carrier concentrations in semiconductors can be predicted by combining the defect formation energetics with thermodynamic modeling of defect and charge carrier equilibria \cite{Goyal2017} We calculated the defect formation energies in Mg$_3$Sb$_2$ using density functional theory (DFT) and the standard supercell approach.\cite{lany2008} Within the supercell approach, the formation energy ($\Delta H_{D,q}$) of a point defect $D$ in charge state $q$ is calculated as:
\begin{equation}\label{defectenthalpy}
\Delta H_{D,q} = (E_{D,q} - E_H) + \sum_{i} n_i\mu_i + qE_F + E_{corr}
\end{equation}
where $E_{H}$ and $E_{D,q}$ represent the total energies of the defect-free, host supercell with no net charge ($E_H$) and the supercell with defect $D$ in charge state $q$, respectively. The chemical potential of element $i$ is denoted by $\mu_i$ and $n_i$ is the number of atoms of element $i$ added ($n_i$\textless 0) or removed ($n_i$\textgreater 0) from the supercell. $E_F$ is the Fermi energy, which represents the reservoir of charge (the Fermi sea). The term $qE_F$ in Eq. \ref{defectenthalpy} is the characteristic energy of exchanging charge between the defect and the reservoir of charge. The supercell approach to calculating defect energetics suffers from artifacts arising due to finite size effects. Additional artifacts are introduced due to the limitations of DFT, most notably, the underestimation of the band gap with standard functionals such as GGA-PBE.\cite{Perdew1996} Various correction schemes are available to correct for the finite size and insufficiently accurate electronic structure artifacts; these corrections, which are briefly discussed in the following paragraphs, are represented by the term $E_{corr}$. Detailed descriptions of the correction scheme can be found in Ref.~\citenum{lany2008}.

The total energies of the supercells we calculate using the generalized gradient approximation of Perdew-Burke-Ernzerhof (PBE) \cite{Perdew1996} within the projector augmented wave (PAW) formalism as implemented in the VASP software package.\cite{Kresse1996} A 3$\times$3$\times$2 supercell containing 90 atoms is used to calculate the defect energetics. The total energies are calculated with a plane-wave energy cutoff of 340 eV and a $\Gamma$-centered 4$\times$4$\times$4 Monkhorst pack $k$-point grid to sample the Brillouin zone. The defect supercells are relaxed following a similar procedure used in Refs. \citenum{ortiz2017,Gorai2018,ohno2017}. Specifically,

The elemental chemical potentials $\mu_i$ are expressed relative to those of the elements in reference elemental phases as $\mu_i=\mu_i^0 + \Delta\mu_i$, where $\mu_i^0$ is the reference chemical potential under standard conditions and $\Delta\mu_i$ is the deviation from the reference. $\Delta\mu_i=0$ corresponds to $i$-rich conditions. For example, $\Delta\mu_\mathrm{Mg}=0$ corresponds to the equilibrium between Mg$_3$Sb$_2$ and Mg-metal conditions (referred to as Mg-rich). Following the FERE approach,\cite{stevanovic2012} the reference chemical potentials ($\mu_i^0$) are fitted to a set of measured formation enthalpies of compounds. The stability Mg$_3$Sb$_2$ relative to decomposition into competing phases determines the region of phase stability and the bounds on the values of $\Delta\mu_i$.

The underestimation of the band gap in DFT is remedied by applying individual valence and conduction band edge shifts (relative to the DFT-computed band edges) as determined from GW quasi-particle energy calculations \cite{lany2008}. We use DFT wave functions as input to the GW calculations. The GW eigen-energies are iterated to self-consistency to remove the dependence of the G$_0$W$_0$ result on the single-particle energies of the initial DFT calculation. The input DFT wave functions are kept constant during the GW calculations, which allows the interpretation of the GW quasi-particle energies in terms of energy shifts relative to the DFT Kohn-Sham energies. Relative to DFT-calculated conduction and valence band edges, the shifts are +0.176 eV (shift to higher energy) and -0.171 eV (shift to lower energy), respectively.

The finite-size corrections that are included in $E_{corr}$, following the methodology in Ref. \citenum{lany2008}, are: (1) image charge correction for charge defects, (2) potential alignment correction for charged defects, (3) band filling correction for shallow defects, and (4) correction of band edges for shallow acceptors/donors. The calculations are organized and the results are analyzed using our software package, pylada-defects, for automation of point defect calculations.\cite{Goyal2016}

Under a given growth condition, the equilibrium $E_F$ and the corresponding free charge carrier concentration is determined by solving charge neutrality. The concentration of donor and acceptor defects are determined using Boltzmann distribution, such that $\left[D_q\right] = N_{s} e^{-\Delta H_{D,q}/\mathrm{k_B}T}$, where $\left[D_q\right]$ is the defect concentration, $N_s$ is the concentration of lattice sites where the defect can be formed, $\mathrm{k_B}$ is the Boltzmann constant, and T is the temperature. At a given T, The concentrations of electrons and holes are functions of $E_F$. To establish charge neutrality, the total positive charges should equal the negative charges. In this equation, $E_F$ is the only free parameter. By solving charge neutrality condition self-consitently, we can determine the equilibrium $E_F$ and the corresponding free carrier concentration.

\section{Results and Discussion}
\begin{figure}[!t]
\centering
\includegraphics[width=0.9\linewidth]{native.pdf}
\caption{\label{native}
Formation energetics of native defects in Mg$_3$Sb$_2$ under (a) Mg-poor, and (b) Mg-rich growth conditions. Subscripts (1) and (2) denote the unique Wyckoff positions of Mg as shown in Figure \ref{structure}. The solid vertical line labeled ``cbm'' is the conduction band minimum. A hypothetical donor dopant is shown with a dotted black line.}
\end{figure}
In previous studies, we have investigated the defect energetics of native point defects \cite{ohno2017} and a suite of candidate $n$-type dopants \cite{Gorai2018}. We briefly revisit these findings in Section \ref{sec:native} to contextualize the results of $n$-type doping with Sc and Y in Section \ref{sec:doping}.

\subsection{Energetics of Native Defects}\label{sec:native}
\begin{figure}[t!]
\centering
\includegraphics[width=0.9\linewidth]{g3doping.pdf}
\caption{\label{scydoping}
Defect energetics of $n$-type doping of Mg$_3$Sb$_2$ by cation substitution under Mg-rich ($\Delta\mu_\mathrm{Mg}$ = 0 eV) and most dopant-rich conditions. The formation energies of native defects and Te doping on Sb site (Te$_\mathrm{Sb}$) are also shown for reference. Doping with group-3 elements: (a) Sc, and (b) Y. The defects corresponding to the extrinsic dopants (Sc$_\mathrm{Mg}$, Y$_\mathrm{Mg}$) are denoted by black dotted lines. Both Sc and Y are found to be more effective $n$-type dopants compared to Te.}
\end{figure}
The formation energetics of native point defects critically influences the dopability of semiconductors, as exemplified in the case of Mg$_3$Sb$_2$. The calculated energetics of native point defects, under Mg-poor and Mg-rich growth conditions, as functions of Fermi energy ($E_F$) in Mg$_3$Sb$_2$ are shown in Figure \ref{native} (adapted from Ref. \citenum{Gorai2018}). Defects plots with positive slopes are donors and negative slopes are acceptors. We found that all native point defects in Mg$_3$Sb$_2$ are shallow; Mg vacancies are the only acceptor defects while all other native defects are donors. More importantly, we found that the growth conditions (Mg-rich \textit{vs.} Mg-poor) have a profound effect on the dopability of Mg$_3$Sb$_2$. Specifically, under Mg-poor conditions, acceptor Mg vacancies are the dominant native defects that hinder $n$-type doping of Mg$_3$Sb$_2$ by compensating or ``killing'' electrons generated by extrinsic donor dopants. This is schematically illustrated in Figure \ref{native}(a), where the electrons generated by the hypothetical donor dopant (dotted black line) is compensated by holes generated by the acceptor Mg vacancies such that the equilibrium $E_\mathrm{F}$ is near mid-gap with low free electron concentration. Commonly, Mg$_3$Sb$_2$ has been synthesized under Mg-poor growth conditions due to the inherently high volatility of Mg; consequently, Mg$_3$Sb$_2$ has been, until recently, synthesized almost exclusively as a $p$-type material \cite{ponnam2013,condron2006,bhardwaj2014}.

In contrast, under Mg-rich growth conditions (Figure \ref{native}b), the formation energy of Mg vacancies are significantly higher and therefore, present in low concentrations. The predominant defects are Mg interstitials that occupy octahedral sites within the [Mg$_2$Sb$_2$]$^{2-}$ layers. Therefore, Mg-rich growth conditions provide an opportunity to extrinsically $n$-type dope  Mg$_3$Sb$_2$; the hypothetical donor dopant (dotted black line) shown in Figure \ref{native}(b) does not suffer from electron compensation.

\subsection{Energetics of Sc and Y Doping}\label{sec:doping}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{doping-comparison.pdf}
\caption{\label{dopingcomparison}
Comparison of predicted free electron concentrations (at growth temperature 900 K) for 15 different $n$-type dopants for Mg$_3$Sb$_2$. The free electron concentration in natively self-doped Mg$_3$Sb$_2$ is denoted by ``nd''. The free electron concentrations for all but Sc and Y doping are extracted from the results in Refs. \citenum{Gorai2018} and \citenum{ohno2017}. Group-3 elements (Sc, Y, La) are more effective compared to chalcogens (Se, Te) that are traditionally used for $n$-type doping of Mg$_3$Sb$_2$.}
\end{figure}
The defect energetics of $n$-type doping of Mg$_3$Sb$_2$ with Sc and Y are shown in Figure \ref{scydoping}. Both Sc and Y are effective Mg-site $n$-type dopants due to the low formation energy of substitutional defects Sc$_\mathrm{Mg(1)}$ and Y$_\mathrm{Mg(1)}$. The formation energy of substitutional Te$_\mathrm{Sb}$ in Figure \ref{scydoping} is reproduced from Refs. \citenum{Gorai2018} and \citenum{ohno2017}. We find that $n$-type doping of Mg$_3$Sb$_2$ with Sc and Y is more effective than with chalcogens. The effectiveness is defined as the formation energy of Sc and Y impurities arouns the conduction band minimum. Lower the defect formation energy, higher the solubility of the impurity and consequently, higher the achievable free electron concentration.

As illustrated in Figure \ref{structure}, there are two unique Wyckoff positions of Mg atoms in Mg$_3$Sb$_2$ that are denoted by Mg(1) and Mg(2), where Mg(1) is the Mg atom residing between the [Mg$_2$Sb$_2$]$^{2-}$ slabs. From Figure \ref{scydoping}, we find that Sc and Y preferentially substitutes Mg(1) over Mg(2). The tendency to preferentially substitute Mg(1) is similar to La doping\cite{Gorai2018} and may be attributed to the large ionic size of Sc$^{3+}$ and Y$^{3+}$ ions that can be spatially accommodated betewen the [Mg$_2$Sb$_2$]$^{2-}$ slabs. The substitutional donor defects Sc$_\mathrm{Mg(1)}$ and Y$_\mathrm{Mg(1)}$ have lower formation energies than the donor Te$_\mathrm{Sb}$ (Figure \ref{scydoping}). Doping with Y, in particular, is more effective than Sc and comparable to La doping\cite{Gorai2018}. The predicted free electron concentrations at 900K under Mg-rich and most Y-rich conditions (as allowed by the phase stability of Mg$_3$Sb$_2$) exceeds 10$^{20}$ cm$^{-3}$.

The high thermoelectric performance of $n$-type Mg$_3$Sb$_2$ is attributed to its highly degenerate conduction band \cite{zhang2017theory} and low lattice thermal conductivity \cite{peng2018}. Heavy doping on the Mg site can have detrimental effects on the intrinsic electron transport properties since the conduction band edge is mainly composed of states derived from Mg($s$) orbitals. The detrimental effects of Mg site doping can be minimized by light doping with an efficitive $n$-type dopant, as experimentally demonstrated for La doping\cite{imasato2018}.

Figure \ref{dopingcomparison} compares the predicted free electron concentrations (at growth temperature 900 K) for a suite of 15 plausible $n$-type dopants in Mg$_3$Sb$_2$. It is evident from Figure \ref{dopingcomparison} that group-3 elements (Sc, Y, La) are far more effective $n$-type dopants compared to chalcogens (Se, Te). With group 3 elements, free electron concentrations exceeding 10$^{20}$ cm$^{-3}$ can be achieved. Recent experimental confirmation\cite{imasato2018} of effective $n$-type doping with La provides further confidence in our prediction about effective $n$-type doping with group-3 elements.

\section{Conclusions}
In this work, we have computationally assessed the effectiveness of $n$-type doping of Mg$_3$Sb$_2$ with group-3 elements such as Sc and Y. We find that both Sc and Y are more effective $n$-type dopants than chalcogens (Se, Te). The predicted doping behaviors of Sc and Y are similar to La; Y is more effective than Sc with predicted free electron concentrations comparable to those achieved by La doping. The higher doping effectiveness of Sc and Y offers greater tunability of the electronic properties in Mg$_3$Sb$_2$-based thermoelectric materials. Also, the relative abundance of group-3 elements (Sc, Y, La) compared to Te offers at additional incentive for considering these dopants. Following the successful experimental confirmation\cite{imasato2018} of our prediction about effective La doping\cite{Gorai2018} of Mg$_3$Sb$_2$-based materials, we call upon the experimentalists to confirm the new predictions presented in this work.\\

\section*{Acknowledgements}
We acknowledge support from NSF DMR program, grant no. 1729594. The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory (NREL).