#  Shape reconstructions by using plasmon resonances

Ming-Hui Ding School of Mathematics, Hunan University, Changsha 410082, China.
Email<EMAIL_ADDRESS>Hongyu Liu Department of Mathematics, City
University of Hong Kong, Kowloon, Hong Kong, China. Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Guang-Hui Zheng School of
Mathematics, Hunan University, Changsha 410082, China. Email:
<EMAIL_ADDRESS>
Abstract

We study the shape reconstruction of an inclusion from the faraway measurement
of the associated electric field. This is an inverse problem of practical
importance in biomedical imaging and is known to be notoriously ill-posed. By
incorporating Drude’s model of the permittivity parameter, we propose a novel
reconstruction scheme by using the plasmon resonance with a significantly
enhanced resonant field. We conduct a delicate sensitivity analysis to
establish a sharp relationship between the sensitivity of the reconstruction
and the plasmon resonance. It is shown that when plasmon resonance occurs, the
sensitivity functional blows up and hence ensures a more robust and effective
construction. Then we combine the Tikhonov regularization with the Laplace
approximation to solve the inverse problem, which is an organic hybridization
of the deterministic and stochastic methods and can quickly calculate the
minimizer while capture the uncertainty of the solution. We conduct extensive
numerical experiments to illustrate the promising features of the proposed
reconstruction scheme.

Keywords: Shape reconstruction; plasmon resonance; sensitivity analysis;
Tikhonov regularization; Laplace approximation

2010 Mathematics Subject Classification: 65N21, 35R30, 78A46, 65H50

##  1 Introduction

Plasmon resonance is the resonant oscillation of conduction electrons at the
interface between negative and positive permittivity material stimulated by
incident field. Plasmonics is revolutionizing many light-based technologies
via electron oscillations in metals. We refer to [32, 10, 46, 49, 43, 11] and
the references cited therein for many striking optical, phononic, biomedical,
diagnostic and therapeutic applications in the physical literature. Recent
studies have revealed the deep and intriguing connection between the plasmon
resonance and the spectral study of the Neumann-Poincaré operator [2, 3, 5,
20, 41, 27, 38, 12]. In addition, there are many theoretical understandings
and conceptual proposals about plasmonic devices.

In [8], by analyzing the imaginary part of the Green function, it is shown
that one can achieve super-resolution and super-focusing by using plasmonic
nanoparticles. In [12, 2, 3], it is shown that the plasmon resonance
concentrates and localises at high-curvature places, which can provide
potential application in super-resolution imaging of plasmon particles. We
would also like to mention in passing some related studies on plasmonic
cloaking [5, 9, 20, 35, 42, 22, 36, 37, 34, 16, 52]. In this paper, we study
the shape reconstruction of an inclusion from the faraway measurement of the
associated electric field. This is an inverse problem of practical importance
in biomedical imaging and is known to be notoriously ill-posed. By
incorporating Drude’s model of the permittivity parameter, we propose a novel
reconstruction scheme by using the plasmon resonance with a significantly
enhanced resonant field.

We next introduce the mathematical formulation of the inverse shape problem
for our study. Let $D\subset\mathbb{R}^{2}$ be a bounded domain with a
connected complement $\mathbb{R}^{2}\backslash\overline{D}$. Given a harmonic
function $H$, we consider the following electrostatic problem:

$\displaystyle\begin{cases}&\nabla\cdot(\varepsilon\nabla u(x))=0\ \ \ \ \ \
\mathrm{in}\ \ \mathbb{R}^{2},\\\ &u(x)-H(x)=O(|x|^{-1})\ \ \mbox{as}\ \
|x|\rightarrow+\infty,\end{cases}$ (1.1)

where

$\displaystyle\varepsilon=\varepsilon_{D}\chi(D)+\varepsilon_{m}\chi(\mathbb{{R}}^{2}\backslash\overline{{D}}),$
(1.2)

and $\chi$ is the characteristic function. (1.1)-(1.2) describes the
transverse electromagnetic propagation in the quasi-static regime. $u$
signifies the transverse electric field and $\varepsilon$ signifies the
permittivity parameter of the medium. Throughout, we shall assume that the
background parameter $\varepsilon_{m}$ is a positive constant, whereas
$\varepsilon_{D}$ is a complex-valued function of the illuminating frequency
and fulfils the Drude’s model. We shall supply more details about the Drude
model in what follows. The shape reconstruction problem can be formulated as
follows:

_Inverse Problem (IP)_ : Identify the shape of the inclusion, namely $\partial
D$, from the measurement data $u^{s}=u-H$ on $\partial\Omega$ with
$D\Subset\Omega$ associated with a fixed incident field $H$. For simplicity,
we take $\Omega$ to be a central ball of radius $r\in\mathbb{R}_{+}$ with $r$
sufficiently large. Hence, the measurement represents the far-field pattern of
the electric field.

The shape reconstruction problem introduced above is severely ill-posed and
highly nonlinear. First, it is well known that due to the diffraction limit,
the far field excited by the object carries information on a scale much larger
than the operating wavelength, while information on a scale smaller than the
operating wavelength is confined near the object itself. In addition, the
scattering information in the quasi-static regime is very weak, and in the
presence of measurement noise, the signal-to-noise ratio in the far field is
low and signal distortion is serious [1]. We also refer to [3, 4, 28, 30, 26,
18, 13, 14, 15, 21, 39, 40, 17, 51] for related studies in the literature on
this inverse shape problem.

In this article, we first perform a shape sensitivity analysis and derive the
shape sensitivity functional with respect to domain perturbation by a delicate
asymptotic analysis. We establish the spectral expansion of the shape
sensitivity functional, from which we can conclude the sharp relationship
between the reconstruction sensitivity and the plasmon resonance. It indicates
that the plasmon resonant field can render a more robust and effective
reconstruction. Moreover, in order to overcome the ill-posedness, we combine
the Tikhonov regularization method with the Laplace approximation (LA) to
solve the inverse problem. This hybrid method is essentially the organic
combination of the deterministic method and stochastic method, which can
rapidly calculate the minimizer (Maximum a posteriori estimation point (MAP))
and capture statistical information of the solution more effectively. To
provide a global view of our study, the major contributions of this work can
be summarised as follows.

  1. 1.

By using the layer-potential perturbation technique, we rigorously derive the
asymptotic expansion of the perturbed far field with respect to the shape
perturbation. Furthermore, we obtain the representation formula of the shape
sensitivity functional.

  2. 2.

Based on the spectral theory of the Neumann-poincaré operator, we establish
the delicate spectral expansion of the shape sensitivity functional. It
indicates that when plasmon resonance occurs, the shape sensitivity can be
improved dramatically.

  3. 3.

Due to the severe ill-posedness of inverse problem, we use plasmon resonance
to enhance the sensitivity, and then combine the Tikhonov regularization
method with the Laplace approximation to solve the inverse problem and
quantify the uncertainty of the solution. Compared with the standard method,
our numerical results show that the proposed method can significantly improve
the accuracy and robustness of the numerical reconstruction.

The rest of the paper is organized as follows. In Section 2, we provide
preliminary knowledge on layer potential operators and plasmon resonance. In
Section 3, we conduct the sensitivity analysis for the perturbed domain, and
derive the spectral expansion of the shape sensitivity functional. In Section
4, we discuss the combination of the Tikhonov regularization method and the
Laplace approximation. Sections 5 and 6 are respectively devoted to numerical
experiments and conclusion.

##  2 Preliminaries

###  2.1 Layer potentials and Neumann-Poincaré operator

We collect a number of preliminary results on the layer potentials, in
particular the Neumann-Poincaré operator for our subsequent use. Throughout
this paper, we consider a domain $D$ with a $C^{2}$ boundary. The $L^{2}$
inner product and the corresponding norm on $\partial D$ are denoted by
$\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ in short, respectively. The single
layer potential $\mathcal{S}_{D}$ and double layer potential $\mathcal{D}_{D}$
associated with $D$ are given by

$\displaystyle\mathcal{S}_{D}[\varphi](x)$ $\displaystyle:=\int_{\partial
D}\Gamma(x,y)\varphi(y)d\sigma(y),\ \ \ x\in\mathbb{R}^{2},$
$\displaystyle\mathcal{D}_{D}[\varphi](x)$ $\displaystyle:=\int_{\partial
D}\frac{\partial\Gamma(x,y)}{\partial\nu(y)}\varphi(y)d\sigma(y),\ \ \
x\in\mathbb{R}^{2}\setminus\partial D,$

where $\varphi\in L^{2}(\partial D)$ is the density function, and the Green
function $\Gamma(x,y)$ to the Laplacian in $\mathbb{R}^{2}$ is given by

$\displaystyle\Gamma(x,y)=\frac{1}{2\pi}\ln|x-y|.$

The notations $u|_{\pm}$ and $\frac{\partial u}{\partial\nu}|_{\pm}$ denote
the traces on $\partial D$ from the outside and inside of $D$, respectively.
The following jump relations hold [4, 5]:

$\displaystyle\frac{\partial\mathcal{S}_{D}[\varphi]}{\partial\nu}\bigg{|}_{\pm}(x)$
$\displaystyle=\left(\pm\frac{1}{2}I+\mathcal{K}^{*}_{D}\right)[\varphi](x),\
\ x\in\partial D,$

where $\mathcal{K}^{*}_{D}$ is known as the Neumann-Poincaré (NP) operator
defined by

$\displaystyle\mathcal{K}_{D}^{*}[\varphi](x)=\frac{1}{2\pi}\int_{\partial
D}\frac{\langle x-y,\nu(x)\rangle}{|x-y|^{2}}\varphi(y)d\sigma(y).$

Next, we recall some useful facts about the NP operator $\mathcal{K}^{*}_{D}$
[5, 20, 41].

######  Lemma 2.1.

(i) The Calder$\acute{o}$n identity holds:
$\mathcal{S}_{D}\mathcal{K}_{D}^{*}=\mathcal{K}_{D}\mathcal{S}_{D}$ on
$H_{0}^{-\frac{1}{2}}$, where $H_{0}^{-\frac{1}{2}}$ is the zero mean subspace
of $H^{-\frac{1}{2}}$;

(ii) The operator $\mathcal{K}_{D}^{*}$ is compact and self-adjoint in the
Hilbert space $H_{0}^{-\frac{1}{2}}$ equipped with the following inner product

$\displaystyle\langle\varphi,\psi\rangle_{\mathcal{H}^{*}(\partial
D)}=-\langle\mathcal{S}_{D}[\psi],\varphi\rangle_{\frac{1}{2},-\frac{1}{2}},$

with $\langle\cdot,\cdot\rangle_{-\frac{1}{2},\frac{1}{2}}$ being the duality
pairing between $H^{-1/2}(\partial D)$ and $H^{1/2}(\partial D)$

(iii) Let $\mathcal{H}^{*}(\partial D)$ be the space
$H_{0}^{-\frac{1}{2}}(\partial D)$ with the new inner product in (ii). Let
$(\lambda_{j},\varphi_{j}),j=0,1,2,...$ be the eigenvalue and normalized
eigenfunction pair of $\mathcal{K}_{D}^{*}$ in $\mathcal{H}^{*}(\partial D)$,
then $\lambda_{j}\in(-\frac{1}{2},\frac{1}{2})$ and $\lambda_{j}\rightarrow 0$
as $j\rightarrow\infty$;

(iv) The following representation formula holds: for any
$\varphi\in\mathcal{H}^{*}(\partial D)$,

$\displaystyle\mathcal{K}^{*}_{D}[\varphi]=\sum_{j=1}^{\infty}\lambda_{j}\langle\varphi,\varphi_{j}\rangle_{\mathcal{H}^{*}(\partial
D)}\varphi_{j}.$ (2.3)

###  2.2 Plasmon resonance

We next briefly discuss the mathematical framework of plasmon resonance. We
first give the form of the solution of equation ($\ref{conducitivity pro})$.
From [1, 5, 20], we have

$\displaystyle u(x)=H(x)+\mathcal{S}_{D}[\phi](x),\ \ \ \mathrm{for}\
x\in\mathbb{{R}}^{2},$

where $\phi\in L_{0}^{2}(\partial D):=\\{\phi\in L^{2}(\partial
D);\int_{\partial D}\phi=0\\}$ satisfies

$\displaystyle(\lambda I-\mathcal{K}^{*}_{D})[\phi](x)=\frac{\partial
H}{\partial\nu}\bigg{|}_{\partial D},\ \ x\in\partial D,$ (2.4)

with $\lambda$ given by

$\displaystyle\lambda=\frac{\varepsilon_{D}+\varepsilon_{m}}{2(\varepsilon_{D}-\varepsilon_{m})}.$
(2.5)

The permittivities of plasmon materials, such as noble metals, are different
from the ordinary materials and may possess negative real parts. In fact, the
electric permittivity $\varepsilon_{D}$ of the plasmon material is changing
with respect to the operating frequency $\omega$. The $\varepsilon_{D}$ can be
described by the Drude’s model (see [47, 6, 24]),

$\displaystyle\varepsilon_{D}=\varepsilon_{D}(\omega)=\varepsilon_{0}\left(1-\frac{\omega_{p}^{2}}{\omega(\omega+i\gamma)}\right),$
(2.6)

where $\varepsilon_{0}$ is the electric permittivity of the vacuum and is
given by $\varepsilon_{0}=9\cdot 10^{-12}F/m$. The plasmon frequency
$\omega_{p}$ and the damping parameter $\gamma$ are strictly positive. When
$\omega<\omega_{p}$ the real part of $\varepsilon_{D}$ can be negative. In
particular, for the gold, the value of the plasmon frequency and the damping
parameter are $\omega_{p}=2\cdot 10^{15}s^{-1}$ and $\gamma=10^{14}s$,
respectively. Let the operating frequency $\omega=6\cdot 10^{14}Hz$ (visible
light frequency), and from Drude’s model (2.6), the electric permittivity of
gold nanoparticle is calculated as
$\varepsilon_{D}\approx(-9.8108+1.8018i)\varepsilon_{0}$.

Now by applying the spectral decomposition of the $\mathcal{K}^{*}_{D}$ to the
integral equation $(\ref{density})$, the density $\phi\in L_{0}^{2}(\partial
D)$ becomes

$\displaystyle\phi=\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle_{\mathcal{H}^{*}(\partial
D)}}{\lambda-\lambda_{j}}\varphi_{j},$ (2.7)

where $\lambda_{j}$ are eigenvalues of $\mathcal{K}^{*}_{D}$ and they satisfy
$|\lambda_{j}|<\frac{1}{2}$. When the real part of $\varepsilon_{D}(\omega)$
is negative, it holds that $|\mathfrak{Re}(\lambda(\omega))<\frac{1}{2}|$,
where and also in what follows $\mathfrak{Re}(\lambda(\omega))$ signifies the
real part of $\lambda(\omega)$. The frequency $\omega$ is called a plasmon
resonance frequency if it satisfies

$\displaystyle\mathfrak{Re}(\lambda(\omega))=\lambda_{j},$ (2.8)

for some $j$, where $\lambda_{j}$ is an eigenvalue of the NP operator
$\mathcal{K}^{*}_{D}$. For this reason, (2.8) is called the plasmon resonance
condition. It is emphasized that only when the imaginary part
$\mathfrak{Im}(\lambda(\omega))=0$ (lossless metal), the resonance condition
(2.8) can imply the density $\phi$ in (2.7) blow up.

In (2.7) the density $\phi$ will be amplified when the plasmon resonance
condition is reached provided that $\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle_{\mathcal{H}^{*}(\partial D)}$ is nonzero.
As a result, the  $j$-th mode of the far field $u-H$ will show a resonant
behavior, and it is called that plasmon resonance occurs (see also [20]). As
an illustration, we consider the cylindrical metal nanorod as an example. When
the external electric field lines are parallel to the circular cross-section
of the cylindrical metal nanorod, the physical process can be described by the
two-dimensional mathematical model (1.1)–(1.2). Next, we set
$\varepsilon_{m}=\varepsilon_{0}$ (vacuum) and notice that, as $D$ is a disk,
the spectrum of $\mathcal{K}_{D}^{*}$ are $\\{0,\frac{1}{2}\\}$. If
$\lambda_{j}=0$, from the plasmon resonance condition (2.8) and the fact that
$\gamma$ is sufficiently small in (2.6), we nearly get
$\mathfrak{Re}(\varepsilon_{D}(\omega))=-\varepsilon_{m}$. This relationship
is called the Fröhlich condition (see [44]). Furthermore, by Drude’s model
(2.6), we can obtain the resonance frequency
$\omega=\sqrt{\frac{\omega_{p}^{2}}{2}-\gamma^{2}}$. When $\gamma=0$ (lossless
metal), the resonance frequency $\omega=\frac{\omega_{p}}{\sqrt{2}}$ or
$f=\frac{\omega_{p}}{2\pi\sqrt{2}}$, which is called the Fröhlich frequency
(see [47]).

##  3 Shape sensitivity analysis

###  3.1 Sensitivity analysis for the perturbed domain

We consider the sensitivity analysis for the shape reconstruction problem with
a perturbed domain, namely evaluating the effect of the domain variations on
the far-field measurement data.

Define the forward operator $\mathcal{F}:X\rightarrow
H^{\frac{1}{2}}(\partial\Omega)$ on a subset $X\subset C^{2}(\partial D)$:

$\displaystyle\mathcal{F}(\partial D)=u^{s}|_{\partial\Omega},\ \ \ x\in X,$

where $u^{s}|_{\partial\Omega}=(u-H)|_{\partial\Omega}$ is the far-field
measurement associated with (1.1) on the boundary $\partial\Omega$.

For a small $\epsilon\in\mathbb{R}_{+}$, we let $\partial D_{\epsilon}$ be an
$\epsilon$ -perturbation of $D$, i.e.,

$\displaystyle\partial D_{\epsilon}:=\\{\tilde{x}=x+\epsilon
h(x)\nu(x),x\in\partial D\\},$

where $h\in C^{1}(\partial D)$, and $\nu$ is the outward unit normal vector to
$\partial D$. The solution $u_{\epsilon}$ to (1.1) with $D_{\epsilon}$ has the
following representation formula

$\displaystyle
u_{\epsilon}=H(x)+\mathcal{S}_{D_{\epsilon}}[\tilde{\phi}_{\epsilon}](x),\ \
x\in\mathbb{R}^{2},$ (3.9)

where the density function $\tilde{\phi}_{\epsilon}$ is the solution to

$\displaystyle(\lambda
I-\mathcal{K}_{D_{\epsilon}}^{*})[\tilde{\phi}_{\epsilon}](\tilde{x})=\frac{\partial
H}{\partial\tilde{\nu}}\bigg{|}_{\partial D_{\epsilon}}\ \ \
\tilde{x}\in\partial D_{\epsilon}.$ (3.10)

Let $\Psi_{\epsilon}$ be the diffeomorphism from $\partial D$ to $\partial
D_{\epsilon}$ given by

$\Psi_{\epsilon}(x)=x+\epsilon h(x)\nu(x),\ \ \ x\in\partial D.$

Moreover, we denote $\tilde{\nu}$ the outward unit normal vector to $\partial
D_{\epsilon}$ and $d\tilde{\sigma}$ the line element of $\partial
D_{\epsilon}$. The following expansions of $\tilde{\nu}$ and $d\tilde{\sigma}$
hold [7]:

$\displaystyle\tilde{\nu}(\tilde{x})$ $\displaystyle=\nu(x)-\epsilon
h^{\prime}(x)T(x)+O(\epsilon^{2}),$ (3.11) $\displaystyle
d\tilde{\sigma}(\tilde{x})$
$\displaystyle=d\sigma(x)-\epsilon\tau(x)h(x)d\sigma(x)+O(\epsilon^{2}).$
(3.12)

Here and throughout the rest of the paper, $\tau(x)$ signifies the curvature
of $\partial D$ at $x$, $T$ is the unit tangential vector to $\partial D$, and
$h^{\prime}$ is the tangential derivative of $h$ on $\partial D$, i.e.,
$h^{\prime}=\frac{\partial h}{\partial T}$, where
$\frac{\partial}{\partial\tilde{\nu}}$ denotes the outward normal derivative
on $\partial D_{\epsilon}$.

In view of (3.9) and (3.10), in order to obtain the asymptotic formula of the
perturbed far field $u_{\epsilon}^{s}|_{\partial\Omega}:=\mathcal{F}(\partial
D_{\epsilon})$, we need to obtain the corresponding asymptotic expansions of
the operators $\mathcal{K}^{*}_{D_{\epsilon}}$ and $\tilde{\phi}_{\epsilon}$.
The following lemmas can be found in [7].

######  Lemma 3.1.

For $\tilde{\phi}\in L^{2}(\partial D_{\epsilon})$, let
$\phi:=\tilde{\phi}\circ\Psi_{\epsilon}$. Then there exists a constant $C$
depending only on the $C^{2}$-norm of $\partial D$ and $\parallel
h\parallel_{C^{1}}$ such that

$\parallel\big{(}\mathcal{K}_{D_{\epsilon}}^{*}[\tilde{\phi}]\big{)}\circ\Psi_{\epsilon}-\mathcal{K}_{D}^{*}[{\phi}]-\epsilon\mathcal{K}_{D}^{(1)}[\phi]\parallel_{L^{2}(\partial
D)}\leq C\epsilon^{2}\parallel\phi\parallel_{L^{2}(\partial D)},$

with the operator $\mathcal{K}_{D}^{(1)}$ defined for any $\phi\in
L^{2}(\partial D)$ by

$\displaystyle\mathcal{K}_{D}^{(1)}[\phi](x)=\mathrm{p.v.}\int_{\partial
D}K_{1}(x,y)\phi(y)d\sigma(y),\ \ \ x\in\partial D,$ (3.13)

where

$\displaystyle K_{1}(x,y)$ $\displaystyle=-2\frac{\langle
x-y,\nu(x)\rangle\langle x-y,h(x)\nu(x)-h(y)\nu(y)\rangle}{|x-y|^{4}}$
$\displaystyle+\frac{\langle
h(x)\nu(x)-h(y)\nu(y),\nu(x)\rangle}{|x-y|^{2}}-\frac{\langle
x-y,\tau(x)h(x)\nu(x)+h^{\prime}(x)T(x)\rangle}{|x-y|^{2}}$
$\displaystyle+\frac{\langle
x-y,\nu(x)\rangle}{|x-y|^{2}}(h(x)\tau(x)-h(y)\nu(y)).$

Here, $\mathrm{p.v.}$ stands for the Cauchy principal value.

In fact, we can rewrite the operator $\mathcal{K}_{D}^{(1)}$ in terms of more
familiar operators as follows,

$\displaystyle\small\mathcal{K}_{D}^{(1)}[\phi]=-\frac{\partial}{\partial
T}\bigg{(}h\frac{\partial\mathcal{S}_{D}[\phi]}{\partial
T}\bigg{)}+\frac{\partial\mathcal{D}_{D}[h\phi]}{\partial\nu}+h\tau\mathcal{K}_{D}^{*}[\phi]-\mathcal{K}_{D}^{*}[h\tau\phi].$
(3.14)

######  Lemma 3.2.

Let $\tilde{\phi}_{\epsilon}=(\lambda
I-\mathcal{K}_{D_{\epsilon}}^{*})^{-1}[\tilde{\nu}\cdot\nabla H]$,
$\phi_{\epsilon}=\tilde{\phi}_{\epsilon}\circ\Psi_{\epsilon}$, and
$\phi=(\lambda I-\mathcal{K}^{*}_{D})^{-1}[\nu\cdot\nabla H]$. Then we have

$\displaystyle\parallel\phi_{\epsilon}-\phi-\epsilon\phi^{(1)}\parallel_{L^{2}(\partial
D)}\leq C\epsilon^{2}\parallel\phi\parallel_{L^{2}(\partial D)},$

where $C$ is a constant depending only on the $C^{2}$-norm of $\partial D$ and
$\parallel h\parallel_{C^{1}}$ and

$\displaystyle\phi^{(1)}=(\lambda
I-\mathcal{K}_{D}^{*})^{-1}\bigg{(}h\langle(\nabla^{2}H)\nu,\nu\rangle-h^{\prime}\langle\nabla
H,T\rangle+\mathcal{K}_{D}^{(1)}\phi\bigg{)}.$ (3.15)

Furthermore, for $\frac{\partial H}{\partial\tilde{\nu}}$ of (3.10), we can
further obtain by using (3.11), that

$\displaystyle\frac{\partial H}{\partial\tilde{\nu}}$ $\displaystyle=\nabla
H(\Psi(x))\cdot\tilde{\nu}(\Psi(x))=\frac{\partial H}{\partial\nu}(x)+\epsilon
G^{(1)}(x)+O(\epsilon^{2}),\ \ \ x\in\partial D$ (3.16)

where

$\displaystyle
G^{(1)}=h(x)\langle\nabla^{2}H(x)\nu(x),\nu(x)\rangle-h^{\prime}(x)\langle\nabla
H(x),T(x)\rangle.$ (3.17)

The asymptotic expansion of $\mathcal{F}(\partial D_{\epsilon})$ under small
perturbations of the boundary $\partial D$ in terms of the asymptotic
parameter $\epsilon$ is given in the following theorem. We would like to point
out that its proof is adapted from those in [25, 7].

######  Theorem 3.3.

Suppose that $\phi,\phi_{d}$ satisfy

$\displaystyle(\lambda I-\mathcal{K}_{D}^{*})[\phi](y)$
$\displaystyle=\frac{\partial H}{\partial\nu}\bigg{|}_{\partial D},$ (3.18)
$\displaystyle(\lambda I-\mathcal{K}_{D})[\phi_{d}](x,y)$
$\displaystyle=\Gamma(x,y),\ \ \ x\in\partial\Omega,\ y\in\partial D,$ (3.19)

the following asymptotic expansion holds:

$\displaystyle
u_{\epsilon}^{s}(x)|_{\partial\Omega}-u^{s}(x)|_{\partial\Omega}=\epsilon\langle
h(y),P(x,y)\rangle_{L^{2}(\partial D)}+O(\epsilon^{2}),\ \ \
x\in\partial\Omega,\ y\in\partial D,$

where

$\displaystyle P(x,y)=\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\bigg{(}\frac{\partial(H(y)+\mathcal{S}_{D}[\phi](y))}{\partial
T(y)}\bigg{)}+\bigg{(}\frac{\partial(\mathcal{D}_{D}[\phi_{d}](x,y)+\Gamma(x,y))}{\partial\nu(y)}\bigg{)}\phi(y).$
(3.20)

###### Proof.

From the form of the solution $(\ref{per solution})$, it holds that

$\mathcal{F}(\partial D_{\epsilon})=\int_{\partial
D_{\epsilon}}\Gamma(x,\tilde{y})\tilde{\phi}_{\epsilon}(\tilde{y})d\tilde{\sigma}(\tilde{y}),\
x\in\partial\Omega,$

where $\tilde{\phi}_{\epsilon}$ is solution to (3.10),
$d\tilde{\sigma}(\tilde{y})$ has an expansion as (3.12), and
$\Gamma(x,\tilde{y})$ has the Taylor expansion as follows

$\Gamma(x,y+\epsilon h(y)\nu(y))=\Gamma(x,y)+\epsilon
h(y)\frac{\partial\Gamma(x,y)}{\partial\nu(y)}+O(\epsilon^{2})\ \ \
y\in\partial D,\ x\in\partial\Omega.$

By virtue of Lemma 3.2, we can obtain that

$\displaystyle\mathcal{F}(\partial D_{\epsilon})-\mathcal{F}(\partial D)=$
$\displaystyle\epsilon\int_{\partial
D}\Gamma(x,y)\phi^{(1)}(y)d\sigma(y)+\epsilon\int_{\partial
D}\bigg{(}\frac{\partial\Gamma(x,y)}{\partial\nu(y)}-\tau(y)\Gamma(x,y)\bigg{)}\phi(y)h(y)d\sigma(y)$
$\displaystyle+O(\epsilon^{2}).$

Next, we first calculate the term $\int_{\partial
D}\Gamma(x,y)\phi^{(1)}(y)d\sigma(y)$. From $(\ref{asy density1})$ and (3.19)
we have

$\displaystyle\small\int_{\partial D}\Gamma(x,y)\phi^{(1)}(y)d\sigma(y)=$
$\displaystyle\int_{\partial D}(\lambda
I-\mathcal{K}_{D})\phi_{d}(x,y)\phi^{(1)}(y)d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial D}\phi_{d}(x,y)(\lambda
I-\mathcal{K}_{D}^{*})\phi^{(1)}(y)d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial
D}\phi_{d}(x,y)\bigg{(}h(y)\langle(\nabla^{2}H(y))\nu(y),\nu(y)\rangle-h^{\prime}(y)\frac{\partial
H(y)}{\partial T(y)}\bigg{)}d\sigma(y)$ $\displaystyle+\int_{\partial
D}\phi_{d}(x,y)\mathcal{K}_{D}^{(1)}[\phi](y)d\sigma(y)$ $\displaystyle=$
$\displaystyle:l_{1}+l_{2}.$

We treat $l_{1},l_{2}$ separately. Since

$\displaystyle\frac{\partial}{\partial
T(y)}\bigg{(}\phi_{d}(x,y)\frac{\partial H(y)}{\partial T(y)}\bigg{)}$
$\displaystyle=\frac{\partial\phi_{d}(x,y)}{\partial T(y)}\frac{\partial
H(y)}{\partial T(y)}+\phi_{d}(x,y)\langle(\nabla^{2}H(y))T(y),T(y)\rangle$
(3.21) $\displaystyle+\phi_{d}(x,y)\tau(y)\frac{\partial
H(y)}{\partial\nu(y)},$

and noting that $H(y)$ is a harmonic function and (3.21), $l_{1}$ has the
following form

$\displaystyle\small l_{1}=$ $\displaystyle\int_{\partial
D}\phi_{d}(x,y)\bigg{(}h(y)\langle(\nabla^{2}H(y))\nu(y),\nu(y)\rangle-h^{\prime}(y)\langle\nabla
H(y),T(y)\rangle\bigg{)}d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial
D}h(y)\bigg{(}\phi_{d}(x,y)\langle(\nabla^{2}H(y))\nu(y),\nu(y)\rangle+\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\frac{\partial H(y)}{\partial T(y)}$
$\displaystyle+\phi_{d}(x,y)\langle(\nabla^{2}H(y))T(y),T(y)\rangle+\phi_{d}(x,y)\tau(y)\frac{\partial
H(y)}{\partial\nu(y)}\bigg{)}d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial
D}h(y)\bigg{(}\frac{\partial\phi_{d}(x,y)}{\partial T(y)}\frac{\partial
H(y)}{\partial T(y)}+\phi_{d}(x,y)\Delta
H(y)+\phi_{d}(x,y)\tau(y)\frac{\partial
H(y)}{\partial\nu(y)}\bigg{)}d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial
D}h(y)\bigg{(}\frac{\partial\phi_{d}(x,y)}{\partial T(y)}\frac{\partial
H(y)}{\partial T(y)}+\phi_{d}(x,y)\tau(y)\frac{\partial
H(y)}{\partial\nu(y)}\bigg{)}d\sigma(y).$

Next, by (3.14), we have

$\displaystyle l_{2}=$ $\displaystyle\int_{\partial
D}\phi_{d}(x,y)\mathcal{K}_{D}^{(1)}[\phi](y)d\sigma(y)$ $\displaystyle=$
$\displaystyle\int_{\partial D}\phi_{d}(x,y)\bigg{[}-\frac{\partial}{\partial
T(y)}\bigg{(}h(y)\frac{\partial\mathcal{S}_{D}[\phi](y)}{\partial
T(y)}\bigg{)}+\frac{\partial}{\partial\nu(y)}\mathcal{D}_{D}[\phi
h](y)+h(y)\tau(y)\mathcal{K}_{D}^{*}[\phi](y)$
$\displaystyle-\mathcal{K}_{D}^{*}[h\tau\phi](y)\bigg{]}d\sigma(y)$
$\displaystyle=$ $\displaystyle\int_{\partial
D}h(y)\bigg{[}\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\frac{\partial(\mathcal{S}_{D}\phi(y))}{\partial
T(y)}+\frac{\partial}{\partial\nu(y)}\mathcal{D}_{D}[\phi_{d}](x,y)\phi(y)+\phi_{d}(x,y)\tau(y)\mathcal{K}_{D}^{*}[\phi](y)$
$\displaystyle-\tau(y)\phi(y)\mathcal{K}_{D}[\phi_{d}](x,y)\bigg{]}d\sigma(y),$

from which, together with the use of Proposition 4.1 in [7], we can further
show that

$\int_{\partial
D}\phi_{d}(x,y)\frac{\partial\mathcal{D}_{D}[h\phi](y)}{\partial\nu(y)}d\sigma(y)=\int_{\partial
D}\frac{\partial\mathcal{D}_{D}[\phi_{d}](x,y)}{\partial\nu(y)}h(y)\phi(y)d\sigma(y).$

Thus we obtain that

$\displaystyle\mathcal{F}(\partial D_{\epsilon})-\mathcal{F}(\partial D)$
$\displaystyle=\int_{\partial
D}h(y)\bigg{[}\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\bigg{(}\frac{\partial(\mathcal{S}_{D}[\phi])(y)}{\partial
T(y)}\bigg{)}+\frac{\partial(\mathcal{D}_{D}[\phi_{d}])(x,y)}{\partial\nu(y)}\phi(y)$
$\displaystyle+\tau(y)\phi_{d}(x,y)\mathcal{K}^{*}_{D}[\phi](y)-\tau(y)\mathcal{K}_{D}[\phi_{d}](x,y)\phi(y)$
$\displaystyle+\phi_{d}(x,y)\tau(y)\frac{\partial
H(y)}{\partial\nu(y)}+\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\frac{\partial H(y)}{\partial T(y)}$
$\displaystyle+\bigg{(}\frac{\partial\Gamma(x,y)}{\partial\nu(y)}-\tau(y)\Gamma(x,y)\bigg{)}\phi(y)\bigg{]}d\sigma(y)+O(\epsilon^{2}).$

From (3.18) and (3.19), we know that

$\displaystyle\mathcal{K}^{*}_{D}[\phi](y)+\frac{\partial
H(y)}{\partial\nu(y)}=\lambda\phi(y),\ \ \
\mathcal{K}_{D}[\phi_{d}](x,y)+\Gamma(x,y)=\lambda\phi_{d}(x,y).$

Thus it follows that

$\displaystyle\mathcal{F}(\partial D_{\epsilon})-\mathcal{F}(\partial D)$
$\displaystyle=\epsilon\int_{\partial
D}h(y)\bigg{[}\frac{\partial\phi_{d}(x,y)}{\partial
T(y)}\bigg{(}\frac{\partial(H(y)+\mathcal{S}_{D}[\phi](y))}{\partial
T(y)}\bigg{)}$
$\displaystyle+\bigg{(}\frac{\partial(\Gamma(x,y)+\mathcal{D}_{D}[\phi_{d}](x,y))}{\partial\nu(y)}\bigg{)}\phi(y)\bigg{]}d\sigma(y)+O(\epsilon^{2}),$

which readily completes the proof. ∎

In order to investigate variations in the measurement resulting from
variations in the shape of the underlying object, we introduce the following
definition of the shape sensitivity functional.

######  Definition 3.1.

The shape sensitivity functional for the far-field measurement
$u^{s}|_{\partial\Omega}$ with respect to the shape of $\partial D$ is defined
as

$\displaystyle SSF(\partial D):=\lim_{\epsilon\rightarrow
0}\frac{u_{\epsilon}^{s}|_{\partial\Omega}-u^{s}|_{\partial\Omega}}{\epsilon}.$

######  Remark 3.1.

From Definition 3.1, it is easy to see that the shape sensitivity function is
actually the shape derivative of the forward operator $\mathcal{F}(\partial
D)$ (cf. [7]). Furthermore, by using Theorem 3.3, the shape sensitivity
function can be rewritten as

$\displaystyle SSF(\partial D)=\langle h,P\rangle_{L^{2}(\partial D)},$

where $P$ is defined in (3.20).

###  3.2 Shape sensitivity analysis and plasmon resonance

In this subsection, we shall derive the spectral representation of the shape
sensitivity functional in this subsection. It indicates that, when the plasmon
resonance occurs, the shape sensitivity functional will be amplified and
exhibits a large peak. Hence, the plasmon resonance can be used to
significantly increase the sensitivity of the far-field measurement with
respect to the shape of the underlying domain. For simplicity, in the
subsequent spectral analysis, we always exclude the essential spectrum $0$
from the spectrum set of NP operator and assume the eigenvalues is simple.

First, we derive the asymptotic formulas for the eigenvalues and
eigenfunctions of the NP operator with respect to the asymptotic parameter
$\epsilon$.

######  Lemma 3.4.

Suppose that $\\{\lambda_{j,\epsilon},\varphi_{j,\epsilon}\\},j=0,1,2,\cdots,$
are the perturbed eigenvalues and eigenfunctions of
$\mathcal{K}^{*}_{D_{\epsilon}}$. Then we have as $\epsilon\rightarrow 0$ that

$\displaystyle\lambda_{j,\epsilon}=\lambda_{j}+\epsilon\lambda_{j,\epsilon}^{(1)}+O(\epsilon^{2}),\quad\varphi_{j,\epsilon}=\varphi_{j}+\epsilon\varphi_{j,\epsilon}^{(1)}+O(\epsilon^{2}),$

where

$\displaystyle\lambda_{j,\epsilon}^{(1)}=\langle\mathcal{K}_{D}^{(1)}\varphi_{j},\varphi_{j}\rangle_{\mathcal{H}^{*}(\partial
D)},\quad\varphi_{j,\epsilon}^{(1)}=\sum_{j\neq
l}\frac{\langle\mathcal{K}_{D}^{(1)}\varphi_{j},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}}{\lambda_{j}-\lambda_{l}}\varphi_{l},$

and $\mathcal{K}_{D}^{(1)}$ is defined by (3.13) or (3.14).

###### Proof.

Since
$\mathcal{K}_{D_{\epsilon}}^{*}[\varphi_{j,\epsilon}]=\lambda_{j,\epsilon}\varphi_{j,\epsilon}$,
and by Lemma 3.1, we can obtain the following result

$(\mathcal{K}_{D}^{*}+\epsilon\mathcal{K}_{D}^{(1)}+O(\epsilon^{2}))(\varphi_{j}+\epsilon\varphi_{j,\epsilon}^{(1)}+O(\epsilon^{2}))=(\lambda_{j}+\epsilon\lambda_{j,\epsilon}^{(1)}+O(\epsilon^{2}))(\varphi_{j}+\epsilon\varphi_{j,\epsilon}^{(1)}+O(\epsilon^{2})),$

which implies that

$\displaystyle\mathcal{K}_{D}^{*}\varphi_{j}$
$\displaystyle=\lambda_{j}\varphi_{j},$
$\displaystyle\mathcal{K}_{D}^{(1)}\varphi_{j}+\mathcal{K}_{D}^{*}\varphi_{j,\epsilon}^{(1)}$
$\displaystyle=\lambda_{j,\epsilon}^{(1)}\varphi_{j}+\lambda_{j}\varphi_{j,\epsilon}^{(1)}.$
(3.22)

By inner producting $\varphi_{l}$ in both sides of (3.22), we further have

$\displaystyle\langle\mathcal{K}_{D}^{(1)}\varphi_{j},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}+\langle\mathcal{K}_{D}^{*}\varphi_{j,\epsilon}^{(1)},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}=\lambda_{j,\epsilon}^{(1)}\langle\varphi_{j},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}+\lambda_{j}\langle\varphi_{j,\epsilon}^{(1)},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}.$

Since the operator $\mathcal{K}_{D}^{*}$ is self-adjoint, we have

$\langle\mathcal{K}_{D}^{*}\varphi_{j,\epsilon}^{(1)},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}=\langle\varphi_{j,\epsilon}^{(1)},\mathcal{K}_{D}^{*}\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}=\lambda_{l}\langle\varphi_{j,\epsilon}^{(1)},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)},$

and thus

$\lambda_{j,\epsilon}^{(1)}=\langle\mathcal{K}_{D}^{(1)}\varphi_{j},\varphi_{j}\rangle_{\mathcal{H}^{*}(\partial
D)}.$

From (3.22), we have

$(\lambda_{j}I-\mathcal{K}_{D}^{*})\varphi_{j,\epsilon}^{(1)}=\mathcal{K}_{D}^{(1)}\varphi_{j}-\lambda_{j,\epsilon}^{(1)}\varphi_{j}.$

Then it is straightforward to see that

$\varphi_{j,\epsilon}^{(1)}=\sum_{j\neq
l}\frac{\langle\mathcal{K}_{D}^{(1)}\varphi_{j},\varphi_{l}\rangle_{\mathcal{H}^{*}(\partial
D)}}{\lambda_{j}-\lambda_{l}}\varphi_{l}.$

The proof is complete. ∎

Let $\text{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))$ be the distance of
$\lambda$ to the spectrum set $\sigma(\mathcal{K}_{D}^{*})$. Next, we are
ready to state our main result in the following theorem.

######  Theorem 3.5.

As $\epsilon\rightarrow 0$, the perturbation
$u^{s}_{\epsilon}|_{\partial\Omega}-u|_{\partial\Omega}$ has the following
spectral expansion:

$\displaystyle
u^{s}_{\epsilon}(x)-u^{s}(x)=\epsilon\mathcal{T}(x){+O\left(\frac{\epsilon^{2}}{\mathrm{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))^{2}}\right)},\
\ \ \ x\in\partial\Omega,$

with

$\displaystyle\mathcal{T}(x)=$
$\displaystyle\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j,\epsilon}^{(1)}\rangle+\langle
G^{(1)},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{Z}(x)}{\lambda-\lambda_{j}}$ (3.23)
$\displaystyle+\sum_{j=1}^{\infty}\frac{\lambda_{j,\epsilon}^{(1)}\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)}{(\lambda-\lambda_{j})^{2}},$

where $G^{(1)}$ is defined by (3.17), $\langle\cdot,\cdot\rangle$ signifies
$\langle\cdot,\cdot\rangle_{\mathcal{H}^{*}(\partial D)}$, and

$\displaystyle\small\mathcal{Z}(x)=\int_{\partial
D}\left(\frac{\partial\Gamma(x,y)}{\partial\nu(y)}h(y)\varphi_{j}(y)+\Gamma(x,y)\varphi_{j,\epsilon}^{(1)}(y)-\tau(y)h(y)\Gamma(x,y)\varphi_{j}(y)\right)d\sigma(y).$

###### Proof.

According to (2.3), $\mathcal{F}(\partial D)$ can be decomposed as follow

$\displaystyle\mathcal{F}(\partial D)$ $\displaystyle=\int_{\partial
D}\Gamma(x,y)(\lambda I-\mathcal{K}_{D}^{*})^{-1}\bigg{[}\frac{\partial
H}{\partial\nu}\bigg{]}(y)d\sigma(y)$
$\displaystyle=\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle}{\lambda-\lambda_{j}}\int_{\partial
D}\Gamma(x,y)\varphi_{j}(y)d\sigma(y)$
$\displaystyle=\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)}{\lambda-\lambda_{j}}.$

In a similar manner, $\mathcal{F}(\partial D_{\epsilon})$ can be decomposed as

$\displaystyle\small\mathcal{F}(\partial D_{\epsilon})$
$\displaystyle=\int_{\partial D_{\epsilon}}\Gamma(x,\tilde{y})(\lambda
I-\mathcal{K}_{D_{\epsilon}}^{*})^{-1}\bigg{[}\frac{\partial
H}{\partial\tilde{\nu}}\bigg{]}(\tilde{y})d\tilde{\sigma}(\tilde{y})$
$\displaystyle=\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\tilde{\nu}},\varphi_{j,\epsilon}\rangle\mathcal{S}_{D_{\epsilon}}[\varphi_{j,\epsilon}](x)}{\lambda-\lambda_{j,\epsilon}}.$

Furthermore, using the perturbation of eigenvalues and eigenfunctions in Lemma
3.4, we can show that

$\displaystyle\small\mathcal{S}_{D_{\epsilon}}[\varphi_{j,\epsilon}](x)=\int_{\partial
D}\Gamma(x,y)\varphi_{j}(y)d\sigma(y)+\epsilon\mathcal{Z}(x)+O(\epsilon^{2}),$
(3.24)

where

$\displaystyle\small\mathcal{Z}(x)=\int_{\partial
D}\bigg{(}\frac{\partial\Gamma(x,y)}{\partial\nu(y)}h(y)\varphi_{j}(y)+\Gamma(x,y)\varphi_{j,\epsilon}^{(1)}(y)-\tau(y)h(y)\Gamma(x,y)\varphi_{j}(y)\bigg{)}d\sigma(y).$

From $(\ref{per H/v})$, $(\ref{per Sd})$ and Lemma 3.4, we deduce that

$\displaystyle\mathcal{F}(\partial D_{\epsilon})=$
$\displaystyle\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\tilde{\nu}},\varphi_{j,\epsilon}\rangle\mathcal{S}_{D_{\epsilon}}[\varphi_{j,\epsilon}](x)}{\lambda-\lambda_{j,\epsilon}}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu}+\epsilon
G^{(1)}+O(\epsilon^{2}),\varphi_{j}+\epsilon\varphi_{j,\epsilon}^{(1)}+O(\epsilon^{2})\rangle(\mathcal{S}_{D}[\varphi_{j}](x)+\epsilon\mathcal{Z}(x)+O(\epsilon^{2}))}{\lambda-\lambda_{j}-\epsilon\lambda_{j,\epsilon}^{(1)}+{O(\frac{\epsilon^{2}}{\lambda-\lambda_{j}})}}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\epsilon\bigg{(}\langle\frac{\partial
H}{\partial\nu},\varphi_{j,\epsilon}^{(1)}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle
G^{(1)},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}]+\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{Z}(x)\bigg{)}+O(\epsilon^{2})}{(\lambda-\lambda_{j})(1-\frac{\epsilon\lambda_{j,\epsilon}^{(1)}}{\lambda-\lambda_{j}}+{O(\frac{\epsilon^{2}}{\lambda-\lambda_{j}}))}}$
$\displaystyle=$
$\displaystyle\bigg{(}\sum_{j=1}^{\infty}\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\epsilon\big{(}\langle\frac{\partial
H}{\partial\nu},\varphi_{j,\epsilon}^{(1)}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle
G^{(1)},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)$
$\displaystyle+\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{Z}(x)\big{)}+O\bigg{(}\epsilon^{2}\bigg{)}\bigg{)}\cdot\bigg{(}\sum_{k=0}^{\infty}\frac{1}{(\lambda-\lambda_{j})}\big{(}\frac{\epsilon\lambda_{j,\epsilon}^{(1)}}{\lambda-\lambda_{j}}+{O(\frac{\epsilon^{2}}{\lambda-\lambda_{j}})\big{)}^{k}}\bigg{)}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)}{\lambda-\lambda_{j}}+\epsilon\bigg{(}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j,\epsilon}^{(1)}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle
G^{(1)},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{Z}(x)}{\lambda-\lambda_{j}}$
$\displaystyle+\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)\lambda_{j,\epsilon}^{(1)}}{(\lambda-\lambda_{j})^{2}}\bigg{)}{+O\left(\frac{\epsilon^{2}}{\text{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))^{2}}\right)}.$

Therefore it follows that

$\displaystyle\mathcal{F}(\partial D_{\epsilon})-\mathcal{F}(\partial D)$
$\displaystyle=\epsilon\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j,\epsilon}^{(1)}\rangle+\langle
G^{(1)},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)+\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{Z}(x)}{\lambda-\lambda_{j}}$
$\displaystyle+\epsilon\sum_{j=1}^{\infty}\frac{\langle\frac{\partial
H}{\partial\nu},\varphi_{j}\rangle\mathcal{S}_{D}[\varphi_{j}](x)\lambda_{j,\epsilon}^{(1)}}{(\lambda-\lambda_{j})^{2}}{+O\left(\frac{\epsilon^{2}}{\text{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))^{2}}\right)}.$

The proof is complete. ∎

######  Remark 3.2.

From (3.23), one readily sees that the $j$-th mode in the expansion formula
contains both $\mathcal{O}\left(\frac{1}{\lambda-\lambda_{j}}\right)$ and
$\mathcal{O}\left(\frac{1}{(\lambda-\lambda_{j})^{2}}\right)$. Thus, for the
sufficiently small loss ($\mathfrak{Im}(\varepsilon_{D})\rightarrow 0$), the
$j$-th mode will exhibit a large peak if the plasmon resonance condition (2.8)
is fulfilled.

From Theorem 3.5, we can straightforwardly obtain the spectral representation
formula of the shape sensitivity function as follows.

######  Corollary 3.6.

If
$\epsilon=o\left(\mathrm{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))^{2}\right)$,
(as $\mathrm{dist}(\lambda,\sigma(\mathcal{K}_{D}^{*}))\rightarrow 0$), the
shape sensitivity function can be represented as

$\displaystyle SSF(\partial D)=\mathcal{T}(x),$

where $\mathcal{T}(x)$ is defined by (3.23).

######  Remark 3.3.

Similar to Remark 3.2, one sees that when
$\lambda=\frac{\varepsilon_{D}+\varepsilon_{m}}{2(\varepsilon_{D}-\varepsilon_{m})}$
is very close to an eigenvalue $\lambda_{j}$ of the NP operator, i.e., the
plasmon resonance occurs, the shape sensitivity function $SSF(\partial D)$
will be amplified dramatically for the $j$-th mode. Hence, the plasmon
resonance can improve the sensitivity of the shape construction. Moreover,
from the Drude’s model (2.6), we can compute the plasmon resonance frequency
for improving the sensitivity.

######  Remark 3.4.

From a numerical point of view, the discrete version of the shape sensitivity
functional is the shape sensitivity matrix. The sensitivity of the shape
reconstruction can be analysed by the singular value decomposition of the
sensitivity matrix. In Section 5, the numerical results shall show that when
the plasmon resonance occurs, the singular values of the sensitivity matrix
increase dramatically. That is the plasmon resonance technology can enhance
the stability of the Gauss-Newton iteration algorithm that we use therein (see
Section 4).

##  4 Tikhonov regularization and Laplace approximation

In this section, we discuss the numerical issues for the shape reconstruction.
First, a numerical implementation requires a parametrization of the boundary
$\partial D$. Here we assume that $\partial D$ is starlike boundary curve with
respect to the origin, i.e. there exists $q\in C^{2}[0,2\pi]$ such that

$\displaystyle\partial D=\\{\vec{q}(t)=q(t)\bigg{(}\begin{array}[]{c}\cos t\\\
\sin t\end{array}\bigg{)},\ t\in[0,2\pi]\\}.$

The admissible set $X=\\{q\in C^{2}([0,2\pi]):q>0\\}$, and the forward
operator maps $X$ into $H^{\frac{1}{2}}(\partial\Omega)$. Without change of
notation we write $\mathcal{F}(q)=u^{s}(x)|_{\partial\Omega}$.

In order to overcome the ill-posedness of the inverse problem, we apply the
Tikhonov regularization method to tackle it. The corresponding variational
functional is given as follows,

$\displaystyle
J[q]:=\frac{1}{2}\parallel\mathcal{F}(q)-u^{s,\delta}\parallel_{2}^{2}+\frac{\mu}{2}\parallel
q\parallel_{L^{2}[0,2\pi]}^{2},$ (4.25)

where,
$u^{s,\delta}=(u_{1}^{s,\delta},u_{2}^{s,\delta},\cdots,u_{n}^{s,\delta})$
signfies the discrete measurement data on $\partial\Omega$. Here,
$\|\cdot\|_{2}$ denotes the $2$-norm in $\mathbb{R}^{n}$, and $\mu$ is the
regularization parameter. Moreover, the measurement data $u^{s,\delta}$ and
the exact data $u^{s}$ satisfies $\|u^{s,\delta}-u^{s}\|_{2}\leq\delta$.

There are many iterative algorithms to solve the variational problem (4.25).
In this paper, we apply the Levenberg-Marquardt method [23, 29] to find the
minimize of $J[q]$, which is essentially a variant of the Gauss-Newton
iteration. Assuming that $q^{*}$ is an approximation of $q$, then the
nonlinear mapping $\mathcal{F}$ in (4.25) can be replaced approximatively by
its linearization around $q^{*}$, i.e.

$\mathcal{F}(q)\approx\mathcal{F}(q^{*})+\mathcal{F}^{\prime}(q^{*})(q-q^{*}).$

The nonlinear inverse problem $\mathcal{F}(q)=u^{s,\delta}$ can then be
converted to a linear inverse problem

$\mathcal{F}^{\prime}(q^{*})(q-q^{*})=u^{s,\delta}-\mathcal{F}(q^{*}).$

Thus, minimizing (4.25) can easily be seen to minimizing

$\displaystyle J[q]:=\frac{1}{2}\parallel\mathcal{F}^{\prime}(q^{*})\delta
q-(u^{s,\delta}-\mathcal{F}(q^{*}))\parallel_{2}^{2}+\frac{\mu}{2}\parallel\delta
q\parallel_{L^{2}[0,2\pi]}^{2},$

where, $\delta q=q-q^{*}$. We formulate the Levenberg-Marquardt iteration as
Algorithm 1.

Algorithm 1 Levenberg-Marquardt algorithm

Input: The regularization parameter $\mu$, the noise level $\delta$ and the
maximum iterations $M$
Output: $q_{k+1}$.
Require:
1: Initialize the solution $q_{0}$;
2: Settle the forward problem and obtain the additional data $u^{s,\delta}$;
3: While $k<M$
4: Compute the forward problem and set
$F_{k}=u^{s,\delta}-\mathcal{F}(q_{k});$
5: Compute the Jacobian matrix $G=\mathcal{F}^{\prime}(q_{k})$;
6: Compute $\delta q_{k}=(G^{T}G+\mu I)^{-1}(G^{T}F_{k});$
7: Update the solution: $q_{k+1}=q_{k}+\delta q_{k};$
8: $k\leftarrow k+1$;
9: Terminate if $\parallel\delta q_{k}\parallel<eps$
10: end

We also decipher the inverse problem from a Bayesian perspective so that we
can capture more statistical information about the solution. From the
classical Bayesian theory, and noticing that the observation error $\xi$ is
assumed to be an independent and identically distributed Gauss random vector
with mean zero and the covariance matrix $B=\delta^{2}I$ (here $I$ is the unit
matrix), we can write the minimization functional as

$\displaystyle J[q]$
$\displaystyle\propto\frac{1}{2\delta^{2}}\parallel\mathcal{F}(q)-u^{\delta}\parallel_{2}^{2}+\frac{\mu}{2\delta^{2}}\parallel
q\parallel_{L^{2}[0,2\pi]}^{2}$
$\displaystyle=\frac{1}{2}\parallel\mathcal{F}(q)-u^{\delta}\parallel_{B}^{2}+\frac{\mu}{2\delta^{2}}\parallel
q\parallel_{L^{2}[0,2\pi]}^{2}$ $\displaystyle=:J_{B}(q),$

where $\parallel\cdot\parallel_{B}$ is a covariance weighted norm given by
$\parallel\cdot\parallel_{B}=\parallel\delta^{-1}I\cdot\parallel_{2}$, and the
minimizer of $J_{B}$ defines the maximum a posteriori estimator

$\displaystyle q_{\mathrm{MAP}}=\arg\min_{q}J_{B}(q).$ (4.26)

The Laplace approximation replaces the complicated posterior with a normal
distribution located at the maximum a posteriori value $q_{\mathrm{MAP}}$. Its
essence is a linearization around the $\mathrm{MAP}$ point $q_{\mathrm{MAP}}$
(cf. [48]). It consists of approximating the posterior measure (or
distribution) by $\omega\approx N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$, where

$\displaystyle
C_{\mathrm{MAP}}=(J^{\prime\prime}_{B})^{-1}=(\frac{\mu}{\delta^{2}}I+\frac{1}{\delta^{2}}G^{T}G)^{-1},$
(4.27)

and $G$ is the Jacobian matrix of the forward operator $\mathcal{F}$ at the
point $q$. Notice that the covariance formula (4.27) only uses the first order
derivatives of $\mathcal{F}$. The implementation of the Laplace approximation
is presented in the following algorithm [31].

Algorithm 2 Laplace approximation (LA) for sampling.

1: Compute $q_{\mathrm{MAP}}$ from $(\ref{MAP})$ by using Algorithm 1, and
$C_{\mathrm{MAP}}$ from $(\ref{C_map})$, respectively;
2: Compute the Cholesky factor $L$ of $C_{\mathrm{MAP}}$, i.e.,
$C_{\mathrm{MAP}}=LL^{T}$;
3: For $j=\\{1,...,N_{e}\\}$, generate $q^{j}=q_{\mathrm{MAP}}+L^{T}z^{j}$,
where $z^{j}\sim N(0,I)$.

In this algorithm, samples generated by algorithm 2 are drawn from
$N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$, and so the ensemble of $N_{e}$
realizations $\\{q_{j}\\}_{j=1}^{N_{e}}$ provides an approximation to
$N(q_{\mathrm{MAP}},C_{\mathrm{MAP}})$ and hence the posterior. Finally, we
use the mean of the sample $\bar{q}=\frac{1}{N_{e}}\sum_{j=1}^{N_{e}}q_{j}$ as
an approximation of $q_{\mathrm{MAP}}$, where the convergence of $\bar{q}$
follows from the strong law of large numbers. From the classical Gaussian
statistic theory, we find that $\bar{q}$ is consistent and the best unbiased
estimate of $q_{\mathrm{MAP}}$.

##  5  Numerical results and discussions

In this section, we present several numerical examples to illustrate the
salient and promising features of the proposed reconstruction scheme.

In all of our numerical examples, the measurement boundary curve
$\partial\Omega$ is given by the circle of radius 3 and centred at the origin,
that is $\partial\Omega=\\{3(\cos t,\sin t),0\leq t\leq 2\pi\\}$. Set the
incident field $H(x)=x_{1}$. We solve the direct problem through the
$\mathrm{Nystr\ddot{o}m}$ method [33], which is discretized with $n=80$ grid
points. Furthermore, in an effort to avoid committing an inverse crime, the
number of collocation points for obtaining the synthetic data was chosen to be
different from the number of collocation points within the inverse solver. In
addition, we approximate the radial function $q(t)$ for unknown interior
boundary curve $\partial D$ by the trigonometric series

$q(t)\approx\sum_{k=0}^{m}a_{k}\cos kt+\sum_{k=1}^{m}b_{k}\sin kt,\ \ \ 0\leq
t\leq 2\pi,$

where $m\in\mathbb{N}$ and the vector
$q=(a_{0},...,a_{m},b_{1},...,b_{m})\in\mathbb{R}^{2m+1}$.

In the iterative process, we choose a circle as the initial guess, which
contains the inclusion $D$. A finite difference method is used to calculate
the Jacobian matrix $G$, and the maximum number of iteration steps is 100. The
number of samples $N_{e}$ is 10000, and we use the following stopping rule

$E_{k}=\parallel q_{k}-q_{k-1}\parallel_{L^{2}}\leq 10^{-5}.$

The noisy measured data are generated by

$u^{s,\delta}=u^{s}(x)+\delta\xi,\ \ x\in\partial\Omega,$

where $u^{s}(x)$ is the exact data, $\delta$ indicates the noise level, and
$\xi$ is the Gaussian random vector with a zero mean and unit standard
deviation.

The accuracy of the approximate solution $\tilde{q}_{N_{e}}$ by the LA
algorithm is characterised by comparing to the exact solution $q(x)$ via the
relative error

$e_{\gamma}=\frac{\parallel\tilde{q}_{N_{e}}-q(x)\parallel_{L^{2}(\partial
D)}}{\parallel q(x)\parallel_{L^{2}(\partial D)}}.$

In particular, we specifically take the permittivity $\varepsilon_{m}=24.8F/m$
(such as alcohol) in the first example, and in all of the examples we set
$\varepsilon_{m}=\varepsilon_{0}$. Moreover, except for Example 5.1, the
spectrum of $\mathcal{K}_{D}^{*}$ needs to calculated numerically, from which
one can then compute the corresponding plasmon resonance frequency. We use the
$\mathrm{Nystr\ddot{o}m}$ method for spectral calculations. However, we would
like to emphasise that in practice, the plasmon resonance frequency can be
measured by special apparatus, when one can observe significant enhancements
in the electric field due to the resonance.

######  Example 5.1.

In this example, we consider the reconstruction of a circle object with

$q(t)=0.5,\ \ \ 0\leq t\leq 2\pi.$

In the iteration algorithm, we choose a circle of radius 0.6 as the initial
guess.

First, we investigate the influence of different $\lambda$ on the
reconstruction effect. The numerical results for Example 5.1 with $\lambda$
taking different values are shown in Figure 5.1 (a). It is well known that the
eigenvalues of $\mathcal{K}_{D}^{*}$ are $\\{0,\frac{1}{2}\\}$ for $D$ being a
disk. If setting $\lambda_{1}=0$, from Section 2.2, we can see that the
Fröhlich condition is satisfied based on the Drude’s model for metal nanorod
and the plasmon resonance occurs (though our theoretical analysis does not
contain the essential spectrum case). In particular, we take the resonance
frequency $\omega_{1}=3.63\cdot 10^{14}Hz$, and according to the experimental
data of gold nanorod in [45] (the data coincide with the Drude model), we can
get $\mathfrak{Re}(\varepsilon_{D}(\omega_{1}))=-24.8$, and the imaginary part
$\mathfrak{Im}(\varepsilon_{D}(\omega_{1}))=0.797$. Furthermore, by (2.5), it
implies that $\lambda\approx 0-8\cdot 10^{-3}i$.

Figure 5.1: Reconstruction the shape for Example 5.1 with 1% noise data and
the regularization parameter $\mu=0.01$. (a) Reconstructions with different
$\lambda$ as well as the exact solution, (b) the relative error $e_{\gamma}$
versus the iterations step $k$, associated to different values of $\lambda$.

Next, we take $\omega_{2}=2.30\cdot 10^{14}Hz$ and based on the experimental
data in [45], we can calculate $\lambda=0.244-6.3\cdot 10^{-3}i$. It is
obvious that when taking $\omega_{2}$, the real part of $\lambda(\omega_{2})$
satisfies $|\mathfrak{Re}(\lambda(\omega_{2}))<\frac{1}{2}|$, which is only an
approximate resonance frequency. As seen in Figure 5.1 (a), a perfect
reconstruction result is obtained when taking the resonance frequency
$\omega_{1}$. When $\lambda=0.244-6.3\cdot 10^{-3}i$, although the
reconstructed result is good, it is still poor compared to the case with the
resonant frequency $\omega_{1}$. The case with $\lambda=-0.98$ indicates a
normal material, i.e. the permittivity $\varepsilon_{D}=8$ is positive
constant, and the reconstruction effect is poorer. The plot in Figure 5.1 (b)
illustrates the relative error $e_{\gamma}$ versus the iteration number with
different $\lambda$. It can be seen that the relative error can quickly reach
the convergence and remain small when $\omega_{1}$ is the plasmon resonance
frequency. When $\lambda=0.244-6.3\cdot 10^{-3}i$, the relative error also
converges quickly, but is larger than that at the plasmon resonance frequency.
In particular, the relative error increases rapidly even up to 27%, and the
inversion is rather imprecise with $\lambda=-0.98$.

Figure 5.2: Reconstruction of the shape for Example 5.1 with 5% noise data and
the regularization parameter $\mu=0.05$. (a) Reconstruction under different
$\lambda$ as well as the exact solution, (b) the accuracy error $e_{\gamma}$
of the number of the iterations $k$, for different values of $\lambda$.

The numerical results for Example 5.1, which has $5\%$ noise in the data, are
shown in Figure 5.2(a) for different $\lambda$. Compared to the $1\%$ noise in
the data shown in Figure 5.1(a), it is clear that as the noise level
increases, the reconstruction effect becomes worse. It is worth noting that
when $\omega$ is a plasmon resonance frequency, the reconstruction is very
good even at larger error levels. For $\lambda=-0.98$ (general materials) the
inversion result deviates from the exact solution quickly as the error level
increases. We also list the relative errors for different noise levels in
Table 1.

Table 1:  Numerical results of Example 5.1 for $e_{\gamma}$ associated with different $\lambda$ and noise level $\delta$ $\lambda$ | $\delta=0.01$ | $\delta=0.05$
---|---|---
$0-8\cdot 10^{-3}i$ | $7.2\cdot 10^{-5}$ | $9.4\cdot 10^{-4}$
$0.24-6.3\cdot 10^{-3}i$ | $0.106$ | $0.574$
$-0.98$ | $0.276$ | $1.546$
|  |

######  Example 5.2.

In this example, we consider the reconstruction of a bean contour with a
radial function,

$q(t)=\frac{4/5+18/25\cos t+3/25\sin 2t}{1+7/10\cos t},\ \ \ 0\leq t\leq
2\pi.$

We choose a circle of radius 0.72 as the initial guess.

Figure 5.3: Reconstruction of the shape in Example 5.2 with 1% noise data and
the regularization parameter $\mu=0.1$. (a) Reconstruction with different
$\lambda$ as well as the exact solution, (b) the relative error $e_{\gamma}$
versus the iteration step $k$ associated with different values of $\lambda$.

The numerical results for Example 5.2 with $\lambda$ taking different values
are shown in Figure 5.3 (a) and the relative error $e_{\gamma}$ versus the
iteration number with different $\lambda$ in Figure 5.3 (b). The real part of
$\lambda=0.16-10^{-6}i$ is the eigenvalue of $\mathcal{K}^{*}_{D}$ in Example
5.2. It represents very accurate approximations of the LA solution to the
exact solution, and in the iteration algorithm, the relative error quickly
reaches a stable state and remains small. In addition, when $\lambda=2$, the
reconstruction is bad. The results show that the relative error becomes more
oscillating as the noise level $\delta$ increases; see Figure 5.4 and Table 2
for $e_{\gamma}$ with different $\lambda$ and $\delta$.

Figure 5.4: Reconstruction of the shape in Example 5.2 with 5% noise data and the regularization parameter $\mu=0.1$. (a) Reconstruction with different $\lambda$ as well as the exact solution, (b) the relative error $e_{\gamma}$ versus iteration steps $k$ associated with different values of $\lambda$. Table 2:  Numerical results of Example 5.2 for $e_{\gamma}$ associated with different $\lambda$ and $\delta$ $\lambda$ | $\delta=0.01$ | $\delta=0.05$
---|---|---
$0.16-10^{-6}i$ | $0.0081$ | $0.1237$
$0.25-10^{-6}i$ | $0.0278$ | $0.2810$
$2$ | $0.5304$ | $0.5804$
|  |

######  Example 5.3.

In this example, we consider that the inclusion is a peanut, and the polar
radius of the peanut is parameterized by

$q(t)=\sqrt{\cos^{2}t+0.26\sin^{2}(t+0.5)},\ \ \ 0\leq t\leq 2\pi.$

In the iteration, we choose a circle of radius 0.78 as the initial guess.

From [19], it is known that the singular value decomposition of the
sensitivity matrix plays a key role in the uncertainty estimation. Let the
singular value decomposition (SVD) of the sensitivity matrix (Jacobian matrix)
$G$ of the forward operator at the true solution $q_{true}$ be denoted as

$\displaystyle\begin{split}G(q_{true})=U\bigg{[}\begin{array}[]{c}\mathrm{\Lambda}\\\
0\end{array}\bigg{]}V^{T},\end{split}$

where $U$ is an $n\times n$ orthogonal matrix, i.e. $U^{T}U=UU^{T}=I_{n}$,
with $U_{1}$ containing the first $2m+1$ columns of $U$ and $U_{2}$ containing
the last $n-(2m+1)$ columns, $U=[U_{1}\ U_{2}]$. The matrix $V$ is an
$(2m+1)\times(2m+1)$ orthogonal matrix, i.e. $V^{T}V=VV^{T}=I_{p}$, and
$v_{i}$ and $u_{i}$ denote the $i$th columns of $V$ and $U$, respectively. The
diagonal matrix $\Lambda=diag(s_{1},...,s_{2m+1})$ with strictly positive
decreasing singular values $s_{i}$, i.e. $s_{1}\geq s_{2}\geq...\geq
s_{2m+1}\geq 0$. Then the estimator $q$ has the following form [19]:

$\displaystyle
q=q_{true}+V\Lambda^{-1}U_{1}^{T}\tilde{\xi}=q_{true}+\sum_{i=1}^{2m+1}\frac{1}{s_{i}}v_{i}u_{i}^{T}\tilde{\xi}.$
(5.28)

From $(\ref{SVD})$, it can be seen that the instability of the inverse problem
is caused by the small singular values. In Example 5.3, the singular values of
the sensitivity matrix are calculated at different values of $\lambda$. In
Figure 5.5(b), when $\lambda(\omega)=0.19-10^{-6}i$ ($\omega$ is a plasmon
resonance frequency for Example 5.3), we can see that all the singular values
of $G$ are lager than $\lambda=0.25-10^{-6}i$ and $\lambda=2$. Consequently,
the numerical solution computed by using the plasmon resonance is in an
excellent agreement with the exact solution (see Figure 5.5(a)). Hence, the
plasmon resonance of the metal nanoparticles can correct the singular value of
the sensitivity matrix and overcome the numerical instability. This numerical
result is in agreement with our analysis in Section 3.2 that the sensitivity
of the far-field data to the shape of the underlying domain can be enhanced at
the resonant frequencies, thus reducing the ill-posedness of the inverse
problem.

Figure 5.5: (a) Reconstruction of the shape in Example 5.3 with 1% noise data
and different $\lambda$ as well as the exact solution, under the same
regularization parameter $\mu=0.05$, (b) the distribution of the singular
values associated to different values of $\lambda$.

(a) $\lambda=0.19-10^{-6}i$

(b) $\lambda=2$

Figure 5.6: The numerical results for Example 5.3 with 1% noise in the data
and 95% confidence interval for different $\lambda$.

Next, we study the variations of the confidence intervals with different
$\mathrm{dist}(\lambda,\sigma(\mathcal{K}^{*}_{D}))$. The confidence interval
can quantify the uncertain information of the solution. The numerical results
for Example 5.3 with different $\lambda$ are shown in Figure 5.6, where the
blue region represents the corresponding 95% confidence region. The comparison
with $\lambda=2$ indicates that the confidence region shrinks at
$\lambda=0.19-10^{-6}i$ (the real part of $\lambda=0.19-10^{-6}i$ is the
eigenvalue of $\mathcal{K}^{*}_{D}$ in Example 5.3). It can be clearly
obtained that as $\mathrm{dist}(\lambda,\sigma(\mathcal{K}^{*}_{D}))$ shrinks,
the accuracy of the inversion can be improved and the sensitivity to the
random error can be reduced.

######  Example 5.4.

In this example, we consider that the inclusion is pear, and the polar radius
of the pear is parameterized by,

$q(t)=18/25+3/20\cos(3t),\ \ \ 0\leq t\leq 2\pi.$

In the iteration, we choose a circle of radius 0.73 as the initial guess.

In Example 5.4, we consider reconstructing a more challenging pear-shaped
inclusion and we can actually reach similar conclusions to previous three
Examples. At the same noise level, $\tilde{q}_{N_{e}}$ coincides well with the
exact solution when $\lambda=0.14-10^{-6}i$ (minimum distance of
$\mathrm{dist}(\lambda,\sigma(\mathcal{K}^{*}_{D}))$, and a steady, fast
convergence of the relative error $e_{\gamma}$ in the iteration algorithm is
shown in Figure 5.7. Moreover, the relative error becomes more oscillating as
the noise level $\delta$ increases; see Table 3 for $e_{\gamma}$ associated
with different $\lambda$ and $\delta$.

Figure 5.7: Reconstruction the shape for Example 5.4 with 1% noise data and regularization parameter $\mu=0.5$. (a) Reconstruction with different $\lambda$ and exact solution, (b) the relative error $e_{\gamma}$ versus iterations steps $k$, for various values of $\lambda$.  Table 3:  Numerical results of Example 5.4 for $e_{\gamma}$ with various $\lambda$ and noise level $\delta$ $\lambda$ | $\delta=0.01$ | $\delta=0.05$
---|---|---
$0.14-10^{-6}i$ | $0.0054$ | $0.0625$
$0.25-10^{-6}i$ | $0.279$ | $0.4850$
$2$ | $0.4206$ | $1.009$
|  |

Figure 5.8: Reconstruction the shape for Example 5.4 with 5% noise data with various regularization parameter $\mu$. Table 4: The relative errors for Example 5.4 with different regularization parameters $\mu$ and $\lambda$, under same noise level $\delta=0.05$. | $\lambda=2$ | $\lambda=0.14-10^{-6}i$
---|---|---
$\mu$ | $0.1$ | $0.8$ | $0.09$ | $0.73$ | $0.83$ | $0.206$ | $0.212$ | $0.320$ | $0.416$ | $0.668$ | $0.5$
$e_{\gamma}$ | $0.881$ | $0.720$ | $1.278$ | $0.548$ | $1.096$ | $1.297$ | $1.324$ | $1.034$ | $0.848$ | $0.668$ | $0.0625$

The choice of the regularization parameters is very important in the
algorithmic calculation. In Example 5.4, we fixed the regularization parameter
$\mu=0.5$ associated to different values of $\lambda$. In order to eliminate
the influence of improper regularization parameters on the inversion results,
the regularization parameters can be considered as a random variable with
uniform distribution, i.e. $\mu\sim U(0,1)$. In fact, there are many other
selection strategies of regularization parameters, such as Morozov s
discrepancy principle, The L-curve method and so on (see [50]). But we only
concern the comparative effectiveness of the plasmon resonance case and the
non-resonance case, and then above simple regularization parameters selection
is used. The more advanced selection strategies will be considered in our
future works. In Figure 5.8, when $\lambda=2$, we obtain four groups of random
numbers from uniform distribution as the regularization parameters, while the
corresponding reconstruction results remain disappointing. Table 4 lists 10
randomly generated regularization parameters with the relative error also
exceeding $50\%$.

##  6 Conclusions

We investigate the inverse problem that utilizing the far-field measurement to
reconstruct the shape of of an inclusion. The plasmon resonance is proposed to
enhance the sensitivity of the reconstruction as well as to reduce the ill-
posedness of the inverse problem. In fact, we derive the representation
formula of shape sensitivity functional by using the asymptotic expansion
method. Based on the asymptotic expansion of the eigenvalues and
eigenfunctions of the Neumann-Poincaré operator, we further derive the
delicate spectral representation of the shape sensitivity functional, which
indicates that the sensitivity is improved greatly as the plasmon resonance
occurs. Moreover, we combine the Tikhonov regularization method with the
Laplace approximation to solve the inverse problem. This hybrid method not
only calculates the minimizer accurately and quickly, but also captures the
statistical information of the solution. Finally, extensive numerical
experiments confirm our theoretical analysis and illustrate the promising and
salient features of the proposed reconstruction scheme.

## Acknowledgments

The work of M Ding and G Zheng were supported by the NSF of China (11301168)
and NSF of Hunan (2020JJ4166). The work of H Liu was supported by the Hong
Kong RGC General Research Fund (projects 12302919, 12301218 and 11300821) and
the NSFC/RGC Joint Research Grant (project N_CityU101/21).

## References

  * [1] H. Ammari, Y. T. Chow, H. Liu, Localized sensitivity analysis at high-curvature boundary points of reconstructing inclusions in transmission problems, SIAM J. Math. Anal., in press, 2022.
  * [2] H. Ammari, Y. Chow, H. Liu, Quantum ergodicity and localization of plasmon resonances, arXiv:2003.03696
  * [3] H. Ammari, Y. Chow, H. Liu, M. Sunkula, Quantum integral systems and concentration of plasmon resonance, arXiv:2109.13008
  * [4] H. Ammari, Y. Chow, K. Liu, J. Zou, Optimal shape design by partial spectral data, SIAM J. Sci. Comput., 37 (2015) B855-B883.
  * [5] H. Ammari, G. Ciraolo, H. Kang, H. Lee, G. Milton, Spectral theory of a Neumann-Poincaré-type operator and analysis of cloaking due to anomalous localized resonance, Arch. Ration. Mech. Anal. 208 (2013) 667-692.
  * [6] H. Ammari, Y. Deng, P. Millien, Surface plasmon resonance of nanoparticles and applications in imaging, Arch. Ration. Mech. Anal., 220 (2016) 109-153.
  * [7] H. Ammari, H. Kang, M. Lim, H. Zribi, The generalized polarization tensors for resolved imaging. Part I: Shape reconstruction of a conductivity inclusion, Math. Comp., 81 (2012) 367-386.
  * [8] H. Ammari, P. Millien, M. Ruiz, H. Zhang, Mathematical analysis of plasmonic nanoparticles: the scalar case, Arch. Ration. Mech. Anal., 224(2) (2017) 597-658.
  * [9] K. Ando, H. Kang, H. Liu, Plasmon resonance with finite frequencies: A validation of the quasi-static approximation for diametrically small inclusions, SIAM J. Appl. Math., 76 (2016) 731-749.
  * [10] J. Anker, W. Hall, O. Lyandres, N. Shah, J. Zhao, R. Van Duyne, Biosensing with plasmonic nanosensors, and Applications, Cambridge University Press, New York, 2010.
  * [11] G. Baffou, C. Girard, R. Quidant, Mapping heat origin in plasmonic structures, Phys. Rev. Lett, 104 (2010) 136805.
  * [12] E. Blåsten, H. Li, H. Liu, Y. Wang, Localization and geometrization in plasmon resonances and geometric structures of Neumann-Poincaré eigenfunctions, ESAIM: Math. Model. Numer. Anal., 54 (2020) 957-976.
  * [13] E. Blåsten, H. Liu, On corners scattering stably and stable shape determination by a single far-field pattern, Indiana Univ. Math. J., 70 (2021), no. 3, 907–947.
  * [14] E. Blåsten, H. Liu, Scattering by curvatures, radiationless sources, transmission eigenfunctions, and inverse scattering problems, SIAM J. Math. Anal. 53 (2021), no. 4, 3801–3837.
  * [15] Blåsten, H. Liu, Recovering piecewise constant refractive indices by a single far-field pattern, Inverse Problems, 36 (2020), no. 8, 085005, 16 pp.
  * [16] G. Bouchitté, B. Schweizer, Cloaking of small objects by anomalous localized resonance, Quart. J. Mech. Appl. Math., 63 (2010) 437-463.
  * [17] X. Cao, H. Diao, H. Liu, J. Zou, On nodal and generalized singular structures of Laplacian eigenfunctions and applications to inverse scattering problems, J. Math. Pures Appl. (9), 143 (2020), 116–161.
  * [18] Y. T. Chow, Y. Deng, Y. He, H. Liu, X. Wang, Surface-localized transmission eigenstates, super-resolution imaging, and pseudo surface plasmon modes, SIAM J. Imaging Sci., 14 (2021), no. 3, 946–975.
  * [19] A. Cintrón-Arias, H. Banks, A. Capaldi, A. Lloyd, A sensitivity matrix methodology for inverse problem formulation, J. Inverse Ill-Pose. P., 17 (2009) 1-20.
  * [20] Y. Deng, H. Liu, G. Zheng, Mathematical analysis of plasmon resonances for curved nanorods, J. Math. Pures Appl. (9), 153 (2021), 248–280.
  * [21] H. Diao, X. Cao, H. Liu, On the geometric structures of transmission eigenfunctions with a conductive boundary condition and applications, Comm. Partial Differential Equations, 46 (2021), no. 4, 630–679.
  * [22] Y. Deng, H. Li, H. Liu, Analysis of surface polariton resonance for nanoparticles in elastic system, SIAM J. Math. Anal., 52 (2020) 1786-1805.
  * [23] A. Doicu, T. Trautmann, F. Schreier, Numerical regularization for atmospheric inverse problems, Springer Science $\&$ Business Media, 2010.
  * [24] X. Fang, Y. Deng, J. Li, Plasmon resonance and heat generation in nanostructures, Math. Method. Appl. Sci., 38(18) (2015) 4663-4672.
  * [25] T. Feng, H. Kang, H. Lee, Construction of GPT-vanishing structures using shape derivative, J. Comput. Math., 35(005) (2017) 569-585.
  * [26] Y. Gao, H. Liu, X. Wang, K. Zhang, On an artificial neural network for inverse scattering problems, J. Comput. Phys. 448 (2022), 110771.
  * [27] D. Grieser, The plasmonic eigenvalue problem, Rev. Math. Phys. 26 (2014) 1450005.
  * [28] M. Hanke, Recent progress in electrical impedance tomography, Inverse Probl., 19 (2003), S65–S90.
  * [29] M. Hanke, A regularizing Levenberg-Marquardt scheme with applications to inverse groundwater filtration problems, Inverse Probl., 13 (1997) 79.
  * [30] M. Hintermüller, A. Laurain, I. Yousept, Shape sensitivities for an inverse problem in magnetic induction tomography based on the eddy current model, Inverse Probl., 31 (2015) 065006.
  * [31] M. Iglesias, K. Law, A. Stuart, Evaluation of Gaussian approximations for data assimilation in reservoir models, Comput. Geosci., 17 (2013) 851-885.
  * [32] P. Jain, K. Lee, I. El-Sayed, M. El-Sayed, Calculated absorption and scattering properties of gold nanoparticles of different size, shape, and composition: applications in biomedical imaging and biomedicine, J. Phys. Chem. B, 110 (2006) 7238-7248.
  * [33] R. Kress, Linear Integral Equations, 2nd edition, Springer, 1999.
  * [34] H. Li, J. Li, H. Liu, On quasi-static cloaking due to anomalous localized resonance in $\mathbb{R}^{3}$, SIAM J. Appl. Math., 75 (2015), no. 3, 1245–1260.
  * [35] H. Li, H. Liu, On anomalous localized resonance and plasmonic cloaking beyond the quasistatic limit, Proc. Roy. Soc. A, 474 (2018).
  * [36] H. Li, J. Li, H. Liu, On novel elastic structures inducing polariton resonances with finite frequencies and cloaking due to anomalous localized resonances, J. Math. Pures Appl. (9), 120 (2018) 195-219.
  * [37] H. Li, S. Li, H. Liu, X. Wang, Analysis of electromagnetic scattering from plasmonic inclusions beyond the quasi-static approximation and applications, ESAIM: Math. Model. Numer. Anal., 53 (2019) 1351-1371.
  * [38] S. Link, M. El-Sayed, Shape and size dependence of radiative, non-radiative and photothermal properties of gold nanocrystals, Int. Rev. Phys. Chem., 19 (2000) 409-453.
  * [39] H. Liu, C. H. Tsou, Stable determination of polygonal inclusions in Calderón’s problem by a single partial boundary measurement, Inverse Problems, 36 (2020), no. 8, 085010, 23 pp.
  * [40] H. Liu, C. H. Tsou and W. Yang, On Calderón’s inverse inclusion problem with smooth shapes by a single partial boundary measurement, Inverse Problems, 37 (2021), no. 5, Paper No. 055005, 18 pp.
  * [41] I. Mayergoyz, D. Fredkin, Z. Zhang, Electrostatic (plasmon) resonances in nanoparticles, Phys. Rev. B, 72 (2005) 155412.
  * [42] G. Milton, N. Nicorovici, On the cloaking effects associated with anomalous localized resonance, Proc. R. Soc. A, 462 (2006) 3027-3059.
  * [43] J. Nam, C. Thaxton, C. Mirkin, Nanoparticle-based bio-bar codes for the ultrasensitive detection of proteins, Science 301 (2003) 1884-886.
  * [44] D. Nicholls, X. Tong, A high-order perturbation of surfaces algorithm for the simulation of localized surface plasmon resonances in two dimensions, J. Sci. Comput. 76 (2018) 1370-1395.
  * [45] M. Ordal, L. Long, R. Bell, S. Bell, R. Bell, R. Alexander, C. Ward, Optical properties of the metals al, co, cu, au, fe, pb, ni, pd, pt, ag, ti, and w in the infrared and far infrared, Appl. Opt., 22 (1983) 1099-1119.
  * [46] G. Raschke et al. , Biomolecular recognition based on single gold nanoparticle light scattering , Nano Lett., 3(7) (2003) 935-938.
  * [47] D. Sarid, W. Challener, Modern Introduction to Surface Plasmons: Theory, Mathematical Modeling, Nat. Mater., 7 (2008) 442-453.
  * [48] C. Schillings, B. Sprungk, P. Wacker, On the Convergence of the Laplace Approximation and Noise-Level-Robustness of Laplace-based Monte Carlo Methods for Bayesian Inverse Problems, Numerical Mathematik, 145 (2020) 915-971.
  * [49] S. Schultz, D. Smith, J. Mock, D. Schultz, Single-target molecule detection with nonbleaching multicolor optical immunolabels, Proc. Natl Acad. Sci. USA, 97 (2000) 996-1001.
  * [50] C. Vogel, Computational Methods for Inverse Problems, SIAM, 2002.
  * [51] W. Yin, W. Yang and H. Liu, A neural network scheme for recovering scattering obstacles with limited phaseless far-field data, J. Comput. Phys., 417 (2020), 109594, 18 pp.
  * [52] G. Zheng, Mathematical analysis of plasmonic resonance for 2-D photonic crystal, J. Differential Equations, 266 (2019), no. 8, 5095–5117.