Shape reconstructions by using plasmon resonances

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.


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 [10,11,36,47,50,53] 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,12,22,23,31,42,45]. 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 [2,3,12], 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, 16, 22, 25, 38-41, 46, 56]. 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 ⊂ R 2 be a bounded domain with a connected complement R 2 \D. Given a harmonic function H, we consider the following electrostatic problem: where ε = ε D χ(D) + ε m χ(R 2 \D), (1.2) and χ 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 ε signifies the permittivity parameter of the medium. Throughout, we shall assume that the background parameter ε m is a positive constant, whereas ε 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 ∂D, from the measurement data u s = u − H on ∂Ω with D Ω associated with a fixed incident field H. For simplicity, we take Ω to be a central ball of radius r ∈ 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, 13-15, 17, 18, 20, 21, 24, 27, 30, 32, 34, 43, 44, 55] 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. 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. 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. 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.

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 ∂D are denoted by ·, · and · in short, respectively. The single layer potential S D and double layer potential D D associated with D are given by where ϕ ∈ L 2 (∂D) is the density function, and the Green function Γ(x, y) to the Laplacian in R 2 is given by The notations u| ± and ∂u ∂ν | ± denote the traces on ∂D from the outside and inside of D, respectively. The following jump relations hold [4,5]: where K * D is known as the Neumann-Poincaré (NP) operator defined by Next, we recall some useful facts about the NP operator K * D [5,22,45].  (iv) The following representation formula holds: for any ϕ ∈ H * (∂D),

Plasmon resonance
We next briefly discuss the mathematical framework of plasmon resonance. We first give the form of the solution of equation (1.1). From [1,5,22], we have where φ ∈ L 2 0 (∂D) := {φ ∈ L 2 (∂D); ∂D φ = 0} satisfies with λ given by . (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 ε D of the plasmon material is changing with respect to the operating frequency ω. The ε D can be described by the Drude's model (see [6,28,51]), where ε 0 is the electric permittivity of the vacuum and is given by ε 0 = 9 · 10 −12 F/m. The plasmon frequency ω p and the damping parameter γ are strictly positive. When ω < ω p the real part of ε D can be negative. In particular, for the gold, the value of the plasmon frequency and the damping parameter are ω p = 2 · 10 15 s −1 and γ = 10 14 s, respectively. Let the operating frequency ω = 6 · 10 14 Hz (visible light frequency), and from Drude's model (2.6), the electric permittivity of gold nanoparticle is calculated as ε D ≈ (−9.8108 + 1.8018i)ε 0 . Now by applying the spectral decomposition of the K * D to the integral equation (2.4), the density φ ∈ L 2 0 (∂D) becomes where λ j are eigenvalues of K * D and they satisfy |λ j | < 1 2 . When the real part of ε D (ω) is negative, it holds that |Re(λ(ω)) < 1 2 |, where and also in what follows Re(λ(ω)) signifies the real part of λ(ω). The frequency ω is called a plasmon resonance frequency if it satisfies for some j, where λ j is an eigenvalue of the NP operator K * D . For this reason, (2.8) is called the plasmon resonance condition. It is emphasized that only when the imaginary part Im(λ(ω)) = 0 (lossless metal), the resonance condition (2.8) can imply the density φ in (2.7) blow up.
In (2.7) the density φ will be amplified when the plasmon resonance condition is reached provided that ∂H ∂ν , ϕ j H * (∂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 [22]). 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 ε m = ε 0 (vacuum) and notice that, as D is a disk, the spectrum of K * D are {0, 1 2 }. If λ j = 0, from the plasmon resonance condition (2.8) and the fact that γ is sufficiently small in (2.6), we nearly get Re(ε D (ω)) = −ε m . This relationship is called the Fröhlich condition (see [48]). Furthermore, by Drude's model (2.6), we can obtain the resonance frequency ω = ω 2 p 2 − γ 2 . When γ = 0 (lossless metal), the resonance frequency ω = ωp √ 2 or f = ωp 2π √ 2 , which is called the Fröhlich frequency (see [51]).

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 F : X → H 1 2 (∂Ω) on a subset X ⊂ C 2 (∂D): where u s | ∂Ω = (u − H)| ∂Ω is the far-field measurement associated with (1.1) on the boundary ∂Ω. For a small ∈ R + , we let ∂D be an -perturbation of D, i.e., where h ∈ C 1 (∂D), and ν is the outward unit normal vector to ∂D. The solution u to (1.1) with D has the following representation formula where the density functionφ is the solution to Let Ψ be the diffeomorphism from ∂D to ∂D given by Moreover, we denoteν the outward unit normal vector to ∂D and dσ the line element of ∂D .
The following expansions ofν and dσ hold [7]: Here and throughout the rest of the paper, τ (x) signifies the curvature of ∂D at x, T is the unit tangential vector to ∂D, and h is the tangential derivative of h on ∂D, i.e., h = ∂h ∂T , where ∂ ∂ν denotes the outward normal derivative on ∂D . In view of (3.9) and (3.10), in order to obtain the asymptotic formula of the perturbed far field u s | ∂Ω := F(∂D ), we need to obtain the corresponding asymptotic expansions of the operators K * D andφ . The following lemmas can be found in [7].
Here, p.v. stands for the Cauchy principal value.
In fact, we can rewrite the operator K D in terms of more familiar operators as follows, (3.14) where C is a constant depending only on the C 2 -norm of ∂D and h C 1 and Furthermore, for ∂H ∂ν of (3.10), we can further obtain by using (3.11), that The asymptotic expansion of F(∂D ) under small perturbations of the boundary ∂D in terms of the asymptotic parameter is given in the following theorem. We would like to point out that its proof is adapted from those in [7,29].
the following asymptotic expansion holds: Proof. From the form of the solution (3.9), it holds that whereφ is solution to (3.10), dσ(ỹ) has an expansion as (3.12), and Γ(x,ỹ) has the Taylor expansion as follows By virtue of Lemma 3.2, we can obtain that We treat l 1 , l 2 separately. Since and noting that H(y) is a harmonic function and (3.21), l 1 has the following form Next, by (3.14), we have from which, together with the use of Proposition 4.1 in [7], we can further show that Thus we obtain that From (3.18) and (3.19), we know that Thus it follows that 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 | ∂Ω with respect to the shape of ∂D is defined as 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 F(∂D) (cf. [7]). Furthermore, by using Theorem 3.3, the shape sensitivity function can be rewritten as where P is defined in (3.20).

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 .
Proof. Since K * D [ϕ j, ] = λ j, ϕ j, , and by Lemma 3.1, we can obtain the following result which implies that By inner producting ϕ l in both sides of (3.22), we further have Since the operator K * D is self-adjoint, we have and thus λ (1) j, ϕ j . Then it is straightforward to see that The proof is complete.
From (3.16), (3.24) and Lemma 3.4, we deduce that Therefore it follows that The proof is complete. (λ−λ j ) 2 . Thus, for the sufficiently small loss (Im(ε D ) → 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 = o dist(λ, σ(K * D )) 2 , (as dist(λ, σ(K * D )) → 0), the shape sensitivity function can be represented as where T (x) is defined by (3.23). Remark 3.3. Similar to Remark 3.2, one sees that when λ = ε D +εm 2(ε D −εm) is very close to an eigenvalue λ j of the NP operator, i.e., the plasmon resonance occurs, the shape sensitivity function SSF (∂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).

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 ∂D. Here we assume that ∂D is starlike boundary curve with respect to the origin, i.e. there exists q ∈ C 2 [0, 2π] such that The admissible set X = {q ∈ C 2 ([0, 2π]) : q > 0}, and the forward operator maps X into H 1 2 (∂Ω). Without change of notation we write F(q) = u s (x)| ∂Ω .
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, where, u s,δ = (u s,δ 1 , u s,δ 2 , · · · , u s,δ n ) signfies the discrete measurement data on ∂Ω. Here, · 2 denotes the 2-norm in R n , and µ is the regularization parameter. Moreover, the measurement data u s,δ and the exact data u s satisfies u s,δ − u s 2 ≤ δ. There are many iterative algorithms to solve the variational problem (4.25). In this paper, we apply the Levenberg-Marquardt method [26,33] 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 F in (4.25) can be replaced approximatively by its linearization around q * , i.e.
The nonlinear inverse problem F(q) = u s,δ can then be converted to a linear inverse problem Thus, minimizing (4.25) can easily be seen to minimizing where, δq = q − q * . We formulate the Levenberg-Marquardt iteration as Algorithm 1.

Algorithm 1 Levenberg-Marquardt algorithm
Input: The regularization parameter µ, the noise level δ 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,δ ; 3: While k < M 4: Compute the forward problem and set F k = u s,δ − F(q k ); 5: Compute the Jacobian matrix G = F (q k ); 6: Compute δq k = (G T G + µI) −1 (G T F k ); 7: Update the solution: q k+1 = q k + δq k ; 8: k ← k + 1; 9: Terminate if δq k < 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 ξ is assumed to be an independent and identically distributed Gauss random vector with mean zero and the covariance matrix B = δ 2 I (here I is the unit matrix), we can write the minimization functional as where · B is a covariance weighted norm given by · B = δ −1 I· 2 , and the minimizer of J B defines the maximum a posteriori estimator The Laplace approximation replaces the complicated posterior with a normal distribution located at the maximum a posteriori value q MAP . Its essence is a linearization around the MAP point q MAP (cf. [52]). It consists of approximating the posterior measure (or distribution) by ω ≈ N (q MAP , C MAP ), where 27) and G is the Jacobian matrix of the forward operator F at the point q. Notice that the covariance formula (4.27) only uses the first order derivatives of F. The implementation of the Laplace approximation is presented in the following algorithm [35]. In this algorithm, samples generated by algorithm 2 are drawn from N (q MAP , C MAP ), and so the ensemble of N e realizations {q j } Ne j=1 provides an approximation to N (q MAP , C MAP ) and hence the posterior. Finally, we use the mean of the sampleq = 1 Ne Ne j=1 q j as an approximation of q MAP , where the convergence ofq follows from the strong law of large numbers. From the classical Gaussian statistic theory, we find thatq is consistent and the best unbiased estimate of q MAP .

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 ∂Ω is given by the circle of radius 3 and centred at the origin, that is ∂Ω = {3(cos t, sin t), 0 ≤ t ≤ 2π}. Set the incident field H(x) = x 1 . We solve the direct problem through the Nyström method [37], 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 ∂D by the trigonometric series where m ∈ N and the vector q = (a 0 , ..., a m , b 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 The noisy measured data are generated by u s,δ = u s (x) + δξ, x ∈ ∂Ω, where u s (x) is the exact data, δ indicates the noise level, and ξ is the Gaussian random vector with a zero mean and unit standard deviation.
The accuracy of the approximate solutionq Ne by the LA algorithm is characterised by comparing to the exact solution q(x) via the relative error In particular, we specifically take the permittivity ε m = 24.8F/m (such as alcohol) in the first example, and in all of the examples we set ε m = ε 0 . Moreover, except for Example 5.1, the spectrum of K * D needs to calculated numerically, from which one can then compute the corresponding plasmon resonance frequency. We use the Nyströ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 ≤ t ≤ 2π.
In the iteration algorithm, we choose a circle of radius 0.6 as the initial guess.
First, we investigate the influence of different λ on the reconstruction effect. The numerical results for Example 5.1 with λ taking different values are shown in Figure 5.1 (a). It is well known that the eigenvalues of K * D are {0, 1 2 } for D being a disk. If setting λ 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 ω 1 = 3.63 · 10 14 Hz, and according to the experimental data of gold nanorod in [49] (the data coincide with the Drude model), we can get Re(ε D (ω 1 )) = −24.8, and the imaginary part Im(ε D (ω 1 )) = 0.797. Furthermore, by (2.5), it implies that λ ≈ 0 − 8 · 10 −3 i.
Next, we take ω 2 = 2.30 · 10 14 Hz and based on the experimental data in [49], we can calculate λ = 0.244 − 6.3 · 10 −3 i. It is obvious that when taking ω 2 , the real part of λ(ω 2 ) satisfies |Re(λ(ω 2 )) < 1 2 |, which is only an approximate resonance frequency. As seen in Figure 5   It can be seen that the relative error can quickly reach the convergence and remain small when ω 1 is the plasmon resonance frequency. When λ = 0.244 − 6.3 · 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 λ = −0.98. The numerical results for Example 5.1, which has 5% noise in the data, are shown in Figure  5.2(a) for different λ. 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 ω is a plasmon resonance frequency, the reconstruction is very good even at larger error levels. For λ = −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.
We choose a circle of radius 0.72 as the initial guess.
The numerical results for Example 5.2 with λ taking different values are shown in Figure 5.3 (a) and the relative error e γ versus the iteration number with different λ in Figure 5.3 (b). The    Table 2 for e γ with different λ and δ.
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) = cos 2 t + 0.26 sin 2 (t + 0.5), 0 ≤ t ≤ 2π.
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 where U is an n × n orthogonal matrix, i.e. U T U = U U 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) × (2m + 1) orthogonal matrix, i.e. V T V = V V T = I p , and v i and u i denote the ith columns of V and U , respectively. The diagonal matrix Λ = diag(s 1 , ..., s 2m+1 ) with strictly positive decreasing singular values s i , i.e. s 1 ≥ s 2 ≥ ... ≥ s 2m+1 ≥ 0. Then the estimator q has the following form [19]: 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. Next, we study the variations of the confidence intervals with different dist(λ, σ(K * D )). The confidence interval can quantify the uncertain information of the solution. The numerical results for Example 5.3 with different λ are shown in Figure 5.6, where the blue region represents the corresponding 95% confidence region. The comparison with λ = 2 indicates that the confidence region shrinks at λ = 0.19 − 10 −6 i (the real part of λ = 0.19 − 10 −6 i is the eigenvalue of K * D in Example 5.3). It can be clearly obtained that as dist(λ, σ(K * D )) shrinks, the accuracy of the inversion can be improved and the sensitivity to the random error can be reduced.
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,q Ne coincides well with the exact solution when λ = 0.14 − 10 −6 i (minimum distance of dist(λ, σ(K * D )), and a steady, fast convergence of the relative error e γ in the iteration algorithm is shown in Figure  5.7. Moreover, the relative error becomes more oscillating as the noise level δ increases; see Table  3 for e γ associated with different λ and δ.
The choice of the regularization parameters is very important in the algorithmic calculation. In Example 5.4, we fixed the regularization parameter µ = 0.5 associated to different values of λ. 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. µ ∼ U (0, 1). In fact, there are many other selection strategies of regularization parameters, such as Morozovs discrepancy principle, The L-curve method and so on (see [54]). 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 λ = 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%.

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.