# Numerical reconstruction of electromagnetic inclusions in three dimensions

Gang Bao, Junshan Lin, Séraphin Mefire

► To cite this version:

Gang Bao, Junshan Lin, Séraphin Mefire. Numerical reconstruction of electromagnetic inclusions in three dimensions. SIAM Journal on Imaging Sciences, 2014, 7 (1), pp.558-577. 10.1137/130937640. hal-01286235

HAL Id: hal-01286235

https://hal.science/hal-01286235

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical reconstruction of electromagnetic inclusions in three dimensions

Gang Bao*, Junshan Lin†, and Séraphin M. Mefire‡

## Abstract

This paper is concerned with the reconstruction of electromagnetic inclusions in a three-dimensional bounded domain by boundary measurements. An accurate and stable reconstruction method is presented to solve the associated inverse problem. The approach consists of two main steps: (1). At low frequency, a Multiple SIgnal Classification (MUSIC) algorithm is used to obtain the locations of inclusions, which serve as an initial guess for reconstructions at higher frequencies; (2). A continuation method based on multiple frequency data is then applied to recover the shapes of inclusions accurately. Numerical examples are provided to illustrate the effectiveness of the approach.

**Keywords.** inverse problems, Maxwell equations, shape optimization

**MS Codes.** 65N21, 65N30, 78A25

## 1 Introduction

Electromagnetic wave has found significant applications in civil and military engineering. There has been long interest and extensive study on inverse electromagnetic scattering theory in different contexts, see for example [3, 17, 19]. In this paper, we consider the reconstruction of inhomogeneities embedded in a background medium by electromagnetic measurements on the boundary of the object. The problem arises naturally in applications such as medical imaging, non-destructive testing of materials, geophysics, etc.

Let $\varepsilon$ and $\mu$ be the electric permittivity and magnetic permeability of the medium respectively in a bounded domain $\Omega \subset \mathbb{R}^3$. Denote the permittivity and permeability of the vacuum as $\varepsilon_0$ and $\mu_0$. Then $\varepsilon$ and $\mu$ can be rewritten as $\varepsilon = \varepsilon_r\varepsilon_0$ and $\mu = \mu_r\mu_0$, where $\varepsilon_r$ and $\mu_r$ are the relative permittivity and permeability respectively. In this work, we are interested in the case when the permittivity $\varepsilon_r$ is piecewise constant inside the domain $\Omega$. To be more precise, let $D_j$ ($j = 1, 2, \dots, j_0$) be the disconnected inclusions inside the domain $\Omega$ such that $D_i \cap D_j = \emptyset$ if $i \neq j$, and $D = \cup_{j=1}^{j_0} D_j$ be the union of inclusions. Then the relative permittivity takes the following form:

$$ \varepsilon_r(x) = \begin{cases} \varepsilon_1 & \text{if } x \in D, \\ \varepsilon_2 & \text{if } x \in \Omega \setminus \overline{D}, \end{cases} $$

<footnote>*Department of Mathematics, Zhejiang University, Hangzhou, China; and Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA. Email: bao@math.msu.edu.</footnote>

<footnote>†Department of Mathematics and Statistics, Auburn University, Auburn, AL 36849, USA. Email: jzl0097@auburn.edu.</footnote>

<footnote>‡Institut Elie Cartan de Lorraine, CNRS UMR 7502, Fac. des Sciences et Tech. de Nancy, Université de Lorraine, 54506 Vandœuvre-lès-Nancy Cedex, France. Email: seraphin.mefire@univ-lorraine.fr.</footnote>
where $\varepsilon_1$ and $\varepsilon_2$ are permittivities for inhomogeneities and background medium respectively. Throughout the paper, we also restrict our attention to the nonmagnetic medium by assuming that $\mu_r \equiv 1$ in $\Omega$.

For a time harmonic electromagnetic wave (with $e^{-i\omega t}$ dependence), the electric field $E$ and magnetic field $H$ satisfy the Maxwell's equations in the bounded domain $\Omega$:

$$ \left\{ \begin{array}{l} \nabla \times E = i\omega\mu_0 H, \\ \nabla \times H = -i\omega\varepsilon_r\varepsilon_0 E, \end{array} \right. $$

where $\omega$ is the operating frequency of the electromagnetic wave. By taking the curl of the first equation, and substituting into the second equation, the equation for the electric field $E$ becomes

$$ \nabla \times \nabla \times E - k^2 \varepsilon_r E = 0 \quad \text{in } \Omega. \qquad (1.1) $$

Here $k = \omega\sqrt{\varepsilon_0\mu_0}$ denotes the wavenumber. Let $\nu$ be the outward unit normal on the boundary $\partial\Omega$ of the domain. Assume that some electric current is applied on the boundary, then the following boundary condition for the electric field is prescribed:

$$ E \times \nu = g \quad \text{on } \partial\Omega. \qquad (1.2) $$

From boundary measurements of the magnetic field $H \times \nu$, or equivalently $\nabla \times E \times \nu$, the inverse problem is to reconstruct the inclusions $D_j$, $j = 1, \dots, j_0$, or their boundaries $\partial D_j$. For clarity of exposition, let us denote the boundaries as $\Gamma := \cup_{j=1}^{j_0} \partial D_j$. The problem considered here has piecewise constant medium, where the reconstruction of the boundaries of inclusions is of interest. The readers are referred to [7] for a related inverse medium scattering problem in two dimensions, and [8, 9, 10] for the reconstruction of a continuous medium by electromagnetic measurements in three dimensions.

The goal of this paper is to study a two-stage method to reconstruct the inclusions in an accurate and stable way. In the first stage, we employ a Multiple Signal Classification (MUSIC) algorithm at low frequency to obtain the location of each inclusion $D_j$ ($j = 1, 2, \dots, j_0$). The locations then serve as an initial guess for the reconstruction process at higher frequencies. The MUSIC approach is based upon an asymptotic expansion of the magnetic field $H$ on the boundary and a perturbation formula at low frequency. Such a perturbation formula is initiated by Ammari et al for the inverse conductivity problem [4], and is further developed in [5] for the three-dimensional Maxwell's equations. From the perturbation formula, a MUSIC algorithm is proposed and studied in details in [6] to locate the electromagnetic inclusions.

In the second stage, in order to obtain the boundaries of inclusions, the nonlinear inverse problem is then recast as a shape optimization problem of solving $D$ such that the resulting magnetic field matches the boundary measurement. We propose to use multiple frequency data to solve this nonlinear optimization problem. Indeed, recent studies for related inverse scattering problems suggest that a recursive linearization (continuation) approach with multiple frequency data overcomes the difficulty of reaching some local minimum, and is promising to lead the solution to the global minimum of the cost functional [8, 9, 11, 12, 16]. In this paper, we extend the continuation approach for the considered nonlinear optimization problem when multiple frequency data are available. The optimization scheme uses the locations of inclusions obtained in the first stage as an initial guess, i.e., by guessing each $D_j$ as a sphere with a small radius centered at the obtained location $z_j$. Then the scheme marches from low to high frequencies. At a fixed frequency, using the solution obtained at the previous frequency as an initial guess, a steepest descent direction $V$ for the defined cost functional is derived on the surface $\Gamma$, and a level set approach is employed to simulate the dynamic evolution of the surface with given vector field $V$.
The rest of the paper is organized as follows. Section 2 begins with the perturbation formula and the MUSIC algorithm at low frequency to locate the inclusions. In Section 3, the reconstruction by the continuation approach with multiple frequency data is proposed. We also discuss the steepest descent direction at a fixed frequency for the cost functional and the dynamic evolution of the surface with such a direction by the level set method. Several numerical examples are presented in Section 4 to demonstrate the effectiveness of the method, and the paper is concluded with some general remarks in Section 5.

## 2 A MUSIC algorithm at low frequency

Assume that the typical size of each inclusion $D_j$ ($j = 1, 2, \dots, j_0$) is $\rho$. Moreover, $D_j$ is represented by $z_j + \rho U_j$, where $z_j$ denotes its location, and $U_j$ is a bounded domain containing the origin and has smooth boundary. The perturbation formula is based on the assumption that

$$ \begin{cases} d_0 \le |z_j - z_l| & \forall j \ne l, \\ d_0 \le \operatorname{dist}(z_j, \partial\Omega) & \forall j, \end{cases} $$

for some positive constant $d_0$. In addition, the wavenumber $k$ is such that the natural weak formulation of the Maxwell equations for the homogeneous medium

$$ \left\{ \begin{array}{ll} \nabla \times \nabla \times E_0 - k^2 E_0 = 0 & \text{in } \Omega, \\ E_0 \times \nu = g & \text{on } \partial\Omega, \end{array} \right. \qquad (2.1) $$

has a unique solution $E_0$. Let $W$ be any smooth solution of the Maxwell equation

$$ \nabla \times \nabla \times W - k^2 W = 0 \text{ in } \bar{\Omega}, $$

then there exists a constant $\rho_0$ that for $0 < \rho < \rho_0$, the following asymptotic formula holds (see [4]):

$$ \int_{\partial\Omega} \nabla \times E \times \nu \cdot W ds - \int_{\partial\Omega} \nabla \times W \times \nu \cdot (\nu \times (E \times \nu)) ds \\ = \rho^3 k^2 \sum_{j=1}^{j_0} \left(\frac{\varepsilon_2}{\varepsilon_1} - 1\right) \left(\mathbf{M}_j\left(\frac{\varepsilon_2}{\varepsilon_1}\right) E_0(z_j)\right) \cdot W(z_j) + O(\rho^4), \qquad (2.2) $$

with $E$ the unique weak solution of the boundary value problem (1.1) - (1.2).

Here the constant $\rho_0$ depends on the shapes $U_j$, $\Omega$, permittivity values $\varepsilon_1$ and $\varepsilon_2$, the wavenumber $k$ and $d_0$, but is otherwise independent of $g$, $W$, and of the points $z_j$, $j = 1, \dots, j_0$. Also, $\mathbf{M}_j$ is the polarization tensor (see [3] for detailed properties), where each entry is given by

$$ m_{pq}^j \left( \frac{\varepsilon_2}{\varepsilon_1} \right) = \frac{\varepsilon_1}{\varepsilon_2} \int_{U_j} \frac{\partial \phi_p}{\partial x_q} dx, \quad 1 \le p, q \le 3, $$

with the scalar potential $\phi_p$ satisfying

$$ \left\{ \begin{array}{ll} \Delta\phi_p = 0 & \text{in } \mathbb{R}^3 \setminus \overline{U_j}, \\ \Delta\phi_p = 0 & \text{in } U_j, \\ \phi_p^+ - \phi_p^- = 0 & \text{on } \partial U_j, \\ \frac{\varepsilon_2}{\varepsilon_1} \left(\frac{\partial\phi_p}{\partial\nu}\right)^+ - \left(\frac{\partial\phi_p}{\partial\nu}\right)^- = 0 & \text{on } \partial U_j, \\ \phi_p(x) - x_p \to 0 & \text{as } |x| \to \infty, \end{array} \right. $$
where the outward unit normal to $\partial(z_j + \rho U_j)$, the boundary of $z_j + \rho U_j$, is also denoted by $\nu$, and the superscripts $+$, $-$ indicate the limiting values as $\partial(z_j + \rho U_j)$ is approached from outside $z_j + \rho U_j$, and from inside $z_j + \rho U_j$. Also, the coordinates of $x \in \mathbb{R}^3$ are represented by $x_q$, $1 \le q \le 3$.

Now we introduce the MUSIC algorithm based on the asymptotic formula (2.2). The method is essentially an approach to characterize the range of an operator and is closely related to the linear sampling method [18, 22]. For simplicity of presentation, we assume that $\Omega$ is a unit sphere. For $p = 1, 2, \dots, m$, let

$$E_{0,p}(x) = \hat{\theta}_{p}^{\perp} e^{ik\hat{\theta}_{p} \cdot x},$$

where $\hat{\theta}_p$ is the incidence direction, and $\hat{\theta}_p^{\perp}$ is the polarization direction. It is clear that $E_{0,p}$ is the solution of the Maxwell's equation for the homogeneous medium (2.1) with current

$$g_p(x) = \hat{\theta}_p^{\perp} e^{ik\hat{\theta}_p \cdot x} \times \nu$$

on the boundary $\partial\Omega$.

Let $E_p$ be the solution of the Maxwell equation (1.1) with the applied current $g = g_p$. Then for each test vector function

$$W_q(x) = \hat{\theta}_q^{\perp} e^{ik\hat{\theta}_q \cdot x}, \quad q = 1, 2, \dots, m,$$

the left-hand side of the asymptotic formula (2.2), which we denote as $A_{pq}$, can be evaluated with vector fields $g_p$, $\nabla \times E_p \times \nu$ and $W_q$ as follows:

$$A_{pq} = \int_{\partial\Omega} \nabla \times E_p \times \nu \cdot W_q ds - \int_{\partial\Omega} \nabla \times W_q \times \nu \cdot (\nu \times g_p) ds. \quad (2.3)$$

With such particular choices of vector fields $g_p$ and $W_q$, the right-hand side of (2.2) becomes

$$\rho^3 k^2 \sum_{j=1}^{j_0} \left( \frac{\varepsilon_2}{\varepsilon_1} - 1 \right) \left( \mathbf{M}_j \left( \frac{\varepsilon_2}{\varepsilon_1} \right) \hat{\theta}_p^{\perp} \right) \cdot \hat{\theta}_q^{\perp} e^{ik(\hat{\theta}_p + \hat{\theta}_q) \cdot z_j} + O(\rho^4).$$

Let $A$ be the $m \times m$ matrix with each entry given by (2.3), which can be evaluated in practice with given data. If the high order terms in the asymptotic formula (2.2) are neglected, it follows that

$$A_{pq} = \rho^3 k^2 \sum_{j=1}^{j_0} \left( \frac{\varepsilon_2}{\varepsilon_1} - 1 \right) \left( \mathbf{M}_j \left( \frac{\varepsilon_2}{\varepsilon_1} \right) \hat{\theta}_p^{\perp} \right) \cdot \hat{\theta}_q^{\perp} e^{ik(\hat{\theta}_p + \hat{\theta}_q) \cdot z_j}.$$

Let $x_0 \in \mathbb{R}^3$ be some constant vector such that $x_0 \cdot \hat{\theta}_p^{\perp} \neq 0$ for each $p$. For each $x \in \Omega$, define the vector

$$v(x) = [x_0 \cdot \hat{\theta}_1^{\perp} e^{ik\hat{\theta}_1 \cdot x}, x_0 \cdot \hat{\theta}_2^{\perp} e^{ik\hat{\theta}_2 \cdot x}, \dots, x_0 \cdot \hat{\theta}_m^{\perp} e^{ik\hat{\theta}_m \cdot x}]. \quad (2.4)$$

Referring to [1], it can be shown that there exists some integer $m_0$ that for $m \ge m_0$, the following claim is true:

$$v(x) \in \text{Range}(AA^*) \text{ if and only if } x \in \{z_1, z_2, \dots, z_{j_0}\}.$$

The singular-value decomposition of $A$ allows for localization of inclusions. In fact, the significant singular-values of $A$ determine the number of detectable inclusions. Typically, there are $j_0$ detectable inclusions when there exist $3j_0$ singular-values of $A$ that are not close to 0. The corresponding left singular-vectors of $A$ denoted by $v_1, v_2, \dots, v_{3j_0}$, allow us to locate these inclusions. Let $V = [v_1, v_2, \dots, v_{3j_0}]$, then $V\overline{V}^T$ is the projection onto the signal space of $A$. Thus the projection onto the null space is
$$P = I - V\overline{V}^T,$$

where $I$ is the $m \times m$ identity matrix, with $m > 3j_0$. The locations of inclusions can now be identified with the function

$$f(x) = \frac{1}{||Pv(x)||_2}, \qquad (2.5)$$

wherein the vector $v(x)$ is given by (2.4), and $||\cdot||_2$ is the usual $L^2$ norm applied here to a vector of $m$ components. Indeed, from the above claim, $x \in \{z_1, z_2, \dots, z_{j_0}\}$ if and only if $Pv(x) = 0$. Therefore, if we plot the function $f(x)$ for each $x \in \Omega$, then the points with high peaks represent the locations of inclusions. The obtained locations serve as an initial guess for the shape reconstruction by optimization at higher frequencies, which is discussed in detail in next section.

# 3 Reconstruction of inclusions at higher frequencies

## 3.1 The continuation approach with multiple frequency data

Suppose that multiple frequency measurements $\{ H_{k_n} \times \nu | n = 1, 2, \dots, N \}$ of the tangential trace of the magnetic field are available on the boundary for a set of wavenumbers $\{ k_n | n = 1, 2, \dots, N \}$, wherein $k_n > k_m$ if $n > m$. Assume that $z_j, j = 1, \dots, j_0$, are the locations of inclusions obtained by the MUSIC algorithm at $k = k_1$, then the shape reconstruction starts with an initial guess where each inclusion $D_j$ is a sphere centered at $z_j$ with a small radius.

To recover the shapes of inclusions at higher frequencies, we propose a continuation approach that marches from low to high wavenumbers. At each fixed wavenumber $k = k_n$ ($n > 1$), the inverse problem is recasted as a shape optimization problem. In addition, the solution at lower wavenumber $k = k_{n-1}$ is used as an initial guess, and a steepest descent method is applied to solve the optimization problem. More precisely, the continuation approach can be summarized as follows:

**FOR** $n = 2, 3, 4, \dots, N$

Let $k = k_n$,

Use the reconstruction at $k = k_{n-1}$ as an initial guess,

Solve the optimization problem at $k = k_n$ by the steepest descent method.

**END**

It is expected that with such a continuation process, the main features of the surface are recovered at low frequencies and the small-scale features are captured at higher frequencies. On the other hand, the approach leads to a convergent reconstruction without falling into some local minimum of the cost functional.

Next we focus on the reconstruction approach at each fixed wavenumber. In particular, the steepest descent direction for the cost functional and its calculation by an adjoint approach are discussed. The evolution of the inclusion surface by the chosen steepest descent direction is simulated by the level set approach, which is also briefly reviewed for the purpose of completeness.

## 3.2 Cost functional and descent directions at fixed frequency

At fixed wavenumber $k$, we recast the inverse problem of recovering the boundaries of inhomogeneities as a shape optimization problem. For a given permittivity $\epsilon_r$, define the forward
map

$$ \mathcal{M}: \varepsilon_r \to \nabla \times E \times \nu|_{\partial\Omega}, \qquad (3.1) $$

where E is the solution of (1.1) with $\varepsilon_r$. Denote the measurement of the tangential trace of the magnetic field on the boundary $\partial\Omega$ by $H^{meas}$. Here we scale the measurement by a factor of $i\omega\mu_0$. Hence $H^{meas} = \nabla \times E^t \times \nu|_{\partial\Omega}$ when no noise is present and $E^t$ is the solution of (1.1) with true permittivity. The residual on the boundary is given by

$$ \mathcal{R}(\varepsilon_r) = \mathcal{M}(\varepsilon_r) - H^{meas}, $$

and the corresponding least-squares cost functional is

$$ \mathcal{F}(\varepsilon_r) = \frac{1}{2} \|\mathcal{M}(\varepsilon_r) - H^{meas}\|_{(L^2(\partial\Omega))^3}^2. \qquad (3.2) $$

Given a vector field $V$ on $\Gamma$ such that $V(x) \in \mathbb{R}^3$ for each point $x \in \Gamma$, we define the surface $\Gamma^\delta$ as the perturbation of $\Gamma$ in the direction $V(x)$ with a step length $\delta$:

$$ \Gamma^\delta := \{x + \delta V \mid x \in \Gamma\}. $$

Let $\delta\varepsilon_r$ be the change of permittivity due to the perturbation. Then it can be shown that for any test function $u \in L^2(\Omega)$, the inner product

$$ (\delta\varepsilon_r, u) := \int_{\Omega} \delta\varepsilon_r(x)\bar{u}(x) \, dx = \int_{\text{symdiff}(D, D^\delta)} \delta\varepsilon_r(x)\bar{u}(x) \, dx. $$

Here $\bar{u}$ represents the conjugate of $u$, $D$ and $D^\delta$ are the corresponding inclusions sets bounded by the surfaces $\Gamma$ and $\Gamma^\delta$ respectively, and the symmetric difference of two sets $D$ and $D^\delta$ is given by

$$ \text{symdiff}(D, D^\delta) = (D \cup D^\delta)\setminus(D \cap D^\delta). $$

In view of the fact that the permittivity values for inhomogeneities and background medium are $\varepsilon_1$ and $\varepsilon_2$ respectively, then for an infinitesimal $\delta$, the above inner product can be simplified as

$$ (\delta\varepsilon_r, u) = \delta \int_{\Gamma} (\varepsilon_1 - \varepsilon_2)(V \cdot \nu)\bar{u}(x) \, dx, \qquad (3.3) $$

where $\nu$ is the outward unit normal on $\Gamma$.

Let $\varepsilon_r^\delta (=\varepsilon_r + \delta\varepsilon_r)$ be the permittivity associated with the perturbed interface $\Gamma^\delta$. The cost functional $\mathcal{F}(\varepsilon_r)$ admits the following expansion:

$$ \mathcal{F}(\varepsilon_r + \delta\varepsilon_r) = \mathcal{F}(\varepsilon_r) + \Re \left( (\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{meas}), \delta\varepsilon_r \right) + O(||\delta\varepsilon_r||^2), \qquad (3.4) $$

where $(\mathcal{M}')^*$ is the adjoint of the Fréchet derivative of the forward mapping $\mathcal{M}$, and $\Re$ denotes the real part of a complex number. By substituting (3.3) into (3.4) and neglecting the high order terms, it follows that

$$ \mathcal{F}(\varepsilon_r + \delta\varepsilon_r) \approx \mathcal{F}(\varepsilon_r) + \delta \Re \int_{\Gamma} (\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{meas}) \overline{(\varepsilon_1 - \varepsilon_2)} (V \cdot \nu) \, ds. $$

Therefore,

$$ \mathcal{F}'(\varepsilon_r) = \lim_{\delta \to 0} \frac{\mathcal{F}(\varepsilon_r + \delta\varepsilon_r) - \mathcal{F}(\varepsilon_r)}{\delta} = \Re \int_{\Gamma} (\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{meas}) \overline{(\varepsilon_1 - \varepsilon_2)} (V \cdot \nu) \, ds. $$
We observe that the vector field $V$ such that

$$V \cdot \nu = -\Re [(\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{\text{meas}})(\overline{\varepsilon_1 - \varepsilon_2})]$$

is a steepest descent direction, and the cost functional $\mathcal{F}(\varepsilon_r^\delta) < \mathcal{F}(\varepsilon_r)$ if $\delta \in (0, t_0]$ for some small positive constant $t_0$. Note that the tangential component of $V(x)$ does not contribute to the deformation of the shape $\Gamma$. We choose the vector field directed to the normal direction on the interface $\Gamma$ in the following way:

$$V = -\Re [(\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{\text{meas}})(\overline{\varepsilon_1 - \varepsilon_2})] \nu. \quad (3.5)$$

To calculate the adjoint of the Fréchet derivative $(\mathcal{M}')^*$, we employ the adjoint state approach, which is stated in the following Theorem.

**Theorem 3.1** Let $E$ be the solution of (1.1)-(1.2), and $\Phi$ be the solution of the following boundary value problem:

$$\left\{ \begin{array}{ll} \nabla \times \nabla \times \Phi - k^2 \bar{\varepsilon}_r \Phi = 0 & \text{in } \Omega, \\ \nu \times (\Phi \times \nu) = \mathcal{M}(\varepsilon_r) - H^{\text{meas}} & \text{on } \partial\Omega. \end{array} \right. $$

$$\text{Then } (\mathcal{M}')^*(\mathcal{M}(\varepsilon_r) - H^{\text{meas}}) = -k^2 \bar{E} \cdot \Phi.$$

By Theorem 3.1, the steepest descent direction $V(x)$ defined by (3.5) can be calculated via the following formula:

$$V = k^2 \Re [(\bar{E} \cdot \Phi)(\overline{\varepsilon_1 - \varepsilon_2})] \nu, \quad (3.6)$$

where $E$ and $\Phi$ are the solutions of the forward and adjoint problems respectively.

**Proof of Theorem 3.1** Let $\delta\varepsilon_r$ be a small perturbation of permittivity. Then the linear term of the perturbation in the magnetic field is $\mathcal{M}'\delta\varepsilon = \nabla \times \delta E \times \nu$, where the perturbation $\delta E$ satisfy the following boundary value problem:

$$\left\{ \begin{array}{ll} \nabla \times \nabla \times \delta E - k^2 \varepsilon_r \delta E = k^2 \delta \varepsilon E & \text{in } \Omega, \\ \delta E \times \nu = 0 & \text{on } \partial \Omega. \end{array} \right. \quad (3.7)$$

By using the relation $\nabla \cdot ((\nabla \times \delta E) \times \bar{\Phi}) = \nabla \times (\nabla \times \delta E) \cdot \bar{\Phi} - (\nabla \times \delta E) \cdot (\nabla \times \bar{\Phi})$, and the Green's formula,

$$\begin{aligned} \langle \mathcal{M}'\delta\varepsilon, \Phi \rangle_{\partial\Omega} &= \int_{\partial\Omega} (\nabla \times \delta E) \times \nu \cdot \bar{\Phi} \, ds \\ &= \int_{\Omega} -\nabla \times (\nabla \times \delta E) \cdot \bar{\Phi} + (\nabla \times \delta E) \cdot (\nabla \times \bar{\Phi}) \, dx \\ &= \int_{\Omega} (-k^2 \varepsilon_r \delta E - k^2 \delta \varepsilon E) \cdot \bar{\Phi} + (\nabla \times \delta E) \cdot (\nabla \times \bar{\Phi}) \, dx, \end{aligned}$$

where we have used the first equation of (3.7) in the last equality. Using the Green's formula one more time, it follows that

$$\langle \mathcal{M}'\delta\varepsilon, \Phi \rangle_{\partial\Omega} = \int_{\Omega} (-k^2 \varepsilon_r \delta E - k^2 \delta \varepsilon E) \cdot \bar{\Phi} + \delta E \cdot \nabla \times (\nabla \times \bar{\Phi}) dx + \int_{\partial\Omega} (\nabla \times \bar{\Phi}) \times \delta E \cdot \nu ds.$$

Since $\nabla \times \nabla \times \Phi - k^2 \bar{\varepsilon}_r \Phi = 0$ and $\delta E \times \nu = 0$, the above formula can be simplified as

$$\langle \mathcal{M}'\delta\varepsilon, \Phi \rangle_{\partial\Omega} = \int_{\Omega} \delta\varepsilon (-k^2) \overline{E \cdot \Phi} dx. \quad (3.8)$$

On the other hand, note that $\mathcal{M}'\delta\varepsilon$ is perpendicular to the normal direction $\nu$, hence

$$\langle \mathcal{M}'\delta\varepsilon, \mathcal{M}(\varepsilon_r) - H^{\text{meas}} \rangle_{\partial\Omega} = \langle \mathcal{M}'\delta\varepsilon, \Phi \rangle_{\partial\Omega}. \quad (3.9)$$

The proof is complete by combining (3.8) and (3.9).
### 3.3 Regularization strategy

In the presence of noisy measurements, suitable regularization strategy has to be employed to achieve a stable solution. We regularize by penalizing the area of the surface, i.e., by adding a term $\alpha \int_{\Gamma} 1 ds$ to the least-squares cost functional (3.2), where $\alpha > 0$ is the regularization parameter. Such a regularization method has been used widely in image processing to shorten the total curve length of the interfaces, see for example [25]. It is closely related to the total variation (TV) regularization method [15]. We also refer to [21] for an overview of other regularization techniques for inverse problems.

The regularization term $\alpha \int_{\Gamma} 1 ds$ used to penalize the area of the surface leads to an additional term in the steepest descent direction for the cost functional. Indeed, by a direct calculation as in [29, 20], it can be shown that

$$ \lim_{\delta \to 0} \left( \frac{1}{\delta} \int_{\Gamma_{\delta}} 1 ds - \int_{\Gamma} 1 ds \right) = \int_{\Gamma} \kappa (V \cdot \nu) ds, \qquad (3.10) $$

where $\kappa$ is the mean curvature of $\Gamma$ given by $\kappa = \nabla \cdot \nu$. Therefore, combining (3.6) and (3.10), the steepest descent direction for the cost functional with the regularization term becomes

$$ V = k^2 \Re [ (\bar{E} \cdot \Phi)(\overline{\varepsilon_1 - \varepsilon_2}) ] \nu - \alpha \kappa \nu. \qquad (3.11) $$

### 3.4 Shape evolution by level set functions

With the steepest descent direction for the cost functional $\mathcal{F}(\varepsilon_r)$ specified by (3.11), we employ the level set approach to simulate the dynamic evolution of the interface $\Gamma$. The level set method has been successfully applied to inverse problems involving shape reconstructions, see [13, 14, 15, 20, 28] and references therein for detailed discussions. For completeness, we briefly review the method and its application for the dynamic evolution of the interface $\Gamma$ with chosen velocity $V$.

Introduce an artificial time variable $t$. Let $\Gamma_t$ be the evolution of the inclusions boundary with velocity $V$, and $D_t$ be the domain bounded by $\Gamma_t$. Define the level set function $\varphi : \mathbb{R}^3 \times (0, \infty) \to \mathbb{R}$ such that

$$ \varphi(x, t) = \begin{cases} < 0, & x \in D_t \\ > 0, & x \in \Omega \setminus D_t. \end{cases} $$

The zero level set of $\varphi(x, t)$ represents the interface $\Gamma_t$ at time $t$. The topological change of the interface can be handled automatically with the level set function $\varphi$, and its dynamic evolution is captured by solving the following Hamilton-Jacobi equation:

$$ \frac{\partial \varphi(x, t)}{\partial t} + V(x, t) \cdot \nabla \varphi(x, t) = 0. \qquad (3.12) $$

In the above equation, $V(x, t)$ is the velocity of the interface evolution given by (3.11). The level set equation may be reformulated as

$$ \frac{\partial \varphi(x, t)}{\partial t} + v(x, t) |\nabla \varphi(x, t)| = \alpha \kappa(x, t) |\nabla \varphi(x, t)|, \qquad (3.13) $$

if we notice that the outward unit normal $\nu = \frac{\nabla \varphi}{|\nabla \varphi|}$. Here, the scalar function

$$ v(x, t) = k^2 \Re [ (\bar{E} \cdot \Phi)(\overline{\varepsilon_1 - \varepsilon_2}) ] . \qquad (3.14) $$
Typically, in order to solve (3.13), the velocity $v(x, t)$ needs to be extended to the entire domain. This can be accomplished by evaluating (3.14) in the whole domain $\Omega$.

To discretize the level set equation, the Euler method may be applied in the time direction. This gives rise to

$$ \frac{\varphi^{n+1}(x) - \varphi^n(x)}{\Delta t} + v^n(x)|\nabla\varphi^n(x)| = \alpha\kappa^n(x)|\nabla\varphi^n(x)|, $$

where $\Delta t$ is the time step length. Whereas in the space direction, the first-order accurate upwind scheme, or ENO, WENO schemes with higher-order spatial accuracy may be applied for the convection term $v^n(x)|\nabla\varphi^n(x)|$, and the center difference scheme can be used for the diffusion term $\alpha\kappa^n|\nabla\varphi^n(x)|$. We refer the reader to [27] that addresses this important topic in detail.

Note that only the zero level set of $\varphi$, which determines the interface $\Gamma$, is critical, while the value of $\varphi$ itself is not that important. In general, the interface is more accurately detected if $\varphi$ is a signed distance function. That is

$$ \varphi(x, t) = d(x, \overline{D_t}) - d(x, \overline{\Omega \setminus D_t}). $$

The distance between $x$ and a closed set $D$ is defined as

$$ d(x, D) = \min_{y \in D} |x - y|. $$

The signed distance function $\varphi$ satisfies the equation

$$ |\nabla \varphi| = 1. $$

Typically, after the evolution of the interface by solving the Hamilton-Jacobi for a few time steps, $\varphi$ is not a signed distance function any more. Usually a reinitialization approach for the level set function is applied periodically such that $\varphi$ becomes a signed distance function. For more details about the construction of the signed distance function, we refer to the Chapter 7 of [27].

# 4 Numerical examples

Several numerical examples are presented in this section to illustrate the efficiency and robustness of the proposed method. In all examples, we assume that the relative permittivity values of inhomogeneities and background medium are $\varepsilon_1 = 2.5$ and $\varepsilon_2 = 2.0$, respectively. For simplicity, the domain $\Omega$ is a unit ball centered at origin in the examples. However, we note that the proposed algorithm is suitable for general domains with piecewise smooth boundaries. To simulate the forward model (1.1)-(1.2) and to obtain the solution for the adjoint problem defined in Theorem 3.1, we employ the second-order edge element method [23, 24, 26].

In all examples, it is assumed that a 5% uniformly distributed noise is added to the simulation data. We also set a uniform frequency spacing in the continuation algorithm such that $k_{n+1} - k_n = 0.3$ for $n = 2, 3, \dots, N-1$. In addition, the relative residual on the boundary for the obtained reconstruction is defined by

$$ r = \frac{\|\mathcal{M}(\varepsilon_r) - H^{\text{meas}}\|_{(L^2(\partial\Omega))^3}}{\|H^{\text{meas}}\|_{(L^2(\partial\Omega))^3}}, $$

where $\mathcal{M}$ is the forward map given in (3.1).
**Example 1** We consider a bean-type inclusion D shown in Figure 1 (left). The surface of the inclusion is described by the equation

$$ \frac{\bar{x}_1^2}{R^2} + \frac{(\bar{x}_2 - \alpha_1 R \cos(\pi \bar{x}_1 / R))^2}{R^2(1 - \alpha_2 \cos(\pi \bar{x}_1 / R))} + \frac{\bar{x}_3^2}{R^2(1 - \alpha_3 \cos(\pi \bar{x}_1 / R))} = 1, \quad (4.1) $$

with $\bar{x}_1 = x_1 - c_1$, $\bar{x}_2 = x_2 - c_2$, $\bar{x}_3 = x_3 - c_3$.

In fact, in the above equation, $(c_1, c_2, c_3)$ represents the location of D, and $2R$ denotes the largest dimension of the object. Here we set $c_1 = 0.2$, $c_2 = 0.1$, $c_3 = 0.05$, and $R = 0.4$. The parameters $\alpha_1, \alpha_2$ and $\alpha_3$ are chosen such that $\alpha_1 = 0.15$, $\alpha_2 = 0.3$, and $\alpha_3 = 0.2$.

The location of the inclusion is obtained by the MUSIC algorithm at lowest frequency $k_1 = 1.0$. It is obvious that with chosen frequency, the non-dimensionalized parameter $k\rho < 1$ in (2.2), thus the asymptotic expansion is valid. The identification function $1/||Pv(x)||_2$ is evaluated in the whole domain $\Omega$. Figure 1 (right) shows the surface plot of $1/||Pv(x)||_2 = 650$, which represents the location of the inclusion D. The function values $1/||Pv(x)||_2$ on the cross-sectional planes $x_1 = c_1$, $x_2 = c_2$ and $x_3 = c_3$ are also displayed in Figure 2. It is observed that the identification function obtains high peaks near $(c_1, c_2, c_3)$.

To reconstruct the shape of the inclusion D, the continuation approach is employed, where the wavenumber starts at $k_2 = 1.2$ and stops at $k_N = 5.7$. As mentioned at the beginning of this section, a uniform frequency spacing is used such that $k_{n+1} - k_n = 0.3$ for $n \ge 2$. Figure 3 plots the relative residual during the iteration process when the wavenumber $k = 2.4, 3.0$, and $4.5$, respectively. In order to demonstrate the convergence the continuation approach at higher frequencies, the evolution of the reconstruction process at different wavenumbers is displayed in Figure 4. It is seen that the inclusion is successfully reconstructed at the largest wavenumber $k = 5.7$. To be more precise, we plot the curves of true inclusion and numerical reconstruction on the cross-sectional planes $x_3 = -0.2, 0.05$ and $0.35$ in Figure 5. The boundary of D is accurately obtained even though the noise is present.

**Example 2** Assume that there are two bean-type inclusions $D_1$ and $D_2$ as shown in Figure 6 (left), where the surface of each is described by the equation (4.1). Let $c_1 = c_3 = 0$, $\tilde{c}_2 = -0.4$, and $\hat{c}_2 = 0.4$. The coordinates $(c_1, \tilde{c}_2, c_3)$ and $(c_1, \hat{c}_2, c_3)$ represent the locations of $D_1$ and $D_2$, i.e., in the equation (4.1), $\bar{x}_1 = x_1 - c_1$, $\bar{x}_2 = x_2 - \tilde{c}_2$, $\bar{x}_3 = x_3 - c_3$ and $\bar{x}_1 = x_1 - c_1$, $\bar{x}_2 = x_2 - \hat{c}_2$, $\bar{x}_3 = x_3 - c_3$ for $D_1$ and $D_2$ respectively.

The locations of $D_1$ and $D_2$ are obtained by the MUSIC algorithm at lowest frequency $k_1 = 1.0$. We display the surface plot of $1/||Pv(x)||_2 = 250$ in Figure 6 (right), where it is seen that two inclusions are successfully identified. Such a resolution beyond the Rayleigh limit has also been reported in [2]. For clearness, the function values $1/||Pv(x)||_2$ on the cross-sectional planes $x_1 = c_1$, $x_2 = \tilde{c}_2$, $x_2 = \hat{c}_2$ and $x_3 = c_3$ are displayed in Figure 7. It is observed that the function values display high peaks near $(c_1, \tilde{c}_2, c_3)$ and $(c_1, \hat{c}_2, c_3)$.

Figure 8 plots the relative residual during the iteration process at specified wavenumbers, and Figure 9 demonstrates the evolution of reconstruction process. To investigate the accuracy of the proposed approach, the curves of true inclusions and numerical reconstructions on the cross-sectional planes $x_3 = -0.15, 0$ and $0.15$ are plotted for comparison (see Figure 10).

**Example 3** In this example, we consider the case where there are three inclusions inside $\Omega$ as shown in Figure 11 (left). The inhomogeneities consist of two ellipsoids $D_1$ and $D_2$ centered at $(c_{1,1}, c_{1,2}, c_{1,3})$ and $(c_{2,1}, c_{2,2}, c_{2,3})$ respectively, and a ball $D_3$ centered at $(c_{3,1}, c_{3,2}, c_{3,3})$, where $(c_{1,1}, c_{1,2}, c_{1,3}) = (0.3, -0.35, 0.3)$, $(c_{2,1}, c_{2,2}, c_{2,3}) = (-0.35, -0.4, -0.2)$,
and $(c_{3,1}, c_{3,2}, c_{3,3}) = (-0.15, 0.35, 0.1)$. In addition, $D_1$ and $D_2$ have the same size, where the semi-axes lengths are 0.32, 0.22, and 0.22. The radius of the inclusion $D_3$ is 0.25.

At lowest frequency $k_1 = 1.0$, the identification function $1/||Pv(x)||_2$ is evaluated in the whole domain $\Omega$ to obtain the locations of $D_1, D_2$ and $D_3$. Figure 11 (right) is surface plot of $1/||Pv(x)||_2 = 400$. We display the function values $1/||Pv(x)||_2$ on the cross-sectional planes $x_1 = c_{i,1}, x_2 = c_{i,2}$ and $x_3 = c_{i,3}$ ($i = 1, 2, 3$) in Figure 12. From the figure, it is seen that the identification function obtains high peaks near the centers of inclusions.

The continuation approach is applied to reconstruct these inclusions, where the relative residual on the boundary during the iteration process for $k = 1.5, 2.1$ and $3.6$ is plotted in Figure 13. The shapes of numerical reconstructions at different wavenumbers are shown in Figure 14, which confirms the convergence of the algorithm. The boundaries of $D_1, D_2$ and $D_3$ are accurately obtained at the highest wavenumber $k = 5.7$, as can be observed from the plots on the cross-sectional planes $x_3 = -0.2, 0.1$ and $0.3$ (see Figure 15).

# 5 Discussions

We have studied a two-stage reconstruction method to solve the inverse problem for the three-dimensional Maxwell's equations. The locations of inhomogeneities are obtained by the MUSIC algorithm at low frequency. Based on knowledge of their locations, the shapes are recovered by a continuation approach with multiple frequency data. The proposed method is accurate and stable, as demonstrated by the numerical examples. One of the challenges along this research direction is to develop fast solvers for the Maxwell's equation at high frequencies. Indeed, the linear system resulting from the discretization by edge element methods becomes more ill-conditioned as frequency increases. More efficient preconditioners need to be developed to accelerate the linear solver. The fast boundary integral equation method may also be designed to solve the Maxwell's equations in a more efficient way.

<figure><img src="image_2.png" /><figcaption>Figure 1: True inclusion (left) and its localization by the MUSIC algorithm (right).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Values of $1/||Pv(x)||_2$ on the cross-sectional planes. (a): $x_1 = c_1$; (b): $x_2 = c_2$; (c): $x_3 = c_3$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: The relative residual during the iteration process for $k = 2.4$, 3.0 and 4.5, respectively.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the reconstruction at $k = 1.2, 1.8, 3.0, 4.2, 4.8, 5.7$ (left to right, top to bottom).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Reconstruction on the cross-sectional planes $x_3 = -0.2$ (left), $x_3 = 0.05$ (middle) and $x_3 = 0.35$ (right). The solid line represents the boundary of true inclusion, and the dotted line is the boundary of reconstruction.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: True inclusions (left) and localization by the MUSIC algorithm (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Values of $1/||Pv(x)||_2$ on the cross-sectional planes. (a): $x_1 = c_1$; (b): $x_2 = \tilde{c}_2$; (c): $x_2 = \hat{c}_2$; (d): $x_3 = c_3$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: The relative residual during the iteration process for $k = 1.5$, $2.1$ and $4.2$, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution of the reconstruction at $k = 1.2, 1.8, 3.0, 4.2, 4.8, 5.7$ (left to right, top to bottom).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Reconstruction on the cross-sectional planes $x_3 = -0.15$ (left), $x_3 = 0$ (middle) and $x_3 = 0.15$ (right). The solid line represents the boundary of true inclusion, and the dotted line is the boundary of reconstruction.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: True inclusions (left) and localization by the MUSIC algorithm (right).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Values of $1/||Pv(x)||_2$ on the cross-sectional planes. (a): $x_1 = c_{1,1}$; (b): $x_2 = c_{1,2}$; (c): $x_3 = c_{1,3}$; (d): $x_1 = c_{2,1}$; (e): $x_2 = c_{2,2}$; (f): $x_3 = c_{2,3}$; (g): $x_1 = c_{3,1}$; (h): $x_2 = c_{3,2}$; (i): $x_3 = c_{3,3}$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: The relative residual during the iteration process for $k = 1.5, 2.1$ and $3.6$, respectively.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Evolution of the reconstruction at $k = 1.2, 2.4, 4.2, 5.7$ (left to right, top to bottom).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Reconstruction on the cross-sectional planes $x_3 = -0.2$ (left), $x_3 = 0.1$ (middle) and $x_3 = 0.3$ (right). The solid line represents the boundary of true inclusion, and the dotted line is the boundary of reconstruction.</figcaption></figure>
# Acknowledgements

G. Bao's research was supported in part by the NSF grants DMS-0908325, DMS-0968360, DMS-1211292, the ONR grant N00014-12-1-0319, a Key Project of the Major Research Plan of NSFC (No. 91130004), and a special research grant from Zhejiang University. S.M. Mefire's research was partly supported by the IMMELEXT Project of the University of Henri Poincaré — Nancy-I and of the Regional Council of Lorraine in France.

# References

[1] H. Ammari, E. Iakovleva, and D. Lesselier, *A MUSIC algorithm for locating small inclusions buried in a half-space from the scattering amplitude at a fixed frequency*, Multiscale Model. Sim., **3** (2005), 597-628.

[2] H. Ammari, E. Iakovleva, D. Lesselier, and G. Perrusson, *MUSIC-type electromagnetic imaging of a collection of small three-dimensional inclusions*, SIAM J. Sci. Comput., **29** (2007), 674-709.

[3] H. Ammari and H. Kang, *Reconstruction of Small Inhomogeneities from Boundary Measurements*, Lecture Notes in Mathematics, vol. 1846, Springer-verlag, Berlin, 2004.

[4] H. Ammari, S. Moskow, and M. Vogelius, *Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter*, ESAIM Control Optim. Calc. Var., **9** (2003), 49-66.

[5] H. Ammari, M. Vogelius, and D. Volkov, *Asymptoic formulas for perturbations in electromagnetic fields due to the presence of imperfections of small diameter II. the full Maxwell equations*, J. Math. Pures. Appl., **80** (2001), 769-814.

[6] M. Asch and S. M. Mefire, *Numerical localization of electromagnetic imperfections from a perturbation formula in three dimensions*, J. Comput. Math., **26** (2008), 149-195.

[7] G. Bao, S. Hou, and P. Li, *Inverse scattering by a continuation method with initial guesses from a direct imaging algorithm*, J. Comput. Phy., **227** (2007), 755-762.

[8] G. Bao and P. Li, *Inverse medium scattering for three-dimensional time harmonic Maxwell's equations*, Inverse Problems, **20** (2004), L1-L7.

[9] G. Bao and P. Li, *Inverse medium scattering problems for electromagnetic waves*, SIAM J. Appl. Math., **65** (2005), 2049-2066.

[10] G. Bao and P. Li, *Numerical solution of an inverse medium scattering problem for Maxwell's equations at fixed frequency*, J. Comput. Phys., **228** (2009), 4638-4648.

[11] G. Bao and J. Lin, *Imaging of local surface displacement on an infinite ground plane: the multiple frequency case*, SIAM J. Appl. Math., **71** (2011), 1733-1752.

[12] G. Bao and F. Triki, *Error estimates for the recursive linearization for solving inverse medium problems*, J. Comput. Math., **28** (2010), 725-744.

[13] M. Burger, *A level set method for inverse problems*, Inverse Problems, **17** (2001), 1327-155.
[14] M. Burger and S. Osher, *A survey on level set methods for inverse problems and optimal design*, Eur. J. Appl. Math., **16** (2005), 263-301.

[15] T. Chan and X. Tai, *Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients*, J. Comput. Phys., **193** (2003), 40-66.

[16] Y. Chen, *Inverse scattering via Heisenberg's uncertainty principle*, Inverse Problems, **13** (1997), 1-13.

[17] F. Cakoni and D. Colton, *Qualitative Methods in Inverse Scattering Theory*, Springer-Verlag, Berlin, 2006.

[18] D. Colton and A. Kirsch, *A simple method for solving inverse scattering problems in the resonance region*, Inverse Problems, **12** (1996), 383-393.

[19] D. Colton and R. Kress, *Inverse Acoustic and Electromagnetic Scattering Theory*, Applied Mathematical Sciences, vol. 93, Springer-Verlag, Berlin, 1998.

[20] O. Dorn and D. Lesselier, *Level set method for inverse scattering*, Inverse Problems, **22** (2006), 67-131.

[21] H. W. Engl, M. Hanke, and A. Neubauer, *Regularization of Inverse Problems*, Mathematics and Its Application, Kluwer Academic Publishers, 1996.

[22] A. Kirsch, *Characterization of the shape of a scattering obstacle using the spectral data of the far-field operator*, Inverse Problems, **14** (1998), 1489-1512.

[23] J. F. Lee, D. K. Sun, and Z. J. Cendes, *Tangential vector finite element for electromagnetic field computation*, IEEE. Trans. on Magnetics, **27** (1991), 4032-4035.

[24] P. Monk, *Finite Element Methods for Maxwell's Equations*, Oxford University Press, Oxford, U.K., 2003.

[25] D. Mumford and J. Shah, *Optimal approximation by piecewise smooth functions and associated variational problems*, Commun. Pure Appl. Math., **42** (1989), 577-685.

[26] J.-C. Nédélec, *Mixed finite elements in $\mathbb{R}^3$*, Numer. Math., **35** (1980), 315-341.

[27] S. Osher and R. Fedkiw, *Level Set Methods and Dynamic Implicit Surfaces*, Applied Mathematical Sciences, vol. 153, Springer, New York, 2003.

[28] F. Santosa, *A level set approach for inverse problems involving obstacles*, ESAIM Control Optim. Calculus Variations, **1** (1996), 17-33.

[29] J. Sokolowski and J. Zolesio, *Introduction to Shape Optimization: Shape Sensitivity Analysis*, Springer-Verlag, Berlin, 1992.