# Mathematical framework for multi-frequency identification of thin insulating
and small conductive inhomogeneities

Habib Ammari1, Jin Keun Seo2 and Tingting Zhang2 1 Department of Mathematics,
ETH Zürich, Rämistrasse 101, 8092 Zürich, Switzerland.  2 Department of
Computational Science and Engineering, Yonsei University, Seoul 120-749 Korea
habib.ammari<EMAIL_ADDRESS>
###### Abstract

We are aiming to identify the thin insulating inhomogeneities and small
conductive inhomogeneities inside an electrically conducting medium by using
multi-frequency electrical impedance tomography (mfEIT). The thin insulating
inhomogeneities are considered in the form of tubular neighborhood of a curve
and small conductive inhomogeneities are regarded as circular disks. Taking
advantage of the frequency dependent behavior of insulating objects, we give a
rigorous derivation of the potential along thin insulating objects at various
frequencies. Asymptotic formula is given to analyze relationship between
inhomogeneities and boundary potential at different frequencies. In numerical
simulations, spectroscopic images are provided to visualize the reconstructed
admittivity at various frequencies. For the view of both kinds of
inhomogeneities, an integrated reconstructed image based on principle
component analysis (PCA) is provided. Phantom experiments are performed by
using Swisstom EIT-Pioneer Set.

##  1 Introduction

Multi-frequency electrical impedance tomography (mfEIT) is a noninvasive
method to provide the images of the conductivity and permittivity at various
frequencies ranging from 1 kHz to 1 MH inside an electrically conducting
object [12, 13, 15, 25, 29]. Recently, multi-frequency imaging techniques have
been paid considerable attention due to their advantages in probing thin
insulating inhomogeneities with their widths. Their potential applications
include biomedical imaging to probe biological tissues comprising insulating
cell membranes [10, 11, 14, 18, 23, 26] and non-destructive testing (NDT) to
probe insulating cracks [17, 19, 20, 22, 30].

In mfEIT, we inject ac currents using surface electrodes to produce the time-
harmonic electrical field inside the imaging object. The induced electrical
field is determined by the distribution of effective conductivity and
permittivity, the geometry of the object, electrode positions, the applied
frequency, and others [26]. In the presence of thin insulating inhomogeneities
or insulating membranes, the time-harmonic electrical potential near these
inhomogeneities varies significantly with the applied frequency. In
particular, the potential jump across the thin insulating inhomogeneities
changes a lot with frequency, because currents pass through a thin insulating
inhomogeneity as frequency increases. Those variations with frequency are
conveyed to the boundary voltages [22, 23, 30]. Multi-frequency EIT (mfEIT)
can take advantage of this frequency-dependent behavior to measure tissue
anomalies and thin insulating inhomogeneities with their thickness.

One challenging issue is how to link the spectroscopic information obtained
from mfEIT to the structural information of an imaging object containing thin
insulating inhomogeneities. We need to describe the role of the thin
insulating inhomogeneities by understanding the frequency-dependent interplays
between the real and imaginary parts of the complex potentials due to the
change of the refraction along the thin insulating inhomogeneities as
frequency varies.

This paper deals with this challenging issue of the spectroscopic mfEIT by
means of rigorous mathematical analysis with both numerical and experimental
validations. We start with the simplest model of a linear-shaped thin
insulating inhomogeneity with a uniform thickness, in order to provide a basis
of the rigorous expression on the frequency-dependent behavior of the
potential due to the influence of thin insulating inhomogeneities. Through a
rigorous asymptotic analysis, we provide the behavior of the time-harmonic
potential $u^{\omega}$ across the linear inhomogeneity with respect to angular
frequency $\omega$. To be precise, the frequency-dependent variation of the
exterior normal derivative $\frac{\partial u^{\omega}}{\partial\nu}$ on the
boundary of the inhomogeneity can be approximately described in terms of the
thickness ($2\delta$), the contrast ratio between the background admittivity
$\sigma_{b}+i\omega\epsilon_{b}$ and the inhomogeneity permittivity
$\epsilon_{c}$:

$\left|\frac{\partial
u^{\omega}}{\partial\nu}\right|~{}~{}\approx~{}~{}\frac{\omega}{2\delta}~{}\frac{\left|[u^{\omega}]\right|}{|(\sigma_{b}+i\omega\epsilon_{b})/\epsilon_{c}|}.$

Assuming the above approximation for a more general shape of thin insulators,
we could get better understanding on spectroscopic admittivity images being
reconstructed from standard mfEIT algorithm. In a special case of inclusions
(highly conducting disks and linear insulating segments) inside the imaging
object, we derive a formula for detecting both inhomogeneities. Finally,
various numerical simulations are performed at multiple frequencies to give
spectroscopic reconstructed images. Using spectroscopic images, an integrated
image based on principle component analysis (PCA) is presented. In addition,
we take use of Swisstom EIT-Pioneer set to carry out phantom experiments at a
frequency range from 50kHz to 250kHz.

##  2 Mathematical model

For rigorous analysis, we use the simplified two-dimensional model by
considering axially symmetric cylindrical sections under the assumption that
the out-of-plane current density is negligible in an imaging slice. We assume
a two dimensional electrically conducting domain $\Omega$ with its connected
$C^{2}$-boundary $\partial\Omega$. We denote the conductivity distribution of
the domain by $\sigma$ and the permittivity distribution by $\epsilon$. Inside
$\Omega$, there exist thin insulating objects
${\mathcal{C}}_{k},~{}k=1,2,\ldots,N_{C}$ and small conductive objects
$D_{k},~{}k=1,2,\ldots,N_{D}$. Let $D=\cup_{k=1}^{N_{D}}D_{k}$ and
${\mathcal{C}}=\cup_{k=1}^{N_{C}}{\mathcal{C}}_{k}$ denote the collections of
the conductive objects and thin insulating objects, respectively. Since the
conductivity $\sigma$ and permittivity $\epsilon$ change abruptly across the
thin objects and conductive objects, we denote

$\gamma^{\omega}(x)=\left\\{\begin{array}[]{ll}\gamma^{\omega}_{c}=\sigma_{c}+i\omega\epsilon_{c}&\quad\mbox{for}~{}x\in\mathcal{C},\\\
\gamma^{\omega}_{d}=\sigma_{d}+i\omega\epsilon_{d}&\quad\mbox{for}~{}x\in
D,\\\
\gamma^{\omega}_{b}=\sigma_{b}+i\omega\epsilon_{b}&\quad\mbox{for}~{}x\in\Omega\backslash(D\cup\mathcal{C}),\\\
\end{array}\right.$ (1)

Because the thin objects ${\mathcal{C}}_{k}$ are highly insulating, we
consider the following extreme contrast case:

$\sigma_{c}/\sigma_{b}\approx 0.$

In the frequency range below 1MHz ($\frac{\omega}{2\pi}\leq 10^{6}$), we
inject a sinusoidal current $g(x)\sin(\omega t)$ at $x\in\partial\Omega$ where
$g$ is the magnitude of the current density on $\partial\Omega$ and $g\in
H^{-1/2}_{\diamond}(\partial\Omega)$. Here
$H^{-1/2}_{\diamond}(\partial\Omega):=\\{\phi\in
H^{-1/2}(\partial\Omega)~{}:~{}\langle\phi,1\rangle=0\\}$ with
$\langle\,,\,\rangle$ being the duality pair between
$H^{-1/2}(\partial\Omega)$ and $H^{1/2}(\partial\Omega)$. The injected current
produces the time-harmonic potential $u^{\omega}$ in $\Omega$ which is
dictated by

$\left\\{\begin{array}[]{ll}\nabla\cdot\left(\gamma^{\omega}\nabla
u^{\omega}\right)=0&\quad\mbox{in}~{}\Omega,\\\ \gamma^{\omega}\frac{\partial
u^{\omega}}{\partial\nu}=g&\quad\mbox{on}~{}\partial\Omega,\end{array}\right.$
(2)

where $\gamma^{\omega}=\sigma+i\omega\epsilon$ is the admittivity
distribution, ${\nu}$ is the outward unit normal vector on $\partial\Omega$,
and $\frac{\partial}{\partial\nu}$ is the normal derivative. Setting
$\int_{\partial\Omega}u^{\omega}ds=0$, we can obtain a unique solution
$u^{\omega}$ to (2) from the Lax-Milgram theorem. Hence, we can define the
Neumann-to-Dirichlet map
$\Lambda_{\omega}:H^{-1/2}_{\diamond}(\partial\Omega)\to
H^{1/2}_{\diamond}(\partial\Omega)$ by
$\Lambda_{\omega}(g)=u^{\omega}|_{\partial\Omega}$. Using $N_{E}-$channel
multi-frequency EIT system, we may inject $N_{E}$ number of linearly
independent currents at several angular frequencies
$\omega_{1},\ldots,\omega_{N_{\omega}}$ and measure the induced corresponding
boundary voltages. We collect these current-voltage data
$\\{\Lambda_{\omega_{j}}(g_{k})~{}:~{}~{}k=1,\ldots,N_{E},~{}j=1,\ldots,N_{\omega}\\}$
at various frequencies ranging from 10Hz to 1MHz. The inverse problem is to
identify the thin insulators ${\mathcal{C}}_{k}$ and small conductive objects
$D_{k}$ from measured current-voltage data in multi-frequency EIT system.

To carry out detailed analysis, we will restrict our considerations to
geometric structures of ${\mathcal{C}}$ and $D$. We assume that each thin
insulating inhomogeneity ${\mathcal{C}}_{k}$ has a uniform thickness of
$\delta_{k}$ and is a neighborhood of a $C^{3}-$smooth open curve
$\mathcal{L}_{k}$:

$\mathcal{C}_{k}=\\{x+h\nu_{x}~{}:~{}x\in\mathcal{L}_{k},~{}-\delta_{k}<h<\delta_{k}\\},~{}~{}~{}(k=1,2,\ldots,N_{C}),$
(3)

where $\nu_{x}$ is the unit normal vector at $x$ on the curve
$\mathcal{L}_{k}$. The thickness to the length ratio is assumed to be very
small, that is, $\delta_{k}\approx 0$. We also assume that each small
conductive object has the form

$D_{k}:=z_{k}+\delta_{D}B_{k},~{}~{}~{}(k=1,2,\ldots,N_{D}),$ (4)

where $z_{k}$ is the center of $D_{k}$, $B_{k}$ is a bounded smooth reference
domain centered at $(0,0)$ and $\delta_{D}$ is related to the diameter of
$D_{k}$. We assume that ${\mathcal{C}}_{k}$ and $D_{k}$ are well separated
from each other as well as from the boundary $\partial\Omega$, i.e., the
separation distance is much larger than the characteristic size of the
conductors or thickness of the insulators. Note that in the non-resolved case,
where the distance between small conductors or insulators is of order of the
characteristic size of $D_{k}$ or thickness of ${\mathcal{C}}_{k}$, the
conductors and the insulators can not be determined separately. Only
equivalent targets can be imaged using boundary measurements (see [4]).
Therefore, throughout this paper, we assume that there exists a positive
constant $d_{0}>0$ [1, 2] such that

$\begin{array}[]{c}\inf_{k\neq
k^{\prime}}\mbox{dist}(D_{k},D_{k^{\prime}})\geq d_{0},~{}~{}\inf_{k\neq
k^{\prime}}\mbox{dist}(\mathcal{C}_{k},\mathcal{C}_{k^{\prime}})\geq d_{0},\\\
\mbox{dist}({\mathcal{C}},\partial\Omega)\geq
d_{0},~{}~{}\mbox{dist}(D,\partial\Omega)\geq
d_{0},~{}~{}\mbox{dist}_{k,j}({\mathcal{C}}_{k},D_{j})\geq d_{0}.\end{array}$
(5)

##  3 Methods

In this part, we will focus on rigorous analysis of the frequency-dependent
behaviors of the complex potential $u^{\omega}$ around the thin insulating
inhomogeneities. We will derive an explicit formula for detecting positions of
conductive inhomogeneities and thin insulating inhomogeneities by using
asymptotic expansions of $u^{\omega}$. The explicit formula depends on the
operating frequency $\omega$ and the insulator thickness $\delta_{k}$.

To start with, since each thin object $\mathcal{C}_{k}$ is highly non-
conductive, there is a noticeable potential jump along the thin insulating
objects [1, 23, 30]. For ease of notation, we define
exterior($+$)/interior($-$) normal derivative on the boundary of
$\mathcal{C}_{k}$ as follows:

$\begin{array}[]{lll}\frac{\partial
u^{\omega}}{\partial\nu}(x-\delta_{k}\nu_{x})|_{\pm}&=&\lim_{s\rightarrow
0^{+}}\frac{\partial u^{\omega}}{\partial\nu}(x-\delta_{k}\nu_{x}\mp
s\nu_{x})\\\ \frac{\partial
u^{\omega}}{\partial\nu}(x+\delta_{k}\nu_{x})|_{\pm}&=&\lim_{s\rightarrow
0^{+}}\frac{\partial u^{\omega}}{\partial\nu}(x+\delta_{k}\nu_{x}\pm
s\nu_{x})\end{array}\quad~{}\quad(x\in\mathcal{L}_{k}).$

Denote by $[u^{\omega}]_{k}$ and $\left[\frac{\partial
u^{\omega}}{\partial\nu}\right]_{k}$ the jump of the potential and the jump of
its normal derivative across the boundary of the thin insulating object
${\mathcal{C}}_{k}$, respectively:

$\begin{array}[]{lll}[u^{\omega}]_{k}(x)&:=&u^{\omega}(x+\delta_{k}\nu_{x})-u^{\omega}(x-\delta_{k}\nu_{x}),\\\
\left[\frac{\partial u^{\omega}}{\partial\nu}\right]_{k}(x)&:=&\frac{\partial
u^{\omega}}{\partial\nu}(x+\delta_{k}\nu_{x})|_{+}-\frac{\partial
u^{\omega}}{\partial\nu}(x-\delta_{k}\nu_{x})|_{+},\end{array}\quad~{}\quad(x\in\mathcal{L}_{k}).$
(6)

Then, based on the above notations, we have the following theorem for jump
conditions across thin insulating inhomogeneities:

######  Theorem 3.1.

Let ${\mathcal{C}}_{k}=\\{x:|x_{1}|\leq\ell_{k},|x_{2}|\leq\delta_{k}\\}$ be a
tubular neighbourhood of segment
${\mathcal{L}}_{k}=\\{x:|x_{1}|\leq\ell_{k},x_{2}=0\\}$ as shown in Figure 1.
Let
$\mathcal{L}_{k}^{\diamond}=\\{x\in\mathcal{L}_{k}:~{}dist(x,\partial\mathcal{L}_{k})>c_{0}\\}$
for $c_{0}>0$ and
${\mathcal{C}}_{k}^{\diamond}=\\{x+h\nu_{x}:~{}x\in\mathcal{L}_{k}^{\diamond},~{}-\delta_{k}<h<\delta_{k}\\},~{}(k=1,2,\cdots,N_{C})$.
For $\delta_{k}\approx 0$, the jump of potential and the jump of normal
derivative across ${\mathcal{C}}_{k}^{\diamond}$ can be approximated by

$\begin{array}[]{lll}[u^{\omega}]_{k}(x)&=&\frac{2\delta_{k}}{\lambda_{c}(\omega)}\frac{\partial
u^{\omega}}{\partial\nu}(x-\delta_{k}\nu_{x})|_{+}+O(\delta_{k}^{2}\ln\delta_{k})\\\
\left[\frac{\partial
u^{\omega}}{\partial\nu}\right]_{k}(x)&=&O(\delta_{k}\ln\delta_{k})\end{array}\quad~{}\quad(x\in\mathcal{L}_{k}^{\diamond}),$
(7)

with

$\lambda_{c}(\omega)=\frac{\sigma_{c}+i\omega\epsilon_{c}}{\sigma_{b}+i\omega\epsilon_{b}}.$
(8)

Figure 1: Linear object model.

###### Proof.

The potential $u^{\omega}$ in (2) can be expressed as

$\displaystyle u^{\omega}(x)$ $\displaystyle=$ $\displaystyle
H_{k}(x)+\int_{\partial{\mathcal{C}}_{k}}\Gamma(x,x^{\prime})\phi_{k}(x^{\prime})ds_{x^{\prime}}\quad\mbox{in
}~{}\Omega,$

where $\Gamma(x,x^{\prime}):=-\frac{1}{2\pi}\ln|x-x^{\prime}|$ is the
fundamental solution of Laplacian in two dimensions, $H_{k}$ is a harmonic
function in a neighborhood of ${\mathcal{C}}_{k}$, and

$\phi_{k}=\left(\lambda_{c}-1\right)\frac{\partial
u^{\omega}}{\partial\nu}\bigg{|}_{-}=\left(1-\frac{1}{\lambda_{c}}\right)\frac{\partial
u^{\omega}}{\partial\nu}\bigg{|}_{+}\quad\mbox{on
}~{}\partial{\mathcal{C}}_{k}.$ (9)

To be precise, $H_{k}$ can be expressed in the following form:

$H_{k}=-\frac{1}{\gamma_{b}^{\omega}}{\mathcal{S}}_{\Omega}g+{\mathcal{D}}_{\Omega}(\Lambda_{\omega}(g))+\sum_{j=1}^{N_{D}}{\mathcal{S}}_{D_{j}}\psi_{D_{j}}+\sum_{j\neq
k}{\mathcal{S}}_{{\mathcal{C}}_{j}}\phi_{j},$

where
$\psi_{D_{j}}=\left(\frac{\gamma_{d}^{\omega}}{\gamma_{b}^{\omega}}-1\right)\frac{\partial
u^{\omega}}{\partial\nu}\bigg{|}_{-}$ on ${\partial D_{j}}$,
${\mathcal{S}}_{\Omega}$, ${\mathcal{S}}_{D_{j}}$,
${\mathcal{S}}_{{\mathcal{C}}_{j}}$ are the single layer potentials over the
domain $\Omega,~{}D_{j},~{}{\mathcal{C}}_{j}$ respectively and
${\mathcal{D}}_{\Omega}$ is the double layer potential over the domain
$\Omega$; see [2, 3, 26].

From the transmission conditions across ${\mathcal{C}}_{k}$, $\phi_{k}$
satisfies

$\left(\lambda
I-{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\right)\phi_{k}=\frac{\partial
H_{k}}{\partial\nu}\quad\mbox{on }~{}\partial{\mathcal{C}}_{k},$ (10)

where $\lambda=\frac{\lambda_{c}+1}{2(\lambda_{c}-1)}$ and
${\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi$ is given by

$\mathcal{K}^{*}_{{\mathcal{C}}_{k}}\phi(x):=\frac{1}{2\pi}\int_{\partial{{\mathcal{C}}_{k}}}\frac{\langle
x^{\prime}-x,\nu_{x}\rangle}{|x-x^{\prime}|^{2}}\phi(x^{\prime})ds_{x^{\prime}},\quad
x\in\partial{{\mathcal{C}}_{k}}.$ (11)

Now, we compute $\left[\frac{\partial u^{\omega}}{\partial\nu}\right]_{k}(x)$
for $x=(s,0)\in{\mathcal{L}}_{k}^{\diamond}$. From (10), we have

$\begin{array}[]{lll}&&\left(\lambda
I-{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\right)\phi_{k}(s,\delta_{k})-\left(\lambda
I-{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\right)\phi_{k}(s,-\delta_{k})\\\
&&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=\int_{-\delta_{k}}^{\delta_{k}}\frac{d^{2}H_{k}(s,t)}{dt^{2}}dt\quad\mbox{for
}~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond}.\end{array}$ (12)

For convenience, let
$\partial{{\mathcal{C}}_{k}}=\partial{{\mathcal{C}}_{k}}^{\sharp}\cup{\mathcal{L}}_{k}^{+}\cup{\mathcal{L}}_{k}^{-}$
where
$\partial{{\mathcal{C}}_{k}}^{\sharp}=\\{x\in\partial{{\mathcal{C}}_{k}}:|x_{1}|>\ell_{k}-c_{0}/2\\}$
and
${\mathcal{L}}_{k}^{\pm}=\\{x\in\partial{{\mathcal{C}}_{k}}:|x_{1}|\leq\ell_{k}-c_{0}/2,x_{2}=\pm\delta_{k}\\}$
as described in Figure 2.

Figure 2: $\partial{{\mathcal{C}}_{k}}$ can be divided as three parts; green
color stands for $\partial{{\mathcal{C}}_{k}}^{\sharp}$, red color is
${\mathcal{L}}_{k}^{+}$ and blue color is ${\mathcal{L}}_{k}^{-}$.

Then, when $x=(s,-\delta_{k})$,
${\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(x)$ can be written as

$\displaystyle
2\pi{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(s,-\delta_{k})$
$\displaystyle=$
$\displaystyle\int_{\partial{{\mathcal{C}}_{k}}^{\sharp}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}$
$\displaystyle+\int_{{\mathcal{L}}_{k}^{+}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}+\int_{{\mathcal{L}}_{k}^{-}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}.$

Since ${\mathcal{L}}_{k}$ is a line, $x^{\prime}-(s,-\delta_{k})$ is
perpendicular to $\nu_{x}$ for $x^{\prime}\in{\mathcal{L}}_{k}^{-}$ and
$\langle x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle=2\delta_{k}$ for
$x^{\prime}\in{\mathcal{L}}_{k}^{+}$, which directly yield

$\displaystyle
2\pi{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(s,-\delta_{k})$
$\displaystyle=$
$\displaystyle\int_{\partial{{\mathcal{C}}_{k}}^{\sharp}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}+\int_{{\mathcal{L}}_{k}^{+}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}$
$\displaystyle=$
$\displaystyle\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}\phi_{k}(s^{\prime})ds^{\prime}+\int_{\partial{{\mathcal{C}}_{k}}^{\sharp}}\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}.$

Analogously, it follows that

$\displaystyle
2\pi{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(s,\delta_{k})$
$\displaystyle=$
$\displaystyle\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}\phi_{k}(s^{\prime})ds^{\prime}+\int_{\partial{{\mathcal{C}}_{k}}^{\sharp}}\frac{\langle
x^{\prime}-(s,\delta_{k}),\nu_{x}\rangle}{|(s,\delta_{k})-x^{\prime}|^{2}}\phi_{k}(x^{\prime})dx^{\prime}.$

since $x^{\prime}-(s,\delta_{k})$ is perpendicular to $\nu_{x}$ for
$x^{\prime}\in{\mathcal{L}}_{k}^{+}$ and $\langle
x^{\prime}-(s,\delta_{k}),\nu_{x}\rangle=2\delta_{k}$ for
$x^{\prime}\in{\mathcal{L}}_{k}^{-}$. Combining formulas (LABEL:Eq:Kminus) and
(LABEL:Eq:Kplus) yields

${\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(s,\delta_{k})-{\mathcal{K}}_{{\mathcal{C}}_{k}}^{*}\phi_{k}(s,-\delta_{k})=\mathcal{T}_{k}\phi^{\sharp}_{k}(s)+E(s)\quad\mbox{for
}~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond},$ (15)

where

$\displaystyle\mathcal{T}_{k}\phi^{\sharp}_{k}(s)=\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}\phi_{k}^{\sharp}(s^{\prime})ds^{\prime},$
(16)
$\displaystyle\phi_{k}^{\sharp}(s)=\phi_{k}(s,\delta_{k})-\phi_{k}(s,-\delta_{k})=\left(1-\frac{1}{\lambda_{c}}\right)\left[\frac{\partial
u^{\omega}}{\partial\nu}\right]_{k}(s,0),$ (17)

and $E(s)$ for $(s,0)\in{\mathcal{L}}_{k}^{\diamond}$ is given by

$E(s)=\frac{1}{2\pi}\int_{\partial{{\mathcal{C}}_{k}}^{\sharp}}\left\\{\frac{\langle
x^{\prime}-(s,\delta_{k}),\nu_{x}\rangle}{|(s,\delta_{k})-x^{\prime}|^{2}}-\frac{\langle
x^{\prime}-(s,-\delta_{k}),\nu_{x}\rangle}{|(s,-\delta_{k})-x^{\prime}|^{2}}\right\\}\phi_{k}(x^{\prime})ds_{x^{\prime}}.$

Since $H(x)$ is harmonic in a neighborhood of ${\mathcal{C}}_{k}$,

$\sup_{(s,0)\in{\mathcal{L}}_{k}^{\diamond}}\left|\int_{-\delta_{k}}^{\delta_{k}}\frac{d^{2}H_{k}(s,t)}{dt^{2}}dt\right|\leq
2\delta_{k}\|\nabla\nabla H_{k}\ \|_{C^{2}({\mathcal{C}}_{k})}=O(\delta_{k}).$
(18)

Moreover, the term $E(s)$ can be estimated by the mean-value theorem. Using
the fact that
$c_{0}/2\leq\mbox{dist}({\mathcal{L}}_{k}^{\diamond},\partial{{\mathcal{C}}_{k}}^{\sharp})$
and $|\langle x^{\prime}-(s,\pm\delta_{k}),\nu_{x}\rangle|\leq 2\delta_{k}$,
we get

$\sup_{(s,0)\in{\mathcal{L}}_{k}^{\diamond}}\left|E(s)\right|\leq\frac{16}{c_{0}^{2}}\|\phi_{k}\|_{L^{1}({\mathcal{C}}^{\sharp}_{k})}\delta_{k}=O(\delta_{k}).$
(19)

Therefore, it follows from (12)-(19) that

$(\lambda
I-\mathcal{T}_{k})\phi_{k}^{\sharp}(s)=\int_{-\delta_{k}}^{\delta_{k}}\frac{d^{2}H_{k}(s,t)}{dt^{2}}dt+E(s)\quad\mbox{for
}~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond}.$ (20)

Note that

$\displaystyle\mathcal{T}_{k}\phi_{k}^{\sharp}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}(\phi_{k}^{\sharp}(s^{\prime})-\phi_{k}^{\sharp}(s))ds^{\prime},$
(21)
$\displaystyle+\phi_{k}^{\sharp}(s)\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}ds^{\prime}.$

Direct computation yields

$\displaystyle\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}ds^{\prime},$
$\displaystyle=\frac{1}{2\pi}\left(\arctan\frac{\ell_{k}-c_{0}/2-s}{2\delta_{k}}-\arctan\frac{-\ell_{k}+c_{0}/2-s}{2\delta_{k}}\right),$
$\displaystyle=\frac{1}{2}+\left(-\frac{2}{\ell_{k}-c_{0}/2-s}+\frac{2}{-\ell_{k}+c_{0}/2-s}\right)\delta_{k}+O(\delta_{k}^{3}),$
$\displaystyle=\frac{1}{2}-\frac{4(\ell_{k}-c_{0}/2)}{(\ell_{k}-c_{0}/2)^{2}-s^{2}}\delta_{k}+O(\delta_{k}^{3}).$

For $s\in{\mathcal{L}}^{\diamond}_{k}$, we have

$\displaystyle(\ell_{k}-c_{0}/2)^{2}-s^{2}\geq(\ell_{k}-c_{0}/2)^{2}-(\ell_{k}-c_{0})^{2}=c_{0}(\ell_{k}-3c_{0}/4).$

Therefore, we obtain that

$\displaystyle\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}ds^{\prime}=\frac{1}{2}+O(\delta_{k}).$
(22)

Since $\frac{\partial\gamma^{\omega}}{\partial x_{1}}=0$ in a neighborhood of
${\mathcal{C}}^{\diamond}_{k}$, $\frac{\partial u^{\omega}}{\partial x_{1}}$
is a weak solution of $\nabla\cdot(\gamma^{\omega}\frac{\partial
u^{\omega}}{\partial x_{1}})=0$ near the region
${\mathcal{C}}^{\diamond}_{k}$. Hence, $\phi_{k}^{\sharp}(s)$ is
differentiable for $(s,0)\in{\mathcal{L}}_{k}^{\diamond}$, and
$\phi_{k}^{\sharp}(s)-\phi_{k}^{\sharp}(s^{\prime})$ can be estimated by

$|\phi_{k}^{\sharp}(s)-\phi_{k}^{\sharp}(s^{\prime})|\leq\left\|\frac{d}{ds}\phi_{k}^{\sharp}(s)\right\|_{L^{\infty}({\mathcal{L}}_{k}^{\diamond})}|s-s^{\prime}|\quad\mbox{for
}~{}(s^{\prime},0),~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond}.$ (23)

From the above estimate (23), we have

$\displaystyle\left|\frac{1}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}(\phi_{k}^{\sharp}(s^{\prime})-\phi_{k}^{\sharp}(s))ds^{\prime}\right|$
(24)
$\displaystyle\quad\leq\frac{\left\|\frac{d}{ds}\phi_{k}^{\sharp}(s)\right\|_{L^{\infty}({\mathcal{L}}_{k}^{\diamond})}}{2\pi}\int_{-\ell_{k}+c_{0}/2}^{\ell_{k}-c_{0}/2}\frac{2\delta_{k}}{|s-s^{\prime}|^{2}+(2\delta_{k})^{2}}|s-s^{\prime}|ds^{\prime}=-\left\|\frac{d}{ds}\phi_{k}^{\sharp}(s)\right\|_{L^{\infty}({\mathcal{L}}_{k}^{\diamond})}\frac{2\delta_{k}}{\pi}\ln
2\delta_{k}$
$\displaystyle\quad\quad+\frac{\delta_{k}\left\|\frac{d}{ds}\phi_{k}^{\sharp}(s)\right\|_{L^{\infty}({\mathcal{L}}_{k}^{\diamond})}}{2\pi}\ln\left[((\ell_{k}-\frac{c_{0}}{2}-s)^{2}+4\delta_{k}^{2})((\ell_{k}-\frac{c_{0}}{2}+s)^{2}+4\delta_{k}^{2})\right]$
$\displaystyle\quad=O(\delta_{k}\ln\delta_{k}).$

From (22) and (24), $\mathcal{T}_{k}\phi_{k}^{\sharp}(s)$ in (21) can be
estimated by

$\mathcal{T}_{k}\phi_{k}^{\sharp}(s)=\frac{1}{2}\phi_{k}^{\sharp}(s)+O(\delta_{k}\ln\delta_{k}).$
(25)

Combining (17), (20) and (25), we arrive at

$\displaystyle\left(\lambda-\frac{1}{2}\right)\phi_{k}^{\sharp}(s)=O(\delta_{k}\ln\delta_{k}),$
$\displaystyle\left[\frac{\partial
u^{\omega}}{\partial\nu}\right]_{k}(s,0)=\frac{\lambda_{c}}{\lambda_{c}-1}\phi_{k}^{\sharp}(s)=O(\delta_{k}\ln\delta_{k})\quad\mbox{for
}~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond}.$ (26)

It easily follows from the above analysis and transmission condition that

$\displaystyle[u^{\omega}]_{k}(s)$ $\displaystyle=$ $\displaystyle
u^{\omega}(s,\delta_{k})-u^{\omega}(s,-\delta_{k})=\int_{-\delta_{k}}^{\delta_{k}}\frac{d}{dt}u^{\omega}(s,t)dt$
$\displaystyle=$ $\displaystyle 2\delta_{k}\frac{\partial u^{\omega}}{\partial
x_{2}}(s,-\delta_{k})|_{-}+\int_{-\delta_{k}}^{\delta_{k}}\frac{d}{dt}u^{\omega}(s,t)-\frac{\partial
u^{\omega}}{\partial x_{2}}(s,-\delta_{k})|_{-}dt$

where

$\displaystyle|\int_{-\delta_{k}}^{\delta_{k}}\frac{d}{dt}u^{\omega}(s,t)-\frac{\partial
u^{\omega}}{\partial x_{2}}(s,-\delta_{k})|_{-}dt|$ $\displaystyle\leq$
$\displaystyle\int_{-\delta_{k}}^{\delta_{k}}\left|\frac{d}{dt}u^{\omega}(s,t)-\frac{\partial
u^{\omega}}{\partial x_{2}}(s,-\delta_{k})|_{-}\right|dt\leq
C\delta_{k}^{2}\ln\delta_{k}.$

Then we have

$[u^{\omega}]_{k}(s)=\frac{2\delta_{k}}{\lambda_{c}(\omega)}\frac{\partial
u^{\omega}}{\partial\nu}(s,-\delta_{k})|_{+}+O(\delta_{k}^{2}\ln\delta_{k})\quad\mbox{for
}~{}(s,0)\in{\mathcal{L}}_{k}^{\diamond}.$ (27)

∎

######  Remark 3.2.

The approximation formulas given in (7) can be proved for special cases of
thin insulators by using layer potential techniques. Unfortunately, the proof
of (7) in general still remains a challenging issue due to technical
difficulties in obtaining a uniform estimate for the Hessian matrix of
$u^{\omega}$ in $\mathcal{C}_{k}^{\diamond}$ with respect to $\delta_{k}$. For
the numerical proof of the jump conditions in (7), we refer to [30].

###  3.1 Effective zero-thickness model

It is very difficult to numerically solve $u^{\omega}$ in (2) due to the
potential discontinuity across the thin insulating inhomogeneities. One way to
deal with thin objects is to treat them as lower dimensional interfaces. In
our case, the thin insulating objects can be considered as curves. Based on
the above Theorem 3.1, we can describe an effective zero-thickness model by
imposing the jump conditions of $[u^{\omega}]_{k}$ and $\left[\frac{\partial
u^{\omega}}{\partial\nu}\right]_{k}$ on the curves $\mathcal{L}_{k}$. This
means that the potential $u^{\omega}$ can be approximated by the corresponding
potential $\widetilde{u}^{\omega}$ satisfying the following effective zero-
thickness model:

$\left\\{\begin{array}[]{l}\nabla\cdot((\gamma^{\omega}_{b}+(\gamma^{\omega}_{d}-\gamma^{\omega}_{b})\chi_{D})\nabla\widetilde{u}^{\omega})=0~{}~{}\quad\textrm{in}\quad\Omega\setminus\cup_{k=1}^{N_{C}}\mathcal{L}_{k},\\\
\left[\frac{\partial}{\partial\nu}\widetilde{u}^{\omega}\right]_{\mathcal{L}_{k}}=0,~{}~{}~{}k=1,2,\ldots,N_{C},\\\
\left[\widetilde{u}^{\omega}\right]_{\mathcal{L}_{k}}=2\delta_{k}\frac{1}{\lambda_{c}(\omega)}\frac{\partial\widetilde{u}^{\omega}}{\partial\nu}|_{+},~{}~{}~{}k=1,2,\ldots,N_{C},\\\
\gamma^{\omega}_{b}\frac{\partial\widetilde{u}^{\omega}}{\partial\nu}=g\quad\mbox{on}~{}\partial\Omega,\end{array}\right.$
(28)

where $\chi_{D}$ is the characteristic function of $D$ and

$\begin{array}[]{lll}[\widetilde{u}^{\omega}(x)]_{\mathcal{L}_{k}}&:=&\lim_{s\rightarrow
0^{+}}\left(\widetilde{u}^{\omega}(x+s\nu_{x})-\widetilde{u}^{\omega}(x-s\nu_{x})\right)\\\
\left[\frac{\partial\widetilde{u}^{\omega}}{\partial\nu}(x)\right]_{\mathcal{L}_{k}}&:=&\lim_{s\rightarrow
0^{+}}\left(\frac{\partial\widetilde{u}^{\omega}}{\partial\nu}(x+s\nu_{x})-\frac{\partial\widetilde{u}^{\omega}}{\partial\nu}(x-s\nu_{x})\right)\end{array}~{}\quad(x\in\mathcal{L}_{k}).$
(29)

Since $u^{\omega}\approx\tilde{u}^{\omega}$ in
$\\{x\in\Omega~{}:\mbox{dist}(x,\partial\Omega)<\frac{d_{0}}{2}\\}$, the
forward model (2) and the effective zero-thickness model (28) have basically
the same Neumann-to-Dirichlet data in terms of the inverse problem. From now
on, let $u^{\omega}$ denote the solution of (28) for ease of notation. From
the above zero-thickness model, the boundary condition along curve
$\mathcal{L}_{k}$ depends on thickness $\delta_{k}$ of thin insulating objects
as well as the value of $\lambda_{c}(\omega)$ which is related with injected
current frequency $\omega$. From the above Theorem 3.1, the jump of the
potential across ${\mathcal{C}}_{k}$ depends on angular frequency $\omega$ as
well as the thickness $\delta_{k}$. At high frequencies, the admittivity ratio
$|\lambda_{c}(\omega)|$ is away from 0, and there will be no potential jump
when $\delta_{k}\approx 0$. Whereas, at low frequencies,
$|\lambda_{c}(\omega)|$ is close to zero since $\sigma_{c}/\sigma_{b}\approx
0$ and potential jump happens across ${\mathcal{C}}_{k}$. For easy
understanding, changes of potential distribution near thin insulating objects
are given in Figure 3.

Figure 3: Changes of potential distribution near thin insulating objects with
frequencies: (a) 10Hz, (b) 10kHz, (c) 500kHz.

Based on the above analysis, we consider separately the following two cases
[21]:

  * •

High-frequency case: $\delta_{k}\approx 0$ and
$0<c_{1}\leq|\lambda_{c}(\omega)|$.

  * •

Low-frequency case: $|\lambda_{c}(\omega)|\approx 0$ and $\delta_{k}\approx 0$
with $|\lambda_{c}(\omega)|^{-1}\delta_{k}\approx\beta$ and $0<\beta<\infty$.

In the upcoming sections, asymptotic formulas connecting the boundary
potential and inhomogeneities are used to derive an explicit formula for
identification of thin insulating objects and small conductive inhomogeneities
inside $\Omega$.

###  3.2 High-frequency case: $\delta_{k}\approx 0$ and
$0<c_{1}\leq|\lambda_{c}(\omega)|$

In high-frequency case, we suppose that the injected current frequency
$\omega$ is not that low, so that $|\lambda_{c}(\omega)|$ is away from zero.
When the thickness $\delta_{k}$ goes to zero, the potential jump along
$\mathcal{C}_{k}$ also goes to zero according to approximation formula (7) as
well as Figure 3. Therefore, the proposed problem can be regarded as
traditional impedance boundary value problem and the influence of thin
insulating inhomogeneities on the high-frequency current-voltage data is very
weak. In this case, the following boundary voltage asymptotic expansion holds
at high-frequencies. For detailed analysis and similar proof, one may refer to
[1, 7, 8, 9] in which the proof can be immediately extended to the complex
valued equation.

######  Theorem 3.3 (Asymptotic expansion at high-frequencies).

Let $\lambda_{c}(\omega)$ and $\delta_{k}$ satisfy the conditions stated in
high-frequency case. Assume that $u^{\omega}$ is a solution of the effective
zero-thickness model (28) and $u_{0}$ is the solution of equation (2) with
$\gamma^{\omega}=\gamma^{\omega}_{b}$. Assume that all $\mathcal{L}_{k}$ are
line segments with endpoints $P_{k},Q_{k}$. For $x\in\partial\Omega$, when the
injection current frequency is high, the perturbations of voltage potential
$u^{\omega}$ due to small conductive objects $D_{k}$ and thin insulating
objects ${\mathcal{C}}_{k}$ can be expressed as

$\displaystyle\left(-\frac{1}{2}{I}+\mathcal{K}_{\Omega}\right)[u^{\omega}-u_{0}](x)$
(30) $\displaystyle=$
$\displaystyle-\sum_{k=1}^{N_{C}}\int_{\mathcal{L}_{k}}\delta_{k}A_{k}(x^{\prime},\lambda_{c}(\omega))\nabla
u_{0}(x^{\prime})\cdot\nabla\Gamma(x,x^{\prime})ds_{x^{\prime}}$
$\displaystyle-\delta_{D}^{2}\sum_{k=1}^{N_{D}}\nabla\Gamma(x,z_{k})\cdot
M(\lambda_{d}(\omega),B_{k})\nabla
u_{0}(z_{k})+O(\delta_{k}^{2})+O(\delta_{D}^{3}),$

where $A_{k}(x,\lambda_{c}(\omega))$ is a $2\times 2$ symmetric matrix whose
eigenvectors are unit normal vector $\nu_{k}(x)$ and unit tangential vector
$\tau_{k}(x)$ to ${\mathcal{L}}_{k}$ and the corresponding eigenvalues are
$2(1-\frac{1}{\lambda_{c}(\omega)})$ and
$2\left(\lambda_{c}(\omega)-1\right)$, respectively. Moreover,
$M(\lambda_{d}(\omega),B_{k})$ is polarization tensor given by

$M_{ij}:=\int_{\partial
B_{k}}y^{j}(\lambda_{d}(\omega)I-\mathcal{K}^{*}_{B_{k}})^{-1}(\nu_{x}\cdot\nabla
x^{i})(y)ds_{y},\quad i,j=1,2,$ (31)

with

$\lambda_{d}(\omega)=\frac{(\sigma_{d}+\sigma_{b})+i\omega(\epsilon_{d}+\epsilon_{b})}{2((\sigma_{d}-\sigma_{b})-i\omega(\epsilon_{d}-\epsilon_{b}))}.$
(32)

######  Remark 3.4.

In Theorem 3.3, the matrix $A_{k}$ is only related to the geometry of the thin
insulators and admittivity ratio $\lambda_{c}(\omega)$. Since $\sigma_{c}$ is
very small, the ratio $\lambda_{c}(\omega)$ changes much with respect to the
operating frequency $\omega$. On the other hand, polarization tensor
$M(\lambda_{d}(\omega),B_{k})$ is related with the geometry of small
conductors and admittivity distribution. When $D_{k}$ is a circular disk,
$M(\lambda_{d}(\omega),B_{k})$ can be explicitly written as
$M(\lambda_{d}(\omega),B_{k})=\frac{|B_{k}|}{\lambda_{d}(\omega)}I$. Since
$\sigma_{d}$ is realtively large compared with the background conductivity
$\sigma_{b}$ and $\sigma_{d}-\sigma_{b}$ is away from zero,
$\lambda_{d}(\omega)$ varies little with respect to frequency. Thus conductive
objects are insensitive to the boundary measurements.

Theorem 3.3 shows that the measured boundary data are influenced by insulating
objects and conductive objects since the first term on the right-hand side of
formula (30) is only related with thin insulators while the second term is
only related with small conductors. Depending on the magnitude of
$\omega,\delta_{k}$ and $\delta_{D}$, the dominative term on right-hand side
of formula (30) may be alternative. To see the effect of $\omega,\delta_{k}$
and $\delta_{D}$ on the measured boundary data more clearly, we need further
analysis on the expansion formula in Theorem 3.3.

To avoid any confusion, we will adopt the following notations:
$x=(x_{1},x_{2})$ denotes a point in $\mathbb{R}^{2}$ and
$\mathbf{x}=x\cdot(1,i)$ will be the corresponding point in ${\mathbb{C}}$.
Let $g=a\cdot\nu$ on $\partial\Omega$ where $a$ is a unit vector in
$\mathbb{R}^{2}$. Similarly, ${\bf a}=a\cdot(1,i)$, the center $z_{k}$ of each
small conductor $D_{k}$ can be expressed as
${\mathbf{z}_{k}}=z_{k}\cdot(1,i)$, the endpoints of line segment
${\mathcal{L}}_{k}$ can be written as
${\mathbf{P}_{k}}=P_{k}\cdot(1,i),~{}\hbox{and}~{}~{}{\mathbf{Q}_{k}}=Q_{k}\cdot(1,i).$
Then we have the following result.

######  Theorem 3.5 (Identification of thin insulators and small conductors).

Let $\lambda_{c}(\omega)$ and $\delta_{k}$ satisfy the conditions stated in
high-frequency case. Assume that all the thin inhomogeneities are line
segments and all the small conductive objects are disks. Then
$(-\frac{1}{2}I+\mathcal{K}_{\Omega})[u^{\omega}-u_{0}]$ on the boundary
$\partial\Omega$ can be expressed as

$\displaystyle\Re\left\\{(-\frac{1}{2}I+\mathcal{K}_{\Omega})[u^{\omega}-u_{0}]({\bf
x})\right\\}=\Re\\{G^{\Re}({\mathbf{x}})\\}+O(\delta_{k}^{2})+O(\delta_{D}^{3}),$
(33)
$\displaystyle\Im\left\\{(-\frac{1}{2}{I}+\mathcal{K}_{\Omega})[u^{\omega}-u_{0}]({\bf
x})\right\\}=\Re\left\\{G^{\Im}(x)\right\\}+O(\delta_{k}^{2})+O(\delta_{D}^{3}),$
(34)

where $G^{\Re}$ and $G^{\Im}$ are meromorphic functions:

$\displaystyle\frac{dG^{\Re}({\mathbf{x}})}{d{\mathbf{x}}}=\sum_{k=1}^{N_{C}}{\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})}\left(\frac{1}{{\mathbf{x}}-\mathbf{Q}_{k}}-\frac{1}{{\mathbf{x}}-\mathbf{P}_{k}}\right)-\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Re}(\omega,\delta_{D})}\frac{1}{({\mathbf{x}}-\mathbf{z}_{k})^{2}}$
(35)
$\displaystyle\frac{dG^{\Im}({\mathbf{x}})}{d{\mathbf{x}}}=\sum_{k=1}^{N_{C}}{\mathfrak{C}_{k}^{\Im}(\omega,\delta_{k})}\left(\frac{1}{{\mathbf{x}}-\mathbf{Q}_{k}}-\frac{1}{{\mathbf{x}}-\mathbf{P}_{k}}\right)-\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Im}(\omega,\delta_{D})}\frac{1}{({\mathbf{x}}-\mathbf{z}_{k})^{2}}$
(36)

and

$\displaystyle\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})$
$\displaystyle=\frac{\delta_{k}}{\pi}\left(\Re\\{(\lambda_{c}(\omega)-1)\\}a_{\tau_{k}}+i\Re\\{(1-\frac{1}{\lambda_{c}(\omega)})\\}a_{\nu_{k}}\right)$
(37) $\displaystyle\mathfrak{C}_{k}^{\Im}(\omega,\delta_{k})$
$\displaystyle=\frac{\delta_{k}}{\pi}\left(\Im\\{\lambda_{c}(\omega)-1\\}a_{\tau_{k}}+i\Im\\{1-\frac{1}{\lambda_{c}(\omega)}\\}a_{\nu_{k}}\right)$
(38) $\displaystyle{\mathfrak{D}}_{k}^{\Re}(\omega,\delta_{D})$
$\displaystyle=-\Re\left\\{\frac{|B_{k}|\delta_{D}^{2}}{2\pi\lambda_{d}(\omega)}\right\\}\mathbf{a},\quad\mathfrak{D}_{k}^{\Im}(\omega,\delta_{D})=-\Im\left\\{\frac{|B_{k}|\delta_{D}^{2}}{2\pi\lambda_{d}(\omega)}\right\\}\mathbf{a}.$
(39)

Here, $a_{\nu_{k}}=a\cdot\nu_{k},a_{\tau_{k}}=a\cdot\tau_{k}$.

###### Proof.

Since $B_{k}$ is a disk, the formula (31) gives
$M(\lambda_{d}(\omega),B_{k})=\frac{|B_{k}|}{\lambda_{d}(\omega)}I$. Hence,
the formula (30) in Theorem 3.3 can be expressed as

$\left(\frac{1}{2}{I}+\mathcal{K}_{\Omega}\right)[u^{\omega}-u_{0}](x)=\Phi(x)+O(\delta_{k}^{2})+O(\delta_{D}^{3})\quad\quad(x\in\partial\Omega),$
(40)

where $\Phi$ is

$\Phi(x)=-\sum_{k=1}^{N_{C}}\delta_{k}\int_{\mathcal{L}_{k}}(A_{k}~{}a)\cdot\nabla\Gamma(x,x^{\prime})ds_{X^{\prime}}-\frac{\delta_{D}^{2}}{2\pi}\sum_{k=1}^{N_{D}}\frac{|B_{k}|}{\lambda_{d}(\omega)}\frac{x-z_{k}}{|x-z_{k}|^{2}}\cdot
a.$ (41)

We use $a=a_{\nu_{k}}\nu_{k}+a_{\tau_{k}}\tau_{k}$ to get

$\displaystyle\Phi(x)$ $\displaystyle=$
$\displaystyle-\frac{1}{2\pi}\sum_{k=1}^{N_{C}}\delta_{k}\int_{\mathcal{L}_{k}}\left(2(1-\frac{1}{\lambda_{c}(\omega)})a_{\nu_{k}}\nu_{k}+2(\lambda_{c}(\omega)-1)a_{\tau_{k}}\tau_{k}\right)\cdot\frac{x-x^{\prime}}{|x-x^{\prime}|^{2}}ds_{x^{\prime}}$
(42)
$\displaystyle-\frac{\delta_{D}^{2}}{2\pi\lambda_{d}(\omega)}\sum_{k=1}^{N_{D}}|B_{k}|\frac{x-z_{k}}{|x-z_{k}|^{2}}\cdot
a\quad\quad\quad\quad(x\in\partial\Omega).$

From now on, we shall identify $\mathbb{R}^{2}$ with the complex plane
$\mathbb{C}$ and use similar ideas as those in [1, 6]. Since
$\lambda_{c}(\omega)$, $\lambda_{d}(\omega)$ as well as $u^{\omega}$ are
complex, we will consider real and imaginary parts of $\Phi(x)$ separately.

The real part of $\Phi(x)$ for $x\in\partial\Omega$ can be expressed as

$\displaystyle\Re\\{\Phi({\mathbf{x}})\\}=\Re\left\\{-\frac{1}{2\pi}\sum_{k=1}^{N_{C}}\delta_{k}\int_{\mathcal{L}_{k}}{\frac{\mbox{\boldmath$\xi$}}{\bf
x-\bf x^{\prime}}}ds_{\bf
x^{\prime}}-\Re\left\\{\frac{\delta_{D}^{2}}{2\pi\lambda_{d}(\omega)}\right\\}\sum_{k=1}^{N_{D}}|B_{k}|\frac{\bf
a}{{\mathbf{x}}-\mathbf{z}_{k}}\right\\},$ (43)

where $\xi$$=\xi\cdot(1,i)$ and $\xi$ is

$\displaystyle\xi=\Re\\{2(1-\frac{1}{\lambda_{c}(\omega)})a_{\nu_{k}}\nu_{k}+2(\lambda_{c}(\omega)-1)a_{\tau_{k}}\tau_{k}\\}.$
(44)

Since $\mathcal{L}_{k}$ is the segment with endpoints $P_{k},Q_{k}$, it can be
written as $P_{k}+t(Q_{k}-P_{k}),0\leq t\leq 1$. Therefore, $\mathcal{L}_{k}$
has its unit tangent vector
${\mbox{\boldmath$\tau$}_{k}}={\frac{\mathbf{Q}_{k}-\mathbf{P}_{k}}{|\mathbf{P}_{k}-\mathbf{Q}_{k}|}}$
and its unit normal vector
${\mbox{\boldmath$\nu$}_{k}}=i{\frac{\mathbf{Q}_{k}-\mathbf{P}_{k}}{|\mathbf{P}_{k}-\mathbf{Q}_{k}|}}$
in $\mathbb{C}$. Hence, the integral term in (43) can be written as

$\displaystyle\int_{\mathcal{L}_{k}}{\bf\frac{\mbox{\boldmath$\xi$}}{x-x^{\prime}}}ds_{\bf
x^{\prime}}$ $\displaystyle=$
$\displaystyle{|\mathbf{Q}_{k}-\mathbf{P}_{k}|}\int_{0}^{1}\frac{\mbox{\boldmath$\xi$}}{{({\mathbf{x}}-\mathbf{P}_{k})}-t{(\mathbf{Q}_{k}-\mathbf{P}_{k})}}dt$
$\displaystyle=$
$\displaystyle\frac{\mbox{\boldmath$\xi$}|\mathbf{Q}_{k}-\mathbf{P}_{k}|}{\mathbf{Q}_{k}-\mathbf{P}_{k}}\ln{\frac{{\mathbf{x}}-\mathbf{P}_{k}}{{\mathbf{x}}-\mathbf{Q}_{k}}}.$

From (44), we have

$\displaystyle\frac{\mbox{\boldmath$\xi$}|\mathbf{Q}_{k}-\mathbf{P}_{k}|}{\mathbf{Q}_{k}-\mathbf{P}_{k}}$
$\displaystyle=\frac{|\mathbf{Q}_{k}-\mathbf{P}_{k}|}{\mathbf{Q}_{k}-\mathbf{P}_{k}}\left(2\Re\\{1-\frac{1}{\lambda_{c}(\omega)}\\}a_{\nu_{k}}{\frac{i(\mathbf{Q}_{k}-\mathbf{P}_{k})}{|\mathbf{P}_{k}-\mathbf{Q}_{k}|}}+2\Re\\{\lambda_{c}(\omega)-1\\}a_{\tau_{k}}{\bf\frac{Q_{k}-P_{k}}{|P_{k}-Q_{k}|}}\right)$
$\displaystyle\quad=2\Re\\{\lambda_{c}(\omega)-1\\}a_{\tau_{k}}+i2\Re\\{1-\frac{1}{\lambda_{c}(\omega)}\\}a_{\nu_{k}}.$

Therefore, (43) can be simplified as

$\displaystyle\Re\\{\Phi({\bf
x})\\}=\Re\left\\{\sum_{k=1}^{N_{C}}{\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})}\ln{\frac{{\mathbf{x}}-\mathbf{Q}_{k}}{{\mathbf{x}}-\mathbf{P}_{k}}}+\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Re}(\omega,\delta_{D})}\frac{1}{{\mathbf{x}}-\mathbf{z}_{k}}\right\\},$
(45)

where $\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})$ and
$\mathfrak{D}_{k}^{\Re}(\omega,\delta_{D})$ are the quantities defined in (37)
and (39). From (45), the real part of $\Phi$ can be viewed as the real part of
the meromorphic function $G^{\Re}({\mathbf{x}})$ given by

$\displaystyle
G^{\Re}({\mathbf{x}}):=\sum_{k=1}^{N_{C}}{\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})}\ln{\frac{{\mathbf{x}}-\mathbf{Q}_{k}}{{\mathbf{x}}-\mathbf{P}_{k}}}+\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Re}(\omega,\delta_{D})}\frac{1}{{\mathbf{x}}-\mathbf{z}_{k}}.$
(46)

Since $G^{\Re}({\mathbf{x}})$ is holomorphic except at the points
$\mathbf{P}_{k},\mathbf{Q}_{k},\mathbf{z}_{k}$ and the segments
$\mathbf{P}_{k}\mathbf{Q}_{k}$ (see, for instance, [27]), it has complex
derivative near $\partial\Omega$ in the complex plane:

$\displaystyle\frac{dG^{\Re}({\mathbf{x}})}{d{\mathbf{x}}}=\sum_{k=1}^{N_{C}}{\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})}\left(\frac{1}{{\mathbf{x}}-\mathbf{Q}_{k}}-\frac{1}{{\mathbf{x}}-\mathbf{P}_{k}}\right)-\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Re}(\omega,\delta_{D})}\frac{1}{({\mathbf{x}}-\mathbf{z}_{k})^{2}}.$
(47)

A similar argument applies for the imaginary part of $\Phi({\mathbf{x}})$. ∎

The followings remark on Theorem 3.5 is in due.

######  Remark 3.6.

According to Theorem 3.5, both $G^{\Re}({\mathbf{x}})$ and
$G^{\Im}({\mathbf{x}})$ can be viewed as known quantities from the knowledge
of $(-\frac{1}{2}I+\mathcal{K}_{\Omega})[u^{\omega}-u_{0}]$ on
$\partial\Omega$. This theorem states that $\frac{dG^{\Re}({\bf x})}{d{\bf
x}}$ is a meromorphic function in $\mathbb{C}$ with simple poles at the
endpoints $\mathbf{P}_{k},~{}\mathbf{Q}_{k}$ of the segments $\mathcal{L}_{k}$
and poles of order $2$ at the center $\mathbf{z}_{k}$ of $D_{k}$. Hence, the
residues of $\frac{dG^{\Re}({\bf x})}{d{\bf x}}$ at the endpoints are given by

$\displaystyle\mbox{Res}\left(\frac{dG^{\Re}({\bf x})}{d{\bf x}},{\bf
Q_{k}}\right)=\mathfrak{C}_{k}^{\Re}(\omega,\delta_{k})=-\mbox{Res}\left(\frac{dG^{\Re}({\bf
x})}{d{\bf x}},{\bf P_{k}}\right).$ (48)

The information of the center of $D_{k}$ is contained in the following
function

$\displaystyle w({\bf x})$ $\displaystyle:=$
$\displaystyle\sum_{k=1}^{N_{D}}{\mathfrak{D}_{k}^{\Re}}(\omega,\delta_{D})\frac{1}{({\mathbf{x}}-\mathbf{z}_{k})^{2}}.$
(49)

Then the function $\frac{w^{\prime}({\bf x})}{w({\bf x})}$ will have simple
poles at the poles of $w({\mathbf{x}})$. Hence, these center points can be
identified from boundary measurements [16].

The above analysis shows that the effect of thin insulators on the boundary
data highly depends on the frequency, while the effect of small conductive
objects does not depend on the operating frequency that much. This relation
leads to the assertion that we can detect the small conductors when the
frequency is very high and both thin insulators and small conductors when the
frequency decreases. Numerical simulations in the later part of this paper
will illustrate this important observation. We can similarly analyze
$\mathfrak{C}_{k}^{\Im}(\omega,\delta_{k})$ and
$\mathfrak{D}_{k}^{\Im}(\omega,\delta_{D})$. For the low frequency case, we
have the following results.

###  3.3 Low-frequency case: $|\lambda_{c}(\omega)|\approx 0$ and
$\delta_{k}\approx 0$ with $|\lambda_{c}(\omega)|^{-1}\delta_{k}\approx\beta$
and $0<\beta<\infty$

In the low-frequency case, the admittivity contrast $\lambda_{c}(\omega)$ is
getting close to zero. As thickness $\delta_{k}$ goes to zero, we suppose that

$\frac{1}{|\lambda_{c}(\omega)|}\delta_{k}\approx\beta.$ (50)

Then by approximation (7), the potential jump along each thin insulator can
not be ignored. According to [1, 5, 9, 21], we have the following asymptotic
expansion formula of the potential $u^{\omega}$ for low frequency current.

######  Theorem 3.7 (Asymptotic expansion at low frequencies).

In the low-frequency regime, we have the following asymptotic formula for the
boundary perturbations of the potential $u^{\omega}$:

$\displaystyle\indent\displaystyle\left(-\frac{1}{2}I+K_{\Omega}\right)[u^{\omega}-u_{0}](x)=-\delta_{D}^{2}\sum_{k=1}^{N_{D}}\nabla\Gamma(x,z_{k})\cdot
M(\lambda_{d},B_{k})\nabla u_{0}(z_{k})$
$\displaystyle\indent\displaystyle+\sum_{k=1}^{N_{C}}\int_{\mathcal{L}_{k}}\frac{\partial\Gamma(x,x^{\prime})}{\partial\nu_{x^{\prime}}}[u^{\omega}]_{k}(x^{\prime})dx^{\prime}+O(\delta^{2}_{k})+O(\delta_{D}^{3}).$
(51)

In this case, since the potential jump
$[u^{\omega}]_{k}(x)=\frac{2\delta_{k}}{\lambda_{c}(\omega)}\frac{\partial
u^{\omega}}{\partial\nu}(x-\delta_{k}\nu_{x})|_{+}$ along $\mathcal{L}_{k}$ is
very large and could not be ignored, the effect of small conductors on the
perturbations of the boundary voltage is hidden by the insulators. Although we
cannot write (3.7) in an explicit way, we know that it is related with the
endpoints as well as the potential jump along the thin insulators. When
multiple thin insulators are well separated from each other, we can always
image them from boundary measurements. However, small conductors at low
frequencies are invisible since the thin insulators will dominate the boundary
measurements. Note that (50) gives the optimal range of operating frequencies
to use for optimal imaging of the thin insulators.

###  3.4 Spectroscopic analysis

Based on the above analysis in low- and high-frequency regimes, we
mathematically derived the frequency dependency of the current-voltage data in
a rigorous way. The current-voltage data is mainly affected by the outermost
thin insulators when the frequency is low, whereas the data mainly depends on
the small conductors when the frequency is high. With this observation, we can
detect the outermost thin insulators at low frequencies. As the frequency
increases, the small conductors become gradually visible whereas thin
insulators fade out (thicker insulator fades out at higher frequency than
thinner one). Hence, multi-frequency EIT system allows to probe these
frequency dependent behavior. On the other hand, we could not give a clear
range of the frequency to best identify thin objects or small conductors since
the penetrating frequency is not only depending on thickness $\delta_{k}$ but
also relying on admittivity contrast between $\gamma_{c}$ and $\gamma_{b}$,
geometry, relative size between $\delta_{k}$ and domain $\Omega$, position
distribution, etc. We plan to deal with this challenging issue in the future.

##  4 Numerical simulations

In this section, we will verify our mathematical analysis through various
numerical simulations by using multi-frequency Electrical Impedance Tomography
(mfEIT). The rough procedure of mfEIT using the standard sensitivity matrix
[14, 26] is as follows:

  1. 1.

$N_{E}$ electrodes
${\mathcal{E}}_{1},{\mathcal{E}}_{2},\cdots,{\mathcal{E}}_{N_{E}}$ are
attached on the boundary of $\Omega$ with a unform distance between adjacent
electrodes. Inject current between all adjacent pair of electrodes
${\mathcal{E}}_{k},{\mathcal{E}}_{k+1}$ at various angular frequencies
($\omega_{1},\omega_{2},\cdots,\omega_{N_{\omega}}$).

  2. 2.

Solve the forward problem using finite element method (FEM) to compute the
potential $u^{\omega}_{k}$ due to $k$-th injection current and collect
simulated boundary voltage data
$V_{\omega}=(V_{\omega}^{1,1},\cdots,V_{\omega}^{N_{E},N_{E}})$, where
$V^{k,j}$ is the $j$th boundary voltage subject to the $k$th injection
current, which is given by
$V^{j,k}_{\omega}=\int_{\Omega}\gamma^{\omega}\nabla u_{k}^{\omega}\cdot\nabla
u_{j}^{\omega}dx$. (In real experiment, we do not use
$V^{k-1,k}_{\omega},V^{k,k}_{\omega},V^{k,k+1}_{\omega}$ due to the unknown
contact impedance.)

  3. 3.

Discretizing the domain $\Omega$ into $M$ elements or pixels as
$\Omega=\sum_{k=1}^{M}{\mathcal{T}}_{k}$, compute the standard sensitivity
matrix (see [14, 26]). We use the standard linearized method to reconstruct
the admittivity images
$\delta\gamma_{\omega}=\\{\delta\gamma^{\omega}_{{\mathcal{T}}_{1}},\delta\gamma^{\omega}_{{\mathcal{T}}_{2}},\cdots,\delta\gamma^{\omega}_{{\mathcal{T}}_{M}}\\}^{T}$
from the boundary data $V_{\omega_{j}}$. Here, the reference homogeneous
admittivity was subtracted.

The numerical simulations are performed on a unit disk
$\Omega=\\{(x,y):x^{2}+y^{2}<1\\}$ with 16 electrodes equally spaced around
its circumference. Inside the unit disk, we consider the following numerical
models with two conductive objects and several thin insulating objects inside
as shown in Figure 4.

(a)(b)(c) Figure 4: Two different numerical models: (a) Two thin insulating
inhomogeneities ${\mathcal{C}}_{1},~{}{\mathcal{C}}_{2}$ with same thickness
and two small conductive objects and (b) two small conductive objects $D_{1}$
and $D_{2}$, and a rectangular insulating inhomogeneity encircling $D_{2}$;
(c) Admittivity distribution for thin insulating inhomogeneities
${\mathcal{C}}_{1},~{}{\mathcal{C}}_{2}$ and small conductive objects
$D_{1},~{}D_{2}$.

In Figure 4 (a), two thin insulating inhomogeneities with the same thickness
and two conductive objects are placed inside a homogenous domain so that we
can investigate the effect of current frequencies on the reconstructed images.
Thereafter, in Figure 4 (b), we will take use of a simulation model with two
conductive objects encircled by a rectangular insulating object
${\mathcal{C}}_{1}$. The thickness of the three thin insulating objects is
$5*10^{-4}$ in Figure 4 (a, b). The admittivity distribution is a piecewise
constant in each subdomain as shown in Figure 4 (c). A wide range of current
frequencies is applied and the resulting boundary potential are collected. The
selected spectroscopic images of admittivity distribution are shown in Figure
5 where the admittivity distribution is reconstructed at frequencies $\omega=$
10Hz, 1kHz, 50kHz, 150kHz, 250kHz and 500kHz.

###  4.1 Multi-frequency images

Figure 5 shows the reconstructed spectroscopic images when the simulation
model is Figure 4 (a); When current is injected at low frequencies
($\frac{\omega}{2\pi}=$ 10Hz, 1kHz), two small conductive objects are
invisible since the low frequency currents are blocked by the thin insulating
inhomogeneities and the boundary potential is mostly influenced by the thin
insulating objects.

Figure 5: Reconstructed admittivity distribution when current is injected at
frequencies $\omega=$ 10Hz, 1kHz, 50kHz, 150kHz, 250kHz, 500kHz. The second
row contains reconstructed images for normalized conductivity $\sigma$ (S/m),
and the third row is reconstructed images for normalized $\omega\epsilon$
(S/m).

As frequency increases, the currents at frequencies 50kHz and 150kHz begin to
penetrate the thin insulating objects and conductive objects start appearing
in the reconstructed images. Finally, when the currents are applied at
frequencies 250kHz and 500kHz, both conductive objects are visible in the
reconstructed images because the potential is mainly affected by the
conductive objects $D_{1},~{}D_{2}$.

Figure 6: Reconstructed admittivity distribution when current is injected at
frequencies $\omega=$ 10Hz, 1kHz, 50kHz, 150kHz, 250kHz, 500kHz. The second
row contains reconstructed images for normalized conductivity $\sigma$ (S/m),
and the third row is reconstructed images for normalized $\omega\epsilon$
(S/m).

As to the second numerical model in Figure 4 (b), we illustrate the
reconstructed spectroscopic images in Figure 6. Similarly, conductive object
surrounded by the rectangular insulating object is not visible at low
frequencies 10kHz and 1kHz because the currents can not pass through the thin
insulating objects. On the other hand, when frequency increases from 50kHz to
250kHz, the insulating rectangular ’wall’ starts to fade out because the
current begins to penetrate insulating object and the boundary voltage is
influenced by the internal conductive object. When the frequency is very high
(500kHz), the insulating rectangular object ${\mathcal{C}}_{1}$ disappears and
only conductive objects are visible in the reconstructed image. It is due to
the fact that the data are mainly influenced by conductive objects. From the
point of view of inverse problem, the spectroscopic images can give the
spectroscopic information for the thin insulating objects
${\mathcal{C}}_{1},{\mathcal{C}}_{2}$ by providing how much the boundary
potential is affected by thin insulating inhomogeneities at various
frequencies. The influence of thin insulating objects on the boundary
measurements decreases as the frequency increase. The dominant influence at
low frequencies is from insulating objects while the boundary potential is
dominated by conductive objects at high frequencies.

###  4.2 Fusion of multi-frequency images

Now we are considering to construct an integrated image based on an
investigation of principal component analysis (PCA) where the information of
both thin insulating objects and small conductive objects can be extracted
from the integrated image. PCA is a useful tool in extracting the dominant
features (principal components) from a set of reconstructed images at various
frequencies [24, 28]. In order to implement this approach, we first obtain the
reconstructed images at a broad range of frequencies
$\omega=\omega_{1},\omega_{2},\cdots,\omega_{N_{\omega}}$ and each
reconstructed admittivity. Then we represent each image by an $M\times 1$
column vector
${\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega_{j}}},j=1,2,\cdots,N_{\omega}$
where
${\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega_{j}}}=\\{\delta\gamma^{\omega_{j}}_{{\mathcal{T}}_{1}},\delta\gamma^{\omega_{j}}_{{\mathcal{T}}_{2}},\cdots,\delta\gamma^{\omega_{j}}_{{\mathcal{T}}_{M}}\\}^{T}$.
Then average admittivity at $\omega$ can be defined as

$\overline{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega}}}=\frac{1}{N_{\omega}}\sum_{j=1}^{N_{\omega}}{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega_{j}}}.$
(52)

We define the mean oscillation of
${\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega}}$ by the
vector

$\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{j}}}={\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega_{j}}}-\overline{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{\omega}}}.$
(53)

Since the real part of admittivity is $\sigma$ and its imaginary part is
$\omega\epsilon$, we resolve the set of admittivities and deal with the real
and imaginary part separately.

Let $\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{j}}}$ be
the real part of
$\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\gamma$}}_{j}}}$. This set
of vectors is then subject to principal component analysis which tries to find
a set of $N_{\omega}$ orthogonal vectors and their associated eigenvalues.
Both information can best describe the distribution of the admittivity. The
vectors and scalars are the eigenvectors and eigenvalues of the covariance
matrix

$\mathbf{C_{\sigma}}=\frac{1}{N_{\omega}}\sum_{j=1}^{N_{\omega}}\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{j}}}~{}\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{j}}}^{T}=\mbox{\boldmath{$\hat{\delta\sigma}~{}\hat{\delta\sigma}^{T}$}},$
(54)

where
$\mbox{\boldmath{$\hat{\delta\sigma}$}}=[\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{1}}}~{}\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{2}}}~{}\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{3}}}\cdots\hat{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{N_{\omega}}}}]$
is of size $M\times N_{\omega}$ and the covariance matrix
$\mathbf{C_{\sigma}}$ is of size $M\times M$ with $M$ the pixel number of the
admittivity distribution. Using the singular value decomposition, the
covariance matrix $\mathbf{C_{\sigma}}$ can be written as

$\mathbf{C_{\sigma}}=\mbox{\boldmath{$\hat{\delta\sigma}~{}\hat{\delta\sigma}^{T}$}}={\mbox{\boldmath{$U$}}}{\mbox{\boldmath{$\Lambda$}}}{\mbox{\boldmath{$U$}}^{T}}=\sum_{i=1}^{M}\lambda_{i}\mathbf{u}_{i}\mathbf{u}_{i}^{T},$
(55)

where $\lambda_{i}$ is the eigenvalue of
$\hat{\delta\sigma}\hat{\delta\sigma}^{T}$ with
$\lambda_{1}\geq\lambda_{2},\cdots\geq\lambda_{M}$ and $\mathbf{u}_{i}$ is the
corresponding eigenvector. Note that in practice, the number of frequencies
($N_{\omega}$) is much smaller than the pixel number ($M$) of the
reconstructed images. Therefore, we consider $N(<N_{\omega})$ leading
eigenvectors of the covariance matrix that are chosen as those with the
largest associated eigenvalues, where $N$ depends on the signal-to-noise ratio
(SNR). Then the principal components for any admittivity image
$\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega_{j}},j=1,\cdots,N_{\omega}$,
can be written as

$p_{i}={\bf
u}_{i}^{T}(\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega_{j}}-\overline{{\mbox{\boldmath{$\delta$}}}{\mbox{\boldmath{$\sigma$}}_{\omega}}}).$
(56)

for $i=1,2,\cdots,N$. Here $p_{i}$ represents the data projected into the
$N-$dimensional space of eigenvectors. In order to enhance the image details,
we rewrite $\hat{\delta\sigma}$ in terms of singular value decomposition

$\mbox{\boldmath{$\hat{\delta\sigma}$}}=\sum_{i=1}^{N_{\omega}}\sqrt{\lambda_{i}}{\mbox{\boldmath{$u$}}_{i}}{\mbox{\boldmath{$v$}}_{i}}^{T},$
(57)

where ${\mbox{\boldmath{$v$}}_{i}}$ is the eigenvectors of
$\mathbf{C_{\sigma}}^{T}$ and
${\mbox{\boldmath{$u$}}_{i}}=\mbox{\boldmath{$\hat{\delta\sigma}$}}{\mbox{\boldmath{$v$}}_{i}}$.
We construct the following matrix from the principal components of
$\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega_{j}}$

$\tilde{\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma_{\omega}$}}}=\sum_{i=1}^{N}{\mbox{\boldmath{$u$}}_{i}}{\mbox{\boldmath{$v$}}_{i}}^{T},$
(58)

where $\tilde{\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma_{\omega}$}}}$
is a matrix of size $M\times N_{\omega}$. The new admittivity distribution can
be obtained by

${\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega}}=\frac{1}{N_{\omega}}\sum_{j=1}^{N_{\omega}}\tilde{\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega,j}},\textsection$
(59)

where
$\tilde{\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma$}}_{\omega,j}}$ is
the $j-$th column of
$\tilde{\mbox{\boldmath{$\delta$}}\mbox{\boldmath{$\sigma_{\omega}$}}}$.
Similarly, we can obtain the imaginary part of integrated admittivity
distribution. In our case, only $N=2$ eigenvectors corresponding to the two
largest eigenvalues are chosen to obtain the new images. The images
corresponding the configurations in Figure 4 (a,b) are shown respectively in
Figures 7 and 8. The integrated images show both on the thin insulating and
the small conductive objects.

Figure 7: Integrated image of real and imaginary part of normalized
admittivity distribution for numerical model in Figure 4 (a). Figure 8:
Integrated image of real and imaginary part of normalized admittivity
distribution for numerical model in Figure 4 (b).

##  5 Phantom experiments

In this section, we present phantom experiments by using 32-channel mfEIT
system (EIT-Pioneer Set made by Swisstom, Switzerland) to illustrate the
frequency dependent behavior of the reconstructed images. The available
injection current frequency of the Swisstom EIT-Pioneer Set is set between 50
kHz and 250 kHz.

Figure 9: Configuration of phantom experiments by using Swisstom EIT-Pioneer
Set.

Figure 9 shows the configuration of phantom experiment. We use a cylindrical
tank with 360 mm in diameter and 32 equally-spaced electrodes are attached.
The tank is filled with agar-gelatin mixture. Inside the phantom, there are
two conductive objects with a diameter of 30 mm and four very thin kitchen
wrap with 1 $\mu$m thickness. One conductive object is encircled by four very
thin kitchen wrap. We injected current of 1 mA at various frequencies of 50
kHz, 100 kHz, 150 kHz, 180kHz, 200 kHz and 250 kHz. Figure 10 (a) presents the
reconstructed images at six different frequencies. The insulating wraps appear
as a solid object from 50 kHz to 100 kHz since currents cannot penetrate the
insulating wrap. As we expected by numerical experiments in the previous
section, the insulating wrap starts to fade out at frequency 150 kHz and
totally disappear at 250kHz. All the experimental results shown in Figures 10
(a) are consistent with the numerical simulations in the previous section.

Figure 10: (a) Spectroscopic images of normalized admittivity distribution
from low frequencies to high frequencies; (b) Fusion of multi-frequency
images.

As for the fusion of multi-frequency images, the experiments are conducted at
a wide range of frequencies from 50kHz to 250kHz with a step size 10kHz since
the minimum frequency step of Swisstom EIT-Pioneer Set is fixed at 10kHz. The
integrated image is shown in Figure 10 (b) where both two conductive objects
and the surrounding insulating wrap are visible.

##  6 Conclusion

In this work, we have provided a rigorous mathematical formula for the jump of
potential and normal derivative across the thin linear-shaped insulating
inhomogeneities based on layer potential techniques. The potential jump is
related with the thickness of insulating objects as well as the current
frequencies. Based on these jump conditions, we have developed two asymptotic
expansions for current-voltage data perturbations due to thin insulators and
small conductors at various frequencies. Using these two asymptotic
expansions, we have mathematically shown that at high frequencies we can
visualize the small conductors, while at low frequencies we can only get the
information of thin insulators. Based on this mathematical analysis, we
conclude that multiple frequencies help us to handle the spectroscopy behavior
of the current-voltage data with respect to thin insulators and small
conductors. When the frequency increase from very low to very high, we can
continuously observe the images of thin insulators (low frequency), both thin
insulators and small conductors (not too low, not too high frequency), and
only small insulators (high frequency). The mathematical results are supported
by a variety of numerical illustrations and phantom experiments.

## Acknowledgements

Ammari was supported by the ERC Advanced Grant Project MULTIMOD–267184. Seo
and Zhang were supported by the National Research Foundation of Korea (NRF)
grant funded by the Korean government (MEST) (No. 2011-0028868,
2012R1A2A1A03670512)

## References

## References

  * [1] H. Ammari, E. Beretta and E. Francini, _Reconstruction of thin conductivity imperfections, II. the case of multiple segments_ , Appl. Anal., 85(2006), pp. 87-105.
  * [2] H. Ammari and H. Kang, _Reconstruction of small inhomogeneities from boundary measurements_ , Lecture Notes in Math. 1846, Springer-Verlag, Berlin, (2004).
  * [3] H. Ammari and H. Kang, _Polarization and moment tensors with applications to inverse problems and effective medium theory_ , Appl. Math.Sci., Springer-Verlag, New York, (2007).
  * [4] H. Ammari, H. Kang, E. Kim and M. Lim, _Reconstruction of closely spaced small inclusions_ , SIAM J. Numer. Anal., 42(2005), pp. 2408-2428.
  * [5] H. Ammari, and J. K. Seo, _An accurate formula for the reconstruction of conductivity inhomogeneities_ , Adv. Appl. Math., 30(2003), pp. 679-705.
  * [6] L. Baratchart, J. Leblond, F. Mandrea and E. B. Saff, _How can the meromorphic approximation help to solve some 2D inverse problems for the Laplacian?_ , Inverse Problems, 15(1999) , pp. 79-90.
  * [7] E. Beretta, E. Francini and M. S. Vogelius, _Asymptotic formulas for steady state voltage potentials in the presence of thin inhomogeneities. a rigorous error analysis_ , J. Math. Pures Appl., 82(2003), pp. 1277-1301.
  * [8] E. Beretta, A. Mukherjee and M. S. Vogelius, _Asymptotic formulas for steady state voltage potentials in the presence of conductivity imperfections of small area_ , Z. Angew. Math. Phys., 52(2001), pp. 543-572.
  * [9] A. Friedman and M. S. Vogelius, _Identifcation of small inhomogeneities of extreme conductivity by boundary measurements: a theorem on continuous dependence_ , Arch. Rational Mech. Anal., 105(1989), pp. 299-326.
  * [10] H. Fricke, _A mathematical treatment of the electrical conductivity of colloids and cell suspensions_ , J. Gen. Physiol., 6(1924), pp. 375-384.
  * [11] S. Gabriel, R. W. Lau and C. Gabriel, _The dielectric properties of biological tissues: II. measurements in the frequency range 10Hz to 20GHz_ , Phys. Med. Biol., 41(1996), pp. 2251-2269.
  * [12] H. Griffiths, H.T.L. Leung and R. L. Williams, _Imaging the complex impedance of the thorax_ , Clin. Phys. Physiol. Meas., 13 A(1992), pp. 77-81.
  * [13] R. Halter, A. Hartov and K. D. Paulsen, _A Broadband High-Frequency Electrical Impedance Tomography System for Breast Imaging_ , IEEE Trans. Biomed. Eng., 55(2008), pp. 650-59.
  * [14] D. Holder, _Electrical impedance tomography: methods, history and applications_ , Bristol, UK, IOP Publishing, (2005).
  * [15] H. Jain, D. Isaacson, P.M. Edic and J.C. Newell, _Electrical impedance tomography of complex conductivity distributions with noncircular boundary_ , IEEE Trans. Biomed. Eng., 44 (1997), pp. 1051-60.
  * [16] H. Kang and H. Lee, _Identifcation of simple poles via boundary measurements and an application of EIT_ , Inverse Problems, 20(2004), pp. 1853-1863.
  * [17] D. M. McCann and M. C. Forde, _Review of ndt methods in the assessment of concrete and masonry structures_ , NDT&E Int., 34(2001), pp. 71-84.
  * [18] M. Pavlin and D. Miklav$\check{\mbox{c}}$i$\check{\mbox{c}}$, _Effective conductivity of a suspension of permeabilized cells: a theoretical analysis_ , Biophys. J., 85(2003), pp. 719-729.
  * [19] M. Pour-Ghaz and J. Weiss, _Application of frequency selective circuits for crack detection in concrete elements_ , J. ASTM. Int., 8(2011), pp. 1–11.
  * [20] M. Pour-Ghaz, M. Niemuth and J. Weiss, _Use of electrical impedance spectroscopy and conductive surface films to detect cracking and damage in cement based materials,_ Structural Health Monitoring Technologies - Part I. Glisic, B., Ed.; American Concrete Institute Special Publication: Farmington Hills , MI, USA, (2012).
  * [21] A. Khelifi and H. Zribi, _Asymptotic expansions for the voltage potentials with two-dimensional and three-dimensional thin interfaces_ , Math. Meth. Appl. Sci., 34(2011), pp. 2274-2290.
  * [22] S. Kim, J. K. Seo, and T. Ha, A nondestructive evaluation method for concrete voids: frequency differential electrical impedance scanning, SIAM J. Appl. Math., 69(2009), pp. 1759-1771.
  * [23] S. Kim, E. J. Lee, E. J. Woo and J. K. Seo, _Asymptotic analysis of the membrane structure to sensitivity of frequency-difference electrical impedance tomography_ , Inverse Problems, 28(2012).
  * [24] A. Sophian, Gui Yun Tian, David Taylor and John Rudlin, _A feature extraction technique based on principal component analysis for pulsed eddy current NDT_ , NDT&E International 36(2003), pp. 37-41.
  * [25] J.K. Seo, J. Lee, S.W. Kim, H. Zribi, and E. J. Woo E J, _Frequency-difference electrical impedance tomography (fdEIT): algorithm development and feasibility study_ , Physiol. Meas., 29 (2008), pp. 929-944.
  * [26] J. K. Seo and E. J. Woo, _Nonlinear inverse problems in imaging_ , John Wiley & Sons, 2012.
  * [27] E. M. Stein and R. Shakarchi, _Complex analysis_ , Princeton University Press, 2003.
  * [28] M. A. Turk and A. P. Pentland, _Face recognition using eigenfaces_ , In Proceedings of the Computer Society Conference on Computer Vision and Pattern Recognition, (1991), pp. 586-591.
  * [29] A. J. Wilson , P. Milnes, A. R. Waterworth, R. H. Smallwood and B. H. Brown, _Mk3.5: a modular, multi-frequency successor to the Mk3a EIS/EIT system_ , Physiol. Meas., (2001), pp. 49-54.
  * [30] T. Zhang, L. Zhou, H. Ammari and J. K. Seo, _Electrical impedance spectroscopy-based defect sensing technique in estimating cracks_ , Sensors, 15(2015), pp. 10909-10922.