Application of the second domain derivative in inverse electromagnetic scattering

We consider the inverse scattering problem of reconstructing a perfect conductor from the far field pattern of a scattered time harmonic electromagnetic wave generated by one incident plane wave. In order to apply iterative regularization schemes for the severely ill-posed problem the first and the second domain derivative of the far field pattern with respect to variations of the domain are established. Characterizations of the derivatives by boundary value problems allow for an application of second degree regularization methods to the inverse problem. A numerical implementation based on integral equations is presented and its performance is illustrated by a selection of examples.


Introduction
A challenging class of inverse problems in scattering theory is the identification of scattering objects by the knowledge of far field patterns of scattered waves (see [3]). We have to distinguish the inverse scattering problems, if the response to any or at least to many incident fields is known, or if data are given only for one or a few incident fields. In this work we are going to consider the extreme situation of the reconstruction of the shape of a perfect conductor just from the knowledge of the far field pattern of one scattered time harmonic electromagnetic wave. * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Derivative based iterative regularization schemes are known to be suitable numerical approaches for this class of problems (see the discussion in section 5.4 in [3] and references cited therein). Thus, we focus on linearization of the far field pattern with respect to variations of the shape of the scattering object. The derivative is given by the far field pattern of the so called domain derivative of the scattered wave. These domain derivatives are well established for most of the usually considered boundary value problems (see [14] and references cited therein). Furthermore, in case of acoustic scattering problems several numerical implementations are documented (see [6,11] for three dimensional examples). Presumably due to the computational effort, there are only a few results for the full vector valued electromagnetic inverse scattering problem (see [9,15,21]). These approaches are based on boundary integral equations for the electromagnetic scattering problem and the first domain derivative. We will extend on these results, mainly by showing the existence and a characterization of the second domain derivative. This gives rise to an application of second degree regularization schemes (see [13,16]).
After this introduction, in section 2, we collect some notations for the scattering problem and describe its weak formulation for later use. In section 3 the first domain derivative of the scattering problem is discussed. Although the derivative has already been established (see [19,24]), we will present it in some detail in preparation of the following investigations for the second derivative. Moreover, the characterization of the domain derivative is used to show a condition, which ensures that a constant variation cannot lead to a vanishing domain derivative. This illuminates the challenging question on injectivity of the domain derivative operator. With these preparations we devote section 4 to the second domain derivative. It is shown that such a derivative exists and can be characterized again by an electromagnetic boundary value problem. Finally, based on these characterizations we explain and discuss in section 5 the regularized Halley method applied to the inverse problem and present its numerical performance by some examples.
Some results of this paper, e.g. theorem 4.5, are part of one of the authors' PhD thesis [8].

The scattering problem
Let us assume a bounded scattering obstacle D ⊂ R 3 with smooth boundary and simply connected complement R 3 \D. The object is surrounded by a homogeneous, linear, isotropic and non-conductive medium with electric permittivity ε 0 > 0 and magnetic permeability μ 0 > 0, for instance vacuum. At frequency ω > 0, the time harmonic Maxwell system for the electric field E and the magnetic field H then reads as curl E − ikH = 0, curl H + ikE = 0, (2.1) with wave number k = ω √ ε 0 μ 0 (see [3]). Given an incident plane wave, E i (x) = pe ikd·x , H i (x) = (d × p)e ikd·x for x ∈ R 3 with polarisation p ∈ C 3 and direction d ∈ S 2 satisfying p · d = 0 the scatterer gives rise to a radiating scattered field (E s , H s ), a solution of the Maxwell system (2.1) in R 3 \D, which satisfies the Silver-Müller radiation condition The interaction of the perfect conductor D with the incident wave can be formulated as a boundary condition for the total field E = E s + E i and is given by where ν denotes the outwards directed normal vector to ∂D.
The following investigations require a variational formulation of the scattering problem. Thus, we choose R > 0 large enough such that D ⊂ B R (0), where B R (0) denotes the open ball of radius R centered in the origin, and introduce the bounded domain Ω = B R (0)\D. In order to derive the weak formulation, let (E, H) be a pair of reasonable smooth solutions of the scattering problem and let V denote a test function with ν × V = 0 an ∂D. By partial integration and the Maxwell system (2.1) we arrive at (2.2) To ensure that a solution E s = E − E i of (2.2) can be extended to a radiating solution of the Maxwell system in R 3 \D we have introduced on the artificial boundary ∂B R (0) the Calderon operator Λ, which maps ν × ϕ onto ν × H s , where (E s , H s ) denote the unique radiating solution of Equation (2.2) is considered in the Sobolev space H(curl, Ω) = {E ∈ L 2 (Ω, C 3 ) : curl E ∈ L 2 (Ω, C 3 )}. The boundary integrals on ∂D and ∂B R (0) exists in the sense of the dual pairing ·, · ∂B R (0) between the range spaces H − 1 2 (Div, ∂Ω) and and we incorporate the boundary condition by the closed subspace Since, the Calderon operator is extendable to a bounded operator Λ : H − 1 2 (Div, ∂B R (0)) → H − 1 2 (Div, ∂B R (0)) (see [22]), we finally can define the bounded sesquilinear form A : H pc (Ω) × H pc (Ω) → C and the antilinear map : H pc (Ω) → C such that (2.2) reads as A weak solution of the scattering problem is then given by a function E ∈ H pc (Ω) such that A(E, V) = (V) holds for all V ∈ H pc (Ω). Assuming D to be a Lipschitz domain it is known that for any ∈ H pc (Ω) * there exists a unique solution E ∈ H pc (Ω) of (2.3) for all V ∈ H pc (Ω), and it exists c > 0 such that E H(curl,Ω) c H * pc (see [22, theorem 10.7]). Due to the radiation condition, the scattered field E s in R 3 \D has the asymptotic behavior E ∞ is called the (electric) far field pattern and is an analytic tangential vector field on the unit sphere S 2 . This motivates the definition of the non-linear boundary to far field operator F, which maps the boundary ∂D onto the far field pattern of E s , i.e., Of course, F depends also on the incident field (E i , H i ) and the wave number k, which we assume to be fixed and known. The domain of F is given by a class of admissible boundaries, for which there is a unique solution of the scattering problem. Thus, the inverse obstacle problem under consideration is given by the inversion of equation (2.4). It is known that the far field pattern uniquely determines the solution of the scattering problem, but, nevertheless, the whole inverse obstacle problem is severely ill-posed. For some more details on inverse electromagnetic scattering we refer to [3], where, for instance, uniqueness of the inverse problem is shown in the sense that if for a fixed wave number k the far field patterns of two objects for all incident plane waves coincide, the scattering objects must be identical. Such a result is not known in case of just one incident field.

Linearization of the inverse problem
In solving for the nonlinear equation (2.4), obviously a linearization is useful. Thus, a derivative of the far field pattern with respect to variations of the boundary of the scattering obstacle D is of specific interest. This leads to the concept of a domain derivative, which is well established in electromagnetic scattering (see [4,12,14,19,24]). For convenience to the reader and in preparation of the next section we present the variational approach in some detail, following very closely [12], where penetrable scattering objects are considered.
A perturbation of the scatterer is described by a vector field h ∈ C 1 (R 3 , R 3 ) with compact support. Given a set D ⊂ R 3 , we denote by D h the corresponding perturbed set is a diffeomorphism. Throughout, we assume that a perturbation h does not change the artificial boundary ∂B R (0). Thus, without loss of generality we have h ∈ Let E ∈ H pc (Ω) be the weak solution of the scattering problem (2.2) and let E h ∈ H pc (Ω h ) denote the solution of the scattering problem with respect to a perturbed scatterer D h , i.e.
for all V h ∈ H pc (Ω h ). Note the same right-hand sides of the weak formulations (2.2) and (3.1), since the boundary integral on the artificial boundary does not change. According to different sets of definition, E h has to be transformed. We use the curl conserving transformation E h → E h , given by where J ϕ denotes the Jacobian matrix of ϕ (see [2, section 3.9]). Then where curl ∼ denotes the curl operator with respect to the untransformed coordinates, and from By the transformation we arrive for for all V ∈ H pc (Ω). We define the bounded sesquilinear form Straight forward calculations yield the asymptotic behavior for h → 0. As a first step we show continuity with respect to the perturbation h.
We calculate and by (3.3) and (3.4) we obtain Looking closely at the linearizations of the coefficients in the weak formulation (3.2), we can prove differentiability. Theorem 3.2. Let E ∈ H pc (Ω) be the solution of (2.2) and E h ∈ H pc (Ω) of (3.2). Then there exists a function W ∈ H pc (Ω), depending linearly on h ∈ C 1 0 (B R (0), R 3 ), such that Proof. We define W ∈ H pc (Ω) as the solution of We introduce the notation V ν = V · ν and V τ = ν × (V × ν) for the decomposition V = V ν ν + V τ of a vector field on ∂D. Furthermore, we introduce the surface gradient Grad ∂D : H Curl ∂D u = Grad ∂D u × ν, and the surface divergence Div ∂D : H − 1 2 (Div, ∂D) → H − 1 2 (∂D), which is defined for a smooth function V by Div ∂D V = divV − ν · J V ν. Note that the surface divergence satisfies Div ∂D (V × ν) = curl V · ν (3.5) and is coupled by the duality With these notations a representation of the domain derivative can be shown.
Proof. By the regularity of ∂D we have E ∈ H 1 (Ω, C 3 ) and therefore E = W − J h E − J E h ∈ L 2 (Ω, C 3 ) (see [1]). Some basic vector calculus shows (3.6) which in particular implies curl E ∈ L 2 (Ω, C 3 ) and therefore E ∈ H(curl, Ω). Additionally, by ν × W = 0 on ∂D, we find and H ν = 0, which follows by Maxwell's equations and (3.5), we conclude the stated boundary condition for E .
It remains to show, that E is a radiating solution to Maxwell's equations, which will be achieved by showing A(E , V) = 0 for any V ∈ H pc (Ω). We have Using again (3.6), we find Since divcurl = 0, we obtain by the divergence theorem Note that no boundary integrals on ∂B R (0) occur, since h is compactly supported in B R (0). The first term vanishes since ν curl V = Div ∂D (V × ν) = 0 for V ∈ H pc (Ω) and for the second term we compute Thus, we have A(E , V) = 0, which finishes the proof.
Before establishing the second domain derivative in the next section let us consider the linearization of the operator F, which by the previous result is given by its Fréchet derivative F [∂D]h = E ∞ , if we specify a linear space of admissible boundaries. Note, that the domain derivative E of E in R 3 \D also allows for a consideration of near field data instead of the far field pattern. In general solving an ill-posed nonlinear equation by iterative regularization schemes based on its derivative requires some additional conditions on the operator F. A quite general one is the tangential cone condition, which can be described by the existence of a con- holds. It ensures to some extent the equivalence of local ill-posedness of a nonlinear equation and illposedness of its linearization (see [17]). To our knowledge the validity of such a condition is an essential open problem in any inverse obstacle scattering problem so far. Here, we just remark on injectivity of F , a necessary consequence from the cone condition, which is a severe problem in itself.

Corollary 3.4. Let the wave number −k 2 be no eigenvalue of the Laplace-Beltrami operator
Assuming h ν = 0, a rotation and taking the surface divergence yields Since −k 2 is no eigenvalue of the Laplace-Beltrami operator, we obtain E ν = 0 on ∂D. Furthermore, by (3.7) we have ν × (H × ν) = 0 on ∂D. Applying the Stratton-Chu representation of E i in D and of the radiating solution E s in R 3 \D (see [18]) we obtain from vanishing boundary values ν × E = 0 and ν × (H × ν) = 0 on ∂D of the total field the contradiction for any x ∈ D. Thus we conclude h ν = 0.
Excluding eigenvalues of the Laplace-Beltrami operator seems to be necessary for injectivity of F . This can be seen from scattering by a ball B ρ (0) of radius ρ > 0, as it was already observed by H Haddar and R Kress in [7]. Since, if we consider an incident field n for a positive integer n ∈ N and m ∈ {−n, . . . , n} (see [18]), some calculations show for the total field for n ∈ N, i.e., if −k 2 is an eigenvalue of the Laplace-Beltrami operator on ∂B ρ (0), there are nontrivial incident fields leading to a vanishing boundary condition (3.7) and therefore to a vanishing domain derivative.
Based on theorem 3.3 we can characterize the adjoint operator (F ) * , which is of specific interest for iterative regularization schemes, and adds to comparable results for the exterior Dirichlet problem in acoustic scattering (see [10]). We specify the investigations to the case of star shaped domains, which also will be considered in the numerical tests below. Without changing notation we consider the operator F : r → E ∞ for r ∈ C 2 (S 2 ) and a star shaped parameterized boundary ∂D = {y = r(ŷ)ŷ :ŷ ∈ S 2 }. Analogously, a variation is given by h(y) =h(ŷ)ŷ.

where E A , H A denote the total fields of the scattering problem with incident field given by the Herglotz wave function
Proof. We introduce the notation E s (y; a,x) for a solution of the scattering problem with incident field E i (y; a,x) = a e ikx·y and H i (y; a,x) =x × E i (y; a,x) with directionx and polarization a⊥x. From the integral representation of the far field pattern of the radiating solution E (see [3]) we compute for anyx ∈ S 2 andx⊥a ∈ C 3 . Applying Green's vector formula in the exterior of D together with the Silver-Müller radiation conditions implies Using this representation with polarizations given by a tangential field A ∈ L 2 t (S 2 , C 3 ) and substituting the boundary condition from theorem 3.3 yields where we have denoted as before the solution of the actual scattering problem by E = E(·; p, d) and H = H(·; p, d). If ∂D is star shaped, parameterized by y = r(ŷ)ŷ, and the variation is given by h(y) =h(ŷ)ŷ we can calculate h ν explicitly and arrive at which shows the assertion.

The second domain derivative
We continue in proving the scattered wave to be twice differentiable with respect to the boundary. If we use two small perturbations h 1 , h 2 ∈ C 1 (R 3 , R 3 ) with compact support in B R (0) to perturb the boundary, we arrive at which is not symmetric with respect to the variations h 1 and h 2 . But we expect a second derivative to be symmetric, see [5,chapter VIII.12]. The perturbation becomes symmetric, if we replace h 1 by h 1 • ϕ −1 2 . This motivates our goal: finding a radiating solution of Maxwell's equations E , depending bilinearly on h 1 and h 2 , being symmetric in h 1 and h 2 , and satisfying denotes the domain derivative with respect to the variation h at the scatterer D.
with (E 1 ) 2 being the domain derivative with respect to the variation h 2 of the domain derivative with respect to the variation h 1 . The second term E h is the domain derivative from theorem 3.3 with respect to the variation h = J h 1 h 2 .
We are going to prove that the second domain derivative is given by (4.1) and present a characterization of E as a radiating solution to Maxwell's equations. Similar to the first derivative, we start by showing existence of the material derivative of the material derivative. Let where we introduced the abbreviation A i for the symmetric matrix for all V ∈ H pc (Ω). In the next lemma, we provide the linearization of the new matrices.

Proof. The linearizations follow from (3.3), (3.4) and the Taylor expansion of the coefficients
As a first step we prove that the material derivative W 1 depends continuously on perturbations h 2 .
Proof. Let h 2 (V) denote the right-hand side of (4.2) with i = 1 and let h 2 ,h 1 denote the righthand side of (4.3). Recall the notation A h 2 for the sesquilinear form, such that the left-hand side of (4.3) is given by A h 2 ( W 1,h 2 , V). Then we have Adding and subtracting the integral With lemma 4.1 and theorem 3.
As before, we consider the linearizations and prove differentiability.
. Then there exists a function W 1 ∈ H pc (Ω), depending linearly and continuously on h 2 ∈ C 1 , such that Proof. Motivated by (4.4) we define W 1 ∈ H pc (Ω) as the solution of As before, we consider the difference . We add and subtract the following integrals i.e. we consider This leads to the estimate for all V ∈ H(curl, Ω). Again by a perturbation argument, we conclude Since W 1 ∈ H pc (Ω) is the material derivative with respect to h 2 of the material derivative with respect to h 1 , it contains by linearity the domain derivative with respect to h 2 of the domain derivative with respect to h 1 , which we denoted by (E 1 ) 2 before. To calculate it, we consider the formal Taylor expansion With the decomposition and, furthermore, with is motivated. Similarly to the first domain derivative, we need higher regularity of the solution and therefore higher regularity of the boundary, to ensure the Ansatz to be well defined.
Proof. See appendix A.
In order to give a characterization of the second domain derivative E = (E 1 ) 2 − E h with h = J h 1 h 2 , we need to introduce the symmetric curvature operator R : ∂D → R 3×3 , which acts on the tangential plane and is given by R(x) = J ν (x), x ∈ ∂D. Furthermore we define the mean curvature κ : ∂D → R by κ = 1 2 div ν. Note, that these definitions require differentiable extensions of the normal vector field ν in a neighborhood of ∂D which is constant in the direction of ν, see [23]. We state the main result of this paper.
since R is acting on the tangential plane. With W 1 ∈ H pc (Ω) the boundary values of (E 1 ) 2 are We use the decomposition of the material derivative W i = E i + J h i + J E h i , for i = 1, 2 to find similarly as before As seen before, we have From the boundary condition ν × E = 0 on ∂D we conclude We gather some identities, Combining and substituting these into ν × E = ν × (E 1 ) 2 + ν × E h yields the boundary condition For any vector field F, we have on the boundary ∂D (see (5.4.50) in [23]), and For the tangential part of the curl we get (see theorem 2.5.20 in [23]). With equation (4.5) we conclude Furthermore, by div(H ) = 0 and H ν = 0, we have With equation (4.7) we obtain Thus, we arrive at By (4.6) we obtain

Furthermore, it holds
By the product rule we finally arrive at a symmetric characterization, i.e., From (4.5) and from (4.6) we see ∂H τ ∂ν = −RH τ , and ∂E ν ∂ν = −2κE ν and conclude We use the boundary condition as stated in the theorem.
We do not claim that the characterization is the most elegant or shortest way to describe the boundary condition of theorem 4.4. But it shows its symmetry with respect to h 1 and h 2 . Note that the boundary condition of the second domain derivatives requires both the solution (E, H) as well as the first domain derivatives (E i , H i ) to be sufficiently smooth in order to be well posed.

The regularized Halley method
Recall the boundary to far field operator defined by F(∂D) = E ∞ ∈ L 2 t (S 2 ). From the previous section, we know F to be twice differentiable where the derivatives are given by the far field patterns of the domain derivatives, i.e.
t (S 2 ), we apply Newton type methods. Choosing a starting guess ∂D 0 a classical Newton step consists in solving the linear equation Due to ill-posedness the equation (5.1) has to be regularized in order to ensure solvability. Applying Tikhonov regularization, we consider the uniquely solvable equation with some chosen regularization parameter α 1 > 0. Thus the regularized Newton scheme uses an update of the boundary by ∂D i+1 = ∂D i h . For more details on iterative regularization methods we refer to [17].
For the second degree method, we modify this approach and use h just as a predictor to linearize the quadratic approximation F(∂D h ) ≈ F(∂D) + F again with some regularization parameter α 2 > 0. Then the so called Halley method is given by an update of the boundary by ∂D i+1 = ∂D i h . To obtain a regularization scheme for the full non-linear problem it is known that we have to add a stopping condition. Therefore, we stop the iteration if the relative residual F(∂D i ) − E ∞ L 2 (S 2 ) / E ∞ L 2 (S 2 ) falls below a chosen threshold, which depends on the noise level. This completes the regularized second degree method, which is sometimes called regularized Halley method. The method is introduced in [13,16], where regularizing properties are shown providing assumptions similar to those used for the regularized Newton method. Mainly, the tangential cone condition, already mentioned in section 3, is required.
However, we consider a numerical implementation of the method, which requires the choice of a set Y of admissible boundaries as an open set of a normed space X . Then, the domain derivatives become Frechet derivatives. We have chosen Y to be the set of star shaped domains with center in the origin and boundary of class C ∞ , discretized in the same way as in [8,9] by spherical harmonics Y m n in the following way: first, we identify the boundary ∂D ∈ Y by the positive smooth function r : S 2 → R, such every x ∈ ∂D is given in spherical coordinates by x = r(d)d for some d ∈ S 2 , i.e. we choose the open set Y = {r ∈ C ∞ (S 2 ) : r > 0} in the space X = C ∞ (S 2 ) as the domain of F. To discretize Y, we choose the finite dimensional subspace X N ⊂ X , using the real and imaginary part of spherical harmonics Y m n up to the degree N ∈ N, i.e., we arrive at the representation matrix for the discretized operator The discretization of the identity I is given by the identity matrix I (N+1) 2 . We observed a better performance of our scheme by using a different penalty matrix J instead, which punishes the curvature of the boundary. Such a matrix J is for example given by the diagonal matrix with entries (J) kk = 1 + λ(k), k = 1, . . . , (N + 1) 2 . Here, λ(k) is the corresponding eigenvalue of the spherical harmonic Y m n with respect to the Laplace Beltrami operator Δ S 2 = Div S 2 Grad S 2 , associated to the kth basis element of B. For the predictor h = (α h , β h ) ∈ R (N+1) 2 we solve the discretized version of equation (5.2) with I replaced by J. In general, a solution of this equation is complex-valued, so we discard the imaginary part.
Let E ∞ (x; h, h) denote the far field pattern of the second domain derivative with respect to the perturbations h 1 = h and h 2 = h, evaluated atx ∈ S 2 . Then the representation matrix for the discretized operator

Full discretization requires the numerical evaluation of F(∂D), F [∂D] and F [∂D]
. Looking closely at the boundary conditions for the first and second domain derivative, we identify the traces of the solutions (E, H) and (E , H ) and some terms involving surface differential operators to these traces. We therefore chose an integral equation approach for the full discretization. Our implementations were carried out in the open source Galerkin boundary element methods library BEMPP (https://bempp.com). For an overview of the library, see [26]. We will briefly present the tools needed to formulate the scattering from a perfect conductor as an integral equation. Let Φ(x, y) = 1 4π e ik|x−y| |x−y| , x = y, denote the fundamental solution of the three-dimensional Helmholtz equation Δu + k 2 u = 0. We define the electric potential and the magnetic potential which are bounded operators from H − 1 2 (Div, ∂D) to H loc (curl 2 , R 3 \∂D). For any radiating solution E ∈ H loc (curl 2 , R 3 \D) of Maxwell's equations, we have the Stratton-Chu repre- The potentials satisfy the following jump conditions on the boundary ∂D, By the mean of the interior and exterior traces of the potentials, we arrive at the electric boundary operator E and magnetic boundary operator H, both bounded linear operators from H − 1 2 (Div, ∂D) to H − 1 2 (Div, ∂D). These operators satisfy Let E be a radiating solution to Maxwell's equations, satisfying a Dirichlet boundary condition γ t E = −F for some right-hand side F, in our case the scattered field E s with F = γ t E i or the domain derivatives E , E with the right hand sides presented in theorems 3.3 and 4.5. In each case, we make the Ansatz E(x) = −Eλ(x), x ∈ R 3 \D for some density λ ∈ H − 1 2 (Div, ∂D). Then, by the trace and γ t E = E we arrive at the indirect electric field equation (EFIE) Assuming k to be no interior eigenvalue of D, the EFIE is uniquely solvable for any right-hand side (see [2]). The major challenge arises from calculating the boundary conditions for the domain derivatives. Recall the boundary condition Numerically calculating the boundary condition requires access to the discrete version of the surface gradient Grad, the rotation operator R, defined by Rγ T ϕ = γ t ϕ, the magnetic trace H τ = (ν × (H × ν)), and the normal component of the electric field E ν . Since basic surface differential operators like surface gradient, surface divergence, and Laplace-Beltrami operator are available in BEMPP, we represent the boundary terms as follows.
We have to calculate discrete products of discretizations for the product h ν E ν and h ν H τ . From (3.5) we conclude E ν = − 1 ik Div ∂D (H × ν), i.e. we can calculate the normal component of E by applying the surface divergence to H × ν. Considering we see, that the negative dual pairing − γ t ϕ, γ T ψ ∂D between H − 1 2 (Div, ∂D) and its dual space H − 1 2 (Curl, ∂D) can be seen as the weak formulation for the rotation operator R.
Since we use the Ansatz E s = −Eλ, the tangential trace of the electric field is given by For the discrete product f · d g = i α i φ i of two functions f and g in a chosen basis of functions φ i , we calculate the L 2 projection of the product onto the bases functions φ i , i.e. we solve the linear system Note that we have to choose a basis of scalar functions for the product h ν E ν and a basis of vector valued functions for the product h ν H τ . For details on the above described implementations and the code of the actual implementations of the first domain derivative and its use to solve an inverse problem, we refer to [8,9] and the tutorials on the homepage of BEMPP (https://bempp.com). For details on how to solve the EFIE with BEMPP, see also [25].
Lets consider now the boundary condition for the second domain derivative E , given by Note, that we have formulated the boundary condition in a way, we can use the same tools as before. We only have to consider additionally a discretization of the curvature operator R and of the mean curvature κ. The discrete scalar product of two functions is another special case of the discrete product · d described above. Recall the definitions R = J ν and κ = 1 2 div ν.
We have ∂ν ∂ν = 0 and R acts only on the tangential plane. Since R = R , we arrive at which is in every component a discrete product of functions. Having calculated each component, we use again L 2 projections to calculate RH τ . For the mean curvature κ, we use the relation −Δ ∂D x i = 2κν i , i ∈ {1, 2, 3}, of κ and the Laplace-Beltrami operator (see equation (2.5.212) in [23] with u = x i ), to calculate The left-hand side can again be implemented by using the discrete product of functions · d and applications of the surface gradient and the surface divergence.

Numerical examples
Before we consider some numerical examples let us summarize on one step of the iteration schemes using the above discretized operators. First for a Newton step: In all presented experiments, we have chosen N = 4, i.e. we used 25 basis functions to describe our reconstructions. For the considered objects, choosing more basis functions did not seem to be worth the greatly increasing computational effort. Now, since we know how to realize the boundary conditions, we present actual reconstructions using the second degree method. We ran reconstructions for exact and also noisy data. As in [3,8,9,15] we consider the following shapes: • A rounded cuboid, implicitly given by Note, that the characterizations of the second domain derivative requires a smooth boundary. We therefore chose the rounded cuboid to challenge our reconstructions with an object close to the non-smooth cuboid. Additionally, we show the cushion as an example for a non-convex object with positive and negative curvature κ. In order to cancel any positive effects due to symmetry, we applied the translation x → x + k 8 (1, 1, 1) such that the center of the rounded cuboid and the cushion-shaped object does not coincide with the center of the star shaped reconstructions.
Furthermore, we consider the plane wave, given by (1 , 2 , 3) ∈ S 2 such that the direction of the plane wave does not coincide with the symmetries of the considered shapes. We also ran experiments where we have considered different star shaped objects, wave numbers, incident directions, or real or complex valued polarizations. All of these led to comparable results as in the presented examples. Additionally, similar to other iterative approaches in acoustic or electromagnetic inverse scattering the Halley method confirmed the known observation that there is some improvement if reconstructions from some far field patterns generated for instance by three or four incident directions are averaged in each iteration step (see [11,15,21]).
In the presented experiments, the wave number is k = 1.4. Note that the diameter of the rounded cuboid and the wavelength 2π/k are of the same magnitude. We calculate the exact data E ∞ = F(∂D) by picking M = 168 evaluation points on the unit sphere S 2 , i.e. E ∞ ∈ C 3×168 . Note, that due to the offset of the exact data, we avoid an inverse crime, since the exact data is being calculated by using meshes unrelated to those used in the reconstructions. We additionally use finer meshes for the calculation of the exact data. For the presented examples of reconstructions, we have used grids with 800 elements and 1200 edges which results in 1200 degrees of freedom. To generate the exact data, we used a grid with 1152 elements and 1728 edges which results in 1728 degrees of freedom. On a machine with 32 CPU cores one iteration took on average about 15 min, while one iteration of the Halley method took on average about 90 min. Note that these computational times can almost surely be improved, since the code is not optimized with respect to a fast performance.
In the case of noisy data of level δ 0, we multiply every element of the far field matrix E ∞ ∈ C 3×168 by a random complex number 1 + δλ 1 e 2πiλ 2 , where λ 1 , λ 2 are uniformly distributed random numbers on (0, 1). In our experiments with noisy data, we have chosen a noise level δ = 0.2 which results in 10% relative error in comparison with the exact data. In each experiment, besides the one presented in figure 5, we have chosen the unit ball as initial guess. To find an appropriate choice for the regularization parameters α 1 , α 2 , we tested a range of parameters α 1 , α 2 and observed in the case of noise free data successful reconstructions for a wide range of parameters. In the case of noisy data, the choice of α 1 = 300 and α 2 = 150 led to reasonable reconstructions in all of our experiments. Considerably lower values led sometimes to deteriorated reconstructions.
To observe stability of the reconstructions with respect to random noise, we ran each experiment with noisy data ten times. In order to judge on the quality of the reconstructions, we considered the relative L 2 (S 2 ) error of the parametrizations of the exact object and the reconstruction. At the above mentioned noise level and with our chosen regularization parameters, one can hardly see which reconstruction is the worst and which one is the best. Considering the parametrizations of the best and the worst reconstructions we observed about 5% difference. For our figures, we picked neither the best nor the worst case to present.
The stopping rule applies, if the residual F(∂D i ) − E ∞ L 2 (S 2 ) / E ∞ L 2 (S 2 ) falls below τ times the noise level, where τ satisfies 2 τ > 1. Since there is no strategy for choosing an optimal τ , we have chosen τ = 1.5 for the final reconstruction in figure 2. With τ = 1.5, one can hardly observe any differences in the reconstructions of the regularized Newton and Halley method. However, by choosing smaller τ , we did observe qualitative differences in the reconstructions. This requires some fine tuning of the parameter τ . To present this behavior, we have chosen τ = 1.13 to pick the final reconstructions in figures 3 and 4.
In most of our experiments we observed a behavior of the residual as shown in figure 1 for the rounded cuboid and the cushion shaped object. It occurs a significant decrease in the first few iteration steps which then slows down rapidly. As seen for the rounded cuboid usually the second degree method shows a faster decrease in the first one to three iteration steps. But finally both methods, the iterative regularized Newton method and the second degree method, lead to similar residual errors in case of noise free data as well as in case of noisy data. Note that there was a case, where the second degree method shows only in the first iteration an improvement in comparison to the Newton type method.
In comparing the reconstructions of both approaches some differences are remarkable. Both the Halley as well as the Newton scheme lead to reasonable, similar reconstructions, as shown in figure 2. The only difference are a few less iteration steps required in the Halley methods, which, of course, are payed for by more computational effort for the second derivative in each step. But, additionally, we observed a more stable performance with respect to the choice of the regularization parameter. The range of possible parameters α 1 within the Newton scheme leading to reasonable results is significantly smaller then for the second degree approach. To be more specific: using α 1 /2 instead of α 1 within the Newton scheme led most of the time to deteriorated shapes after some iterations, whereas the using of α 2 /2 instead of α 2 within the Halley method almost always still produced reasonable reconstructions.
A more significant effect in the reconstructions occurs in case of noisy data. Due to regularization the iteration has to be stopped if the residual error is becoming to small. It is seen also in the acoustic case that the iterated shapes start to deteriorate if a receding cusp is developed. This effect is regularized slightly by choosing the matrix J instead of I in the Tikhonov regularization, but if it occurs, in general the iteration process can not compensate on it and the iteration must be stopped. Of course, by a larger regularization parameter we can avoid the effect but then reconstructions become worse, close to the initial guess. Here we observed the main advantage of considering the second degree method, since it turned out that the method reaches frequently the stopping level before such cusps occur. This can be seen from the figures 3 and 4 showing reconstructions from noisy data with the iterative regularized Newton method and with the Halley method.
Finally, we consider the performance of the schemes in the case, where we know that the tangential cone condition fails. Thus we consider a ball, where the radius is chosen such that the wave number is an eigenvalue of the Laplace-Beltrami operator. From corollary 3.4 we have seen that injectivity of F does not hold if we illuminate the ball by an incident field generated by a vector surface harmonics. Especially F h = 0 if h ν is constant. Thus the iteration schemes cannot just expand or shrink the size of the ball. Exactly this was observed. Using such a ball as an initial guess both the Newton as well as the second degree scheme slow down and do not reach as good reconstructions as in non critical cases (see figure 5).
As a conclusion from these numerical investigations we can state a slightly more stable performance of the Halley method compared to an iterative regularized Newton approach as it was already observed for the acoustic case (see [13]), but by the prize of a higher computational effort in each iteration step. Additionally, the last observations in case of non injective domain derivatives confirm that further research is required in understanding the performance of iterative regularization schemes in inverse obstacle scattering.
Let us consider now With the previous calculations, we get

Recall the definition of the symmetric matrix
Considering vector fields E, V, h and a symmetric matrix A, elementary calculus yields the following identities: With these identities, we conclude We combine (5.6), (5.7) and (5.5) to get By the divergence theorem, we have since h 1 , h 2 are compactly supported in Ω. Application of the partial integration formula leads to In the last integral, we only need to consider the tangential component of (A 1 curlE) × h 2 , which is given by We use again (3.5) to conclude ν curl V = −Div ∂D (ν × V) = 0, since V ∈ H pc (Ω), which yields Finally, we have shown A second application of the divergence theorem, together with the boundary condition ν × E = ν × V = 0 on ∂D leads to Applying (5.8) and (5.9), we have We consider the term curl ((E × h 1 ) × h 2 ). Elementary calculations yields which leads to We apply again the divergence theorem to conclude The second integral occurs to
We finally conclude for all V ∈ H pc (Ω), i.e., (E 1 ) 2 is a radiating solution to homogeneous Maxwell's equations.