Sensitivity analysis for 3D Maxwell's equations and its use in the resolution of an inverse medium problem at fixed frequency

ABSTRACT

This paper deals with the reconstruction of small-amplitude perturbations in the electric properties (permittivity and conductivity) of a medium from boundary measurements of the electric field at a fixed frequency. The underlying model is the three-dimensional time-harmonic Maxwell equations in the electric field. Sensitivity analysis with respect to the parameters is performed, and explicit relations between the boundary measurements and the characteristics of the perturbations are found from an appropriate integral equation and extensive numerical simulations in 3D. The resulting non-iterative algorithm allows to retrieve efficiently the centre and volume of the perturbations in various situations from the simple sphere to a realistic model of the human head.

KEYWORDS:

• Inverse medium problem
• 3D Maxwell equations
• small-amplitude inhomogeneities
• Gâteaux derivative
• integral equation
• edge finite elements
• computer science

AMS CLASSIFICATIONS:

• 65N21
• 65R32
• 78M10
• 65D18

1. Introduction

The study of dielectric properties of biological tissues or materials is of great interest in medical or industrial applications. The dielectric behaviour of a tissue and its interaction with electromagnetic fields are able to describe and provide information about its characteristics and composition. This information can be used to develop new non-invasive modalities in many practical applications of electric fields in agriculture, bioengineering and medical diagnosis. Dielectric properties of biological tissues are frequency-dependent or dispersive, and experimental investigation has shown that they vary with respect to low or high frequencies of the applied electric field or current (e.g. [1,2]). Among the different imaging modalities based on electric fields, we may cite Electrical Impedance Tomography (EIT) which operates at frequencies between 10 and 100 kHz and Microwave Imaging between 300 MHz and 300 GHz. The principle of EIT is to provide the electrical permittivity and conductivity inside a body from simultaneous measurements of electrical currents and potentials at the boundary. With regard to medical applications, microwave imaging is studied with the aim of detecting and monitoring cerebrovascular accidents (or strokes). Indeed, strokes result in variations of the dielectric properties of the affected tissues, and experimental research has found that the contrast of dielectric parameters between the abnormal and normal cerebral tissues can be imaged within the microwave spectrum (at frequencies of the order of 1 GHz). New devices based on these properties are currently designed and studied [3,4]. Microwave breast imaging offers also a promising alternative method to mammography [5]. Compared to other medical imaging technologies such as Magnetic Resonance Imaging (MRI) and Computerized Tomography (CT-scan), EIT and microwave imaging have a low resolution due to the ill-posed nature of the image reconstruction problem. There is a lot of interest (low-cost device, harmless procedure, etc.) in finding ways to improve their resolution and these modalities are areas of active research.

From a mathematical point of view, one has to deal with the theoretical and numerical study of an inverse medium problem. The goal is to retrieve the complex refractive index of a medium, namely the electric permittivity (real part) and conductivity (imaginary part) from boundary measurements at a fixed frequency. This inverse problem is severely ill-posed (e.g. [6]). Indeed, coefficients of elliptic problems (like the conductivity equation) in a bounded domain are uniquely determined by the entire (scalar or vector) Dirichlet-to-Neumann map on the whole boundary of the domain which is in general not available in practical applications. The fundamental example of parameter reconstruction is Calderón's inverse conductivity problem [7]. The theoretical and numerical study of the EIT inverse problem has also been extensively addressed in the last two decades. We refer for instance to [8–11] and references therein.

In this paper, we focus on the inverse medium problem associated with the time-harmonic 3D Maxwell equations formulated in the electric field with a possible application in microwave imaging. For uniqueness and stability results for the Maxwell system from total or partial data, we refer the reader for instance to the works of Ola et al. [12], Caro et al. [13,14], Kenig et al. [15] and references therein. In practice, only partial information on the (vector) Dirichlet-to-Neumann map is available. The challenging issues are thus to provide numerical methods for reconstructing the dielectric properties of a medium from a finite number of boundary measurements of the electric field. A classical way consists in formulating the inverse problem as the minimization of a cost function representing the difference between the measured and predicted fields. To solve the minimization problem, a gradient-based algorithm is currently used. For instance, Beilina et al. (e.g. [16,17]) have developed an adaptive finite element method based on a posteriori estimates for the simultaneous reconstruction of the real-valued electric permittivity and magnetic permeability functions of the 3D Maxwell's system. De Buhan and Darbas [18] have combined the quasi-Newton BFGS method and an iterative process (called the Adaptive Eigenspace Inversion) for determining the complex dielectric permittivity of a medium with 2D numerical validations. Another way to express the inverse medium problem is to search small anomalies in the electric parameters on a known background. We can cite the significant results of Ammari et al. (e.g. [19–21]) who have derived small-volume expansions of the electromagnetic field, the volume of the imperfections being the asymptotic parameter. This yields constructive numerical methods for the localization of small-volume electromagnetic defects from measurements on a part of the boundary (see e.g. [22] for 3D numerical results). This asymptotic approach has also been combined with an exact controllability method for retrieving small-amplitude perturbations in the permeability of a medium [23]. In this case, the time-dependent Maxwell equations and dynamic boundary measurements are considered. The performance of the reconstruction method in 2D has been addressed in [24]. In the present work, our aim is to detect and identify small-amplitude perturbations in the dielectric parameters of a medium from (time-independent) boundary measurements of the electric field. In this sense, our work falls into the previous class of approaches. We propose to investigate the problem from a different point of view. Our reconstruction method is based on explicit relations between sensitivity (with respect to the physical parameters) and characteristics of the imperfections. Sensitivity is the derivation of a given quantity (cost functional, physical field, etc.) with respect to the parameters. Sensitivity gives an interesting tool for understanding the impact of local changes in interior parameters on the observed boundary measurements at the surface of the studied object. For instance, sensitivity information has been recently used to answer concrete clinic questions in EEG (electroencephalography) for neonates [25] or to design resolution-based discretizations of the conductivity space in EIT [26]. In practical applications, the parameters are in general discretized, for instance by P0 or P1 finite elements. This leads to a Jacobian matrix which can be used to find the critical points of some least-square functional [27,28]. Another approach that is widely studied, is the topological sensitivity (or shape sensitivity) with respect to the shape of a perturbation in the parameters (see e.g. [29,30] for EIT). The topological sensitivity is also used for the detection and shape identification of scatterers (e.g [31–33] for acoustics and electromagnetism respectively). Data are in this case measurements of the far-field pattern of the scattered field.

In the present paper, the aim is to investigate the impact of small-amplitude perturbations in the dielectric parameters of a medium on boundary measurements of the electric field. The novelty lies in proposing a rigourous sensitivity analysis of the electric field with respect to the variations of the electric permittivity and conductivity, noticing that the perturbation in the measurements is proportional to sensitivity for small-amplitudes of the parameters. We address both theoretical and numerical aspects. The sensitivity analysis is the first step for developing a new non- iterative inversion algorithm that allows to determine the location and volume of small-amplitude anomalies from boundary measurements of the perturbed electric field. To our knowledge, it's the first time that such a sensitivity analysis with respect to parameters (and not to the shape) is realized for solving an inverse electromagnetic medium problem.

The remainder of the paper is organized as follows. In Section 2, we present the forward problem under consideration and define the functional setting. In Section 3, we propose a theoretical sensitivity analysis of the electric field with respect to the electric permittivity and/or conductivity of a medium which is illustrated by numerical simulations. Section 4 is devoted to the sensitivity analysis in the case of a constant background which allows to write the sensitivity boundary data as the solution of an integral equation. In Section 5, we explain how to use the previous results for solving an inverse medium problem. The localization procedure is described in Section 6, and various three-dimensional numerical results are reported to illustrate the method. Finally, we give some conclusions and perspectives in the last section.

2. The forward problem

2.1. Time-harmonic Maxwell's equations

Let Ω denote a bounded and simply connected domain in ${\mathbb{R}}^3$ of Lipschitz boundary $\mathrm{\Gamma }=\mathrm{\partial \Omega }$. The unit outward normal to Ω is denoted by $\mathbf{n}$. We introduce the vector spaces $\begin{array}{rl}H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)& =\left\{\mathbf{u}\in L^2(\mathrm{\Omega }{)}^3;\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{u}\in L^2(\mathrm{\Omega }{)}^3\right\},\\ Y\left(\mathrm{\Gamma }\right)& =\left\{\mathbf{f}\in H^{-1/2}(\mathrm{\Gamma }{)}^3;\mathrm{\exists }\mathbf{u}\in H(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l})|\mathbf{u}\times \mathbf{n}=\mathbf{f}\right\}.\end{array}$

We are interested in time-harmonic Maxwell's equations with Neumann boundary condition: $\left\{\begin{array}{rlrl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}-k^2\kappa \mathbf{E}=\mathbf{F}& & \mathrm{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}\times \mathbf{n}=\mathbf{g}& & \mathrm{on}\mathrm{\Gamma }.\end{array}\right.\mathcal{M}$ Here, $\mathbf{E}$ denotes the electric field intensity in Ω, and the fields $\mathbf{F}$ and $\mathbf{g}$ are given source terms in, respectively, $L^2(\mathrm{\Omega }{)}^3$ and $Y\left(\mathrm{\Gamma }\right)$. The number $k:=\omega \sqrt{{\mu }_0{\epsilon }_0}$ is the wavenumber with ω the wave angular frequency, ${\mu }_0$ and ${\epsilon }_0$ respectively the magnetic permeability and electric permittivity in vacuum. We assume that the magnetic permeability in Ω is equal to ${\mu }_0$. Let ϵ and σ denote, respectively, the electric permittivity and conductivity in Ω. The refractive index κ of the medium in Ω is defined by $\kappa \left(\mathbf{x}\right)=\frac{1}{{\epsilon }_0}\left(\epsilon \left(\mathbf{x}\right)+i\frac{\sigma \left(\mathbf{x}\right)}{\omega }\right),\mathbf{x}\in \mathrm{\Omega }.$

We assume that Ω is decomposed into P connected Lipschitz subdomains, denoted by ${\mathrm{\Omega }}_p$ for $1\le p\le P$, such that $\overline{\mathrm{\Omega }}=\bigcup _{p=1}^P{\overline{\mathrm{\Omega }}}_p\mathrm{and}{\mathrm{\Omega }}_p\cap {\mathrm{\Omega }}_q=\mathrm{\varnothing },\mathrm{\forall }p\ne q.$ Moreover, the following assumptions are made on the parameter κ: $\left\{\begin{array}{rl}& \mathrm{\forall }1\le p\le P,{\kappa ∣}_{{\mathrm{\Omega }}_p}\in H^3\left({\mathrm{\Omega }}_p\right),\\ & \mathrm{\exists }{\alpha }_R>0,\mathcal{R}\left(\kappa \right)\ge {\alpha }_R\text{in}\mathrm{\Omega },\\ & \mathrm{\exists }{\alpha }_I>0,\mathrm{\forall }1\le p\le P\text{either}\mathcal{I}\left({\kappa ∣}_{{\mathrm{\Omega }}_p}\right)\ge {\alpha }_I\text{or}\mathcal{I}\left({\kappa ∣}_{{\mathrm{\Omega }}_p}\right)=0,\\ & \mathcal{I}\left(\kappa \right)\not\equiv 0\text{in}\mathrm{\Omega }.\end{array}\right.{\mathcal{H}}_{\kappa }$

The variational formulation of ($\mathcal{M}$) is $\left\{\begin{array}{rl}& \mathrm{Find}\mathbf{E}\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E},\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa \mathbf{E},\mathit{\phi })=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right),\end{array}\right.{\mathcal{M}}_v$ where $(\cdot ,\cdot )$ denotes the dot-product in $L^2(\mathrm{\Omega }{)}^3$ and $⟨\cdot ,\cdot {⟩}_{\mathrm{\Gamma }}$ denotes the dual product from $Y\left(\mathrm{\Gamma }\right)$ to $Y(\mathrm{\Gamma }{)}^{\prime }$. Here, ${(\cdot )}_T$ denotes the extension to $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ of the map: $\mathbf{v}\mapsto \mathbf{n}\times ({\mathbf{v}∣}_{\mathrm{\Gamma }}\times \mathbf{n}),$ classically defined on $\mathcal{C}(\overline{\mathrm{\Omega }}{)}^3$ (see [34, Theorem 3.31]).

Theorem 2.1

Under the above assumptions, the problem (${\mathcal{M}}_v$) admits a unique solution $\mathbf{E}$ in $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ depending continuously on $\mathbf{F}$ and $\mathbf{g}$.

Sketch of the proof: The proof is adapted from [34]. The main ingredients are a judicious Helmholtz decomposition of $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ and the Fredholm alternative.

2.2. Regularity of the variational solution

In a setting where the subdomains ${\mathrm{\Omega }}_p$ are only Lipschitz, the solution of (${\mathcal{M}}_v$) may have poor regularity. Indeed, singularities are likely to occur at the corners and edges of the subdomains, and the solution $\mathbf{E}$ does not belong, in general, to $H^1({\mathrm{\Omega }}_p{)}^3$ (see [35]). In the present paper and with regard to the (biomedical) applications that we have in mind, we do not deal with these questions of singularities. Therefore, we assume from now on that Ω as well as all subdomains ${\mathrm{\Omega }}_p$ are at least of class ${\mathcal{C}}^{1,1}$. This assumption allows to prove regularity results for the solution of (${\mathcal{M}}_v$) and the sensitivity equation that will be presented hereafter.

According to the partition of Ω into P subdomains $({\mathrm{\Omega }}_p{)}_p$, we introduce the spaces of piecewise smooth functions: for s>0, let ${\mathrm{PH}}^s\left(\mathrm{\Omega }\right)=\left\{v\in L^2\left(\mathrm{\Omega }\right);{v∣}_{{\mathrm{\Omega }}_p}\in H^s\left({\mathrm{\Omega }}_p\right)\mathrm{\forall }1\le p\le P\right\}.$ We adopt the notation ${\mathrm{PH}}^s(\mathrm{\Omega }{)}^3$ to denote spaces of piecewise smooth vector fields. We further introduce the classical space $H\left(\mathrm{div}\right)=\left\{\mathbf{u}\in L^2(\mathrm{\Omega }{)}^3;\mathrm{div}\mathbf{u}\in L^2(\mathrm{\Omega })\right\}$ as well as the trace space $H^s\left({\mathrm{div}}_{\mathrm{\Gamma }}\right)=\left\{\mathbf{f}\in H_t^s\left(\mathrm{\Gamma }\right);{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{f}\in H^s\left(\mathrm{\Gamma }\right)\right\}$ where ${\mathrm{div}}_{\mathrm{\Gamma }}$ denotes the surface divergence operator defined on the subspace $H_t^s\left(\mathrm{\Gamma }\right)$ of $H^s(\mathrm{\Gamma }{)}^3$ of tangential fields $\mathbf{f}$. A rigorous definition of ${\mathrm{div}}_{\mathrm{\Gamma }}$ can be found in [34,36]. The following regularity result can be deduced from [35]:

Theorem 2.2

Let Ω as well as all subdomains ${\mathrm{\Omega }}_p$ be of class ${\mathcal{C}}^{1,1}$. Assume that the source term $\mathbf{F}$ belongs to $H\left(\mathrm{div}\right)$ and satisfies $\mathbf{F}\cdot \mathbf{n}\in H^{1/2}\left(\mathrm{\Gamma }\right)$. Assume further that $\mathbf{g}\in H^{1/2}\left({\mathrm{div}}_{\mathrm{\Gamma }}\right)$. Let κ be a piecewise constant function with respect to the partition of Ω that satisfies the assumptions of Theorem 2.1. Then, the solution of (${\mathcal{M}}_v$) belongs to ${\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$ and satisfies (1a) $\begin{array}{rl}& -k^2\mathrm{div}\left(\kappa \mathbf{E}\right)=\mathrm{div}\mathbf{F}\text{in}\mathrm{\Omega },\end{array}$(1a) (1b) $\begin{array}{rl}& k^2\kappa \mathbf{E}\cdot \mathbf{n}=-\mathbf{F}\cdot \mathbf{n}+{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{g}\text{on}\mathrm{\Gamma }.\end{array}$(1b)

Proof.

Let $\mathbf{E}\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ be the solution of (${\mathcal{M}}_v$). Since $\mathbf{E}$ satisfies ($\mathcal{M}$) in the distributional sense, we get (1a(1a) $\begin{array}{rl}& -k^2\mathrm{div}\left(\kappa \mathbf{E}\right)=\mathrm{div}\mathbf{F}\text{in}\mathrm{\Omega },\end{array}$(1a) ) where the right-hand side belongs to $L^2\left(\mathrm{\Omega }\right)$. We then deduce (1b(1b) $\begin{array}{rl}& k^2\kappa \mathbf{E}\cdot \mathbf{n}=-\mathbf{F}\cdot \mathbf{n}+{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{g}\text{on}\mathrm{\Gamma }.\end{array}$(1b) ) from Green's formula and (${\mathcal{M}}_v$). Now, let $p\in H^1\left(\mathrm{\Omega }\right)/\mathbb{R}$ be the unique solution of the Neumann problem $\begin{array}{rl}& -\mathrm{div}\left(\kappa \mathrm{\nabla }p\right)=\mathrm{div}\mathbf{F}\text{in}\mathrm{\Omega },\\ & \kappa {\mathrm{\partial }}_{n}p=-\mathbf{F}\cdot \mathbf{n}+{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{g}\text{on}\mathrm{\Gamma }.\end{array}$ According to the assumptions on the data $\mathbf{F}$ and $\mathbf{g}$ and due to the regularity of Ω and its subdomains, the scalar potential p belongs to ${\mathrm{PH}}^2\left(\mathrm{\Omega }\right)$. Then, let ${\mathbf{E}}_0=\mathbf{E}-\left(1/k^2\right)\mathrm{\nabla }p$. ${\mathbf{E}}_0$ obviously belongs to $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ and is divergence free in the sense that $\mathrm{div}\left(\kappa {\mathbf{E}}_0\right)=0$. It also satisfies the homogeneous boundary condition $\kappa {\mathbf{E}}_0\cdot \mathbf{n}=0$ on Γ. Therefore, we can apply Theorem 3.5 of [35], and deduce that the field ${\mathbf{E}}_0$ admits a decomposition ${\mathbf{E}}_0={\mathbf{E}}_{0,R}+\mathrm{\nabla }{p}_0$ where ${\mathbf{E}}_{0,R}\in {\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$ and $p_0\in H^1\left(\mathrm{\Omega }\right)/\mathbb{R}$ is the unique solution of a Neumann problem with the right-hand side in $L^2\left(\mathrm{\Omega }\right)$ and homogeneous boundary condition. Again, $p_0$ belongs to ${\mathrm{PH}}^2\left(\mathrm{\Omega }\right)$ in the present setting of regular subdomains. This shows that ${\mathbf{E}}_0$ belongs to ${\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$ and implies $\mathbf{E}\in {\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$ due to the regularity of the scalar potential p.

In the case of regular data and a constant parameter κ, a stronger regularity result can be obtained for $\mathbf{E}$:

Theorem 2.3

Let Ω be of class ${\mathcal{C}}^{2,1}$. Let κ be a constant such that $\mathcal{R}\left(\kappa \right)\ge 0$ and $\mathcal{I}\left(\kappa \right)\ge 0$. Assume that $\mathbf{F}\in H^1(\mathrm{\Omega }{)}^3$ with $\mathrm{div}\mathbf{F}\in H^1\left(\mathrm{\Omega }\right)$ and $\mathbf{F}\cdot \mathbf{n}\in H^{3/2}\left(\mathrm{\Gamma }\right)$. Assume further that $\mathbf{g}\in H^{3/2}\left({\mathrm{div}}_{\mathrm{\Gamma }}\right)$. Then, $\mathbf{E}\in H^2(\mathrm{\Omega }{)}^3$.

Proof.

The proof is based on a result from [37]: if Ω is of class ${\mathcal{C}}^{m,1}$ for $m\in \mathbb{N}$, the spaces $\left\{\mathbf{v}\in L^2(\mathrm{\Omega }{)}^3;\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{v}\in H^{m-1}(\mathrm{\Omega }{)}^3;\mathrm{div}\mathbf{v}\in H^{m-1}\left(\mathrm{\Omega }\right);\mathbf{v}\times \mathbf{n}\in H^{m-1/2}(\mathrm{\Gamma }{)}^3\right\}$ and $\left\{\mathbf{v}\in L^2(\mathrm{\Omega }{)}^3;\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{v}\in H^{m-1}(\mathrm{\Omega }{)}^3;\mathrm{div}\mathbf{v}\in H^{m-1}\left(\mathrm{\Omega }\right);\mathbf{v}\cdot \mathbf{n}\in H^{m-1/2}\left(\mathrm{\Gamma }\right)\right\}$ are both continuously imbedded in $H^m(\mathrm{\Omega }{)}^3$.

Now, let $\mathbf{w}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}$. According to ($\mathcal{M}$) and the regularity result of Theorem 2.2, we have $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{w}=k^2\kappa \mathbf{E}+F\in H^1(\mathrm{\Omega }{)}^3$ as well as $\mathrm{div}\mathbf{w}=0\text{in}\mathrm{\Omega }$ and $\mathbf{w}\times \mathbf{n}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}\times \mathbf{n}=\mathbf{g}\in H^{3/2}\left(\mathrm{\Gamma }\right)$. Thus, $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}=\mathbf{w}\in H^2(\mathrm{\Omega }{)}^3$. Furthermore, the regularity assumptions on $\mathbf{F}$ and $\mathbf{g}$ ensure that $\mathrm{div}\mathbf{E}\in H^1\left(\mathrm{\Omega }\right)$ and $\mathbf{E}\cdot \mathbf{n}\in H^{3/2}\left(\mathrm{\Gamma }\right)$. Therefore, $\mathbf{E}\in H^2(\mathrm{\Omega }{)}^3$.

3. Sensitivity analysis with respect to a perturbation of electric parameters

Sensitivity analysis determines how the solution of a problem varies when a slight perturbation is induced in some of its physical parameters. Here, we are interested in the sensitivity analysis of the electric field with respect to the electrical permittivity and/or conductivity. Mathematically, it may be described rigorously by the Gâteaux derivative (see for example [38]).

Definition 3.1

Let $F:X\to Y$ be an application between two Banach spaces X and Y. Let $U\subset X$ be an open set. The Gâteaux derivative of F at $\tau \in U$ in the direction $\varrho \in X$ is defined as $D_{\varrho }F\left(\tau \right)=\underset{h\to 0}{lim}\frac{F(\tau +h\varrho )-F\left(\tau \right)}{h}$ if the limit exists. If it exists for any direction $\varrho \in X$ and if the application $\varrho \mapsto D_{\varrho }F\left(\tau \right)$ is linear and continuous from X to Y, then we say that F is Gâteaux differentiable at τ.

3.1. Sensitivity equation

We define the space of parameters $\mathcal{P}=\left\{(\epsilon ,\sigma )\in L^{\mathrm{\infty }}(\mathrm{\Omega }{)}^2;\kappa \in {\mathrm{PH}}^3(\mathrm{\Omega })\right\}$ which is a Banach space, equipped with the norm ${∥(\epsilon ,\sigma )∥}_{\mathcal{P}}=max({∥\epsilon ∥}_{L^{\mathrm{\infty }}},{∥\sigma ∥}_{L^{\mathrm{\infty }}}),\mathrm{\forall }(\epsilon ,\sigma )\in \mathcal{P}.$ We define the open set of admissible parameters ${\mathcal{P}}_{\mathrm{adm}}=\left\{(\epsilon ,\sigma )\in \mathcal{P};{\epsilon }_{\min }<\epsilon <{\epsilon }_{\max }\mathrm{and}{\sigma }_{\min }<\sigma <{\sigma }_{\max }\mathrm{in}\mathrm{\Omega }\right\}$ where $0<{\epsilon }_{\min }<{\epsilon }_{\max }$ and $0<{\sigma }_{\min }<{\sigma }_{\max }$ are real constants. From Theorem 2.1, we deduce that, for any $\tau =(\epsilon ,\sigma )\in {\mathcal{P}}_{\mathrm{adm}}$, the problem (${\mathcal{M}}_v$) admits a unique solution, denoted by $\mathbf{E}(\cdot ,\tau )$.

Theorem 3.2

Let $\tau \in {\mathcal{P}}_{\mathrm{adm}}$. Let $h_0>0$ be such that $\tau +h\varrho \in {\mathcal{P}}_{\mathrm{adm}}$ for any $h\in [-h_0,h_0]$ and $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })\in \mathcal{P}$. Then $\mathbf{E}(\cdot ,\tau )$ is Gâteaux differentiable at τ in the direction ϱ. Moreover, its derivative $D_{\varrho }\mathbf{E}(\cdot ,\tau )$ is the unique solution of the following variational problem: $\begin{array}{rl}& \mathrm{Find}{\mathbf{E}}^1\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right)\mathbf{E},\mathit{\phi }\right),\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right).\mathcal{S}\end{array}$

To simplify the writing of the proof of this result, we introduce the following notation, for any couple of positive (or null) reals a and b: $a\lesssim b⟺(\mathrm{\exists }C>0|a\le Cb),$ where C is a constant independent of a and b.

Proof.

Let $h\in [-h_0,h_0]\setminus \left\{0\right\}$ and $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })\in \mathcal{P}$. Let ${\mathbf{E}}_h=\mathbf{E}(\cdot ,\tau +h\varrho )$.

The field $\mathbf{E}$ is the unique solution of (2) $\begin{array}{rl}& \mathrm{Find}\mathbf{E}\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E},\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-\frac{k^2}{{\epsilon }_0}\left(\left(\epsilon +i\frac{\sigma }{\omega }\right)\mathbf{E},\mathit{\phi }\right)\\ & \phantom{-}=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right),\end{array}$(2) whereas the field ${\mathbf{E}}_h$ is the unique solution of (3) $\begin{array}{rl}& \mathrm{Find}{\mathbf{E}}_h\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-\frac{k^2}{{\epsilon }_0}\left(\left((\epsilon +h{\varrho }_{\epsilon })+i\frac{\sigma +h{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right)\\ & \phantom{-}=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right).\end{array}$(3) Let $\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$. We compute the difference between (3(3) $\begin{array}{rl}& \mathrm{Find}{\mathbf{E}}_h\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-\frac{k^2}{{\epsilon }_0}\left(\left((\epsilon +h{\varrho }_{\epsilon })+i\frac{\sigma +h{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right)\\ & \phantom{-}=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right).\end{array}$(3) ) and (2(2) $\begin{array}{rl}& \mathrm{Find}\mathbf{E}\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E},\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-\frac{k^2}{{\epsilon }_0}\left(\left(\epsilon +i\frac{\sigma }{\omega }\right)\mathbf{E},\mathit{\phi }\right)\\ & \phantom{-}=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right),\end{array}$(2) ) and we divide by h to find (4) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_h^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right),$(4) where ${\mathbf{E}}_h^1=\left({\mathbf{E}}_h-\mathbf{E}\right)/h$.

We now compute the difference between (4(4) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_h^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right),$(4) ) and ($\mathcal{S}$) to obtain (5) $\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right({\mathbf{E}}_h^1-{\mathbf{E}}^1),\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2\left(\kappa \right({\mathbf{E}}_h^1-{\mathbf{E}}^1),\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right)({\mathbf{E}}_h-\mathbf{E}),\mathit{\phi }\right).$(5) We note ${\tilde{\mathbf{E}}}_h={\mathbf{E}}_h^1-{\mathbf{E}}^1$ and ${\tilde{\mathbf{F}}}_h=\left(k^2/{\epsilon }_0\right)({\varrho }_{\epsilon }+i{\varrho }_{\sigma }/\omega )({\mathbf{E}}_h-\mathbf{E})$. We get (6) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\tilde{\mathbf{E}}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\tilde{\mathbf{E}}}_h,\mathit{\phi })=({\tilde{\mathbf{F}}}_h,\mathit{\phi }).$(6)

As ${\varrho }_{\epsilon }$ and ${\varrho }_{\sigma }$ are in $L^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$, we have ${\tilde{\mathbf{F}}}_h\in L^2(\mathrm{\Omega }{)}^3$. Then (6(6) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\tilde{\mathbf{E}}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\tilde{\mathbf{E}}}_h,\mathit{\phi })=({\tilde{\mathbf{F}}}_h,\mathit{\phi }).$(6) ) can be seen as the variational formulation of Maxwell's equations with a homogeneous Neumann boundary condition and the source term ${\tilde{\mathbf{F}}}_h$. From Theorem 2.1, we deduce that ${\tilde{\mathbf{E}}}_h\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ is the unique field satisfying (6(6) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\tilde{\mathbf{E}}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\tilde{\mathbf{E}}}_h,\mathit{\phi })=({\tilde{\mathbf{F}}}_h,\mathit{\phi }).$(6) ) for all $\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$. Moreover, we know that ${∥{\tilde{\mathbf{E}}}_h∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}\lesssim {∥{\tilde{\mathbf{F}}}_h∥}_0\lesssim {∥{\mathbf{E}}_h-\mathbf{E}∥}_0\lesssim {∥{\mathbf{E}}_h-\mathbf{E}∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}.$ We now use the definition of ${\mathbf{E}}_h^1$ to get $|h{|}^{-1}{∥{\mathbf{E}}_h-\mathbf{E}∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}={∥{\mathbf{E}}_h^1∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}.$ Furthermore, ${\mathbf{E}}_h^1$ satisfies (4(4) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_h^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right),$(4) ) for all $\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ and we have ${∥{\mathbf{E}}_h^1∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}\lesssim {∥{\mathbf{E}}_h∥}_0\lesssim {∥{\mathbf{E}}_h∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}\lesssim {∥\mathbf{F}∥}_0+{∥\mathbf{g}∥}_{H^{-1/2}}$ since ${\mathbf{E}}_h$ is the unique solution of (3(3) $\begin{array}{rl}& \mathrm{Find}{\mathbf{E}}_h\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)\text{such that}\\ & (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_h,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-\frac{k^2}{{\epsilon }_0}\left(\left((\epsilon +h{\varrho }_{\epsilon })+i\frac{\sigma +h{\varrho }_{\sigma }}{\omega }\right){\mathbf{E}}_h,\mathit{\phi }\right)\\ & \phantom{-}=(\mathbf{F},\mathit{\phi })+⟨\mathbf{g},{\mathit{\phi }}_T{⟩}_{\mathrm{\Gamma }},\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right).\end{array}$(3) ). Combining these inequalities, we get ${∥{\tilde{\mathbf{E}}}_h∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}\lesssim |h|.$ Thus ${\mathbf{E}}_h^1$ converges to ${\mathbf{E}}^1$ in $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$.

In order to prove the linearity of the application $\varrho \mapsto D_{\varrho }\mathbf{E}(\cdot ,\tau )$, let $\varrho =\lambda {\varrho }_1+{\varrho }_2$ with $\lambda \in \mathbb{C}$ and ${\varrho }_j=({\varrho }_{j,1},{\varrho }_{j,2})\in \mathcal{P},\mathrm{\forall }j\in \{1,2\}$. For $j\in \{1,2\}$, we set ${\mathbf{E}}_j^1:=D_{{\varrho }_j}\mathbf{E}(\cdot ,\tau )$. Thus ${\mathbf{E}}_j^1$ solves $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_j^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_j^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{j,1}+i\frac{{\varrho }_{j,2}}{\omega }\right)\mathbf{E},\mathit{\phi }\right),\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right).$ Let ${\mathbf{E}}^1=\lambda {\mathbf{E}}_1^1+{\mathbf{E}}_2^1$. By linearity, we have $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}\left(\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right)\mathbf{E},\mathit{\phi }\right),\mathrm{\forall }\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right),$ where $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })=(\lambda {\varrho }_{1,1}+{\varrho }_{2,1},\lambda {\varrho }_{1,2}+{\varrho }_{2,2})$. Then ${\mathbf{E}}^1$ is solution of the problem satisfied by $D_{\varrho }\mathbf{E}(\cdot ,\tau )$. From the uniqueness of the solution, we deduce that $D_{\lambda {\varrho }_1+{\varrho }_2}\mathbf{E}(\cdot ,\tau )=\lambda D_{{\varrho }_1}\mathbf{E}(\cdot ,\tau )+D_{{\varrho }_2}\mathbf{E}(\cdot ,\tau ).$

We obtain that $D_{\varrho }\mathbf{E}(\cdot ,\tau )$ is solution of ($\mathcal{S}$). Moreover we have ${∥D_{\varrho }\mathbf{E}(\cdot ,\tau )∥}_{H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)}\lesssim {∥\left({\varrho }_{\epsilon }+i\frac{{\varrho }_{\sigma }}{\omega }\right)\mathbf{E}∥}_0\lesssim {∥\varrho ∥}_{\mathcal{P}}.$ Thus, the application $\varrho \mapsto D_{\varrho }\mathbf{E}(\cdot ,\tau )$ is linear and continuous from $\mathcal{P}$ to $H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$.

3.2. Regularity of the solution to the sensitivity equation

The derivative ${\mathbf{E}}^1=D_{\varrho }\mathbf{E}(\cdot ,\tau )$ of $\mathbf{E}$ with respect to the parameter $\tau =(\epsilon ,\sigma )$ in the direction $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })$ is solution of the following boundary value problem: (7) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-k^2\kappa {\mathbf{E}}^1=k^2\chi \mathbf{E},\text{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on}\mathrm{\Gamma }\end{array}$(7) where $\chi =\left(1/{\epsilon }_0\right)({\varrho }_{\epsilon }+i({\varrho }_{\sigma }/\omega \left)\right)$.

Theorem 3.3

Let $\chi \in W^{1,\mathrm{\infty }}\left(\mathrm{\Omega }\right)$. Under the assumptions of Theorem 2.2, the solution of ($\mathcal{S}$) belongs to ${\mathrm{PH}}^1\left(\mathrm{\Omega }\right)$ and satisfies (8) $\begin{array}{rl}& -\mathrm{div}\left(\kappa {\mathbf{E}}^1\right)=\mathrm{div}\left(\chi \mathbf{E}\right)\text{in}\mathrm{\Omega },\end{array}$(8) (9) $\begin{array}{rl}& \kappa {\mathbf{E}}^1\cdot \mathbf{n}=-\chi \mathbf{E}\cdot \mathbf{n}\text{on}\mathrm{\Gamma }.\end{array}$(9)

Proof.

Let ${\mathbf{E}}^1\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$ be the solution of ($\mathcal{S}$). The regularity assumption on χ implies that $\mathrm{div}\left(\chi \mathbf{E}\right)$ belongs to $L^2\left(\mathrm{\Omega }\right)$, and (8(8) $\begin{array}{rl}& -\mathrm{div}\left(\kappa {\mathbf{E}}^1\right)=\mathrm{div}\left(\chi \mathbf{E}\right)\text{in}\mathrm{\Omega },\end{array}$(8) ) follows immediately from (7(7) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-k^2\kappa {\mathbf{E}}^1=k^2\chi \mathbf{E},\text{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on}\mathrm{\Gamma }\end{array}$(7) ). The second identity (9(9) $\begin{array}{rl}& \kappa {\mathbf{E}}^1\cdot \mathbf{n}=-\chi \mathbf{E}\cdot \mathbf{n}\text{on}\mathrm{\Gamma }.\end{array}$(9) ) can be obtained as in Theorem 2.2. Since $\mathbf{E}$ belongs to ${\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$, its normal trace on Γ is an element of $H^{1/2}\left(\mathrm{\Gamma }\right)$. The same arguments as in Theorem 2.2 then yield ${\mathbf{E}}^1\in {\mathrm{PH}}^1(\mathrm{\Omega }{)}^3$.

As for the solution of (${\mathcal{M}}_v$), we get more regularity in the case of a constant function κ.

Theorem 3.4

Let $\chi \in W^{2,\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ and κ a constant. Under the assumptions of Theorem 2.3, the solution of ($\mathcal{S}$) belongs to $H^2(\mathrm{\Omega }{)}^3$.

Proof.

Under the given assumptions, the solution $\mathbf{E}$ of (${\mathcal{M}}_v$) belongs to $H^2(\mathrm{\Omega }{)}^3$ according to Theorem 2.3. Together with the regularity of χ, we thus infer from (8(8) $\begin{array}{rl}& -\mathrm{div}\left(\kappa {\mathbf{E}}^1\right)=\mathrm{div}\left(\chi \mathbf{E}\right)\text{in}\mathrm{\Omega },\end{array}$(8) ) and (9(9) $\begin{array}{rl}& \kappa {\mathbf{E}}^1\cdot \mathbf{n}=-\chi \mathbf{E}\cdot \mathbf{n}\text{on}\mathrm{\Gamma }.\end{array}$(9) ) that $\mathrm{div}{\mathbf{E}}^1\in H^1\left(\mathrm{\Omega }\right)$ and ${\mathbf{E}}^1\cdot \mathbf{n}\in H^{3/2}\left(\mathrm{\Gamma }\right)$. As in the proof of Theorem 2.3, $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}$ can be shown to belong to $H^2(\mathrm{\Omega }{)}^3$. The regularity result follows from [37].

3.3. Some properties of the sensitivity

In this section, we prove some properties of the sensitivity that are directly linked to the linearity of the Gâteaux derivative.

Proposition 3.5

Let $\tau \in {\mathcal{P}}_{\mathrm{adm}}$. Let $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })\in \mathcal{P}$. Then we have $D_{\varrho }\mathbf{E}(\cdot ,\tau )=D_{({\varrho }_{\epsilon },0)}\mathbf{E}(\cdot ,\tau )+D_{(0,{\varrho }_{\sigma })}\mathbf{E}(\cdot ,\tau ).$

Proposition 3.6

Let $\tau \in {\mathcal{P}}_{\mathrm{adm}}$. Let $\varrho \in {\mathrm{PH}}^3\left(\mathrm{\Omega }\right)$. We set ${\mathbf{E}}_{\epsilon }^1:=D_{(\varrho ,0)}\mathbf{E}(\cdot ,\tau )$ and ${\mathbf{E}}_{\sigma }^1:=D_{(0,\varrho )}\mathbf{E}(\cdot ,\tau )$. Then ${\mathbf{E}}_{\epsilon }^1=-i\omega {\mathbf{E}}_{\sigma }^1.$

Proof.

Let ${\mathbf{E}}^1={\mathbf{E}}_{\epsilon }^1+i\omega {\mathbf{E}}_{\sigma }^1$. The result will be proved if we show that ${\mathbf{E}}^1=0$. To this end, let $\mathit{\phi }\in H\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\right)$. We have (10) $\begin{array}{rl}& (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}^1,\mathit{\phi })=(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\epsilon }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\epsilon }^1,\mathit{\phi })\\ & \phantom{=}+i\omega \left((\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\sigma }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\sigma }^1,\mathit{\phi })\right).\end{array}$(10) We then apply Theorem 3.2 to ${\mathbf{E}}_{\epsilon }^1$ and ${\mathbf{E}}_{\sigma }^1$ to find (11) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\epsilon }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\epsilon }^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}(\varrho \mathbf{E},\mathit{\phi })$(11) and (12) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\sigma }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\sigma }^1,\mathit{\phi })=i\frac{k^2}{{\epsilon }_0\omega }(\varrho \mathbf{E},\mathit{\phi }).$(12) We now inject (11(11) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\epsilon }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\epsilon }^1,\mathit{\phi })=\frac{k^2}{{\epsilon }_0}(\varrho \mathbf{E},\mathit{\phi })$(11) ) and (12(12) $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\sigma }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\sigma }^1,\mathit{\phi })=i\frac{k^2}{{\epsilon }_0\omega }(\varrho \mathbf{E},\mathit{\phi }).$(12) ) in (10(10) $\begin{array}{rl}& (\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}^1,\mathit{\phi })=(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\epsilon }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\epsilon }^1,\mathit{\phi })\\ & \phantom{=}+i\omega \left((\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\sigma }^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2(\kappa {\mathbf{E}}_{\sigma }^1,\mathit{\phi })\right).\end{array}$(10) ) to find $(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1,\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathit{\phi })-k^2\left(\kappa {\mathbf{E}}^1\mathit{\phi }\right)=\left(\frac{k^2}{{\epsilon }_0}-\omega \frac{k^2}{{\epsilon }_0\omega }\right)(\varrho \mathbf{E},\mathit{\phi })=0.$ Then ${\mathbf{E}}^1$ is solution of (${\mathcal{M}}_v$) with $\mathbf{F}=0$ and $\mathbf{g}=0$. By the uniqueness of the solution, we find that ${\mathbf{E}}^1=0$.

Remark 3.1

For large frequencies ω, Proposition 3.6 thus implies that the derivatives of the electric field with respect to the parameters ϵ and σ are not of the same order whenever the directions in which the derivatives are taken have comparable norms of order $\mathcal{O}\left(1\right)$. This statement suggests to study sensitivity with respect to the permittivity in a direction of order $1/\omega$. We refer to Figure 3 for an illustration.

We are interested in studying how the location of a perturbation affects the electric field. Therefore, we focus on derivatives in the direction of characteristic functions of the perturbations' supports. The numerical results of Section 3.4 show that in this case the sensitivity is localized and illustrate the following proposition.

Proposition 3.7

Let $(P_j{)}_{1\le j\le N}$ be a collection of $N\in {\mathbb{N}}^{\ast }$ subsets of Ω such that $P_{j_1}\cap P_{j_2}=\mathrm{\varnothing },\mathrm{\forall }{j}_1\ne j_2.$ For all $1\le j\le N,$ we denote by ${\varrho }_j$ the indicator function of $P_j$. Let ϱ be the indicator function of $\bigcup _{j=1}^N{P}_j$. Then we have $D_{(\varrho ,0)}\mathbf{E}(\cdot ,\tau )=\sum _{j=1}^N{D}_{({\varrho }_j,0)}\mathbf{E}(\cdot ,\tau )\mathrm{and}{D}_{(0,\varrho )}\mathbf{E}(\cdot ,\tau )=\sum _{j=1}^N{D}_{(0,{\varrho }_j)}\mathbf{E}(\cdot ,\tau ),$ for any $\tau \in {\mathcal{P}}_{\mathrm{adm}}$.

3.4. Numerical results and comments

We implemented the numerical solver for 3D Maxwell's equations with FreeFem++ (see [39]). Our test domain Ω is the unit ball of ${\mathbb{R}}^3$. We consider a tetrahedral mesh ${\mathcal{T}}_h$. For any $T\in {\mathcal{T}}_h$, let $h_T$ be its diameter. Then $h=\underset{T\in {\mathcal{T}}_h}{max}{h}_T$ is the mesh parameter of ${\mathcal{T}}_h$. For any h, we denote by $N_e$ the number of edges. Edge finite elements of order 1 (see [34,40]) are used to approximate the respective solutions of the problem (${\mathcal{M}}_v$) and of the sensitivity Equation ($\mathcal{S}$).

We consider that Ω is filled with a homogeneous medium of constant electrical permittivity $\epsilon ={10}^{-8}\mathrm{F}\cdot {\mathrm{m}}^{-1}$ and conductivity $\sigma =0.33\mathrm{S}\cdot {\mathrm{m}}^{-1}$ at the fixed frequency $\omega ={10}^6\mathrm{Hz}$. The mesh characteristics are h=0.12 and $N_e=167402$. The sensitivity ${\mathbf{E}}^1$ of the electric field in a given direction $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })$ is computed as the solution of Equation ($\mathcal{S}$).

First, we compare the modulus of the sensitivity with respect to a perturbation either of the conductivity or the permittivity (see Figure 1, left and right). This perturbation is modelled by a sphere $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ of radius $\alpha =0.1$, centred at ${\mathbf{x}}_0=(-0.8,0,0)$. The respective directions are $\varrho =({\mathbf{1}}_B/\omega ,0)$ for the permittivity and $\varrho =(0,{\mathbf{1}}_B)$ for the conductivity. This result illustrates Proposition 3.6 which implies $|{\mathbf{E}}_{\epsilon }^1|=|{\mathbf{E}}_{\sigma }^1|$. In the sequel, we consider a perturbation of the conductivity only. In the bottom of Figure 1, the perturbation is placed at a different position. The simulation indicates how the position of the inhomogeneity affects sensitivity. In particular, it shows that the sensitivity is localized to a surface area close to the inhomogeneity.

Figure 1. Sensitivity of the electric field on the surface with respect to the parameter (top left: inhomogeneity in the conductivity, top right: inhomogeneity in the permittivity) and the position of the inhomogeneity (top: centred at ${\mathbf{x}}_0=(-0.8,0,0)$ on the x-axis, bottom: centred at ${\mathbf{x}}_0=(0,-0.8,0)$ on the y-axis).

In Figure 2, we report the sensitivity corresponding to a spherical inhomogeneity centred at ${\mathbf{x}}_0=(-0.55,0,0)$ for different volumes. Compared to the perturbation at ${\mathbf{x}}_0=(-0.8,0,0)$ (see Figure 1, left), we observe that a deeper inhomogeneity leads to more spreaded surfacic perturbations. Moreover, increasing the inhomogeneity's size does not change the shape, but increases the amplitude of the sensitivity (see Figure 2, right).

Figure 2. Sensitivity of the electric field on the surface with respect to volume of the inhomogeneity (left: $\alpha =0.1$, right: $\alpha =0.3$).

Figure 3. Sensitivity of the electric field on the surface corresponding to two disjoint inhomogeneities.

In Figure 3, we present the sensitivity corresponding to two spherical inhomogeneities: one centred at ${\mathbf{x}}_0=(-0.85,0,0)$ of radius $\alpha =0.1$ and the other centred at $(0,-0.7,0)$ of radius $\alpha =0.2$. We retrieve two surfacic perturbations corresponding to each inhomogeneity in agreement with Proposition 3.7.

The numerical results of Figures 1–3 emphasize that sensitivity analysis provides information about those surface areas on which the electrical field is affected by small variations in the electric parameters of the medium. More precisely, the values of the sensitivity ${\mathbf{E}}^1$ give insights about the inhomogeneities' locations and sizes. We will see in Section 5 that it is a useful tool for solving the inverse problem of reconstructing the support of a perturbation in the permittivity and/or conductivity from boundary data. Indeed, the solution of the sensitivity Equation ($\mathcal{S}$) is linked to the perturbed electric field in the following way. Let $\tau :=(\epsilon ,\sigma )\in {\mathcal{P}}_{\mathrm{adm}}$ and $\varrho =({\varrho }_{\epsilon },{\varrho }_{\sigma })\in \mathcal{P}$, both fixed. As suggested in [41], a first order Taylor expansion of the solution ${\mathbf{E}}_p:=\mathbf{E}(\cdot ,\tau +h\varrho )$ of the perturbed problem with parameters $\tau +h\varrho$ for small-amplitudes of order h, $0<h\ll 1$, yields (13) $\mathbf{E}(\cdot ,\tau +h\varrho )-\mathbf{E}(\cdot ,\tau )\approx h{D}_{\varrho }\mathbf{E}(\cdot ,\tau ).$(13) In other words, for small values of h, the boundary data (measurements) $({\mathbf{E}}_p-\mathbf{E})\times \mathbf{n}$ have the same behaviour as the Gâteaux derivative of the electric field $\mathbf{E}$ in the direction ϱ.

4. Sensitivity analysis in the case of a constant background

In this section, we identify the tangential trace of the Gâteaux derivative ${\mathbf{E}}^1$ as the solution of a boundary integral equation. Estimates of the right-hand side of this equation exhibit some relations between the sensitivity ${\mathbf{E}}^1$ and geometric characteristics of the perturbation.

In the sequel, we assume that the material parameters ϵ and σ of the unperturbed background medium are positive constants. In order to simplify the notations, we introduce the complex-valued wavenumber ξ which is defined by (14) ${\xi }^2:=k^2\kappa ={\omega }^2{\mu }_0\left(\epsilon +i\frac{\sigma }{\omega }\right),$(14) where $\epsilon >0$ and $\sigma >0$.

Let $\mathbf{E}$ be the solution of (${\mathcal{M}}_v$) associated with ξ. Let $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ be the sphere of radius $\alpha >0$ and centre ${\mathbf{x}}_0\in \mathrm{\Omega }$. We assume that the distance between B and the boundary Γ is at least equal to a given value $\beta >0$ and we choose a neighbourhood ${\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$ of Γ such that ${\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)\cap B=\mathrm{\varnothing }$.

The perturbation occurs in the domain B and will be described by a function $f:\mathbf{x}\mapsto {\omega }^2{\mu }_0\left({\varrho }_{\epsilon }\left(\mathbf{x}\right)+i\frac{{\varrho }_{\sigma }\left(\mathbf{x}\right)}{\omega }\right),$ where $({\varrho }_{\epsilon },{\varrho }_{\sigma })\in \mathcal{P}$. We assume that f is regular, $f\in W^{2,\mathrm{\infty }}\left(\mathrm{\Omega }\right)$, and that $\mathrm{supp}\left(f\right)\subset B$.

According to Section 3 (see (7(7) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-k^2\kappa {\mathbf{E}}^1=k^2\chi \mathbf{E},\text{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on}\mathrm{\Gamma }\end{array}$(7) )), the Gâteaux derivative of $\mathbf{E}$ with respect to the (constant) parameters $\tau =(\epsilon ,\sigma )$ is solution of the boundary value problem, (15) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-{\xi }^2{\mathbf{E}}^1=f\mathbf{E},\text{in }\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on }\mathrm{\Gamma }.\end{array}$(15)

4.1. An integral equation

In order to state the integral equation for the tangential trace $\mathbf{a}={\mathbf{E}}^1\times \mathbf{n}$ on Γ, we introduce in the sequel suitable integral operators. Let Φ denote the fundamental solution in ${\mathbb{R}}^3$ of the Helmholtz equation with complex wavenumber ξ, $-\left(\mathrm{\Delta \Phi }+{\xi }^2\mathrm{\Phi }\right)={\delta }_0,$ satisfying the outgoing Sommerfeld condition as $|\mathbf{x}|\to \mathrm{\infty }$. Function $\mathrm{\Phi }\left(\mathbf{x}\right)$ is given by $\mathrm{\Phi }:\mathbf{x}\mapsto \frac{1}{4\pi }\frac{{\mathrm{e}}^{i\xi |\mathbf{x}|}}{|\mathbf{x}|},\mathbf{x}\ne 0.$ Define the space of continuous tangential fields on Γ, $\mathcal{T}\left(\mathrm{\Gamma }\right)=\left\{\mathbf{a}\in {\mathcal{C}}^0(\mathrm{\Gamma }{)}^3;\mathbf{a}\cdot \mathbf{n}=0\right\}.$ For $\mathbf{z}\in {\mathbb{R}}^3\setminus \mathrm{\Gamma }$, the vector potential $\mathbf{A}\left(\mathbf{z}\right)$ with density $\mathbf{a}\in \mathcal{T}\left(\mathrm{\Gamma }\right)$ is defined by (16) $\mathbf{A}\left(\mathbf{z}\right)={\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}.$(16) For the bounded domain $\mathrm{\Omega }\subset {\mathbb{R}}^3$, we denote by ${\mathbf{A}}^{-}$ the restriction of $\mathbf{A}$ to Ω. Similarly, ${\mathbf{A}}^{+}$ is the restriction of $\mathbf{A}$ to the exterior of Ω, ${\mathbf{A}}^{+}={\mathbf{A}}_{{\mathbb{R}}^3\setminus \overline{\mathrm{\Omega }}}$. The following theorem from [36] describes the behaviour of $\mathbf{A}$ on the boundary Γ.

Theorem 4.1

Assume that Ω is a domain of class ${\mathcal{C}}^2$ and let $\mathbf{a}\in \mathcal{T}\left(\mathrm{\Gamma }\right)$. Then, $\mathbf{A}$ is continuous across Γ, i.e. (17) $\mathbf{A}\left(\mathbf{z}\right)={\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\mathrm{d}s\left(\mathbf{x}\right),\mathrm{\forall }\mathbf{z}\in {\mathbb{R}}^3.$(17) Furthermore, the Neumann trace satisfies the jump condition (18) $\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma }:\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{\pm }\left(\mathbf{z}\right)\times \mathbf{n}={\int }_{\mathrm{\Gamma }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_z\left(\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\right)\times {\mathbf{n}}_{\mathbf{z}}\mathrm{d}s\left(\mathbf{x}\right)\mp \frac{1}{2}\mathbf{a}\left(\mathbf{z}\right),$(18) and the following relation holds true uniformly for all $\mathbf{z}\in \mathrm{\Gamma }:$ (19) $\underset{h\to 0^{+}}{lim}\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{+}(\mathbf{z}+h{\mathbf{n}}_z)-\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{-}(\mathbf{z}-h{\mathbf{n}}_z)\right)\times {\mathbf{n}}_{\mathbf{z}}=0.$(19)

We next introduce the magnetic dipole operator $\mathcal{M}$ which is defined for $\mathbf{a}\in \mathcal{T}\left(\mathrm{\Gamma }\right)$ by (20) $\mathcal{M}\mathbf{a}\left(\mathbf{z}\right)=2{\int }_{\mathrm{\Gamma }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{z}}\left(\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\right)\times {\mathbf{n}}_{\mathbf{z}}\mathrm{d}s\left(\mathbf{x}\right),\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma }.$(20) Finally, we introduce the fundamental solution of the Maxwell equations that can be derived from Φ in the following way: (21) $G:\mathbf{x}\mapsto -\mathrm{\Phi }\left(\mathbf{x}\right)\mathbb{I}+\frac{1}{{\xi }^2}{D}^2\mathrm{\Phi }\left(\mathbf{x}\right).$(21) Here, $\mathbb{I}\in {\mathcal{M}}_3\left(\mathbb{R}\right)$ is the identity matrix, and $D^2\mathrm{\Phi }\left(\mathbf{x}\right)$ denotes the Hessian of Φ at the point $\mathbf{x}$. G can be shown to solve the following equation in ${\mathbb{R}}^3$: (22) $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}G-{\xi }^{2}G=-{\delta }_0\mathbb{I},$(22) where the curl of the matrix valued function G has to be understood column wise [42]. Since Φ satisfies the outgoing radiation condition, G satisfies the following Silver-Müller condition [36]: $\underset{|\mathbf{x}|\to \mathrm{\infty }}{lim}|\mathbf{x}|\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}G\times \frac{\mathbf{x}}{|\mathbf{x}|}-i\xi G\right)=0.$ Then, we are able to state the following theorem.

Theorem 4.2

Let ${\mathbf{E}}^1$ be the solution of (15(15) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-{\xi }^2{\mathbf{E}}^1=f\mathbf{E},\text{in }\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on }\mathrm{\Gamma }.\end{array}$(15) ) for ξ given by (14(14) ${\xi }^2:=k^2\kappa ={\omega }^2{\mu }_0\left(\epsilon +i\frac{\sigma }{\omega }\right),$(14) ) with $\epsilon >0$ and $\sigma >0$. For $\mathbf{z}\in \mathrm{\Gamma }$, define $\mathbf{T}\left(\mathbf{z}\right)$ by (23) $\mathbf{T}\left(\mathbf{z}\right)=-2\left({\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z})f\left(\mathbf{x}\right)\mathbf{E}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}\right)\times \mathbf{n}.$(23) Under the regularity assumptions of Theorem 3.4, the tangential trace $\mathbf{a}={\mathbf{E}}^1\times \mathbf{n}$ is solution of the following integral equation on Γ: (24) $(\mathcal{I}-\mathcal{M})\mathbf{a}=\mathbf{T},$(24) where $\mathcal{I}$ denotes the identity operator.

The proof of Theorem 4.2 is adapted from [20] where an asymptotic expansion of the perturbed field is obtained in terms of the (small) radius of the perturbation. Notice that in the present study of sensitivity, the analysis is simplified and no asymptotic parameter occurs. In other words, we do not assume that the perturbation is small in size, but only in amplitude in order to connect the derivative to the perturbed field (13(13) $\mathbf{E}(\cdot ,\tau +h\varrho )-\mathbf{E}(\cdot ,\tau )\approx h{D}_{\varrho }\mathbf{E}(\cdot ,\tau ).$(13) ).

Proof.

Let $\mathbf{z}\in {\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$. According to Theorem 3.4, the solution of the sensitivity equation ($\mathcal{S}$) belongs to $H^2(\mathrm{\Omega }{)}^3\hookrightarrow {\mathcal{C}}^0(\overline{\mathrm{\Omega }}{)}^3$, and thus the duality product $<{\delta }_0\mathbb{I},{\mathbf{E}}^1>={\mathbf{E}}^1\left(\mathbf{z}\right)$ is well defined. From (22(22) $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}G-{\xi }^{2}G=-{\delta }_0\mathbb{I},$(22) ), we get (25) $-{\mathbf{E}}^1\left(\mathbf{z}\right)={\int }_{\mathrm{\Omega }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{x}}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{x}}G(\mathbf{x}-\mathbf{z}){\mathbf{E}}^1\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}-{\xi }^2{\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z}){\mathbf{E}}^1\left(\mathbf{x}\right)\mathrm{d}\mathbf{x},$(25) where the integrals have been to understood as duality products in $H^s\left(\mathrm{\Omega }\right)$ for appropriate values of s.

The following partial integration formula holds true for matrix valued functions $A:\mathrm{\Omega }\to {\mathbb{C}}^{3\times 3}$ and $B:\mathrm{\Omega }\to {\mathbb{C}}^{3\times p},p\in {\mathbb{N}}^{\ast }$: (26) ${\int }_{\mathrm{\Omega }}{}^{t}A\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}B-{\int }_{\mathrm{\Omega }}{}^t\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}A\right)B={\int }_{\mathrm{\Gamma }}{}^t(A\times \mathbf{n})B,$(26) where the vector product $A\times \mathbf{n}$ is taken column wise. Applying (26(26) ${\int }_{\mathrm{\Omega }}{}^{t}A\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}B-{\int }_{\mathrm{\Omega }}{}^t\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}A\right)B={\int }_{\mathrm{\Gamma }}{}^t(A\times \mathbf{n})B,$(26) ) twice to the first term on the right-hand side of (25(25) $-{\mathbf{E}}^1\left(\mathbf{z}\right)={\int }_{\mathrm{\Omega }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{x}}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{x}}G(\mathbf{x}-\mathbf{z}){\mathbf{E}}^1\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}-{\xi }^2{\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z}){\mathbf{E}}^1\left(\mathbf{x}\right)\mathrm{d}\mathbf{x},$(25) ) yields ${\mathbf{E}}^1\left(\mathbf{z}\right)={\int }_{\mathrm{\Gamma }}{}^t\left({\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{x}}G\right(\mathbf{x}-\mathbf{z})\times \mathbf{n}){\mathbf{E}}^1\left(\mathbf{x}\right)\mathrm{d}s\left(\mathbf{x}\right)-{\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z})f\left(\mathbf{x}\right)\mathbf{E}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}$ taking into account the symmetry of either G and $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}G$ as well as the strong formulation of the sensitivity equation (15(15) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1-{\xi }^2{\mathbf{E}}^1=f\mathbf{E},\text{in }\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}^1\times \mathbf{n}=0,\text{on }\mathrm{\Gamma }.\end{array}$(15) ). The remaining boundary integral on the right-hand side can be written in terms of the tangential trace of ${\mathbf{E}}^1$ taking into account that $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\left(D^2\mathrm{\Phi }\right)=0$ in the definition of G. Finally, the following identity holds true for any $\mathbf{z}$ in the neighbourhood ${\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$ of the boundary: (27) ${\mathbf{E}}^1\left(\mathbf{z}\right)-{\int }_{\mathrm{\Gamma }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{z}}\left(\mathrm{\Phi }\right(\mathbf{x}-\mathbf{z}\left)\right({\mathbf{E}}^1\left(\mathbf{x}\right)\times \mathbf{n}\left)\right)\mathrm{d}s\left(\mathbf{x}\right)=-{\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z})f\left(\mathbf{x}\right)\mathbf{E}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}.$(27) Notice that the volume integral on Ω is well defined since $\mathrm{supp}\left(f\right)\subset B$ and $\mathbf{z}\in {\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$. Hence, $\mathbf{x}\ne \mathbf{z}$ for any $\mathbf{x}\in B$ and G is regular on the integration domain.

Now, for $\mathbf{z}\in {\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$, let ${\mathbf{z}}_{\mathrm{\Gamma }}$ denote the projection of $\mathbf{z}$ onto Γ. Notice that ${\mathbf{z}}_{\mathrm{\Gamma }}$ is well defined since Γ is regular and ${\mathcal{V}}_{\beta }\left(\mathrm{\Gamma }\right)$ can be assumed to be sufficiently small. Taking the vector product in the identity (27(27) ${\mathbf{E}}^1\left(\mathbf{z}\right)-{\int }_{\mathrm{\Gamma }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_{\mathbf{z}}\left(\mathrm{\Phi }\right(\mathbf{x}-\mathbf{z}\left)\right({\mathbf{E}}^1\left(\mathbf{x}\right)\times \mathbf{n}\left)\right)\mathrm{d}s\left(\mathbf{x}\right)=-{\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z})f\left(\mathbf{x}\right)\mathbf{E}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}.$(27) ) with the vector ${\mathbf{n}}_{\mathbf{z}}={\mathbf{n}}_{{\mathbf{z}}_{\mathrm{\Gamma }}}$ and passing to the limit as $\mathbf{z}\to {\mathbf{z}}_{\mathrm{\Gamma }}\in \mathrm{\Gamma }$ yields the integral equation (24(24) $(\mathcal{I}-\mathcal{M})\mathbf{a}=\mathbf{T},$(24) ) on Γ according to (18(18) $\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma }:\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{\pm }\left(\mathbf{z}\right)\times \mathbf{n}={\int }_{\mathrm{\Gamma }}{\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}}_z\left(\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\right)\times {\mathbf{n}}_{\mathbf{z}}\mathrm{d}s\left(\mathbf{x}\right)\mp \frac{1}{2}\mathbf{a}\left(\mathbf{z}\right),$(18) ).

In order to obtain an estimate of ${\mathbf{E}}^1\times \mathbf{n}$, we analyse the operator involved in the integral equation (24(24) $(\mathcal{I}-\mathcal{M})\mathbf{a}=\mathbf{T},$(24) ). We introduce the following normed spaces of tangential fields with surface divergence: ${\mathcal{T}}_d\left(\mathrm{\Gamma }\right)=\left\{\mathbf{a}\in \mathcal{T}\left(\mathrm{\Gamma }\right);{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{a}\in {\mathcal{C}}^0\left(\mathrm{\Gamma }\right)\right\}$ and ${\mathcal{T}}_d^{0,\alpha }\left(\mathrm{\Gamma }\right)=\left\{\mathbf{a}\in \mathcal{T}\left(\mathrm{\Gamma }\right);\mathbf{a}\in {\mathcal{C}}^{0,\alpha }\left(\mathrm{\Gamma }\right);{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{a}\in {\mathcal{C}}^{0,\alpha }\left(\mathrm{\Gamma }\right)\right\}$ equipped with the respective graph norms.

Theorem 4.3

Under the assumptions of Theorem 4.2, the operator $\mathcal{I}-\mathcal{M}$ is bijective on the space ${\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$ and has a bounded inverse.

Proof.

According to [36, Theorems 6.15 and 6.16], we can state that the operator $\mathcal{M}:{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)\to {\mathcal{T}}_d^{0,\alpha }\left(\mathrm{\Gamma }\right)$ is continuous, whereas ${\mathcal{T}}_d^{0,\alpha }\left(\mathrm{\Gamma }\right)$ is compactly embedded in ${\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$. Hence, $\mathcal{M}$ is a compact operator on ${\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$.

We prove that $\mathcal{I}-\mathcal{M}$ is injective. To this end, let $\mathbf{a}\in {\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$ and define the vector field ${\mathbf{E}}_{\mathbf{a}}$ for $\mathbf{z}\in {\mathbb{R}}^3\setminus \mathrm{\Gamma }$ by ${\mathbf{E}}_{\mathbf{a}}\left(\mathbf{z}\right)=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\mathrm{d}s\left(\mathbf{x}\right)=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{A}\left(\mathbf{z}\right).$ Here $\mathbf{A}$ is the vector potential with density $\mathbf{a}$ introduced in (17(17) $\mathbf{A}\left(\mathbf{z}\right)={\int }_{\mathrm{\Gamma }}\mathrm{\Phi }(\mathbf{x}-\mathbf{z})\mathbf{a}\left(\mathbf{x}\right)\mathrm{d}s\left(\mathbf{x}\right),\mathrm{\forall }\mathbf{z}\in {\mathbb{R}}^3.$(17) ). We have $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathbf{a}}^{\pm }-{\xi }^2{\mathbf{E}}_{\mathbf{a}}^{\pm }=0\text{in}{\mathrm{\Omega }}^{\pm }$ as well as (28) ${\mathbf{E}}_{\mathbf{a}}^{\pm }\times \mathbf{n}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{\pm }\times \mathbf{n}=\frac{1}{2}\left(\mathcal{M}\mp \mathcal{I}\right)\mathbf{a},\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma }.$(28) Now, let $\mathbf{a}\in {\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$ such that $(\mathcal{I}-\mathcal{M})\mathbf{a}=0$. On the exterior domain ${\mathrm{\Omega }}^{+}={\mathbb{R}}^3\setminus \overline{\mathrm{\Omega }}$, ${\mathbf{E}}_{\mathbf{a}}^{+}$ is solution of the exterior Maxwell problem with homogeneous Dirichlet boundary condition. In addition, ${\mathbf{E}}^{+}$ satisfies the outgoing radiation condition at infinity as does the fundamental solution Φ. Consequently, ${\mathbf{E}}_{\mathbf{a}}^{+}\equiv 0$ on ${\mathrm{\Omega }}^{+}$ due to the uniqueness of the solution to the exterior Maxwell problem.

Next, applying identity (19(19) $\underset{h\to 0^{+}}{lim}\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{+}(\mathbf{z}+h{\mathbf{n}}_z)-\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{-}(\mathbf{z}-h{\mathbf{n}}_z)\right)\times {\mathbf{n}}_{\mathbf{z}}=0.$(19) ) to $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathbf{a}}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{A}$, we get $\underset{h\to 0^{+}}{lim}\left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathbf{a}}^{+}(\mathbf{z}+h{\mathbf{n}}_{\mathbf{z}})-\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathbf{a}}^{-}(\mathbf{z}-h{\mathbf{n}}_{\mathbf{z}})\right)\times {\mathbf{n}}_{\mathbf{z}}=0.$ But ${\mathbf{E}}_{\mathbf{a}}^{+}$ vanishes on ${\mathrm{\Omega }}^{+}$ and therefore, $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathbf{a}}^{-}\times \mathbf{n}=0\text{on}\mathrm{\Gamma }.$ The field ${\mathbf{E}}_{\mathbf{a}}^{-}$ is thus solution of the interior Maxwell problem with homogeneous Neumann boundary condition. According to the properties of the constant parameters $\epsilon >0$ and $\sigma >0$ in the definition of the wave number ξ, the only solution to the interior Maxwell problem is ${\mathbf{E}}_{\mathbf{a}}^{-}\equiv 0$ (see Theorem 2.1).

From the homogeneous integral equation $(\mathcal{I}-\mathcal{M})\mathbf{a}=0$, we deduce $\mathcal{M}\mathbf{a}=\mathbf{a}$. Together with the identity (28(28) ${\mathbf{E}}_{\mathbf{a}}^{\pm }\times \mathbf{n}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{A}}^{\pm }\times \mathbf{n}=\frac{1}{2}\left(\mathcal{M}\mp \mathcal{I}\right)\mathbf{a},\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma }.$(28) ), this yields ${\mathbf{E}}_{\mathbf{a}}^{-}\times \mathbf{n}=\frac{1}{2}(\mathcal{M}+\mathcal{I})\mathbf{a}=\mathbf{a}$ and thus $\mathbf{a}=0$ on Γ. This proves that the operator $\mathcal{I}-\mathcal{M}$ is injective. Since $\mathcal{M}$ has been shown to be compact on ${\mathcal{T}}_d\left(\mathrm{\Gamma }\right)$, we deduce from the Fredholm alternative that $\mathcal{I}-\mathcal{M}$ is bijective. Finally, $\mathcal{I}-\mathcal{M}$ has a bounded inverse according to the inverse (or open) mapping theorem.

Theorems 4.2 and 4.3 imply that the tangential trace $\mathbf{a}={\mathbf{E}}^1\times \mathbf{n}$, solution to the integral equation $(\mathcal{I}-\mathcal{M})\mathbf{a}=\mathbf{T},$ can be estimated by (29) $\parallel \mathbf{a}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}\le C_{\mathcal{M}}\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}.$(29)

Remark 4.1

The norm of the inverse operator in (29(29) $\parallel \mathbf{a}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}\le C_{\mathcal{M}}\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}.$(29) ) actually depends on the wavenumber ξ. This may be seen from a thorough analysis of the eigenvalues of the operator $\mathcal{M}$ when Ω is a sphere. In this case, an exact analytical expression of the eigenvalues can be obtained in function of Ricatti–Bessel and Ricatti–Hankel functions [43]. This study allows us to numerically observe the behaviour of the spectrum of the operator $(\mathcal{I}-\mathcal{M})$ and its inverse with respect to the wavenumber ξ. The operator $(\mathcal{I}-\mathcal{M})$ is called the MFIE (Magnetic Field Integral Equation) operator [42].

4.2. Estimates for $\mathbf{T}$

In this section, we will prove some estimates of the functional $\mathbf{T}$ on the boundary that are at the origin of the localization algorithm described in Section 6. The proof is mainly based on the following estimates of the fundamental solution G.

Lemma 4.4

Let the wavenumber $\xi \in \mathbb{C}$ be defined as in (14(14) ${\xi }^2:=k^2\kappa ={\omega }^2{\mu }_0\left(\epsilon +i\frac{\sigma }{\omega }\right),$(14) ) and let G be the associated fundamental solution of the Maxwell equations as in (21(21) $G:\mathbf{x}\mapsto -\mathrm{\Phi }\left(\mathbf{x}\right)\mathbb{I}+\frac{1}{{\xi }^2}{D}^2\mathrm{\Phi }\left(\mathbf{x}\right).$(21) ). There is a polynomial p of degree 3 with positive coefficients depending on $|\xi |$ satisfying $p\left(0\right)=0$ such that (30) $\left|G_{j\ell }\left(\mathbf{x}\right)\right|\le p\left(\frac{1}{|\mathbf{x}|}\right)\mathrm{\forall }1\le j,\ell \le 3,\mathrm{\forall }\mathbf{x}\ne 0.$(30) Similarly, there is a polynomial $q\in {\mathbb{P}}_4\left(\mathbb{R}\right)$ of degree 4 with positive coefficients depending on $|\xi |$ satisfying $q\left(0\right)=0,$ such that (31) $\left|{\mathrm{\partial }}_m{G}_{j\ell }\left(\mathbf{x}\right)\right|\le q\left(\frac{1}{|\mathbf{x}|}\right)\mathrm{\forall }1\le j,\ell ,m\le 3.$(31)

Proof.

Let $1\le j,\ell \le 3$ and $\mathbf{x}\ne 0$. We have $G_{j\ell }\left(\mathbf{x}\right)=-\frac{1}{4\pi }\frac{{\mathrm{e}}^{i\xi |\mathbf{x}|}}{|\mathbf{x}|}{\delta }_{j\ell }+\frac{1}{{\xi }^2}{\mathrm{\partial }}_j{\mathrm{\partial }}_{\ell }\mathrm{\Phi }\left(\mathbf{x}\right).$ A straightforward computation of the derivatives yields $G_{j\ell }\left(\mathbf{x}\right)=\frac{{\mathrm{e}}^{i\xi |\mathbf{x}|}}{4\pi }\left(\left(-\frac{1}{|\mathbf{x}|}+\frac{i}{\xi |\mathbf{x}{|}^2}-\frac{1}{{\xi }^2|\mathbf{x}{|}^3}\right){\delta }_{j\ell }-\left(\frac{1}{|\mathbf{x}|}+\frac{3i}{\xi |\mathbf{x}{|}^2}-\frac{3}{{\xi }^2|\mathbf{x}{|}^3}\right)\frac{x_j{x}_{\ell }}{|\mathbf{x}{|}^2}\right)$ which can be estimated by $\left|G_{j\ell }\left(\mathbf{x}\right)\right|\le \frac{1}{4\pi }\left(\frac{2}{|\mathbf{x}|}+\frac{4}{|\xi ||\mathbf{x}{|}^2}+\frac{4}{|\xi {|}^2|\mathbf{x}{|}^3}\right).$ This yields (30(30) $\left|G_{j\ell }\left(\mathbf{x}\right)\right|\le p\left(\frac{1}{|\mathbf{x}|}\right)\mathrm{\forall }1\le j,\ell \le 3,\mathrm{\forall }\mathbf{x}\ne 0.$(30) ) where the coefficients of the polynomial p depend on the wavenumber ξ. Estimate (31(31) $\left|{\mathrm{\partial }}_m{G}_{j\ell }\left(\mathbf{x}\right)\right|\le q\left(\frac{1}{|\mathbf{x}|}\right)\mathrm{\forall }1\le j,\ell ,m\le 3.$(31) ) follows in the same way.

Theorem 4.5

Let $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ be the sphere of radius $\alpha >0$ and centre ${\mathbf{x}}_0\in \mathrm{\Omega }$ and assume that $\mathrm{dist}(B,\mathrm{\Gamma })\ge \beta >0$. Then there is a constant $C=C(\beta ,f,\mathbf{E})>0$ such that for any $\mathbf{z}\in \mathrm{\Gamma },$ (32) $\begin{array}{rl}|\mathbf{T}\left(\mathbf{z}\right)|& \le Cp\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right)\mathrm{and}\end{array}$(32) (33) $\begin{array}{rl}\left|{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{z}\right)\right|& \le Cq\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right),\end{array}$(33) where p and q are the polynomials from Lemma 4.4.

Proof.

We recall that $\mathbf{T}\left(\mathbf{z}\right)=-2\left({\int }_{\mathrm{\Omega }}G(\mathbf{x}-\mathbf{z})f\left(\mathbf{x}\right)\mathbf{E}\left(\mathbf{x}\right)\mathrm{d}\mathbf{x}\right)\times \mathbf{n}.$ Since the background parameters ϵ and σ are constant, the solution $\mathbf{E}$ of problem (${\mathcal{M}}_v$) belongs to $H^2\left(\mathrm{\Omega }\right)\hookrightarrow {\mathcal{C}}^0\left(\overline{\mathrm{\Omega }}\right)$ and $\parallel \mathbf{E}{\parallel }_{\mathrm{\infty },B}$ is well defined. Taking into account that $\mathrm{supp}\left(f\right)\subset B$, we get $|\mathbf{T}\left(\mathbf{z}\right)|\le 2\parallel f{\parallel }_{\mathrm{\infty },B}\underset{1\le j,\ell \le 3}{max}\parallel G(\cdot -\mathbf{z}){\parallel }_{\mathrm{\infty },B}\parallel \mathbf{E}{\parallel }_{\mathrm{\infty },B}\mathrm{vol}\left(B\right).$ Now, let $\mathbf{x}\in B$ and $\mathbf{z}\in \mathrm{\Gamma }$. According to the assumptions on B, we have $|\mathbf{x}-\mathbf{z}|\ge \beta >0$. Moreover, since $|{\mathbf{x}}_0-\mathbf{x}|\le \alpha$, we get $|{\mathbf{x}}_0-\mathbf{z}|\le \alpha +|\mathbf{x}-\mathbf{z}|\le |\mathbf{x}-\mathbf{z}|\left(\frac{\alpha }{|\mathbf{x}-\mathbf{z}|}+1\right)\le |\mathbf{x}-\mathbf{z}|\left(\frac{\alpha }{\beta }+1\right).$ The constant $\left(\alpha /\beta \right)+1$ can be majored for all possible values of α by $\left(\mathrm{diam}\left(\mathrm{\Omega }\right)/2\beta \right)+1$. Noticing that the polynomial p in Lemma 4.4 has positive coefficients allows to write $p\left(\frac{1}{|\mathbf{x}-\mathbf{z}|}\right)\le Cp\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)$ and completes the proof of the estimate (32(32) $\begin{array}{rl}|\mathbf{T}\left(\mathbf{z}\right)|& \le Cp\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right)\mathrm{and}\end{array}$(32) ).

In order to obtain the estimate for the surface divergence of $\mathbf{T}$, we notice that ${\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{z}\right)=\mathbf{n}\cdot \left(\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{T}\right)\left(\mathbf{z}\right),$ for any $\mathbf{z}\in \mathrm{\Gamma }$. But the computation of $\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{T}$ involves the first order derivatives of the fundamental solution $G(\cdot -\mathbf{z})$ which are estimated with the help of the polynomial q (see (31(31) $\left|{\mathrm{\partial }}_m{G}_{j\ell }\left(\mathbf{x}\right)\right|\le q\left(\frac{1}{|\mathbf{x}|}\right)\mathrm{\forall }1\le j,\ell ,m\le 3.$(31) )). This yields (33(33) $\begin{array}{rl}\left|{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{z}\right)\right|& \le Cq\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right),\end{array}$(33) ) noticing again that q has positive coefficients.

4.3. Estimates for ${\mathbf{E}}^1\times \mathbf{n}$

We deduce from the previous section the following estimates that yield relations between the tangential trace of the sensitivity and characteristics of the perturbation located in the ball B. As before, let $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ be the ball of radius α and centre ${\mathbf{x}}_0\in \mathrm{\Omega }$. We assume that $\mathrm{dist}(B,\mathrm{\Gamma })\ge \beta$ for a fixed constant $\beta >0$. We denote by $\hat{\mathbf{x}}$ the projection of ${\mathbf{x}}_0$ on the boundary Γ.

Proposition 4.6

Let ${\mathbf{E}}^1$ denote the solution of the sensitivity Equation ($\mathcal{S}$). Under the assumptions of Theorem 4.5, we have (34) $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{0,\mathrm{\Gamma }}\le C\left(p\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)+q\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)\right)\mathrm{vol}\left(\mathrm{B}\right),$(34) where p and q are the polynomials of Lemma 4.4, and C>0 is a constant independent from α and d.

Proof.

First notice that $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{0,\mathrm{\Gamma }}\le \mathrm{meas}(\mathrm{\Gamma }{)}^{1/2}\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{\mathrm{\infty },\mathrm{\Gamma }}$ since ${\mathbf{E}}^1$ is continuous on $\overline{\mathrm{\Omega }}$ according to the regularity assumptions. We next recall that $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{\mathrm{\infty },\mathrm{\Omega }}\le C_{\mathcal{M}}\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)},$ where $C_{\mathcal{M}}$ denotes the norm of $(\mathcal{I}-\mathcal{M}{)}^{-1}$. According to Theorem 4.5, $|\mathbf{T}\left(\mathbf{z}\right)|$ (resp. $|{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{z}\right)|$) can be estimated for any $\mathbf{z}\in \mathrm{\Gamma }$ by the polynomial p (resp. q) which has positive coefficients and satisfies $p\left(0\right)=0$ (resp. $q\left(0\right)=0$). Notice further that $1/|{\mathbf{x}}_0-\mathbf{z}|$ takes its maximum value for $\mathbf{z}=\hat{\mathbf{x}}$, and so do $p\left(1/|{\mathbf{x}}_0-\mathbf{z}|\right)$ and $q\left(1/|{\mathbf{x}}_0-\mathbf{z}|\right)$. Then, (34(34) $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{0,\mathrm{\Gamma }}\le C\left(p\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)+q\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)\right)\mathrm{vol}\left(\mathrm{B}\right),$(34) ) follows from (32(32) $\begin{array}{rl}|\mathbf{T}\left(\mathbf{z}\right)|& \le Cp\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right)\mathrm{and}\end{array}$(32) ) and (33(33) $\begin{array}{rl}\left|{\mathrm{div}}_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{z}\right)\right|& \le Cq\left(\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}\right)\mathrm{vol}\left(B\right),\end{array}$(33) ).

The next estimate states that ${\mathbf{E}}^1\times \mathbf{n}$ behaves similar to $\mathbf{T}$ on the boundary Γ:

Proposition 4.7

Under the assumptions of Theorem 4.5, there are constants $\tilde{C}>0$ and c>0 independent from $\mathbf{z}$ and the ball $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ such that for any $\mathbf{z}\in \mathrm{\Gamma },$ (35) $|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|\le \frac{\tilde{C}}{|{\mathbf{x}}_0-\mathbf{z}|}+c.$(35)

Proof.

Let $\mathbf{a}={\mathbf{E}}^1\times \mathbf{n}=(\mathcal{I}-\mathcal{M}{)}^{-1}\mathbf{T}$. The following identity can be easily verified: $(\mathcal{I}-\mathcal{M}{)}^{-1}\mathbf{T}=\mathbf{T}+(\mathcal{I}-\mathcal{M}{)}^{-1}\mathcal{M}\mathbf{T}.$ Now, define the constant $K_{\mathcal{M}}=\parallel (\mathcal{I}-\mathcal{M}{)}^{-1}\mathcal{M}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}$. One gets $\frac{1}{K_{\mathcal{M}}}\left|\mathbf{a}\left(\mathbf{z}\right)\right|\le \frac{1}{K_{\mathcal{M}}}\left|\mathbf{T}\left(\mathbf{z}\right)\right|+\frac{1}{K_{\mathcal{M}}}\parallel (\mathcal{I}-\mathcal{M}{)}^{-1}\mathcal{M}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}$ or, equivalently, $\left|\mathbf{a}\left(\mathbf{z}\right)\right|\le \left|\mathbf{T}\left(\mathbf{z}\right)\right|+K_{\mathcal{M}}\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}.$ An estimate for $|\mathbf{T}\left(\mathbf{z}\right)|$ has been obtained in Theorem 4.5. Here, we only keep the dominating terms. Since $|{\mathbf{x}}_0-\mathbf{z}|\ge \beta$, we have $\frac{1}{|{\mathbf{x}}_0-\mathbf{z}{|}^n}\le \frac{1}{{\beta }^{n-1}}\frac{1}{|{\mathbf{x}}_0-\mathbf{z}|}$ for any $n\ge 1$ which yields the first term on the right-hand side of (35(35) $|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|\le \frac{\tilde{C}}{|{\mathbf{x}}_0-\mathbf{z}|}+c.$(35) ).

In order to get an estimate of $\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}$, we state as in the proof of Proposition 4.6 that $p\left(1/|{\mathbf{x}}_0-\mathbf{z}|\right)$ and $q\left(1/|{\mathbf{x}}_0-\mathbf{z}|\right)$ reach their maximum values at $\mathbf{z}=\hat{\mathbf{x}}$. Therefore, we have $\parallel \mathbf{T}{\parallel }_{{\mathcal{T}}_d\left(\mathrm{\Gamma }\right)}\le \frac{C_1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}$ and the right-hand side of the above inequality can be majored by the constant $c=\left(C_1/\beta \right)>0$ independently from the ball B.

5. The sensitivity analysis for solving an inverse problem

The inverse medium problem, that we are interested in, is to localize inhomogeneities in the electrical parameters of the medium from total or partial boundary data on Γ for a given (boundary) source term, at a fixed frequency ω. The setting is similar to the ones in [18,23,24]. Our inverse method is based on the information obtained by the sensitivity analysis.

5.1. An inverse medium problem

We assume that Ω is filled with a medium of electrical permittivity and conductivity (36) ${\epsilon }_p=\epsilon +a_{\epsilon }{\varrho }_{\epsilon }\text{and}{\sigma }_p=\sigma +a_{\sigma }{\varrho }_{\sigma },$(36) where ${\varrho }_{\epsilon }$ and ${\varrho }_{\sigma }$ are the characteristic functions of a perturbation in the homogeneous background parameters ϵ and σ.

Let ${\mathbf{E}}_{\mathit{\eta }}$ be a plane wave of direction $\mathit{\eta }\in {\mathbb{R}}^3$, ${\mathbf{E}}_{\mathit{\eta }}\left(\mathbf{x}\right)={\mathit{\eta }}^{\perp }{\mathrm{e}}^{i\mathit{\eta }\cdot \mathbf{x}}.$ Here, ${\mathit{\eta }}^{\perp }$ is a unit vector orthogonal to $\mathit{\eta }$. ${\mathbf{E}}_{\mathit{\eta }}$ is acting as a boundary source term for the Neumann trace. Notice that other source terms could have been considered. In the absence of inhomogeneities, the electric field $\mathbf{E}$ is solution to (37) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}-k^2\frac{1}{{\epsilon }_0}\left(\epsilon +i\frac{\sigma }{\omega }\right)\mathbf{E}=0\mathrm{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{E}\times \mathbf{n}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathit{\eta }}\times \mathbf{n}\mathrm{on}\mathrm{\Gamma }.\end{array}$(37) Next, consider the electric field ${\mathbf{E}}_p$ in the presence of the inhomogeneities and subject to the same boundary data. ${\mathbf{E}}_p$ is solution to the perturbed problem (38) $\begin{array}{rl}& \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_p-k^2\frac{1}{{\epsilon }_0}\left({\epsilon }_p+i\frac{{\sigma }_p}{\omega }\right){\mathbf{E}}_p=0\mathrm{in}\mathrm{\Omega },\\ & \mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_p\times \mathbf{n}=\mathbf{c}\mathbf{u}\mathbf{r}\mathbf{l}{\mathbf{E}}_{\mathit{\eta }}\times \mathbf{n}\mathrm{on}\mathrm{\Gamma }.\end{array}$(38) We focus on perturbations of small amplitude and simple geometries (sphere, ellipsoid, etc. ). The inverse problem consists in retrieving their centres and volumes from boundary data $({\mathbf{E}}_p-\mathbf{E})\times \mathbf{n}$. According to Taylor expansion (13(13) $\mathbf{E}(\cdot ,\tau +h\varrho )-\mathbf{E}(\cdot ,\tau )\approx h{D}_{\varrho }\mathbf{E}(\cdot ,\tau ).$(13) ), the boundary data are related to the sensitivity data ${\mathbf{E}}^1\times \mathbf{n}$ of the electric field with respect to these perturbations. The estimates of Section 4, completed by a numerical study, allow to find three explicit relations between the data and the characteristics of the inhomogeneities. These relations and their link to the theoretical estimates are presented in Section 5.2. They have been validated numerically for a large number of configurations and the results of this verification are presented in Section 5.3.

5.2. Explicit relations between data and inhomogeneities

We infer from the results of Section 4.3 that the boundary sensitivity data ${\mathbf{E}}^1\times \mathbf{n}$ behave approximately as follows: $|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|\approx \frac{C}{|\mathbf{z}-{\mathbf{x}}_0|}+c,\mathrm{\forall }\mathbf{z}\in \mathrm{\Gamma },R1$ where C and c are (unknown) positive constants. In other words, the modulus of the data on the boundary takes its maximum value at the projection $\hat{\mathbf{x}}$ of the perturbation's centre ${\mathbf{x}}_0$ and should be small far away from $\hat{\mathbf{x}}$. This allows to retrieve the position of the projection $\hat{\mathbf{x}}$ from the modulus of the data (see Figure 4).

Figure 4. Illustration of the relation ($R1$). The part of the boundary Γ where the largest values of the modulus of the data are reached.

Next, we aim to reconstruct the depth of the perturbation. Together with the projection $\hat{\mathbf{x}}$ obtained in the previous step, this yields the centre ${\mathbf{x}}_0$. To this end, let $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ be the ball of radius α and centre ${\mathbf{x}}_0$ and denote by $d=|{\mathbf{x}}_0-\hat{\mathbf{x}}|$ the distance of the centre of the perturbation to its projection on the boundary. Let ${\mathbf{E}}_d^1$ be the sensitivity in the direction ${\varrho }_d=(0,1_B)$. We introduce the set (39) ${\mathrm{\Gamma }}_{\vartheta }\left(d\right)=\left\{\mathbf{z}\in \mathrm{\Gamma };|({\mathbf{E}}_d^1\times \mathbf{n})\left(\mathbf{z}\right)|\ge \vartheta {∥{\mathbf{E}}_d^1\times \mathbf{n}∥}_{\mathrm{\infty },\mathrm{\Gamma }}\right\},$(39) where $0<\vartheta <1$ is a fixed threshold. The following relation has been obtained from numerical simulations: $\frac{\mathrm{meas}\left({\mathrm{\Gamma }}_{\vartheta }\right(d\left)\right)}{\mathrm{meas}\left(\mathrm{\Gamma }\right)}\approx \frac{1}{1+{\mathrm{e}}^{p_{\vartheta }\left(d\right)}}R2,$ where $p_{\vartheta }$ is a (known) polynomial function of degree 4 that is independent from the radius α of the perturbation (see Figure 6). Since the left-hand side of relation ($R2$) can be computed from the boundary data ${\mathbf{E}}_d^1\times \mathbf{n}$, relation ($R2$) allows to compute the depth d by inversion of the function $1/(1+{\mathrm{e}}^{p_{\vartheta }\left(d\right)})$.

Remark 5.1

Relation ($R2$) could be interpreted in the following probabilistic way. Assume that the boundary point $\mathbf{z}\in \mathrm{\Gamma }$ is chosen randomly following a uniform distribution. Then, the term (40) $\frac{|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|}{\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{\mathrm{\infty },\mathrm{\Gamma }}}$(40) may be interpreted as a random variable X that follows a probabilistic law described by a density function $f_X(\mathbf{z};d)$ depending on the depth d. Consequently, the left-hand side of relation ($R2$) is given by $\frac{\mathrm{meas}\left({\mathrm{\Gamma }}_{\vartheta }\right(d\left)\right)}{\mathrm{meas}\left(\mathrm{\Gamma }\right)}\approx \mathbb{P}(X\ge \vartheta )=1-F_X(\vartheta ;d)$ where $F_X(\vartheta ;d)$ is the cumulative distribution function associated with the density $f_X(\mathbf{z};d)$. Now, we infer from Section 4.3, that on the one hand, $|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|$ takes its maximum at the point $\mathbf{z}=\hat{\mathbf{x}}$, and, on the other, $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{\mathrm{\infty },\mathrm{\Gamma }}$ behaves roughly speaking as $1/|{\mathbf{x}}_0-\hat{\mathbf{x}}|=1/d$. If we assume that the random variable X follows a logistic law (which is consistent with the theoretical results), $X\sim \frac{{\mathrm{e}}^{-x/d}}{d(1+{\mathrm{e}}^{-x/d}{)}^2}$ we get $\frac{\mathrm{meas}\left({\mathrm{\Gamma }}_{\vartheta }\right(d\left)\right)}{\mathrm{meas}\left(\mathrm{\Gamma }\right)}=1-\frac{1}{1+{\mathrm{e}}^{-\vartheta /d}}=\frac{1}{1+{\mathrm{e}}^{\vartheta /d}}.$ We may notice that this behaviour fits qualitatively with the numerical observations (see Figure 6). For a better concordance with the numerical results, however, the quantity $1-F_X$ has been fitted with the help of a polynomial function $p_{\vartheta }$ such that $1-F_X(\vartheta ;d)\approx \frac{1}{1+{\mathrm{e}}^{p_{\vartheta }\left(d\right)}}.$ This yields relation ($R2$). Notice also that both the numerator and the denominator in (40(40) $\frac{|({\mathbf{E}}^1\times \mathbf{n})\left(\mathbf{z}\right)|}{\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{\mathrm{\infty },\mathrm{\Gamma }}}$(40) ) depend linearly on the volume of the perturbation. Consequently, relation ($R2$) should not behave on $\mathrm{vol}\left(B\right)$ which is confirmed by the numerical results.

Finally, we aim to obtain the volume of the perturbation. To this end, we recall that according to (34(34) $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{0,\mathrm{\Gamma }}\le C\left(p\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)+q\left(\frac{1}{|{\mathbf{x}}_0-\hat{\mathbf{x}}|}\right)\right)\mathrm{vol}\left(\mathrm{B}\right),$(34) ), the $L^2$-norm of the boundary data ${\mathbf{E}}^1\times \mathbf{n}$ is related to $\mathrm{vol}\left(B\right)$ by a linear relation with a constant depending on $1/|{\mathbf{x}}_0-\hat{\mathbf{x}}|=1/d$. This constant can be fitted numerically and leads to the following relation between the data and the volume: $\parallel {\mathbf{E}}^1\times \mathbf{n}{\parallel }_{0,\mathrm{\Gamma }}\approx {\mathrm{e}}^{p\left(d\right)}\mathrm{vol}\left(B\right),R3$ where this time p is a polynomial function of degree 2.

The relations ($R1$)–($R3$) have been obtained from estimates on the right-hand side of an appropriate integral equation with the sensitivity as unknown. Their precise formulation is based on the numerical fitting of polynomial parameters from a data base. To the best of our knowledge, this point of view has not yet been adopted in literature. Integral operators have been used in [44,45] to develop asymptotic expansions that allow to retrieve information about the localization and shape of small-volume perturbations in the parameters of both the conductivity and Helmholtz equation.

5.3. Numerical verification of the relations ($R1$), ($R2$), and ($R3$)

The numerical verification of the above explicit relations has been done in the case where the computational domain is the unit sphere. Let us consider a single spherical perturbation $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ of radius $\alpha >0$ and centre ${\mathbf{x}}_0\in \mathrm{\Omega }$. We study the Gâteaux derivative of the electric field in the direction $\varrho =(0,{\mathbf{1}}_B)$ for sample values of ${\mathbf{x}}_0$ and α. For each couple $({\mathbf{x}}_0,\alpha )$, we compute the tangential trace ${\mathbf{E}}^1\times \mathbf{n}$ on Γ where ${\mathbf{E}}^1$ is the solution of the sensitivity Equation ($\mathcal{S}$).

The physical parameters are the same as in Section 3.4. We fix h=0.14 and $N_e=114457$. In the sequel, we present some illustrations of the three relations, and explain in which way the polynomial functions of relations ($R2$) and ($R3$) have been obtained.

5.3.1. Projection of the perturbation's centre

Property ($R1$) is illustrated in Figure 5. We report the modulus of the sensitivity ${\mathbf{E}}^1\times \mathbf{n}$ at the boundary in the direction $\varrho =(0,{\mathbf{1}}_B)$ for a perturbation of the conductivity centred at ${\mathbf{x}}_0=(-0.85,0,0)$ with radius $\alpha =0.1$. In order to improve the readability of the image, we use an equirectangular projection and create an image of ratio 2:1. The coordinates $(x,y)$ of each pixel are mapped to $(\theta ,\varphi )\in [0,2\pi [\times [-(\pi /2),(\pi /2\left)\right]$. Then, the colour of the pixel corresponds to the value of the function at coordinates $(\theta ,\varphi )$ on the sphere. We observe that the trace ${\mathbf{E}}^1\times \mathbf{n}$ is localized in a neighbourhood of the point $\hat{\mathbf{x}}=(-1,0,0)$ which is the projection of the perturbation's centre ${\mathbf{x}}_0$ on Γ. The same localization property is observed for any tested couple $({\mathbf{x}}_0,\alpha )$.

Figure 5. Illustration of the relation ($R1$). Top: modulus of the trace $|{\mathbf{E}}^1\times \mathbf{n}|$ on the unfolded sphere. Bottom: an amplitude peak appears at the point $\hat{\mathbf{x}}$ (view from side).

Figure 6. Illustration of relation ($R2$). Left: fitting of the numerical data to a curve of shape $1/(1+{\mathrm{e}}^{p_{\vartheta }\left(d\right)})$, $\alpha =0.1,\vartheta =0.2$. Right: illustration of the independence of the relation over the perturbation's size α.

5.3.2. Depth of the perturbation

In order to verify relation ($R2$) and determine numerically the coefficients of the involved polynomial function $p_{\vartheta }$, we fix $\hat{\mathbf{x}}=(-1,0,0)$ and a radius α, and consider different centres ${\mathbf{x}}_0\left(d\right)=(1-d)\hat{\mathbf{x}},d\in ]\alpha ,1[$. For each sample value d, we compute numerically the ratio $\mathrm{meas}\left({\mathrm{\Gamma }}_{\vartheta }\right(d\left)\right)/\mathrm{meas}\left(\mathrm{\Gamma }\right)$ from the boundary data ${\mathbf{E}}_d^1\times \mathbf{n}$ associated with the direction $\varrho =(0,{\mathbf{1}}_B)$ where $B=B_{\alpha }\left({\mathbf{x}}_0\right(d\left)\right)$. We get the distribution function given by ($R2$) where the polynomial function $p_{\vartheta }$ has been chosen to fit the data (see Figure 6, left). Next, we perform the same test for different radii α. It turns out that the behaviour of the plotted curve is independent from α (Figure 6, right). Therefore, the same polynomial $p_{\vartheta }$ allows to retrieve the depth d by inverting the relation ($R2$) independently from the (unknown) radius α.

5.3.3. Volume of the perturbation

Finally, let us study relation ($R3$). In Figure 7 (left) we plot the $L^2$-norm of the boundary data ${\mathbf{E}}^1\times \mathbf{n}$ in terms of the volume of the perturbation B with fixed centre ${\mathbf{x}}_0=(0,0,0)$ and different values of $\alpha \in [0,1[$, where ${\mathbf{E}}^1$ has been computed in the direction $\varrho =(0,{\mathbf{1}}_B)$. This agrees with the statement ($R3$) if we neglect the constant term in the affine relation. However, the linearity constant in ($R3$) is likely to depend on the depth d. We thus check numerically the value of the constant for different depths d. These data fit to a relation of exponential shape, (41) $K\left(d\right)={\mathrm{e}}^{p\left(d\right)}$(41) with a given polynomial p of degree 2.

Figure 7. Illustration of the linear relation ($R3$) (left). Evolution of the linearity constant with respect to the depth d (right).

Figure 8. Example of a thresholded modulus of the trace ${\mathbf{E}}^1\times \mathbf{n}$. The surfacic perturbations of largest amplitudes are shown in white.

6. The localization algorithm

We develop a reconstruction algorithm of interior perturbations from the knowledge of the tangential trace ${\mathbf{E}}^1\times \mathbf{n}$ of the sensitivity (of the electric field). The algorithm is based on the relations ($R1$), ($R2$), and ($R3$) of the former section. Notice that the proposed algorithm could easily be applied to the case where the input data is the tangential trace of the perturbed field according to Taylor expansion (13(13) $\mathbf{E}(\cdot ,\tau +h\varrho )-\mathbf{E}(\cdot ,\tau )\approx h{D}_{\varrho }\mathbf{E}(\cdot ,\tau ).$(13) ) which is valid for small-amplitudes.

6.1. Database generation

A first step of the inversion algorithm consists in simulating a large number of possible spherical perturbations with same projection $\hat{x}$ and defined by their depth d and their volume. We compute the corresponding boundary data ${\mathbf{E}}^1\times \mathbf{n}$. By varying the inhomogeneity's parameters, we are able to estimate the coefficients of the polynomial functions in ($R2$) and ($R3$). These coefficients are stored and will be used in the resolution phase. The database is generated with a given mesh ${\mathcal{M}}_{\mathrm{data}}$. It is important to notice that this preliminary step is required only once for a given computational domain Ω. It is described in Algorithm 1.

6.2. The inversion procedure

The algorithm is presented hereafter in the case of one single perturbation in the conductivity. The aim is to retrieve the following parameters from given (synthetic) boundary data ${\mathbf{E}}^1\times \mathbf{n}$:

• $\hat{\mathbf{x}}$, the projection on the sphere of ${\mathbf{x}}_0$, the centre of the interior perturbation,

• $d=|{\mathbf{x}}_0-\hat{\mathbf{x}}|$, the depth of the centre ${\mathbf{x}}_0$,

• the volume $v_B$ of the perturbation which yields the radius α in the case of a spherical perturbation.

6.3. Numerical simulations

6.3.1. Generation of synthetic boundary data

In the absence of measurements, we generate discrete synthetic boundary data in the following way.

• Fix a perturbation $B=B_{\alpha }\left({\mathbf{x}}_0\right)$ with given centre ${\mathbf{x}}_0$ and radius α.

• For any direction $\mathit{\eta }$ in a given set $N_m$ of m directions

  • Compute the boundary data ${\mathbf{E}}^1\times \mathbf{n}$ where ${\mathbf{E}}^1$ is the sensitivity in the direction $\varrho =(0,{\mathbf{1}}_B)$ for the source term ${\mathbf{E}}_{\eta }$.

  • Add some noise to the data (Section 6.3.3 only).

• Take $\mathbf{a}=({\mathbf{E}}^1\times \mathbf{n})+\mathrm{noise}$ as input data for the inversion algorithm.

The following sets of directions are used in the sequel: $\begin{array}{rl}{N}_1& =\left\{\right(0,1,0\left)\right\},\\ N_6& =N_1\cup \left\{\right(1,0,0),(0,0,1),(-1,0,0),(0,-1,0),(0,0,-1\left)\right\},\\ N_{14}& =N_6\cup \left\{\right(1,1,1),(-1,1,1),(1,-1,1),(1,1,-1),\\ & \phantom{=N_6\cup }(-1,-1,1),(-1,1,-1),(1,-1,-1),(-1,-1,-1)\}.\end{array}$ In order to avoid an inverse crime (in the sense of [36, p. 133]), the boundary data are computed on a tetrahedral mesh ${\mathcal{M}}_{\mathrm{inv}}$ of size h=0.12 ($N_e=167402$ edges) which is different from the mesh used to generate the database. Notice also that we focus on perturbations of the conductivity parameter. We refer to Section 6.3.4 where both parameters, ϵ and σ, undergo a perturbation.

6.3.2. Spherical perturbation

We first apply Algorithm 2 in a homogeneous background medium containing a spherical perturbation in the conductivity. We keep the physical settings defined in Section 3.4. In the case of multiple directions, we apply the algorithm on each direction $\mathit{\eta }\in N_m$ and compute the mean of the results.

We first test our algorithm with a spherical perturbation centred at ${\mathbf{x}}_0=(-0.7,0,0)$ and of radius 0.2. The parameters to retrieve are then $\hat{\mathbf{x}}=(-1,0,0)$ (or, in spherical coordinates, $\hat{\mathbf{x}}=(\pi ,0)$), d=0.3 and $\alpha =0.2$. We report the results in Table 1. The approximations of $\hat{\mathbf{x}}$, d and α are respectively denoted by ${\hat{\mathbf{x}}}_h$, $d_h$ and ${\alpha }_h$. We choose the Euclidian norm of the spherical coordinates of the points $\hat{\mathbf{x}}$ and ${\hat{\mathbf{x}}}_h$ to compute the projection error.

Table 1. Numerical resolution of the inverse problem in a homogeneous medium. One spherical perturbation in the conductivity, of centre ${\mathbf{x}}_0=(-0.7,0,0)$ (i.e. $\hat{\mathbf{x}}=(\pi ,0)$) and radius $\alpha =0.2$.

Incidences	$N_1$	$N_6$	$N_{14}$
${\hat{\mathbf{x}}}_h$	(3.156,0.018)	(3.152,0.015)	(3.151,0.009)
${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	7.502e−03	5.713e−03	4.173e−03
$d_h$	0.316	0.311	0.313
${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	5.459e−02	3.728e−02	4.193e−02
${\alpha }_h$	0.214	0.214	0.210
${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$	6.943e−02	7.158e−02	5.037e−02

We retrieve the projection centre $\hat{\mathbf{x}}$ with a very good accuracy (less than 1% error). The approximation error on $\hat{\mathbf{x}}$ decreases with respect to the number of incident waves. The other perturbation's characteristics d and α are well approximated, too (about 4% to 7% error). Their approximation does not really depend on the number of waves. In order to keep a reasonable number of computations (resp. measurements in the context of biomedical applications that we have in mind), we decide to work from now on with the set $N_6$ of six different incident waves (unless specified otherwise).

We now test a spherical perturbation which is centred at a different point. We keep the parameters d=0.3 and $\alpha =0.2$, and consider the centre ${\mathbf{x}}_0=(0,-0.7,0)$ (which corresponds to $\hat{\mathbf{x}}=(0,-1,0)$. In Table 2, we compare the errors obtained with the two configurations. The relative errors of the two approximations are comparable. Changing the position of the perturbation does not affect the quality of the approximation. From now on, we thus will consider perturbations centred on the x-axis (unless indicated otherwise).

Table 2. Comparison of the results for different centres. Relative errors for a spherical perturbation of radius $\alpha =0.2$ in the conductivity.

Perturbation's centre	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
${\mathbf{x}}_0=(-0.7,0,0)$	5.713e−03	3.728e−02	7.158e−02
${\mathbf{x}}_0=(0,-0.7,0)$	5.551e−03	5.528e−02	9.906e−02

We next apply the algorithm for different depths and volumes and report the errors in Table 3. It may be stated that the accuracy of the reconstruction does not depend significantly on the depth or volume of the perturbation.

Table 3. Comparison of the results for different depths and radii. Relative errors for a spherical perturbation in the conductivity.

Perturbation's centre	r	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
${\mathbf{x}}_0=(-0.7,0,0)$	0.2	5.713e−03	3.728e−02	7.158e−02
${\mathbf{x}}_0=(-0.4,0,0)$	0.2	9.658e−03	6.661e−03	1.354e−02
${\mathbf{x}}_0=(-0.4,0,0)$	0.4	1.110e−02	9.791e−03	1.057e−02

6.3.3. Noisy data

An important feature in numerical reconstruction is noise robustness. Noisy synthetic data are generated in the following way. Let ${\left(z_k\right)}_k={(x_k+i{y}_k)}_k$ be the degrees of freedom of the sensitivity ${\mathbf{E}}^1$. We add noise independently in the real and imaginary parts. Let ${\left(n_k\right)}_k$ be a vector of real random numbers following a normal law $\mathcal{N}(0,\sigma )$. Then, we generate additive noise in the real part with the following formula: $x_k^{\mathrm{noise}}=x_k+\frac{M_x-m_x}{2}{n}_k+\frac{m_x+M_x}{2},$ where ${\displaystyle m_x=\underset{k}{min}{x}_k}$ and ${\displaystyle M_x=\underset{k}{max}{x}_k}$. The imaginary part is jittered in the same way.

In Table 4, we compare the results obtained by our algorithm using noisy or non-noisy data (with a single incident wave of direction $(0,1,0)$). As expected, noise affects the reconstruction. A numerical observation of the noisy boundary data allows to get a deeper insight on the effect of noise and to propose a strategy for noise reduction. Indeed, in a normally distributed noise vector, a few coefficients may have large values. Consequently, there can be a small number of points on the boundary located far away from the projection $\hat{x}$, and at which the modulus $|{\mathbf{E}}^1\times \mathbf{n}|$ is greater or equal than at the points around the projection $\hat{\mathbf{x}}$ (see Figure 9). This causes problems to the first step of Algorithm 2. A solution is to apply a post-processing which consists in identifying outliers and deciding whether they should be retained or rejected. As before, we retrieve the points where the largest values of the modulus are reached. An outlier detection is applied to remove only those points which are not in the neighbourhood of the principal peak (i.e. centred at $\hat{\mathbf{x}}$). For instance, the deviation around the median can be used (see [47]). By applying this simple algorithm twice, we reduce the number of outliers significantly. In Table 5, we report the reconstruction results obtained from the noisy data that have undergone the outlier detection, compared to results from non-noisy data. Up to 2% of noise, the results are not affected by the noise, and the projection $\hat{x}$ is retrieved with a similar precision even for 5% of noise. At 10% of noise, the projection is found with 12% error, but the precision of the other parameters is no longer significant. This first study of noisy data indicates that the detection of outliers is an interesting and promising approach, and different noise reduction methods could be tested to improve the results. This was however beyond the scope of this paper.

Figure 9. Example of the modulus of a noised trace (here with 2% of noise). Some points with a value bigger than the central peak can be observed.

Table 4. Reconstruction of one spherical perturbation centred at ${\mathbf{x}}_0=(-0.7,0,0)$ of radius r=0.2 from non-noisy or noisy data.

	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
Non-noisy data	7.502e−03	5.459e−02	7.474e−02
Noisy data ($\sigma =2\mathrm{\%}$)	2.835e−01	4.945e−01	1.243e−01

Table 5. Reconstruction of one spherical perturbationcentred at ${\mathbf{x}}_0=(-0.7,0,0)$ of radius r=0.2 from non-noisy or noisy data (after the detection of outliers).

	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
Non-noisy data	7.502e−03	5.459e−02	7.474e−02
Denoised data ($\sigma =1\mathrm{\%}$)	3.942e−03	4.587e−02	6.175e−02
Denoised data ($\sigma =2\mathrm{\%}$)	4.826e−03	4.541e−02	6.286e−02
Denoised data ($\sigma =5\mathrm{\%}$)	4.642e−03	2.244e−01	1.711e−01
Denoised data ($\sigma =10\mathrm{\%}$)	1.241e−02	9.144e−01	5.082e−01

6.3.4. Simultaneous reconstruction of both parameters

We test our algorithm for determining a perturbation in the conductivity and the permittivity, located in the same sphere of centre ${\mathbf{x}}_0=(-0.7,0,0)$ and of radius $\alpha =0.2$. This happens for instance when a stroke occurs in the brain. In Table 6, we compare the results with the ones obtained where only the conductivity is perturbed. It has to be noted that only perturbations in the conductivity have been used to generate the database. The good accuracy on the approximated centre and depth (i.e. the localization error) is preserved. We observe a slight loss of precision for the radius: for an exact radius of $\alpha =0.2$, we find ${\alpha }_h=0.242$ instead of ${\alpha }_h=0.214$ if only the conductivity undergoes a perturbation. This is due to the fact that the $L^2$-norm of the tangential trace (used in $\left(R3\right)$) is bigger if both parameters are perturbed. A solution could be to construct a database relative to perturbations on both parameters. Nevertheless, the present study shows that the reconstruction is satisfactory even if the database does not take into account information about which parameter is perturbed.

Table 6. Perturbation in both the conductivity and the permittivity, centred at ${\mathbf{x}}_0=(-0.7,0,0)$ of radius $\alpha =0.2$.

Perturbed parameter (s)	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
Conductivity	5.713e−03	3.728e−02	7.158e−02
Conductivity & permittivity	5.666e−03	3.813e−02	2.080e−01

6.3.5. Reconstruction of two perturbations

As stated in Proposition 3.7, in the case of n>1 disjoint perturbations, the total sensitivity can be separated into n sensitivities, each representing one perturbation. We use this property to handle the configuration with two or more inhomogeneities. To this end, we proceed as follows. The entry data is the trace ${\mathbf{E}}^1\times \mathbf{n}$ which is containing n surfacic perturbations (see Figure 3). The first step is to detect different amplitude clusterings. This can be achieved by applying a detection algorithm such as DBSCAN which yields the connected components of the thresholded data. DBSCAN is used in data mining (see [48]) and has the advantage that the a priori knowledge of n is not needed. Next, for each connected component, we build an artificial piecewise trace of the sensitivity: it equals the values of the original trace on the region of the connected component, and is zero otherwise. Finally, we use this new trace as the entry of our localization procedure. In Table 7, we show the results for two disjoint spherical perturbations in the conductivity. One is centred at ${\mathbf{x}}_{0,1}=(-0.18,0.41,0)$ with radius ${\alpha }_1=0.35$ and the other at ${\mathbf{x}}_{0,2}=(0,-0.7,0)$ with ${\alpha }_2=0.2$. Both perturbations are very well localized, and the volume of the biggest one is better approximated. This procedure will work if the perturbations give raise to well-separated projections on the boundary.

Table 7. Reconstruction of two disjoint perturbations in the conductivity.

Perturbation's centre	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
${\mathbf{x}}_{0,1}=(-0.18,0.41,0)$	8.556e−03	2.829e−04	3.599e−02
${\mathbf{x}}_{0,2}=(0,-0.7,0)$	2.031e−03	9.596e−02	1.334e−01

6.3.6. Ellipsoidal perturbation

In real life applications, the shape of the perturbation is in general not known. In this subsection, we want to retrieve an ellipsoidal perturbation. It is centred at ${\mathbf{x}}_0=(-0.4,0,0)$, of x-radius 0.2, y-radius 0.4 and z-radius 0.2 which yields a volume v=6.702e$-02$. We recall that only spherical perturbations have been considered to generate the database. We report the results in Table 8 and compare them with the results obtained in the case of a sphere with the same volume. This reconstruction is illustrated in Figure 10. The ellipsoidal shape does not really affect the precision of the approximations. This may indicate that the polynomials in relations ($R2$) and ($R3$) are independent from the shape of perturbation.

Figure 10. Localization of a perturbation of ellipsoidal shape. Left: induced ellipsoidal perturbation. Right: retrieved perturbation (ball of equivalent volume).

Table 8. Reconstruction of an ellipsoidal perturbation in the conductivity, centred at ${\mathbf{x}}_0=(-0.4,0,0)$ of volume $v=6.702\mathrm{e}-02$.

Perturbation's shape	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|v-v_h|}{\left|v\right|}}$
Sphere	9.482e−03	8.522e−03	1.344e−02
Ellipsoid	2.725e−04	2.993e−02	3.742e−01

6.3.7. Three-layer head spherical model

In the biomedical applications we have in mind, e.g. for the diagnostic of strokes, it is important to take into account the heterogeneity of the medium. A classical spherical head model is commonly used in the literature. This model is built of three concentric spheres representing brain, skull and scalp (see Figure 11). The parameters of this model are described in Table 9. We simulate a spherical inhomogeneity in the conductivity of the brain layer, centred at ${\mathbf{x}}_0=(-0.57,0,0)$ of radius $\alpha =0.2$. The errors are reported in Table 10 and compared to the reconstruction results of the same perturbation in a homogeneous background. The approximations are of the same order. The theory has been developed in a case of a homogeneous background but the localization algorithm still offers good results in more realistic configurations. We emphasize that the polynomials in the database were generated with the piecewise constant background conductivity.

Figure 11. Three-layer head spherical model.

Table 9. Parameters of the three-layers head model.

Layer	Brain (${\mathrm{\Omega }}_1$)	Skull (${\mathrm{\Omega }}_2$)	Scalp (${\mathrm{\Omega }}_3$)
Definition	$B_{0.87}\left(0\right)$	$B_{0.92}\left(0\right)\setminus {\mathrm{\Omega }}_1$	$\mathrm{\Omega }\setminus ({\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2)$
Permittivity (F · m^−1)	$8.854\cdot {10}^{-10}$	$3.542\cdot {10}^{-10}$	$8.854\cdot {10}^{-11}$
Conductivity (S · m^−1)	0.33	0.04	0.33

Table 10. Three-layer spherical head model. Reconstruction of a spherical perturbation in the conductivity centred at ${\mathbf{x}}_0=(-0.57,0,0)$ of radius $\alpha =0.2$.

Medium	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
Homogeneous	3.616e−03	5.432e−02	2.845e−02
Heterogeneous	1.164e−02	3.326e−02	8.000e−02

6.3.8. A realistic head model

Finally, we consider a realistic head mesh. We use the Colin27 adult brain atlas (version 2, see [49,50]). As shown in Figure 12, the tetrahedrons are smaller in some regions, to reflect the complexity of the brain. The elements are of size between 4.67e−5 and 2e−2 m, leading to a total number of $N_t=425224$ tetrahedrons and $N_e=499136$ edges. We generate a new database with this mesh in the case of one perturbation in a homogeneous background. We refer to Section 3.4 for the values of the physical parameters. Data used for the inverse problem are generated with another head mesh, slightly different, of same mesh size but with $N_e=535921$ edges. The errors are reported in Table 11. In Figure 13, we compare graphically the expected perturbation and the approximated one. These results show that the algorithm is not limited to the academic case of the unit ball, but can also be applied on more realistic geometries.

Figure 12. Realistic head model. Left: view on the boundary. Right: cut in the middle of the x-axis.

Figure 13. Realistic head mesh. Localization of a spherical perturbation in the conductivity. Left: expected perturbation. Right: result of the algorithm. The brain region, where the tetrahedrons are smaller, is shown in for more readibility.

Table 11. Realistic head mesh. Reconstruction of a spherical perturbation in the conductivity.

Domain geometry	${\displaystyle \frac{|\hat{\mathbf{x}}-{\hat{\mathbf{x}}}_h|}{\left|\hat{\mathbf{x}}\right|}}$	${\displaystyle \frac{|d-d_h|}{\left|d\right|}}$	${\displaystyle \frac{|\alpha -{\alpha }_h|}{\left|\alpha \right|}}$
Head	7.007e−03	2.708e−03	7.198e−03

An interesting fact is that the computation time is still reasonable with this mesh, which is finer than the ball we used for the other tests. In all the tests we present here, the localization procedure is achieved in a few seconds on a personal computer (quad-core processor clocked at 2.5 GHz, with 4 Gio of RAM).

7. Conclusion and future works

In this paper, we have proposed a new and efficient algorithm for localizing small-amplitude perturbations in the electric parameters of a medium from boundary field measurements at a fixed frequency. The approach is based on a rigourous sensitivity analysis of the electric field with respect to the variations of the permittivity and conductivity. Sensitivity is proportional to the boundary measurements of the physical field for perturbations of arbitrary shape but small amplitudes. We have proved, both theoretically and numerically, that the trace of the sensitivity of the electric field (on the boundary of the domain) contains relevant information on the perturbations in the medium. Up to our knowledge, this is the first time that this kind of sensitivity analysis for the 3D Maxwell equations is used in the reconstruction of parameter perturbations. From an integral equation in a homogeneous background and extensive numerical simulations, we have obtained explicit relations between the sensitivity and some characteristics (centre and volume) of the perturbations. These relations lead to a constructive algorithm for determining the centre, the depth and the volume of inhomogeneities in the permittivity and/or the conductivity of a medium. Its implementation makes use of diverse tools from scientific computing as, for example, 3D finite element discretization, geometric algorithms, or outlier detection for noise reduction. A large variety of three-dimensional numerical results attests the efficiency of the method: one or two perturbations with different locations and sizes, perturbation in the permittivity and/or conductivity, spherical or ellipsoidal defects, non-noisy or noisy data, a constant or piecewise constant background. The study of a general variable background would be an interesting perspective for future work. Furthermore, in view of biomedical applications (e.g. microwave imaging), we have also provided simulations in the case where the computational domain is the head. Two types of head models have been considered: the classical three-layer spherical model and a realistic one. In each configuration, the perturbations are localized with a very good accuracy, and information on its volume are obtained, too. This non-iterative localization procedure would be an interesting initial guess for gradient-based descent algorithms intended to retrieve the physical coefficient values in the perturbations. This is part of ongoing work.

All the results have been obtained at a fixed frequency. Considering a non-ionizing electromagnetic radiation, a possible application of our work is to discriminate between healthy and abnormal brain tissues. This can be achieved with a single frequency and changing the frequency does not provide a priori more information, excepted if the frequency dependence of the effective electric properties of the tissue is known. In this case, a multifrequency approach could be interesting. Recent works on the multifrequency electrical impedance tomography (mfEIT) have been addressed [51,52].

Disclosure statement

No potential conflict of interest was reported by the authors.