Optimization approach for the simultaneous reconstruction of the dielectric permittivity and magnetic permeability functions from limited observations

We consider the inverse problem of the simultaneous reconstruction of the dielectric permittivity and magnetic permeability functions of the Maxwell's system in 3D with limited boundary observations of the electric field. The theoretical stability for the problem is provided by the Carleman estimates. For the numerical computations the problem is formulated as an optimization problem and hybrid finite element/difference method is used to solve the parameter identification problem.


Introduction
This paper is focused on the numerical reconstruction of the dielectric permittivity coefficient ε(x) and the magnetic permeability coefficient µ(x) for Maxwell's system basing our observation on a single measurement data of the electric field E(x, t). That means that we use boundary measurements of E(x, t) which are generated by a single direction of a plane wave. For the numerical reconstructions of the dielectric permittivity coefficient ε(x) and the magnetic permeability coefficient µ(x), we consider a similar hybrid finite element/difference method (FE/FDM) as was developed in [3].
In the literature, the stability results for Maxwell's equations have been proposed by Dirichlet-to-Neumann map or by Carleman estimates. For results using Dirichlet-to-Neumann map with infinitely many boundary observations, see [9,12,13,31,41,46]. For 1 results with a finite number of observations several Carleman estimates have been derived in [10,21,30,34,35,36]. In order to solve the coefficient inverse problem of Maxwell's equations numerically, we consider only those theoretical results that involve finite number of observations.
Since real observations are generally corrupted by noise, it is important to verify whether close observations lead to close estimations of the coefficients. The results of this paper give such conditions on the observations for estimating the parameters µ(x) and ε(x) of the Maxwell's system (1) from the observations of E(x, t) on the boundary of the domain. In particular, we obtain a stability inequality of the form for f (s) such that lim s→0 f (s) = 0, which links the distance between two sets of coefficients with the distance between two sets of boundary observations of the electric field E(x, t). Such stability inequalities lead to the uniqueness of the coefficients (µ(x), ε(x)) given the observation "Observation(E(x, t))" on a small neighborhood of the boundary of the domain of interest. They are also useful for the numerical reconstruction of the coefficients using noise-free observations [27]. Our main theoretical result concerns stability inequality which gives an estimate of the norm of two coefficients ε and µ in terms of observation of only the electric field E(x, t) on the boundary of the domain. This implies directly a uniqueness result. In the domain of the inverse problems, associated to reconstructing the Maxwell's coefficients from finite number of observations, except the reference [30], and up to our knowledge, there exists no result involving only one component E(x, t) (or the component H(x, t)).
In this work, the minimization problem of reconstructing functions ε(x) and µ(x) is reformulated as the problem of finding a stationary point of a Lagrangian involving a forward equation (the state equation), a backward equation (the adjoint equation) and two equations expressing that the gradients with respect to the coefficients ε(x) and µ(x) vanish. Moreover, in our work the forward and adjoint problems are given by time-dependent Maxwell's equations for the electric field. This means that we have observations of the electric field in space and in time which provides us better reconstruction of both coefficients. In order to get the computed values of ε(x) and µ(x), we arrange an iterative process by solving in each step the forward and backward equations and updating the coefficients ε(x) and µ(x) at every step of our iterations.
Recall, that in our optimization procedure the forward and adjoint problems are given by the time-dependent Maxwell's system for the electric field. For the numerical solution of the Maxwell equations, different formulations are available. One can consider, for example, the edge elements of Nédélec [40], the node-based first-order formulation of Lee and Madsen [33], the node-based curl-curl formulation with divergence condition of Paulsen and Lynch [42] or the interior-penalty discontinuous Galerkin FEM [26]. In this work, for the discretization of the Maxwell's equations, we use stabilized domain decomposition method of [3] with divergence condition of Paulsen and Lynch [42] which removes spurious solutions when local mesh refinement is applied and material discontinuities are not too big [5]. Numerical tests of [3] in the case of coefficient inverse problems (CIPs), similar to our of recent experimental work [7], show that these spurious solutions will not appear.
It is well known that edge elements are the most satisfactory from a theoretical point of view [37] since they automatically satisfy the divergence free condition. However, in the case of time-dependent computations they are less attractive. First, they are time consuming because the solution of a resulted linear system is required at every time iteration. Second, in [22] was shown that in the case of triangular or tetrahedral edge elements, the entries of the diagonal matrix resulting from mass-lumping are not necessarily strictly positive, and thus, explicit time stepping cannot be used in general. The method of [3] where nodal elements were used leads to a fully explicit and efficient scheme when mass-lumping is applied [29,22]. This method is efficiently implemented in the software package WavES [51] in C++ using PETSc [43] and MPI (message passing interface) and is convenient for our simulations.
We note, that we reconstruct functions ε(x) and µ(x) simultaneously. Up to our knowledge, there are few studies (e.g. [25]) that consider the simultaneous recovery of ε(x) and µ(x) within the spatio-temporal Maxwell's equations. In the papers [7,3], which use the same optimization approach, the coefficient µ is assumed to have a known and a constant value, i.e. µ = 1, which means that the medium is non-magnetic. However, it is well known that in many situations, we need to deal with magnetic materials (e.g. metamaterials [48], low-loss materials [45]).
Potential applications of our algorithm are in reconstructing the electromagnetic parameters in nanocomposites or artificial materials [45,47,48], imaging of defects and their sizes in a non-destructive testing of materials and in photonic crystals [15], measurement of the moisture content [14] and drying processes [38], for example.
Our numerical simulations show that we are able to accurately reconstruct simultaneously contrasts for both functions ε(x) and µ(x) as well as their locations. In our future work, similarly with [7,3,6], we are planning also to reconstruct shapes of the inclusions using a posteriori error estimates in the Tikhonov functional or in the reconstructed coefficients and based on them an adaptive finite element method.
The paper is organized as follows. In Section 2 we state the theoretical results of recovering ε(x) and µ(x) from the limited boundary measurements of E(x, t). Section 3 is devoted to presenting the numerical method used in this article. In Section 4, we give detailed information about our discrete numerical method and we outline the algorithm for the solution of our inverse problem. The numerical results are presented in Section 5. Discussions and conclusions are given in Section 6.

Theoretical results
Let us consider a bounded domain Ω ⊂ R 3 with a smooth boundary ∂Ω, and define by Ω T := Ω × (0, T ), ∂Ω T := ∂Ω × (0, T ), where T is a strictly positive constant. The electromagnetic equations in an inhomogeneous isotropic case in the bounded domain Ω ⊂ R d , d = 2, 3 with boundary ∂Ω, are described by the first order system of partial differential equations where E(x, t), H(x, t), D(x, t), B(x, t) are three-dimensional vector-valued functions of the time t and the space variable x = (x 1 , x 2 , x 3 ), and correspond to the electric and magnetic fields and the electric and magnetic inductions, respectively. The dielectric permittivity, ε(x) > 0 and the magnetic permeability, µ(x) > 0, depend on x ∈ Ω, ν = ν(x) denotes the unit outward normal vector to ∂Ω.
Our goal is to reconstruct the coefficients ε(x) and µ(x) in the system (1) with appropriate initial conditions E 0 and H 0 on the electric and magnetic inductions, using only a finite number of observations of the electric field E(x, t) on the boundary ∂Ω of the domain Ω.
We base our theoretical approach on the work [10]. In this work the authors obtained a stability inequality for the dielectric permittivity ε(x) and the magnetic permeability µ(x) involving boundary observations of both magnetic induction B(x, t) and electric induction D(x, t) (see theorem 1 of [10]). With the same assumptions, as those in [10], we can reformulate this stability result using only the observations of the electric field E(x, t) (or correspondingly, using only the observations of the magnetic field H(x, t)): Assume that the functions µ(x) and ε(x) in C 2 (Ω), x ∈ Ω obey for some ε 0 > 0 and µ 0 > 0. Next, for simplicity, we introduce some similar notations to the ones introduced in [10]. Let us pick x 0 ∈ R 3 \Ω, set c(x) = (µ(x)ε(x)) −1 for x ∈ Ω, c 0 = (µ 0 ε 0 ) −1 and assume that the following condition holds for some ρ ∈ (0, c 0 ) This technical condition is claimed by the weight function ψ 0 = |x − x 0 | 2 that is used to create the Carleman estimate established to prove the theorem 2.1 in [10]. In other terms, (3) arises from the classical pseudo-convexity condition. Another standard hypothesis is that the coefficients ε and µ are known in a neighborhood of the boundary of Ω Next, we define ω = Ω ∩ O where O is some neighborhood of ∂Ω in R 3 . Further let M 0 > 0 and two given functions µ ♯ , λ ♯ belong to C 2 (ω). Now we can define the admissible set of unknown coefficients µ and ε as Λ ω (M 0 ) = (µ, ε) verifying (2); (µ, ε) C 2 (Ω) ≤ M 0 and (µ, ε) = (µ ♯ , ε ♯ ) in ω .
Further, for the identification of (µ, ε), imposing (as will appear in the sequel) that (B(x, t), D(x, t)) are observed twice, we consider two sets of initial data (D k 0 , B k 0 ), k = 1, 2 such that, and define the 12 × 6 matrix where x ∈ Ω. Then we can write that ( . The extension of the time interval to (-T,T) corresponds to a technical point in the proof of the stability inequality (see lemma 3 in [10]). Now we recall the main theoretical result of the paper [10]: Under some hypothesis on T and choosing initial conditions (B k 0 , D k 0 ), k = 1, 2 verifying some additional assumptions, then there are two constants C > 0 and κ ∈ (0, 1), depending on Ω, ω, T , M and M 0 , such that we have: Under the same assumptions and considering the definition of the electric induction D i and the magnetic induction B i for i = 1, 2 we can write in the neighborhood ω of the boundary ∂Ω the following relations: Since ε i = ε ♯ , for i = 1, 2 in ω, we can write It is straightforward to verify that We define by the Hilbert space equipped with the norm .Then from (4) and since B k . Thus we can deduce our theorem in the following form for some M > 0. Then there are two constants C > 0 and κ ∈ (0, 1), depending on Ω, ω, T , M and M 0 , such that we have: 3 Statement of the forward and inverse problems

The mathematical model
Below for any vector function u ∈ R 3 our notations u ∈ L 2 (Ω) or u ∈ H k (Ω), k = 1, 2, mean that every component of the vector function u belongs to this space. Next, we decompose Ω into two subregions, Ω FEM and Ω FDM such that Ω = Ω FEM ∪Ω FDM , Ω FEM ∩ Ω FDM = ∅ and ∂Ω FEM ⊂ ∂Ω FDM , for an illustration of the domain decomposition, see figure 1. In Ω FEM we use finite elements and in Ω FDM we will use finite difference method with first order absorbing boundary conditions [24]. The boundary ∂Ω is such that  By eliminating B and D from (1) we obtain the model problem for the electric field E with the perfectly conducting boundary conditions at the boundary ∂Ω is as follows: Here we assume that We note that a similar equation can be derived also for H. For numerical solution of (6)- (9) in Ω FDM we use the finite difference method on a structured mesh with constant coefficients ε = ε ♯ = 1 and µ = µ ♯ = 1.
In Ω FEM , we use finite elements on a sequence of unstructured meshes K h = {K}, with elements K consisting of triangles in R 2 and tetrahedra in R 3 satisfying maximal angle condition [11].
In this work, for the discretization of the Maxwell's equations we use stabilized domain decomposition method of [3] and consider Maxwell's system in a convex geometry without reentrant corners where material discontinuities are not too big [5]. Since in our numerical simulations the relative permeability ε and relative permittivity µ does not vary much, such assumptions about the coefficients are natural. Moreover, recent experimental works [7,8,49] show that material discontinuities in real applications of reconstruction of dielectrics (or refractive indices) of unknown targets placed in the air or ground using microwave imaging technology are not big. For application where the coefficients are assumed to be smooth see [20] and [52]. [20] considers breast imaging and [52] soil moisture imaging. In our work we treat discontinuities of coefficients ǫ and µ in a similar way as in formula (3.25) of section 3.5 in [17]. See also [4] for the case when only ǫ is unknown. Further, in our computations all materials with values of ε > 10 are treated as metals and we call ε as "apparent" or "effective" dielectric constant", see [7,8,49] for more information and explanation.
To stabilize the finite element solution using standard piecewise continuous functions, we enforce the divergence condition (7) and add a Coulomb-type gauge condition [1,39] Let S T := ∂ 1 Ω × (0, T ) where ∂ 1 Ω is the backscattering side of the domain Ω with the time domain observations, and define by S 1,1 := ∂ 1 Ω × (0, t 1 ], S 1,2 := ∂ 1 Ω × (t 1 , T ), We use the Neumann boundary conditions at the left and right hand sides of a domain Ω (recall, that Ω = Ω FEM ∪ Ω FDM and in Ω FDM our coefficients ε(x) = µ(x) = 1) and first order absorbing boundary conditions [24] at the rest of the boundaries. Application of Neumann boundary conditions allows us to assume infinite structure in lateral directions and thus, allows us to consider the CIPs of the reconstruction of unknown parameters ε and µ in a waveguide. Further, we assume homogeneous initial conditions.
It was demonstrated numerically in [3] that the solution of the problem (14) approximates well the solution of the original Maxwell's equations for µ = 1, s = 1. The energy estimate of Theorem 4.1 of [3] implies an uniqueness result for the forward problem (14) with ε ≥ 1, µ = 1. Similar theorem can be proven for the case when both coefficients are unknown, and it can be done in a future research.
We assume that our coefficients ε (x) , µ(x) of equation (14) are such that The values of constants d 1 , d 2 in (15) are chosen from real life experiments similarly with [7,8,48,49] and we assume that we know them a priori. We consider the following Inverse Problem (IP) Suppose that the coefficients ε (x) and µ(x) satisfies (15) such that numbers d 1 , d 2 > 1 are given. Assume that the functions ε (x) , µ(x) are unknown in the domain Ω\Ω FDM . Determine the functions ε (x) , µ(x) for x ∈ Ω\Ω FDM , assuming that the following functionẼ (x, t) is known A priori knowledge of an upper and lower bounds of functions ε (x) and µ(x) corresponds well with the inverse problems concept about the availability of a priori information for an ill-posed problem [2,23,50]. In applications, the assumption ε (x) = µ(x) = 1 for x ∈ Ω FDM means that the functions ε (x) and µ(x) have a known constant value outside of the medium of interest Ω\Ω FDM . The functionẼ (x, t) models time dependent measurements of the electric wave field at the backscattering boundary ∂ 1 Ω of the domain of interest. In practice, measurements are performed on a number of detectors, see [8,49].

Optimization method
We reformulate our inverse problem as an optimization problem, where we seek for two functions, the permittivity ε(x) and permeability µ(x), which result in a solution of equations (14) with best fit to time and space domain observationsẼ, measured at a finite number of observation points on ∂ 1 Ω. Our goal is to minimize the Tikhonov functional whereẼ is the observed electric field, E satisfies the equations (14) and thus depends on ε and µ, ε 0 is the initial guess for ε and µ 0 is the initial guess for µ, and γ i , i = 1, 2 are the regularization parameters. Here z δ (t) is a cut-off function, which is introduced to ensure that the compatibility conditions at Ω T ∩ {t = T } for the adjoint problem (22) are satisfied, and δ > 0 is a small number. We choose a function z δ such that Next, we introduce the following spaces of real valued vector functions To solve the minimization problem, we introduce the Lagrangian where u = (E, λ, ε, µ) ∈ U 1 , and search for a stationary point with respect to u satisfying ∀ū = (Ē,λ,ε,μ) ∈ U 1 L ′ (u;ū) = 0, where L ′ (u; ·) is the Jacobian of L at u. We assume that λ (x, T ) = ∂ t λ (x, T ) = 0 and seek to impose such conditions on the function λ that L (E, λ, ε, µ) := L (u) = F (E, ε, µ) . Next, we use the fact that λ(x, T ) = ∂λ ∂t (x, T ) = 0 and E(x, 0) = ∂E ∂t (x, 0) = 0, as well as µ = ε = 1 on ∂Ω, together with boundary conditions ∂ n E = 0 and ∂ n λ = 0 on S 3 . The equation (17) expresses that for allū, Further, we obtain two equations that express that the gradients with respect to ε and µ vanish: The equation (18) is the weak formulation of the state equation (14) and the equation (19) is the weak formulation of the following adjoint problem We note that the adjoint problem is solved backward in time.

Numerical method 4.1 Finite element discretization
We discretize Ω FEM × (0, T ) denoting by K h = {K} a partition of the domain Ω FEM into tetrahedra K (h = h(x) being a mesh function, defined as h| K = h K , representing the local diameter of the elements), and we let J k be a partition of (0, T ) into time intervals J = (t k−1 , t k ] of uniform length τ = t k − t k−1 . We assume also a minimal angle condition on the K h [11].
To formulate the finite element method, we define the finite element spaces V h , W E h and W λ h . First we introduce the finite element trial space W E h for every component of the electric field E defined by where P 1 (K) and P 1 (J) denote the set of linear functions on K and J, respectively. We also introduce the finite element test space W λ h defined by Hence, the finite element spaces W E h and W λ h consist of continuous piecewise linear functions in space and time, which satisfy certain homogeneous initial and first order absorbing boundary conditions.
To approximate functions µ(x) and ε(x) we will use the space of piecewise constant functions V h ⊂ L 2 (Ω), where P 0 (K) is the piecewise constant function on K.
Next, we define Usually dim U h < ∞ and U h ⊂ U 1 as a set and we consider U h as a discrete analogue of the space U 1 . We introduce the same norm in U h as the one in U 0 , · U h := · U 0 . This means that in finite dimensional spaces all norms are equivalent and in our computations we compute coefficients in the space V h . The finite element method now reads: Find u h ∈ U h , such that

Fully discrete scheme
We expand E and λ in terms of the standard continuous piecewise linear functions {ϕ i (x)} M i=1 in space and {ψ k (t)} N k=1 in time and substitute them into (14) and (22) to obtain the following system of linear equations: with initial conditions : Here, M is the block mass matrix in space, K is the block stiffness matrix corresponding to the rotation term, C, D are the stiffness matrices corresponding to the divergence terms, S k is the load vector at time level t k , E k and λ k denote the nodal values of E(·, t k ) and λ(·, t k ), respectively, τ is the time step. Let us define the mapping F K for the reference elementê such that F K (ê) = e and letφ be the piecewise linear local basis function on the reference elementê such that ϕ • F K =φ. Then the explicit formulas for the entries in system (23) at each element e can be given as: where (·, ·) e denotes the L 2 (e) scalar product.
To obtain an explicit scheme, we approximate M with the lumped mass matrix M L (for further details, see [18]). Next, we multiply (23) with (M L ) −1 and get the following explicit method: Finally, for reconstructing ε(x) and µ(x) we can use a gradient-based method with an appropriate initial guess values ε 0 and µ 0 . The discrete versions of the gradients with respect to coefficients ε and µ in (20) and (21), respectively, take the form: Here, λ h and E h are computed values of the adjoint and forward problems using explicit scheme (24), and ε h , µ h are approximated values of the computed coefficients.

The algorithm
In this algorithm we iteratively update approximations ε m h and µ m h of the function ε h and µ h , respectively, where m is the number of iteration in our optimization procedure. We denote are computed by solving the state and adjoint problems with ε := ε m h and µ := µ m h .

Algorithm
Step 0. Choose the mesh K h in Ω and time partition J of the time interval (0, T ) . Start with the initial approximations ε 0 h = ε 0 and µ 0 h = µ 0 and compute the sequences of ε m h , µ m h via the following steps: Step ( 14) and adjoint (22) problems on K h and J.
Step 2. Update the coefficient ε h := ε m+1 h and µ h := µ m+1 h on K h and J using the conjugate gradient method , where α i , i = 1, 2, are step-sizes in the gradient update [44] and

Numerical Studies
In this section we present numerical simulations of the reconstruction of two unknown functions ε(x) and µ(x) inside a domain Ω FEM using the algorithm of section 4.3. These functions are known inside Ω FDM and are set to be ε(x) = µ(x) = 1. The goal of our numerical tests is to reconstruct two magnetic metallic targets of figure 2 with µ = 2.0. We note that when metallic targets are presented then our model problem (14) is invalid, see discussion about it [8,49]. This is one of the discrepancies between our mathematical model (14) and the simulated backscattering data. We refer to [49] for the description of other discrepancies in a similar case. However, one can treat metallic targets as dielectrics with large dielectric constants and it was shown computationally using experimental data in [8,32,49]. Similarly with [8,32,49] we call these large dielectric constants as apparent or effective dielectric constants and choose values for them in the interval ε (metallic target) ∈ (10, 30) .
For the relative magnetic permeability we choose values for them on the interval µ ∈ [1, 3], see [48]. In our studies, we initialize only one component E 2 of the electrical field E = (E 1 , E 2 , E 3 ) as the boundary condition in (14) on S T ( see (27)). Initial conditions are set to be zero. In all computations we used modification of the stabilized domain decomposition method of [3] which was implemented using the software package WavES [51] with two non-constant functions ε(x) and µ(x).
The computational geometry Ω is split into two geometries, Ω FEM and Ω FDM such that Ω = Ω FEM ∪ Ω FDM , see figure 1. Next, we introduce dimensionless spatial variables x ′ = x/ (1m) and obtain that the domain Ω FEM is transformed into dimensionless computational domain The dimensionless size of our computational domain Ω for the forward problem is The space mesh in Ω FEM and in Ω FDM consists of tetrahedral and cubes, respectively. In the optimization algorithm we choose the mesh size h = 0.1 in our geometries in the hybrid FEM/FDM method, as well as in the overlapping regions between FEM and FDM domains. In all our computational tests, we choose in (10) the penalty factor s = 1 in Ω FEM .
Note that in Ω FDM because of the domain decomposition method and conditions (15), the Maxwell's system transforms to the wave equation We initialize only one component of the electrical field E 2 as a plane wave f (t) in Ω in time T = [0, 1.2] such that while other two components E 1 , E 3 are initialized as zero. Thus, in Ω FDM we solve the problem (26) and in Ω FEM we have to solve Here, ∂Ω FDM I is internal boundary of the domain Ω FDM , and ∂Ω FEM is the boundary of the domain Ω FEM . Similarly, in Ω FDM the adjoint problem (22) transforms to the wave equation Thus, in Ω FDM we solve the problem (28) and in Ω FEM we have to solve We define exact functions ε(x) = 12 and µ(x) = 2 inside two small inclusions, see Figure  2, and µ(x) = ε(x) = 1 at all other points of the computational domain Ω FEM . We choose in our computations the time step τ = 0.003 which satisfies the CFL condition [19] and run computations in time [0, 1.2].
We consider the following test cases for the generation of the backscattering data: i) frequency ω = 21 with 3% additive noise ii) frequency ω = 21 with 10% additive noise iii) frequency ω = 30 with 3% additive noise iv) frequency ω = 30 with 10% additive noise To generate backscattering data at the observation points at S T in each cases i)-iv), we solve the forward problem (14), with function f (t) given by (27) in the time interval t = [0, 1.2] with the exact values of the parameters ε(x) = 12.0, µ(x) = 2 inside scatterers of figure 2, and ε(x) = µ(x) = 1.0 everywhere else in Ω. We avoid the variational crime in our tests since the data were generated on a locally refined mesh where inclusions were presented. However, our optimization algorithm works on a different structured mesh with the same mesh size h = 0.1.
The isosurfaces of the simulated exact solution of the initialized component E 2 (x, t) of the electrical field E(x, t) in the forward problem (14) with ω = 30 at different times are presented in figure 3. Using this figure we observe the backscattering wave field of the component E 2 (x, t).
We start the optimization algorithm with guess values of the parameters ε(x) = 1.0 , µ(x) = 1.0 at all points in Ω. Such choice of the initial guess provides a good reconstruction for both functions ε(x) and µ(x) and corresponds to starting the gradient algorithm from the homogeneous domain, see also [2,3] for a similar choice of initial guess. Using (25) the minimal and maximal values of the functions ε(x) and µ(x) in our computations belongs to the following sets of admissible parameters The solution of the inverse problem needs to be regularized since different coefficients can correspond to similar wave reflection data on ∂ 1 Ω. We regularize the solution of the inverse problem by starting computations with two different regularization parameters γ 1 = 0.01, γ 2 = 0.9 in (16). Our computational studies have shown that such choices for the regularization parameters are optimal in our case. We choose the regularization parameters in a computational efficient way such that the values of regularization parameters give the smallest reconstruction error given by the relative L 2 error e ε = ||ε−ε h || ||ε h || for the reconstructed ε and e µ = ||µ−µ h || ||µ h || for the reconstructed µ. Here, ε, µ are exact values of the coefficients and ε h , µ h are computed ones. We refer to [2,23,28] for different techniques for the choice of regularization parameters. The tolerance θ in our algorithm (section 4.3) is set to θ = 10 −6 . Figure 4 shows a case of backscattering data without presence of the additive noise. Figures 5 and 6 present typical behavior of noisy backscattering data with ω = 21 and ω = 30, respectively. Figure 7 presents a comparison between computed components E 2 and E 3 of the backscattering data with 10 % additive noise for both frequencies ω = 21 and ω = 30. Figure 8 presents the differences in backscattering data between 3% and 10% additive noise for both considered frequencies, ω = 21 on the left and ω = 30 on the right in figure 8.
To get images of figures 9 -15, we use a post-processing procedure. Suppose that functions ε n (x) and µ l (x) are our reconstructions obtained by algorithm of section 4.3 where n and l are number of iterations in gradient method when we have stopped to compute ε(x) and µ(x). Then to get images in figures 9 -15, we set

Discussion and Conclusion
In this work we have used time dependent backscattering data to simultaneously reconstruct both coefficients, ε(x) and µ(x), in the Maxwell's system as well as their locations. In order to do that we have used optimization approach which was similar to the method used in [4]. We tested our algorithm with two different noise levels (3% and 10% of additive noise) and with two different frequencies (ω = 21 and ω = 30, see (27)). The bigger noise level (10%) seemed to produce artefact in reconstructing ε with frequency ω = 30, see figure 15. However, we are able to reconstruct functions ε(x) and µ(x) with contrasts that are within the limits of (29). An important observation is that in our computations, we are able to obtain large contrasts for dielectric function ε(x) which allows us to conclude that we are          able to reconstruct metallic targets. At the same time, the contrast for the function µ(x) is within limits of (29). We could reconstruct size on z-direction for ε, however, size for µ(x) should be still improved. In our future research, we are planning to refine the obtained images through the adaptive finite element method in order to get better shapes and sizes of the inclusions. In [7,4,6] it was shown that this method is powerful tool for the reconstruction of heterogeneous targets, their locations and shapes accurately.
We note that our algorithm of simultaneous reconstruction of both parameters as well as its adaptive version can also be applied for the case when usual edge elements are used for the numerical simulation of the forward and adjoint problems in step 1 of our algorithm, see [16,17,25] for finite element analysis and discretization in this case. This as well as comparison between two different discretization techniques for the solution of our CIP can be considered as a challenge fot the future research.