SOME REMARKS ON THE SMALL ELECTROMAGNETIC INHOMOGENEITIES RECONSTRUCTION PROBLEM

. This work considers the problem of recovering small electromagnetic inhomogeneities in a bounded domain Ω ⊂ R 3 , from a single Cauchy data, at a ﬁxed frequency. This problem has been considered by several au- thors, in particular in [4]. In this paper, we revisit this work with the objective of providing another identiﬁcation method and establishing stability results from a single Cauchy data and at a ﬁxed frequency. Our approach is based on the asymptotic expansion of the boundary condition derived in [4] and the extension of the direct algebraic algorithm proposed in [1].

1. Introduction. In this paper, we consider the inverse problem of recovering small electromagnetic inhomogeneities, first in the framework of Helmholtz's equation and then in full Maxwell's equations.
Let Ω ⊂ R 3 be an open bounded domain with a Lipschitz boundary Γ. Suppose that Ω contains m small inhomogeneities D j of the form (1) D j = S j + B j j = 1, · · · , m where S j ∈ Ω are the locations of these inhomogeneities and B j ⊂ R 3 are bounded domains containing the origin. The constant is the common order of magnitude of the diameters of D j supposed to be a sufficiently small positive real number, taken without loss of generality, smaller than 1. Moreover, the inhomogeneities D j are assumed away from the boundary and their centers S j are supposed to be mutually distinct. Then, the electric potential u, in the case of the presence of these inhomogeneities, is the solution of the Helmholtz equation with a given Dirichlet boundary data f where ω > 0 is a given frequency and µ and ρ are respectively the permeability and the permittivity coefficients, defined as Here µ j and ρ j , j = 0, · · · , m, are positive known constants. The inverse problem, we are interested in, consists in reconstructing the inhomogeneities' number m, their centers S j and some characteristics related to and B j using only a single Cauchy data (f, g) := (u |Γ , ∂u ∂ν |Γ ) and at a fixed frequency ω. Here, ν denotes the outward unit normal to Γ.
The main objective of this paper is the development of an effective novel direct algebraic identification process for the reconstruction of the inhomogeneities' parameters using the high order terms of the complete asymptotic expansion of ∂u ∂ν | Γ when → 0 given in [4].
Several methods were developed in the literature for inhomogeneities reconstruction based on the boundary asymptotic expansion of the solution. Indeed, the development of the asymptotic expansion for this type of inverse problems was initially proposed by Friedman and Vogelius in [13] for the recovery of small inhomogeneities with extreme conductivity in Poisson's equation (ω = 0). Later, a similar asymptotic formula for steady-state voltage potential, in the presence of finite well-separated small inhomogeneities, was developed in [7] accompanied with a method based on the choice of the boundary voltage potential and a least-square approach to reconstruct the location and the size of these inhomogeneities. Then, for the reconstruction of electromagnetic inhomogeneities and based on the leading order terms of an asymptotic expansion of the normal derivative developed in [19], the authors in [3] used an approach with inverse Fourier transforms for identifying the order of magnitude and the center of these inhomogeneities using realistic wave sources. Their method exploited a number of boundary measurements obtained using different direction vectors on a sphere. In [5], the authors have adopted a MUSIC-type algorithm, [16], for detecting the number and the location of the corrosive parts in a domain based on the development of an asymptotic expansion of the boundary perturbations in terms of the size of the corrosive parts. In their paper, the algorithm reconstructed some characteristics of the inhomogeneities employing a number of boundary measurements based on the use of multiple frequencies. Furthermore, other identification methods were proposed in [6,8,9,17,18,20]. In these methods, the authors used either multiple boundary measurements and/or multiple frequencies to identify the inhomogeneities parameters.
The main novelty presented in this paper lies in the proposed identification method and the established stability result, using only a single Cauchy data at a fixed frequency. This method is employed priorly in the problem of sources reconstruction [1] and is based on the use of a reciprocity gap functional and specific test functions that lead to a number of algebraic relationships. Specifically, in here, the related developed relationships are utilized in order to recover the parameters of the domain inhomogeneities. Moreover, it generalizes and extends the work in [15] .
The remainder of this paper is organized as follows. In Section 2, we review some basic functions and recall the asymptotic expansion of ∂u ∂ν | Γ used, later, for the development of the identification method. In Section 3, we develop the algebraic relationships related to the inhomogeneities reconstruction. Then, Section 4 shows the corresponding identification process and algorithm employed for the identification of the unknown parameters. Later, an application over Maxwell's equation is adopted in Section 5. Finally, the stability of the inhomogeneities' positions is studied in Section 6. 2. Asymptotic formula. Before developing our identification method to restore the inhomogeneities, we need to state several notations and results essential for achieving the required relationships.
First, let u 0 be the solution of (2) in the absence of inhomogeneities. In other words, defining κ 0 as ( In order to ensure the existence of a solution u 0 and consequently that of u, solution of (2), one assumes that, [19], κ 2 0 is not an eigenvalue of the operator −∆ with Dirichlet boundary conditions. Now, for X = (x, y, z), consider G(X, S) as the Dirichlet Green function of ∆+κ 2 0 in Ω, namely, for all S ∈ Ω, G is the solution of (4) where δ S stands for the Dirac distribution at the point S.
Then, based on the calculations in [4], the asymptotic expansion of ∂u ∂ν over Γ when → 0 is computed as with S j , and B j respectively are the center, the diameter and the support of the inhomogeneities D j defined in (1) and ψ is a function that is evaluated using a single layer potential operator and another function denoted ϕ . Actually, for all X ∈ Γ and each multi-index , the couple (ϕ , ψ ) is the unique solution to ∂ν j where ν j denotes the outward unit normal to ∂B j and the subscripts − and + indicate the limiting values as the point approaches ∂B j from outside and inside B j respectively. Moreover, for any bounded domain B in R 3 and any positive constant κ, S κ B is the single layer potential defined over R 3 as with Θ κ being the fundamental solution of ∆ + κ 2 i.e. Θ κ is explicitly given as For more details about the asymptotic expansion and its coefficients, refer to [4]. Now, we are able to present our direct identification algorithm.
3. Algebraic relationships. The principal aim of this section is to establish algebraic relationships that permit to relate the Cauchy data (f, g) := (u |Γ , ∂u ∂ν |Γ ) to the parameters (m, S j , , B j ) of the inhomogeneities D j . This is accomplished by using the approach utilized in [1]. Indeed, we multiply (3) with special test functions, integrate by parts, use Green's formula and finally employ the asymptotic expansion (5) to reach our objective. As a matter of fact, in order to apply this algebraic approach, we proceed as follows. First, consider the homogenous space (7) H κ0 = {v ∈ H 1 (Ω) : ∆v + κ 2 0 v = 0} and define the following opertor Then, multiplying equation (3) by a function v in H κ0 , integrating over Ω and using Green's formula, we get Furthermore, using the value of ∂u 0 ∂ν in the asymptotic expansion (5) leads to On the other hand, from (4), one has Hence, deriving the latter equation β times with respect to S and replacing it in (9), the algebraic relationships (9) become The idea behind our identification method, as developed in [1], consists in choosing particular test functions v in H κ0 , as v a n (X) = (x + iy) n e iκ0z , X = (x, y, z) ∈ Ω.
Then, replacing v a n by its value in (10), for β = β 1 + β 2 + β 3 , we get To simplify this equation, we take α = β 1 + β 2 to get (11) c a n : Hence, using the latter equation, we can find, modulo 6 , the xy− projections of the inhomogeneities' locations. Before solving the equations (11), we need to know if the projections P a j are mutually distinct. Indeed, one can remark that there is only a finite number of planes containing the origin such that at least two points among S j are projected onto the same point on this plane. So, if a basis is chosen randomly, one is almost sure that the positions S j are projected onto distinct points. Therefore, without loss of generality, we assume that: (H) The projections onto the xy, yz and xz-planes of the points S j are mutually distinct.
Before presenting our identification algorithm we will introduce some notations and definitions that are used throughout this paper.
First, denote by P b j and P c j the projections of S j onto the yz and xz-complex planes respectively. Then, using in (10) the following test functions v b n = (y + iz) n e ik0x and v c n = (x + iz) n e ik0y and since v b n , v c n ∈ H κ0 , the general algebraic relationships are written as where η α,a j is defined in (12), Thus, the number, the 2−dimensional projections and other characteristics of the homogeneities are reconstructed, modulo 6 , based on the resolution of the algebraic relationships (13), whose identification algorithm is obtained in [1] and is recalled briefly in Section 4 for the convenience of the reader. 4. Identification algorithm. The identification process proposed to solve the algebraic relations (13) is attained in two steps. The first step consists in determining the number of the inhomogeneities based on the rank of a Hankel matrix reconstructed using the known measurements c r n . Then, the second step consists in determining the projections P r j of the centers S j by means of the eigenvalues of a companion matrix reconstructed using the calculated number of inhomogeneities.
As a matter of fact, being given an upper boundJ of the number and rewriting the equations (13) in a matrix form, we obtain the desired parameters as mentioned in the following theorems.
Theorem 4.1. LetJ ∈ N * be a known upper bound of 5m. Let H r J be the Hankel matrix defined as where c r n , n = 0, · · · , 2J − 2 are the coefficients defined in (13). Assuming (H) is satisfied, then, the number of the inhomogeneities is obtained as m = 1 5 rank(H r J ) if and only if η 4,r j = 0, ∀j = 1, · · · , m, r = a, b, c.
Proof. The proof is similar to that of [1, Theorem 2.8]. It is based on the decomposition of the Hankel matrix H r J given in the Lemma 2.9 in [1]. We present a sketch of the proof in the following.
Step 1. The Hankel matrix H r J can de decomposed as H r J = A rĪ r (A r ) t where the complex matrix A r is of sizeJ × J and defined as (15) A r = (U 0,r , · · · , U 4,r ) with, for α = 0, · · · , 4, U α,r being the confluentJ × m Vandermonde matrices andĪ r is the multi-diagonal matrix defined as Step 2. From (H), we can check that rank(A r ) t = J andĪ r is a nonsingular matrix if and only if η 4,r j , j = 1, · · · , m, are non null. This implies that with non-zero η 4,r j , I r (A r ) t is surjective. Consequently, we have rank(A rĪ r (A r ) t ) = rank(A r ) and since rank(A r ) = J, we obtain the needed result.
Assuming (H) is satisfied, then, 1. B r admits m eigenvalues of multiplicity 5. 2. The m eigenvalues are the projections P r j of the centers S j . Proof. The proof of this theorem is similar to [1, Theorem 2.10]. We present in the following a sketch of the proof.
Step 2. One can easily derive the relations: Here, the matrix A r is invertible since the projected points P r j are assumed distinct. Moreover, since rank(H r J ) = J, the family (ξ r n ) n=0,...,J−1 is a basis of C J . Therefore, writing explicitly A r T r (A r ) −1 in this basis coincides with B r defined in (16) and thus we obtain the needed results.

5.
Inhomogeneities reconstruction over Maxwell's equation. The methodology used above can be also employed for the case of recovering the parameters of inhomogeneities in Maxwell's equation. Indeed, as in the Helmholtz case, the idea is to reconstruct the number, the centers, and some characteristics related to the diameter and the support of the inhomogeneities D j of form (1), based on the use of the asymptotic expansion of the solution over the boundary when → 0 and on a direct reconstruction method. In here, our work is based on the methods developed in [10,11]. More precisely, Subsection 5.1 is devoted to the statement of the problem. Then, in Subsection 5.2, we recall the key elements needed in Subsection 5.3 to develop the related algebraic reconstruction method. Finally, Subsection 5.4 presents briefly the direct algorithm used to recover the needed parameters.  (1), where each inhomogeneity is of magnetic permeability µ j , electric permittivityρ j and electric conductivity σ j . Then, the electric fieldĒ and the magnetic fieldB satisfy the full Maxwell's equations given as with µ j ,ρ j , and σ j , j = 0, · · · , m, being known positive constants where µ 0 ,ρ 0 and σ 0 are respectively the magnetic permeability, the electric permittivity and the electric conductivity of the background medium.
In here, we consider the time-harmonic Maxwell solutions. These solutions are of the formĒ where ω > 0 is a given frequency and E and B satisfy the equations Therefore, to obtain the equation satisfied by the electric field E, we divide equation (19) by the so-called complex permittivity ρ , defined as, ρ =ρ + i σ ω and then we take the curl of equation (18). Thus, in this case, the electric field E satisfies the system where f is the given tangential component of E.
The goal of the inverse problem is to reconstruct the number of inhomogeneities m, their locations S j , j = 1, · · · , m, and some characteristics related to and B j from a single given Cauchy data (f, g) = ((E × ν) |Γ , (∇ × E × ν) |Γ ) prescribed on a sufficiently regular boundary Γ of Ω, at a fixed frequency.

Asymptotic formula.
As achieved in the Helmholtz case, in order to develop the algebraic direct method used to identify the inhomogeneities, we need to define several notations and functions indispensable for the identification algorithm.
First, let E 0 be the electric field of (20) in the absence of inhomogeneities. In other words, setting κ 0 as Moreover, consider the 3 × 3 matrix-valued function G(X, S) the solution of Here, I denotes the 3 × 3 identity matrix and the operator ∇ X × applies to the matrices, column by column. Then, based on the calculations in [2], the asymptotic expansion of (∇ × E) × ν over Γ when → 0 is computed as where for all j = 1, · · · , m, r j = (r j,1 , r j,2 , r j, with, for a constant c, each polarization tensor M j (c) is a 3 × 3, symmetric, positive definite matrix, [13], associated with the jth inhomogeneity where each element of the matrix M j is given by where δ i j is the Kronecker delta and for 1 ≤ j ≤ 3, φ j is the unique function that satisfies with φ j continuous across ∂B j and lim |y|→+∞ φ j (y) = 0. Here, {e j } 3 j =1 is an orthonormal basis of R 3 , ν j denotes the outward unit normal to ∂B j and the subscripts − and + indicate the limiting values as the point approaches ∂B j from outside and inside B j respectively.

Remark 2.
Note that in this case, the leading order of the boundary perturbations resulting from the presence of these inhomogeneities is of modulo ( 4 ) rather than ( 6 ), which was obtained in the case of Helmholtz equation treated above. In fact, in order to generalize the obtained asymptotic expansion (23), one needs further techniques that we don't have in hand for the moment. Therefore, in this paper, the identification method for Maxwell's equation is developed modulo ( 4 ).

Algebraic relations.
In this subsection, our aim is to achieve the algebraic relationships that permit us to recover the parameters of the inhomogeneities using the test functions chosen in [11].
As mentioned in the Helmholtz case, the algebraic relations (29), reconstruct, modulo 4 , only the xy− projections of the inhomogeneities locations. To determine their 3D positions, we proceed in the same way considering the test functions Then, the general algebraic relationships, in this case, are written as where P a j , P b j and P c j are the projections of S j respectively on the xy−, yz− and xz− complex planes and λ a j is defined in (30), λ b j = 3 e −iκ0xj (r j,2 + ir j,3 ) + κ 0 (p j,2 + ip j,3 ) and λ c j = 3 e −iκ0yj (r j,1 + ir j,3 ) + κ 0 (p j,1 + ip j,3 ) .

Reconstruction algorithm.
The reconstruction algorithm of our developed direct method is based on using [10, Theorem 2] and proceeding using the following steps.
Step 1. Suppose the existence of an upper boundJ of the number of inhomogeneities m. Then, using the Cauchy data pair (f, g) on the boundary Γ, compute c r 0 , c r 1 , · · · , c r 2J−1 . Therefore, according to the proof of [10,Theorem 2], m can be determined as the rank of one of the three Hankel matrices H r J , r = a, b, c, defined in (14) with coefficients c r n defined in (31). The rank of these matrices can be numerically estimated, for example, using the Singular Value Decomposition method with an appropriate threshold, following [14].
Step 2. Solve the linear system H r m C r = (c r m , c r m+1 , · · · , c r 2m−1 ) t . Thus, the 2D projection points P r j of the inhomogeneities' positions are obtained as the m simple eigenvalues of the companion matrix B r , defined in (16) and for which the obtained C r form the last line.
Step 3. The coefficients λ r j and consequently characteristics related to and B j are obtained by solving the systems A r (λ r 1 , λ r 2 , · · · , λ r m ) t = (c r 0 , c r 1 , · · · , c r m−1 ) t where A r is the Vandermonde matrix defined as 6. Stability of the inhomogeneities' centers. In this section, we study the stability of the inhomogeneities reconstruction in Helmholtz equation and then in Maxwell equation. The stability estimates for the position of the inhomogeneities over Helmholtz equation are based on the method proposed in [12]. Indeed, in the latter paper, the authors derived the stability estimates for the reconstruction of monopolar and dipolar sources. Their method, that we employ in here, is based on the choice a suitable test function using the 2D projection of two different configurations of the inhomogeneities. Then, using the reciprocity gap functional and suitable estimation inequalities, we obtain the needed results.  (1) with their centers S j = (x j , y j , z j ), j = 1, · · · , m. As denoted before, we take P a j = x j + iy j and P b j = y j + iz j as the projections of the centers onto xy-and yz-complex planes respectively, and set P a = (P a j ) 1≤j≤m and P b = (P b j ) 1≤j≤m . Then, we introduce the following real coefficients (34) a = min 1≤j,n≤m,j =n P a j − P a n and b = min 1≤j,n≤m,j =n n which henceforth will be called, respectively, the "separability coefficient" of the projected sources P a and P b and set Now, for two points configurations R = (R j ) 1≤j≤m for = 1, 2, we recall the Hausdorff distance between R 1 and R 2 , defined as follows 6.2. Stability in Helmholtz's case. In this subsection, our aim is to establish the stability estimates for the centers of the inhomogeneities in Helmholtz equation following the guidelines in [12]. Indeed, let (λ j, β , S j ) 1≤j≤m for = 1, 2, β = 0, · · · , 4, be two inhomogeneities configurations with S j = (x j , y j , z j ) for = 1, 2 and λ j, β are the coefficients defined in (6). Moreover, let c 1 and c 2 be the following constants: where is defined in (33). Then, we obtain a Hölder stability estimate of the centers of the inhomogeneities as shown in the following theorem. Theorem 6.1. Let, u for = 1, 2 be the solutions of (2) corresponding to the inhomogeneities D j = S j + B j , j = 1, · · · , m characterized by the configurations (λ j, β , S j ) 1≤j≤m , with λ j, 4 = 0. Let (f , g ) := (u |Γ , ∂u ∂ν |Γ ) for = 1, 2, be the corresponding measurements on the boundary Γ. Then, assuming S j ∈ Ω α , defined in (32), the following estimate holds where is defined in (35) and c 1 , c 2 in (36).
To prove Theorem 6.1, we need to show the following results on the Hausdorff distance between the projections onto xy− and yz− complex planes of the configurations S 1 and S 2 . Lemma 6.2. Let, for = 1, 2, P a, = (P a, j ) 1≤j≤m and P b, = (P b, j ) 1≤j≤m , be, respectively, the corresponding projected point sources onto xy-and yz-complex planes of S = (S j ) 1≤j≤m . Under the assumptions of Theorem 6.1, we have , where a , b are defined in (34) and c 1 , c 2 in (36).
Ψ n belongs to the space H k0 (defined in (7)). Therefore, integrating by parts, using the reciprocity gap functional (8) and employing (10), we obtain Now, by taking the difference between the previous two sums, for all n = 1, · · · , m 2 , we get which can be rewritten, using the definitions of R and Φ n , as follows Moreover, using the separability coefficient a defined in (34) and Hölder estimate, we get from (38), where c 1 is defined in (36). Now, one has to estimate Ψ n L 2 (Γ) and ∂Ψ n ∂ν .
Indeed, according to the definition of in (33) and using Cauchy-Schwarz inequality, one obtains where |Γ| is defined in (33). Thus, where c 2 is defined in (36).
Repeating the same procedure by replacing the function Φ n by Φ n , Φ n (x, y) = (x+iy−P a,1 n ) 4 m 2

j=1
(x+iy−P a,2 j ) 5 m 1 j =n (x+iy−P a,1 j ) 5 for n = 1, · · · , m 1 leads to max 1≤n≤m 1 min 1≤j≤m 2 P a,1 n − P a,2 Finally, taking the maximum between (39) and (40), we get the desired first estimate. To prove the second one, it suffices to replace P a,1 j , P a,2 j , by P b,1 j , P b,2 j respectively and consider the following tests functions Now, thanks to Lemma 6.4, we can easily prove Theorem 6.1.
Proof. Proof of Theorem 6.1. Since . Using Lemma 6.4, we get the desired result. remains in the space H k0 . It is rather interesting to mention that, to reach a better identification of the centers of the inhomogeneities , it is desirable to project the point on a plane ( u, v) where the associated separability coefficient is the largest possible, which can also be seen in the dependence of the upper bound in the stability estimate. In practice, to attain such a plane, we can assume, for example, that u = (cos(φ) cos(θ), cos(φ) sin(θ), sin(φ)), v = (sin(φ) cos(θ), sin(φ) sin(θ), − cos(φ)) and then take the pair (φ, θ) ∈ [0, π 2 ] × [0, 2π] that realizes the largest m with the best separability coefficient.
6.3. Stability in Maxwell's case. To study the stability of the inhomogeneities' centers in Maxwell's equation (20), we follow the same methodology as proceeded in Helmholtz case and shown in what follows.
To prove Theorem 6.3, we need, as seen in Helmholtz's case, to show the following results on the Hausdorff distance between the projections onto xy− and yz− complex planes of the configurations S 1 and S 2 . Lemma 6.4. Let, for = 1, 2, P a, = (P a, j ) 1≤j≤m and P b, = (P b, j ) 1≤j≤m , be, respectively, the corresponding projected point sources onto xy-and yz-complex planes of S = (S j ) 1≤j≤m . Under the assumptions of Theorem 6.3, we have where a , b are defined in (34) and c 3 in (41).
Indeed, according to the definition of in (33) and using Cauchy-Schwarz inequality, one obtains Finally, taking the maximum between (45) and (46), we get the desired first estimate. To prove the second one, it suffices to replace P a,1 j , P a,2 j , by P b,1 j , P b,2 j respectively and consider the following tests functions Ψ n (x, y, z) = Φ n (y, z) Now, as in Theorem 6.1, thanks to Lemma 6.4, we can easily prove Theorem 6.3.

Remark 4.
From the stability estimates (37) and (42), one can remark the effect of different parameters on the reconstruction of the centers of the inhomogeneities. Indeed, one can see that the number of inhomogeneities m, their diameter and the separability coefficient between the projections have an impact on the identification process. As a matter of fact, it may be noted that as the number of sources increases and consequently the separability coefficient decreases (due to a fixed domain), the stability estimate is quite low. This is predictable since in the case of a large number of inhomogeneities, the interactions between them inside the domain become strong enough and therefore it will be more difficult to identify their positions. Furthermore, the diameter of the inhomogeneities plays a crucial role in the identification process. As seen in the stability estimates (37) and (42), taking a very small coefficient causes the degradation of the reconstruction method. Hence, one should find a compromise between the need to find more characteristics about the inhomogeneities by taking high orders of the asymptotic expansion and taking a suitable small enough . This causes difficulties in the numerical implementations of the identification process.
7. Conclusion. In this work, we adopted the asymptotic expansions obtained in [3,2] for the solutions of Helmholtz's equation and Maxwell's equation respectively in order to develop a simple direct identification process of small inhomogeneities. Different from prior methods, our method uses a single Cauchy data pair with a single frequency and enables us to reconstruct the number, the centers and some characteristics of the support and the diameter of these inhomogeneities. Finally, to study the effect of various parameters on the reconstruction method, Hölder stability estimates of the locations of the inhomogeneities are obtained.