Ihsen Alouani
email: ihsen.alouani@gmail.com

Smail Niar
email: smail.niar@univ-valenciennes.fr

Mohamed Jemai

Fadi Kurdahi²

Mohamed Abid

Working-Conditions-Aware Fault Injection Technique 1 Working-Conditions-Aware Fault Injection Technique

Keywords: Dependability, fault tolerant, fault injection, voltage scaling

With new integration rates, the circuits sensitivity to environmental and working conditions has increased dramatically. Thus, presenting reliable, less consuming energy and error resilient architectures is being one of the major problems to deal with. Besides, evaluating robustness and effectiveness of the proposed architectures is also an urgent need. In this paper, we present an extension of SimpleScalar simulation tool having the ability to inject faults in a given cache architecture. We tried to focus on transient errors which are due to cosmic rays (soft errors) or to voltage scaling and high temperatures. During the fault injection process, there is no software traps and the target application isn't modified. Transient errors are modeled by a bit flip. The bit flip may occur in any part of the cache arrays and also in the cache blocks' address. In this paper we present an overview on the existing fault injection techniques. Then, we propose a fuzzy system which models the relationship between working conditions (the supply voltage and temperature) and the error probability. The impact of the injected faults on the cache performance as well as the experimental results are also shown. The results show that the impact of faults depends also on the application itself.

Introduction

The electronic systems dependability is becoming a growing issue especially with new process technologies. Besides, evaluating the dependability of a processing architecture is a complex task. The complexity of this task increases with the electronic circuits' performance requirements in power consumption and rapidity. Modeling analytically the phenomena of occurring faults is very difficult because of the complexity and relative randomness of these phenomena. Fault Injection is defined by [START_REF] Bonso | Fault Injection Techniques and Tools for Embedded Systems Reliability Evaluation[END_REF] as the dependability validation technique that is based on the realization of the controlled experiments where the observation of the system behavior in presence of faults, is explicitly induced by the deliberate introduction (injection) of faults into the system. It consists in the accomplishment of controlled experiments where the observation of the system's behavior in presence of faults is induced explicitly by the introduction (injection) of faults in the system. In the literature, there are four categories of fault injectors: Hardware-based, software-based, simulation-based and emulation-based fault injectors: The hardware-based fault injection technique uses specific additional hardware in order to inject faults in the target system's hardware. Depending on the faults and their locations, hardwareimplemented fault injection methods fall into two categories: [START_REF] Ziade | A Survey on Fault Injection Techniques[END_REF] 1) Hardware fault injection with contact: The fault is injected using physical contact with the target system by introducing voltage or current pulses into the circuit. 2) Hardware fault injection without contact: The fault injector has no physical contact with the target system. Instead, an external source produces some physical phenomenon, such as heavy ion radiation and electromagnetic interference, causing spurious currents inside the target chip.

The software-based fault injection technique involves the modification of the software running on the tested system. Several kinds of faults may be injected: register and memory errors, dropped or replicated network packets and erroneous flags. [START_REF] Ziade | A Survey on Fault Injection Techniques[END_REF] The simulation-based fault injection technique is different from the two first ones as it involves the construction of a simulation model of the analyzed system architecture, including a detailed model of the used processor. Finally, the fourth method is the emulation-based fault injection technique: In order to make use of simulation benefits and take into account the effects due to real working conditions. The target system is implemented onto FPGA. The synthesized circuit is connected to a computer which role is to control the fault injection triggers and display the results. In our case, we opted for the simulation-based fault injection approach. The idea is to inject faults in a given architecture at the simulation level in order to measure its effectiveness and its level of accuracy. This method incurs a low cost and does not require any special-purpose hardware. These features make the fault injection process more flexible. In fact, any signal value can be easily corrupted and the results of the corruption are easily observable regardless of the location of the corrupted signal. This flexibility facilitates error modeling and injection. We extended SimpleScalar3.0 Mips processor simulator [START_REF] Austin | SimpleScalar: an infrastructure for computer system modeling[END_REF] in order to make a simulation environment equipped with a fault injector. Presently, our fault injection tool has been applied on different benchmarks of Mibench [START_REF] Guthaus | A free, commercially representative embedded benchmark suite[END_REF] in order to study the cache system behavior under aggressive working conditions. We note that the majority of works already done on fault injection were focusing on modeling faults, injecting them in hardware or software. Generally, choosing the time and the location in wich the fault injection is triggered is made either randomly or in a deterministic manner. The problem is that this choice may not correspond to the fault occurring phenomenon in a real process. For these reasons, we propose in this paper a fuzzy system that predicts the error probability in terms of supply voltage and temperature. The Fuzzy logic control is an alternative to traditional notions of logic. It has emerged as a profitable tool for the controlling and steering of systems and complex industrial processes, as well as for other expert systems and applications. As the nature of the studied phenomena is more fuzzy than deterministic, the fuzzy logic tool corresponds to the features of our problem. This paper is organized as follows. The proposed fault injector and the proposed predictor of error probability are discussed in section 2. Section 3 shows the impact of injected faults on the cache system performance. Finally, section 4 concludes the paper and opens horizons to future work.

Fault Injection Process

A. The Proposed Fault Injector

In this paper, we are focusing on transient errors that are triggered by environmental conditions such as power-line fluctuation, electromagnetic interference, or cosmic radiations (soft errors). Generally, transient errors are modeled by a bit flipping even if it can cause a permanent stuck at 0 or stuck at 1 [START_REF] Ziade | A Survey on Fault Injection Techniques[END_REF]. In terms of triggering time, the transient errors are generally considered as randomly occurring phenomena. The majority of related works in the literature dealing with this type of faults uses either deterministic or random fault injection. The deterministic method consists in injecting faults in preselected times before running. However, the random technique randomly carries out the injection during runtime. In our case, we are modeling transient errors using a bit flip in the cache arrays. The error may occur either in data bloc, tag or address bits. If a bit flip occurs in a data block, it leads automatically to an execution failure if the affected bloc is accessed after fault injection. However, in the case of either tag or address, it can produce an execution failure or simply produce a miss and increase the latency of the cache access. In this work, we propose a pseudo-random process of fault injector triggering. Our fault injection mechanism is based on fault rates of transient errors. These rates are deducted from the percentage of saved supply voltage and the working temperature by a proposed fuzzy system of error probability prediction. In addtion to the triggering time, the second issue our paper proposes to study is the location where errors has to be injected within the cache. In the literature, the fault injection location as well as the fault injection time, is either considered as random location, or deterministic choice, i.e. pre-defined locations [START_REF] Ziade | A Survey on Fault Injection Techniques[END_REF]. In our case, we chose to inject faults in the same address of the cache access. This can be considered as "the worst case" methodology. In fact, choosing a given address to inject fault doesn't lead necessarily and immediately to an error because the address may not be used immediately. However, if we inject a fault in the selected line address to access, it is more likely that this bloc would be needed in the execution.

B. Error Probability Prediction System

The majority of works already done on fault injection were focusing on modeling faults, injecting them in hardware or software. They generally choose either random or deterministic methods to define time and location of the fault injection process. The problem is that this choice may not represent correctly how faults occur in reality. In fact, supposing that faults are occurring randomly may imprecisely model the phenomena behind faults taking place. Faults occur randomly in time, but they are not random phenomena in terms of working conditions. Thus, the whole fault injection tool may be far from reality and in a context of evaluating fault resilient architecture, it can lead to refusing an efficient architecture or accepting a weak one. Hence, the need to build a model that mimes the real phenomenon of occurring faults is being imperative. The figure below represents the fault injection mechanism that we propose to design. The inputs are temperature and saved percentage of supply voltage which is the reduction in Vdd per cent. In fact, in a context of saving energy and reducing power consumption, reducing the supply voltage leads to increasing the cache failure probability. Besides, as shown in [START_REF] Khajeh | TRAM: A Tool for Temperature and Reliability[END_REF], the temperature of the circuit has an important impact on error probability. As there is a great difficulty in modeling analytically the phenomenon of fault occurring, we have though to use a fuzzy system that is able to integrate the interference between working conditions and describe the behavior of fault occurring mechanism in the reality. Fuzzy Logic was initiated in 1965 [START_REF] Zadeh | Fuzzy Sets, Information and Control[END_REF] by Lotfi Zadeh. Basically, Fuzzy Logic (FL) is a multivalued logic that allows intermediate values to be defined between conventional evaluations like true/false, yes/no, high/low…etc. Notions like rather tall or very fast can be formulated mathematically and processed by computers, in order to apply a more human-like way of thinking in the programming of computers [START_REF] Guthaus | A free, commercially representative embedded benchmark suite[END_REF]. Our proposed fuzzy system's role is to estimate the error probability value in terms of the chip working conditions and communicate it with the simulator.

FIG. 1 -Diagram of fault injection process.

Given a set of rules describing the behavior of the error probability taken from [START_REF] Swift | Single-event upset susceptibility testing of the Xilinx Virtex II FPGA[END_REF], we have to build a fault injection approach that models the way error probability evolves while working conditions like temperature and saved power are changing. So, we have to give a model that takes temperature and saved power as inputs and gives the probability of error as an output. In the proposed system, the modeling of the relationship between the probability of errors and the working conditions (Temperature and Voltage) is based on observed behavior of error probability under changed running parameters. Based on these observations we proposed a fuzzy logic system that transforms inputs from numbers to fuzzy sets admitting a number of membership functions. After that, basing on fuzzy rules, the system extracts the output using a given defuzzification method (centroïd in our case). Then the system construction is presented as follows (figure 2). FIG. 2 -General form of the proposed fuzzy system.

In our case, we have two inputs ("temperature" and "saved supply voltage percentage") and one output ("error probability"). The range of temperature we have chosen is between 20°C and 120°C distributed on three fuzzy sets: mild, medium and aggressive. However, for the saved power percentage, we divide the range to 5 fuzzy sets: low, mild, medium, aggressive and destructive. The choice we made about the number of fuzzy sets can be explained by referring to [START_REF] Khajeh | TRAM: A Tool for Temperature and Reliability[END_REF] which shows that, below the nominal Vdd, the impact of the voltage scaling on cache failure probability is more important than the temperature one. Thus, while in our case we do not surpass the nominal supply voltage, the voltage scaling impact has to be more precisely described in the fuzzy system design. The output variable which is the error probability for the SRAM is placed in a range between 10 -6 and 10 -1 and this choice is chosen in harmony with the approximation for SRAM error probabilities maid in [START_REF] Sampson | EnerJ: Approximate Data Types for Safe and General Low-Power Computation[END_REF]. It is consisting of 5 fuzzy sets which are: low, mild, medium, high, aggressive, very aggressive and destructive. The rules we have implemented are shown in the table below. We implemented the fuzzy inference system using Fuzzy Logic Toolbox of Matlab R2007a and we got the following variations of error probability in terms of temperature and saved Vdd percentage. Figure 3 illustrates the exponential growth of the cache failure probability while reducing supply voltage. We suppose in our case that the frequency is maintained constant. However, in reality, changing cycle time leads to different probability of failure [START_REF] Khajeh | TRAM: A Tool for Temperature and Reliability[END_REF]. This means that our predictor targets a precised architecture with given performance and has to be changed with other cases. FIG. 3 -Error probability in cache in terms of saved supply voltage percentage and temperature. The range of error probabilities are obtained for 32-nm technology for a 32kB cache [START_REF] Sampson | EnerJ: Approximate Data Types for Safe and General Low-Power Computation[END_REF].

Increasing the percentage of the gain means a reduction in the supply voltage. Therefore, as shown in the figure 3, varying the saved Vdd percentage from 0% to 50% leads to multiplying the error probability by 10 4 . We also note that increasing the running temperature increases the error probability.

The Impact of Injected Faults on The Cache Performance

In order to measure the impact of error rate on cache performance, we varied error probabilities in a range from 10 -6 to 10 -3 which is in accordance with [START_REF] Sampson | EnerJ: Approximate Data Types for Safe and General Low-Power Computation[END_REF]. We represented its impact on some cache performance as shown in figure 4 below.

FIG. 4-The impact of the injected errors on the first level instruction cache (IL1). The probability varies

from 0 to 0.01. The IL1 size is 16 KB. Basic-math benchmark (158M Instructions) was used.

From this figure we can notice the increase in instruction cache misses from 800 to more than 1200 misses while injecting faults. The increasing number of misses means an overhead in time and power which represents a degradation in cache performance. Figure 4 gives also the progress of the injected faults number while increasing the error probabilities.

After studying the impact of injecting errors on cache behavior in terms of data and misses, we modified the SimpleScalar based simulator sim-wattch in order to get information about the impact of injected errors on the cache power and energy consumption. probability in the Rijndael benchmark from Mibench. We note that injecting faults degrades the processor performance and increases the power consumption. This is due to the raise of the occurring misses' number. In fact, its direct impact is to decrease the executed instructions per cycle. This is due to the lost time in additional misses. In parallel with the IPC decrease, we note that the power consumption of the cache is increasing while rising fault injection probability.

FIG. 5-Variation of IPC and the instruction cache power consumption for the Rijndael benchmark (39M instructions). The IL1 and DL1 size is set to 16KB.

The same behavior is recorded with Basic Math benchmark. In the figure 6, we represent the total power consumption per cycle as well as instruction per cycle evolution in terms of injected fault probability.

FIG. 6-Variation of the IPC and the instruction cache power consumption for the Rijndael benchmark (39M instructions). The IL1 and DL1 size is set to 16KB.

From the figure above we validate the remarks done for Rijndael application and we can note the increase of the total power consumption with the rising number of injected faults. In fact, the additional number of cache misses creates a need to additional power consumed to access next level caches.

After studying the impact of injected faults on power consumption and processor performance, we tried to show the impact of the cache size on the power consumption while the fault injection probability is maintained. The following figure shows the evolution of cache power consumption and IPC in terms of the cache size. We fixed the error probability and varied the first level cache size from 8 kB to 128kB. We note that the IPC is increasing with cache size and this is explained by the decrease of the misses number. At the same time, spite of the diminution of the cache misses, the power consumption increased. In fact, raising the cache size leads to increasing the cache power consumption. So, a designer has to look for the appropriate tradeoffs between reducing power consumption which result in raising the error probability and increasing the cache size to deal with performance degradation which result in increasing power consumption.

1,045E+09

1,045E+09

Conclusion

In this paper, we persent a fault injection environment at the simulation level. It is implemented in C language as an extension of SimpleScalar simulator. It simulates transient faults occurring in a cache architecture in runtime. The work focuses on transient faults because of their increasing impact especially with new transistors integration rates and process technologies. The injected errors are injected during cache access and modeled by a bit flip in the cache arrays. The fault injection process is random in time. However, it is controllable via an error probability that is estimated using a proposed fuzzy logic predictor. The latest, which is the main contribution of this paper, is an inference system that predicts the error probabilitiy using both estimated circuit temperature and saved supply voltage percentage. The current version can be improved in several ways. We are currently working on a possible extension aiming to integrate fault tolerant architectures and implementing a fault injection process in the pipeline stages and the processor's registers. We also intend to generalize the approach by integrating the running frequency as an input to the error probability predictor.

Figure 4

 4 represents the variation of instruction cache power consumption, and instruction per cycle (IPC) while increasing the error

  Variation of data cache power, total power per cycle and IPC in terms of cache size (from 8kB to 128kB) 2-way associative cache. The error probability is fixed and it's equal to 10

		1,92				1,046E+09
		1,915				1,045E+09 1,046E+09
		1,91				1,045E+09	ipc
	FIG. 7-	1,905				icache power
		0	5E-06	1E-05	1,9E-05	8,5E-05
		0,902				1,379E+10
		0,901 0,9012 0,9014 0,9016 0,9018				1,378E+10 1,378E+10 1,378E+10 1,379E+10 1,379E+10 1,379E+10
		0,9008				1,378E+10
						IPC
						total power per
						cycle

-4


Frédérique Le Louër
email: f.lelouer@math.uni-goettingen.de

Fast spectral methods for the shape identification problem of a perfectly conducting obstacle

Keywords: Maxwell equations, perfectly conducting interface, boundary integral equation, spectral method, Fréchet derivative, regularized Newton method

We are concerned with fast methods for the numerical implementation of the direct and inverse scattering problems for a perfectly conducting obstacle. The scattering problem is usually reduced to a single uniquely solvable modified combined-field integral equation (M-CFIE). For the numerical solution of the M-CFIE we propose a new high-order spectral algorithm by transporting this equation on the unit sphere via the Piola transform. The inverse problem is formulated as a nonlinear least squares problem for which the iteratively regularized Gauss-Newton method is applied to recover an approximate solution. Numerical experiments are presented in the special case of star-shaped obstacles.

Introduction

In this paper we are concerned with the numerical aspects of the solution of the direct and inverse electromagnetic scattering problems for a three-dimensional bounded and perfectly conducting obstacle lit by time-harmonic incident plane waves.

We assume the perfect conductor be represented by a bounded domain Ω in R 3 . Let Ω c denote the exterior domain R 3 \Ω and n denote the outer unit normal vector to the boundary Γ. Let κ denote the exterior wavenumber. The propagation of electromagnetic waves are governed by the system of Maxwell equations and the time-harmonic Maxwell system can be reduced to a second order equation for the electric field only. In this case the forward problem is formulated as follows : Given an incident electric wave E inc which is assumed to solve the second order Maxwell equation in the absence of any scatterer, find the electric scattered wave E s solution to the time-harmonic Maxwell equation

curl curl E s -κ 2 E s = 0 in Ω c , (1.1a)
and satisfying the boundary condition,

n × (E s + E inc ) = 0 on Γ. (1.1b)
In addition the scattered field E s has to satisfy the Silver-Müller radiation condition:

lim |x|→+∞ |x| curl E s (x) × x |x| -iκE s (x) = 0. (1.1c)
This scattering problem can be reduced in several different ways to a single modified combinedfield boundary integral equation (M-CFIE). In section 2, we describe, via direct and indirect approaches, the single integral equations that are usually used to solve the problem (1.1a)- (1.1c).

The construction of efficient inverse obstacle scattering algorithm greatly depends on an accurate approximation of the solution to the forward problem. Among all the existing numerical methods to solve integral equations, we focus here on spectral methods. Ganesh and Hawkins already proposed in [START_REF] Ganesh | A spectrally accurate algorithm for electromagnetic scattering in three dimensions[END_REF][START_REF]An efficient surface integral equation method for the time-harmonic Maxwell equations[END_REF][START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF][START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF] several spectrally accurate methods to implement the electric field integral equation (EFIE). One of them [START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF] consists in transporting the EFIE on the unit sphere via a normal transformation acting from the tangent plane to the boundary Γ onto the tangent plane to the unit sphere, so that one only have to seek a solution in terms of tangential vector spherical harmonics. To treat numerically the M-CFIE, we need to implement the hypersingular part of the electromagnetic double layer boundary integral operator. This requires a transformation that move the nullspace of the surface divergence on Γ onto the nullspace of the surface divergence on S 2 . The implementation of a hypersingular integral can then be avoided by involving integration by parts and surface derivatives of the vector spherical harmonics that can be computed analytically. Such a transformation is well known as the Piola transform. In section 3 we give the reformulation of the integral equation in spherical coordinates, based on this approach, and describe an appropriate splitting of the various singularities. The fully discrete high order spectral algorithm and numerical examples are presented in section 4. The Piola transform is already used by Hohage and Le Louër [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF] to implement systems of second-kind boundary integral equations in an inverse dielectric obstacle scattering algorithm. Therefore, section 3 and 4 use notations and contain various results from [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF] that are essential for this paper.

The radiation condition implies that the scattered field E s has an asymptotic behavior of the form

E s (x) = e iκ|x| |x| E ∞ ( x) + O 1 |x| , |x| → ∞,
uniformly in all directions x = x |x| . The far-field pattern E ∞ is a tangential vector function defined on the unit sphere S 2 of R 3 and is always analytic.

Let consider the scattering of m incident plane waves of the form E inc k (x) = p k e iκx•d k where d k , p k ∈ S 2 and d k • p k = 0. We denote by F k the boundary to far-field operator that maps the boundary Γ onto the far-field pattern E ∞ k of the solution to the forward problem (1.1a)-(1.1c) for the incident wave E inc k . The inverse problem is formulated as follows: Given noisy far field measurements E ∞ 1,δ , . . . , E ∞ m,δ obtained from the scattering of the m incident plane waves characterized by the couples of directions and polarizations (d k , p k ) k=1,...,m , solve

F k (Γ) = E ∞ k,δ , for k = 1, . . . , m. (1.2)
Here, the noise level is measured in the L 2 -norm, i.e.

m k=1 ||E ∞ k,δ -E ∞ k || 2 L 2 1 2
< δ, and the error bound δ is assumed to be known. Kress [START_REF] Kress | Electromagnetic waves scattering : Scattering by obstacles[END_REF] proved that a perfect conductor can be uniquely determined from the knowledge of the far-field pattern for all incoming plane waves. However, it remains an open question whether or not the boundary Γ is uniquely determined by only a finite number of incoming plane waves.

To solve the inverse problem (1.2), we reformulate it as a nonlinear equation posed on an open set of parametrized boundaries. Then we apply the iteratively regularized Gauss-Newton method (IRGNM) to the nonlinear equation via first order linearization. The computation of the iterates requires the analysis and an explicit form of the first Fréchet derivative of the parametrized form of the boundary to far-field operator. The first Fréchet derivative is usually characterized as the far-field pattern of the solution to a new exterior boundary value problem. The fast spectral solver of the direct problem is then incorporated in the IRGNM to compute the far-field pattern of the solution and its derivative at each iteration steps. In electromagnetism, Fréchet differentiability was first investigated by Potthast [START_REF] Potthast | Domain derivatives in electromagnetic scattering[END_REF] for the perfect conductor problem, using the boundary integral equation approach. The characterization of the first derivative was then improved by Kress [START_REF] Kress | Electromagnetic waves scattering : Scattering by obstacles[END_REF]. In section 5, we formulate the IRGNM for the inverse problem (1.2). We also recall the main result on the Fréchet differentiability of the boundary to far-field operator and give a characterization of the adjoint operator, following ideas of [START_REF]Iterative Methods in Inverse Obstacle Scattering: Regularization Theory of Linear and Nonlinear Exponentially Ill-Posed Problems[END_REF][START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF], which is needed in the implementation of the utilized regularized Newton method. We present numerical experiments in the special case of star-shaped obstacles.

The solution of the perfect conductor problem

In this paper, we will assume that the boundary Γ of Ω is a smooth and simply connected closed surface, so that Ω is diffeomorphic to a ball.

Notation: We denote by H s (Ω), H s loc (Ω c ) and H s (Γ) the standard (local in the case of the exterior domain) complex valued, Hilbertian Sobolev space of order s ∈ R defined on Ω, Ω c and Γ respectively (with the convention H 0 = L 2 .) Spaces of vector functions will be denoted by boldface letters, thus H s = (H s ) 3 . If P is a differential operator, we write:

H s (P, Ω) = {v ∈ H s (Ω) : Pv ∈ H s (Ω)}, H s loc (P, Ω c ) = {v ∈ H s loc (Ω c ) : Pv ∈ H s loc (Ω c )}.
The space H s (P, Ω) is endowed with the graph norm ||v||

H s (P,Ω) = ||v|| 2 H s (Ω) + ||Pv|| 2 H s (Ω) 1 2 .
This defines in particular the Hilbert spaces H s (curl, Ω) and H s (curl curl, Ω) and the Fréchet spaces H s loc (curl, Ω) and H s loc (curl curl, Ω). When s = 0 we omit the upper index 0. We use the following surface differential operators: The tangential gradient denoted by grad Γ , the surface divergence denoted by div Γ , the tangential vector curl denoted by curl Γ and the surface scalar curl denoted by curl Γ . For their definitions we refer to [15, Appendix A] or [17, pages 68-75]. For s ∈ R, we introduce the Hilbert space

H s div (Γ) = {j ∈ H s (Γ); j • n = 0 and div Γ j ∈ H s (Γ)} , H s curl (Γ) = {j ∈ H s (Γ); j • n = 0 and curl Γ j ∈ H s (Γ)} ,
endowed with the norms

|| • || H s div (Γ) = || • || 2 H s (Γ) + || div Γ • || 2 H s (Γ) 1/2 , || • || H s curl (Γ) = || • || 2 H s (Γ) + || curl Γ •|| 2 H s (Γ) 1/2 .
Recall that for a vector function u ∈ H(curl, Ω) ∩ H(curl curl, Ω), the traces n × u |Γ and

n × curl u |Γ are in H -1 2 div (Γ). The dual space of H -1 2 div (Γ) for the L 2 duality product is H -1 2
curl (Γ) and the exterior product with the normal vector defines a bicontinuous isomorphism between

H -1 2 div (Γ) and H -1 2 curl (Γ); so that (n × u |Γ ) × n and (n × curl u |Γ ) × n are in H -1 2 curl (Γ). Let Φ(κ, z) = e iκ|z| 4π|z|
be the fundamental solution of the Helmholtz equation ∆u + κ 2 u = 0.

For any solution E s to the Maxwell equation (1.1a) that satisfies the radiation condition (1.1c), it holds the Stratton-Chu representation formula:

E s (x) = Γ curl x Φ(κ, x -y) n(y) × E s (y) ds(y) + 1 κ 2 Γ curl curl x Φ(κ, x -y) n(y) × curl E s (y) ds(y).
(2.1)

Conversely, for any j ∈ H Using these results, the forward problem (1.1a)-(1.1c) can be reduced, in several different ways, to a single uniquely solvable modified combined-field integral equation (M-CFIE). We will consider the following two different approaches. The first one can be used to solve electromagnetic scattering problem with general Dirichlet conditions of the form n × E s = f where f ∈ H -1 2 div (Γ) is given. It is based on the layer ansatz :

E s = Dj + iηS Λj, (2.2)
where Λ is a bounded operator from H -1 2 div (Γ) to itself, self-adjoint and elliptic for the bilinear form

(j, m) → Γ j • (n × m) ds (2.3)
and η is a non vanishing real constant. By the jump relations, the field E s given by (2.2) solves the Dirichlet boundary value problem (1.1a)-(1.1c) if the density j solves the following integral equation

(I -M κ -iηC κ Λ)j = 2f on Γ. (2.4)
Here the single layer potential C κ and the double layer potential M κ are defined by

M κ j(x) = - Γ n(x) × curl x {2Φ(κ, x -y)j(y)}ds(y), C κ j(x) = - 1 κ Γ n(x) × curl curl x {2Φ(κ, x -y)j(y)}ds(y) = -κ n(x) × Γ 2Φ(κ, x -y)j(y)ds(y) + 1 κ curl Γ Γ 2Φ(κ, x -y) div Γ j(y)ds(y).
The operator

M κ : H -1 2 div (Γ) → H -1 2 div (Γ) is compact and the operator C κ has a hypersingular kernel but is bounded on H -1 2 div (Γ). The operator Λ is then chosen such that (I -M κ -C κ Λ) is a Fredholm operator of index zero. Kress first proposed a compact regularization Λj = n × S 2
0 j where S 0 is the single layer boundary integral operator of the Laplace equation [START_REF] Colton | Integral Equation Methods in Scattering Theory[END_REF], thus (I + M κ + C κ Λ) is a Fredhom operator of the second kind. One can also use the elliptic and invertible operator Λj = C 0 j which is a variant of the operator C κ [START_REF] Costabel | On the Kleinman-Martin integral equation method for the electromagnetic scattering problem by a dielectric body[END_REF][START_REF] Steinbach | Modified combined field integral equations for electromagnetic scattering[END_REF] defined on H

-1 2 div (Γ) by C 0 j = n × S 0 j + curl Γ S 0 div Γ j.
The far-field pattern can be computed via the integral representation formula E ∞ = F ∞ j where the far-field operator

F ∞ : H -1 2 div (Γ) → L 2 t (S 2
) is defined for x ∈ S 2 by :

F ∞ j( x) = iκ 4π Γ e -iκ x•y x × j(y) ds(y) + iη κ 4π Γ e -iκ x•y x × Λj(y) × x ds(y).
The alternative boundary integral equation approach is based on the Stratton-Chu representation formula (2.1) and the following equality for

x ∈ Ω c 0 = 1 κ 2 Γ curl curl x Φ(κ, x -y) n(y) × curl E inc (y) ds(y) + Γ curl x Φ(κ, x -y) n(y) × E inc (y) ds(y).
Adding the last equation to (2.1) and using the boundary condition of E s , it yields

E s (x) = 1 κ 2 Γ curl curl x Φ(κ, x -y) n(y) × curl E s (y) + E inc (y) ds(y), x ∈ Ω c .
Taking the Neumann and Dirichlet trace of the above equation we obtain

C κ 1 κ n × curl E s + E inc = 2n × E inc , and
(I + M κ ) 1 κ n × curl E s + E inc = 2 1 κ n × curl E inc . Setting m = κ -1 n × curl E s + E inc ,
by a linear combination of the two above equalities we obtain the integral equation

(I + M κ + iηΛC κ )m = 2 1 κ n × curl E inc + iηΛ n × E inc on Γ, (2.5)
where Λ is still an elliptic operator for the bilinear form (2.3).

The far-field pattern can then be computed via the integral representation formula

E ∞ ( x) = κ 4π Γ e -iκ x•y x × m(y) × x ds(y).

Spherical reformulation of the integral equation

We transport the surface integral equations (2.4) and (2.5) on a spherical reference domain and work in the spherical coordinate system. We denote by θ, φ the spherical coordinates of any point x ∈ S 2 , that means we can write

x = ψ(θ, φ) = (sin θ cos φ, sin θ sin φ, cos θ) T , with (θ, φ) ∈ ]0; π[×[0; 2π[ ∪{(0, 0); (0, π)}.
The tangent and the cotangent planes at any point x = ψ(θ, φ) ∈ S 2 is generated by the unit vectors e θ = ∂ψ ∂θ (θ, φ) and e φ = 1 sin θ ∂ψ ∂φ (θ, φ). The triplet ( x, e θ , e φ ) forms a direct orthonormal system. The determinant of the Jacobian of the change of variable is J ψ (θ, φ) = sin θ.

Let q : S 2 → Γ be a parametrization of class C 1 at least. The total derivative [Dq( x)] maps the tangent plane T x to S 2 at the point x onto the tangent plane T q( x) to Γ at the point q( x). The latter is generated by the vectors t 1 ( x) = [Dq( x)]e θ and t 2 ( x) = [Dq( x)]e φ . The determinant J q of the Jacobian of the change of variable q : S 2 →∈ Γ and the normal vector n • q can be computed via the formulas J q = t 1 × t 2 and τ q (n) = t1×t2

Jq . The parametrization q : S 2 → Γ being a diffeomorphism, we set [Dq( x)] -1 = [Dq -1 ] • q( x). The transposed matrix [Dq( x) * ] -1 maps the cotangent plane T *

x to S 2 at the point x onto the cotangent plane T * q( x) to Γ at the point q( x). The latter is generated by the vectors t 1 ( x) = t2(q( x))×n(q( x))

Jq( x)
= [Dq( x) * ] -1 e θ and t 2 ( x) = n(q( x))×t1(q( x))

Jq( x)
= [Dq( x) * ] -1 e φ . To transport the boundary integral equations on the unit sphere, we use the Piola transform of q, which is an invertible operator from H

-1 2 div (Γ) to H -1 2 div (S 2
) defined as follows [15, Lemma 3.1]:

P q : H -1 2 div (Γ) -→ H -1 2 div (S 2 ) j → j s = J q [Dq] -1 (j • q). (3.1)
This allows us to use the following identities stated in [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF]:

(grad Γ u) • q = [Dq * ] -1 grad S 2 (u • q), (curl Γ u) • q = 1 J q [Dq] curl S 2 (u • q), (div Γ v) • q = 1 J q div S 2 J q [Dq] -1 (v • q) , (curl Γ w) • q = 1 J q curl S 2 [Dq * ](w • q) . (3.2)
The construction of the spectral method is based on the combination of the boundary integral operators M κ and C κ with the Piola transform (3.1) and its inverse. In practice, we apply the operator (3.1) to both sides in the boundary integral equations (2.4) and (2.5) and we insert the identity

I H -1 2 div (Γ)
= P q P -1 q between the integral operators and the unknown. We thus replace the operators Λ, M κ and C κ in (2.4) and (2.5) by the following operators defined for any density

j s ∈ H -1 2 div (S 2
) by:

Λ s j s = P -1 q ΛP q j s , M κ j s = P -1 q M κ P q j s = -J q [Dq] -1 S 2 τ q (n) × curl{2Φ(κ, |q(•) -q( y)|)[Dq( y)]j s ( y)}ds( y), C κ j s = P q C κ P q j s = -κ J q [Dq] -1 S 2 τ q (n) × {2Φ(κ, |q(•) -q( y)|)[Dq( y)]j s ( y)}ds( y) + 1 κ curl S 2 S 2
{2Φ(κ, |q(•) -q( y)|) div S 2 j s ( y)}ds( y).

The new unknown will be a tangential vector density in H

H -1 2 div (S 2 ) × H -1 2 div (S 2 ) → C (j s , m s ) → S 2 j s • ( x × m s ) ds (3.3) Proof. Let j s ∈ H -1 2 div (S 2
). Notice that we have (P -1 q j s ) • q = 1 Jq [Dq]j s . In view of the definition of the vectors t 1 and t 2 , we obtain the following equality :

n × P -1 q j s • q = [Dq * ] -1 ( x × j s ) . (3.4)
Now, we set j = P -1 q j s . We have Λj = P -1 q Λ s j s and

Γ Λj • (n × j) ds = S 2 J q ((Λj) • q) • ((n × j) • q) ds = S 2 J q P -1 q Λ s j s • n × P -1 q j s • q ds = S 2 [Dq]Λ s j s • [Dq * ] -1 ( x × j s ) ds = S 2 Λ s j s • ( x × j s ) ds.
It follows that the ellipticity of Λ s for the bilinear form (3.3) results from the ellipticity of Λ for the bilinear form (2.3), and conversely.

Remark 3.2 A suitable choice is to define Λ s such that j s → (Λ s j s ) × x is the duality operator from H -1 2 div (S 2 ) to H -1 2 curl (S 2 ) [17, pp. 208], that is Λ s j s = -∇ S 2 (-∆ S 2 ) -3 2 curl S 2 -curl S 2 (-∆ S 2 ) -1 2 div S 2 .
In practice we will use the operator defined by Λ s j s = x × j s .

The implementation of the operator M κ is already discussed in [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF]. Here, we focus on the operator C κ . We introduce the functions :

S 1 (q; κ, x, y) = 1 2π cos(κ|q( x) -q( y)|), S 2 (q; κ, x, y) = 1 2π    sin(κ|q( x) -q( y)|) |q( x) -q( y)| x = y, κ x = y. and R(q; x, y) = | x -y| |q( x) -q( y)| .
The operator C κ can then be rewritten as

C κ j s ( x) = S 2 R(q; x, y) | x -y| W 1 (q; κ, x, y)j s ( y)ds( y) + i S 2
W 2 (q; κ, x, y)j s ( y)ds( y)

+ 1 κ curl S 2 S 2 R(q; x, y) | x -y| S 1 (q; κ, x, y) div S 2 j s ( y)ds( y) + i κ curl S 2 S 2 S 2 (q; κ, x, y) div S 2 j s ( y)ds( y)
where W 1 (q; κ, x, y) and W 2 (q; κ, x, y) are 3 × 3 matrices given by W 1 (q; κ, x, y) = κS 1 (q; κ, x, y)V(q; x, y), W 2 (q; κ, x, y) = κS 2 (q; κ, x, y)V(q; x, y),

with V(q; x, y) = t 2 ( x) • t 1 ( y) e θ ( x) ⊗ e θ ( y) + t 2 ( x) • t 2 ( y) e θ ( x) ⊗ e φ ( y) -t 1 ( x) • t 1 ( y) e φ ( x) ⊗ e θ ( y) -t 1 ( x) • t 2 ( y) e φ ( x) ⊗ e φ ( y).
Next, we introduce a change of coordinate system in order to move all the singularities in the weakly singular integrals to only one point that is chosen to be the North pole. For x ∈ S 2 we consider an orthogonal transformation T x which maps x onto the North pole denoted by η. If x = x(θ, φ) then T x := P (φ)Q(-θ)P (-φ) where P and Q are defined by

P (φ) =   cos φ -sin φ 0 sin φ cos φ 0 0 0 1   and Q(θ) =   cos θ 0 sin φ 0 1 0 -sin φ 0 cos φ   .
We also introduce an induced linear tranformation T x defined by T x u( y) = u(T -1

x y) and we still denote by T x its bivariate analogue

T x v( y 1 , y 2 ) = v(T -1
x y 1 , T -1 x y 2 ). If we write z = T x y then we have the identity

| x -y| = |T -1 x ( η -z)| = | η -z|.
The boundary integral operator C κ can be rewritten in the form:

C κ j s ( x) = S 2 T x R(q; η, z) | η -z| T x W 1 (q; κ, η, z) + iT x W 2 (q; κ, η, z) T x j s ( z)ds( z) + 1 κ curl S 2 S 2 T x R(q; η, z) | η -z| T x S 1 (q; κ, η, z) + iT x S 2 (q; κ, η, z) T x div S 2 j s ( z)ds( z), (3.5)
and it can be shown that the mappings (θ , φ ) → T x R(q; η, z(θ , φ ))T x S 1 (q; κ, η, z(θ , φ )) and (θ , φ ) → T x R(q; η, z(θ , φ ))T x C 1 (q; κ, η, z(θ , φ )) are smooth. An important point is that the singularity

1 |η-z(θ ,φ )| = 1 2 sin θ 2
is cancelled out by the surface element ds( z) = sin θ dθ dφ .

Finally, setting f s = P q f , the parametrized form of the equation (2.4) is

I H -1 2 div (S 2 ) -M κ -iηC κ Λ s j s = 2f s , on S 2 , (3.6)
and setting u inc 1 = P q n × E inc and u inc 2 = P q n × curl E inc the parametrized form of the equation (2.5) is

I H -1 2 div (S 2 ) + M κ + iηΛ s C κ j s = 2 1 κ u inc 2 + iηΛ s u inc 1 , on S 2 . (3.7)

A high order spectral algorithm

The parametrized boundary integral equations (3.6) and (3.7) are now propicious to obtain spectrally accurate approximations of the solution by combining the spectral method of Ganesh and Graham [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF] and the hybrid spectral method of Ganesh and Hawkins [START_REF]A hybrid high-order algorithm for radar cross section computations[END_REF][START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF]. The numerical scheme is based on the numerical integration formula over the unit sphere for a continuous function

S 2 u( x)ds( x) ≈ 2n+1 ρ=0 n+1 τ =1 µ ρ ν τ u( x(θ τ , φ ρ )), (4.1)
where θ τ = arccos ζ τ where ζ τ , for τ = 1, . . . , n + 1, are the zeros of the Legendre polynomial P 0 n+1 of degree n + 1 and ν τ , for τ = 1, . . . , n + 1, are the corresponding Gauss-Legendre weights and

µ ρ = π n + 1 , φ ρ = ρπ n + 1
, for ρ = 0, . . . , 2n + 1.

For l ∈ N and 0 ≤ j ≤ l, let P j l denote the j-th associated Legendre function of order l. When j = 0 we write P 0 l = P l . The spherical harmonics defined by Y l,j ( x) = (-1)

|j|+j 2 2l + 1 4π (l -|j|!) (l + |j|!) P |j| l (cos θ)e ijφ
for j = -l, . . . , l and l = 0, 1, 2, . . . form a complete orthonormal system in L 2 (S 2 ). The formula (4.1) is exact for the spherical polynomials of order less than or equal to 2n + 1 (see [START_REF]Vector hyperinterpolation on the sphere[END_REF]). This induces the discrete inner product (•,

•) n (ϕ 1 , ϕ 2 ) n = 2n+1 ρ=0 n+1 τ =1 µ ρ ν τ ϕ 1 ( x ρτ )ϕ 2 ( x ρτ ),
on the space P n of all scalar spherical polynomials of degree less than or equal to n. We have (ϕ 1 , ϕ 2 ) S 2 = (ϕ 1 , ϕ 2 ) n for all ϕ 1 , ϕ 2 ∈ P n . We introduce a projection operator n on P n defined by

n u = n l=0 l j=-l (u, Y l,j ) n Y l,j ,
and we set

L n u =   n u 1 n u 2 n u 3   where u = (u 1 , u 2 , u 3 ) T .
The tangential vector spherical harmonics defined by

Y (1) l,j = 1 l(l + 1) grad S 2 Y l,j and Y (2) l,j = 1 l(l + 1) curl S 2 Y l,j
for j = -l, . . . , l and l = 1, 2, . . . form a complete orthonormal system in L 2 t (S 2 ). The formula (4.1) is exact for the vector spherical polynomials of order less than or equal to 2n and we have Y

(k) l,j • Y (k )
l ,j ∈ P 2n+1 for l, l ≤ n, |j| ≤ l, |j | ≤ l and k, k = 1, 2 (see [START_REF]Vector hyperinterpolation on the sphere[END_REF]). This induces the discrete inner product

(v 1 |v 2 ) n = 2n+1 ρ=0 n+1 τ =1 µ ρ ν τ v 1 ( x ρτ ) • v 2 ( x ρτ ), on the subspace T n ⊂ H -1 2 div (S 2
) of finite dimension 2(n + 1) 2 -2 generated by the orthonormal basis of tangential vector spherical harmonics of degree less than or equal to n ∈ N. We introduce a projection operator L n on T n defined by

L n v = 2 i=1 n l=1 l j=-l (v|Y (i) lj ) n Y (i) lj .
In a first step, the operator C κ is approached by

C κ,n j s ( x) = S 2 1 |η -z| L n {T x R(q; η, •)T x W 1 (q; κ, η, •)T x j s (•)}( z)ds( z) + i S 2 L n {T x W 2 (q; κ, η, •)T x j s (•)}( z)ds( z) + 1 κ curl S 2 S 2 1 |η -z| n {T x R(q; x, •)T x S 1 (q; κ, x, •)T x div S 2 j s (•) }( z)ds( z) + i κ curl S 2 S 2 n {T x S 2 (q; κ, x, •)T x div S 2 j s (•) }( z)ds( z).
for some n = an + 1 with fixed a > 1 and n -n > 3 (see [START_REF]A high-order tangential basis algorithm for electromagnetic scattering by curved surfaces[END_REF] and [8, Appendix A]). By the use of the fact that the scalar spherical harmonics are eigenfunctions of the single layer potential on the sphere and additional identities [START_REF]Inverse acoustic and electromagnetic scattering theory[END_REF] we obtain

C κ,n j s ( x) = W κ,n j s ( x) + 1 κ curl S 2 S κ,n div S 2 j s ( x), (4.2)
with

W κ,n j s ( x) = 2n +1 ρ =0 n +1 τ =1 µ ρ ν τ α τ T x R(q; η, z ρ τ )T x W 1 (q; κ, η, z ρ τ )T x j s ( z ρ τ ) + i 2n +1 ρ =0 n +1 τ =1 µ ρ ν τ T x W 2 (q; κ, η, z ρ τ )T x j s ( z ρ τ ), and
S κ,n ϕ s ( x) = 2n +1 ρ =0 n +1 τ =1 µ ρ ν τ α τ T x R(q; η, z ρ τ )T x S 1 (q; κ, η, z ρ τ )T x ϕ s ( z ρ τ ) + i 2n +1 ρ =0 n +1 τ =1 µ ρ ν τ T x S 2 (q; κ, η, z ρ τ )T x ϕ s ( z ρ τ ),
where α τ = n l=0 P l (ζ τ ). We proceed in the same way for the operator M κ as for the weakly singular part of C κ and we denote by M κ,n its approximation. Our fully discrete scheme for (3.6) is as follows: compute j n ∈ T n such that

j n -L n M κ,n j n -iηL n C κ,n Λ s j n = 2L n f s , (4.3)
where Λ s is defined as in remark 3.2. Since any density j n ∈ T n can be written

j n = n l=1 l j=-l α l,j Y (1)
l,j + β l,j Y (2)
l,j , the equation (4.3) is equivalent to a system of 2 × (n + 1) 2 -1 equations for the 2 × (n + 1) 2 -1 unknown coefficients (α l,j ) and (β l,j ) by applying the scalar product

( • |Y (1) l,j ) n and ( • |Y (2)
l,j ) n , for l = 1, . . . , n and j = -l, . . . , l to the equation (4.3). The same remarks hold true for the equation (3.7).

In the remaining of the section we briefly describe the fully discrete scheme for the equation (4.3) and give numerical results. The discrete approximation W κ of the weakly singular part of the operator of C κ is of the form

W κ = W 1,1 W 1,2 W 2,1 W 2,2 ,
where W a,b , for a, b = 1, 2 is a (n + 1) 2 -1 × (n + 1) 2 -1 matrix. The coefficients of W a,b , for 1 ≤ l, l ≤ n, |j| ≤ l and |j | ≤ l are given by

W ljl j a,b = (L n W κ,n Y (b) l,j |Y (a) l j ) n .
The discrete approximation of the operator M κ is obtained in the same way and is denoted by M κ . We refer to the paper of Hohage and Le Louër [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF] for details on their numerical implementation. The discrete approximation of the hypersingular part of the operator C κ is of the form

H κ = 0 0 H 2,1 0 , where T 2,1 is a (n + 1) 2 -1 × (n + 1) 2 -1 matrix. The coefficients of H 2,1 , for 1 ≤ l, l ≤ n,
|j| ≤ l and |j | ≤ l are given by

H ljl j 2,1 = κ -1 (L n curl S 2 S κ,n div S 2 Y (1) l,j |Y (2) l j ) n = -κ -1 (curl S 2 n S κ,n l(l + 1)Y l,j , Y (2)
l j ) n = -κ -1 ( n S κ,n l(l + 1)Y l,j , l (l + 1)Y l j ) n .
where S κ,n is the discrete approximation of the acoustic single layer potential and we refer to the paper of Ganesh and Graham [START_REF] Ganesh | A high-order algorithm for obstacle scattering in three dimensions[END_REF] for details on its numerical implementation. The discrete approximation of the operator

C κ is C κ = W κ + H κ . The discrete approximation of Λ s is of the form Λ n = 0 -I C (n+1) 2 -1 I C (n+1) 2 -1 0 .
Finaly, the discrete approximation of the integral operator associated to (2.4) is I-M κ -iηC κ Λ n , and the discrete approximation of the integral operator associated to (2.5) is

I + M κ + iηΛ n C κ .
The discrete approximation of the right-hand side 2f s = 2P q f of one of the boundary integral equations is the vector 2f , where f = (f 1 , f 2 ) T , whose coefficients are given for k = 1, 2, l = 1, . . . , n and j = -l, . . . , l by

f lj k = (L n P q f |Y (k) lj ) n . The parametrized form of the far-field operator F ∞ is F ∞ j s = F 1 j s + iηF 2 Λ s j s , with (F 1 j s ) ( x) = iκ 4π S 2 e -iκ x•q( y) e θ ( y) • j s ( y) x × t 1 ( y) + e φ ( y) • j s ( y) x × t 2 ( y) ds( y), and
(F 2 m s ) ( x) = κ 4π S 2 e -iκ x•q( y) e θ ( y) • m s ( y) x × t 1 ( y) + e φ ( y) • m s ( y) x × t 2 ( y) × x ds( y).
The discrete approximation of F ∞ evaluated at the 2(n ∞ + 1) 2 Gauss-quadrature points on the unit far sphere is

F ∞ = F 1 F 2 I C 2(n+1) 2 -2 0 0 Λ n with F 1 = F 1,1 F 1,2 , F 2 = F 2,1 F 2,2
where

F a,b , for a, b = 1, 2 is a 6(n ∞ + 1) 2 × ((n + 1) 2 -1) matrix. The coefficients of F a,b , for a, b = 1, 2, 1 ≤ l ≤ n, |j | ≤ l and ρ = 0, . . . , 2n ∞ + 1 and τ = 1, . . . , n ∞ + 1 are given by F ρτ l j a,b = F a Y (b) l j ( x ρτ ).
Table 2 exhibits emphasized fast convergence for the far-field pattern E ∞ for boundaries with parametric representations given in Table 1. As a first test we compute the electric far-field denoted, E ∞ ps , created by an off center point source located inside the perfect conductor :

E inc (x) = grad Φ(κ, |x -s|) × p, s ∈ Ω and p ∈ S 2 .
Table 1: Parametric representation of the perfectly conducting obstacles [START_REF] Andersson | Stellations, compounds, periodicity and the exponential scale[END_REF][START_REF] Onaka | Simple equations giving shapes of various convex polyhedra: The regular polydedra and polyhedra composed of crystallographically low-index planes[END_REF] surface parametric representation sphere q • ψ(θ, φ) = Rψ(θ, φ), R ∈ R * + ; great stellated dodecahedron q • ψ(θ, φ) = 0.33 r(θ, φ)ψ(θ, φ), where r(θ, φ) > 0 solves

p -|r(θ,φ)f (δ,ξ,0)| p + p -|r(θ,φ)f (δ,-ξ,0)| p + p -|r(θ,φ)f (0,δ,ξ)| p + p -|r(θ,φ)f (0,δ,-ξ)| p + p -|r(θ,φ)f (ξ,0,δ)| p + p -|r(θ,φ)f (-ξ,0,δ)| p = 2.5, with f (a, b, c) = a sin θ cos φ + b sin θ sin φ + c cos θ, δ = 5- √ 5 10 , ξ = 5+ √ 5 10 , p ∈ N, p > 2.
Table 2: Numerical examples for the perfect conductor problem In this case the total exterior wave has to vanish so that the far-field pattern of the scattered wave E s is the opposite of the far field pattern of the incident wave:

surface n ||[E ∞ ps ]n -E ∞ exact ||∞ Re[E ∞ pw (d)]n • p Im[E ∞ pw (d)]n
E ∞ exact ( x) = - iκ 4π e -iκ x•s ( x × p).
We choosed s = 0, 0.1

√ 2 , -0.1 √ 2
T and p = (1, 0, 0) T . In the tabulated results we indicate the uniform-norm error (by taking the maximum of errors obtained over 1300 observed directions, i.e n ∞ = 25) :

||[E ∞ ps ] n -E ∞ exact || ∞ = max x∈S 2 [E ∞ ps ] n -E ∞ exact .
As a second task we compute the electric far field denoted, E ∞ pw , created by the scattering of an incident plane wave :

E inc (x) = p e iκx•
Q = {q ∈ H s (Γ ref , R 3 ) : q injective, det(Dq( x)) = 0 for all x ∈ Γ ref }. For s > 2 the set Q is open in H s (Γ ref , R 3 ) and H s (Γ ref , R 3 ) ⊂ C 1 (Γ ref , R 3
). Since now, for q ∈ Q we define Γ q := q(Γ ref ) and denote by n q the exterior unit normal vector to Γ q . More generally we will label all quantities and operators related to the perfect conductor problem for the boundary Γ q by the index q. We denote by F k : Q → L 2 t (S 2 ) the boundary to far field operator that maps the parametrization q of the boundary Γ q onto the far-field pattern E ∞ q,k corresponding to the scattering of the incident wave E inc k by the perfect conductor Γ q . These operators may be combined into one operator F : Q → L 2 t (S 2 ) m , F (q) := (F 1 (q), . . . , F m (q)) T . We also combine the measured far-field patterns into a vector

E ∞ δ := (E ∞ 1,δ , . . . , E ∞ m,δ ) T ∈ L 2 t (S 2
) m such that the inverse problem (1.2) can be written as the following nonlinear equation

F (q) = E ∞ δ . (5.1)
A regularized iterative method can then be applied to the linearized equation

F [q]ξ = E ∞ δ -F [q]
. Here, we use the Iteratively Regularized Gaus-Newton Method (IRGNM). This requires to prove the Fréchet differentiability of the boundary to far-field operator F and an explicit form of the first Fréchet dérivative F [q] and its adjoint F [q] * . Then the iterates of the IRGNM can be computed by

q δ N +1 := argmin q∈H s (S 2 ,R) ||F [q δ N ](q -q δ N ) + F (q δ N ) -E ∞ δ || 2 L 2 + α N ||q -q 0 || 2 H s , (5.2)
Here q 0 = q δ 0 is some initial guess and the regularization parameters are chosen of the form α N = α 0 2 3 N . The updates (∂q) N := q δ N +1 -q δ N are the unique solutions to the linear equations

α N I + F [q δ N ] * F [q δ N ] (∂q) δ N = F [q δ N ] * E ∞ κ,δ -F (q δ N ) + α N q δ 0 -q δ N . (5.3)
For the analysis of the derivative and its adjoint, we restrict ourselves to the case m = 1 since the general case can be reduced to this special case by the obvious formulas

F [q]ξ = (F 1 [q]ξ, . . . , F m [q]ξ) T and F [q] * h = m k=1 F k [q] * h.
The following theorem is a rewriting, for the electric field only, of the theorem established in [START_REF] Kress | Electromagnetic waves scattering : Scattering by obstacles[END_REF] by Kress.

Theorem 5.1 (characterization of F [q]). The mapping F : Q → L 2 t (S 2 ) with s > 2 is Fréchet differentiable at all q ∈ Q for which Γ q is of class C 2 , and the first derivative at q in the direction ξ

∈ H s (Γ ref , R 3 ) is given by F [q]ξ = E ∞ q,ξ
, where E ∞ q,ξ is the far-field pattern of the solution E s q,ξ to the Maxwell equation (1.1a) in Ω c q that satisfies the Silver-Müller radiation condition and the following Dirichlet boundary condition

n q × E s q,ξ = -ξ • q -1 • n q u q × n q - 1 κ 2 curl Γq (ξ • q -1 • n q ) div Γq u q .
(5.4) on Γ q where u q = n q × curl E s q + E inc and E s q is the solution of the scattering problem (1.1a)-(1.1c) satisfying the Dirichlet boundary condition

n q × E s q + E inc = 0 on Γ q .
To define the adjoint of

F [q] : H s (Γ ref ; R 3 ) → L 2 t (S 2
), we interpret the naturally complex Hilbert space L 2 t (S 2 ) as a real Hilbert space with the real-valued inner product Re • , • L 2 t (S 2 ) . For bounded linear operator between complex Hilbert spaces such a reinterpretation of the spaces as real Hilbert spaces does not change the adjoint.

Proposition 5.2 (characterization of the adjoint F

[q] * ). Let E inc h (y) = 1 4π S 2 e -iκ x•y h( x) ds( x), y ∈ R 3
be the vector Herglotz function with kernel h ∈ L 2 t (S 2 ) and E q, h the total wave solution to the scattering problem for the perfectly conducting interface Γ q and the incident wave E inc h . Moreover, let j H s →L 2 denote the embedding operator from H

s (Γ ref , R 3 ) to L 2 (Γ ref , R 3 ). Then F [q] * h = j * H s →L 2 J q n q Re -n q × curl E q, h • u q + 1 κ 2 div Γq n q × curl E q, h • div Γq u q • q .
Proof. It suffices to follow the three-step proof of Proposition 6.3 in [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF].

Remark 5.3 As mentioned in [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF], the implementation of the formulas in Theorem 5.1 and Proposition 5.2 is straightforward using our discretization. Indeed, by the use of the Piola transform it remains to compute some surface derivatives on S 2 and both div S 2 and curl S 2 are diagonal with respect to the chosen bases of spherical harmonics and vector spherical harmonics.

For the numerical experiments, we consider perfectly conducting scatterers that are starshaped with respect to the origin. In this case, it is natural to choose the unit sphere as a reference boundary (Γ ref = S 2 ) and to consider special parametrizations of the form q = Rr where Rr ( x) = r( x) x for all x ∈ S 2 and r is a positive-valued function defined on the unit sphere. The function r is uniquely determined by Γ q . We restrict the set of admissible parametrization to radial parametrization and choose Q star = {r ∈ H s (S 2 , R) : r > 0}. We have R(Q star ) ⊂ Q, and we can define F star : Q star → L 2 t (S 2 ) m by F star := F • R. Then F star is injective if a starshaped interface Γ is uniquely determined by the far field data E ∞ 1 , . . . , E ∞ m . Using [START_REF] Harbrecht | Fast methods for three-dimensional inverse obstacle scattering problems[END_REF]Corollary 4] we obtain that F star [r] * h = j H s →L 2 r 2 Re{. . .} • q where the expression in the curly brackets coincides with that in Proposition 5.2. In practice, we apply the IRGNM to the nonlinear equation F star (r) = E ∞ δ . We refer to [START_REF] Hohage | A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem[END_REF] for details on the inverse scattering algorithm. Our tests concern the reconstruction of the (rounded) stellation of the dodecahedron defined by the implicit equation given in table 1. This shape consists of 20 peaks located at the 20 vertex of the regular dodecahedron. The diameter of the obstacle is roughly 2. To compute the exact far-field data we use the boundary integral equation (2.4) with n = 30. To compute the far-field data at each iteration step we use the boundary integral equation (2.5) with n = 25. This allows us to retrieve the boundary value data of the scattered field that are needed to evaluate the adjoint. To compute the Fréchet derivatives of the boundary to far-field operator we use the boundary integral equation (2.4). In each of the following test, we use a mesh grid of 128 (n = 7) degrees of freedom to evaluate the far-field pattern. The radial functions describing the reconstruction belong to the space of spherical harmonics of order less than or equal to 20. The initial guess is the unit sphere. We choose α 0 = 10 -2 .

Test 1: We reconstruct the obstacle from noisy far-field measurements corresponding to the scattering of 6 incident plane waves from top, bottom, front, back, left and right. The noise level is 5%. The wavenumber κ is chosen such that the diameter of the obstacle is equivalent to the wavelength, that is κ = π. After ten iterations we obtain the result, presented in the top right picture below, which is pretty good. We observe that all the peaks are recovered.

Test 2: We reconstruct the obstacle from noisy far-field measurements corresponding to the scattering of one incident plane wave from bottom and polarized in the (Ox) direction. The other parameters does not change. After twenty iterations, we obtain the result presented in the bottom left picture below. Test 3: We reconstruct the obstacle from noisy far-field measurements corresponding to the scattering of 6 incident plane waves. The noise level is 1%. Here The wavenumber κ is chosen such that the diameter of the obstacle is greater than the wavelength, that is κ = π 2 . After twenty iterations, we obtain the result presented in the bottom right picture below.

  curl x Φ(κ, • -y)j(y) ds(y) and Dj = Γ curl x Φ(κ, • -y)j(y) ds(y) satisfy the Maxwell equation and the Silver-Müller radiation condition.

- 1 2 1 2 1 2

 111 div (S 2 ) obtained by applying the Piola transform (3.1) to the unknown in (2.4) or (2.5).Lemma 3.1The operator Λ is an elliptic operator on H div (Γ) for the bilinear form (2.3) if and only if Λ s is an elliptic operator on H div (S 2 ) for the bilinear form

  d , where d, p ∈ S 2 and d • p = 0.In the tabulated results we indicate the real part and the imaginary part of the polarization component of the electric far field evaluated at the incident direction :[E ∞ pw (d)] n • p. We choosed d = 0, 0, 1T and p = (1, 0, 0) T . In each of the examples we take η = 1.5 Numerical solution of the inverse problem and IRGNMWith the objective of using a regularized Newton-type method to solve numerically the inverse problem (1.2) we first reformulate it as a nonlinear equation posed on an open subset of admissible parametrizations of a Hilbert space. More precisely, we choose some reference domain Ω ref with a simply connected closed boundary Γ ref and consider mappings q : Γ ref → Γ belonging to

  Figure. Reconstruction of the (rounded) great stellated dodecahedron.

Acknowledgment: We would like to thank Olha Ivanyshyn and Rainer Kress for matlab code and for helpful discussions. Financial support of DFG through CRC 755 is gratefully acknowledged.