An adaptive finite element/finite difference domain decomposition method for applications in microwave imaging

A new domain decomposition method for Maxwell's equations in conductive media is presented. Using this method reconstruction algorithms are developed for determination of dielectric permittivity function using time-dependent scattered data of electric field. All reconstruction algorithms are based on optimization approach to find stationary point of the Lagrangian. Adaptive reconstruction algorithms and space mesh refinement indicators are also presented. Our computational tests show qualitative reconstruction of dielectric permittivity function using anatomically realistic breast phantom.


Introduction
In this work are presented reconstructions algorithms for the problem of determination of the spatially distributed dielectric permittivity function in conductive media using scattered time-dependent data of the electric field at the boundary of investigated domain. Such problems are called Coefficient Inverse Problems (CIPs). A CIP for a system of timedependent Maxwell's equations for electric field is a problem about the reconstruction of unknown spatially distributed coefficients of this system from boundary measurements.
One of the most important application of algorithms of this paper is microwave imaging including microwave medical imaging and imaging of improvised explosive devices (IEDs). Potential application of algorithms developed in this work are in breast cancer detection.
In numerical examples of current paper we will focus on microwave medical imaging of realistic breast phantom provided by online repository [59]. In this work we develop simplified version of reconstruction algorithms which allow determine the dielectric permittivity function under the condition that the effective conductivity function is known. Currently we are working on the development of similar algorithms for determination of both spatially distributed functions, dielectric permittivity and conductivity, and we are planning report about obtained results in a near future.
Microwave medical imaging is non-invasive imaging. Thus, it is very attractive addition to the existing imaging technologies like X-ray mammography, ultrasound and MRI imaging. It makes use of the capability of microwaves to differentiate among tissues based on the contrast in their dielectric properties.
In [30] were reported different malign-to-normal tissues contrasts, revealing that malign tumors have a higher water/liquid content, and thus, higher relative permittivity and conductivity values, than normal tissues. The challenge is to accurately estimate the relative permittivity of the internal structures using the information from the backscattered electromagnetic waves of frequencies around 1 GHz collected at several detectors.
Since the 90-s quantitative reconstruction algorithms based on the solution of CIPs for Maxwell's system have been developed to provide images of the complex permittivity function, see [17] for 2D techniques, [15,18,31,38] for 3D techniques in the frequency domain and [49,56] for time domain (TD) techniques.
In all these works microwave medical imaging remained the research field and had little clinical acceptance [37] since the computations are inefficient, take too long time, and produce low contrast values for the inside inclusions. In all the above cited works local gradient-based mathematical algorithms use frequency-dependent measurements which often produce low contrast values of inclusions and miss small cancerous inclusions. Moreover, computations in these algorithms are done often in MATLAB, sometimes requiring around 40 hours for solution of inverse problem.
It is well known that CIPs are ill-posed problems [2,32,53,55]. Development of non-local numerical methods is a main challenge in solution of a such problems. In works [6,7,51,52] was developed and numerically verified new non-local approximately globally convergent method for reconstruction of dielectric permittivity function. The twostage global adaptive optimization method was developed in [6] for reconstruction of the dielectric permittivity function. The two-stage numerical procedure of [6] was verified in several works [7,51,52] on experimental data collected by the microwave scattering facility.
The experimental and numerical tests of above cited works show that developed meth-2 ods provide accurate imaging of all three components of interest in imaging of targets: shapes, locations and refractive indices of non-conductive media. In [38], see also references therein, authors show reconstruction of complex dielectric permittivity function using convexification method and frequency-dependent data. Potential applications of all above cited works are in the detection and characterization of improvised explosive devices (IEDs).
The algorithms of the current work can efficiently and accurately reconstruct the dielectric permittivity function for one concrete frequency using single measurement data generated by a plane wave.
A such plane wave can be generated by a horn antenna as it was done in experimental works [7,51,52]. We are aware that conventional measurement configuration for detection of breast cancer consists of antennas placed on the breast skin [1,18,19,37,49]. In this work we use another measurement set-up: we assume that the breast is placed in a coupling media and then the one component of a time-dependent electric plane wave is initialized at the boundary of this media. Then scattered data is collected at the transmitted boundary. This data is used in reconstruction algorithms developed in this work. Such experimental set-up allows avoid multiply measurements and overdetermination since we are working with data resulted from a single measurement. An additional advantage is that in the case of single measurement data one can use the method of Carleman estimates [33] to prove the uniqueness of reconstruction of dielectric permittivity function.
For numerical solution of Maxwell's equations we have developed finite element/finite difference domain decomposition method ( FE/FD DDM).
This approach combines the flexibility of the finite elements and the efficiency of the finite differences in terms of speed and memory usage as well as fits the best for reconstruction algorithms of this paper. We are unaware of other works which use similar set-up for solution of CIP for time-dependent Maxwell's equations in conductive media solved via FE/FD DDM, and this is the first work on this topic.
An outline of the work is as follows: in section 2 we present the mathematical model and in section 3 we describe the structure of domain decomposition. Section 4 presents reconstruction algorithms including formulation of inverse problem, derivation of finite element and finite difference schemes together with optimization approach for solution of inverse problem. Section 5 shows numerical examples of reconstruction of dielectric permittivity function of anatomically realistic breast phantom at frequency 6 GHz of online repository [59]. Finally, section 6 discusses obtained results and future research.

The mathematical model
Our basic model is given in terms of the electric field E (x,t) = (E 1 , E 2 , E 3 ) (x,t) , x ∈ R 3 changing in the time interval t ∈ (0, T ) under the assumption that the dimensionless relative magnetic permeability of the medium is µ r ≡ 1. We consider the Cauchy problem for the Maxwell equations for electric field E (x,t), further assuming that that the electric volume charges are equal zero, to get the model equation for x ∈ R 3 ,t ∈ (0, T ]. Here, ε r (x) = ε(x)/ε 0 is the dimensionless relative dielectric permittivity and σ (x) is the effective conductivity function, ε 0 , µ 0 are the permittivity and permeability of the free space, respectively, and c = 1/ √ ε 0 µ 0 is the speed of light in free space. We are not able numerically solve the problem (1) in the unbounded domain, and thus we introduce a convex bounded domain Ω ⊂ R 3 with boundary ∂ Ω. For numerical solution of the problem (1), a domain decomposition finite element/finite difference method is developed and summarized in Algorithm 1 of section 3.
A domain decomposition means that we divide the computational domain Ω into two subregions, Ω FEM and Ω FDM such that Ω = Ω FEM ∪ Ω FDM with Ω FEM ⊂ Ω, see Figure 2. Moreover, we will additionally decompose the domain Ω FEM = Ω IN ∪ Ω OUT with Ω IN ⊂ Ω FEM such that functions ε r (x) and σ (x) of equation (1) should be determined only in Ω IN , see Figure 2. When solving the inverse problem IP this assumption allows stable computation of the unknown functions ε r (x) and σ (x) even if they have large discontinuities in Ω FEM .
The communication between Ω FEM and Ω FDM is arranged using a mesh overlapping through a two-element thick layer around Ω FEM , see elements in blue color in Figure 1a),b). This layer consists of triangles in R 2 or tetrahedrons in R 3 for Ω FEM , and of squares in R 2 or cubes in R 3 for Ω FDM .
The key idea with such a domain decomposition is to apply different numerical methods in different computational domains. For the numerical solution of (1) in Ω FDM we use the finite difference method on a structured mesh. In Ω FEM , we use finite elements on a sequence of unstructured meshes K h = {K}, with elements K consisting of tetrahedron's in R 3 satisfying minimal angle condition [34].
We assume in this paper that for some known constants d 1 > 1, d 2 > 0, the functions ε r (x) and σ (x) of equation (1) satisfy Turning to the boundary conditions at ∂ Ω, we use the fact that (2) and (1) imply that since ε r (x) = 1, σ (x) = 0 for x ∈ Ω FDM ∪ Ω OUT , then a well known transformation makes the equations (1) independent on each other in Ω FDM , and thus, in Ω FDM we need to solve the equation We write ∂ Ω = ∂ Ω 1 ∪ ∂ Ω 2 ∪ ∂ Ω 3 , meaning that ∂ Ω 1 and ∂ Ω 2 are the top and bottom sides of the domain Ω, while ∂ Ω 3 is the rest of the boundary. Because of (4), it seems natural to impose first order absorbing boundary condition for the wave equation [22], Here, we denote the outer normal derivative of electrical field on ∂ Ω by ∂ · ∂ n , where n denotes the unit outer normal vector on ∂ Ω.
It is well known that for stable implementation of the finite element solution of Maxwell's equation divergence-free edge elements are the most satisfactory from a theoretical point of view [40,43]. However, the edge elements are less attractive for solution of timedependent problems since a linear system of equations should be solved at every time iteration. In contrast, P1 elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix [20,29]. It is also well known that numerical solution of Maxwell equations using nodal finite elements can be resulted in unstable spurious solutions [41,46]. There are a number of techniques which are available to remove them, see, for example, [26-28, 42, 46].
In the domain decomposition method of this work we use stabilized P1 FE method for the numerical solution of (1) in Ω FEM . Efficiency of usage an explicit P1 finite element scheme is evident for solution of CIPs. In many algorithms which solve electromagnetic CIPs a qualitative collection of experimental measurements is necessary on the boundary of the computational domain to determine the dielectric permittivity function inside it. In this case the numerical solution of time-dependent Maxwell's equations are required in the entire space R 3 , see for example [6,7,11,51,52], and it is efficient to consider Maxwell's equations with constant dielectric permittivity function in a neighborhood of the boundary of the computational domain. An explicit P1 finite element scheme with σ = 0 in (1) is numerically tested for solution of time-dependent Maxwell's system in 2D and 3D in [3]. Convergence analysis of this scheme is presented in [4] and CFL condition is derived in [5]. The scheme of [3] is used for solution of different CIPs for determination of dielectric permittivity function in non-conductive media in time-dependent Maxwell's equations using simulated and experimentally generated data, see [7,11,51,52].
The stabilized model problem considered in this paper is: with functions ε r , σ satisfying conditions (2).

The domain decomposition algorithm
We now describe the domain decomposition method between two domains Ω FEM and Ω FDM where FEM is used for computation of the solution in Ω FEM , and FDM is used in Ω FDM , see Figures 1,2. Overlapping nodes between Ω FDM and Ω FEM are outlined in Figure 2 by green circles (boundary nodes of Ω FEM ) and blue diamonds (inner boundary nodes of Ω FDM ). The communication between two domains Ω FEM and Ω FDM is achieved by overlapping of both meshes across a two-element thick layer around Ω FEM -see Figure 2. The nodes of the computational domain Ω belong to either of the following sets (see Figure 2 By conditions (2) functions ε r = 1 and σ = 0 at the overlapping nodes between Ω FEM and Ω FDM , and thus, the Maxwell's equations will transform to the system of uncoupled acoustic wave equations (4) which leads to the fact that the FEM and FDM discretization schemes coincide on the common structured overlapping layer. In this way we avoid instabilities at interfaces in the domain decomposition algorithm.

Reconstruction algorithms
In this section we develop different optimization algorithms which allow determination of the relative dielectric permittivity function using scattered data of the electric field at the boundary of the investigated domain. In all algorithms we use assumption that the effective conductivity function is known in the investigated domain.
In summary, the main algorithms presented in this section are: • Algorithm 2: The domain decomposition algorithm for efficient solution of forward and adjoint problems used in algorithms 3, 4, 5.
• Algorithm 3: Optimization algorithm for determination of the relative dielectric permittivity function under condition that the effective conductivity function is known.
• Algorithm 4, 5: Adaptive optimization algorithms for determination of the relative dielectric permittivity function. These algorithms use local adaptive mesh refinement based on a new error indicators for improved determination of location, material and sizes of the inclusions to be identified.
Let the boundary ∂ Ω = ∂ Ω out FDM ∪ ∂ Ω in FDM be the outer boundary ∂ Ω out FDM of Ω together with the inner boundary ∂ Ω in FDM of Ω FDM , and ∂ Ω FEM be the boundary of Ω FEM . Let at S T := ∂ Ω out FDM × (0, T ) we have time-dependent backscattering observations. Our coefficient inverse problem will be the following. Inverse Problem (IP) Assume that the functions ε r (x), σ (x) satisfy conditions (2) for known d 1 > 1, d 2 > 0. Let the function ε r be unknown in the domain Ω\(Ω FDM ∪ Ω OUT ). Determine the function ε r (x) for x ∈ Ω\(Ω FDM ∪ Ω OUT ), assuming that the function σ (x) is known in Ω and the following functionẼ (x,t) is measured at S T : The functionẼ (x,t) in (7) represents the time-dependent measurements of all components of the electric wave field E( To solve IP we minimize the corresponding Tikhonov functional and use a Lagrangian approach to do that. We present details of derivation of optimization algorithms in the next section.

Derivation of optimization algorithms
For solution of the IP for Maxwell's system (6) it is natural to minimize the following Tikhonov functional whereẼ is the observed electric field in (7) at the observation points located at ∂ Ω out FDM , δ obs = ∑ δ (∂ Ω out FDM ) is a sum of delta-functions at the observations points located at ∂ Ω out FDM , E satisfies the equations (6) and thus depends on ε r , σ . We denote by ε 0 the initial guess for ε r , and by γ the regularization parameter. Here, z δ is a cut-off function ensuring the compatibility conditions for data, see details in [11].
Let us introduce the following spaces of real valued functions To solve the minimization problem we take into account conditions (2) on the function ε r and introduce the Lagrangian where u = (E, λ , ε r ).
To solve the minimization problem (10) we find a stationary point of the Lagrangian with respect to u satisfying ∀ū = (Ē,λ ,ε r ) where L (u; ·) is the Jacobian of L at u. For solution of the minimization problem (12) we develop conjugate gradient method for reconstruction of parameter ε r . To obtain optimality conditions from (12), we integrate by parts in space and time the Lagrangian (11), assuming that λ (x, T ) = ∂ λ ∂t (x, T ) = 0, ∂ λ ∂t = ∂ λ ∂ n , and impose such conditions on the function λ that L(E, λ , ε r ) := L(u) = J(E, ε r ). Using the facts that λ (x, T ) = ∂ λ ∂t (x, T ) = 0, ∇ · (ελ ) = 0 and σ = 0, ε r = 1 on ∂ Ω, together with initial and boundary conditions of (6), we get following optimality conditions for allū ∈ U 1 , Finally, we obtain the main equation for iterative update ε r in the conjugate gradient algorithm which express that the gradient with respect to ε r vanishes: The equation (13) is the weak formulation of the forward problem (6) and the equation (14) is the weak formulation of the following adjoint problem

The domain decomposition FE/FD method for solution of forward and adjoint problems 4.2.1 Finite element discretization
We denote by where l is the total number of elements K in Ω FEM .
Here, h = h(x) is a piecewise-constant mesh function defined as representing the local diameter of the elements. We also denote by ∂ K h = {∂ K} a partition of the boundary ∂ Ω FEM into boundaries ∂ K of the elements K such that vertices of these elements belong to ∂ Ω FEM . We let J τ be a partition of the time interval (0, T ) into time intervals J = (t k−1 ,t k ] of uniform length τ = T /N for a given number of time steps N. We assume also a minimal angle condition on the K h [10,34]. To formulate the finite element method in Ω for (12) we define the finite element spaces C h , W E 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) denote the set of piecewise-linear functions on K.
To approximate function ε r we define the space of piecewise constant functions C h ⊂ L 2 (Ω), The finite element method for (12) now reads: find 11 The equation (19) expresses discretized versions of optimality conditions given by (13)- (15). To get function ε r via optimality condition (15) we need solutions first of the forward problem (6), and then of the adjoint problem (16). To solve these problems via the domain decomposition method, we decompose the computational domain Ω = Ω FEM ∪ Ω FDM as it is described in section 3. Thus, in Ω FEM we have to solve the following forward problem: Here, g is the solution obtained by the finite difference method in Ω FDM which is saved at ∂ Ω FEM . The equation (19) expresses that the finite element method in Ω FEM for the solution of the forward problem (20) will be: Here, we define (6) in Ω FEM .
To get the discrete scheme for (21) we approximate E h (kτ) by E k h for k = 1, 2, ..., N using the following scheme for k = 1, 2, . . . , Rearranging terms in (22) we get for k = 1, 2, . . . , The adjoint problem in Ω FEM will be the following: The finite element method for the solution of adjoint problem (24) in Ω FEM reads: Find (24) in Ω FEM . We note that the adjoint problem should be solved backwards in time, from time t = T to t = 0. To get the discrete scheme for (25) we approximate λ h (kτ) by λ k h for k = N, N − 1, ..., 1 using the following scheme for k = N − 1, . . . , 1: Multiplying both sides of (26) by τ 2 c 2 /ε rh and rearranging the terms we obtain: We note that usually dimU 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 , where U 0 is defined in (9). From (28) follows that all norms in finite dimensional spaces are equivalent. This allows us in numerical simulations of section 5 compute the discrete function ε rh , which is approximation of ε r (x), in the space C h .

Fully discrete scheme in Ω FEM
In this section we present schemes for computations of the solutions of forward (6) and adjoint (16) problems in Ω FEM . After expanding functions E h (x) and λ h (x) in terms of the standard continuous piecewise linear functions where E h i and λ h i denote unknown coefficients at the mesh point (23) and (27), correspondingly, withλ (x,t) =Ē(x,t) = ∑ M j=1 ϕ j (x), and obtain the system of linear equations for computation of the forward problem (6): Here, M, M 1 , M 2 are the assembled block mass matrices in space, G 1 , G 2 , G 3 are the assembled block matrices in space, F k is the assembled load vector at the time iteration k, E k denote the nodal values of E h (·,t k ), τ is the time step. Now we define the mapping F K for the reference elementK such that F K (K) = K and letφ be the piecewise linear local basis function on the reference elementK such that ϕ • F K =φ. Then, the explicit formulas for the entries in system of equations (29) at each element K can be given as: where (·, ·) K denotes the L 2 (K) scalar product and ∂ K is the part of the boundary of element K which lies at ∂ Ω FEM . For the case of adjoint problem (27) we get the system of linear equations: 14 Here, M, M 1 , M 2 , G 1 , G 2 , G 3 are the assembled block matrices in space with explicit entries given in (30), and P k 1 , P k 2 are assembled load vectors at the time iteration k with explicit entries λ k denote the nodal values of λ h (·,t k ), τ is the time step. Finally, for reconstructing ε r (x) in Ω IN we can use a gradient-based method with an appropriate initial guess values ε 0 . The discrete versions in space of the gradients given in (15), after integrating by parts in space of the third term in the right hand side of (15), have the form ∀x ∈ Ω IN : where ε 0 h is interpolant of ε 0 . We note that because of usage of the domain decomposition method, gradient (33) should be updated only in Ω IN since in Ω FDM and in Ω OUT by condition (2) we have ε r = 1, σ = 0. In (33) E h and λ h are computed values of the forward and adjoint problems using schemes (29), (31), correspondingly, and ε rh is approximate value of the computed relative dielectric permittivity function ε r .

Finite difference formulation
We recall now that from conditions (2) it follows that in Ω FDM the function ε r (x) = 1, σ = 0. This means that in Ω FDM the model problem (6) transforms to the following forward problem for uncoupled system of acoustic wave equations for E = (E 1 , E 2 , E 3 ): where ∂ E FEM ∂ n are known values at ∂ Ω in FDM .
Using standard finite difference discretization of the first equation in (??) in Ω FDM we obtain the following explicit scheme for every component of the solution E of the forward problem (??) with correspondingly discretized absorbing boundary conditions. In equations above, E k l, j,m is the finite difference solution on the time iteration k at the discrete point (l, j, m), τ is the time step, and ∆E k l, j,m is the discrete Laplacian. The adjoint problem in Ω FDM will be: where ∂ λ FEM ∂ n are known values at ∂ Ω in FDM . Similarly with (36) we get the following explicit scheme for the solution of adjoint problem (35) in Ω FDM which we solve backward in time: with corresponding boundary conditions. In equations (34), (36) (·) k l, j,m is the solution on the time iteration k at the discrete point (l, j, m), We note that we use FDM only inside Ω FDM , and thus computed values of ∂ E FEM ∂ n and ∂ λ FEM ∂ n can be approximated and will be known at ∂ Ω in FDM through the finite element solution in Ω FEM , see details in the domain decomposition Algorithm 2.

The domain decomposition algorithm to solve forward and adjoint problems
First we present domain decomposition algorithm for the solution of state and adjoint problems. We note that because of using explicit finite difference scheme in Ω FDM we need to choose time step τ accordingly to the CFL stability condition [4,5,14] such that the whole scheme remains stable.

Algorithm 2
The domain decomposition algorithm to solve forward and adjoint problems 1: Construct the finite element mesh K h in Ω FEM and the finite difference mesh in Ω FDM as well as time partition J τ of the time interval (0, T ). At every time step k we perform the following operations: 2: On the mesh in Ω FDM compute E k+1 , λ k−1 from (34), (36), correspondingly, using absorbing boundary conditions at the outer boundary ∂ Ω, with E k , E k−1 and λ k , λ k+1 known. 3: On the mesh K h in Ω FEM compute E k+1 , λ k−1 using the finite element schemes (29), (31), correspondingly, with E k , E k−1 and λ k , λ k+1 known. 4: Use the values of the functions E k+1 , λ k−1 at nodes ω * overlapping with nodes ω , which are computed using the finite element schemes (29), (31), correspondingly, as a boundary conditions at the inner boundary ∂ Ω in FDM for the finite difference method in Ω FDM . 5: Use the values of the functions E k+1 , λ k−1 at nodes ω o overlapping with nodes ω + , which are computed using the finite difference schemes (34), (36), correspondingly, as a boundary conditions at ∂ Ω FEM for the finite element method in Ω FEM . 6: Apply swap of the solutions for the computed functions E k+1 , λ k−1 . Set k = k + 1 for forward problem and k = k − 1 for adjoint problem and go to step 2.

Reconstruction algorithm for the solution of inverse problem IP
We use conjugate gradient method (CGM) for iterative update of approximation ε r m h of the function ε rh , where m is the number of iteration in the optimization algorithm. We introduce the following function where functions E m h , λ m h are computed by solving the state and adjoint problems with ε r := ε r m h , σ := σ m h .

Adaptive algorithms for solution of the inverse problem IP
Adaptive algorithm allows improvement of already computed relative dielectric permittivity function ε r M h obtained on the initially non-refined mesh in the previous optimization algorithm (Algorithm 3). The idea of the local mesh refinement (note that we need it only in Ω IN ) is that it should be refined in all neighborhoods of all points in the mesh K h where where α is the step-size in the gradient update [47] and Here, d 0 (x) = −g 0 h (x). the function |hε rh | achieves its maximum value, or where |J ε r (ε rh )| achieves its maximal values. These local mesh refinements recommendations are based on a posteriori error estimates for the error |ε r − ε rh | in the reconstructed function ε r ( see the first mesh refinement indicator), and for the error |J(ε r ) − J(ε rh )| in the Tikhonov's functional (see the second mesh refinement indicator), respectively. The proofs of these a posteriori error estimates for arbitrary Tikhonov's functional is given in [1]. A posteriori error for the Tikhonov's functional (8) can be derived using technique of [11], and it is a topic of ongoing research. Assuming that we have proof of these a posteriori error indicators, let us show how to compute them.
Assuming now that solutions E(ε r , σ ), λ (ε r , σ ) are sufficiently stable we can write that the Frechét derivative of the Tikhonov functional is the following function Inserting (15) into (39), we get In the second mesh refinement indicator is used discretized version of (40) computed for approximations (ε rh , σ h ).
• The First Mesh Refinement Indicator Refine the mesh in neighborhoods of those points of K h where the function |hε rh | attains its maximal values. In other words, refine the mesh in such subdomains of K h where Here, β ∈ (0, 1) is a number which should be chosen computationally and h is the mesh function (17) of the finite element mesh K h .
• The Second Mesh Refinement Indicator Refine the mesh in neighborhoods of those points of K h where the function |J ε r (E, ε rh )| attains its maximal values. More precisely, let β ∈ (0, 1) be the tolerance number which should be chosen in computational experiments. Refine the mesh K h in such subdomains where Remarks • 1. We note that in (41) exact values of E(x,t), λ (x,t) are used obtained with the already computed functions (ε rh , σ h ), see (40). However, in our algorithms and in computations we approximate exact values of E(x,t), λ (x,t) by the computed ones E h (x,t), λ h (x,t).
• 2. In both mesh refinement indicators we used the fact that functions ε r , σ are unknown only in Ω IN .
Algorithm 4 Adaptive Algorithm, first version 1: Construct the finite difference mesh in Ω FDM . Choose an initial space-time mesh K h0 × J τ 0 in Ω FEM × [0, T ]. Compute the sequence of ε rk , k > 0, via following steps: 2: Obtain numerical solution ε rk with known function σ k on K h k using the Algorithm 3 (Conjugate Gradient Method). 3: Refine such elements in the mesh K h k where the first mesh refinement indicator is satisfied. Here, the tolerance numbers β k ∈ (0, 1) are chosen by the user. 4: Define a new refined mesh as K h k+1 and construct a new time partition J τ k+1 such that the CFL condition is satisfied. Interpolate ε rk , σ k on a new mesh K h k+1 and perform steps 2-4 on the space-time mesh K h k+1 × J τ k+1 . Stop mesh refinements when ||ε rk − ε rk−1 ||< tol 1 or ||g k h (x)||< tol 2 , where tol i , i = 1, 2 are tolerances chosen by the user.
We define the minimizer of the Tikhonov functional (8) and its approximated finite element solution on k times adaptively refined mesh K h k by ε r and ε rk , correspondingly. In our both mesh refinement recommendations we need compute the functions ε rk on the mesh K h k . To do that we apply Algorithm 3 (conjugate gradient algorithm). We will define by ε rk := ε r M h values obtained at steps 3 of the conjugate gradient algorithm.
Algorithm 5 Adaptive Algorithm, second version 1: Choose an initial space-time mesh K h 0 × J τ 0 in Ω FEM . Compute the sequence ε rk , k > 0 with known σ k , on a refined meshes K h k via following steps: 2: Obtain numerical solutions ε rk on K h k × J τ k using the Algorithm 3 (Conjugate Gradient Method). 3: Refine the mesh K h k at all points where the second mesh refinement indicator is satisfied. Here, indicator g k h is defined in (37). Tolerance number β k ∈ (0, 1) should be chosen in numerical examples. 4: Define a new refined mesh as K h k+1 and construct a new time partition J τ k+1 such that the CFL condition is satisfied. Interpolate ε rk , σ k on a new mesh K h k+1 and perform steps 1-3 on the space-time mesh K h k+1 × J τ k+1 . Stop mesh refinements when ||ε rk − ε rk−1 ||< tol 1 , or ||g k h (x)||< tol 2 , where tol i , i = 1, 2 are tolerances chosen by the user.

Remarks
• 1. First we make comments how to choose the tolerance numbers β k , β k in (42), (43). Their values depend on the concrete values of max ingly. If we will take values of β k , β k which are very close to 1 then we will refine the mesh in very narrow region of the Ω IN , and if we will choose β k , β k ≈ 0 then almost all elements in the finite element mesh will be refined, and thus, we will get global and not local mesh refinement.
• 2. To compute L 2 norms ||ε rk − ε rk−1 ||, in step 3 of adaptive algorithms the reconstruction ε rk−1 is interpolated from the mesh K h k−1 to the mesh K h k .
• 3. The computational mesh is refined only in Ω FEM such that no new nodes are added in the overlapping elements between two domains, Ω FEM and Ω FDM . Thus, the mesh in Ω FDM , where finite diffirence method is used, always remains unchanged.

Numerical examples
In this section, we present numerical simulations of the reconstruction of permittivity function of three-dimensional anatomically realistic breast phantom taken from online repository [59] using an adaptive reconstruction Algorithm 4 of section (4.6). We have tested performance of an adaptive Algorithm 5 and it is slightly more computationally expensive in terms of time compared to the performance of Algorithm 4. Additionally, relative errors in the reconstructions of dielectric permittivity function are slightly smaller for Algorithm 4 and thus, in this section we present results of reconstruction for Algorithm 4.

Description of anatomically realistic data
We have tested our reconstruction algorithm using three-dimensional realistic breast phantom with ID = 012204 provided in the online repository [59]. The phantom comprises the structural heterogeneity of normal breast tissue for realistic dispersive properties of normal breast tissue at 6 GHz reported in [35,36]. The breast phantoms of database [59] are derived using T1-weighted MRIs of patients in prone position. Every phantom presents 3D mesh of cubic voxels of the size 0.5 × 0.5 × 0.5 mm. Tissue types and corresponding media numbers of breast phantoms are taken from [59] and are given in Table 1. Spatial distribution of these media numbers for phantom with ID = 012204 is presented in Figure 5. Figures 5-a)  demonstrate distribution of media numbers on the original coarse mesh consisting of 34 036 992 nodes. Clearly, performing computations on a such big mesh is computationally demanding task, and thus, we have sampled the original mesh. In all our computations we have used the mesh consisting of 63492 nodes as a coarse finite element mesh which was obtained by taking every 8-th node in x 1 , x 2 and x 3 directions of the original mesh. Figures 3-4 shows spatial distribution of dielectric permittivity ε r and effective conductivity σ (S/m) on original and sampled meshes. Testing of our algorithms on other sampled meshes is computationally expensive task, requiring running of programs in parallel infrastructure, and can be considered as a topic for future research.
We note that in all our computations we scaled original values of ε r and σ of database [59] presented in Figures 3-4 and considered weighted versions of these parameters, in order to satisfy conditions (2) as well as for efficient implementation of FE/FD DDM for solution of forward and adjoint problems. Table 1 presents weighted values of ε r and σ used in numerical tests of this section. Thus, in this way we get computational set-up corresponding to the domain decomposition method which was used in Algorithms 2-5.
The following model problem was used in all computations: We initialize a plane wave f (t) = (0, f 2 , 0)(t) for one component E 2 of the electric field (44). The function f 2 (t) represents the single direction of a plane wave which is initialized at ∂ 1 Ω in time t = [0, 3.0] and is defined as The goal of our numerical tests Test 1, Test 2 was to reconstruct weighted dielectric permittivity function ε r shown in Figures 7-a), b). Figures 9-a)-c), 10-a)-c) present simulated solution |E h | in Ω FEM of model problem (44) for Test 1 and Test 2, correspondingly.
To perform computations for solution of inverse problem, we add normally distributed Gaussian noise with mean µ = 0 to simulated electric field at the transmitted boundary ∂ 2 Ω. Then we have smoothed out this data in order to get reasonable reconstructions, see details of data-preprocessing in [51,52]. Computations of forward and inverse problems were done in time T = [0, 3] with equidistant time step τ = 0.006 satisfying to CFL condition. Thus, it took 500 timesteps at every iteration of reconstruction Algorithm 4 to solve forward or adjoint problem. The time interval T = [0, 3] was chosen computationally such that the initialized plane wave could reach the transmitted boundary ∂ 2 Ω in order to obtain meaningful reflections from the object inside the domain Ω FEM . Figures 8-a)-i), 9-a)-c), 10-a)-c) show these reflections in different tests. Experimentally such signals can be produced by a Picosecond Pulse Generator connected with a horn antenna, and scattered time-dependent signals can be measured by a Tektronix real-time oscilloscope, see [51,52] for details of experimental set-up for generation of a plane wave and collecting time-dependent data. For example, in our computational set-up, the experimental time step between two signals can beτ = 6 picoseconds and every signal should be recorded duringT = 3 nanoseconds.
We have chosen following set of admissible parameters for reconstructed function ε r (x) as well as tolerance θ = 10 −5 at step 3 of the conjugate gradient Algorithm 3. Parameters β k in the refined procedure of Algorithm 4 was chosen as the constant β k = 0.8 for all refined meshes K hk . These figures show that largest by amplitude reflections, or transmitted data, are obtained from the second component E 2 of the electric field E. The same observation is obtained in previous works [3,11] where was used a similar computational set-up with a plane wave. However, comparison of all three components was not presented in [11]. Domination of reflections at the transmitted boundary from the E 2 component can be explained by the fact that we initialize only one component of the electric field E = (E 1 , E 2 , E 3 ) as a plane wave f (t) = (0, f 2 , 0)(t) at Γ 1,1 in the model problem (44), and thus, two other components E 1 , E 3 will be smaller by amplitude than the E 2 when we use the explicit scheme (29) for computations. See also theoretical justification of this fact in [50].
Numerical tests of [11] show that the best reconstruction results of the space-dependent function ε r (x) for σ = 0 in Ω are obtained for ω = 40 in (45). Thus, we performed simulations of the forward problem (44) taking σ = 0 for different ω = 40, 60, 80, 100 in (45). It turned out that for chosen computational set-up with final time T = 3 maximal values of  scattered function E 2 are obtained for ω = 40. Thus, we take ω = 40 in (45) in all our tests. We assume that both functions ε r , σ satisfy conditions (2): they are known inside Ω out ∪ Ω FDM and unknown inside Ω IN . The goal of our numerical tests is to reconstruct the function ε r of the domain Ω FEM of Figure 7 under conditions (2) and the additional condition that the function σ (x) of this domain is known. See Table 1 for distribution of ε r , σ in Ω FEM .
The computational set-up for solution of inverse problem is as follows. We generate transmitted data by solving the model problem (44) on three times adaptively refined mesh. In this way we avoid variational crimes when we solve the inverse problem. The transmitted data is collected at receivers located at every point of the transmitted boundary ∂ 2 Ω, and then normally distributed Gaussian noise δ = 3%, 10% with mean µ = 0 is added to this data, see Figures 9-d)-f) -10-d)-f). The next step is data pre-processing: the noisy data is smoothed out, see Figures 9-g)-i) -10-g)-i). Next, to reconstruct ε r we minimize the Tikhonov functional (8). For solution of the minimization problem we introduce Lagrangian and search for a stationary point of it using an adaptive Algorithm 4, see details in section 4.6.
We take the initial approximation ε 0 = 1 at all points of the computational domain what corresponds to starting of our computations from the homogeneous domain. This is done because of previous computational works [11] as well as experimental works of [31,37,49] where was shown that a such choice gives good results of reconstruction of dielectric permittivity function.

Test 1
In this test we present numerical results of reconstruction of ε r when exact values of this function are given in Table 1, see Test 1. Isosurface of the exact function ε r to be reconstructed in this test is shown in Figure 7-a). We note that the exact function ε r has complicated structure. Using Figure 7-a) one can observe that isosurface presents a discontinuous function with a lot of big and small inclusions in the domain Ω FEM . Figures 11-a)-i) show results of the reconstruction on adaptively locally refined meshes when noise level in the data was δ = 10%. We start computations on a coarse mesh K h0 . Figure 11-a)-c) shows that the location of the reconstructed function ε h0 is imaged correctly and the reconstructed isosurface covers the domain where the exact ε r is located. We refer to Table 2 for the reconstruction of the maximal contrast in ε h0 . For improvement of the contrast and shape obtained on a coarse mesh K h0 , we run computations on locally adaptivelly refined meshes. Figures 11-d)-f) show reconstruction obtained on the final two times refined mesh K h2 . Table 2 presents results of reconstructions for ε hk obtained on the refined meshes K hk , k = 0, 1, 2. We observe that with mesh refinements we achieve better contrast for function ε r . Also reconstructed isosurface of this function more precisely covers the domain where the exact ε r is located, compare Figure 11-a) with Figure 11-d). Figures  11-g)-i) show locally adaptively refined mesh K h2 .

Test 2
Since it is quite demanding reconstruct very complicated structure of ε r taken in Test 1, in this test we will reconstruct ε r with exact isosurface as it is presented in the Figure 7-b). Exact values of this function are taken as in fibroconnective/glandular-1 media (see Table  1) inside isosurface of Figure 7-b), and outside of this isosurface all values of ε r = 1. Figures 12-a)-i) show results of the reconstruction on adaptively refined meshes when   Computational time δ = 3% δ = 10% Time (sec) Relative time n Test 1 110. 59 1.779 · 10 −6 71360 Test 2 106. 58 1.714 · 10 −6 69699 Time (sec) Relative time n Test 1 116. 22 1.869 · 10 −6 75052 Test 2 111. 53 1.793 · 10 −6 65359 Table 6: Performance of solution of forward problem (44) in Tests 1 and 2 of section 5 on the mesh K h0 in terms of computational time (in seconds) and relative computational time computed by (47). Here, n is number of the nodes on three times adaptively refined original coarse mesh (consisting of 63492 nodes) which we used for generation of transmitted data.
noise level in the data was δ = 10%. We refer to the Table 4 for reconstruction of the contrast in ε r . Using the Table 4 we now observe that with mesh refinements we achieve slightly higher maximal contrast 9.45 in reconstruction ε h1 compared to the exact one 9. Moreover, on the mesh K h1 for σ = 10% we get more than 8 times smaller relative error in the reconstruction compared to the error obtained on the coarse mesh K h0 . Figures 12-d)-i) show good matching of the reconstructed ε h1 compared with the exact one. Figures 11-j)-l) show locally adaptively refined mesh K h2 .

Performance comparison
All computations were performed on a linux workstation Intel Core i7-9700 CPU with one processor using software package WavES [57] efficiently implemented in C++/PETSc [48]. We have estimated the relative computational time T r of the forward problem using the following formula Here, t is the total computational time of the forward problem on the mesh K hl where l = 0, 1, 2, ... is number of the refined mesh, n is the total number of nodes in the mesh K hl , n t is number of timesteps. We take n t = 500 in all computational tests, see clarification in section 5.2. Computational times (in seconds) for solution of forward problem are presented in Table 6. Using this table we observe that the relative time is approximately the same for all tests and we can take it as T r ≈ 1.8 · 10 −6 . Next, using this relative time we can estimate approximate computational time for solution of forward problem for any mesh consisting of n nodes. For example, if we will take original mesh consisting of n = 34036992 nodes, then computational time will be already t = T r · n t · n = 1.8 · 10 −6 · 500 · 34036992 = 30633 seconds, and this time is not computationally efficient. Clearly, computing of the solution of inverse problem on the sampled mesh allows significantly reduce computational times. We have estimated also the relative computational time T ip r of the solution of inverse 30 problem using the formula Here, t ip is the total computational time to run inverse Algorithm 4 on the mesh K hl where l = 0, 1, 2, ... is number of the refined mesh, nno is the total number of nodes in the mesh K hl , n t is number of timesteps. Computational times (in seconds) for solution of inverse problem for Test 1 and Test 2 are presented in Tables 3,5, respectivelly. Using these tables we observe that computational times are depend on the number of iterations M k in the conjugate gradient method (CGM) and number of the nodes nno in the meshes K hl . We took n t = 500 for all tests and thus, computational times presented in these tables are not depend on number of times steps for different refined meshes. We note, that the number of time steps n t can be chosen adaptively as well. However, we are performing only adaptive mesh refinement in space and not in time. The full space-time adaptive algorithm can be considered as a topic for future research.
Using Table 3 we observe that computational time in Test 1 is around 20 minutes for both noise levels σ = 3% and σ = 10%. On every mesh K hl , l = 0, 1, 2, was performed two iterations CGM , or M K = 2. Thus, the total computational time to obtain final reconstruction in Test 1 is 60 min. Table 5 shows that computational time in Test 2 with noise in data δ = 3% is around 20 minutes for non-refined mesh K h0 , 60 min for one time refined mesh K h1 , and 20 minutes for twice refined mesh K h2 . Thus, the total computational time to obtain final reconstruction in Test 2 is 100 minutes. Computational time in this test is larger than in the previous Test 1 since CGM converged only at 5-th iteration on the one time refined mesh K h1 . However, the total computational time with noise in data δ = 10% is around 60 minutes. This is because the solution was obtained already on the one time refined mesh K h1 . Tables 3, 5 also demonstrate that it takes around 10 minutes to compute solution of inverse problem on the one iteration of the conjugate gradient algorithm.
We note that PETSc supports parallel implementation and thus, current version of code can be extended to the version with parallel implementation such that times reported in Tables 6, 3 and Table 5 can be significantly reduced.

Conclusions
This work describes reconstruction methods for determination of the relative dielectric permittivity function in conductive media using scattered data of the time-dependent electric field at number of detectors placed at the boundary of the investigated domain.
Reconstruction methods use optimization approach where a functional is minimized via a domain decomposition finite element/finite difference method. In an adaptive reconstruction method the space mesh is refined only in the domain where a finite element method is used with a feedback from a posteriori error indicators. Developed adaptive algorithms allow us to obtain correct values and shapes of the dielectric permittivity function to be determined. Convergence and stability analysis of the developed methods is ongoing work and will be presented in forthcoming publication. The algorithms of the current work are designed from previous adaptive algorithms developed in [7,11] which reconstruct the wave speed or the dielectric permittivity function. However, all previous algorithms are developed for non-conductive medium.
Our computational tests show qualitative and quantitative reconstruction of dielectric permittivity function using anatomically realistic breast phantom which capture the heterogeneity of normal breast tissue at frequency 6 GHz taken from online repository [59]. In all tests we used assumption that the conductivity function is known. Currently we are working on algorithms when both dielectric permittivity and conductivity functions can be reconstructed. Results of this work will be presented in our future research.
All computations are performed in real time presented in Tables 3, 5 and 6. Some data (Matlab code to read data of database [59], visualize and produce discretized values of ε r , σ , etc.) used in computations of this work is available for download and testing, see [58]. Additional data (computational FE/FD meshes, transmitted data, C++/PETSc code) can be provided upon request.
In summary, the main features of algorithms of this work are as follows: • Ability to reconstruct shapes, locations and maximal values of dielectric permittivity function of targets in conductive media under the condition that the conductivity of this media is a known function.
• More exact reconstruction of shapes and maximal values of dielectric permittivity function of inclusions because of local adaptive mesh refinement.
• Computational greater efficiency because of usage software package WavES [57] implemented in C++/PETSc [48].  Table 1 for breast phantom of object ID 012204 of database [59]. Table 1 clarifies description of media numbers and corresponding tissue types.   Reconstructions ε h2 ≈ 5 obtained on refined mesh K h2 . g)i): refined mesh K h2 . The noise level in the data is δ = 10%. See Table 2 for obtained contrasts max Ω FEM ε hk , k = 0, 1, 2. For comparison we also present exact isosurface with value corresponding to reconstructed one and outlined by red color.  Table 4 for obtained contrasts max Ω FEM ε hk , k = 0, 1. For comparison we also present exact isosurface of ε r with value corresponding to reconstructed one and outlined by red color.