High-order functional derivatives of the diffracted field according to the permittivity-contrast function
Slimane Arhab Dimitrios Anagnostou Maminirina Joelson




]High-order functional derivatives of the diffracted field according to the permittivity-contrast function
[-0.12in*  Corresponding author: Slimane Arhab (slimane.arhab@univ-avignon.fr).]
[-0.12in1 Université d'Avignon et des Pays de Vaucluse, UMR 1114 EMMAH, 84018 Avignon Cedex, France.]











In this work, we propose to extend an approach to calculate at any order (n), the functional derivative of the diffracted field with respect to the permittivity-contrast function. These derivatives obtained for different orders are used to perform an expansion of the data according to the studied model parameter. Its convergence is discussed throughout some numerical results, obtained in the case where the forward model used to simulate the diffracted field is built in the framework of the volume integral formulation. In particular, we show that taking into account higher order derivatives improve drastically the data-fitting. The numerical application considered consists of a cylindrical object illuminated by an incident field under a TE polarisation (electric component parallel to the invariance axis). 








§ INTRODUCTION


Reconstructing the permittivity-contrast function of an object from the diffracted field is one of the most celebrated inverse problems. It is encountered in research areas such as optical digital tomographic microscopy <cit.> or the characterisation of buried objects in natural environments <cit.>. A widely studied physical configuration concerns an object illuminated by an incident field, produced by an emitter. The resulting diffracted field is then detected on receivers. This diffracted field represents the data of the inverse problem. They can be simulated by running the forward model (solution of the forward problem), with permittivity-contrast function as input. From a mathematical point of view, the forward model often acts as a nonlinear operator between permittivity-contrast function and data. Due to this nonlinear aspect data cannot be inverted directly. Moreover, this inverse problem is known to be ill-posed in the sense of Hadamard <cit.>. To circumvent this difficulty, a priori informations are introduced by different regularization procedures <cit.>. The method we are extending in this paper does not address the ill-posed aspect of inverse problems, but contributes to the iterative inverse methods by introducing a new link between a variation of the permittivity-contrast function and the corresponding variation on the data. It is expressed as a functional expansion performed for any order (n).


This expansion is reduced to its first order approximation in local optimisation methods. This is the case for the gradient method <cit.> which is based on the calculus of the Gâteaux derivative, or the Newton-Kantorovich method <cit.> which involves the calculus of the Fréchet derivative. The latter can be obtained directly by differentiating the forward model <cit.>, or indirectly by applying the adjoint method <cit.>. In this last approach, the Fréchet derivative is constructed by solving the forward problem and its associated adjoint problem. In the latter, one considers the same object subject to a fictitious incident field, obtained by simulating a backpropagation from the receivers of the difference between the reference data and data simulated at each step of the iterative inversion process.


The concept of an adjoint problem has been introduced to solve inverse problems in several research areas. In electromagnetic imaging, it has been used to reconstruct the geometric shape of a cylindrical object <cit.>. This concept has also been introduced for the purposes of acoustic imaging <cit.>, and subsequently generalized within the framework of elastic wave theory <cit.>. In global seismology, inversion schemes to characterize the earth structure also involve the resolution of an adjoint problem <cit.>.
 

A similar approach has been introduced in electromagnetism for a configuration with emitters and receivers positioned in a farfield region, and where a harmonic time dependence is assumed <cit.>. Maxwell's equations are first introduced in the sense of distributions. Then, a calculus involving the reciprocity theorem leads to elegant forms of the Fréchet derivative, performed, respectively, with respect to the permittivity-contrast function and to the geometrical shape of the studied object. Such derivatives are built by solving the forward problem, and the associated reciprocal problems obtained by switching the emitter with each receiver. This principle has been applied in different inversion schemes for applications in detection of burried objects <cit.>, optical profilometry <cit.> and reconstruction of the geometric shape of three dimensional objects in microwave regime <cit.>.


The current work aims to extend this approach to the calculus of the functional derivative for any order. The case of a three-dimensional object illuminated by an incident field in harmonic regime, is considered, where the emitter and receivers are positioned in the nearfield region <cit.>. The paper has the following outline. First, the theoretical approach involving the reciprocity theorem <cit.> is used in section <ref> to calculate the functional derivative of the first order, which gives the expression of the Fréchet derivative with respect to the permittivity-contrast function. Then the main result of this work is to show that, it is possible to apply the same approach on the remainder of the first order functional limited expansion, which gives rise to the second order functional derivative and its associated remainder. By repeating an identical calculus on this new remainder and the following ones, it becomes possible to reach the functional derivative for any order (n), and then to perform a functional expansion of the data with respect to the permittivity-contrast function. The convergence of such an expansion is studied from a numerical point of view in section <ref>, where a two-dimensional configuration of a cylindrical object illuminated by an incident field under a TE polarisation (electric component parallel to the invariance axis), is considered. In this case, data are modelled by an operator built in the framework of the volume integral formulation. The conclusion of this work is summarized in the section <ref>.



§ THEORETICAL APPROACH




 §.§ Maxwell's Equations


We consider an object of finite volume V, delimited by its surface S. We assume that this object is made of a dielectric (i.e. non magnetic), linear and isotropic medium. We treat the general case of an object with an inhomogeneous permittivity, described by a scalar function that depends on the position vector r=∑__j=1^_3 x__j x__j in Cartesian coordinates. We introduce the complex notation, and we work under the assumption of a harmonic time dependence of the form ∝  e^- iω t, where i^^2=-1. For example, the source current that generates the incident field is written : J( r,t)={Ĵ( r)  e^- iω t}, with {} designating the real part. For the sake of clarity, from now on, the factor e^- iω t is systematically omitted and the complex vector amplitude Ĵ( r) is denoted J( r). We do the same for the electric field E and the magnetic field H. Finally, we write Maxwell's equations in the sense of distributions. To do this, we rely on the reference <cit.> which deals with the practical aspect of distributions (For a more thorough discussion the interested reader may wish to consult the reference <cit.>). In this way, we get :


    {[            ∇ ×  E =  i ω μ__0  H;                                 ; ∇ ×  H = - i ω ε  E + δ_ r_ o  J;                                 ;           ∀   r∈ the whole space ].

 
Here the term δ_ r_ o≡δ( r- r_ o) is the point-like Dirac distribution. We still have ε = ε__0 ε_ r, with ε_ r( r) for r∈𝒱, denoting the inhomogeneous relative permittivity and where (ε__0,μ__0) are respectively permittivity and permeability of the vacuum.



 §.§ Calculus Method




Let us now introduce the four physical configurations (a__l), (b__m), (b__l) and (a__m) illustrated in figure Fig (<ref>). Note that they are physically independent of each other. Each one of them is described by the permittivity-function ε_ v( r) (v= a, b and r∈𝒱), by the source current J_w (w=l,m) located at r_w and by the total electric field E_ v,w. The latter is given by the sum of the incident field generated by the source current J_w and the field diffracted by the studied object, namely E_ v,w= E^ i_w+ E^ d_ v,w. Then, the system of coupled equations Eqs (<ref>), can be rewritten for the generic configuration (v__w), in the following form:


    ( v__w), {[                 ∇ ×  E_ v,w =  i ω μ__0  H_ v,w;                                                ;  ∇ ×  H_ v,w = - i ω ε_ v  E_ v,w + δ_ r_w  J_w;                                                ; ∀   r∈ the whole space, with  v= a, b and w=l,m ].


From elementary vector calculus, we obtain for the electromagnetic fields of configurations (a__l) and (b__m) the following relation:


    [ ∇ · [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] = [ E_ a,l · (∇ ×  H_ b,m) -  H_ b,m · (∇ ×  E_ a,l)] -;                                                 [ E_ b,m · (∇ ×  H_ a,l) -  H_ a,l · (∇ ×  E_ b,m)] ]
 

By replacing the curl of the electric and magnetic fields, by their respective expressions given in Eq (<ref>), we obtain:


    ∇ · [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] = - i ω (ε_ b - ε_ a)  E_ a,l ·  E_ b,m + δ_ r_m  E_ a,l ·  J_m - δ_ r_l  E_ b,m ·  J_l
 
 
Each term of the above equation is a distribution, that we can apply on a test function. After some lengthy calculations (see appendix <ref>), this equation is rewritten in the sense of functions, in the form:


    [ ∫_ C [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] ·  e__ρ   d C = - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱 +  E_ a,l( r_m) ·  J_m -;                                                                                                                 E_ b,m( r_l) ·  J_l ]


The left-hand side of the above equation in nil (see appendix <ref>), so we have:


    E_ b,m( r_l) ·  J_l -  E_ a,l( r_m) ·  J_m = - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱


In the absence of the scatterers of permittivities ε_ a and ε_ b, this equation is rewritten only with the incident fields:


    E^ i_m( r_l) ·  J_l -  E^ i_l( r_m) ·  J_m = 0
 

By subtracting (<ref>) from (<ref>), we get:


    E^ d_ b,m( r_l) ·  J_l -  E^ d_ a,l( r_m) ·  J_m = - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱


By replacing in the above equation the configuration (a__l) by the configuration (b__l), we get the following equation which is nothing but the reciprocity theorem between the configurations (b__l) and (b__m):


    E^ d_ b,m( r_l) ·  J_l -  E^ d_ b,l( r_m) ·  J_m = 0


Once this equality is introduced in equation Eq (<ref>), the latter can be rewritten in the form:


    [ E^ d_ b,l( r_m) -  E^ d_ a,l( r_m)] ·  J_m = - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱


To derive the variational form of expression Eq (<ref>), we add the assumption that ε_ b = ε_ a + δε. This implies, for the fields of the four studied configurations, the following relations: E^ d_ b,l( r_m) =  E^ d_ a,l( r_m) + δ E^ d_ a,l( r_m) and E_ b,m( r) =  E_ a,m( r) + δ E_ a,m( r). It follows that the expression Eq (<ref>) can be rewritten as:


    δ E^ d_ a,l( r_m) ·  J_m = - iω∫_𝒱 δε( r)   E_ a,l( r) ·  E_ a,m( r)   d𝒱 - iω∫_𝒱 δε( r)   E_ a,l( r) · δ E_ a,m( r)   d𝒱

 
The electric fields of the above expression are all calculated for the same permittivity ε_ a. Thus, for the sake of clarity, we omit the index (a). Moreover, the introduction into the same expression of the relations: ε = ε__0 ε_ r and χ = ε_ r - 1, leads us to rewrite it according to the permittivity-contrast function χ: 


    {[ δ E^ d_l( r_m) ·  J_m = - iωε__0∫_𝒱 δχ( r)   E_l( r) ·  E_m( r)   d𝒱 - iωε__0∫_𝒱 δχ( r)   E_l( r) · δ E_m( r)   d𝒱;                       ⇔      δ𝒟_lm = ℱ__lm^^(1)(δχ) +  o__lm(||δχ||) ,                         (in a compact form) ].

According to the left-hand side of the above expression, data element δ𝒟_lm is given by the dot product between the variation of the diffracted field δ E^ d_l at r_m (this field is diffracted by the object after interacting with the incident field emitted by the source current J_l located at r_l) and the source current J_m that would be positioned at r_m (This source current radiates an incident field that interacts with the object to give the total electric field E_m( r), for  r∈𝒱). The above calculus was first proposed in farfield configuration <cit.>, in this context the remainder o__lm(||δχ||) has been neglected. In this work, this calculus is adapted to localized emitters and receivers. Now, from the right-hand side of the above expression, and by taking into account the remainder o__lm(||δχ||), it is now possible to extract the functional derivatives up to any order (n). 



 §.§ Functional derivatives of orders (1), (2) and (3)




  ∙ First order functional derivative (Fréchet derivative) This derivative which appears directly in expression Eq (<ref>), is obtained by neglecting the remainder o__lm(||δχ||) that contains all the derivatives of higher orders:


    {[  ℱ__lm^^(1)(δχ) = - iωε__0∫_𝒱 δχ( r)   E_l( r) ·  E_m( r)   d𝒱; Where ℱ__lm^^(1) denotes the first order functional derivative ].



  ∙ Second order functional derivative This second derivative is obtained by writing: o__lm(||δχ||) = ℱ__lm^^(2)(δχ,δχ) +  o__lm(||δχ||^^2). To deduce the explicit form of ℱ__lm^^(2), let us start with the expression of o__lm(||δχ||):


    o__lm(||δχ||) = - iωε__0∫_𝒱 δχ( r)   E_l( r) · δ E_m( r)   d𝒱 = ∫_𝒱 δ E_m( r) · [- iωε__0 δχ( r)  E_l( r)   d𝒱]


The novelty of this work (to the best of the authors' knowledge) is to set [- iωε__0 δχ( r)  E_l( r)   d𝒱] =  J_l( r). Indeed, given its physical dimension, this term can be assimilated to a source current. However, to make the calculus more convenient, we introduce the reduced source current J̅_l( r), such that: J_l( r) = J̅_l( r)  d𝒱. This allows us to rewrite the expression Eq (<ref>) as: 


    o__lm(||δχ||) = - iωε__0∫_𝒱 δχ( r)   E_l( r) · δ E_m( r)   d𝒱 = ∫_𝒱 δ E_m( r) · J̅_l( r)   d𝒱


Based on equation Eq (<ref>), this integral can be transformed into:


    [ ∫_𝒱 δ E_m( r) · J̅_l( r)   d𝒱 = ∫_𝒱 {- iωε__0∫_𝒱 δχ( r')   E_m( r') · E̅_l( r,  r')   d𝒱';                               - iωε__0∫_𝒱 δχ( r')   E_m( r') · δE̅_l( r,  r')   d𝒱'}   d𝒱 ]

 
In this new expression, E̅_l( r,  r') is the reduced electric field calculated at the point r'. It is the result of an interaction between the studied object and a fictitious incident field, generated by the reduced source current J̅_l( r) located at the point r. Eq (<ref>) can also be rewritten in the form:


    [ ∫_𝒱 δ E_m( r) · J̅_l( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r')   E_m( r') · {∫_𝒱E̅_l( r,  r')  d𝒱}   d𝒱';                                - iωε__0∫_𝒱 δχ( r')   E_m( r') · {∫_𝒱δE̅_l( r,  r')  d𝒱}   d𝒱' ]


In the above expression, the integral {∫_𝒱E̅_l( r,  r')  d𝒱}= E_l, r^*_0.4(1)( r') is the resulting electric field, due to an interaction between the studied object and all the fictitious incident fields generated by the set of source currents J_l( r), where r spans the domain 𝒱. Here, the index r^*_0.4(1) has two significations: the symbol * means that we integrate in the domain with respect to the position r, in order to take precisely into account the contributions of all the source currents in the calculation of the total field. The index (1) means that the variation of the permittivity-contrast function is implicitly embedded once in the total electric field E_l, r^*_0.4(1). Then, Eq (<ref>) is rewritten in the form:


    [ ∫_𝒱 δ E_m( r) · J̅_l( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r')   E_m( r') ·  E_l, r^*_0.4(1)( r')   d𝒱';                                - iωε__0∫_𝒱 δχ( r')   E_m( r') · δ E_l, r^*_0.4(1)( r')   d𝒱' ]


To simplify the notation, let's replace r' by r in the right hand side of the above equation. Eqs (<ref>,<ref>) lead to the following result:


    {[                     - iωε__0∫_𝒱 δχ( r)   E_l( r) · δ E_m( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r)   E_m( r) ·  E_l, r^*_0.4(1)( r)   d𝒱;                                                                    - iωε__0∫_𝒱 δχ( r)   E_m( r) · δ E_l, r^*_0.4(1)( r)   d𝒱; ⇔       o__lm(||δχ||) = ℱ__lm^^(2)(δχ,δχ) +  o__lm(||δχ||^^2) ,                                          (in a compact form) ].


With:


    {[ ℱ__lm^^(2)(δχ,δχ) = - iωε__0∫_𝒱 δχ( r)   E_m( r) ·  E_l, r^*_0.4(1)( r)   d𝒱;                   Where ℱ__lm^^(2) is the second order functional derivative ].



  ∙ Third order functional derivative To find the expression of this derivative, we repeat the same calculus scheme, i.e by writing : o__lm(||δχ||^^2) = ℱ__lm^^(3)(δχ,δχ,δχ) +  o__lm(||δχ||^^3). So, we start from the expression of the remainder o__lm(||δχ||^^2): 


    o__lm(||δχ||^^2) = - iωε__0∫_𝒱 δχ( r)   E_m( r) · δ E_l, r^*_0.4(1)( r)   d𝒱 = ∫_𝒱 δ E_l, r^*_0.4(1)( r) · [- iωε__0 δχ( r)  E_m( r)   d𝒱]

 
By following the same procedure one obtains:


    {[   - iωε__0∫_𝒱 δχ( r)   E_m( r) · δ E_l, r^*_0.4(1)( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(1)( r) ·  E_m, r^*_0.4(1)( r)   d𝒱;                                                              - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(1)( r) · δ E_m, r^*_0.4(1)( r)   d𝒱; ⇔       o__lm(||δχ||^^2) = ℱ__lm^^(3)(δχ,δχ,δχ) +  o__lm(||δχ||^^3) ,                                          (in a compact form) ].


With:


    {[ ℱ__lm^^(3)(δχ,δχ,δχ) = - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(1)( r) ·  E_m, r^*_0.4(1)( r)   d𝒱;                                         Where ℱ__lm^^(3) is the third functional derivative ].





 §.§ Functional derivative of order (n)
 The previous results can be generalized to the calculus of the functional derivative for any order n. We treat the two cases where n is even n=2p or odd n=2p+1, with p⩾ 2 a natural number. 



  ∙ Functional derivative of order (2p) Its expression is given by:


    {[  - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p-1)( r) · δ E_m, r^*_0.4(p-1)( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(p-1)( r) ·  E_l, r^*_0.4(p)( r)   d𝒱;                                                                             - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(p-1)( r) · δ E_l, r^*_0.4(p)( r)   d𝒱; ⇔       o__lm(||δχ||^^2p-1) = ℱ__lm^^(2p)(δχ, ... ,δχ_2p times) +  o__lm(||δχ||^^2p) ,                                          (in a compact form) ].


With:


    {[ ℱ__lm^^(2p)(δχ, ... ,δχ_2p times) = - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(p-1)( r) ·  E_l, r^*_0.4(p)( r)   d𝒱;                                               Where ℱ__lm^^(2p) is the functional derivative of order (2p) ].



  ∙ Functional derivative of order (2p+1) It is written in the form: 


    {[          - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(p-1)( r) · δ E_l, r^*_0.4(p)( r)   d𝒱 = - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p)( r) ·  E_m, r^*_0.4(p)( r)   d𝒱;                                                                                   - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p)( r) · δ E_m, r^*_0.4(p)( r)   d𝒱; ⇔       o__lm(||δχ||^^2p) = ℱ__lm^^(2p+1)(δχ, ... ,δχ_2p+1 times) +  o__lm(||δχ||^^2p+1) ,                                          (in a compact form) ].


With:


    {[ ℱ__lm^^(2p+1)(δχ, ... ,δχ_2p+1 times) = - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p)( r) ·  E_m, r^*_0.4(p)( r)   d𝒱;                                             Where ℱ__lm^^(2p+1) is the functional derivative of order (2p+1) ].




The total electric field E_w, r^*_0.4(p)( r) (with w=l,m) is obtained after solving p+1 forward problems. This resolution is illustrated by the following diagram:



The next section as well as appendix <ref> describe how the volume integral formulation is used, for the numerical evaluation of the electric fields involved in the different steps of the above diagram.



§ APPLICATION TO A FORWARD MODEL




 §.§ Study Configuration and Incident Field







  ∙ Configuration of the Study In a fixed Cartesian coordinate system ( O,  x__1,  x__2,  x__3), we consider a dielectric cylinder with circular cross section shape and with O x__2 as its invariance axis Fig (<ref>). This object is described by its permittivity-contrast function χ=χ(x__1, x__3), defined in the domain of volume (𝒱). The latter is discretized with Q=3436 square pixels of side Δ≃1.5  mm and area Δ𝒱≃2.25  mm^_2. We are dealing with the case of TE polarization (electric component parallel to the invariance axis), and we model the electric component of the field, which is parallel to the invariance axis. The total field modelled in the object, and the scattered field simulated on the receivers, are respectively written: E= E(x__1, x__3)  x__2 and E^ d= E^ d(x__1, x__3)  x__2. A harmonic time dependence of the form e^- iω t with pulsation ω is assumed. In this study, we consider M antennas arranged along a circle (Γ) of radius 11  cm, and concentric with the object of radius 5  cm. All these antennas are situated at equal distance from the center of the object. Each one of them plays the role of emitter and receiver but not simultaneously. When an antenna l acts like an emitter of source current J_l (l∈{1, 2,...,  M}), all others act as receivers to detect the scattered field E^ d_m, m∈{1, 2,...,  M}_m≠ l. We work with a number of antennas M=16, for a total of generated data equal to 240.

  ∙ Incident Field A source current J_l of unit amplitude J_l= J_l x__2=1  x__2, positioned at point r_l=x^_l__1 x__1+x^_l__3 x__3 and oriented along the invariance axis O x__2, radiates in two-dimensional free space the following incident field: 

    E^ i_l( r)  x__2 = -1/4 ωμ__0  H^_(1)__0( k | r- r_l|)  x__2


Here, E^ i_l is obtained by solving the scalar Helmholtz equation <cit.>, and k=2π/λ denotes the wavevector modulus of the incident medium. H^_(1)__0 is the first-kind Hankel function of zero-order. The wavelength of the incident field is set to λ=10  cm (i.e. 3 GHz of frequency).  





 §.§ Scattering Operator

All vectorial physical quantities are aligned along the invariance axis. Therefore, the forward problem is recast as a scalar problem which is solved in the plane of incidence. This is done by the volume integral formulation, which builds a scattering operator. The latter consists of the following two coupled integral equations:



  ∙ State Equation Its solution is the total field E_l. It is the result of an interaction between the incident field E^ i_l produced by the source current J_l and the object which is described by its permittivity-contrast function χ.

    E_l( r) =  E^ i_l( r) + ( i k^_2/4)∫_𝒱 H^_(1)__0( k | r- r'|) χ( r')  E_l( r')   d𝒱' ,          r∈𝒱


  ∙ Observation Equation Once the total field E_l is computed using the state equation, it is replaced in the observation equation to model the diffracted field E_l^ d on the receiver, located at r_m.

    E^ d_l( r_m) = ( i k^_2/4)∫_𝒱 H^_(1)__0( k | r_m- r'|) χ( r')  E_l( r')   d𝒱' ,                         r_m∈Γ  and   r_m≠ r_l


  ∙ Symbolic Notation for the Scattering Operator and the Data By varying the position r_l of the source current emitter J_l, l∈{1, 2,..., M}, and the position r_m, m∈{1, 2,..., M}_m≠ l of the receiver, we generate data of the diffracted field, denoted: 𝒟={𝒟_lm= E^ d_l( r_m) ·  J_m, for: l,m∈{1, 2,...,  M}_l≠ m}. The above equations Eqs. (<ref>,<ref>) give a nonlinear link between the data and the permittivity-contrast function, which from a mathematical point of view plays the role of a nonlinear scattering operator 𝒩, such that: 𝒟=𝒩(χ). Its numerical implementation is performed with the method of moments as described in references <cit.>.
 



 §.§ Numerical Study



  ∙ Limited Functional Expansion of the Data The functional derivatives expressed above can now be used to express the variation on the data δ𝒟={δ𝒟_lm=δ E^ d_l( r_m) ·  J_m, for : l,m∈{1, 2,...,  M}_l≠ m} as a limited functional expansion with respect to the permittivity-contrast function χ( r), for r∈𝒱. We are interested in the convergence of this development which is written in the form: 


    δ𝒟 = ∑_k=1^nℱ^^(k)(δχ, ... ,δχ_k fois) +  o(||δχ||^^n)


For odd (n), each element δ𝒟_lm of δ𝒟 is given by its expression in a condensed form:


    {[                                                                                       δ𝒟_lm = δ E^ d_l( r_m) ·  J_m = δ E^ d_l( r_m);                                                                                                                                     ; = ℱ__lm^^(1)(δχ) + ℱ__lm^^(2)(δχ,δχ) + ℱ__lm^^(3)(δχ,δχ,δχ) + ℱ__lm^^(4)(δχ, ... ,δχ_4 times) + ℱ__lm^^(5)(δχ, ... ,δχ_5 times) +...;                                                                                                                                     ;                                    +ℱ__lm^^(2p)(δχ, ... ,δχ_2p times) + ℱ__lm^^(2p+1)(δχ, ... ,δχ_2p+1 times) +  o__lm(||δχ||^^2p+1) ].

The explicit form of the above expansion, and the details of the numerical evaluation of the electric fields involved in its terms are described in appendix <ref>.


  ∙ Convergence Criteria To monitor the convergence of development Eq (<ref>), we introduce the residual error ξ and study its behaviour as a function of the truncation order (n):


    {[                                The residual error on data is the misfit: ξ(n) = Δ𝒟-∑_k=1^nℱ^^(k)(Δχ, ... ,Δχ_k times);                                                                                                                      ;   ℒ__2 norm of ξ expressed as a percentage:  P__ξ(n) = 100 × √(∑_lm=1^ M^^2-M|ξ__lm(n)|^^2/∑_lm=1^ M^^2-M|Δ𝒟__lm|^^2);                                                                                                                      ; ℒ__2 norm of ξ expressed in decibels:  L__ξ(n) = 10 × log__10 (∑_lm=1^ M^^2-M|ξ__lm(n)|^^2/∑_lm=1^ M^^2-M|Δ𝒟__lm|^^2) ].


Note that the exact calculus of the difference on the data is directly obtained with the forward model, namely : Δ𝒟=𝒩(χ+Δχ)-𝒩(χ). Concerning the variation on the contrast, it is numerically evaluated by the finite difference Δχ. 



  ∙ Numerical Results and Interpretation We take as a numerical application, an object of homogeneous dielectric contrast χ__ ref=6.0. At this reference value, we add a difference Δχ giving a new object. Then, this finite difference introduced to modify the initial contrast is estimated with the expression P__χ=100×(Δχ/χ__ ref). Figs (<ref>,<ref>,<ref>) illustrate the convergence in decibels L__ξ of the residual error ξ with respect to the truncation order n, and for different values of P__χ. For a better reading of these figures, we have added levels of value of this residual error, estimated with the expression of the percentage P__ξ. The latter are added to Figs (<ref>,<ref>,<ref>), and are represented by horizontal dashed lines. In Fig (<ref>), we consider the finite differences of the permittivity-contrast corresponding to the values P__χ=1%, 2%, 3%, 4% and 5%. For each case, we note that the convergence is reached for a truncation of order n, less than 15. Moreover, we note that the speed of this convergence decreases with the increase of P__χ. This tendency continues in Fig (<ref>), for finite differences of the contrast corresponding to the values P__χ=6%, 7% et 8%. Indeed, for these three values, the convergence is obtained for n<51. For a finite difference of value P__χ=9%, we note a slow decrease in the residual error L__ξ, suggesting that the convergence occurs for a truncation order n>51. This is precisely what happens when the order of truncation is carried up to n=161 (see Fig (<ref>)). Concerning the value P__χ=10%, the increase of error L__ξ according to the truncation order n shown in Fig (<ref>), is the mark that we are outside the domain of convergence of the functional limited expansion. These simulations were performed, with no particular optimisations, using MATLAB on a personal computer, with a quad-core processor at 2.60GHz and 16 Gb central memory. In particular, for a truncation order n=161, a simulation took 22 minutes of computation time. 







§ CONCLUSION


Previous studies have shown that use of an approach involving the reciprocity theorem, combined with a formulation of Maxwell's equations in the sense of distributions, is an efficient approach to calculate the Fréchet derivative. The latter stands for the functional derivative of the first order of the diffracted field, with respect to the permittivity-contrast function. The remainder, which contains all the derivatives of higher orders, has been neglected. In the present work, we have applied this principle of calculus on the remainder, and have succeeded in obtaining all the functional derivatives of higher orders. This allowed us to express, the variation of the data as a limited functional expansion with respect to the permittivity-contrast function. Finally, we have numerically illustrated the convergence of this development, in the case of a diffraction by a dielectric cylindrical object illuminated by an incident field in TE polarization. Such an expansion can be used to propose inversion methods in research areas involving the problem of permittivity-contrast reconstruction. We expect some improvements in the quality of the solutions as well as the speed of convergence.


Let's start with the equality:

    ∇ · [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] = - i ω (ε_ b - ε_ a)  E_ a,l ·  E_ b,m + δ_ r_m  E_ a,l ·  J_m - δ_ r_l  E_ b,m ·  J_l
 
 
Each term of the above equation being a distribution, we can apply it on a test function. Such a function is denoted u( r) and is chosen so as to satisfy the following conditions:



  ∙ u( r) is a function with bounded support.

  ∙ u( r) is a function of class C^∞.


We add the following additional conditions without loss of generality <cit.>: 


    {[   u( r) = 1   r ∈  C_ i         and   u( r) = 1      r ∈  C; ∇ u( r) ≠ 0   r ∈  F_ i         and   u( r) = 0      r ∈  E ].


Here, C denotes the surface of a sphere of volume C_ i, E the surface of a larger sphere including the first and having the volume E_ i. The term F_ i =  E_ i -  C_ i corresponds to the volume comprised between the surfaces C and E. In practice, it is possible to construct a test function that satisfies these three criteria. We proceed now on the application of each term of the Eq (<ref>) to the test function u( r). Starting with its left-hand side, setting T =  H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m, and applying then the distribution ∇ ·  T on the test function u, it follows:


    <∇ ·  T|u> = -< T|∇ u> = -∫_ F_ i  T · ∇ u   d F_ i = -∫_ F_ i ∇ · ( Tu)   d F_ i +∫_ F_ i (∇ ·  T) u   d F_ i


The integral ∫_ F_ i (∇ ·  T) u   d F_ i is nil because, according to Equation Eq (<ref>), each term of its right-hand side is equal to zero for all r∈ F_ i. On the other hand, the divergence theorem allows us to transform the volume integral -∫_ F_ i ∇ · ( Tu)   d F_ i into surface integral:


    <∇ ·  T|u> = -∫_ F_ i ∇ · ( T u)   d F_ i = +∫_ C  Tu ·  e__ρ   d C -∫_ E  Tu ·  e__ρ   d E


Here e__ρ denotes the unit vector orthogonal to the surfaces C and E. From the third property of the test function u, we have : u( r)=0 ∀  r∈ E and u( r)=1 ∀  r∈ C. This allows us to write the following result for the calculus of the left-hand side of equation Eq (<ref>): 


    <∇ ·  T|u> = +∫_ C  T ·  e__ρ   d C = +∫_ C [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] ·  e__ρ   d C


Let us now apply the right-hand side of equation Eq (<ref>) on the test function u. This gives directly:


    [ <- i ω (ε_ b - ε_ a)  E_ a,l ·  E_ b,m + δ_ r_m  E_ a,l ·  J_m - δ_ r_l  E_ b,m ·  J_l|u> =;                                                                                            ; - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱 +  E_ a,l( r_m) ·  J_m -  E_ b,m( r_l) ·  J_l ]

 
Now, in view of Eqs (<ref>,<ref>), we can rewrite Eq (<ref>) in the sense of functions, in the form:


    [ ∫_ C [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] ·  e__ρ   d C = - i ω∫_𝒱 (ε_ b - ε_ a)  E_ a,l ·  E_ b,m   d𝒱 +  E_ a,l( r_m) ·  J_m -;                                                                                                                 E_ b,m( r_l) ·  J_l ]





To show that the integral ∫_ C [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] ·  e__ρ   d C is nil, let C =  C_∞ be the surface of a sphere with a radius large enough to be in the farfield region. If we consider the configuration (a__l), its total magnetic field H_ a,l is given by the sum between the diffracted field, and the incident field radiated by a point-like source current located at a finite distance from the considered object. In the farfield region, this total field satisfies the Sommerfeld radiation condition lim___ρ→ +∞ρ [∂/∂ρ- ik]  H_ a,l=0, and is written as a product of a radial function by an angular function H_ a,l = e^ ikρ/4πρ  A_ a,l(θ,ϕ), with both expressed in spherical coordinates <cit.>. The magnetic farfields of the configurations (a__l) and (b__m) are given by the following expressions:


    {[ H_ a,l = e^ ikρ/4πρ  A_ a,l ∝ e^ ikρ/ρ  A_ a,l;                                               ; H_ b,m = e^ ikρ/4πρ  A_ b,m ∝ e^ ikρ/ρ  A_ b,m ].


Due to the relation E=1/ωε__0∇× H, the electric field is obtained from the magnetic field. After a few calculations, we have:


    {[ E_ a,l =  ik/4πωε__0 e^ ikρ/ρ  e__ρ ×  A_ a,l ∝ e^ ikρ/ρ  e__ρ ×  A_ a,l;                                                                         ; E_ b,m =  ik/4πωε__0 e^ ikρ/ρ  e__ρ ×  A_ b,m ∝ e^ ikρ/ρ  e__ρ ×  A_ b,m ].


Here e__ρ denotes the radial unit vector. In virtue of expressions Eqs (<ref>,<ref>), the following integral vanishes:


    {[                                      ∫_ C_∞ [ H_ b,m ×  E_ a,l -  H_ a,l ×  E_ b,m] ·  e__ρ   d C;                                                                                                  ; =e^ i2kρ/ρ^2∫_ C_∞ [{( A_ b,m ·  A_ a,l) ( e__ρ ·  e__ρ) - ( A_ b,m ·  e__ρ) ( A_ a,l ·  e__ρ)} -;               {( A_ a,l ·  A_ b,m)( e__ρ ·  e__ρ) - ( A_ a,l ·  e__ρ) ( A_ b,m ·  e__ρ)}]   d C=0 ].





For odd (n), each element δ𝒟_lm of δ𝒟 is given by the following detailed expression:


    [                                                                                                                                                                               δ𝒟_lm = δ E^ d_l( r_m) ·  J_m = δ E^ d_l( r_m);                                                                                                                                                                                                                             ;                                                                                        = - iωε__0∫_𝒱 δχ( r)   E_l( r)   E_m( r)  d𝒱_ℱ__lm^^(1)(δχ)  - iωε__0∫_𝒱 δχ( r)   E_m( r)   E_l, r^*_0.4(1)( r)  d𝒱_ℱ__lm^^(2)(δχ,δχ);                                                                                                                                                                                                                             ;                                  - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(1)( r)   E_m, r^*_0.4(1)( r)  d𝒱_ℱ__lm^^(3)(δχ,δχ,δχ)  - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(1)( r)   E_l, r^*_0.4(2)( r)  d𝒱_ℱ__lm^^(4)(δχ, ... ,δχ_4 times);                                                                                                                                                                                                                             ; - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(2)( r)   E_m, r^*_0.4(2)( r)  d𝒱_ℱ__lm^^(5)(δχ, ... ,δχ_5 times)  ... + ... + ...  - iωε__0∫_𝒱 δχ( r)   E_m, r^*_0.4(p-1)( r)   E_l, r^*_0.4(p)( r)   d𝒱_ℱ__lm^^(2p)(δχ, ... ,δχ_2p times);                                                                                                                                                                                                                             ;                         - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p)( r)   E_m, r^*_0.4(p)( r)   d𝒱_ℱ__lm^^(2p+1)(δχ, ... ,δχ_2p+1 times)  - iωε__0∫_𝒱 δχ( r)   E_l, r^*_0.4(p)( r)  δ E_m, r^*_0.4(p)( r)   d𝒱_ o__lm(||δχ||^^2p+1) ]


The following steps describe how to calculate the total fields E_w( r), E_w, r^*_0.4(1)( r), E_w, r^*_0.4(2)( r),..., E_w, r^*_0.4(p)( r) (with w=l,m):





  ∙ Step 0: In view of the Eq (<ref>), the source current J_w=1 located at r_w radiates the incident field E^ i_w( r) = -1/4 ωμ__0  H^_(1)__0( k | r- r_w|), (r∈𝒱). The resulting total field, solution of the state equation Eq (<ref>) is denoted by E_w( r), (r∈𝒱).



  ∙ Step 1: The domain 𝒱 is meshed with Q pixels of area Δ𝒱 centred on the positions r_q (with q=1,2,..., Q), the source currents - iωε_ o Δχ( r_q)  E_w( r_q) Δ𝒱 (with r_q spanning the domain 𝒱), radiate the total incident field E^ i_w, r^*_0.4(1)( r)= i( k^_2/4) Δ𝒱∑_q=1^ QΔχ( r_q)  E_w( r_q)  H^_(1)__0( k | r- r_q|), (r, r_q∈𝒱). Note that when r= r_q, the value of H^_(1)__0 is obtained by performing the calculus of the average (1/Δ𝒱)∫_Δ𝒱 H^_(1)__0( k | r- r_q|) d r, where Δ𝒱 is the area of the pixel q centred on the position r_q (details on this calculus are presented in <cit.>). As stated in section <ref>, the index r^*_0.4(1) has two significations: the symbol * means that we sum over all positions r_q (q=1,2,..., Q), in order to take into account the contributions of all the source currents in the total incident field. The index (1) means that the variation of the permittivity-contrast function δχ (replaced in this numerical study by the finite difference Δχ), is included once in the expression of the total incident field E^ i_w, r^*_0.4(1). Then, the total electric field E_w, r^*_0.4(1)( r), (r∈𝒱) is obtained by solving the state equation Eq (<ref>), with E^ i_w, r^*_0.4(1), (r∈𝒱) playing the role of incident field.



  ∙ Step 2: In this step we replace E_w by E_w, r^*_0.4(1), and we repeat the same calculus as in the previous step, namely: all source currents - iωε_ o Δχ( r_q)  E_w, r^*_0.4(1)( r_q) Δ𝒱 (with r_q spanning the domain 𝒱), radiate the total incident field E^ i_w, r^*_0.4(2)( r)= i( k^_2/4) Δ𝒱∑_q=1^ QΔχ( r_q)  E_w, r^*_0.4(1)( r_q)  H^_(1)__0( k | r- r_q|), (r, r_q∈𝒱). The index (2) means that the finite difference Δχ is included twice in the expression of the total incident field E^ i_w, r^*_0.4(2) (once explicitly and once implicitly because embedded in the field E_w, r^*_0.4(1)). Finally, the total field E_w, r^*_0.4(2)( r), (r∈𝒱) is obtained in the same way, that is to say by solving the state equation Eq (<ref>), with the term E^ i_w, r^*_0.4(2)( r), (r∈𝒱) as the incident field.

.........................

.........................

.........................


  ∙ Step p: The source currents - iωε_ o Δχ( r_q)  E_w, r^*_0.4(p-1)( r_q) Δ𝒱 (with r_q spanning the domain 𝒱), radiate the total incident field E^ i_w, r^*_0.4(p)( r)= i( k^_2/4) Δ𝒱∑_q=1^ QΔχ( r_q)  E_w, r^*_0.4(p-1)( r_q)  H^_(1)__0( k | r- r_q|), (r, r_q∈𝒱). The index (p) means that the finite difference Δχ is included p times in the expression of the total incident field E^ i_w, r^*_0.4(p) (once explicitly and p-1 times implicitly because embedded in the field E_w, r^*_0.4(p-1)). The total field E_w, r^*_0.4(p)( r), (r∈𝒱) is then obtained by solving the state equation Eq (<ref>), for which E^ i_w, r^*_0.4(p)( r), (r∈𝒱) is used as the incident field.



99

kawata1987optical Kawata, S., Nakamura, O. and Minami, S., “Optical microscope tomography. I. Support constraint," JOSA A, Vol. 4, No. 1, 292–297, 1987. 

nakamura1988optical Nakamura, O., Kawata, S. and Minami, S., “Optical microscope tomography. II. Nonnegative constraint by a gradient-projection method," JOSA A, Vol. 5, No. 4, 554–561, 1988.

cmielewski2007two Cmielewski, O., Tortel, H., Litman, A. and Saillard, M., “A two-step procedure for characterizing obstacles under a rough surface from bistatic measurements," IEEE Transactions on Geoscience and Remote Sensing, Vol. 45, No. 9, 2850–2858, 2007. 

gurbuz2014efficient Gurbuz, T. U., Aslanyurek, B., Karabulut, E. P. and Akduman, I., “An efficient nonlinear imaging approach for dielectric objects buried under a rough surface," IEEE Transactions on Geoscience and Remote Sensing, Vol. 52, No. 5, 3013–3022, 2014. 

hadamard1902 Hadamard, J., “Sur les problèmes aux dérivées partielles et leur signification physique [On partial derivative problems and their physical significance]," Princeton university bulletin 13, 1902.

tikhonov1977 Tikhonov, A. N. and Arsenin, V. Y., Solutions of ill-posed problems , W. H. Winston, Washington D. C., 1977.

kleinman1992modified Kleinman, R. E. and Van den Berg, P. M., “A modified gradient method for two-dimensional problems in tomography," Journal of Computational and Applied Mathematics, Vol. 42, No. 1, 17–35, 1992.

van1997contrast Van den Berg, P. M. and Kleinman, R. E., “A contrast source inversion method," Inverse problems, Vol. 13, No. 6, 1607–1620, 1997.

kantorovich1948newton Kantorovich, L. V., “On Newton's method for functional equations," Dokl. Akad. Nauk SSSR, Vol. 59, No. 7, 1237–1240, 1948.

roger1981newton Roger, A., “Newton-Kantorovitch algorithm applied to an electromagnetic inverse problem," IEEE Transactions on Antennas and Propagation, Vol. 29, No. 2, 232–238, 1981.

roger1980inverse Roger, A. and Maystre, D., “Inverse scattering method in electromagnetic optics: application to diffraction gratings," JOSA, Vol. 70, No. 12, 1483–1495, 1980.

cea1986conception Cea, J., “Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût [Optimal design or shape identification, fast calculation of the directional derivative of the cost function.]," ESAIM: Mathematical Modelling and Numerical Analysis, Vol. 20, No. 3, 371–402, 1986.

litman1998reconstruction Litman, A., Lesselier, D. and Santosa, F., “Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set," Inverse problems, Vol. 14, No. 3, 685-706, 1998.

tarantola1984 Tarantola, A., “Inversion of seismic reflection data in the acoustic approximation," Geophysics, vol. 49, no. 8, 1259–1266, 1984.

tarantola1987 Tarantola, A., Inversion of travel times and seismic waveforms , Springer, Dordrecht, 1987.

tarantola1988 Tarantola, A., “Theoretical background for the inversion of seismic waveforms including elasticity and attenuation," Pure and Applied Geophysics, vol. 128, no. 1, 365–399, 1988.

tromp2008 Tromp, J., Komatitsch, D. and Liu, Q., “Spectral-element and adjoint methods in seismology," Communications in Computational Physics, vol. 3, no. 1, 1–32, 2008.

monteiller2015 Monteiller, V., Chevrot, S., Komatitsch, D. et Wang, Y., “Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM–DSM hybrid method," Geophysical Journal International, vol. 202, no. 2, 811–827, 2015.

roger1982 Roger, A., “Reciprocity theorem applied to the computation of functional derivatives of the scattering matrix," Electromagnetics, vol. 2, no. 1, 69–83, 1982.

arhab2016 Arhab, S. and Soriano, G., “Inverse Wave Scattering of Rough Surfaces with Emitters and Receivers in the Transition Zone," Progress In Electromagnetics Research M, vol. 45, 131–141, 2016.

arhab2011 Arhab, S., Soriano, G., Belkebir, K., Sentenac, A. and Giovannini, H., “Full wave optical profilometry," Journal of Optical Society of America A, vol. 28, no. 4, 576–580, 2011.

arhab2012 Arhab, S., Giovannini, H., Belkebir, K. and Soriano, G., “Full polarization optical profilometry," Journal of Optical Society of America A, vol. 29, no. 8, 1508–1515, 2012.

el2009adjoint El-Shenawee, M., Dorn, O. and Moscoso, M., “An adjoint-field technique for shape reconstruction of 3-D penetrable object immersed in lossy medium," IEEE Transactions on Antennas and Propagation, vol. 57, no. 2, 520–534, 2009.

yaghjian1986overview Yaghjian, A. D., “An Overview of Near-Field Antenna Measurements," IEEE Transactions on Antennas and Propagation, vol. 34, no. 1, 30–45, 1986.

appel2007 Appel, W. and Kowalski, E., Mathematics for physics and physicists , Princeton University Press Princeton, NJ, USA; Oxford, UK, 2007.

rudin1991 Rudin, W., Functional analysis , International series in pure and applied mathematics, McGraw-Hill, Inc., New York, 1991.

abramowitz2012handbook Abramowitz, M. and Stegun, I. A., Handbook of mathematical functions: with formulas, graphs, and mathematical tables , Courier Corporation, 2012.

richmond1965scattering Richmond, J., “Scattering by a dielectric cylinder of arbitrary cross section shape," IEEE Transactions on Antennas and Propagation, vol. 13, no. 3, 334–341, 1965.

lesselier1991buried Lesselier, D. and Duchêne, B., “Buried, 2-D penetrable objects illuminated by line sources: FFT-based iterative computations of the anomalous field," Progress In Electromagnetics Research, vol. 5, 351–389, 1991.

tsang2001scattering Tsang, L., Kong, J. A., Ding, K. H. and Ao, C. O., Scattering of electromagnetic waves: Vol. II. Numerical simulations , Wiley Series on Remote Sensing, Wiley, 2001.

tsang2000scattering Tsang, L., Kong, J. A. and Ding, K. H., Scattering of electromagnetic waves: Vol. I. Theories and applications , Wiley Series on Remote Sensing, Wiley, 2000.



