The Imaging of Small Perturbations in an Anisotropic Media

TL;DR

This paper introduces an asymptotic analysis-based method for imaging small defects in anisotropic media using far field data, combining the MUSIC algorithm and transmission eigenvalue analysis.

Contribution

It develops a novel approach integrating asymptotic expansion, MUSIC algorithm, and eigenvalue analysis for detecting and characterizing small perturbations in anisotropic materials.

Findings

Successful reconstruction of defect locations using MUSIC algorithm.

Analysis of how defects influence transmission eigenvalues.

Numerical validation demonstrating method effectiveness.

Abstract

In this paper, we employ asymptotic analysis to determine information about small volume defects in a known anisotropic scattering medium from far field scattering data. The location of the defects is reconstructed via the MUSIC algorithm from the range of the multi-static response matrix derived from the asymptotic expansion of the far field pattern in the presence of small defects. Since the same data determines the transmission eigenvalues corre- sponding to the perturbed media, we investigate how the presence of the defects changes the transmission eigenvalues and use this information to recover the strength of the small defects. We provide convergence results on transmission eigenvalues as the size of the defects tends to zero as well as derive the first correction term in the asymptotic expansion of the simple transmission eigen- values. Numerical examples are presented to show the viability of our imaging method.

Figures

Tables

Table 1. Table 1: The estimated order of convergence for two eigenvalues where D 0 subscript 𝐷 0 D_{0} is two disks of radius ε 𝜀 \varepsilon centered at ( 0.25 , 0 ) 0.25 0 (0.25,0) and ( − 0.25 , − 0.25 ) 0.25 0.25 (-0.25,-0.25) .

$\epsilon$	${\mathrm{EOC}}_1$	${\mathrm{EOC}}_2$
1/4	$-$	$-$
1/8	2.4423	0.8252
1/16	2.3365	2.0673
1/32	2.0881	2.0001
1/64	2.1705	2.1549

Table 2. Table 2: The estimated order of convergence for two eigenvalues where D 0 subscript 𝐷 0 D_{0} is the disk of radius ε 𝜀 \varepsilon centered at the origin.

$\epsilon$	${\mathrm{EOC}}_1$	${\mathrm{EOC}}_2$
1/4	$-$	$-$
1/8	1.9304	2.2957
1/16	2.1519	2.2278
1/32	2.1161	2.0304
1/64	3.1701	2.5851

Equations

$${A}_{b}(x)=\left\{\begin{array}[]{rl}I&\quad\,x\in\mathbb{R}^{d}\setminus\overline{D}\\ A(x)&\quad\,x\in D\end{array}\right.\quad\text{ and }\quad{n}_{b}(x)=\left\{\begin{array}[]{rl}1&\quad\,x\in\mathbb{R}^{d}\setminus\overline{D}\\ n(x)&\quad\,x\in D.\end{array}\right.$$

$$\nabla \cdot A_b(x)\nabla u_b+k^2{n}_b(x)u_b=0$$

$$\underset{r\to \infty }{\lim }{r}^{\frac{d-1}{2}}\left(\frac{\partial u_b^s}{\partial r}-ik{u}_b^s\right)=0,$$

$$u_b^s(x,\hat{y})=\frac{{\text{e}}^{ik|x|}}{|x{|}^{\frac{d-1}{2}}}\left\{u_b^{\infty }(\hat{x},\hat{y})+\mathcal{O}\left(\frac{1}{|x|}\right)\right\}\text{ as }|x|\to \infty$$

$$|z_i-z_j|\ge c_0>0\text{for all }i\ne j\text{with }i,j=1,2,\dots M\text{mbox}and$$

$$\text{dist}(z_m,\partial D)\ge c_0>0\text{for all }m=1,2,\dots M.$$

$${A}_{\varepsilon}(x)=\left\{\begin{array}[]{rl}A_{m}&\quad\,x\in(z_{m}+\varepsilon B_{m})\\ {A}(x)&\quad\,x\in\mathbb{R}^{d}\setminus\overline{D}_{\varepsilon}\end{array}\right.\quad\mbox{and}\quad{n}_{\varepsilon}(x)=\left\{\begin{array}[]{rl}n_{m}&\quad\,x\in(z_{m}+\varepsilon B_{m})\\ {n}(x)&\quad\,x\in\mathbb{R}^{d}\setminus\overline{D}_{\varepsilon}.\end{array}\right.$$

$$\nabla \cdot A_{\epsilon }(x)\nabla u_{\epsilon }+k^2{n}_{\epsilon }(x)u_{\epsilon }=0$$

$$\underset{r\to \infty }{\lim }{r}^{\frac{d-1}{2}}\left(\frac{\partial u_{\epsilon }^s}{\partial r}-ik{u}_{\epsilon }^s\right)=0.$$

$$\underset{x\in D}{\inf }\underset{|\xi |=1}{\inf }\overline{\xi }\cdot A(x)\xi =A_{min}>0$$

$$\underset{m=1\dots M}{\min }\underset{|\xi |=1}{\inf }\overline{\xi }\cdot A_m\xi =a_{min}>0$$

$$(Fg)(\hat{x})={\int }_{\mathbb{S}}\left[u_{\epsilon }^{\infty }(\hat{x},\hat{y})-u_b^{\infty }(\hat{x},\hat{y})\right]g(\hat{y})d\hat{y}\hat{x}\in \mathbb{S}$$

$$\nabla \cdot A_{\epsilon }\nabla w+k_{\epsilon }^2{n}_{\epsilon }w=0\text{and }\Delta v+k_{\epsilon }^{2}v=0$$

$$w=v\text{and }\frac{\partial w}{\partial {\nu }_{A_{\epsilon }}}=\frac{\partial v}{\partial \nu }$$

$$\nabla \cdot A(x)\nabla \mathbb{G}(\cdot ,z)+k^{2}n(x)\mathbb{G}(\cdot ,z)=-\delta (\cdot -z)\text{ in }{\mathbb{R}}^d$$

$$\underset{r\to \infty }{\lim }{r}^{\frac{d-1}{2}}\left(\frac{\partial \mathbb{G}(\cdot ,z)}{\partial r}-ik\mathbb{G}(\cdot ,z)\right)=0.$$

$${\mathbb{G}}^{\infty }(\hat{x},z)=\gamma u_b(z,-\hat{x})\text{ where }\gamma =\frac{{\text{e}}^{i\pi /4}}{\sqrt{8\pi k}}\text{ in }{\mathbb{R}}^2\text{ and }\gamma =\frac{1}{4\pi }\text{ in }{\mathbb{R}}^3.$$

$$u_{\epsilon }(x,\hat{y})=u_b(x,\hat{y})+\sum _{m=1}^M{k}^2\underset{z_m+\epsilon B_m}{\int }(n_m-n)\mathbb{G}(x,z)u_{\epsilon }(z,\hat{y})dz$$

$$+\underset{z_m+\epsilon B_m}{\int }(A-A_m){\nabla }_z\mathbb{G}(x,z)\cdot \nabla u_{\epsilon }(z,\hat{y})dz.$$

$$\displaystyle u^{\infty}(\hat{x},\hat{y})=\gamma\varepsilon^{d}k^{2}\sum\limits_{m=1}^{M}|B_{m}|\big{(}n_{m}-n\big{)}u_{b}(z_{m},-{\hat{x}})u_{b}(z_{m},\hat{y})$$

$$+\gamma {\epsilon }^d\sum _{m=1}^M{\mathbf{M}}^{(m)}\nabla u_b(z_m,-\hat{x})\cdot \nabla u_b(z_m,\hat{y})+o({\epsilon }^d),$$

$${\mathbf{M}}_{i,j}^{(m)}=e_i\cdot (A_m-A)e_j+\underset{\partial B_m}{\int }\left[\nu (y)\cdot (A_m-A)e_j\right]{\varphi }_i^{+}(y)d{s}_y$$

$$\nabla \cdot A(x)\nabla {\varphi }_i=0\text{ in }{\mathbb{R}}^d\setminus {\overline{B}}_m$$

$$\nabla \cdot A_m\nabla {\varphi }_i=0\text{ in }{B}_m$$

$${\varphi }_i^{-}-{\varphi }_i^{+}=x_i\text{ on }\partial B_m$$

$$\frac{\partial }{\partial {\nu }_{A_m}}{\varphi }_i^{-}-\frac{\partial }{\partial {\nu }_A}{\varphi }_i^{+}=\frac{\partial }{\partial {\nu }_A}{x}_i\text{ on }\partial B_m$$

$${\varphi }_i(x)=\mathcal{O}\left(\frac{1}{|x{|}^{d-1}}\right)$$

$$\displaystyle{\bf F}_{i,j}=\gamma\varepsilon^{d}\sum\limits_{m=1}^{M}k^{2}|B_{m}|\big{(}n_{m}-n\big{)}u_{b}(z_{m},-{\hat{x}}_{i})u_{b}(z_{m},{\hat{y}}_{j})+$$

$$+\gamma {\epsilon }^d\sum _{m=1}^M{\mathbf{M}}^{(m)}\nabla u_b(z_m,-{\hat{x}}_i)\cdot \nabla u_b(z_m,{\hat{y}}_j).$$

$${\bf g}_{z}=\big{(}\,u_{b}(z,-{\hat{x}}_{1}),\,\dots,u_{b}(z,-{\hat{x}}_{N})\,\big{)}^{\top}$$

$${\bf g}_{z,b}=\big{(}b\cdot\nabla u_{b}(z,-\hat{x}_{1}),\,\dots\,,b\cdot\nabla u_{b}(z,-\hat{x}_{N})\big{)}^{\top},$$

$${\mathbf{g}}_{z,(1,b)}={\mathbf{g}}_z+{\mathbf{g}}_{z,b}.$$

$$\mathcal{I}(z)=\left[\sum\limits_{j=r+1}^{N}\left|\big{(}{\bf g}_{z,(1,b)},{\bf w}_{j}\big{)}\right|_{\ell^{2}}^{2}\right]^{-1}<\infty\quad\text{ if and only if }\quad z\in\{z_{m}\,:\,m=1,\dots,M\}.$$

$${\bf F}=\big{[}u^{\infty}_{\varepsilon}(\hat{x}_{i},\hat{x}_{j})-u^{\infty}_{b}(\hat{x}_{i},\hat{x}_{j})\big{]}_{i,j=1}^{N},$$

$$\hat{x}_{i}=\Big{(}\cos\left({2\pi(i-1)}/{N}\right),\sin\left({2\pi(i-1)}/{N}\right)\Big{)}\quad\text{ for }\quad i=1,\dots,N.$$

$$\underset{D}{\int }{A}_{\epsilon }\nabla w\cdot \nabla {\overline{\phi }}_1-\nabla v\cdot \nabla {\overline{\phi }}_2-k_{\epsilon }^2(n_{\epsilon }w{\overline{\phi }}_1-v{\overline{\phi }}_2)dx=0\text{for all}({\phi }_1,{\phi }_2)\in X(D).$$

Full text

Keywords: Inhomogeneous media, anisotropic media, inverse scattering, MUSIC, transmission eigenvalues, asymptotic methods.

1 Introduction

The imaging of anisotropic media from scattering data is a challenging problem mainly due to the non-uniqueness issue [13]. Yet, in many applications in medical imaging and non-destructive testing, the scattering media exhibit anisotropic properties in the interaction with probing waves. The so-called qualitative methods in inverse scattering [5] provide imaging techniques to obtain information on changes in material properties of a known anisotropic media. This work concerns the imaging of small volume (possibly anisotropic) perturbations of a known anisotropic inhomogeneous media in acoustic wave propagation (for the case of ${\mathbb{R}}^{3}$) or specially polarized electromagnetic wave propagation (for the case of ${\mathbb{R}}^{2}$). Combining asymptotic analysis with MUSIC and the related transmission eigenvalue problem we derive a range test for the location of small perturbations and computable formulas that provide information about the strength (involving the contrast and geometrical features) of the small perturbation. There is a vast literature on the MUSIC algorithm for a variety of scattering problems [1], [2], [21] and we recall here its formulation for the anisotropic inhomogeneous media. The asymptotic analysis of the transmission eigenvalue problem for isotropic media is studied in [10] and [11]. One of the main contributions of this study is the asymptotic analysis of the transmission eigenvalue problem for anisotropic media with the first order correction term for the perturbation of the eigenvalues. Note that the transmission eigenvalue problem is non-linear and non-selfadjoint, and the mathematical structure of this problem for anisotropic media is different from the isotropic case. In addition, we show how to use the asymptotic expansion for the perturbation of transmission eigenvalues together with the MUSIC algorithm to image small volume perturbations of anisotropic media.

Let us now precisely formulate the problem under consideration. To this end let $D\subset\mathbb{R}^{d}$ (for $d=2$ or 3) be a bounded domain with piecewise smooth boundary which denotes the support of the anisotropic media to be tested. The real valued symmetric matrix $A(x)\in{C}^{1}(D,\mathbb{R}^{d\times d})$ with smooth entries and the smooth function $n\in{C}^{1}(D)$ represent the constitutive parameters for the unperturbed (“healthy”) anisotropic media. Without loss of generality we assume that outside the scatterer $D$ the background media has refractive index scaled to one, i.e. $A(x)=I$ and $n(x)=1$ in $x\in\mathbb{R}^{d}\setminus\overline{D}$, where $I$ denotes the identity matrix. We define

[TABLE]

Now the scattering of a time harmonic incident plane wave $\text{e}^{\text{i}kx\cdot\hat{y}}$ with incident direction $\hat{y}\in\mathbb{S}$ by the unperturbed media (i.e. without defects) is mathematically formulated as: find $u_{b}\in H^{1}_{loc}(\mathbb{R}^{d})$ with $u_{b}=u^{s}_{b}+\text{e}^{\text{i}kx\cdot\hat{y}}$ such that

[TABLE]

where $\mathbb{S}$ denotes the unit circle/sphere, $r=|x|$, and the Sommerfeld radiation condition (2) is satisfied uniformly with respect to $\hat{x}=x/|x|$. Here $u_{b}$ is the total field in the background (including the homogeneous part and the media of compact support $\overline{D}$) and $u^{s}_{b}$ is the scattered field due to the region $D$. Recall that the scattered radiating field $u_{b}^{s}(\cdot,\hat{y})$, which depends on the incident direction $\hat{y}$, has the following asymptotic expansion [12]

[TABLE]

where $\hat{x}:=x/|x|$, and $u_{b}^{\infty}(\hat{x},\hat{y})$, which depends on the incident direction $\hat{y}$ and observation direction $\hat{x}$, is the corresponding far field pattern. Now we consider the small defective regions that are given by $z_{m}+\varepsilon B_{m}$ where $B_{m}$ is a smooth deformation of a ball centered at the origin. Let $A_{m}$ and $n_{m}$ be constant constitutive parameters for the defective regions given by $z_{m}+\varepsilon B_{m}$ and assume that

[TABLE]

[TABLE]

The union of the defective regions is denoted by $D_{\varepsilon}=\bigcup\limits_{m=1}^{M}(z_{m}+\varepsilon B_{m})$ and we let

[TABLE]

The scattering problem for the media with the defective region $D_{\varepsilon}$ now reads: find $u_{\varepsilon}\in H^{1}_{loc}(\mathbb{R}^{d})$ with $u_{\varepsilon}=u_{\varepsilon}^{s}+\text{e}^{\text{i}kx\cdot\hat{y}}$ such that

[TABLE]

Similarly since $u^{s}_{\varepsilon}$ is a radiating solution to the Helmholtz equation in $\mathbb{R}^{d}\setminus\overline{D}$, it assumes a similar asymptotic expansion as (3), and we denote by $u_{\varepsilon}^{\infty}(\hat{x},\hat{y})$ its corresponding far field pattern. In this study we assume that the media is non-absorbing, and $\inf_{x\in D}n(x)=n_{0}>0$, $n_{m}>0$, and

[TABLE]

For later use let us denote

[TABLE]

The inverse problem we consider here is to determine the location $\{z_{m}\}_{m=1,M}$ of the perturbations and information about $A_{m}$ and $n_{m}$ from knowledge of $u_{\varepsilon}^{\infty}(\hat{x},\hat{y})$ for several $\hat{x},\hat{y}\in{\mathbb{S}}$, provided that $A_{b}(x)$ and $n_{b}(x)$ are known.

In general, the support $D_{\epsilon}$ of the defects can be determined from the far field operator

[TABLE]

via the factorization method [8]. In addition, it is well-known [5] that the far field operator $F$ determines the real transmission eigenvalues which are defined below.

Definition 1.1**.**

Transmission eigenvalues are the values $k_{\varepsilon}\in\mathbb{C}$ for which there is a non-trivial solution $(w,v)\in H^{1}(D)\times H^{1}(D)$ of

[TABLE]

These transmission eigenvalues can be used to obtain information about $A_{\epsilon}$ and $n_{\epsilon}$ [5]. In this paper we will make use of the small volume feature of the defects and use asymptotic analysis to determine the locations $z_{m}$, $m=1\dots M$ of the small inhomogeneities. Then, based on the perturbation formulas of the transmission eigenvalues, we determine geometric and physical information about these small defects via polarization tensors in the asymptotic formulas.

2 Asymptotic Formulas and the MUSIC Algorithm

To avoid technical difficulties with asymptotic expansions, without loss of generality we assume that the anisotropic media is homogeneous, i.e. the matrix $A$ and the scalar $n$ are constant. We derive the multi-static response matrix by exploiting the fact that the each of the defective regions has small volume as in [21], which will be used to reconstruct the defective regions. The multi-static response matrix can be seen as the discrete version of the far field operator $F$ given by (8). To this end we first recall $\mathbb{G}(\cdot,\cdot)$ the Green’s function for the background layered media, i.e. the solution of

[TABLE]

Let $\mathbb{G}^{\infty}(\cdot\,,z)\in L^{2}(\mathbb{S})$ be it’s far field pattern. Since $A(x)$ is a symmetric constant positive definite matrix for $x\in D$ and $n(x)$ is a positive constant for $x\in D$ we have by Theorem 5.1 in [8] that

[TABLE]

It can be shown (see [2]) by using Green’s identities and the Sommerfeld radiation condition that $u_{\varepsilon}$ satisfies the Lippmann-Schwinger representation formula given by

[TABLE]

By linearity it is clear that the scattered field $u^{s}=u_{\varepsilon}^{s}-u_{b}^{s}$ is due to the defective regions $D_{\varepsilon}$ and from (12) an asymptotic expansion for $u^{\infty}(\hat{x},\hat{y})$ can be obtained by combining the asymptotic results from [2] and [17] together with (11). Thus we obtain

[TABLE]

where the polarization tensor ${\bf M}^{(m)}$ is given by

[TABLE]

with $e_{i}$ being the $i$th basis vector in $\mathbb{R}^{d}$ and $\phi_{i}$ is the solution to

[TABLE]

We now wish to use the leading term in (13) to determine the location of the defective regions $z_{m}+\varepsilon B_{m}$. To this end assume that there are $N$ incident and observation directions given by $\hat{y}_{j},\hat{x}_{i}\in{\mathbb{S}}$ for $i,j=1,2,\dots N$. Now we define the multi-static response matrix ${\bf F}\in\mathbb{C}^{N\times N}$ given by

[TABLE]

Remark 2.1**.**

The asymptotic expansion in (13) as well as the multi-static response matrix ${\bf F}$ given by (14) can be constructed for the more general case of an inhomogeneous background (i.e. where $A(x)$ is a matrix valued function and $n(x)$ is a scalar function in $D$). In this case $A$ and $n$ are replaced by $A(z_{m})$ and $n(z_{m})$. A general mixed reciprocity for inhomogeneous media is proven in [6].

The range of multi-static response matrix ${\bf F}$ determines the location of small inhomogeneities assuming that the data is collected at sufficiently many directions (see [1],[14] and references therein). In particular, we define the vectors ${\bf g}_{z}\in\mathbb{C}^{N}$ and ${\bf g}_{z,b}\in\mathbb{C}^{N}$ for any point $z\in\mathbb{R}^{d}$ and $b\neq 0\in\mathbb{C}^{d}$ by

[TABLE]

[TABLE]

and let

[TABLE]

Then the following range test can be proven (see [1],[14] and references therein), which essentially says that $z\in\{z_{m}\,:\,m=1,\dots,M\}$ if and only if ${\bf g}_{z,(1,b)}$ is in the range of ${\bf F}{\bf F}^{*}$.

Theorem 2.1**.**

Let ${\bf w}_{j}$ be the $j$-th orthonormal eigenvector of ${\bf F}{\bf F}^{*}$ and let r={\mbox{\em Rank}}\big{(}{\bf F}{\bf F}^{*}\big{)}. Assume that the set $S=\{{\hat{x}}_{i}\,:\,i\in\mathbb{N}\}$ is dense in $\mathbb{S}$ such that any analytic function that vanishes on $S$ also vanishes on $\mathbb{S}$. If $z\in D$ then there is a number $N_{0}\in\mathbb{N}$ such that for all $N\geq N_{0}$ we have that

[TABLE]

Numerical Validation of the MUSIC Algorithm

We present here the numerical implementation of the MUSIC algorithm in the ${\mathbb{R}}^{2}$ case. To this end, we use simulated far-field data to reconstruct the defects in a square scatterer. The simulated data comes from solving the direct scattering problems (1)-(2) and (4)-(5) using a cubic finite element method with a perfectly matched layer. From this we will have the approximated scattered fields $u^{s}_{\varepsilon}(\cdot\,,\hat{x})$ and $u^{s}_{b}(\cdot\,,\hat{x})$. The multi-static response matrix will be defined as

[TABLE]

where the far-field patterns are given by the solutions of the direct problems using the finite element method. In the following we use $N$ different directions on the unit circle given by

[TABLE]

In all examples we take $D=[-2,2]^{2}$ and we fix the wave number $k=1$.

We want to illustrate the performance of the MUSIC algorithm in reconstructing the defective regions $z_{m}+\varepsilon B_{m}$ inside $D$. We give examples with random noise added to the simulated data for $u^{\infty}_{\varepsilon}(\hat{x}_{i},\hat{x}_{j})$. The random noise level is given by $\delta$ where the noise is added to the far-field data $u^{\infty}_{\varepsilon}(\hat{x}_{i},\hat{x}_{j})+\delta E_{i,j}$ and the random matrix $E$ is such that $\|E\|_{2}=1$. Since we have that $A$, $\,n$ and $D$ are known for non-destructive testing we can assume that the far-field pattern $u^{\infty}_{b}(\hat{x}_{i},\hat{x}_{j})$ is computed from solving the scattering problem for the known background. Reconstruction examples are presented in Figure 1, Figure 2 and Figure 3 where the configuration and reconstruction parameters are explained in the respective labels.

3 Convergence of the Transmission Eigenvalue Spectrum

Recall that the multi-static data used for the MUSIC algorithm can also determine the transmission eigenvalues corresponding to the perturbed media. Having reconstructed the location of the small defects, we would like to obtain information about the strength of the perturbations $A_{m}$ and $n_{m}$ from the transmission eigenvalues. To this end, we investigate how the small defects affect the transmission eigenvalues. For the convergence analysis of the transmission eigenvalues we assume more regularity on the coefficients of the unperturbed (without defects) media, i.e. they are given by the symmetric matrix $A(x)\in{C}^{2}(D,\mathbb{R}^{d\times d})$ and $n(x)\in{C}^{1}(D)$. We start by showing that the eigenvalues $k_{\varepsilon}$ for the perturbed media converge to the eigenvalues for the unperturbed media as $\varepsilon\rightarrow 0$. Then we derive an asymptotic formula with correction term of the first order that can be used to obtain more information about the small defects. To analyze (9)-(10) we define the variational space X(D):=\Big{\{}(w,v):\;w,v\in H^{1}(D)\,|\,w-v\in H^{1}_{0}(D)\Big{\}} equipped with the $H^{1}(D)\times H^{1}(D)$ inner product. It is clear that the variational form of (9)-(10) is given by

[TABLE]

For convenience we define the bounded sesquilinear forms

[TABLE]

[TABLE]

[TABLE]

Therefore we have that (9)-(10) can be written as for all $(\varphi_{1},\varphi_{2})\in X(D)$

[TABLE]

Let us define by ${\bf A}_{\epsilon},\,{\bf B}_{\epsilon}$ and ${\bf C}:X(D)\to X(D)$ the bounded linear operators defined from \mathcal{A}_{\varepsilon}\big{(}\cdot\,;\cdot\big{)}, \mathcal{B}_{\varepsilon}\big{(}\cdot\,;\cdot\big{)} and \mathcal{C}\big{(}\cdot\,;\cdot\big{)} by means of the Riesz representation theorem. The unperturbed media corresponds to $\varepsilon=0$, where $A_{0}:=A$, $n_{0}:=n$ and $k_{0}:=k$. It can be shown using $\mathbb{T}$-coercivity that if $A_{min}$ and $a_{min}>1$ that ${\bf A}_{\epsilon}$ is invertible with the norm of the inverse independent of $\varepsilon\geq 0$. To this end we consider the isomorphism $\mathbb{T}(w,v)=(w,-v+2w):X(D)\mapsto X(D)$ (it is easy to check that $\mathbb{T}=\mathbb{T}^{-1}$). Then

[TABLE]

[TABLE]

and by Young’s inequality we obtain that

[TABLE]

where we let $\alpha=\min\{A_{min},\,a_{min}\}$. Therefore we have proven that \mathcal{A}_{\varepsilon}\big{(}(w,v);\mathbb{T}(w,v)\big{)} is coercive provided that $\delta\in(1/\alpha,\,1)$, implying ${\bf A}_{\epsilon}$ is invertible for $\epsilon\geq 0$. Similar arguments hold for $A_{max}$ and $a_{max}<1$ where $A_{min}$ is replaced by $A_{max}$ in \mathcal{A}_{\varepsilon}\big{(}\cdot\,;\cdot\big{)} and \mathcal{C}\big{(}\cdot\,;\cdot\big{)} with $\mathbb{T}(w,v)=(w-2v,-v)$. It is clear that in either case that ${\bf B}_{\epsilon}$ and ${\bf C}$ are compact operators by appealing to the compact embedding of $H^{1}(D)$ in $L^{2}(D)$. Now by (15) it is clear that $(w,v)$ are eigenfunctions corresponding to the eigenvalue $k_{\varepsilon}$ provided that

[TABLE]

Let us denote the eigenvalue parameter $\tau_{\varepsilon}=k^{2}_{\varepsilon}$ and define ${\bf T}_{\varepsilon}:X(D)\to X(D)$

[TABLE]

We can now rephrase (15) as a non-linear eigenvalue problem

[TABLE]

Note that it is clear that ${\bf T}_{\varepsilon}(\tau)$ depends analytically on $\tau$ in any subset of the complex plane that does not include the origin.

Convergence of the Spectrum

In this section, we study the convergence of ${\bf T}_{\varepsilon}(\tau)$ in the operator norm to the unperturbed operator ${\bf T}_{0}(\tau)$ and then use results from [22] to prove convergence for the transmission eigenvalues and eigenfunctions. To this end, notice that ${\bf B}_{\varepsilon}$ and ${\bf C}$ are compact operators and that $\|{\bf A}_{\varepsilon}^{-1}\|$ is uniformly bounded with respect to $\epsilon$, so we can conclude that ${\bf T}_{\varepsilon}(\tau)$ is compact for all $\varepsilon\geq 0$. Hence the convergence of ${\bf T}_{\varepsilon}(\tau)$ would then imply the convergence of the transmission eigenvalues. We start by studying the convergence of the operator ${\bf B}_{\varepsilon}$ to ${\bf B}_{0}$.

Theorem 3.1**.**

${\bf B}_{\epsilon}\rightarrow{\bf B}_{0}$* in the operator norm. Moreover for some $\alpha\in(0,1)$ we have that $\left\|{\bf B}_{\varepsilon}-{\bf B}_{0}\right\|\leq C\varepsilon^{\alpha}\,\text{ in }\,\mathbb{R}^{d}$ for some C independent of $\varepsilon$, $d=2,3$. *

Proof.

By definition we have that

[TABLE]

Therefore, we have that \left\|\big{(}{\bf B}_{\varepsilon}-{\bf B}_{0}\big{)}(w,v)\right\|_{X(D)}\leq||(n_{\varepsilon}-n)w||_{L^{2}(D_{\varepsilon})}. Now since $w\in H^{1}(D)$ we have from Sobolev’s embedding in $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ that $w\in L^{p}(D)$ for some $p\geq 2$ (see e.g. [4] for embedding results). We then conclude that $|w|^{2}\in L^{p/2}(D)$. Now let $q$ be defined by $\frac{1}{p/2}+\frac{1}{q}=1$ notice that $\frac{1}{q}=\frac{p-2}{p}$. Therefore by using the duality between $L^{p/2}(D)$ and $L^{q}(D)$ along with Sobolev’s embedding we have that

[TABLE]

Hence, we have that

[TABLE]

where the constant $C$ incorporates the norm of the contrasts but is independent of $\varepsilon$. Now for the $\mathbb{R}^{2}$ for any choice of $p>2$ we have that $\frac{1}{q}<1$ giving the result. For the case in $\mathbb{R}^{3}$ we can choose $p<6$ giving that $\frac{1}{q}<2/3$ and therefore $d/2q<1$, which gives the result in $\mathbb{R}^{3}$. ∎

We are now interested in the convergence of ${\bf A}_{\varepsilon}^{-1}{\bf B}_{\varepsilon}$ and ${\bf A}_{\varepsilon}^{-1}{\bf C}$ as $\varepsilon$ tends to zero. Recall that ${\bf A}_{\varepsilon}^{-1}$ exists as a bounded linear operator for all $\varepsilon\geq 0$ where the norm of ${\bf A}_{\varepsilon}^{-1}$ is uniformly bounded with respect to $\varepsilon$. To study the convergence of ${\bf A}_{\varepsilon}^{-1}{\bf B}_{\varepsilon}$ and ${\bf A}_{\varepsilon}^{-1}{\bf C}$ we first need some regularity results pertaining to ${\bf B}_{0}$ and ${\bf C}$. Notice that by the variational definition of ${\bf B}_{0}$ we have that for any $(f,g)\in X(D)$ if we denote ${\bf B}_{0}(f,g)=(w,v)$ then

[TABLE]

Therefore by elliptic regularity we have that $w$ and $v$ are in $H^{3}_{loc}(D)$ provided that $n$ is continuously differentiable, and for any $\Omega\subset D$

[TABLE]

Next, as for the operator ${\bf C}$ we have that for any $(f,g)\in X(D)$ if we denote ${\bf C}(f,g)=(w,v)$ then

[TABLE]

and we have the elliptic regularity estimates for any $\Omega\subset D$

[TABLE]

Theorem 3.2**.**

We have that

[TABLE]

in the operator norm as $\varepsilon\rightarrow 0$.

Proof.

Consider the pair $(w_{\varepsilon},v_{\varepsilon})$ and $(w,v)$ in $X(D)$ defined by $(w_{\varepsilon},v_{\varepsilon})={\bf A}_{\varepsilon}^{-1}(f,g)$ and $(w,v)={\bf A}_{0}^{-1}(f,g)$ for any $(f,g)\in X(D)$. By definition we have that

[TABLE]

whence using the $\mathbb{T}$-coercivity we conclude that

[TABLE]

Next we have that ${\bf A}_{0}^{-1}{\bf B}_{0}(f,g)=(w,v)$ due to the variational form of ${\bf A}_{0}$ satisfies

[TABLE]

Recalling ${\bf B}_{0}(f,g)=(p,q)$ we have $\|p\|_{H^{3}(\Omega)}+\|q\|_{H^{3}(\Omega)}\leq C\left(\|f\|_{H^{1}(D)}+\|g\|_{H^{1}(D)}\right)$ where $\Omega\subset D$ and by elliptic regularity given any $\Omega^{\prime}\subset\Omega\subset D$ we have that

[TABLE]

Fixing $\Omega^{\prime}$ and $\Omega$ such that $D_{\varepsilon}\subset\Omega^{\prime}\subset\Omega\subset D$ for all $\varepsilon$ sufficiently small and using that $H^{3}(\Omega^{\prime})\subset C^{1}(\Omega^{\prime})$ we have the following estimates

[TABLE]

Now appealing to the continuity of the embedding of $H^{3}(\Omega^{\prime})$ into $C^{1}(\Omega^{\prime})$ and the regularity estimate we have that

[TABLE]

Using that

[TABLE]

along with the uniform boundedness of $||{\bf A}_{\varepsilon}^{-1}||$ and the norm convergence of ${\bf B}_{\varepsilon}$ to ${\bf B}_{0}$ implies that ${\bf A}_{\varepsilon}^{-1}{\bf B}_{\varepsilon}\rightarrow{\bf A}_{0}^{-1}{\bf B}_{0}$ in norm. The same arguments work for showing that ${\bf A}_{\varepsilon}^{-1}{\bf C}\rightarrow{\bf A}_{0}^{-1}{\bf C}$ in norm, which ends the proof. ∎

Corollary 3.1**.**

Let the operators ${\bf A}_{\varepsilon}$, ${\bf A}_{0}$, ${\bf B}_{\varepsilon}$, ${\bf B}_{0}$ and ${\bf C}$ be defined by the variational forms given above. Then we have that for $d=2,3$

[TABLE]

and \left\|{\bf A}_{\varepsilon}^{-1}\big{(}{\bf B}_{\varepsilon}-{\bf B}_{0}\big{)}\right\|=\mathcal{O}(\varepsilon^{\alpha}) for some $\alpha\in(0,1)$ .

Combining the above results we have:

Theorem 3.3**.**

Let the operator ${\bf T}_{\varepsilon}(\tau)$ be as defined in (17) and $\tau\in U$ with $U$ being any bounded subset of $\mathbb{C}$ with zero not a limit point of $U$. Then we have that

[TABLE]

Moreover if $n_{\varepsilon}=n$ for all $\varepsilon\geq 0$ then we have that

[TABLE]

Having proven the convergence of the operator ${\bf T}_{\varepsilon}(\tau)$ we are ready to study the convergence of the real transmission eigenvalues using the abstract result from [22].

Lemma 3.1**.**

Let $\tau$ be a non-linear eigenvalue of ${\bf T}_{0}$ and assume that ${\bf T}_{0}(\tau)$ and ${\bf T}_{\varepsilon}(\tau)$ are both meromorphic in some region $U$ of $\mathbb{C}$ containing $\tau$. Also assume that ${\bf T}_{\varepsilon}(\tau)\rightarrow{\bf T}_{0}(\tau)$ in the operator norm. Then for any ball around $\tau$ there exists a $\varepsilon_{0}>0$ such that ${\bf T}_{\varepsilon}$ has a non-linear eigenvalue in the ball for all $\varepsilon<\varepsilon_{0}$. Conversely if $\tau_{\varepsilon}$ is a sequence of non-linear eigenvalues of ${\bf T}_{\varepsilon}$ that converges as $\varepsilon\rightarrow 0$, then the limit $\tau$ is a non-linear eigenvalue of ${\bf T}_{0}$.

By Theorem 3.3 we have that ${\bf T}_{\varepsilon}(\tau)\rightarrow{\bf T}_{0}(\tau)$ in the operator norm in any in region $U$ of $\mathbb{C}\setminus\{0\}$ and from the definition of the operator ${\bf T}_{\varepsilon}(\tau)$ we have that it depends analytically on $\tau$ in any subset of the complex plane that does not include the origin. Finally to conclude the convergence of the eigenvalues, we need bounds on the eigenvalues independent of $\epsilon$. The existence of real transmission eigenvalues and monotonicity property with respect to the refractive index are proven in [7] and [9]. The monotonicity property implies $\epsilon$-independent bounds on these real transmission eigenvalues since $A_{\epsilon}$ and $n_{\epsilon}$ are bounded above and below uniformly with respect to $\epsilon>0$ (more specifically such bounds can be obtained by modifying the proof of Theorem 2.6 and Theorem 2.10 in [9] in a similar way as in the proof of Corollary 2.6 in [7].)

4 Asymptotic Formula for the Transmission Eigenvalues

Having proven the convergence of the transmission eigenvalues, we now want to obtain an asymptotic formula for the real transmission eigenvalues. To this end, we need to construct an appropriate corrector that will give an explicit formula for the first term in the asymptotic expansion for the transmission eigenvalues. For technical reasons that have to do with the rate of convergence of ${\bf B}_{\epsilon}$ to ${\bf B}_{0}$ (which will be explained later on) we derive this corrector for the case when there is no contrast in the lower term, i.e. $n_{\epsilon}=n$. To avoid technicalities in the presentation, the corrector will be derived for a homogeneous anisotropic media and the results can be generalized for an inhomogeneous media as in [11]. Hence in this section we again assume that the coefficients $A$ and $n$ are constant in $D$.

Correction for the Operator ${\bf A}_{\varepsilon}^{-1}-{\bf A}_{0}^{-1}$

Consider the pair $(w_{\varepsilon},v_{\varepsilon})$ and $(w,v)$ in $X(D)$ defined by

[TABLE]

and we assume that $w$ is a smooth function. Without loss of generality in the following we perform the calculations only for one inhomogeneity. For multiple inhomogeneity one simply sum the correctors. To this end, assume that the defective region is of the form $\varepsilon B$ where $B$ is centered at the origin with constant matrix $A_{1}$ being the constitutive parameter. We make the scaling $y=x/\varepsilon$ and $\widetilde{D}=\frac{1}{\varepsilon}D$ and let $w_{\varepsilon}^{(1)}(y)\in H^{1}_{0}(\widetilde{D})$ be the unique solution to

[TABLE]

with $\widetilde{A}=A_{1}\,\chi_{B}+A(1-\chi_{B})$.

Theorem 4.1**.**

Assume that $(w_{\varepsilon},v_{\varepsilon})$ and $(w,v)$ are defined by (20) with $w$ being a smooth function, then we have that

[TABLE]

Proof.

Recall that $x=\varepsilon y$ and we define the error functions in $X(\widetilde{D})$ \big{(}note that w_{\varepsilon}^{(1)}(y)\in H^{1}_{0}(\widetilde{D})$$\big{)}

[TABLE]

Now let $(\varphi_{1},\varphi_{2})\in{X(\tilde{D})}$ and define the sesquilinear form

[TABLE]

Using (20) we have that

[TABLE]

Using integration by parts and (21) gives that

[TABLE]

Recall that $w$ is smooth, therefore $\nabla_{x}\cdot(A-A_{1})\nabla_{x}w(\varepsilon y)$ is bounded in $B$. Also notice that by Taylor’s expansion we have that the term $(\nabla_{x}w(\varepsilon y)-\nabla_{x}w(0))=\mathcal{O}(\varepsilon)$. Therefore, we can conclude that there is a constant $C$ independent of $\varepsilon$ such that

[TABLE]

Using the $\mathbb{T}$-coercivity of the sesquilinear form \tilde{\mathcal{A}}_{\varepsilon}\big{(}\cdot\,;\cdot\big{)} in $X(\widetilde{D})$ gives that

[TABLE]

and the result follows from scaling. ∎

Notice that from (21) we have that $\|w^{(1)}_{\varepsilon}(y)\|_{H^{1}(\widetilde{D})}$ is bounded independently of $\varepsilon$ by the Lax-Milgram lemma. Therefore by scaling we have that $\|w^{(1)}_{\varepsilon}(x/\varepsilon)\|_{H^{1}(D)}\leq C\varepsilon^{d/2-1}$ with $C$ independent of $\varepsilon$, which gives the following result.

Corollary 4.1**.**

Assume that $(w_{\varepsilon},v_{\varepsilon})$ and $(w,v)$ are defined by (20) with $w$ being a smooth function then we have that

[TABLE]

Notice that the corrector $w^{(1)}_{\varepsilon}(y)$ depends on $\varepsilon$, hence we now wish to construct a corrector that is independent of the small parameter $\varepsilon$. To this end, we define the function $w^{(1)}(y)\in H^{1}(\mathbb{R}^{d})$ such that for all $\varphi\in H^{1}(\mathbb{R}^{d})$

[TABLE]

Note that the variational problem (25) implies that

[TABLE]

We now have $\left|w^{(1)}\right|\rightarrow 0$ as $|y|\rightarrow\infty$, exponentially fast [12]. This gives that $\nabla_{y}w^{(1)}$ decays faster than the gradient of a solution to Laplace’s equation, therefore

[TABLE]

Theorem 4.2**.**

Let $w^{(1)}_{\varepsilon}$ and $w^{(1)}$ be defined as the solutions to (21) and (25) respectively, then we have that

[TABLE]

Proof.

Let $u_{\varepsilon}=w^{(1)}_{\varepsilon}(y)-w^{(1)}(y)$, there exists a constant $\alpha>0$ such that

[TABLE]

Notice that the variational form (25) implies that

[TABLE]

Therefore integration by parts gives that

[TABLE]

Now by using the boundary value problem for $w^{(1)}$ we have that

[TABLE]

[TABLE]

By the scaling we have that

[TABLE]

Since

[TABLE]

we can conclude that

[TABLE]

which gives the result by scaling the norm back to the domain $D$. ∎

By appealing to the triangle inequality we have the following result.

Corollary 4.2**.**

Let $w^{(1)}$ be the solutions to (25), also assume that $(w_{\varepsilon},v_{\varepsilon})$ and $(w,v)$ are defined by (20) with $w$ being a smooth function then we have that

[TABLE]

The arguments used in this section carry over to the case of multiple inhomogeneities. Indeed, for multiple inhomogeneities centered at $z_{m}$ with anisotropic material parameter $A_{m}$ we have that by using translation and summing over a finite number of inhomogeneities gives that the corrector takes the form

[TABLE]

where $w^{(1)}_{m}(x/\varepsilon)$ is the solution to

[TABLE]

for all $\varphi\in H^{1}(\mathbb{R}^{d})$ with $\widetilde{A}_{m}=A_{m}\,\chi_{B_{m}}+A(1-\chi_{B_{m}})$. The convergence results in this section still hold for $w(0)w^{(1)}(x/\varepsilon)$ replaced by $\widetilde{w}^{(1)}(x/\varepsilon)$.

Asymptotic Formulas

Finally we have all the ingredients to give an asymptotic formula for the transmission eigenvalues using the results in [22]. Note that we have assumed that contrast in the defect is only in the matrix valued material parameter (i.e. $n_{\varepsilon}=n$ for all $\varepsilon>0$), and we still take $A$ and $A_{m}$ constant matrices. Under this assumption we have that the operator ${\bf T}_{\varepsilon}(\tau)={\bf A}_{\varepsilon}^{-1}{\bf B}_{0}+\frac{1}{\tau}{\bf A}_{\varepsilon}^{-1}{\bf C}$ converges in the operator norm.

We now recall Theorem 4.1 of [22] which is a generalization of Osborn’s Theorem (see [20] for Osborn’s result) to nonlinear eigenvalue problems.

Theorem 4.3**.**

Let $X$ be a Hilbert space and ${\bf T}_{\varepsilon}(\tau):X\to X$ be a compact operator valued functions of $\tau$ which are analytic in a region $U$ of the complex plane, such that $\left\|{\bf T}_{\varepsilon}(\tau)-{\bf T}_{0}(\tau)\right\|\rightarrow 0$ for all $\tau\in U$. Now assume that $\tau$ is a simple nonlinear eigenvalue of ${\bf T}_{0}(\tau)$ with normalized eigenfunction $\phi$. Then if

[TABLE]

we have that

[TABLE]

[TABLE]

with $\tau_{\varepsilon}$ is a nonlinear eigenvalue for ${\bf T}_{\varepsilon}(\tau)$.

Theorem 4.3 only holds for simple eigenvalues. Notice that we have established the order of convergence of the operator defined by the transmission eigenvalue problem. In particular, the results in the previous section (see equation (24)) gives that

[TABLE]

We now consider the point wise convergence for the adjoint operator.

Lemma 4.1**.**

Let $(w_{\tau},v_{\tau})\in X(D)$ be the smooth eigenfunction corresponding to the eigenvalue $\tau$ of the operator ${\bf T}_{0}(\tau)$, then we have that

[TABLE]

Proof.

Notice that ${\bf T}^{*}_{\varepsilon}(\tau)={\bf B}_{0}{\bf A}_{\varepsilon}^{-1}+\frac{1}{\tau}{\bf C}{\bf A}_{\varepsilon}^{-1}$ where we define $(w,v)={\bf A}_{0}^{-1}(w_{\tau},v_{\tau})$ and $(w_{\varepsilon},v_{\varepsilon})={\bf A}_{\varepsilon}^{-1}(w_{\tau},v_{\tau})$. Now for any $(\varphi_{1},\varphi_{2})\in X(D)$

[TABLE]

Since the sesqulinear form $\mathcal{B}_{0}$ only has $L^{2}(D)$ terms, we have that

[TABLE]

By rescaling the $L^{2}$ norm in equation (23) gives that

[TABLE]

therefore $\left\|{\bf B}_{0}({\bf A}_{\varepsilon}^{-1}-{\bf A}_{0}^{-1})(w_{\tau},v_{\tau})\right\|_{X(D)}=\mathcal{O}(\varepsilon^{d/2+1})$. A similar argument gives that $\left\|{\bf C}({\bf A}_{\varepsilon}^{-1}-{\bf A}_{0}^{-1})(w_{\tau},v_{\tau})\right\|_{X(D)}=\mathcal{O}(\varepsilon^{d/2+1})$, proving that claim. ∎

Remark 4.1**.**

This result shows why the case where $n_{\varepsilon}\neq n$ can not be handled by this analytic framework. In particular, the rate of convergence in Theorem 3.1 for ${\bf B}_{0}-{\bf B}_{\varepsilon}$ is not fast enough to provide an improved convergence rate for ${\bf T}^{*}_{\varepsilon}(\tau)-{\bf T}^{*}_{0}(\tau)$ which is necessary to apply Theorem 4.3.

We have just shown that the remainder term for the non-linear eigenvalue corrector formula is of the order $\varepsilon^{d+1}$. To construct an asymptotic formula for the transmission eigenvalues we need to construct an asymptotic formula for

[TABLE]

where $(w_{\tau},v_{\tau})$ are the eigenfunctions for $\varepsilon=0$. By equation (15) we have that ${\bf B}_{0}(w,v)+\frac{1}{\tau}{\bf C}(w,v)=\frac{1}{\tau}{\bf A}_{0}(w,v)$. Since the operator ${\bf A}_{\varepsilon}$ is self-adjoint for all $\varepsilon\geq 0$ the definition of ${\bf T}_{\varepsilon}(\tau)$ in (17) gives that

[TABLE]

[TABLE]

This gives that we only need to construct an asymptotic formula for

[TABLE]

We now derive an asymptotic formula for ${\bf A}_{\varepsilon}^{-1}-{\bf A}_{0}^{-1}$ with respect to the sesquilinear form \mathcal{A}_{0}\big{(}\cdot\,;\cdot\big{)} which is given in the following result.

Theorem 4.4**.**

Let $(w_{\tau},v_{\tau})$ be the eigenfunctions for $\varepsilon=0$ with transmission eigenvalue $\tau$ and define $(w,v)={\bf A}_{0}^{-1}(w_{\tau},v_{\tau})$, then we have that

[TABLE]

Proof.

We will prove the result for a single defect centered at the origin then by using translation and summing a finite number of such inhomogeneities, the asymptotic result follows. Letting $(w_{\varepsilon},v_{\varepsilon})={\bf A}_{\varepsilon}^{-1}(w_{\tau},v_{\tau})$, we have that

[TABLE]

Recall, that by elliptic regularity we have the for any $\Omega$ such that $\varepsilon B\subset\Omega\subset D$ the eigenfunctions are in $C^{1}(\Omega)$. Using this along with the support of $A-A_{\varepsilon}$ and Corollary 4.2 we can now estimate the first term

[TABLE]

We now consider the second term of (27) which is given by

[TABLE]

where we have used Taylor’s expansion about the origin to estimate the first integral. Now by the divergence theorem we have that the volume integral involving the eigenfunction and the corrector is given by

[TABLE]

Now by rescaling the second integral for $x=\varepsilon y$ and using a Taylor’s expansion we have that integration is given by

[TABLE]

proving the result. ∎

Now we have all we need for an asymptotic formula for simple transmission eigenvalues. Notice that $\frac{d}{d\tau}{\bf T}_{0}(\tau)=-\frac{1}{\tau^{2}}{\bf A}_{0}^{-1}{\bf C}$, therefore we have that

[TABLE]

For convenience let the constant

[TABLE]

Therefore we have that simple transmission eigenvalues have the expansion.

Theorem 4.5**.**

Let $(w_{\tau},v_{\tau})$ be the eigenfunctions for $\varepsilon=0$ with simple transmission eigenvalue $\tau$ and define $(w,v)={\bf A}_{0}^{-1}(w_{\tau},v_{\tau})$, then we have that

[TABLE]

where $q_{m}$ is given by (28) and

[TABLE]

Numerical Validation of the Asymptotic Formula

The asymptotic formula given in Theorem 4.5 can potentially be used to determine the strength of the small defective region(s). Notice that the MUSIC algorithm gives the location of the defect(s) and recall that the transmission eigenvalues for the perturbed media $\tau_{\varepsilon}$ can be measured from the same scattering data needed for MUSIC but for a range of wave numbers $k$ (see [3], [5], [15]), whereas the transmission eigenvalues $\tau$ and eigenfunctions $(w_{\tau},v_{\tau})$ for the unperturbed media can be computed since $A$ and $n$ are assumed to be known. In particular, denoting by ${\mathbf{F}}_{\epsilon}:=\left(u_{\epsilon}(\hat{x}_{i},x_{j},k)\right)_{i,j=1..N}$, the far field matrix due to the inhomogeneity $D$ with perturbation $D_{\epsilon}$, where we indicate its dependence on $k$, to determine the $\tau_{\epsilon}:=k_{\epsilon}^{2}$ we solve the regularized equation

[TABLE]

for ${\mathbf{g}}(k)$ for a range of $k$. The transmission eigenvalues are those values of $k$ for which $\|{\mathbf{g}}(k)\|_{\ell^{2}}$ blows up. To compute the transmission eigenvalues for the unperturbed media we use a continuous finite element method with the eigenvalue searching technique described in [15], [18] and [19]). In order to use the asymptotic formula in Theorem 4.5 one also needs the functions $(w,v)={\bf A}_{0}^{-1}(w_{\tau},v_{\tau})$ which can be solved for (e.g. using the FEM) since $A$ and $(w_{\tau},v_{\tau})$ are known. Having identified the location of the defect(s) from the MUSIC algorithm (i.e. the points $z_{m}$ are known) one can determine the strength of the defect(s) which is given by $(A_{m}-A)|B_{m}|$ and $q_{m}$ from the knowledge of two transmission eigenvalues. Notice, that the strength of the defect(s) only depend on the constitutive coefficients and geometry of the defect(s).

We first consider a few examples to illustrate the convergence of the transmission eigenvalues as $\varepsilon\rightarrow 0$ in $\mathbb{R}^{2}$. To do so, we denote the transmission eigenvalues for the unperturbed media by $k^{2}_{j}$ and the first transmission eigenvalue for the perturbed media by $k^{2}_{j}(\varepsilon)$. To test our asymptotic formula we will check the order of convergence for two transmission eigenvalues. We compute the error and estimated order of convergence by

[TABLE]

In our calculations we see that the order of convergence seems to be approximately second order which is what is predicted by Theorem 4.5

Example 1. Here we let $D=[-1,1]^{2}$ where $D_{0}$ is given by two disks of radius $\varepsilon$ centered at $(0.25,0)$ and $(-0.25,-0.25)$. For this case we take $n=n_{\varepsilon}=1$ for all $\varepsilon$ with

[TABLE]

and $A_{\varepsilon}=2I$. Below in Table 1 we show estimated order of convergence for two transmission eigenvalues.

Example 2. Here we let ${\displaystyle D=\left\{(x_{1},x_{2})\in\mathbb{R}^{2}\,:\,\frac{x_{1}^{2}}{4}+x_{2}^{2}<1\right\}}$ where $D_{0}$ is the disk centered at the origin of radius $\varepsilon$. For this case we take $n=n_{\varepsilon}=1$ for all $\varepsilon$ with $A=10I$ and $A_{\varepsilon}=2I$. Below in Table 2 we show estimated order of convergence for two transmission eigenvalues.

Next we provide an example on how to use the asymptotic formula in Theorem 4.5 to obtain information about the strength of the perturbation. For the case of a homogeneous isotropic unit disc $D$ with $A=\alpha I$ with $\alpha$ a positive constant and $n=1$ we have that the radially symmetric eigenfunctions corresponding to the eigenvalues $k^{2}$ are given by

[TABLE]

Here $\text{J}_{0}$ is the first kind Bessel function of order zero.Using the variational formulation, the solution to ${\bf A}_{0}(w,v)=(w_{k},v_{k})$ can be show to be

[TABLE]

where the constants $c_{1}$ and $c_{2}$ satisfy

[TABLE]

where here $\text{I}_{0}$ is the third kind Bessel function of order zero. This implies that the corrector term in Theorem 4.5 is known up to the weighted contrast $(A_{m}-A)|B_{m}|$ and ‘polarization’ constant $q_{m}$. Assuming that two transmission eigenvalues for the unperturbed media are known and the corresponding two eigenvalues for the perturbed media is computed via the far-field data, by ignoring the $o(\varepsilon^{d})$ term in the asymptotic formula given in Theorem 4.5 one obtains a 2$\times$2 linear system of equations to determine the weighted contrast and polarization constant. For proof of concept, we consider the case where $D$ is the unit disk with $D_{0}$ being the disk of radius $\varepsilon$ centered at $(0.25,0)$ where coefficients are taken be $A=10I$ and $A_{\varepsilon}=2I$. We assume that for this configuration the center of the disc is reconstructed using MUSIC as shown in Figure 4.

Using the first two radially symmetric eigenvalues and functions we wish to determine the contrast for the particular case of $\varepsilon=1/2$. In this example the contrast is given by $-8$ and solving the 2$\times$2 linear system derived from the asymptotic formula recovers a contrast of $-7.3465$ if the exact location is used in the formula. Using the reconstructed center $(0.255,0.005)$ for the case without noise in Figure 4 we obtain that the contrast is $-7.1222$, and using the reconstructed center $(0.267,0.017)$ for the case with $10\%$ noise in Figure 4 we obtain that the contrast is $-6.6147$. This preliminary example shows that one can determine information about the location and material properties of the small defects from a knowledge of the far-field data. Of course further investigation is needed to numerically validate our imaging method.

Acknowledgments

The research of F. Cakoni is supported in part by AFOSR grant FA9550-17-1-0147, NSF Grant DMS-1602802 and Simons Foundation Award 392261. The research of S. Moskow is supported in part by NSF Grant DMS-1411721.