A direct sampling method for simultaneously recovering inhomogeneous inclusions of different nature

In this work, we propose a novel parameter reconstruction method in which we decouple measurements from one (or at most two) pair(s) of Cauchy data and locate two different types of inhomogeneities in the model. Let us consider an open bounded domain $\Omega$ in $\mathbb{R}^{d}$ ($d=2$, $3$) with a smooth boundary $\partial\Omega$, and the following elliptic PDE:

$$\begin{cases}-\nabla\cdot(\sigma\nabla u)+Vu=0\quad&\text{in}\quad\Omega\,,\\ \frac{\partial u}{\partial\nu}=f\quad&\text{on}\quad\partial\Omega\,,\end{cases}$$

(1.1)

where the coefficients $\sigma$, $V\in L^{\infty}(\Omega)$ represent two unknown physical inclusions in the physical ranges $c<\sigma<C$ and $-C<V<C$ for some $c>0$, $C<\infty$. Let $\sigma_{0}$ and $V_{0}$ be the respective coefficients describing the homogeneous background medium $u_{0}$. We assume two physical inclusions are in the interior of the domain, i.e., $\text{supp}(\sigma-\sigma_{0})$, $\text{supp}(V-V_{0})$. Our goal is to simultaneously identify and reconstruct these two inclusions, i.e., $\text{supp}(\sigma-\sigma_{0})$ and $\text{supp}(V-V_{0})$, using the data $u$ measured on the boundary corresponding to a boundary influx $f$. We like to point out that our proposed method can be appropriately generalized to handle other types of boundary conditions that may arise in real applications, e.g. the Robin boundary condition, although this work focuses only on a Neumann boundary condition (cf. (1.1)).

Inverse problems of the elliptic system (1.1) may arise from a wide range of applications, such as medical imaging, geophysical prospecting, nano-optics, and nondestructive testing; see, e.g., [19, 33, 37, 41] and the references therein. The solution $u$ and two coefficients $\sigma$ and $V$ may represent different physical state and parameters in different applications. For instance, in the diffusion-based optical tomography [5], $u$, $\sigma$ and $V$ represent the photon density, diffusion and absorption coefficients, respectively; Identification of locations of inhomogeneities of $\sigma$ and $V$ helps determine the distribution of different types of tissues. The model (1.1) can also represent the inverse electromagnetic scattering problem. Under the transverse electric symmetry, the three-dimensional full Maxwell equations may be reduced to (1.1), where $\sigma$ and $-V$ stand for the permeability and permittivity of the media [39]. The system (1.1) is also adopted in the ultrasound medical imaging, where $\sigma$ and $V$ represent the volumetric mass density and bulk modulus, respectively, while $u$ describes the acoustic pressure [2]. For the convenience of descriptions, we shall often call $\sigma$ and $V$ as conductivity and potential throughout this work.

The uniqueness and simultaneous identifiability for the elliptic inverse problem (1.1)) have been widely investigated. In particular, a negative result was proved in [6], that is, no uniqueness for the simultaneous reconstruction of $\sigma$ and $V$ when both coefficients are smooth. For piecewise constant $\sigma$ and piecewise analytic $V$, the uniqueness and simultaneous identifiability were established in [23] for real-valued coefficients, as long as all possible Neumann-to-Dirichlet data are available. This uniqueness result also hints why there are many reasonable numerical results for simultaneous reconstructions, even though there is still no general uniqueness result.

During the recent two decades, many efficient numerical methods were proposed for the inverse problem (1.1). Minimizing a least-squared functional with appropriate regularizations is a very popular methodology in many applications, along with iterative methods; see, e.g., [7, 17, 18, 29, 40]. Usually, a locally convergent Newton-type method is employed. However, an iterative scheme may be trapped often in local optima, owing to high ill-posedness and high non-convexity of the objective functional. Moreover, the high dimension of the optimization problem also hinders the performance of this type of algorithms. Therefore there is a significant interest to develop some alternative numerical methods, that are fast, computationally cheap and robust against noisy data, to provide a reasonable initial guess for these iterative methods. On the other hand, some rough estimates of the inhomogeneous inclusions directly from the measurement data may be sufficient for many practical applications.

Motivated by these two important applications, many non-iterative schemes were developed for a large class of inverse problems for parameter identifications. Most of those methods are sampling-type, which rely on an appropriately designed functional that is expected to attain relatively large values inside the inhomogeneity. These include linear sampling method [15], singular source method [34], and factorization method [26], etc. Recently, MUSIC-type method using the multistatic response matrix (MSR) [4, 27], algorithms based on the topological derivative [8], and the reverse time migration [10] were also developed for the purpose. We refer to several recent monographs [9, 11, 28, 35] for more developments in this direction. Nevertheless, to the best of our knowledge, there seems to exist little development of sampling type methods for simultaneously reconstructing two different types of inhomogeneities.

In this work, we make the first effort to develop a new sampling type method, a direct sampling method (DSM), for simultaneously identifying and recovering multiple inhomogeneous inclusions corresponding to two different physical parameters. In particular, a specific attempt is made to ensure that the method can apply to the important scenarios where very limited data is available, e.g., only noisy data collected at one or two measurement event. DSMs have been developed recently through a series of efforts, e.g., [12, 13, 14, 24, 31, 34], for recovering the inhomogeneous media, first for the wave type inverse problems, and then for the non-wave inverse problems. This family of direct sampling methods construct an index function that leverages upon an almost orthogonality property between the family of fundamental’s functions of the forward problem and a particular family of probing functions under a properly selected Sobolev duality product. All the existing DSMs were designed for the cases when there are only inhomogeneous inclusions of same physical nature. In this work, we make the first attempt to design DSMs for simultaneously recovering multiple inhomogeneous inclusions of two very different physical parameters. These inverse scattering problems are much more ill-posed and challenging than those associated with inclusions of same physical nature. A natural mathematical and technical issue is how to identify which inclusions come from one physical parameter, not from the other; and how to locate and separate the multiple inclusions corresponding to one parameter from those corresponding to the other. We shall make use of an important observation that the near field or scattered data satisfies a fundamental property that it can be approximated as a combination of Green’s functions and their gradients at a set of discrete points. With this observation, we shall develop two separate families of probing functions, namely, the monopole and dipole probing functions, which enable us to construct two separate index functions for decoupling the multiple inhomogeneous inclusions associated with one physical parameter from those associated with the other parameter. In order for this decoupling to function effectively, we introduce a new and key concept, the mutually almost orthogonality property, between the family of fundamental functions and their gradients, and two families of monopole and dipole probing functions. Furthermore, we take advantage of an additional parameter, namely the probing direction of the dipole probing function, and an appropriate boundary influx to decouple the multiple inhomogeneous inclusions of one parameter from those of the other parameter. As we will see, the new method is computationally cheap and numerically stable, and works quite satisfactorily, as demonstrated in section 6 by several typical challenging numerical examples with very limited data available, e.g., only noisy data collected at one or two measurement event. The outputs generated by the new method can serve as reasonable approximations for many important applications where general rough locations and shapes of inhomogeneous inclusions are sufficient, or as a quick and stable initial guess of some expensive nonlinear optimization approaches when more accurate reconstructions are needed.

The rest of our work is as follows. We address in section 2 the general principles of DSMs, including the fundamental property and the new mutually almost orthogonality property. We then show in section 3 that the fundamental property holds for our inverse problem in many cases that we encounter in practice. We propose in section 4 two index functions for the reconstruction process and discuss their properties, including an alternative characterization. In section 5, we derive some explicit representations of the probing and index functions in some special sampling domains and discuss the mutually almost orthogonality properties in those cases. We will also address some appropriate boundary influxes to further decouple the monopole and the dipole effects in the measurement. Numerical experiments are conducted in section 6 to illustrate the effectiveness of the new method.

We briefly explain in this section some general observations that motivate our direct sampling method with coupled measurement. The development of our DSM hinges on a basic fact that our measurement data can be approximated by a sum of Green’s functions of the homogeneous equation and their gradients. With this in mind, along with an appropriate choice of the Sobolev duality product, those Green’s functions and their gradients located at different sampling points are respectively nearly orthogonal with two properly selected families of probing functions. These two families of probing functions are monopole-type and dipole-type functions, and couple well with the Green’s function and its gradient respectively. This is a very important property to our new method, and will be called the mutually almost orthogonality property, namely, the Green’s functions interact well only with monopole probing functions, while the gradient of Green’s functions interact well solely with dipole probing functions. This allows us to decouple the monopole and the dipole effects. Moreover, different types of boundary influxes and probing directions can be chosen to maximize the decoupling effect.

To be more precise, we aim to make use of the following two properties to develop an effective and robust direct sampling method:

1. 1.

   (Fundamental property) The boundary data, i.e., $u-u_{0}$ on $\partial\Omega$, of the model (1.1) can be represented approximately by a sum of Green’s functions of the homogeneous medium and their gradients:

   $$(u-u_{0})(x)\approx\sum_{j=1}^{n}c_{j}\,G_{q_{j}}(x)+\sum_{i=1}^{m}a_{i}\,d_{i}\cdot\nabla G_{p_{i}}(x)\,,\quad x\in\partial\Omega$$

   for some choices of coefficients $\{c_{j}\}_{j=1}^{n}\in\mathbb{C}$, $\{(a_{i},d_{i})\}_{i=1}^{n}\in\mathbb{C}\times\mathbb{S}^{d-1}$, and the sets of discrete points $\{q_{j}\}_{j=1}^{n}\in\text{supp}(V-V_{0})$, $\{p_{i}\}_{i=1}^{m}\in\text{supp}(\sigma-\sigma_{0})$.

2. 2.

   (Mutually almost orthogonality property) There are two sets of probing functions, namely $\{\zeta_{x}\}_{x\in\Omega}$ representing a family of monopole probing functions at sources $x\in\Omega$, and $\{\eta_{x,d}\}_{x\in\Omega,d\in\mathbb{S}^{d-1}}$ representing a family of dipole probing functions at sources $x\in\Omega$ and dipole directions $d\in\mathbb{S}^{d-1}$, such that the following four kernels

   	$\displaystyle(x,z)$	$\displaystyle\mapsto$	$\displaystyle{K_{1}}(x,z):=\frac{(\zeta_{x},G_{z})_{\text{mo}}}{C_{\text{mo}}(x)}\,,$
   	$\displaystyle(x,z,d_{z})$	$\displaystyle\mapsto$	$\displaystyle{K_{2,d_{z}}}(x,z):=\frac{(\zeta_{x},d_{z}\cdot\nabla G_{z})_{\text{mo}}}{C_{\text{mo}}(x)}\,,$
   	$\displaystyle(x,z,d_{x})$	$\displaystyle\mapsto$	$\displaystyle{K_{3,d_{x}}}(x,z):=\frac{(\eta_{x,d_{x}},G_{z})_{\text{di}}}{C_{\text{di}}(x,d_{x})}\,,$
   	$\displaystyle(x,z,d_{x},d_{z})$	$\displaystyle\mapsto$	$\displaystyle{K_{4,d_{x},d_{z}}}(x,z):=\frac{(\eta_{x,d_{x}},d_{z}\cdot\nabla G_{z})_{\text{di}}}{C_{\text{di}}(x,d_{x})}\,$

   have the following properties, under two appropriate couplings $(\cdot,\cdot)_{\text{mo}}$, $(\cdot,\cdot)_{\text{di}}$ and weights $C_{\text{mo}}(x)$ for $x\in\Omega$ and $C_{\text{di}}(x,d)$ for $x\in\Omega$, $d\in\mathbb{S}^{d-1}$:

   	$\displaystyle{K_{1}}(x,z)$		is of large magnitude if $x$ is close to $z$, and is small otherwise,
   	$\displaystyle{K_{2,d_{z}}}(x,z)$		is relatively small,
   	$\displaystyle{K_{3,d_{x}}}(x,z)$		is relatively small,
   	$\displaystyle{K_{4,d_{x},d_{z}}}(x,z)$		is of large magnitude if $x\approx z$ and $d_{x}\approx d_{z}$, and is small otherwise.

The above mutually almost orthogonality property means that the two families of probing functions, i.e., monopole and dipole probing functions, interact well with only the Green’s functions and their gradients respectively. This is a very important property that allows us to decouple the monopole and dipole effects in the measurement data.

With the above definitions and the fundamental property, we can define two index functions

$\displaystyle I_{\text{mo}}(x):=\frac{(\zeta_{x},u-u_{0})_{\text{mo}}}{C_{\text{mo}}(x)}\,\quad\text{and}\quad I_{\text{di}}(x,d_{x}):=\frac{(\eta_{x,d_{x}},u-u_{0})_{\text{di}}}{C_{\text{di}}(x,d_{x})}\,,$

(2.1)

which have the approximations

	$\displaystyle I_{\text{mo}}(x)$	$\displaystyle\approx$	$\displaystyle\sum_{j=1}^{n}c_{j}K_{1}(x,q_{j})+\sum_{i=1}^{m}a_{i}K_{2,d_{i}}(x,p_{i})\,,$
	$\displaystyle I_{\text{di}}(x,d_{x})$	$\displaystyle\approx$	$\displaystyle\sum_{j=1}^{n}c_{j}K_{3,d_{x}}(x,q_{j})+\sum_{i=1}^{m}a_{i}K_{4,d_{x},d_{i}}(x,p_{i})\,.$

From the above, we can see from the mutually almost orthogonality property that the index function $I_{\text{mo}}(x)$ has a large magnitude if $x$ is close to one of the points $\{q_{j}\}_{j=1}^{m}$ inside the potential inclusions, i.e., supp$(V-V_{0})$, and is small otherwise. Meanwhile, the index function $I_{\text{di}}(x,d_{x})$ has a large magnitude if $x$ is close to one of the points $\{p_{i}\}_{i=1}^{n}$ inside the conductivity inclusions, i.e., supp$(\sigma-\sigma_{0})$, as well as $d_{x}\approx d_{i}$ for such $i$, and is small otherwise. Therefore, this decouples the effect of Green’s functions and their gradients in the near field or scattered data with the help of monopole and dipole probing functions, thanks to the mutually almost orthogonality property. In order to maximize such a decoupling effect, different types of boundary influxes and probing directions are also analysed. The above properties and strategies for decoupling will be addressed in further detail in the rest of the work.

Under the settings above, two index functions in (2.1) give rise to our new Direct Sampling Method:

Given the measurement data $u-u_{0}$ on $\partial\Omega$, and a set of discrete sampling points $x\in\Omega$,

(i) evaluate $I_{\text{mo}}$ to recover the potential inclusions, i.e., supp$(V-V_{0})$;

(ii) evaluate $I_{\text{di}}$ to recover the conductivity inclusions, i.e., supp$(\sigma-\sigma_{0})$.

In this section, we aim to verify the fundamental property introduced in section 2 for some typical cases that we encounter in real applications. In particular, we intend to derive an approximation of the measurements as a combination of the Green’s functions of the homogeneous medium and their gradients when $\sigma$ is either smooth or piecewise constant.

Associated with the model (1.1), the incident field $u_{0}$ from the homogeneous background satisfies

$$\begin{cases}-\nabla\cdot(\sigma_{0}\nabla u_{0})+V_{0}u_{0}=0\quad&\text{in}\quad\Omega\,,\\ \frac{\partial u_{0}}{\partial\nu}=f\quad&\text{on}\quad\partial\Omega\,.\end{cases}$$

(3.1)

Combining the systems (1.1) and (3.1), we readily see

$$\begin{cases}-\Delta(u-u_{0})+\frac{V_{0}}{\sigma_{0}}(u-u_{0})=\frac{1}{\sigma_{0}}[\nabla\cdot((\sigma-\sigma_{0})\nabla u)-(V-V_{0})u]\quad&\text{in}\quad\Omega\,,\\ \frac{\partial(u-u_{0})}{\partial\nu}=0&\text{on}\quad\partial\Omega\,.\end{cases}$$

(3.2)

If $V_{0}\neq 0$, we consider the Green’s function $G_{x}$ for $x\in\Omega$ satisfying

$$-\Delta G_{x}+\frac{V_{0}}{\sigma_{0}}G_{x}=\delta_{x}\quad\text{in}\quad\Omega\,,\qquad\frac{\partial G_{x}}{\partial\nu}=0\quad\text{on}\quad\partial\Omega\,.$$

(3.3)

Then the difference $u-u_{0}$ can be represented by

$$(u-u_{0})(x)=\frac{1}{\sigma_{0}}\int_{\Omega}\bigg{[}\nabla_{y}\cdot((\sigma-\sigma_{0})\nabla_{y}u)-(V-V_{0})u\bigg{]}G_{x}dy\,.$$

(3.4)

On the other hand, if $V_{0}=0$, we consider the following Green’s function $G_{x}$ for $x\in\Omega$ instead:

$$-\Delta G_{x}=\delta_{x}\quad\text{in}\quad\Omega\,,\qquad\frac{\partial G_{x}}{\partial\nu}=-\frac{1}{|\partial\Omega|}\quad\text{on}\quad\partial\Omega\,,\qquad\int_{\partial\Omega}G_{x}ds=0\,.$$

(3.5)

Then we can obtain a similar representation to (3.4).

From now on, we shall consider only the following two typical cases: either $\sigma$ is smooth or piecewise constant. First for the case when $\sigma\in C^{1}(\Omega)$, by writing $D:=\text{supp}(\sigma-\sigma_{0})\bigcup\text{supp}(V-V_{0})\Subset\Omega$, we can readily derive from (3.4) by the divergence theorem that

$$\begin{split}(u-u_{0})(x)=&\,\frac{1}{\sigma_{0}}\bigg{[}\int_{\partial\Omega}(\sigma-\sigma_{0})G_{x}\frac{\partial u}{\partial\nu}ds(y)-\int_{\Omega}(\sigma-\sigma_{0})\nabla_{y}u\cdot\nabla_{y}G_{x}dy-\int_{\Omega}(V-V_{0})uG_{x}dy\bigg{]}\\ =&-\frac{1}{\sigma_{0}}\bigg{[}\int_{\Omega}(\sigma-\sigma_{0})\nabla_{y}u\cdot\nabla_{y}G_{x}dy+\int_{\Omega}(V-V_{0})uG_{x}dy\bigg{]}\,.\end{split}$$

(3.6)

Next, we consider the case when $\sigma$ is piecewise constant. We assume that $D=\overline{\cup_{i=1}^{m}\Omega_{i}}$, where $\Omega_{i}$ are open subsets of $\Omega$ with smooth boundary such that $\Omega_{i}\bigcap\Omega_{j}=\emptyset$, and that $\sigma=\sigma_{i}$ in $\Omega_{i}$ for some constant $\sigma_{i}$. And we further write $\Omega_{0}=\bar{\Omega}\setminus D$ for simplicity. Then for all $\phi\in H^{1}({\Omega})$, we derive from (1.1) that

$$\begin{split}0&=\sum_{i=0}^{m}\bigg{(}\int_{\Omega_{i}}\sigma_{i}\nabla u\cdot\nabla\phi dy\bigg{)}-\int_{\partial\Omega}\sigma_{0}f\phi ds(y)+\int_{\Omega}Vu\phi dy\\ &=\sum_{i=0}^{m}\bigg{[}\int_{\Omega_{i}}\bigg{(}-\sigma_{i}\Delta u+Vu\bigg{)}\phi dy\bigg{]}+\sum_{i=1}^{m}\bigg{[}\int_{\partial\Omega_{i}}\bigg{(}\sigma_{i}\frac{\partial u^{-}}{\partial\nu}-\sigma_{0}\frac{\partial u^{+}}{\partial\nu}\bigg{)}\phi ds(y)\bigg{]}\,.\\ \end{split}$$

(3.7)

Noticing that the normal derivative of $u$ has a jump across $\partial\Omega_{i}$, we get for $v:=\sigma u$ from (3.7) that

$$-\Delta v+\frac{V_{0}}{\sigma_{0}}v=(\frac{\sigma}{\sigma_{0}}V_{0}-V)u\quad\text{in}\quad\Omega\setminus(\cup_{i=1}^{m}\partial\Omega_{i})\,;\quad\frac{\partial v^{+}}{\partial\nu}|_{\partial\Omega_{i}}=\frac{\partial v^{-}}{\partial\nu}|_{\partial\Omega_{i}}~{}~{}\text{on}~{}~{}\partial\Omega_{i}\,,\quad\frac{\partial v}{\partial\nu}=\frac{f}{\sigma_{0}}~{}~{}\text{on}~{}~{}\partial\Omega\,,$$

where we have chosen the normal vector to point towards $\Omega_{0}$ on each $\partial\Omega_{i}$, and will write the jump of any function $w$ across the boundary $\partial\Omega_{i}$ as $[w]:=w^{+}-w^{-}$. The above equation readily implies the equation for $\gamma:=\sigma u-\sigma_{0}u_{0}$

$$\begin{cases}-\Delta\gamma+\frac{V_{0}}{\sigma_{0}}\gamma=-(V-V_{0})u+(\sigma-\sigma_{0})\frac{V_{0}}{\sigma_{0}}u&\text{in}\quad\Omega\setminus(\cup_{i=1}^{m}\partial\Omega_{i})\,,\\ \frac{\partial\gamma}{\partial\nu}=0\quad\text{on}\quad\partial\Omega\,,\qquad\frac{\partial\gamma^{+}}{\partial\nu}|_{\partial\Omega_{i}}=\frac{\partial\gamma^{-}}{\partial\nu}|_{\partial\Omega_{i}}&\text{on}\quad\partial\Omega_{i}\,.\end{cases}$$

(3.8)

For any $x\in\Omega_{0}$, we can easily write

$$\begin{split}&\sum_{i=1}^{m}\int_{\partial\Omega_{i}}[\gamma]\frac{\partial G_{x}}{\partial\nu}ds(y)\\ =&\sum_{i=1}^{m}\int_{\partial\Omega_{i}}\bigg{(}\gamma^{+}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{+}}{\partial\nu}G_{x}\bigg{)}ds(y)-\sum_{i=1}^{m}\int_{\partial\Omega_{i}}\bigg{(}\gamma^{-}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{-}}{\partial\nu}G_{x}\bigg{)}ds(y)\\ =&\bigg{[}\sum_{i=1}^{m}\int_{\partial\Omega_{i}}\bigg{(}\gamma^{+}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{+}}{\partial\nu}G_{x}\bigg{)}ds(y)-\int_{\partial\Omega}\bigg{(}\gamma\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma}{\partial\nu}G_{x}\bigg{)}ds(y)\bigg{]}\\ &-\bigg{[}\sum_{i=1}^{m}\int_{\partial\Omega_{i}}\bigg{(}\gamma^{-}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{-}}{\partial\nu}G_{x}\bigg{)}ds(y)\bigg{]}\,.\end{split}$$

(3.9)

Applying the Green’s formula in $\Omega_{0}$ to the first part of the above difference, we obtain

$$\begin{split}&\sum_{i=1}^{m}\bigg{[}\int_{\partial\Omega_{i}}\bigg{(}\gamma^{+}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{+}}{\partial\nu}G_{x}\bigg{)}ds(y)\bigg{]}-\int_{\partial\Omega}\bigg{(}\gamma\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma}{\partial\nu}G_{x}\bigg{)}ds(y)\\ =&\int_{\Omega_{0}}\bigg{(}\Delta\gamma G_{x}-\gamma\Delta G_{x}\bigg{)}dy=\int_{\Omega_{0}}\bigg{[}(V-V_{0})uG_{x}+\sigma_{0}(u-u_{0})\bigg{]}dy\,.\end{split}$$

(3.10)

Meanwhile, for the second part of the difference in (LABEL:pc1), we notice the following for each $\Omega_{i}$:

$$\begin{split}&\int_{\partial\Omega_{i}}\bigg{(}\gamma^{-}\frac{\partial G_{x}}{\partial\nu}-\frac{\partial\gamma^{-}}{\partial\nu}G_{x}\bigg{)}ds(y)=\int_{\Omega_{i}}\bigg{[}-(V-V_{0})uG_{x}+(\sigma_{i}-\sigma_{0})u\Delta G_{x}\bigg{]}dy\\ =&-\int_{\Omega_{i}}\bigg{[}(V-V_{0})uG_{x}+(\sigma_{i}-\sigma_{0})\nabla u\cdot\nabla G_{x}\bigg{]}dy+\int_{\partial\Omega_{i}}(\sigma_{i}-\sigma_{0})u^{-}\frac{\partial G_{x}}{\partial\nu}ds(y)\,.\end{split}$$

(3.11)

Combining (LABEL:pc1)-(LABEL:pc3), we come to the difference of the potentials

		$\displaystyle(u-u_{0})(x)$		(3.12)
	$\displaystyle=$	$\displaystyle-\frac{1}{\sigma_{0}}\bigg{\{}\int_{\Omega}(V-V_{0})uG_{x}dy+\sum_{i=1}^{m}\bigg{[}\int_{\Omega_{i}}(\sigma_{i}-\sigma_{0})\nabla u\cdot\nabla G_{x}dy+\int_{\partial\Omega_{i}}([\gamma]-(\sigma_{i}-\sigma_{0})u^{-})\frac{\partial G_{x}}{\partial\nu}ds(y)\bigg{]}\bigg{\}}\,.$

Now using some appropriate numerical quadrature rule, we can easily see from the expressions (3.6) and (3.12) that the boundary data or the scattered field can be approximated by

$$(u-u_{0})(x)\approx\sum_{j=1}^{n}c_{j}G_{q_{j}}(x)+\sum_{i=1}^{m}a_{i}d_{i}\cdot\nabla G_{p_{i}}(x)\,,\quad x\in\partial\Omega$$

(3.13)

for some coefficients $a_{i}\in\mathbb{C}$, $c_{j}\in\mathbb{C}$, $d_{i}\in\mathbb{S}^{d-1}$, and some quadrature points $p_{i}\in\text{supp}(\sigma-\sigma_{0})$ and $q_{j}\in\text{supp}(V-V_{0})$. We have therefore verified the fundamental property introduced in section 2.

In this section, we aim at obtaining some explicit expressions of our choices of probing functions in some special domains for more efficient numerical computation. With the same technique, We can also obtain explicit expressions of kernels $K_{i}$ introduced in (4.10) and (4.11) in those cases, which help us understand more precisely the behaviour of those kernels, and verify the mutually almost orthogonality properties.

For the notational sake, we shall write $k^{2}:={V_{0}}/{\sigma_{0}}$ from now on. The Poincare-Steklov operator plays an essential role in our subsequent analysis. We define the Neumann-to-Dirichlet map (NtD) as $\Lambda f=g$, where $f$ and $g$ satisfy the equations

$$\begin{cases}-\Delta\Phi+k^{2}\Phi=0\quad&\text{in}\quad\Omega\,,\\ \frac{\partial}{\partial\nu}\Phi={f}\quad&\text{on}\quad\partial\Omega\,,\\ \Phi=g\quad&\text{on}\quad\partial\Omega\,.\end{cases}$$

(5.1)

We recall that $\Lambda:H^{-\frac{1}{2}}(\partial\Omega)\rightarrow H^{\frac{1}{2}}(\partial\Omega)$ is a compact self-adjoint operator when we restrict ourselves to $L^{2}(\partial\Omega)$. Therefore, there exists a complete orthonormal basis consists of eigenfunctions of $\Lambda$. We notice that, in some special cases, this set of eigenfunctions coincides with the set of eigenfunctions of the surface Laplacian $\Delta_{\partial\Omega}$. This helps us to write both probing functions and the $H^{\gamma}(\partial\Omega)$ semi-inner product defined in (4.7) explicitly via Fourier coefficients with respect to the same orthonormal basis. In this section, we will focus on one such case, that is when $\partial\Omega=R\mathbb{S}^{d-1}$ for some $R>0$ and $d\geq 2$, which is a typical geometric shape used in many applications.

We like to point out that, although the two sets of eigenfunctions differ in general, they are comparable to each other based on the following observation: if we denote $\|\xi\|^{2}_{g(x)}:=\langle\xi,\,g^{-1}(x)\xi\rangle$, the dual norm of $\xi$ under the metric $g(x)$ on the surface, then the principle symbol of $\Delta_{\partial\Omega}$ is $\|\xi\|_{g(x)}^{2}$, while that of $\Lambda$ is $\|\xi\|_{g(x)}^{-1}$ (Proposition 8.53, [38]). With this, via an application of the generalized Weyl’s law, we can obtain a precise comparison of the pointwise asymptotic average squared density between the two sets of eigenfunctions. In fact, one readily checks that the volume of the variety coming from the two Hamiltonians $\{\xi:\|\xi\|_{g(x)}^{2}=1\}$ and $\{\xi:\|\xi\|_{g(x)}^{-1}=1\}$ are in fact the same, and the generalized Weyl’s law will therefore render us that the two sets of eigenfunctions have the same pointwise asymptotic average squared density in some sense mathematically. We skip the details of this argument for the sake of exposition, and focus only on the case $\partial\Omega=R\mathbb{S}^{d-1}$ for some $R>0$, when the two sets of eigenfunctions coincide.