Gang Bao
email: bao@math.msu.edu

Junshan Lin

Séraphin M Mefire
email: seraphin.mefire@univ-lorraine.fr.

Numerical reconstruction of electromagnetic inclusions in three dimensions

Keywords: inverse problems, Maxwell equations, shape optimization MS Codes. 65N21, 65N30, 78A25

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.

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 [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF][START_REF] Cakoni | Qualitative Methods in Inverse Scattering Theory[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. 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 ε and µ be the electric permittivity and magnetic permeability of the medium respectively in a bounded domain Ω ⊂ R 3 . Denote the permittivity and permeability of the vacuum as ε 0 and µ 0 . Then ε and µ can be rewritten as ε = ε r ε 0 and µ = µ r µ 0 , where ε r and µ r are the relative permittivity and permeability respectively. In this work, we are interested in the case when the permittivity ε r is piecewise constant inside the domain Ω. To be more precise, let D j (j = 1, 2, • • • , j 0 ) be the disconnected inclusions inside the domain Ω such that D i ∩ D j = ∅ if i = j, and D = ∪ j 0 j=1 D j be the union of inclusions. Then the relative permittivity takes the following form:

ε r (x) = ε 1 if x ∈ D, ε 2 if x ∈ Ω\D,
where ε 1 and ε 2 are permittivities for inhomogeneities and background medium respectively. Throughout the paper, we also restrict out attention to the nonmagnetic medium by assuming that µ r ≡ 1 in Ω. For a time harmonic electromagnetic wave (with e -iωt dependence), the electric field E and magnetic field H satisfy the Maxwell's equations in the bounded domain Ω:

∇ × E = iωµ 0 H, ∇ × H = -iωε r ε 0 E,
where ω 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 ∇ × ∇ × E -k 2 ε r E = 0 in Ω. ( 1.1)
Here k = ω √ ε 0 µ 0 denotes the wavenumber. Let ν be the outward unit normal on the boundary ∂Ω 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 × ν = g on ∂Ω. (1.2)
From boundary measurements of the magnetic field H × ν, or equivalently ∇ × E × ν, the inverse problem is to reconstruct the inclusions D j , j = 1, ..., j 0 , or their boundaries ∂D j . For clarity of exposition, let us denote the boundaries as Γ := ∪ j 0 j=1 ∂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 [START_REF] Bao | Inverse scattering by a continuation method with initial guesses from a direct imaging algorithm[END_REF] for a related inverse medium scattering problem in two dimensions, and [START_REF] Bao | Inverse medium scattering for three-dimensional time harmonic Maxwell's equations[END_REF][START_REF] Bao | Inverse medium scattering problems for electromagnetic waves[END_REF][START_REF] Bao | Numerical solution of an inverse medium scattering problem for Maxwell's equations at fixed frequency[END_REF] 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 (MU-SIC) algorithm at low frequency to obtain the location of each inclusion D j (j = 1, 2, • • • , 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 [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF], and is further developed in [START_REF] Ammari | Asympotoic formulas for perturbations in electromagnetic fields due to the presence of imperfections of small diamerter II. the full Maxwell equations[END_REF] for the three-dimensional Maxwell's equations. From the perturbation formula, a MUSIC algorithm is proposed and studied in details in [START_REF] Asch | Numerical localization of electromagnetic imperfections from a perturbation formula in three dimensions[END_REF] to locate the electromagnetic inclusions.

In the second stage, in order to obtain the boundaries of inclusions, the nonlinear inverse problem is then recasted 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 [START_REF] Bao | Inverse medium scattering for three-dimensional time harmonic Maxwell's equations[END_REF][START_REF] Bao | Inverse medium scattering problems for electromagnetic waves[END_REF][START_REF] Bao | Imaging of local surface displacement on an infinite ground plane: the multiple frequency case[END_REF][START_REF] Bao | Error estimates for the recursive linearization for solving inverse medium problems[END_REF][START_REF] Chen | Inverse scattering via Heisenberg's uncertainty principle[END_REF]. 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 Γ, 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.

A MUSIC algorithm at low frequency

Assume that the typical size of each inclusion D j (j = 1, 2, • • • , j 0 ) is ρ. Moreover, D j is represented by z j +ρ 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

d 0 ≤ |z j -z l | ∀ j = l, d 0 ≤ dist(z j , ∂Ω) ∀ j ,
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

∇ × ∇ × E 0 -k 2 E 0 = 0 in Ω, E 0 × ν = g on ∂Ω, (2.1)
has a unique solution E 0 . Let W be any smooth solution of the Maxwell equation

∇ × ∇ × W -k 2 W = 0 in Ω,
then there exists a constant ρ 0 that for 0 < ρ < ρ 0 , the following asymptotic formula holds (see [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF]):

∂Ω ∇ × E × ν • W ds - ∂Ω ∇ × W × ν • (ν × (E × ν))ds = ρ 3 k 2 j 0 j=1 ε 2 ε 1 -1 M j ( ε 2 ε 1 )E 0 (z j ) • W (z j ) + O(ρ 4 ), (2.2)
with E the unique weak solution of the boundary value problem (1.1) -(1.2).

Here the constant ρ 0 depends on the shapes U j , Ω, permittivity values ε 1 and ε 2 , the wavenumber k and d 0 , but is otherwise independent of g, W , and of the points z j , j = 1, ..., j 0 . Also, M j is the polarization tensor (see [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF] for detailed properties), where each entry is given by

m j pq ε 2 ε 1 = ε 1 ε 2 U j ∂φ p ∂x q dx, 1 ≤ p, q ≤ 3,
with the scalar potential φ p satisfying

             ∆φ p = 0 in R 3 \U j , ∆φ p = 0 in U j , φ + p -φ - p = 0 on ∂U j , ε 2 ε 1 ∂φp ∂ν + -∂φp ∂ν - = 0 on ∂U j , φ p (x) -x p → 0 as |x| → ∞,
where the outward unit normal to ∂(z j + ρU j ), the boundary of z j + ρU j , is also denoted by ν, and the superscripts +,indicate the limiting values as ∂(z j + ρU j ) is approached from outside z j + ρU j , and from inside z j + ρU j . Also, the coordinates of x ∈ R 3 are represented by x q , 1 ≤ q ≤ 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 [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Kirsch | Characterization of the shape of a scattering obstacle using the spectral data of the far-field operator[END_REF]. For simplicity of presentation, we assume that Ω is a unit sphere. For p = 1, 2,

• • • , m, let E 0,p (x) = θ⊥ p e ik θp•x ,
where θp is the incidence direction, and θ⊥ p 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) = θ⊥ p e ik θp•x × ν on the boundary ∂Ω.
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) = θ⊥ q e ik θq•x , q = 1, 2, • • • , 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 , ∇ × E p × ν and W q as follows:

A pq = ∂Ω ∇ × E p × ν • W q ds - ∂Ω ∇ × W q × ν • (ν × g p )ds.
(

With such particular choices of vector fields g p and W q , the right-hand side of (2.2) becomes

ρ 3 k 2 j 0 j=1 ε 2 ε 1 -1 M j ε 2 ε 1 θ⊥ p • θ⊥ q e ik( θp+ θq)•zj + O(ρ 4 ).
Let A be the m × 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 = ρ 3 k 2 j 0 j=1 ε 2 ε 1 -1 M j ε 2 ε 1 θ⊥ p • θ⊥ q e ik( θp+ θq)•zj .
Let x 0 ∈ R 3 be some constant vector such that x 0 • θ⊥ p = 0 for each p. For each x ∈ Ω, define the vector

v(x) = [x 0 • θ⊥ 1 e ik θ1 •x , x 0 • θ⊥ 2 e ik θ2 •x , • • • , x 0 • θ⊥ m e ik θm•x ]. (2.4)
Referring to [START_REF] Ammari | A MUSIC algorithm for locating small inclusions buried in a half-space from the scattering amplitude at a fixed frequency[END_REF], it can be shown that there exists some integer m 0 that for m ≥ m 0 , the following claim is true:

v(x) ∈ Range(AA * ) if and only if x ∈ {z 1 , z 2 , • • • , 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 , • • • , v 3j 0 , allow us to locate these inclusions. Let V = [v 1 , v 2 , • • • , v 3j 0 ], then V V
T is the projection onto the signal space of A. Thus the projection onto the null space is

P = I -V V T ,
where I is the m × m identity matrix, with m > 3j 0 . The locations of inclusions can now be identified with the function

f (x) = 1 ||P v(x)|| 2 , (2.5)
wherein the vector v(x) is given by (2.4), and . 2 is the usual L 2 norm applied here to a vector of m components. Indeed, from the above claim, x ∈ {z 1 , z 2 , • • • , z j 0 } if and only if P v(x) = 0. Therefore, if we plot the function f (x) for each x ∈ Ω, 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

| n = 1, 2, • • • , N }, wherein k n > k m if n > m.
Assume that z j , j = 1, • • • , 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, • • • , 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.

Cost functional and descent directions at fixed frequency

At fixed wavenumer k, we recast the inverse problem of recovering the boundaries of inhomogeneities as a shape optimization problem. For a given permittivity ε r , define the forward map

M : ε r → ∇ × E × ν| ∂Ω , (3.1)
where E is the solution of (1.1) with ε r . Denote the measurement of the tangential trace of the magnetic field on the boundary ∂Ω by H meas . Here we scale the measurement by a factor of iωµ 0 . Hence H meas = ∇ × E t × ν| ∂Ω 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

R(ε r ) = M(ε r ) -H meas ,
and the corresponding least-squares cost functional is

F (ε r ) = 1 2 ||M(ε r ) -H meas || 2 (L 2 (∂Ω)) 3 . (3.2)
Given a vector field V on Γ such that V (x) ∈ R 3 for each point x ∈ Γ, we define the surface Γ δ as the perturbation of Γ in the direction V (x) with a step length δ:

Γ δ := {x + δ V | x ∈ Γ}.
Let δε r be the change of permittivity due to the perturbation. Then it can be shown that for any test function u ∈ L 2 (Ω), the inner product

(δε r , u) := Ω δε r (x)ū(x) dx = symdiff(D,D δ ) δε r (x)ū(x) dx.
Here ū represents the conjugate of u, D and D δ are the corresponding inclusions sets bounded by the surfaces Γ and Γ δ respectively, and the symmetric difference of two sets D and D δ is given by

symdiff(D, D δ ) = (D ∪ D δ )\(D ∩ D δ ).
In view of the fact that the permittivity values for inhomogeneities and background medium are ε 1 and ε 2 respectively, then for an infinitesimal δ, the above inner product can be simplified as

(δε r , u) = δ Γ (ε 1 -ε 2 )(V • ν)ū(x) dx, (3.3)
where ν is the outward unit normal on Γ. Let ε δ r (= ε r + δε r ) be the permittivity associated with the perturbed interface Γ δ . The cost functional F (ε r ) admits the following expansion:

F (ε r + δε r ) = F (ε r ) + ℜ ( (M ′ ) * (M(ε r ) -H meas ), δε r ) + O(||δε r || 2 ), (3.4)
where (M ′ ) * is the adjoint of the Fréchet derivative of the forward mapping M, and ℜ denotes the real part of a complex number. By substituting (3.3) into (3.4) and neglecting the high order terms, it follows that

F (ε r + δε r ) ≈ F (ε r ) + δ ℜ Γ (M ′ ) * (M(ε r ) -H meas )(ε 1 -ε 2 )(V • ν) ds.
Therefore,

F ′ (ε r ) = lim δ→0 F (ε r + δε r ) -F (ε r ) δ = ℜ Γ (M ′ ) * (M(ε r ) -H meas )(ε 1 -ε 2 )(V • ν) ds.
We observe that the vector field V such that

V • ν = -ℜ [ (M ′ ) * (M(ε r ) -H meas )(ε 1 -ε 2 ) ]
is a steepest descent direction, and the cost functional F (ε δ r ) < F (ε r ) if δ ∈ (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 Γ. We choose the vector field directed to the normal direction on the interface Γ in the following way:

V = -ℜ [ (M ′ ) * (M(ε r ) -H meas )(ε 1 -ε 2 ) ] ν. (3.5)
To calculate the adjoint of the Fréchet derivative (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 Φ be the solution of the following boundary value problem:

∇ × ∇ × Φ -k 2 εr Φ = 0 in Ω, ν × (Φ × ν) = M(ε r ) -H meas on ∂Ω. Then (M ′ ) * (M(ε r ) -H meas ) = -k 2 Ē • Φ.
By Theorem 3.1, the steepest descent direction V (x) defined by (3.5) can be calculated via the following formula:

V = k 2 ℜ [ ( Ē • Φ)(ε 1 -ε 2 ) ] ν, (3.6)
where E and Φ are the solutions of the forward and adjoint problems respectively.

Proof of Theorem 3.1 Let δε r be a small perturbation of permittivity. Then the linear term of the perturbation in the magnetic field is M ′ δε = ∇ ×δE ×ν, where the perturbation δE satisfy the following boundary value problem:

∇ × ∇ × δE -k 2 ε r δE = k 2 δεE in Ω, δE × ν = 0 on ∂Ω. (3.7)
By using the relation ∇

• ((∇ × δE) × Φ) = ∇ × (∇ × δE) • Φ -(∇ × δE) • (∇ × Φ)
, and the Green's formula,

< M ′ δε, Φ > ∂Ω = ∂Ω (∇ × δE) × ν • Φ ds = Ω -∇ × (∇ × δE) • Φ + (∇ × δE) • (∇ × Φ) dx = Ω (-k 2 ε r δE -k 2 δεE) • Φ + (∇ × δE) • (∇ × Φ) dx,
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

< M ′ δε, Φ > ∂Ω = Ω (-k 2 ε r δE -k 2 δεE) • Φ + δE • ∇ × (∇ × Φ)dx + ∂Ω (∇ × Φ) × δE • νds.
Since ∇ × ∇ × Φk 2 εr Φ = 0 and δE × ν = 0, the above formula can be simplified as

< M ′ δε, Φ > ∂Ω = Ω δε (-k 2 ) Ē • Φdx. (3.8)
On the other hand, note that M ′ δε is perpendicular to the normal direction ν, hence

< M ′ δε, M(ε r ) -H meas > ∂Ω =< M ′ δε, Φ > ∂Ω . (3.9)
The proof is complete by combining (3.8) and (3.9).

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 α Γ 1 ds to the least-squares cost functional (3.2), where α > 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 [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problems[END_REF]. It is closely related to the total variation (TV) regularization method [START_REF] Chan | Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients[END_REF]. We also refer to [START_REF] Engl | Regularization of Inverse Problems, Mathematics and Its Application[END_REF] for an overview of other regularization techniques for inverse problems. The regularization term α Γ 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 [START_REF] Sokolowski | Introduction to Shape Optimization: Shape Sensitivity Analysis[END_REF][START_REF] Dorn | Level set method for inverse scattering[END_REF], it can be shown that lim

δ→0 1 δ Γ δ 1 ds - Γ 1 ds = Γ κ (V • ν) ds, (3.10)
where κ is the mean curvature of Γ given by κ = ∇ • ν. Therefore, combining (3.6) and (3.10), the steepest descent direction for the cost functional with the regularization term becomes

V = k 2 ℜ [ ( Ē • Φ)(ε 1 -ε 2 ) ] ν -ακν. (3.11)

Shape evolution by level set functions

With the steepest descent direction for the cost functional F (ε r ) specified by (3.11), we employ the level set approach to simulate the dynamic evolution of the interface Γ. The level set method has been successfully applied to inverse problems involving shape reconstructions, see [START_REF] Burger | A level set method for inverse problems[END_REF][START_REF] Burger | A survey on level set methods for inverse problems and optimal design[END_REF][START_REF] Chan | Level set and total variation regularization for elliptic inverse problems with discontinuous coefficients[END_REF][START_REF] Dorn | Level set method for inverse scattering[END_REF][START_REF] Santosa | A level set approach for inverse problems involving obstacles[END_REF] and references therein for detailed discussions. For completeness, we briefly review the method and its application for the dynamic evolution of the interface Γ with chosen velocity V . Introduce an artificial time variable t. Let Γ t be the evolution of the inclusions boundary with velocity V , and D t be the domain bounded by Γ t . Define the level set function ϕ :

R 3 × (0, ∞) → R such that ϕ(x, t) = < 0, x ∈ D t > 0, x ∈ Ω\D t .
The zero level set of ϕ(x, t) represents the interface Γ t at time t. The topological change of the interface can be handled automatically with the level set function ϕ, and its dynamic evolution is captured by solving the following Hamilton-Jacobi equation:

∂ϕ(x, t) ∂t + V (x, t) • ∇ϕ(x, t) = 0. (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

∂ϕ(x, t) ∂t + v(x, t)|∇ϕ(x, t)| = ακ(x, t)|∇ϕ(x, t)|, (3.13)
if we notice that the outward unit normal ν = ∇ϕ |∇ϕ| . Here, the scalar function

v(x, t) = k 2 ℜ [ ( Ē • Φ)(ε 1 -ε 2 ) ]. (3.14)
Typically, in order to solve (3.13), the velocity v(x, t) needs to extended to the entire domain. This can be accomplished by evaluating (3.14) in the whole domain Ω. To discretize the level set equation, the Euler method may be applied in the time direction. This gives rise to

ϕ n+1 (x) -ϕ n (x) ∆t + v n (x)|∇ϕ n (x)| = ακ n (x)|∇ϕ n (x)|,
where ∆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)|∇ϕ n (x)|, and the center difference scheme can be used for the diffusion term ακ n |∇ϕ n (x)|. We refer the reader to [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] that addresses this important topic in detail.

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

ϕ(x, t) = d(x, D t ) -d(x, Ω\D t ).
The distance between x and a closed set D is defined as

d(x, D) = min y∈D |x -y|.
The signed distance function ϕ satisfies the equation

|∇ϕ| = 1.
Typically, after the evolution of the interface by solving the Hamilton-Jacobi for a few time steps, ϕ is not a signed distance function any more. Usually a reinitialization approach for the level set function is applied periodically such that ϕ becomes a signed distance function.

For more details about the construction of the signed distance function, we refer to the Chapter 7 of [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF].

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 ε 1 = 2.5 and ε 2 = 2.0, respectively. For simplicity, the domain Ω 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 [START_REF] Lee | Tangential vector finite element for electromagnetic field compuation[END_REF][START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF][START_REF] Nédélec | Mixed finite elements in R 3[END_REF].

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+1k n = 0.3 for n = 2, 3, • • • , N -1. In addition, the relative residual on the boundary for the obtained reconstruction is defined by

r = ||M(ε r ) -H meas || (L 2 (∂Ω)) 3 ||H meas || (L 2 (∂Ω)) 3 ,
where 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

x2 1 R 2 + (x 2 -α 1 R cos(πx 1 /R)) 2 R 2 (1 -α 2 cos(πx 1 /R)) + x2 3 R 2 (1 -α 3 cos(πx 1 /R)) = 1, (4.1)
with x1 = x 1 -c 1 , x2 = x 2 -c 2 , x3 = 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 α 1 , α 2 and α 3 are chosen such that α 1 = 0.15, α 2 = 0.3, and α 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ρ < 1 in (2.2), thus the asymptotic expansion is valid. The identification function 1/||P v(x)|| 2 is evaluated in the whole domain Ω. Figure 1 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+1k n = 0.3 for n ≥ 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. 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/||P v(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 [START_REF] Ammari | MUSIC-type electromagnetic imaging of a collection of small three-dimensional inclusions[END_REF]. For clearness, the function values 1/||P v(x)|| 2 on the cross-sectional planes x 1 = c 1 , x 2 = c2 , x 2 = ĉ2 and x 3 = c 3 are displayed in Figure 7. It is observed that the function values display high peaks near (c 1 , c2 , c 3 ) and (c 1 , ĉ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). At lowest frequency k 1 = 1.0, the identification function 1/||P v(x)|| 2 is evaluated in the whole domain Ω to obtain the locations of D 1 , D 2 and D 3 . Figure 11 (right) is surface plot of 1/||P v(x)|| 2 = 400. We display the function values 1/||P v(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).

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.

  (right) shows the surface plot of 1/||P v(x)|| 2 = 650, which represents the location of the inclusion D. The function values 1/||P v(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 ).

Example 2

 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, c2 = -0.4, and ĉ2 = 0.4. The coordinates (c 1 , c2 , c 3 ) and (c 1 , ĉ2 , c 3 ) represent the locations of D 1 and D 2 , i.e., in the equation (4.1), x1 = x 1c 1 , x2 = x 2 -c2 , x3 = x 3c 3 and x1 = x 1c 1 , x2 = x 2 -ĉ2 , x3 = x 3c 3 for D 1 and D 2 respectively.

Example 3

 3 In this example, we consider the case where there are three inclusions inside Ω as shown in Figure11(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.

Figure 1 :

 1 Figure 1: True inclusion (left) and its localization by the MUSIC algorithm (right).

Figure 2 :Figure 3 :

 23 Figure 2: Values of 1/||P v(x)|| 2 on the cross-sectional planes. (a): x 1 = c 1 ; (b): x 2 = c 2 ; (c): x 3 = c 3 .

Figure 4 :

 4 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).

Figure 5 :

 5 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.

Figure 6 :

 6 Figure 6: True inclusions (left) and localization by the MUSIC algorithm (right).

Figure 7 :

 7 Figure 7: Values of 1/||P v(x)|| 2 on the cross-sectional planes. (a): x 1 = c 1 ; (b): x 2 = c2 ; (c): x 2 = ĉ2 ; (d): x 3 = c 3 .

Figure 8 :

 8 Figure 8: The relative residual during the iteration process for k = 1.5, 2.1 and 4.2, respectively.

Figure 9 :

 9 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).

Figure 10 :

 10 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.

Figure 11 :

 11 Figure 11: True inclusions (left) and localization by the MUSIC algorithm (right).

Figure 12 :

 12 Figure 12: Values of 1/||P v(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 .

Figure 13 :

 13 Figure 13: The relative residual during the iteration process for k = 1.5, 2.1 and 3.6, respectively.

Figure 14 :

 14 Figure 14: Evolution of the reconstruction at k = 1.2, 2.4, 4.2, 5.7 (left to right, top to bottom).

Figure 15 :

 15 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.

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.