Keywords: Calcul matriciel, Calcul parallèl et distribué Physics of waves, Numerical methods, Matrix computation, Large-scale parallel and distributed computing

Titre : Diffraction électromagnétique par des réseaux et des surfaces rugueuses aléatoires: Mise en oeuvre de méthodes hautement efficaces pour la résolution de systèmes aux valeurs propres et de problèmes aux conditions initiales Mots clés : Physique des ondes, Méthodes numériques,

List of Figures

 [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF]to calculate 3 eigenvalues with largest magnitude, the result is compared with IRAM (16) . . . . . . . . . . . . . . . . . . . 7.9 M IRAM ns [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] to calculate 3 eigenvalues with smallest magnitude, the result is compared with IRAM (16) . . . . . . . . . . . . . . . . . . . 7.10 SS-RR method with Center = 12, Radius = 4, 364 eigenvalues are calculated, average residual is 3.53698402454027090 ⇥ 10 6 . . . . . . . . . . . 7.11 SS-RR method with Center = 6, Radius = 2, 136 eigenvalues are calculated, average residual is 5.78977861683748897 ⇥ 10 6 . . . . . . . . . . . 7.12 SS-RR method with Center = 3, Radius = 1, 95 eigenvalues are calculated, average residual is 2.18255949228685363 ⇥ 10 4 . . . . . . . . . . . 7.13 SS-RR method with Center = 1.5, Radius = 0.5, 72 eigenvalues are calculated, average residual is 2.59612185434167264 ⇥ 10 6 . . . . . . . . 7.14 SS-RR method with Center = 0.75, Radius = 0.25, 57 eigenvalues are calculated, average residual is 1.45660089462261519 ⇥ 10 6 . . . . . . . . 7.15 SS-RR method with Center = 0.3, Radius = 0.2, 95 eigenvalues are calculated, average residual is 2.59080105213133636 ⇥ 10 4 . . . . . . . . 7.16 SS-RR method with Center = 0.055, Radius = 0.045, 80 eigenvalues are calculated, average residual is 1.40601714058967572 ⇥ 10 5 . . . . . . . . 7.17 SS-RR method with Center = 0.0055, Radius = 0.0045, 45 eigenvalues are calculated, average residual is 2.61507800766467365 ⇥ 10 6 . . . . . . 7.18 SS-RR method with Center = 0.00055, Radius = 0.00045, [START_REF] Tsang | Wave scattering with the UV multilevel partitioning method: 2. three-dimensional problem of nonpenetrable surface scattering[END_REF]

Introduction

Electromagnetic wave scattering is an active and interdisciplinary area of research with myriad practical applications in fields ranging from atomic physics to optics, medical imaging, geoscience and remote sensing [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF][START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Stratton | Electromagnetic Theory. International series in physics[END_REF][START_REF] Ishimaru | Electromagnetic wave propagation, radiation and scattering[END_REF][START_REF] Palmer | Di↵raction grating handbook[END_REF][START_REF] Petit | Electromagnetic Theory of Gratings[END_REF]. In particular, the subject of wave scattering by gratings [START_REF] Palmer | Di↵raction grating handbook[END_REF][START_REF] Petit | Electromagnetic Theory of Gratings[END_REF][START_REF] Petit | Sur la di↵raction d'une onde plane par un rseau infiniment conducteur[END_REF][START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface[END_REF][START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Van Den Berg | The Rayleigh hypothesis in the theory of reflection by a grating[END_REF][START_REF] Keller | Singularities and Rayleigh's hypothesis for di↵raction gratings[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF] and rough surfaces [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF][START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF][START_REF] Bourlier | Azimuthal harmonic coe cients of the microwave backscattering from a non-Gaussian ocean surface with the first-order SSA model[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF][START_REF] Voronovich | Wave scattering from rough surfaces[END_REF][START_REF] Warnick | Numerical simulation methods for rough surface scattering[END_REF][START_REF] Saillard | Rigorous solutions for electromagnetic scattering from rough surfaces[END_REF][START_REF] Pak | Numerical simulations and backscattering enhancement of electromagentic waves from two-dimensional dielectric random rough surfaces with the sparse-matrix canonical grid method[END_REF][START_REF] Wagner | Monte Carlo simulations of electromagentic scattering from two-dimensional random rough surfaces[END_REF][START_REF] Jandhyala | Fast algorithm for the analysis of scattering by dielectric rough surfaces[END_REF][START_REF] Torrungrueng | Numerical studies of backscattering enhancement of electromagnetic waves from two-dimensional random rough surfaces with the forward-backward novel spectral acceleration method[END_REF][START_REF] Soriano | Scattering of electromagentic waves from twodimensional rough surfaces with an impedance approximation[END_REF][START_REF] Tsang | Wave scattering with the UV multilevel partitioning method: 2. three-dimensional problem of nonpenetrable surface scattering[END_REF] presents great theoretical challenges due to the large number of degrees of freedom in these systems and a need to include multiple scattering e↵ects accurately. In the past several decades, considerable theoretical progress has been made in elucidating and understanding the scattering processes involved in such problems. Diagrammatic techniques and e↵ective medium theories remain essential for analytical studies; however, rapid advances in computer technology have opened new doors for researchers with the full power of Monte Carlo simulations in the numerical analysis of random media scattering [START_REF] Warnick | Numerical simulation methods for rough surface scattering[END_REF][START_REF] Saillard | Rigorous solutions for electromagnetic scattering from rough surfaces[END_REF][START_REF] Pak | Numerical simulations and backscattering enhancement of electromagentic waves from two-dimensional dielectric random rough surfaces with the sparse-matrix canonical grid method[END_REF][START_REF] Wagner | Monte Carlo simulations of electromagentic scattering from two-dimensional random rough surfaces[END_REF][START_REF] Jandhyala | Fast algorithm for the analysis of scattering by dielectric rough surfaces[END_REF][START_REF] Torrungrueng | Numerical studies of backscattering enhancement of electromagnetic waves from two-dimensional random rough surfaces with the forward-backward novel spectral acceleration method[END_REF][START_REF] Soriano | Scattering of electromagentic waves from twodimensional rough surfaces with an impedance approximation[END_REF][START_REF] Tsang | Wave scattering with the UV multilevel partitioning method: 2. three-dimensional problem of nonpenetrable surface scattering[END_REF]. Numerical simulations allow us to solve Maxwell's equations without the limitations of analytical approximations, whose regimes of validity are often di cult to assess [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF][START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF][START_REF] Bourlier | Azimuthal harmonic coe cients of the microwave backscattering from a non-Gaussian ocean surface with the first-order SSA model[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF][START_REF] Voronovich | Wave scattering from rough surfaces[END_REF][START_REF] Born | Principles of optics Electromagnetic theory of propagation, Interference and di↵raction of light[END_REF].

In this thesis, we study the electromagnetic di↵raction by gratings and random rough surfaces. The C-method is an exact method developed for this aim. It is based on Maxwells equations under covariant form written in a nonorthogonal coordinate system [START_REF] Post | Formal structure of electromagnetic[END_REF][START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF][START_REF] Chandezon | Application d'une nouvelle méthode de résolution des équations de Maxwell à l'étude de la propagation des ondes électromagnétiques dans les guides périodiques[END_REF]. Discretizing the Maxwells equations under the non-orthogonal coordinate system and separating variables lead to solving the eigenvalue problem of the high dimension, dense and non-symmetric scattering matrix. All the eigenvalues and eigenvectors of the scattering matrix are needed. The scattered field is expanded as a linear combination of eigensolutions satisfying the outgoing wave condition. The boundary conditions allow the di↵raction amplitudes to be determined. This method has been used for analyzing gratings used in optics ,waveguides [START_REF] Dusséaux | Etude de transformateur plan-E dans un système de coordonnées non orthogonales[END_REF][START_REF] Dusséaux | Analysis of rectangular waveguide Hplane junctions in nonorthogonal coordinate system[END_REF][START_REF] Dusséaux | Analyse de composants plan-E symétriques en guides d'onde à section rectangulaire[END_REF][START_REF] Dusséaux | Telegraphist's equations for rectangular waveguides and analysis in nonorthogonal coordinates[END_REF] and rough surfaces [START_REF] Benali | A new theory for scattering of electromagentic waves from conducting or dielectric rough surfaces[END_REF][START_REF] Dusséaux | Scattering of a plane wave by one-dimensional dielectric rough surfaces-study of the field in a nonorthogonal coordinate system[END_REF][START_REF] Baudier | Scattering of a plane wave by one-dimensional dielectric random rough surfaces -study with the curvilinear coordinate method[END_REF][START_REF] Braham | Scattering of electromagnetic waves from two-dimensional perfectly conducting random rough surfaces -study with the curvilinear coordinate method[END_REF][START_REF] Dusséaux | Implementation and validation of the curvilinear coordinate method for the scattering of electromagnetic waves from two-dimensional dielectric random rough surfaces[END_REF][START_REF] Braham | The curvilinear coordinate method associated with the short-coupling-range approximation for the study of scattering from onedimensional random rough surfaces[END_REF][START_REF] Dusséaux | Study of backscatter signature for seedbed surface evolution under rainfall -influence of radar precision[END_REF][START_REF] Edee | Beam implementation in a nonorthogonal coordinate system: application to the scattering from random rough surfaces[END_REF][START_REF] Dusséaux | Eigenvalue system for the scattering from rough surfaces-saving in computation time by a physical approach[END_REF][START_REF] Prémel | Computation of a quasi-static field induced by two long straight parallel wires in a conductor with a rough surface[END_REF][START_REF] Prémel | Generalization of the second order vector potential formulation for arbitrary non-orthogonal curvilinear coordinates systems from the covariant form of Maxwell's equations[END_REF] We focus on the numerical aspect of the C-method, trying to develop an e cient implementation of this exact method. Iterative eigensolvers, such as Krylov subspace methods or Jacobi-Davidson methods [START_REF] Bai | Templates for the solution of algebraic eigenvalue problems[END_REF] have been developed to deal with large-scale eigenvalue problems. However, they have the possibility of missing some eigenvalues. So the standard iterative methods are ine↵ective for the C-method because all the eigenvalues and eigenvectors are needed. In contrast, the QR algorithm, which is based on similarity transformations, calculates all the eigenvalues and eigenvectors with very little danger and only with a warning of missing some eigensolutions. We propose a specifically designed parallel QR algorithm for the C-method to solve the eigenvalue problem.

This parallel QR algorithm is a variant of QR algorithm based on three techniques:

early shift, parallel bulge chasing [START_REF] Braman | The multishift QR algorithm. part I: Maintaining well-focused shifts and level 3 performance[END_REF][START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF] and parallel aggressive early deflation (AED) [START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF][START_REF] Braman | The multishift QR algorithm. part II: Aggressive early deflation[END_REF]. We propose the "early shift" for the scattering matrix according to the property we have observed. That is the C-method and the physical interpretation behind the Cmethod provides very good approximations of some eigenvalues before any calculations.

The "early shift" provides the possibility of quick deflation. We mixed the "early shift", Wilkinson's shift and exceptional shift together to accelerate the convergence. Especially, we use the "early shift" to have quick deflation of the approximated eigenvalues of the scattering matrix. They provide the possibility of quick deflation. For the bulge chasing, instead of only a single bulge, containing two shifts, a chain of several tightly coupled bulges, each containing two shifts, is chased in the course of one multishift QR iteration. This idea and the delay-and-accumulate technique [START_REF] Braman | The multishift QR algorithm. part I: Maintaining well-focused shifts and level 3 performance[END_REF][START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF] allow performing most of the computational work in terms of matrix-matrix multiplications to benefit from level 3 Basic Linear Algebra Subprograms (BLAS, the level 3 contains matrixmatrix operations) [70]. Aggressive early deflation is a QR algorithm deflation strategy that takes advantage of matrix perturbations outside of the subdiagonal entries of the Hessenberg QR iterate. It identifies and deflates converged eigenvalues long before the classic small-subdiagonal strategy would. Aggressive early deflation can significantly enhance the convergence of the QR algorithm.

We also propose a spectral projection method to solve the eigenvalue problem e ciently.

We propose a global eigensolver by a combination of the SS method (Sakuria and Sugiura method, proposed by Sakuria and Sugiura in [START_REF] Sakuria | A projection method for generalized eigenvalue problem using numerical integration[END_REF]) and MIRAMns (Multiple Implicitly Restarted Arnolid Method with nested subspaces, proposed by S. A. Shahzadeh Fazeli et al. in [START_REF] Shahzadeh Fazeli | A key to choose subspace size in implicitly restarted arnoldi method[END_REF]). This proposed global eigensolver allows us to calculate a large number (or all) of the eigenvalues of a general matrix. According to our experiments presented in chapter 7, this method has the advantage of having very good scalability compared to the QR algorithm. This promising method can be continued in future work.

The original C-method is not very e cient when we are dealing with multilayer gratings.

We want to find other solutions. Especially, we want to explore the potential parallelization of multilayer gratings. So, we propose a new version of C-method which leads to a di↵erential system with initial conditions for gratings. We show that this new version of C-method can be used to study multilayer gratings with homogeneous medium. We

show that this formulation is an interesting tool for analyzing perfectly conducting or dielectric gratings with deep grooves. The proposed method allows analyzing the complex phenomenon of incident energy absorption. We apply this method to multilayer gratings with an arbitrary number of interfaces. We show how to combine the local scattering matrix to obtain the global one. We validate our method by comparing experiments results with that from published paper. We show that this new version of C-method has very good accuracy as well as a nature of two level parallel property [START_REF] Pan | Curvilinear coordinate method as an initial value problem: application to gratings[END_REF][START_REF] Pan | The c-method as an initial value problem: application to multilayer gratings[END_REF]. This new version of C-method is an attractive alternative to analyze multilayered grating having parallel or non-parallel interfaces.

During my Ph.D study, I worked in three laboratories: Laboratoire Atmosphères, Milieux, Observations Spatiales (LATMOS), Maison de la Simulation (MDLS) and Laboratoire Parallélisme, Réseaux, Systèmes, Modélisation (PRiSM). LATMOS is a laboratory specializing on the fundamental physical and chemical processes of atmospheres, environments, and spatial observations. MDLS is a laboratory specializing on scientific computing and simulations using HPC. PRiSM is a laboratory specializing on computer science.

This thesis is structured as follows. In chapter 2, we present the electromagnetic field theory fundamentals. In chapter 3, we present the scattering problem by gratings and by random rough surfaces. In chapter 4, we present the curvilinear coordinate method.

In chapter 5, we propose the parallel QR algorithm for the C-method. In chapter 6, we present an implementation of our parallel QR algorithm and present the results of the numerical experiments. In chapter 7, we propose an alternative to QR algorithm for solving the eigenvalue problem. In chapter 8, we propose a new version of C-method, which is the C-method as an initial value problem.

Chapter 2

Electromagnetic field theory fundamentals

Maxwell's equations

The general theory of electromagnetic phenomena is based on Maxwell's equations, which constitute a set of four coupled first order vector partial di↵erential equations relating the space and time changes of electric and magnetic fields to their scalar source densities (divergence) and vector source densities (curl) [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF]. Maxwell's equations are usually formulated in di↵erential form (i.e., as relationships between quantities at the same point in space and at the same instant in time) or in integral form where, at a given instant, the relations of the fields with their source are considered over an extensive region of space [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF]. The two formulations are related by the divergence and Stokes' theorems.

For stationary media, Maxwell's equations in di↵erential and integral forms are:

Di↵erential form of Maxwell's equation r • D(r, t) = ⇢(r, t) (Gauss 0 law) (2.1) r • B(r, t) = 0 (Gauss 0 law f or magnetic f ields) (2.2) r ⇥ Ẽ(r, t) = @ B(r, t) @t (F araday 0 s law) (2.3) r ⇥ H(r, t) = J(r, t) + @ D(r, t) @t (Generalized Ampere 0 s law) (2.4)
where r is the Nabla operator, r • D designates the divergence of D and r ⇥ Ẽ, the curl of Ẽ.

Integral form of Maxwell's equations

I S D(r, t) • ds = Q T (t) (Gauss 0 law) (2.5) I S B(r, t) • ds = 0 (Gauss 0 law f or magnetic f ields) (2.6) I Ẽ(r, t) • d l = Z S @ B(r, t) @t • ds (F araday 0 s law) (2.7) I H(r, t) • d l = Z S ( J(r, t) + @ D(r, t) @t ) • ds (Generalized Ampere 0 s law) (2.8)
where S is any fixed open surface and is the associated boundary curve.

Maxwell's equations, involve only macroscopic electromagnetic fields and, explicitly, only macroscopic densities of free-charge, ⇢(r, t), which are free to move within the medium, giving rise to the free-current densities, J(r, t). The e↵ect of the macroscopic charges and current densities bound to the medium's molecules is implicitly included in the auxiliary magnitudes D and H which are related to the electric and magnetic fields, Ẽ and B by the so-called constitutive equations that describe the behavior of the medium.

In general, the quantities in these equations are arbitrary functions of the position r and the time t. The definition and units of these quantities are: The equations (2.1)-(2.4) or (2.5)-(2.8) as a whole are associated with the name of Maxwell's equations because he was responsible for completing them, adding to Ampere's original equation, r ⇥ H(r, t) = J(r, t), the displacement current density term or, in short, the displacement current, @ D @t , as an additional vector source for the field H. This term has the same dimensions as the free current density but its nature is different because no free charge movement is involved. Its inclusion in Maxwell's equation is fundamental to predict the existence of electromagnetic wave which can propagate through empty space at the constant velocity of light c. The concept of displacement is also fundamental to deduce from equation (2.4) the principle of charge conservation by means of the continuity equation:

Ẽ = electric f ield intensity (volt/meter, V • m 1 ) B = magnetic f lux density (weber/square meter, W b • m 2 ) D = electric f lux density (coulomb/square meter, C • m 2 ) H = magnetic f ield density (ampere/meter, A • m 1 ) ⇢ = free electric charge density (coulomb/cubic meter, C • m 3 ) Q T = net f
r • J = @⇢ @t (2.9)
or, in integral form:

I J • ds = dQ t dt (2.10)
With these equations, Maxwell showed not only that the electric and magnetic fields are interrelated but also that they are in fact two aspects of a single concept, the electromagnetic field.

Maxwell's equations together with the Lorenz's force constitute the basic mathematical formulation of the physical laws that at a macroscopic level explain and predict all the electromagnetic phenomena which basically comprise the remote interaction of charges and currents taking place via the electric and/or magnetic fields that they produce.

In applications, Maxwell's equations have to be complemented by appropriate initial and boundary conditions. The initial conditions involve values or derivatives of the fields at t = 0, while the boundary conditions involve the values or derivatives of the fields on the boundary of the spatial region of interest. Usually, we consider the initial conditions as a form of boundary conditions and refer to the solution of Maxwell's equation, with all these conditions, as a boundary-value problem.

Next, we briefly describe the physical meaning of Maxwell's equations.

Gauss' law is a direct mathematical consequence of Coulomb's law, which states that the interaction force between electric charges depends on the distance r, between them, as r 2 . According to Gauss's law, the divergence of the vector field D is the volume density of free electric charges which are the sources or sinks of the field D, i.e. the lines of D begin on positive charges (⇢ > 0) and end on negative (⇢ < 0). In its integral form, Gauss' law relates the flux of the vector D through a closed surface S to the total free charge within that surface.

Gauss' law for magnetic fields states that the B field does not have scalar sources, i.e., it is divergenceless or solenoidal. This is because no free magnetic charges or monopoles have been found in nature which would be the magnetic analogues of electric charges for Ẽ. Hence, there are no sources or sinks where the field lines of B start or finish, i.e., the field lines of B are closed. In its integral form, this fact indicates that the flux of the B field through any closed surface S is null.

Faraday's law establishes that a time-varying B field produces a non conservative electric field whose field lines are closed. In its integral form, Faraday's law states that the time variation of the magnetic flux ( R B • ds) through any surface S bounded by an arbitrary closed loop , induces an electromotive force given by the integral of the tangential component of the induced electric field around . The line integration over the contour must be consistent with the direction of the surface vector ds according to the righthand rule. The minus sign in the equations of the law represents the feature by which the induced electric field, when it acts on charges, would produce an induced current that opposes the change in the magnetic flux (Lenz's law).

Ampere's generalized law, constitutes another connection, di↵erent from Faraday's law, between Ẽ and B. It states that the vector sources of the magnetic field may be free currents, J, and/or displacement currents, @ D @t . Thus, the displacement current performs, as a vector source of H, a similar role to that played by @ B @t as a source of Ẽ. In its integral form, the left-hand side of the generalized Ampere's law equation represents the integral of the magnetic field tangential components along an arbitrary closed loop and the right-hand side is the sum of the flux, through any surface S bounded by a closed loop , of both currents: the free current J and the displacement current @ D @t .

Constitutive equations

In the vacuum, Maxwell's equations can be written without using the artificial fields D and H, as

r • Ẽ(r, t) = ⇢ all " 0 (r, t) (2.11) r • B(r, t) = 0 (2.12) r ⇥ Ẽ(r, t) = @ B(r, t) @t (2.13) r ⇥ B(r, t) = µ 0 Jall (r, t) + µ 0 " 0 @ Ẽ(r, t) @t (2.14)
where " 0 = 10 9 36⇡ (farad/meter, F • m 1 ) and µ 0 = 4⇡10 7 (henry/meter, H • m 1 ) are two constants called electric permittivity and magnetic permeability of free space, respectively. The subscript all indicates that all kinds of charges (free and bound) must be individually included in ⇢ and J. These equations are, within the limits of classical electromagnetic theory, absolutely general. Nevertheless, in order to make it possible to study the interaction between an electromagnetic field and a medium and to take into account the discrete nature of matter, it is necessary to develop macroscopic models to obtain Maxwell's macroscopic equations, in which only macroscopic quantities are used and in which only the densities of free charges and currents explicitly appear as sources of the fields. To this end, the atomic and molecular physical properties, which fluctuate greatly over atomic distances, are averaged over microscopically large volume elements, v, so that these contain a large number of molecules but at the same time are macroscopically small enough to represent accurate spatial dependence at a macroscopic scale. As a result of this average, the properties of matter related to atomic and molecular charges and currents are described by the macroscopic parameters, electric permittivity ", magnetic permeability µ, and electrical conductivity . These parameters, called constitutive parameters, are in general smoothed point functions. The derivation of the constitutive parameters of a medium from its microscopic properties is, in general, an involved process that may require complex models of molecules as well as quantum and statistical theory to describe their collective behavior. Fortunately, in most of the practical situations, it is possible to achieve good results using simplified microscopic models.

To define the electric permittivity and describe the behavior of the electric field in the presence of matter, we must introduce a new macroscopic field quantity, P (C • m 2 ), called electric polarization vector, such that

D = " 0 Ẽ + P (2.15)
and defined as the average dipole moment per unit volume

P = lim v!0 P N v n=1 pn v (2.16)
where N is the number of molecules per unit volume and the numerator is the vector sum of the individual dipolar moments, pn , if atoms and molecules contained in a macroscopically infinitesimal volume v. For many materials, called linear isotropic media, P can be considered collinear and proportional to the electric field applied. Thus we have

P = " 0 e Ẽ (2.17)
where the dimensionless parameter e , called the electric susceptibility of the medium, describes the capability of a dielectric to be polarized. Equation (2.15) can be written in a more compact form as

D = (1 + e )" 0 Ẽ (2.18) so that D = " 0 " r Ẽ = " Ẽ (2.19)
where

" r = 1 + e (2.

20)

and

" = " 0 " r (2.21)
are the relative permittivity and the permittivity of the medium, respectively.

To define the magnetic permeability and describe the behavior of the magnetic field in the presence of magnetic materials, we must introduce another macroscopic field quantity, called magnetization vector M (A • m 1 ), such that

H = B µ 0 M (2.22)
where M is defined, in a similar way to that of the electric polarization vector, as the average magnetic dipole moment per unit volume

M = lim v!0 P N v n=1 mn v (2.23)
where N is the number of atomic current elements per unit volume and the numerator is the vector sum of the individual magnetic moments, mn contained in a macroscopically infinitesimal volume v.

In general, M is a function of the history of B or H, which is expressed by the hysteresis curve. Nevertheless, many magnetic media can be considered isotropic and linear, such that

M = m H (2.24)
where m is the dimensionless magnetic susceptibility magnitude, being negative and small for diamagnetic, positive and small for paramagnet, and positive and large for ferromagnet. Thus

H = 1 (1 + m )µ 0 B = 1 µ B (2.25)
where

µ r = (1 + m ) (2.26) and µ = µ r µ 0 (2.27)
are the relative magnetic permeability and the permeability of the medium, respectively.

In a vacuum, or free space, " r = 1, µ r = 1, and therefore the fields vectors D and Ẽ, as well as B and H, are related by Ohm's law may not be applicable. For most metals is a scalar with a magnitude that depends on the temperature and that, at room temperature, has a very high value of the order of 10 7 mho • m 1 . Then very often metals are considered as perfect conductors with an infinite conductivity.

D = " 0 Ẽ (2.28) B = µ 0 H (2.
The relations between macroscopic quantities, (2.15), (2.25) and (2.30), are called constitutive relations. Depending on the characteristics of the constitutive macroscopic parameters ", µ and , which are associated with the macroscopic response of atoms and molecules in medium, this medium can be classified as:

• Inhomogeneous or homogeneous: according to whether or not the constitutive parameter of interest is a function of the position, " = "(r), µ = µ(r), = (r).

• Anisotropic or isotropic: according to whether or not the response of the medium depends on the orientation of the field. In isotropic media all magnitudes of interest are parallel, i.e., Ẽ and D, Ẽ and Jc , B and H. In anisotropic materials the constitutive parameter of interest is a tensor.

• Non linear or linear: according to whether or not the constitutive parameters depend on the magnitude of the applied fields. For instance "(E), (E) and µ(H).

= (t).

• Dispersive: according to whether or not, for time-harmonic fields, the constitutive parameters depend on the frequency, F T ("

) = F T (")(!), F T (µ) = F T (µ)(!), F T ( ) =
F T ( )(!), here, F T () represents the Fourier transformation. The materials in which these parameters are functions of the frequency are called dispersive.

• Magnetic medium: if µ 6 = µ 0 . Otherwise the medium is called nonmagnetic because its only significant reaction to the electromagnetic field is polarization.

Fortunately, in many cases the medium in which the electromagnetic field exists can be considered homogeneous, linear, isotropic, non dispersive and non magnetic. Indeed, this assumption is not very restrictive since many electromagnetic phenomena can be studied using this simplification.

Equations (2.15) and (2.22) are simplified. In practice, we must write [START_REF] Jackson | Classical Electrodynamics[END_REF]:

D = D( Ẽ, B) (2.31)
H = H( Ẽ, B) (2.32)
In this thesis, we only consider:

• Non magnetic medium for which,

H(r, t) = 1 µ 0 B(r, t) (2.33)
• Linear, isotropic, homogeneous and time invariant medium with respect to electrical properties,

D(r, t) = " 0 " r (t) ⇤ Ẽ(r, t) = " 0 Z t 0 <t " r (t t 0 ) Ẽ(r, t 0 )dt 0 (2.34)
where "(t) is the impulse dielectric permittivity.

Boundary conditions

As is evident from the Maxwell's equation, in general the fields Ẽ, B, D and H are discontinuous at points where ", µ and also are. Hence the field vectors will be discontinuous at a boundary between two media with di↵erent constitutive parameters.

The integral form of Maxwell's equations can be used to determine the relations, called boundary conditions, of the normal and tangential components of the fields at the interface between two regions with di↵erent constitutive parameters ", µ and where surface density of sources may exist along the boundary.

The boundary condition for D can be calculated using a very thin, small pill-box that crosses the interface of the two media. Applying the divergence theorem [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF] to (2.1) and we have:

I D • ds = Z Base1 D1 • ds + Z Curved surf ace D • ds + Z Base2 D2 • ds = Z ⇢dv (2.35)
where D1 denotes the value of D in medium 1, and D2 the value in medium 2. Since both bases of the pillbox can be made as small as we like, the total outward flux of D over them is (D n1 D n2 )ds = ( D1 D2 ) • nds, where these D n are the normal drawn from medium 2 to medium 1 and n is the unit normal vector. At the limit, by taking a shallow enough pillbox, we can disregard the flux over the curved surface, whereupon the sources of D reduce to the density of surface free charge ⇢ s on the interface,

n • ( D1 D2 ) = ⇢ s (2.36)
Hence the normal component of D changes discontinuously across the interface by an amount equal to the free charge surface density ⇢ s on surface boundary.

Similarly the boundary condition for B can be established using the Gauss' law for magnetic fields. Since the magnetic field is solenoidal, it follows that the normal components of B are continuous across the interface between two media, n • ( B1 B2 ) = 0 (2.37)

The behavior of the tangential components of Ẽ can be determined using an infinitesimal rectangular loop at the interface which has sides of length dh, normal to the interface, and sides of length dl parallel to it. From the integral form of the Faraday's law and defining t as the unit tangent vector parallel to the direction of integration on the upper side of the loop, we have:

( Ẽ1 • t Ẽ2 • t)dl + contributions of sides dh = @ B @t • ds (2.38)
In the limit, as dh ! 0, the area ds = dhdl bounded by the loop approaches zero and, since B is finite, the flux of B vanishes. Hence ( Ẽ1 Ẽ2 ) • t = 0 and we conclude that the tangential components of Ẽ are continuous across the interface between two media.

In term of the normal n to the boundary, this can be written as:

n ⇥ ( Ẽ1 Ẽ2 ) = 0 (2.39)
where the symbol ⇥ designate the cross product.

Analogously, using the same infinitesimal rectangular loop, it can be deduced from the generalized Ampere's law that

( H1 • t H2 • t)dl + contributions of sides dh = ( @ D @t + J) • ds (2.40)
where, since D is finite, its flux vanishes. Nevertheless, the flux of the surface current can have a non-zero value when the integration loop is reduced to zero, if the conductivity of the medium 2, and consequently Js , is finite. This requires the surface to be a perfect conductor. Thus, n ⇥ ( H1 H2 ) = Js (2.41) the tangential component of H is discontinuous by the amount of surface current density Js . For finite conductivity, the tangential magnetial field is continuous across the boundary.

A summary of the boundary conditions are given for the general case and for the case when the medium 2 is a perfect conductor:

General boundary conditions n ⇥ ( Ẽ1 Ẽ2 ) = 0 (2.42) n ⇥ ( H1 H2 ) = Js (2.43) n • ( D1 D2 ) = ⇢ s (2.44) n • ( B1 B2 ) = 0 (2.45)
where Js and ⇢ s are potential surface density of charge or current.

Boundary conditions when the medium 2 is a perfect conductor ( 2 ! 1)

n ⇥ Ẽ1 = 0 (2.46) n ⇥ H1 = Js (2.47) n • D1 = ⇢ s (2.48) n • B1 = 0 (2.49)

The conservation of energy

Poynting's theorem represents the electromagnetic energy-conservation law. To derive the theorem, let us calculate the divergence of the vector field Ẽ ⇥ H in a homogeneous, linear and isotropic finite region V bounded by a closed surface S. If we assume that V contains power sources generating currents J, then, from Maxwell's equations, we get:

r • ( Ẽ ⇥ H) = H • r ⇥ Ẽ Ẽ • r ⇥ H = H • @ Ẽ @t Ẽ • @ D @t Ẽ • ( Ẽ + J) (2.50)
where J represents the source current density distribution which is the primary origin of the electromagenetic fields, while the induced conduction current density is written as Jc = Ẽ.

As the medium is assumed to be linear and no dispersive, the derivatives with respect to time can be written as

Ẽ • @ D @t = " Ẽ • @ Ẽ @t = @ @t ( "E 2 2 ) = @ @t ( Ẽ • D 2 ) (2.51) H • @ B @t = µ H • @ H @t = @ @t ( µH 2 2 ) = @ @t ( B • H 2 ) (2.52)
By introducing the equations (2.51) and (2.52) into (2.50), integrating over the volume V , applying the divergence theorem, and then rearranging terms, we have

Z V J • Ẽdv = @ @t Z V 1 2 ( Ẽ • D + B • H)dv Z V E 2 dv I S ( Ẽ ⇥ H) • ds (2.53)
To interpret this result we accept that

U ev = Ẽ • D 2 (2.54)
and

U mv = B • H 2 (2.55)
represent, as a generalization of their expression for static fields, the instantaneous electric energy density, U ev , and magnetic energy density, U mv , stored in the respective fields. Let us recall the empirical Lorenz force equation, which gives the electromagnetic force density, f (in N • m 3 ), acting on a volume charge density ⇢ moving at a velocity u (in m • s 1 ) in a region where an electromagnetic filed exists,

f = ⇢( Ẽ + ũ ⇥ B) = ⇢ Ẽ + J ⇥ B (2.56)
where J = ⇢ũ is the current density in terms of the mean drift velocity of the particles, which is independent of any random velocity due to collisions. The total force F exerted on a volume of charge is calculated by integrating f in this volume. For a single particle with charge q the Lorentz force is:

F = q( Ẽ + ũ ⇥ B) (2.57)
The work done by the electromagnetic field that acting on a volume density ⇢ inside a volume dv during a time interval dt is

dW = f • ũdtdv = ⇢( Ẽ + ũ ⇥ B) • ũdtdv = ⇢ Ẽ • ũdtdv = Ẽ • Jdtdv (2.58)
This work is transformed into heat. The corresponding power density

P v (in W • m 3 )
that the electromagnetic field supplies to the charge distribution is:

P v = dP dv = dW dtdv = Ẽ ⇥ J (2.

Time-harmonic electromagnetic fields

A particular case of great interest is on in which the sources vary sinusoidally in time.

In linear media, the time-harmonic dependence of the source gives rise to fields which, once having reached the steady state, also vary sinusoidally in time. However, timeharmonic analysis is important not only because many electromagnetic systems operate with signals that are practically harmonic, but also because arbitrary periodic time functions can be expanded into Fourier series of harmonic sinusoidal components while transient nonperiodic functions can be expressed as Fourier integrals. Thus, since the Maxwell's equations are linear di↵erential equations, the total fields can be synthesized from its Fourier components.

Analytically, the time-harmonic variation is expressed using the complex exponential notation based on Euler's formula, where it is understood that the physical fields are obtained by taking the real part, whereas their imaginary part is discarded. For example, an electric field with time-harmonic dependence given by cos(!t + '), where ! is the angular frequency, is expressed as

Ẽ = Re( Ẽe j!t ) = 1 2 ( Ẽe j!t + ( Ẽe j!t ) ⇤ ) = Ẽ0 cos(!t + ') (2.62)
where Ẽ is the complex phasor,

Ẽ = Ẽ0 e j' (2.63)
of amplitude E 0 and phase ', which will in general be a function of angular frequency and coordinates. The asterisk ⇤ indicates the complex conjugate, and Re() represents the real part of what is in the brackets.

Assuming e j!t time dependence, we can get the phasor form or time-harmonic form of Maxwell's equations simply by changing the operator @ @t to the factor j! and eliminating the factor e j!t . Maxwell's equations in di↵erential and integral forms for time-harmonic fields are given below.

Di↵erential form of Maxwell's equation for time-harmonic fields

r • D = ⇢ (Gauss 0 law) (2.64)
r • B = 0 (Gauss 0 law f or magnetic f ields) (2.65)

r ⇥ Ẽ = j! B (F araday 0 s law) (2.66)
r ⇥ H = J + j! D (Generalized Ampere 0 s law) (2.67)

Integral form of Maxwell's equation for time-harmonic fields

I S D • ds = Q T (Gauss 0 law) (2.68) I S B • ds = 0 (Gauss 0 law f or magnetic f ields) (2.69) I Ẽ • d l = j! Z S B • ds (F araday 0 s law) (2.70) I H • d l = Z S ( J + j! D) • ds (Generalized Ampere 0 s law) (2.71)
For the linear, homogeneous and invariant electrical medium, the constitutive relations become:

H(r, f ) = 1 µ 0 B(r, f ) (2.72) D(r, f ) = " c Ẽ(r, f ) (2.73)
with

" c = " 0 "r (f ), "r (f ) = F T (" r ) (2.74)
where "(f ) depends on the frequency for a dispersive medium. We can write "r = " 0 r +j" 00 r where " 00 r < 0 for f > 0. "r is the relative complex permittivity.

Similar process occurs in magnetic and conducting media, and, within a given frequency range, there may be a phase shift between Ẽ and Jc or between B and H which, at the macroscopic level, is reflected in the corresponding complex constitutive parameters c = 0 + j 00 and µ c = µ 0 + jµ 00 .

For a medium with complex permittivity, the complex phasor form of the displacement current is:

j! D = j!" c Ẽ = !" 00 Ẽ + j!" 0 Ẽ (2.75)
with " 0 = " 0 " 0 c and " 00 = " 0 " 00 c . While the sum, of the displacement and conduction current, called total induced current, Ji , is

Ji = Ẽ + j!" c Ẽ = ( + !" 00 ) Ẽ + j!" 0 Ẽ = Jd + Jr (2.76)
where Jd , called the dissipative current,

Jd = ( + !" 00 ) Ẽ (2.77)
in phase with the electric field, is the real part of the induced current Ji while Jr , called the reactive current,

Jr = j!" 0 Ẽ (2.78)
is the imaginary part of the induced current which is in phase quadrature with the electric field. The dissipative current can be expressed in a more compact form as

Jd = e Ẽ (2.79)
where e is the e↵ective or equivalent conductivity e = + !" 00 (2.80) which includes the ohmic losses due to and the damping losses due to !" 00 . Thus the induced current, can be written as

Ji = e Ẽ + j!" 0 Ẽ = ec Ẽ (2.81)
where ec is the complex e↵ective conductivity, defined as

ec = e + j!" 0 (2.82)
Thus a medium with conductivity ec and null permittivity is formally equivalent to one with conductivity and permittivity, and " c , respectively.

For harmonic signals the boundary conditions of the normal and tangential components of the fields at the interface between two regions with di↵erent constitutive parameters ", µ and , become:

General boundary conditions n ⇥ ( Ẽ1 Ẽ2 ) = 0 (2.83) n ⇥ ( H1 H2 ) = Js (2.84) n • ( D1 D2 ) = ⇢ s (2.85) n • ( B1 B2 ) = 0 (2.86)
Boundary conditions when the medium 2 is a perfect conductor ( 2 ! 1)

n ⇥ Ẽ1 = 0 (2.87) n ⇥ H1 = Js (2.88) n • D1 = ⇢ s (2.89) n • B1 = 0 (2.90)
In formulating the conservation energy equation for time-harmonic fields, it is convenient to find, first, the time-average Poynting vector over a period, i.e. the time-average power passing through a unit area perpendicular to the directon of P. We have:

Ẽ = Re( Ẽe j!t ) = 1 2 ( Ẽe j!t + ( Ẽe j!t ) ⇤ ) (2.91) H = Re( He j!t ) = 1 2 ( He j!t + ( He j!t ) ⇤ ) (2.92)
Thus, the instantaneous Poynting vector can be written as:

P = Ẽ ⇥ H = Re( Ẽe j!t ) ⇥ Re( Hj!t ) = 1 2 Re( Ẽ ⇥ H⇤ + Ẽ ⇥ He 2j!t ) (2.93)
The time-average value of the instantaneous Poynting vector can be calculated integrating the above equation over period, i.e.,

Pav = 1 T Z T 0 Pdt = 1 2T Z T 0 Re( Ẽ⇥ H⇤ + Ẽ⇥ He 2j!t )dt = 1 2 Re( Ẽ⇥ H⇤ ) = Re( Pc ) (2.94)
since the time average of Ẽ ⇥ He 2j!t vanishes. The magnitude,

Pc = 1 2 Ẽ ⇥ H⇤ (2.95)
is termed as the complex Poynting vector. Thus the time-average of the Poynting vector is equal to half of the real part of the complex Poynting vector.

Plane wave and propagation equations

In this thesis, we work with the time-harmonic electromagnetic fields. The time dependence of the plane wave is e j!t , it will be omitted in the calculus. To exhibit the propagation equations, we apply the following mathematical formula to Maxwell's equations [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF]:

r ⇥ (r ⇥ Ṽ ) = r(r • Ṽ ) r 2 Ṽ (2.96)
Then the propagation equations of Ẽ and H can be written as:

r 2 Ẽ + " c µ 0 ! 2 Ẽ = 1 " c r⇢ + j!µ 0 J (2.97) r 2 H + " c µ 0 ! 2 H = r ⇥ J (2.98)
If the medium contains neither free-charge nor free-current (⇢ = 0, J = 0), then the propagation equations become:

r 2 Ẽ + " c µ 0 ! 2 Ẽ = 0 (2.99) r 2 H + " c µ 0 ! 2 H = 0 (2.100)
These are the Helmholtz equations [START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Stratton | Electromagnetic Theory. International series in physics[END_REF][START_REF] Ishimaru | Electromagnetic wave propagation, radiation and scattering[END_REF]. They have a particular solution in the following form:

Ẽ(r, t) = Ẽ0 e j k•r (2.101) H(r, t) = H0 e j k•r (2.102)
where Ẽ0 and H0 are independent of r.

The vector k is the wave vector of propagation medium. We have the following equation:

k2 = " 0 "r µ 0 ! 2 (2.103)
If the medium is transparent, then "r and k are real. We have that:

k = 2⇡ (2.104)
where is the wavelength. From the two equations above, we have:

" 0 "r µ 0 v 2 = 1 (2.105)
where v is a constant velocity characterizing the propagation medium. In particular, if the medium is the vacuum, then we have:

" 0 "r c 2 = 1 (2.106)
where c is the velocity of light. For a lossless medium with optical index ⌫ = p "r , the constant v is:

v = c ⌫ (2.107)
Equation (2.101) and (2.102) express a monochromatic plane wave with the propagation direction given by the wave vector k. For the plane wave described by k • r = Constant, the amplitudes of the components of the wave are independent of the position and remain constant.

From Maxwell's equations, we know that:

Ẽ = H ⇥ k " 0 "r ! (2.108)
With the definition of impedance of medium: Z = q µ 0 " 0 "r , we have:

Ẽ = Z H ⇥ k k (2.109)
For the wave propagation, if we consider the Cartesian coordinate (Oxyz), with orthogonal basis (ũ x , ũy , ũz ). An incident monochromatic plane wave propagates in the space constituted of two media that are separated by an interface. The incident wave vector k0 is located in the plane (xOz). The direction of propagation of the incident wave is represented by the angle ✓ 0 with respect to the Oz axis. The polarization of a plane wave is then determined based on the curve that is going to describe the electric field Ẽ in a wave plane. This polarization is in general elliptical and can be decomposed into a combination of two linear polarizations: horizontal and vertical.

Horizontal polarization corresponds to the case where the electric field Ẽ is perpendicular to the plane of incidence formed by the couple of vectors ( k0 , ũz ) where ũz is the unit vector of Oz axis. It is also called the transverse electric polarization (denoted TE, h or S). We shall call polarization E parallel and will be denoted by E // because the Ẽ field is parallel to the plane (xOy). The situation is similar with vertical polarization where

H replaces Ẽ.
The following table gives convention of notation, the notations in the same column represent the same polarization.

Horizontal Polarization Vertical Polarization

h v

E // H //
TE TM S P

The wave vector k0 can be represented by the incidence zenith angle ✓ 0 and the azimuth angle ' 0 :

k0 = ↵ 0 ũx + 0 ũy 0 ũz (2.110) with 8 > > > > > > > > > < > > > > > > > > > : ↵ 0 = k 0 sin ✓ 0 cos ' 0 0 = k 0 cos ✓ 0 0 = k 0 sin ✓ 0 sin ' 0 k 0 = 2⇡ (2.111)
In fact, all the solutions of the wave propagation problem can be expressed as a combination of elementary plane waves with di↵erent amplitudes and wave vectors. We will discuss this in details in the next chapter.

Conclusion

In this chapter, we present the electromagnetic field theory fundamentals. Especially, we present the Maxwell's equations and the interaction of an electromagnetic field with an object.

Despite their apparent simplicity, Maxwell's equations are in general not easy to solve.

In fact, even in the most favorable situation of homogeneous, linear and isotropic media, there are not many problems of interest that can be analytically solved except for those presenting a high degree of geometrical symmetry. Moreover, the frequency range of scientific and technological interest can vary by many orders of magnitude, expanding from frequency value of zero (or very low) to roughly 10 with both negative permittivity (" < 0)) and negative permeability (µ < 0). These media are called double-negative metamaterials and, owing to their unusual electromagnetic properties, they present many potential technological applications.

Another important factor to study the interaction of an electromagnetic field with an object is the electric size of the body, i.e., the relationship between the wavelength and the body size, which can also vary by several orders of magnitude. All these circumstances make it in general necessary to use analytical, semi-analytical or numerical methods appropriate to each situation. In particular, numerical methods are fundamental for simulating and solving complex problems that do not admit analytical solutions.

In the next chapter, we will present the scattering by gratings and by random rough surfaces.

Chapter 3

Scattering by gratings and by random rough surfaces

In this thesis, the main aim is to study the di↵raction by gratings and the scattering by random rough surfaces illuminated by a electromagnetic plane wave. In this chapter, I will introduce the fundamental theory about these two aspects.

The theory of di↵raction gratings

Introduction

Di↵raction gratings are optical components used to separate light into its component wavelengths. Di↵raction gratings are used in spectroscopy, or for integration into spectrophotometers or monochromators. Di↵raction gratings consist of a series of closely packed grooves that have been engraved or etched into the grating surface. Di↵raction gratings can be either transmissive or reflective. As light transmits through or reflects o↵ a grating, the grooves cause the light to di↵ract, dispersing the light into its component wavelengths [START_REF] Palmer | Di↵raction grating handbook[END_REF].

For practical applications, gratings generally have ridges or rulings on their surface rather than dark lines. Such gratings can be either transmissive or reflective. Gratings which modulate the phase rather than the amplitude of the incident light are also produced, frequently using holography [START_REF] Palmer | Di↵raction grating handbook[END_REF].

The grating equation

When monochromatic light is incident on a grating surface (i.e. a periodic surface), it is di↵racted into discrete directions. We can picture each grating groove as being a very small, slit-shaped source of di↵racted light. The light di↵racted by each groove combines to form set of di↵racted wavefronts. The usefulness of grating depends on the fact that there exists a unique set of discrete angles along which, for a given spacing D between grooves, the di↵racted light from each facet is in phase with the light di↵racted from any other facet, leading to constructive interface [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF].

Di↵raction by a grating can be visualized from the geometry in the figure 3 The geometry path di↵erence between light from adjacent grooves is seen to be D sin ✓ 0 + D sin ✓ n . The principle of constructive interface dictates that only when this di↵erence equals the wavelength of the light, or some integral multiple thereof, will the light from adjacent grooves be in the phase (leading to constructive interface). All other angles the wavelets originating from the groove facets will interface destructively.

These relationships are expressed by the equation:

1 < sin ✓ n = sin ✓ 0 + n D < 1 (3.1)
which governs the angular locations of the principal intensity maxima when light of wavelength is di↵racted from a gratings of groove spacing D. Here n is the di↵raction order (or spectral order), which is an integer. For a wavelength , all values of n for which | sin ✓ 0 +n /D| < 1 correspond to propagating (rather than evanescent) di↵raction orders. The special case n = 0 leads to the law of reflection

✓ n = ✓ 0 .
It is sometimes convenient to write the grating equation as

Gn = sin ✓ 0 + sin ✓ n (3.2)
where G = 1/D is the groove frequency or groove density, more commonly called "groove per millimeter".

Equation 3.1 and 3.2 are the common forms of the grating equation, but their validity is restricted to cases in which the incident and di↵racted rays lie in a plane which is perpendicular to the grooves (at the center of grating). The majority of grating systems fall within this category, which is called classical di↵raction. If the incident light beam is not perpendicular to the grooves, the grating equation must be modified:

Gn = cos "( sin ✓ 0 + sin✓ n ) (3.3)
Here, " is the angle between the incident light path and the plane perpendicular to the groove at the grating center. In geometries, for which " 6 = 0, the di↵racted spectra lie on a cone rather than in a plane, so such cases are termed conical di↵raction.

For a grating of groove spacing D, there is a purely mathematical relationship between the wavelength and the angles of incidence and di↵raction. In a given spectral order n, the di↵erent wavelength of polychromatic wavefronts incident at angle ✓ 0 are separated in angle

✓ n ( ) = sin 1 ( n D + sin ✓ 0 ) (3.4)
When n = 0, the grating acts as a mirror, and the wavelength are not separated (✓ n = ✓ 0 for all ), this is called specular reflection or simply the zeroth order.

A special but common case is that in which the light is di↵racted back towards the direction from which it came (i.e. ✓ 0 = ✓ n ), this is called the Littrow configuration, for which the grating equation becomes:

n = 2d sin ✓ 0 (3.5)

Di↵raction orders

Generally several integers n will satisfy the grating equation. We call each of these values a di↵raction order [START_REF] Palmer | Di↵raction grating handbook[END_REF].

For a particular groove spacing D, wavelength and incidence angle ✓ 0 , the grating equation 3.1 is generally satisfied by more than one di↵raction angle ✓ n . In fact, subject to restrictions discussed below, there will be several discrete angles at which the condition for constructive interference is satisfied. The physical significance of this is that the constructive reinforcement of wavelets di↵racted by successive grooves merely requires that each ray be retarded (or advanced) in phase with every other, this phase di↵erence must therefore correspond to a real distance (path di↵erence) which equals an integral multiple of the wavelength. This happens, for example, when the path di↵erence is one wavelength, in which case, we speak of the positive first di↵raction order (n = 1) or the negative first di↵raction order (n = 1), depending on whether the rays are advanced or retarded as we move from groove to groove.

The grating equation reveals that only those spectral orders for which | sin ✓ 0 +n /D| < 1 can exist. This restriction prevents light of wavelength from being di↵racted in more than finite number of orders. Specular reflection (n = 0) is always possible. In most cases, the grating equation allows light of wavelength to be di↵racted into both negative and positive orders as well. Explicitly, spectra of all orders n exist for which,

( 1 sin ✓ 0 )D < n < (1 sin ✓ 0 )D, n is an integer (3.6)
For /D ⌧ 1, a large number of di↵racted orders will exist.

The most troublesome aspect of multiple order behavior is that successive spectral overlap. It is evident from the grating equation that light of wavelength di↵racted by a grating along direction ✓ n will be accompanied by integral fraction /2, /3, etc. That is for any grating instrument configuration, the light of wavelength di↵racted in the n = 1 order will coincide with the light of wavelength /2 in the n = 2 order, etc. This superposition of wavelengths, which would lead to ambiguous spectroscopic data, is inherent in the grating equation itself and must be prevented by suitable filtering (called order sorting), since the detector cannot generally distinguish between light of di↵erent wavelengths incident on it.

Rayleigh expansion

Suppose the interface of the two media is described by the function z = a(x, y), outside the deformation, the di↵racted field ( Ẽ, H) could be represented by the so-called Rayleigh expansion. For example, in the medium 1, when z > max a(x, y), the electric field Ẽ and the magnetic field H in E // could be represented with the help of the particular solutions as in equation (2.101) and (2.102). These are called the Rayleigh expansions [START_REF] Petit | Electromagnetic Theory of Gratings[END_REF]. In the discrete case, they can be represented as a linear combination of elementary plane waves. In the continuous case, they can be represented as a integral of the elementary plane waves. The Rayleigh expansion is only valid outside the modulated zone (i.e. z > max a(x, y) or z < min a(x, y) ) [START_REF] Petit | Sur la di↵raction d'une onde plane par un rseau infiniment conducteur[END_REF][START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface[END_REF][START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Van Den Berg | The Rayleigh hypothesis in the theory of reflection by a grating[END_REF][START_REF] Keller | Singularities and Rayleigh's hypothesis for di↵raction gratings[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF].

For example, if we consider only the one-dimensional interface for simplicity. In the Cartesian referential Oxyz, the grating is represented by a periodic cylindrical surface y = a(x). This surface separates the air (medium 1) from the medium with a real or complex refractive index (medium 2). The grating of period D is illuminated by a monochromatic plane wave under the incidence ✓ 0 . The incident wave vector lies in the xOy plane. The letter m denotes indi↵erently the upper medium (m = 1) or the lower medium (m = 2). Henceforth, n (m) , Z (m) and k (m) indicate the optical index, the impedence and the wave number of medium m. In this case, when y > max(a(x)) and y < min(a(x)), the di↵racted field can be represented by a combination of elementary plane waves, the Rayleigh expansion:

8 > < > : F (m) c (x, y) = P n (c (m+) n exp( j↵ n x)exp( j (m) n y) + c (m ) n exp( j↵ n x)exp(j (m) n y)) G (m) c (x, y) = P n (m) n k (m) (c (m+) n exp( j↵ n x)exp( j (m) n y) + c (m ) n exp( j↵ n x)exp(j (m) n y)) (3.7)
The subscript (c) denotes the Cartesian components of electromagnetic field. In

E // polarization, F (m) c (x, y) = E (m) z (x, y), G (m) c = Z (m) H (m)
z (x, y) and in H // polarization,

F (m) c (x, y) = Z (m) H (m) z (x, y), G (m) c (x, y) = E (m)
x (x, y). Superscripts (+) and ( ) denotes a plane wave moving in direction along the y-axis and inverse the y-axis, respectively. The propagation coe cients of the n-th order di↵raction are presented by

↵ n and (m) n
with the relation

↵ 2 n + ( (m) n ) 2 = k (m)2 (3.8)
where Im( is imaginary or complex.

(m) n ) < 0 and ↵ n = k (1) sin✓ 0 + n 2⇡ D . (↵ n , (m)

The scattering matrix S

With the discrete version of Rayleigh expansion, we can define the scattering matrix(Smatrix) to relate the amplitudes of outgoing plane waves to those of incoming waves.

Take the one-dimensional case for simplicity, we have:

0 @ c (1+) c (2 ) 1 A = S 0 @ c (1 ) c (2+) 1 A (3.9)
here we use c (m±) to represent a vector containing the scattering amplitudes c

(m±) n

. For a perfectly conduction surface, the scattering matrix is given by:

c (1+) = Sc (1 ) (3.10)
Figure (3.2) illustrates the link between incoming and outcoming plane waves. n /k (m) . According to the definition of the complex Poynting vector (2.95), for an incidence angle ✓ 0 , the e ciency ✏ (m) n is given by [START_REF] Petit | Electromagnetic Theory of Gratings[END_REF]:

✏ (m) n = ⌫ (m) cos ✓ (m) n cos ✓ 0 |c (m) n | 2 (3.11)
For a lossless grating, the sum of e ciencies is equal to 1 according with the conservation of power [START_REF] Petit | Electromagnetic Theory of Gratings[END_REF]. For the two media case, considering the propagation plane waves, we have:

X n ✏ (1) n + X n ✏ (2) n = 1 (3.12)
3.2 Scattering from random rough surfaces

Introduction

The problem of electromagnetic scattering from random rough surfaces has aroused the interest of physicists and engineers for many years because of its wide range of applications in optics, material science, communications, oceanography and remote sensing.

The three classical analytical methods commonly used in random rough surfaces scattering are the small-perturbation method [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF], the Kirchho↵ method [START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF][START_REF] Bourlier | Azimuthal harmonic coe cients of the microwave backscattering from a non-Gaussian ocean surface with the first-order SSA model[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] and the small slope approximation [START_REF] Voronovich | Wave scattering from rough surfaces[END_REF].

The electromagnetic analysis of rough surfaces with parameters close to the incident wavelength requires a rigorous formalism. Numerical methods based on Monte Carlo simulations are available for the study of electromagnetic wave scattering from onedimensional and two-dimensional random rough surfaces. In the frequency domain, the boundary integral method can be used to analyze the scattering problem by rough surfaces. In this case, the electric or magnetic field integral equation is converted into matrix equations using the method of moments (MoM) [START_REF] Warnick | Numerical simulation methods for rough surface scattering[END_REF][START_REF] Saillard | Rigorous solutions for electromagnetic scattering from rough surfaces[END_REF]. The average mesh length determines the accuracy of the MoM solution. The number of unknowns N is proportional to surface area in square wavelength. Several fast methods have been proposed to reduce the CPU time [START_REF] Pak | Numerical simulations and backscattering enhancement of electromagentic waves from two-dimensional dielectric random rough surfaces with the sparse-matrix canonical grid method[END_REF][START_REF] Wagner | Monte Carlo simulations of electromagentic scattering from two-dimensional random rough surfaces[END_REF][START_REF] Jandhyala | Fast algorithm for the analysis of scattering by dielectric rough surfaces[END_REF][START_REF] Torrungrueng | Numerical studies of backscattering enhancement of electromagnetic waves from two-dimensional random rough surfaces with the forward-backward novel spectral acceleration method[END_REF][START_REF] Soriano | Scattering of electromagentic waves from twodimensional rough surfaces with an impedance approximation[END_REF][START_REF] Tsang | Wave scattering with the UV multilevel partitioning method: 2. three-dimensional problem of nonpenetrable surface scattering[END_REF] and lead to a computational e ciency of O(N log N ) . These methods are fully capable of describing the field scattered from surfaces of large size. Other numerical approaches are also suggested and for topical reviews, see [START_REF] Warnick | Numerical simulation methods for rough surface scattering[END_REF] and [START_REF] Saillard | Rigorous solutions for electromagnetic scattering from rough surfaces[END_REF]. Exact methods require solutions for many realizations of twodimensional rough surfaces. The Monte-Carlo is used and the average scattered power is estimated over results of several surface realizations.

Random rough surface generation

In order to study the di↵raction phenomenon of random rough surfaces, we use a numerical solution method based on exact electromagnetic model. This method requires the inputs describing surfaces to be analyzed. These inputs are numerical representation of the actual surfaces.

If we take the method of small perturbations (SPM) for example [START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF], which is an approximate analytical method, we only need the geometric parameters (statistics) which characterize the surface to be analyzed (standard deviation of heights, length correlation ...) to study its average response to electromagnetic excitation. In our case, these parameters will not intervene directly in the electromagnetic treatment. They will be used during a preliminary step in the numeric generation of surfaces that we want to analyze. These numeric profiles will be the inputs of the electromagnetic model that we have developed.

In a space based on the Cartesian orthonormal (Oxyz), a surface whose generating line is based on a function y = a(x) is called as cylindrical surface or one-dimensional surface.

This surface is invariant along the z direction. For the two-dimensional surface, we use the equation z = a(x, y) to represent the surface function.

In this thesis, we consider the one-dimensional or two-dimensional bounded supported random surfaces. They are expected from random process, verifying some assumptions that we explain later in this chapter. The randomness of these surfaces requires a statistical study to characterize it. To better understand the interaction between electromagnetic waves and the rough surface will require a good description of them.

In signal theory, a random process represents the evolution of a random variable with time or space. It is symbolized by a random function depending on time and/or on space and on a parameter W reflecting the randomness. A random process depending only on the space (time-independent) is noted as ⇠(r, W). For a given value W 0 of W, we obtain a realization of the random process ⇠(r, W 0 ). This realisation is deterministic.

If we vary the random parameter W, we get a set of realisations from the same random process.

The statistical description of a random spatial process, whether in the one-dimensional or two-dimensional case, can be done by studying its spatial fluctuations depending on the position. For a given position random functions ⇠(x 0 , W) for the one-dimensional case, and ⇠(x 0 , y 0 , W) for the two-dimensional case, are random variables.

The distribution function of the random variable ⇠(r 0 , W) is defined as:

F ⇠ (h) = P rob(⇠  h) (3.13)
For a real random variable, the derivative of its distribution function (if it exists) gives the density of probability P ⇠ (h). The density of probability verifies the following properties:

P rob(⇠ 2 [a, b]) = Z b a P ⇠ (h)dh (3.14) 8h, P ⇠ (h) > 0, Z +1 1 P ⇠ (h)dh = 1 (3.15)
In practice, this function can be estimated by the normalized histogram of the values taken by the random variable ⇠(x, W) or (x, y, W) for a given position. In the Gaussian case, the probability density is determined entirely by two parameters which are the statistical mean of the random variable which we denote as m ⇠ and the standard deviation ⇠ . ⇠ measures the dispersion of the values of the random variable around the mean value m ⇠ .

m ⇠ = E[⇠] = Z +1 1 hP ⇠ (h)dh (3.16) 2 ⇠ = E[⇠ 2 ] E[⇠] 2 (3.17) E[⇠ 2 ] = Z +1 1 h 2 P ⇠ (h)dh (3.18)
where E denotes the expectation which provides the statistical mean of the considered random variable. The analytical expression of a Gaussian probability density is given by:

P ⇠ (h) = 1 ⇠ p 2⇡ e 1 2 ( h m ⇠ ⇠ ) 2 (3.19)
Now we consider a two-dimensional Gaussian random variable (⇠ 1 , ⇠ 2 ), with statistical average (m ⇠ 1 , m ⇠ 2 ) and standard deviation ( ⇠ 1 , ⇠ 2 ). The associated joint probability density is:

P ⇠ 1 ,⇠ 2 (h 1 , h 2 ) = 1 2⇡ ⇠ 1 ⇠ 2 p 1 ⇢ 2 c exp( (h 1 m ⇠ 1 ) 2 2 ⇠ 1 2⇢c(h 1 m ⇠ 1 )(h 2 m ⇠ 2 ) ⇠ 1 ⇠ 2 + (h 2 m ⇠ 2 ) 2 2 ⇠ 2 2(1 ⇢ 2 c ) ) (3.20)
with ⇢ c the correlation coe cient:

⇢ c = E[⇠ 1 ⇠ 2 ] m ⇠ 1 m ⇠ 2 ⇠ 1 ⇠ 2 (3.21)
To analyze a random process ⇠(r, W), we can consider a study of a single realisation (W fixed). This allows us to know its spatial moments. In the general case, these moments may depend on the realisation, that is to say the random W. The other way to do is to look at the values of the random variable associated with the process for a given position r0 and several realizations of the process ⇠(r, W) (family of realisations). In this case, the process is described by using these statistical moments (statistical average and higher order moments).

Knowing the probability density of the random variable ⇠(r, W) for a given r0 , the statistical moment of order n, associated to this random variable is the expectation of order n. It is defined for a continuous random variable by:

m ⇠ n = E[⇠ n ( r0 , W)] = Z +1 1 h n P ⇠ (h)dh (3.22)
To evaluate the correlation that may exist between two values taken by the random variable ⇠(r, W ) , at two di↵erent points r0 and r0 + r, we compute its statistical autocorrelation function defined by:

R( r0 , r0 + r) = Z +1 1 Z +1 1 hh 0 P (⇠ r0 ,⇠ r0 +r ) (h, h 0 )dhdh 0 (3.23) P ⇠ r0 ,⇠ r0 +r
(h, h 0 ) is the joint probability. This function depends on the position of two points, that is to say, ( r0 , r + r0 ).

If we work now on a single realisation, we have access to spatial moments. If we denote the extent of ⇠(r, W 0 ), then the spatial average of ⇠(r, W 0 ) can be estimated by:

⇠(r, W) = lim !+1 1 Z 2 2 ⇠(r, W 0 )dr (3.24)
The spatial autocorrelation is defined as:

C ⇠,⇠ (r, W 0 ) = ⇠( r0 , W 0 )⇠( r0 + r, W 0 ) = lim !+1 1 Z 2 2 ⇠( r0 , W 0 )⇠( r0 + r, W 0 )d r0 (3.25)
Without prior assumptions about the random processes, the spatial moments depend on the realization. We will see what are the assumptions that allow the equality between spatial moments and statistical moments up to a certain order n. In general , the mean and spatial autocorrelation depend on the realisation. This means that these are random variables. If the spatial moments are independent of the random hazard until the order 2, then the spatial process is ergodic to order 2. Concerning the statistical moments, the calculated statistical average depends, in general, on the position r and the autocorrelation function of the two positions r0 and r + r0 . If now the autocorrelation function depends only on the distance between the two points and that the statistical average is constant when we change position, then the process is said to be stationary to order 2. If these last two properties (ergodicity and stationarity to order 2) are made simultaneously satisfied, then Birko↵'s theorem allows to say that the statistical moments are equal to the spatial moments up to order 2. We can write:

m = E[⇠( r0 , W)] = ⇠(r, W 0 ) = Constant (3.26) R ⇠⇠ (r) = C ⇠⇠ (r) (3.27)
For the surface, the random variable a(r, W) defines the height of the surface at all the points r. The simulated surfaces are assumed to satisfy the two properties above.

Moreover, we assume that the surface satisfies the following conditions:

• The probability density of the surface height P (h) is Gaussian with mean E[a(r, W)] = 0 and standard deviation a . The density is given by:

P (h) = 1 a p 2⇡ exp( h 2 2 2 ) (3.28)
• The autocorrelation function is given by:

R(x, y) = 8 > < > : 2 a exp( x 2 l 2 x ), one dimensional case 2 a exp( x 2 l 2 x y 2 l 2 y
), two dimensional case

(3.29)
If the surface is isotropic, we have l x = l y = l.

We generate the random surface based on the principle of linear filtration of Gaussian white noise. In fact, the formula is as follows:

a(r) = (h f ⇤ B)(r) = Z +1 1 h f (r r0 )B( r0 )d r0 (3.30)
where B is the entry signal of the filter, it is related to a Gaussian white noise characterized by the Gaussian probability density and the autocorrelation function:

R BB (r) = 2 a (r) (3.31)
with (r) the Dirac distribution. From the equation (3.30), we have the formula for R aa (r):

R aa (r) = (C h f ⇤ R BB )(r) (3.32)
where C h f (r) is the spatial autocorrelation of the impulse response of filter:

C h f (r) = Z +1 1 h f (r)h f (r r0 )d r0 (3.33)
So we obtain that

R aa (r) = 2 a C h f (r) (3.34)
Now suppose that Ĥf (↵, ) is the Fourier transformation (FT) of the impulse response:

Ĥf (↵, ) = F T [h f (x, y)] (3.35)
Then we have

| Ĥf (↵, )| 2 = F T [C h f (x, y)] (3.36)
The formula of filtration gives:

F T [R aa (r)] = 2 a | Ĥf (↵, )| 2 (3.37)
Given the function R aa (r) and suppose that Ĥf (↵, ) = | Ĥf (↵, )|, we can calculate the impulse response:

h f (x, y) = F T 1 [ Ĥf (↵, )] = F T 1 [ p F T [R aa (x, y)]] (3.38)
In particular, for the two-dimensional isotropic Gaussian surface, the impulse response is

h f (x, y) = 2 l p ⇡ exp( 2( x 2 + y 2 l ) 2 ) (3.39)
To implement this method, we need to discrete version of the formula. Suppose x and y are the length of the step along the direction ũx and ũy . We have:

x = L x N x , y = L y N y (3.40)
and

a(x i , y j ) = x y X p X q h f (U p , V q )B(U p x i , V q y j ) (3.41)
with

x i = i x, y j = j y, U p = p x, V q = q y (3.42)

A = L x L y is the area of the generated surface. Thereafter, L = L x = L y and N x = N y = N e and N 2
e is the number of samples.

Beam of elementary plane waves

For the two dimensional interface, suppose the interface separating the air from a dielectric medium and described by the function z = a(x, y). For the E // polarization, the component E z will be zero and for the H // polarization, the component H z will be zero.

Nonzero components of Ẽ Nonzero components of H

E // E x , E y H x , H y , H z H // E x , E y , E z H x , H y
Outside the deformation, the di↵racted field ( Ẽ, H) could be represented by the Rayleigh integral when z > max a(x, y) or z < min a(x, y). The electric field Ẽ and the magnetic field H in E // could be represented as follows:

Ẽ(1) (x, y, z) = 1 4⇡ 2 Z +1 1 Z +1 1 C (1,E // ) (↵, ) h(↵, )e j k(1) (↵, )•r d↵d (3.43) Z (1) H(1) (x, y, z) = 1 4⇡ 2 Z +1 1 Z +1 1 C (1,E // ) (↵, )( k(1) (↵, ) k (1) ⇥ h(↵, ))e j k(1) (↵, )•r d↵d (3.44)
where h is the unit polarization vector:

h(↵, ) = p ↵ 2 + 2 ũx + ↵ p ↵ 2 + 2 ũy (3.45)
and

k(1) = ↵ũ x + ũy + (1) ũz (3.46) with ↵ 2 + 2 + ( (1) ) 2 = (k (1) ) 2 , Im( (1) )  0, Re( (1) ) 0. If ↵ 2 + 2 > (k (1)
) 2 , the constant of propagation is pure imaginary and corresponds to the evanescent wave. (1) is real and corresponds to the propagation wave. C (1,E // ) (↵, ) is the amplitude of the elementary wave e j k(1) (↵, )•r .

If ↵ 2 + 2  (k (1) ) 2 ,
For the H // , we have For the one-dimensional case, suppose the interface is described by y = a(x) which is invariant in the direction Oz. If the surface is illuminated by a plane wave with wave vector k0 which is in the plane Oxy. As in the two-dimensional case, we have the following table:

Z (1) H(1) (x, y, z) = 1 4⇡ 2 Z +1 1 Z +1 1 C (1,H // ) (↵, ) h(↵, )e j k(↵, )•r d↵d (3.47) Ẽ(1) (x, y, z) = 1 4⇡ 2 Z +1 1 Z +1 1 C (1,H // ) (↵, )( k(1) (↵, ) k (1) ⇥ h(↵, ))e j k(1) (↵, )•r d↵d (3.
Nonzero components of Ẽ Nonzero components of H E // E z H x , H y H // E x , E y H z
For E // , the electromagnetic fields are:

Ẽ(1) (x, y) = 1 2⇡ Z +1 1 C (1,E // ) (↵)e j k(1) (↵)•r d↵ũ z (3.49) Z (1) H(1) (x, y) = 1 2⇡ Z 1 1 C (1,E // ) (↵)( k(1) (↵) k (1) ⇥ ũz )e j k(1) (↵)•r d↵ (3.50) with k(1) (↵) = ↵ũ x + (1) ũy (3.51)
and

↵ 2 + ( (1) ) 2 = (k (1) ) 2 , Im( (1) )  0, Re( (1) ) 0 (3.52)
For the H // , the fields are:

Z (1) H(1) (x, y) = 1 2⇡ Z +1 1 C (1,H // ) (↵)e j k(1) (↵)•r d↵ũ z (3.53) Ẽ(1) (x, y) = 1 2⇡ Z 1 1 C (1,H // ) (↵)( k(1) (↵) k (1) ⇥ ũz )e j k(1) (↵)•r d↵ (3.54)

Scattering patterns

In the far-field zone, the Rayleigh expansion (3.43) and (3.44) is reduced to the only contribution of the propagation waves. For E // , the method of stationary phase [START_REF] Born | Principles of optics Electromagnetic theory of propagation, Interference and di↵raction of light[END_REF] leads to the asymptotic fields at the point M (r, ✓, '):

Ẽ(1) far (r, ✓, ') = C (1,E // ) (k (1) sin ✓ cos ', k (1) sin ✓ sin ') cos ✓ e jk (1) r r e j ⇡ 2 ũ' (3.55)
Z (1) H( 1) far (r, ✓, ') = C (1,E // ) (k (1) sin ✓ cos ', k (1) sin ✓ sin ') cos ✓ e jk {E // , H // }), the normalized bistatic scattering coe cient (ba) is defined as follows:

(ba) (✓, ') = 1

P (a) 0 dP (ba) d⌦ = |C (1,ba) (k (1) sin ✓ cos ', k (1) sin ✓ sin ') cos ✓| 2 2 L 2 cos ✓ 0 (3.57)
where dP (ba) d⌦ is the power scattered per unit solid angle d⌦ = sin ✓d✓d' with

dP (ba) = 1 2 Re( Ẽ(ba) far ⇥ H(ba)⇤ far dSũ r ) (3.58)
The symbol ⇤ designates the complex conjugate. dS is the element surface with dS = r 2 d⌦. The unit vectors ũr , ũ✓ , ũ' are drawn in the direction of increasing r, ✓ and ' such as to constitute a right-hand base system. P (a) 0

is the flux of incident power through the modulated region with L the modulated length (see equations 2.61 and 2.94 ):

P (a) 0 = 1 2 Z +L/2 L/2 Z L/2 L/2 Re( Ẽ(a) i ⇥ H(a)⇤ i dxdyũ z ) (3.59)
For random rough surface, the average bistatic scattering coe cient is defined: (1,ba) (k (1) sin ✓ cos ', k (1) sin ✓ sin ') cos ✓| 2 ] (3.60)

E[ (ba) (✓, ')] = 1 2 L 2 cos 2 ✓ cos ✓ 0 E[|C
For infinite extension surfaces, the three classical analytical methods (First-order perturbation method, first-order small slope approximation and the Kirchho↵ method) lead to closed-form formula for the average bistatic scattering coe cient [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF][START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF][START_REF] Voronovich | Wave scattering from rough surfaces[END_REF]. Exact methods require solutions for many realizations of two-dimensional rough surfaces. The Monte Carlo technique is applied to estimate the average bi-static coe cient from results over N R di↵erent realizations [START_REF] Warnick | Numerical simulation methods for rough surface scattering[END_REF].

E[ (ba) (✓, ')] = 1 N R N R X j=1 (ba) j (✓, ') (3.61)
Some authors prefer to use the radar cross section which is 4⇡ cos

✓ 0 L 2 E[ (ba) (✓, ')].

Conclusion

In this chapter, we present the theory of di↵raction gratings and the theory of scattering from random rough surfaces. We recall that the Rayleigh expansion is only valid outside the modulated zone. This is why we can not simply use Rayleigh expansion in the modulated zone and the C-method is needed here. For analyzing gratings, we present the concept of matrix S, see equation (3.9), this concept gives rise to chapter 8, the new version of C-method: C-method as an initial value problem. We also define the average bistatic coe cient that is a quantity measured in remote sensing and in optics.

In chapter 6, we compare the average bistatic coe cient estimated with the C-method and experimental data.

In the next chapter, we present the C-method and we show how this method leads to eigenvalue problem.

This formalism has been extensively used in the theory of grating: grating with finite conductivity in conical incidence, multilayer gratings with parallel or non parallel interface, bi-crossed gratings [START_REF] Popov | Conical di↵raction mounting. generalization of a rigourous di↵erential method[END_REF][START_REF] Elston | Polarization conversion from di↵raction gratings[END_REF][START_REF] Chandezon | Multicoated gratings, a di↵erential formalism applicable in the entire optical region[END_REF][START_REF] Li | Multilayer-coated di↵raction gratings: di↵erential method of Chandezon et al revisited[END_REF][START_REF] Cotter | Scattering-matrix approach to multilayer di↵raction[END_REF][START_REF] Granet | Scattering by a periodically corrugated dielectric layer with non-identical faces[END_REF][START_REF] Preist | Periodic multilayer gratings of arbitrary shape[END_REF][START_REF] Li | Some topics in extending the C-method to multilayer-coated gratings of di↵erent profiles[END_REF][START_REF] Granet | Di↵raction par des surfaces bipériodiques: résolution en coordonnées non-orthogonales[END_REF][START_REF] Granet | Analysis of di↵raction by surface-relief crossed gratings with use of the Chandezon method: Application to multilayer crossed gratings[END_REF]. The Maxwell's equations in covariant form leads to new pertubation methods [START_REF] Dusséaux | New pertubation theory of di↵raction gratings and its application to the study of ghosts[END_REF] and to models of two roughness levels [START_REF] Dusséaux | Model with two roughness levels for di↵raction gratings: the generalized Rayleigh expansion[END_REF]. For these works, the medium is linear, homogeneous and isotropic. E. Popov and M. Nevière have extended the C-method to the grating containing materials with nonzero susceptibility In most studies, the grating surface is described as a function and the study is done in the translation coordinate system. Plumey et al [START_REF] Plumey | Coordinate transformation method as applied to asymmetric gratings with vertical facets[END_REF] have extended the C-method to study gratings that are not described by functions. Granet et al studied gratings given by parametric equations [START_REF] Granet | La méthode des coordonnées curvilignes appliquée à la di↵raction par des réseaux dont le profil est donné par des équations paramétriques: application à la di↵raction par un réseau cycloidal[END_REF]. This formalism has given rise to some works in waveguide for the study of waveguide bends and power divider [START_REF] Dusséaux | Etude de transformateur plan-E dans un système de coordonnées non orthogonales[END_REF][START_REF] Dusséaux | Analysis of rectangular waveguide Hplane junctions in nonorthogonal coordinate system[END_REF][START_REF] Dusséaux | Analyse de composants plan-E symétriques en guides d'onde à section rectangulaire[END_REF][START_REF] Dusséaux | Telegraphist's equations for rectangular waveguides and analysis in nonorthogonal coordinates[END_REF]. The C-method is also an e cient theoretical tool for analyzing rough surfaces illuminated by a plane wave [START_REF] Benali | A new theory for scattering of electromagentic waves from conducting or dielectric rough surfaces[END_REF][START_REF] Dusséaux | Scattering of a plane wave by one-dimensional dielectric rough surfaces-study of the field in a nonorthogonal coordinate system[END_REF][START_REF] Baudier | Scattering of a plane wave by one-dimensional dielectric random rough surfaces -study with the curvilinear coordinate method[END_REF][START_REF] Braham | Scattering of electromagnetic waves from two-dimensional perfectly conducting random rough surfaces -study with the curvilinear coordinate method[END_REF][START_REF] Dusséaux | Implementation and validation of the curvilinear coordinate method for the scattering of electromagnetic waves from two-dimensional dielectric random rough surfaces[END_REF][START_REF] Braham | The curvilinear coordinate method associated with the short-coupling-range approximation for the study of scattering from onedimensional random rough surfaces[END_REF][START_REF] Dusséaux | Study of backscatter signature for seedbed surface evolution under rainfall -influence of radar precision[END_REF] or a electromagnetic beam [START_REF] Edee | Beam implementation in a nonorthogonal coordinate system: application to the scattering from random rough surfaces[END_REF][START_REF] Dusséaux | Eigenvalue system for the scattering from rough surfaces-saving in computation time by a physical approach[END_REF]. Recently, D. Prémel et al have implemented an original formulation based on the field-potential vectors and applied to the domain of low frequencies [START_REF] Prémel | Generalization of the second order vector potential formulation for arbitrary non-orthogonal curvilinear coordinates systems from the covariant form of Maxwell's equations[END_REF][START_REF] Caire | Semi-analytical computation of a quasi-static field induced by a 3D eddy current probe scanning a 2D layered conductor with parallel rough interfaces[END_REF].

The Maxwell's equations in covariant form and the translation system

We will derive the Maxwell's equation in covariant form in this subsection. For the onedimensional case, we have the surface function y = a(x). We consider the translation system:

8 > > > > > < > > > > > : x 0 = x y 0 = u = y a(x) z 0 = z (4.1)
Then the transformation matrix is:

A i i 0 = @x i @x 0 i = 0 B B B @ A x x 0 A x y 0 A x z 0 A y x 0 A y y 0 A y z 0 A z x 0 A z y 0 A z z 0 1 C C C A = 0 B B B @ 1 0 0 da dx 1 0 0 0 1 1 C C C A (4.2) A i 0 i = @x i 0 @x i = 0 B B B @ 1 0 0 da dx 1 0 0 0 1 1 C C C A (4.
3)

The covariant basis vectors can be expressed from the basis vectors of the Cartesian coordinate system (ũ x , ũy , ũz ):

8 > > > > > < > > > > > : ũx 0 = ũx + da dx ũy ũy 0 = ũy ũz 0 = ũz (4.4)
and the contravariant basis vectors can be written as:

8 > > > > > < > > > > > : ũx 0 = ũx ũy 0 = da dx ũx + ũy ũz 0 = ũz (4.5)
The covariant and contravariant metric tensors are [START_REF] Stratton | Electromagnetic Theory. International series in physics[END_REF]:

g i 0 j 0 = X i,j A i i 0 A j j 0 g ij = 0 B B B @ 1 + ( da dx ) 2 da dx 0 da dx 1 0 0 0 1 1 C C C A (4.6) g i 0 j 0 = X i,j A i 0 i A j 0 j g ij = 0 B B B @ 1 da dx 0 da dx 1 + ( da dx ) 2 0 0 0 1 1 C C C A (4.7)
So in the new coordinate system, the covariant components of a vector ṽ can be written as,

8 > > > > > < > > > > > : v x 0 = v x + da dx v y v y 0 = v y v z 0 = v z (4.8)
and the contravariant components as,

8 > > > > > < > > > > > : v x 0 = v x v y 0 = da dx v x + v y v z 0 = v z (4.9)
The covariant components v y 0 and v z 0 become identified with Cartesian ones v y and v z .

Moreover, the covariant component v x 0 and v y 0 are parallel to the interface given by u = 0 (i.e. y = a(x)).

By a similar procedure, for the two-dimensional surface z = a(x, y), we have by using the translation system (x 0 , y 0 , z 0 ) = (x, y, z a(x, y)):

A i i 0 = 0 B B B @ 1 0 0 0 1 0 @a @x @a @y 1 1 C C C A (4.10) A i 0 i = 0 B B B @ 1 0 0 0 1 0 @a @x @a @y 1 1 C C C A (4.11) g i 0 j 0 = 0 B B B @
1 + ( @a @x ) 2 @a @x @a @y @a @x @a @x @a @y 1 + ( @a @y ) 2 @a @y @a @x @a @y

1 1 C C C A (4.
12)

g i 0 j 0 = 0 B B B @
1 0 @a @x 0 1 @a @y @a @x @a @y 1 + ( @a @x ) 2 + ( @a @y )

2 1 C C C A (4.13)
So, in the translation coordinate system, the covariant components of a vector ṽ can be written as, 8 > > > > > < > > > > > :

v x 0 = v x + @a @x v z v y 0 = v y + @a @y v z v z 0 = v z (4.14)
and the contravariant component as,

8 > > > > > < > > > > > : v x 0 = v x v y 0 = v y v z 0 = v z @a @x v x @a @y v y (4.15)
The covariant component v z 0 is simply the vertical component v z . Moreover, the covariant components v x 0 and v y 0 are parallel to the surface coordinate z 0 = 0 (i.e. z = a(x, y)).

From the Ostrogradsky theorem and the Stokes' theorem expressed in a non orthogonal coordinate system (x i 0 , x j 0 , x k 0 ), the Maxwell's equations for time-harmonic fields can be written as [START_REF] Stratton | Electromagnetic Theory. International series in physics[END_REF]:

1 p g 0 X i 0 @ @x i 0 ( p g 0 B x i 0 ) = 0 (4.16) 1 p g 0 ( @E x k 0 @x j 0 @E x j 0 @x k 0 ) = j!B x i 0 (4.17) 1 p g 0 X i 0 @ @x i 0 ( p g 0 D x i 0 ) = 0 (4.18) 1 p g 0 ( @H x k 0 @x j 0 @H x j 0 @x k 0 ) = j!D x i 0 (4.19)
where g 0 = det(g i 0 j 0 ). Here, we assume that there is no current density and no charge density. For a linear, homogeneous, isotropic and non magnetic medium, the constitutive relations for time-harmonic fields can be written as:

D x i 0 = " c E x i 0 = " c 3 X j 0 =1
g i 0 j 0 E x j 0 (4.20)

B x i 0 = µ 0 H x i 0 = µ 0 3 X j 0 =0 g i 0 j 0 E x j 0 (4.21)
So the equations (4.17) and (4.19) associated with the constitutive relations (4.20) and (4.21) can be written as:

@E x k 0 @x j 0 @E x j 0 @x k 0 = j!µ 0 X j 0 g i 0 j 0 H x j 0 (4.22)
and @H x k 0 @x j 0 @H x j 0 @x k 0 = j!" c X j 0 g i 0 j 0 E x j 0 (4.23) more specifically, equations (4.22) and (4.23) are equivalent to the following six equations [START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF]:

@E z 0 @y 0 @E y 0 @z 0 = j!µ 0 (g x 0 x 0 H x 0 + g x 0 y 0 H y 0 + g x 0 z 0 H z 0 ) (4.24) @E z 0 @x 0 @E x 0 @z 0 = j!µ 0 (g y 0 x 0 H x 0 + g y 0 y 0 H y 0 + g y 0 z 0 H z 0 ) (4.25) @E y 0 @x 0 @E x 0 @y 0 = j!µ 0 (g z 0 x 0 H x 0 + g z 0 y 0 H y 0 + g z 0 z 0 H z 0 ) (4.26) @H z 0 @y 0 @H y 0 @z 0 = j!" c (g x 0 x 0 E x 0 + g x 0 y 0 E y 0 + g x 0 z 0 E z 0 ) (4.27) @H z 0 @x 0 @H x 0 @z 0 = j!" c (g y 0 x 0 E x 0 + g y 0 y 0 E y 0 + g y 0 z 0 E z 0 ) (4.28)
@H y 0 @x 0 @H x 0 @y 0 = j!" c (g z 0 x 0 E x 0 + g z 0 y 0 E y 0 + g z 0 z 0 E z 0 ) (4.29)

Formulation for one-dimensional case

We first consider the one-dimensional case. As we have described before, a surface by equation y = a(x) separates two di↵erent media. It is illuminated by a monochromatic plane wave with wavelength under incident angle ✓ 0 . The incident wave vector k0 is defined by the incident angle ✓ 0 .

k0 = ↵ 0 ũx + 0 ũy (4.30)
with ↵ 0 = k sin ✓ 0 , 0 = k cos ✓ 0 . The surface could be periodic or non periodic. Here we consider periodic surfaces or periodic random surfaces.

We represent the vector function by its complex vector function and omit its timedependence factor exp(j!t). So for the horizontal (E // ) polarization and vertical (H // ) polarization,

F 0 (x, y) = exp( j↵ 0 x + j 0 y) = 8 > < > : E 0,z (x, y) E // Z 1 H 0,z (x, y) H // (4.31)
and

Z 1 H = k0 k ^Ẽ (4.32)
where

Z 1 = 120⇡ • Ohms.
The reflected and transmitted plane waves can be written in a similar form. But, for rough surface, we have, in addition to the incident, reflected and transmitted plane waves, a scattered field F (x, y) because of the deformation. The problem consists in working out the scattered field within the two media. The rough surface here is generated by simulation.

Equations (4.24) to (4.29) enable us to write Maxwell's equations associated with the constitutive relations:

8 > > > > > < > > > > > : j k 1 @F (x,u) @u = j k 1 b(x) @F (x,u) @x + c(x)G(x, u) j k 1 @G(x,u) @u = 1 k 2 1 @
@x (c(x) @F (x,u) @x ) + ⌫ 2 F (x, u)

+ j k 1 @ @x (b(x)G(x, u)) (4.33) with b(x) = da dx 1 + ( da dx ) 2 , c(x) = 1 1 + ( da dx ) 2
and in medium (m),

F (x, u) = F m (x, u), G(x, u) = G m (x, u), ⌫ = ⌫ m , m = 1, 2. In E // polarization, F (x, u) = E z 0 (x, u), G(x, u) = Z 1 H x 0 (x, u). In H // polarization, F (x, u) = Z 1 ⌫ H z 0 (x, u), G(x, u) = ⌫E x 0 (x, u).
System (4.33) can be written in the form as follows:

j k 1 @ (x, u) @u = L (x, u) (4.34)
with

L = 0 @ j k 1 b(x) @• @x c(x)• 1 k 2 1 @ @x (c(x) @• @x ) + ⌫ 2 • j k 1 @b(x)• @x 1 A (4.35) and (x, u) = 0 @ F (x, u) G(x, u) 1 A (4.36)
We separate the variables by writing (x, u) = '(x)(u), then we get:

j k 1 (u) d(u) du = L'(x) '(x) = r = constant (4.37)
So we conclude that (u) = Cexp( jk 1 ru), L'(x) = r'(x) and (x, u) = Aexp( jk 1 ru)'(x).

If the function a(x) is a period function with D its period, then one has,

8 > > > > > > < > > > > > > : a(x) = P m a m exp( j2⇡mx/D) f (x) = P m f m exp( j↵ m x) g(x) = P m g m exp( j↵ m x) (4.38) with '(x) = 0 @ f (x) g(x) 1 A and ↵ m = k 1 sin✓ 0 + m 2⇡ D .
Under this function decomposition, the eigenproblem L'(x) = r'(x) has a matrix form as follows:

[L]' = r ' (4.39)

with

[L] = 0 @ [L ff ] [L fg ] [L gf ] [L gg ] 1 A and ' = 0 @ f g 1 A (4.40)
where

[L ff ] = [C][ Ȧ][↵], [L fg ] = [C], [L gf ] = ⌫ 2 [I] [↵][C][↵], [L gg ] = [↵][C][ Ȧ], [ Ȧ] pq = ȧp q = (p q) 2⇡ d a p q , [I] pq = p q , [C] = ([I]+[ Ȧ][ Ȧ]) 1 , ↵p = ↵p k 1 , [↵] pq = p q ↵p , ( f ) p = f p , (g) p = g p , 8(p, q) 2 Z 2 .
Equations (4.39) and (4.40) give an eigenvalue system of infinite dimension. In a numerical computation, one can truncate it to a finite order problem with a truncation order M . Theoretically, increasing the truncation order M will increase the precision of results as well as increase the computational time. We want to ensure a certain precision and also keep M relatively small. For the lossless medium, i.e. the medium with optical index ⌫ real, the power balance criterion (3.12) is checked to see if the truncation order M is large enough for a certain precision. The new system is similar to the original one except that now, we have M  p, q  M . By solving the truncated eigenvalue problem, with the eigenvalues r n and eigenvectors 'n , 1 6 = n 6 = 2M + 1 one gets:

8 > > > > > > > > > > < > > > > > > > > > > : F n (x, u) = f n (x)exp( jk 1 r n u) = P M mM f mn exp( j↵ m x)exp( jk 1 r n u) G n (x, u) = g n (x)exp( jk 1 r n u) = P M mM g mn exp( j↵ m x)exp( jk 1 r n u) (4.41)
So we are left with the eigenproblem of order 4M + 2. The signs of the real and imaginary parts of the eigenvalues r n define the nature of the wave corresponding to the elementary wavefunction. In particular, the associated expression represent an outgoing wave propagating with no attenuation if Re(r n ) > 0 and Im(r n ) = 0. For an evanescent wave, Im(r n ) < 0. Finally, the field scattered in the air can be represented as a linear combination of all the solutions that verifies the outgoing conditions.

(i) (x, u) = 2M +1 X n=1 C (i) n (i) n (x, u), i = 1, 2 (4.42)
and the amplitudes

C (i)
n are determined by solving the boundary conditions at u = 0 (i.e., at y = a(x)). The boundary conditions stipulate the continuity of the electric and magnetic components parallel to the surface. These components are (H x 0 , E z 0 ) in E // polarization and (E x 0 , H z 0 ) in H // polarization.

Formulation for the two-dimensional case

Now we consider the two-dimensional case. We assume the surface z = a(x, y) is illuminated by a monochromatic plane wave with wavelength . a(x, y) is a local function with L denotes the deformation length with respect to the Ox and Oy axis. For the formulation applied to crossed gratings, we refer the readers to [START_REF] Granet | Di↵raction par des surfaces bipériodiques: résolution en coordonnées non-orthogonales[END_REF][START_REF] Granet | Analysis of di↵raction by surface-relief crossed gratings with use of the Chandezon method: Application to multilayer crossed gratings[END_REF]. The incident wave vector k0 is defined by the zenith angle ✓ 0 and the azimuth angle ' 0 .

k0 = ↵ 0 ũx + 0 ũy 0 ũz (4.43)
with

↵ 0 = k sin ✓ 0 cos ' 0 , 0 = k sin ✓ 0 sin ' 0 , 0 = k cos ✓ 0 (4.44) and k = 2⇡ (4.45)
For the E // polarization, the incident field can be expressed as:

Ẽ0 (x, y, z) = he j k0 r and Z H0 = k0 k ⇥ Ẽ0 (4.46)
For the H // polarization, the incident field can be expressed as:

Z H0 (x, y, z) = he j k0 r and Z H0 = k0 k ⇥ Ẽ0 (4.47)
Here

h = sin ' 0 ũx + cos ' 0 ũy (4.48)
and

r = xũ x + yũ y + zũ z (4.49)
We want to know the scattered field, but it cannot be expressed by the Rayleigh integral (3.43) in the modulated zone if the perturbation amplitude is too large. We can obtain an expression of field that is valid over the surface by solving Maxwell's equation in the translation coordinate system:

8 > > > > > <
> > > > > :

x 0 = x y 0 = y z 0 = z a(x, y) (4.50)
In a source-free medium, from equations (4.24)-(4.29), we can obtain that the longitudinal components E z 0 and ZH z 0 obey to the propagation equation [START_REF] Granet | Analysis of di↵raction by surface-relief crossed gratings with use of the Chandezon method: Application to multilayer crossed gratings[END_REF]: @ @z 0 (g x 0 z 0 @ @x 0 + @g x 0 z 0 @x 0 ) @ @z 0 (g y 0 z 0 @ @y 0 + @g y 0 z 0 @y 0 ) + jkg z 0 z 0 @ 0 @z 0 = @ 2 @x 02 + @ 2 @y 02 + k 2 and (x 0 , y 0 , z 0 ) = E z 0 (x 0 , y 0 , z 0 ) or ZH z 0 (x 0 , y 0 , z 0 ). And g x 0 z 0 , g y 0 z 0 and g z 0 z 0 are elements of metric tensor which depend on the derivatives of function a(x 0 , y 0 ) with respect to x 0 and y 0 . From equation (4.13), we have:

8 > > > > > < > > > > > : g x 0 z 0 = @a @x 0 g y 0 z 0 = @a @y 0 g z 0 z 0 = 1 + ( @a @x 0 ) 2 + ( @a @y 0 ) 2 (4.53)
Again from (4.24)-(4.29), we obtain expression of components E x 0 , E y 0 , H x 0 and H y 0 in terms of longitudinal components E z 0 and ZH z 0 only.

@ 2 E x 0 @z 02 + k 2 E x 0 = @ 2 E z 0 @x 0 @z 0 k 2 g x 0 z 0 E z 0 jkg y 0 z 0 @ZH z 0 @z 0 jk @ZH z 0 @y 0 (4.54) @ 2 E y 0 @z 02 + k 2 E y 0 = @ 2 E z 0 @y 0 @z 0 k 2 g y 0 z 0 E z 0 jkg x 0 z 0 @ZH z 0 @z 0 + jk @ZH z 0 @x 0 (4.55)
@ 2 ZH x 0 @z 02 + k 2 ZH x 0 = @ 2 ZH z 0 @x 0 @z 0 k 2 g x 0 z 0 ZH z 0 jkg y 0 z 0 @E z 0 @z 0 jk @E z 0 @y 0 (4.56) @ 2 ZH y 0 @z 02 + k 2 ZH y 0 = @ 2 ZH z 0 @y 0 @z 0 k 2 g y 0 z 0 ZH z 0 jkg x 0 z 0 @E z 0 @z 0 + jk @E z 0 @x 0 (4.57)

The covariant components E x 0 and E y 0 are parallel to the interface. Consequently, for instance, for a perfectly conducting surface, we have E x 0 = E y 0 = 0 at z 0 = 0. We need to solve the propagation equation (4.51).

To solve equation (4.51), we use a Fourier transform with respect to x 0 and y 0 , then the equations take the following form Here, K⇤ L is the convolution product of two Fourier transforms K(↵, , z 0 ) and L(↵, , z 0 ). Now, convolution products are approximated as follows:

@ @z 0 [j↵(ĝ x 0 z 0 ⇤ ˆ ) + jĝ x 0 z 0 ⇤ (↵ ˆ ) + j (ĝ y 0 z 0 ⇤ ˆ ) + jĝ y 0 z 0 ⇤ ( ˆ )] + jkĝ z 0 z 0 ⇤ @ ˆ 0 @z 0 = 2 ˆ (4.58) with j k @ ˆ @z 0 = ˆ 0 (4.
( K ⇤ L)(↵, , z 0 ) = 1 4⇡ 2 Z 1 1 Z 1 1 K(↵ 0 , 0 , z 0 ) L(↵ ↵ 0 , 0 )d↵ 0 d 0 ⇡ ↵ 2 4⇡ 2 X p X q K(↵ p , q , z 0 ) L(↵ ↵ p , q ) (4.61)
where

↵ p = k sin ✓ 0 cos ' 0 + p ↵, q = k sin ✓ 0 sin ' 0 + q ↵ (4.62)
and ↵ = = 2⇡ L is the spectral resolution. Using this approximation and applying the point matching method at discrete values (↵ s , t ) to equation (4.58), we obtain

j k @ @z 0 ( X p,q ( ↵ s k ĝx 0 z 0 s p,t q + ĝx 0 z 0 s p,t q ↵ a k + t k ĝy 0 z 0 s p,t q + ĝy 0 z 0 s p,t q q k ) ˆ (↵ p , q , z 0 )) + j k @ @z 0 ( X p,q ĝz 0 z 0 s p,t q ˆ 0 (↵ p , q , z 0 )) = 2 st k 2 ˆ (↵ s , t , z 0 ) (4.63) j k @ ˆ (↵ s , t , z 0 ) @z 0 = ˆ 0 (↵ s , t , z 0 ) (4.64) with ĝx 0 z 0 p,q = ↵ 2 4⇡ 2 ĝx 0 z 0 (↵ p , q ) (4.65) ĝy 0 z 0 p,q = ↵ 2 4⇡ 2
ĝy 0 z 0 (↵ p , q ) (4.66)

ĝz 0 z 0 p,q = pq + X u,v ĝx 0 z 0 p u,q v ĝx 0 z 0 u,v + X u,v ĝy 0 z 0 p u,q v ĝy 0 z 0 u,v (4.67)
Equation (4.63) can be written in matrix form according to the sample theorem [START_REF] Charbit | Eléments de théorie du signal: les signaux aléatoires[END_REF], the elementary wave functions ˆ mn (↵, , z 0 ) and ˆ 0 mn (↵, , z 0 ) can be constructed from mn and 0 mn ˆ mn (↵, , z 0 ) =exp( jkr mn z 0 )

j k [L l ] @ @z 0 0 @ ~ ~ 0 1 A = [L r ] 0 @ ~ ~ 0 1 A (4.
⇥ s=M X s= M t=M X t= M mn (↵ s , t )sinc( ⇡ ↵ (↵ ↵ s ))sinc( ⇡ ↵ ( t )) (4.71) ˆ 0 mn (↵, , z 0 ) =exp( jkr mn z 0 ) ⇥ s=M X s= M t=M X t= M 0 mn (↵ s , t )sinc( ⇡ ↵ (↵ ↵ s ))sinc( ⇡ ↵ ( t )) (4.72)
Finally, the Fourier transform of Oz-component is defined as a linear combination of M s eigensolutions satisfying the outgoing wave condition:

ˆ d (↵, , z 0 ) = X (m,n)2Ds A mn ˆ mn (↵, , z 0 ) (4.73) ˆ 0 d (↵, , z 0 ) = X (m,n)2Ds
A mn ˆ 0 mn (↵, , z 0 ) (4.74)

Substituting E z 0 = 0 and applying the same procedure in the spectral domain, we obtain the Fourier transforms of horizontal polarized transverse components:

ˆ (ha) dT (↵, , z 0 ) = X (m,n)2Ds A (ha) mn ˆ (ha)
T,mn (↵, )exp( jkr mn z 0 ) (4.75)

with

ˆ (ha) dT (↵, , z 0 ) = 0 B B B B B B @ Ê(ha) dx 0 (↵, ) Ê(ha) dy 0 (↵, ) Z Ĥ(ha) dx 0 (↵, ) Z Ĥ(ha) dy 0 (↵, ) 1 C C C C C C A and ˆ (ha) T,mn (↵, ) = 0 B B B B B B @ Ê(ha) x 0 ,mn (↵, ) Ê(ha) y 0 ,mn (↵, ) Z Ĥ(ha) x 0 ,mn (↵, ) Z Ĥ(ha) y 0 ,mn (↵, ) 1 C C C C C C A (4.76)
According to the sampling theorem [START_REF] Charbit | Eléments de théorie du signal: les signaux aléatoires[END_REF], we write

ˆ (ha) T,mn (↵, ) = M X s= M M X t= M ˆ (ha) T,mn (↵ s , t )sinc( ⇡ ↵ (↵ ↵ s ))sinc( ⇡ ↵ ( t )) (4.77)
where

Ê(ha) x 0 ,mn (↵ s , t ) = k 2 M X p= M M X q= M ĝy 0 z 0 s p,t q 0 mn (↵ p , q ) k t mn (↵ s , t ) (4.78) Ê(ha) y 0 ,mn (↵ s , t ) = k 2 M X p= M M X q= M ĝx 0 z 0 s p,t q 0 mn (↵ p , q ) + k↵ s mn (↵ s , t ) (4.79)
Z Ĥ(ha) x 0 ,mn (↵ s , t ) = k 2 M X p= M M X q= M
ĝx 0 z 0 s p,t q mn (↵ p , q ) k↵ s 0 mn (↵ s , t ) (4.80)

Z Ĥ(ha) y 0 ,mn (↵ s , t ) = k 2 M X p= M M X q= M
ĝy 0 z 0 s p,t q mn (↵ p , q ) k t 0 mn (↵ s , t ) (4.81)

Taking H z 0 = 0 and substituting Ê(ha) by Z Ĥ(va) and Z Ĥ(ha) by Ê(va) , we obtain the vertical components of magnetic and electric fields.

The scattering amplitudes A (ha) mn and A (va) mn are found by solving the boundary conditions.

So, for an incident wave in (a) polarization, we can write

E (1,ha) dx 0 (x 0 , y 0 , z 0 ) + E (1,va) dx 0 (x 0 , y 0 , z 0 ) E (2,ha) dx 0 (x 0 , y 0 , z 0 ) E (2,va) dx 0 (x 0 , y 0 , z 0 ) = (E (a) 0x 0 (x 0 , y 0 , z 0 ) + ⇢ (a) r E (a) rx 0 (x 0 , y 0 , z 0 )) ⇢ (a) t E (a) tx 0 (x 0 , y 0 , z 0 ) (4.82) E (1,ha) dy 0 (x 0 , y 0 , z 0 ) + E (1,va) dy 0 (x 0 , y 0 , z 0 ) E (2,ha) dy 0 (x 0 , y 0 , z 0 ) E (2,va) dy 0 (x 0 , y 0 , z 0 ) = (E (a) 0y 0 (x 0 , y 0 , z 0 ) + ⇢ (a) r E (a) ry 0 (x 0 , y 0 , z 0 )) ⇢ (a) t E (a)
ty 0 (x 0 , y 0 , z 0 ) (4.83)

H (1,ha) dx 0 (x 0 , y 0 , z 0 ) + H (1,va) dx 0 (x 0 , y 0 , z 0 ) H (2,ha) dx 0 (x 0 , y 0 , z 0 ) H (2,va) dx 0 (x 0 , y 0 , z 0 ) = (H (a) 0x 0 (x 0 , y 0 , z 0 ) + ⇢ (a) r H (a) rx 0 (x 0 , y 0 , z 0 )) ⇢ (a) t H (a)
tx 0 (x 0 , y 0 , z 0 ) (4.84)

H (1,ha) dy 0 (x 0 , y 0 , z 0 ) + H (1,va) dy 0 (x 0 , y 0 , z 0 ) H (2,ha) dy 0 (x 0 , y 0 , z 0 ) H (2,va) dy 0 (x 0 , y 0 , z 0 ) = (H (a) 0y 0 (x 0 , y 0 , z 0 ) + ⇢ (a) r H (a)
ry 0 (x 0 , y 0 , z 0 )) ⇢ are the Fresnel reflection and transmission coe cients. After a Fourier transform, the point matching method is applied, then a 4M s -dimensional matrix system is obtained, the inversion of which leads to scattering amplitude A (ha) mn and A (va) mn .

These scattering amplitudes lead to the bistatic coe cients as defined in equation (3.60).

Conclusion

In this chapter, we show the C-method and how this C-method leads to eigenvalue problem. We present the formulations for both one-dimensional case and two-dimensional case. The computational time of the C-method is a key topic of our research.

The computational time of the C-method is mainly spent on the computation of eigenvalues and eigenvectors. We give a figure here to show the computation time of the C-method. The figure 4.1 shows the computational time of numerical experiment of one realisation. The perfectly conducting surface we consider is of 64 square wavelengths and a = and l x = l y = 1.41 . The incident angle is chosen as ✓ 0 = 30 and ' 0 = 0 .

The computational time varies as N 3 where N = 2(2M + 1) 2 is the order of the considered matrix. In fact, from figure 4.1, we can see that we have approximately the computational time t = a(2(2M + 1) 2 ) ↵ , with ↵ ⇡ 3.1 a good fit. In terms of computational time, the C-method is not competitive with respect to fast integral method whose complexity is O(N log N ) [START_REF] Pak | Numerical simulations and backscattering enhancement of electromagentic waves from two-dimensional dielectric random rough surfaces with the sparse-matrix canonical grid method[END_REF][START_REF] Wagner | Monte Carlo simulations of electromagentic scattering from two-dimensional random rough surfaces[END_REF][START_REF] Jandhyala | Fast algorithm for the analysis of scattering by dielectric rough surfaces[END_REF][START_REF] Torrungrueng | Numerical studies of backscattering enhancement of electromagnetic waves from two-dimensional random rough surfaces with the forward-backward novel spectral acceleration method[END_REF][START_REF] Soriano | Scattering of electromagentic waves from twodimensional rough surfaces with an impedance approximation[END_REF][START_REF] Tsang | Wave scattering with the UV multilevel partitioning method: 2. three-dimensional problem of nonpenetrable surface scattering[END_REF]. The computational time is a weak point of the C-method, in particular, for analyzing rough random surfaces insofar as the average scattered intensity is estimated over results of several surface realizations (Monte-Carlo method). However, the strength of the C-method is that it leads to the eigensolutions of the scattering problem. It is an accurate method and it can be used as a reference for the analytical methods [START_REF] Elfouhaily | A critical survey of approximate scattering wave theories from random rough surfaces[END_REF].

In the next chapter, we propose a parallel QR algorithm adapted to the C-method for reducing the computational time. The proposed method keeps the strength of C-method and improves the weak point of C-method. [START_REF] Bai | Templates for the solution of algebraic eigenvalue problems[END_REF] have been developed to deal with large-scale eigenvalue problems. However, they have the possibility of missing some eigenvalues. So these iterative methods are ine↵ective for the C-method because all the eigenvalues and eigenvectors are needed. In contrast, the QR algorithm, which is based on similarity transformations, calculates all the eigenvalues and eigenvectors with very little danger, and only with a warning of missing some eigensolutions. We propose a specifically designed parallel QR algorithm for the C-method to save the computation time. Three techniques are used in the implementation of parallel QR algorithm: early shift, parallel bulge chasing and parallel aggressive early deflation (AED). The early shifts are introduced in parallel algorithm to give approximation of a part of the eigenvalues of the matrix. The early shifts are based on physical interpretation and observation based on the C-method and they are specifically designed and first introduced in our work. They provide the possibility of quick deflation. For the bulge chasing, instead of only a single bulge, containing two shifts, a chain of several tightly coupled bulges, each containing two shifts, is chased in the course of one multishift QR iteration. As described in [START_REF] Bai | On a block implementation of hessenberg multishift qr iteration[END_REF], this idea allows performing most of the computational work in terms of matrix-matrix multiplications to benefit from highly e cient level 3 BLAS. The idea of AED allows to detect converged eigenvalues much earlier than conventional deflation strategies. We will first present QR sequential algorithm and the shift strategy to accelerate the convergence, then we present all the parallel techniques for this specifically designed QR algorithm.

The basic QR algorithm

The QR algorithm computes a Schur decomposition of a matrix. It is certainly one of the most important algorithms in eigenvalue computations. As QR seems to be the only method that can provide us all the eigenvalues and eigenvectors, we choose to use the QR algorithm.

The QR algorithm consists of two separate stages. First, by means of a similarity transformation, the original matrix is transformed in a finite number of steps to Hessenberg form. This first stage of algorithm prepares its second stage, the actual QR iterations that are applied to the Hessenberg matrix. The overall complexity (number of floating points) of the algorithm is O(N 3 ) where the matrix A is assumed to be of the order

N ⇥ N .
We start with a basic iteration, given by algorithm 1. We notice that: Algorithm 1 Basic QR algorithm Input: A 2 C N ⇥N Output: An upper triangular matrix T and a unitary matrix U such that A = UT U ⇤ is the Schur decomposition of A.

A k 1 = Q k R k 4: A k = R k Q k 5: U k = U k 1 Q k 6: end for 7: Set T = A 1 and U = U 1 A k = R k Q k = Q ⇤ k A k 1 Q k (5.1)
and hence A k and A k 1 are unitary similar. The matrix sequence {A k } converges (under certain assumptions) towards an upper triangular matrix [START_REF] Golub | Matrix Computations[END_REF]. Let us assume that the eigenvalues are pairwise di↵erent in magnitude and we can therefore number the eigenvalues such that

| 1 | > | 2 | > ... > | N |.
Then the elements of A k below the diagonal converge to zero like [START_REF] Golub | Matrix Computations[END_REF]:

|a (k) ij | = O(| i j | k ), i > j (5.2)
From equation (5.1), we have:

A k = Q ⇤ k A k 1 Q k = Q ⇤ k Q ⇤ k 1 A k 2 Q k 1 Q k = Q ⇤ k ...Q ⇤ 1 A 0 Q 1 ...Q k (5.3)
With the same assumption on the eigenvalues, A k tends to an upper triangular matrix and U k = Q 1 ...Q k converges to the matrix of Schur vectors.

The convergence of the basic QR algorithm is slow and expensive. We want to:

• find a matrix structure that is preserved by the QR algorithm and that lowers the cost of a single iteration step.

• improve the convergence properties of the algorithm.

The desired matrix structure is a Hessenberg matrix: a matrix H is a Hessenberg matrix if its elements below the lower o↵-diagonal are zero, h ij = 0 for i > j+1. The Hessenberg form is preserved by the QR algorithm and this form can lower the cost of a single iteration step [START_REF] Golub | Matrix Computations[END_REF].

There are several means of Hessenberg reduction such as Gram-Schmidt transformation, Householder reduction and Givens rotations [START_REF] Golub | Matrix Computations[END_REF]. An E cient parallel algorithm for this Hessenberg reduction is implemented in the ScaLAPACK [80] (Scalable Linear Algebra PACKage) routine PZGEHRD. So, we will focus on the iterative part that comes after this Hessenberg reduction and try to improve the convergence properties of the algorithm.

QR algorithm with shift

We will show how the convergence of the Hessenberg QR algorithm can be improved dramatically by introducing spectral shifts into the algorithm.

Lemma 1. Let H be an irreducible Hessenberg matrix, i.e., h i+1,i 6 = 0 for all i = 1, ..., N

1. Let H = QR be the factorization of H. Then for the diagonal elements of R, we have

|r kk | > 0, for all k < N. Thus, if H is singular then r NN = 0.
This lemma gives the motivation of shift strategy to speed up the convergence of the QR algorithm. To see this, assume that is an eigenvalue of the irreducible Hessenberg matrix H. We perform:

Algorithm 2 The single shift QR algorithm (one iteration)

1: H I = QR 2: H = RQ + I
We can see that H is similar to H:

H = Q ⇤ (H I)Q + I = Q ⇤ HQ (5.4)
By the lemma, we have:

H I = QR, with R NN = 0 (5.5)
So,

H = RQ + I = 0 @ H 1 h 1 0 1 A (5.6)
So if we apply a QR step with a perfect shift to a Hessenberg matrix, the eigenvalue drops out. We then have a deflation, i.e. we can proceed the algorithm with a smaller matrix H 1 of size (N 1) ⇥ (N 1).

For the single shift, when the item h N 1,N 1 is O(h N,N 1 ), the convergence could be slow even the Rayleigh quotient shift gives a very good approximation (e.g. h N,N 1 is very small). In practice, the double shift QR algorithm is very commonly used for real matrices and can be extended to complex matrices [START_REF] Mark | A parallel eigenvalue routine for complex Hessenberg matrices[END_REF]. The algorithm is characterized by a "bulge chasing" procedure.

Suppose 1 and 2 are two shifts of the Hessenberg matrix H. The algorithm proceeds as follows:

1. Calculate the first column of the shift polynomial

v = (H 1 I)(H 2 I)e 1 = 0 B B B B B B B B B B B B B @ ⇤ ⇤ ⇤ 0 . . . 0 1 C C C C C C C C C C C C C A (5.7)
2. Construct a 3 ⇥ 3 Householder transformation Q 1 such that the second and third entries of v are transformed to zero. The similarity transformation gives the updated matrix H 1 :

H 1 = Q ⇤ 1 HQ 1 = 0 B B B B B B B B B B B B B B B B @ ⇤ ⇤ ⇤ ⇤ ⇤ ⇤ ••• X X X ⇤ ⇤ ⇤ ••• X X X ⇤ ⇤ ⇤ ••• X X X ⇤ ⇤ ⇤ ••• 0 0 0 ⇤ ⇤ ⇤ ••• 0 0 0 0 ⇤ ⇤ ••• . . . . . . . . . . . . . . . . . . 1 C C C C C C C C C C C C C C C C A (5.8)
The Hessenberg structure is damaged by the bulge that we denote with symbol "X".

3. Construct a 3 ⇥ 3 Householder transformation Q 2 such that the third and fourth entries of the first column of H 1 reduce to zero. The similarity transformation gives the updated matrix H 2 : The bulge will be chased to vanish at the bottom right corner and lead to zeros, thus deflations. For example, H 3 will be look like as follows:

H 2 = Q ⇤ 2 H 1 Q 2 = 0 B B B B B B B B B B B B B B B B @ ⇤ ⇤ ⇤ ⇤ ⇤ ⇤ ••• ⇤ ⇤ ⇤ ⇤ ⇤ ⇤ ••• 0 X X X ⇤ ⇤ ••• 0 X X X ⇤ ⇤ ••• 0 X X X ⇤ ⇤ ••• 0 0 0 0 ⇤ ⇤ ••• . . . . . . . . . . . . . . . . . . 1 C C C C C C C C C C C C C C C C A (
H 3 = Q ⇤ 3 H 2 Q 3 = 0 B B B B B B B B B B B B B B B B @ ⇤ ⇤ ⇤ ⇤ ⇤ ⇤ ••• ⇤ ⇤ ⇤ ⇤ ⇤ ⇤ ••• 0 ⇤ ⇤ ⇤ ⇤ ⇤ ••• 0 0 X X X ⇤ ••• 0 0 X X X ⇤ ••• 0 0 X X X ⇤ ••• . . . . . . . . . . . . . . . . . . 1 C C C C C C C C C C C C C C C C A
(5.10)

Early shift

Some of the eigenvalues of the scattering matrix can be approximated as follows: for the one-dimensional case,

r (m) n = ± q (⌫ (m) ) 2 ↵2 n (5.11)
with ↵n = ↵ n /k (1) , and for the two-dimensional case,

r (m) pq = ± q (⌫ (m) ) 2 ↵2 p ˜ 2 q
(5.12)

with ↵p = ↵ p /k (1) , ˜ q = q /k (1) . The values r

(m) n or r (m)
pq constitute very good approximations when the index n or p, q are small relative to the matrix size. In general if ⌫ (m) is real or has a very small imaginary part, for the n that satisfies ↵2 n < Re(⌫ (m) ) 2 , or for the pair (p, q) that satisfies ↵2 p + ˜ 2 q < Re(⌫ (m) ) 2 , the approximation can be quite accurate. For a lossless medium (m), a real value r We take the one-dimensional case as an example to show why r (m) n can be a good approximation. In fact, we can check that the functions is just the set of eigenvalues of (4.37). And the problem under consideration (4.39) is just a truncated form of the infinite dimensional eigenvalue problem (4.37). It is therefore attempting to take advantage of this analytical solution to represent the solution as a linear combination of F± , however, if we operate in such a way, we can show that this method is equivalent to the well known Rayleigh expansion method, which leads to a numerical failure [START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF].

F± (x, u) = exp( jk (m) ↵ n x ± jk (m) r (m) n u) (5.
We therefore propose to use (5.11) and (5.4) as shifts. Moreover, according to our observation, we use only these approximations when ⌫ (m) is real or has a very small imaginary part, and n satisfies ↵2 n < Re(⌫ (m) ) 2 , or for the pair (p, q) that satisfies ↵2 p + ˜ 2 q < Re(⌫ (m) ) 2 . The approximations come in pair and each used pair will be used only once to create a 3 ⇥ 3 bulge.

We include an example for the one-dimensional case, where M = 15, ⌫ = 1, = 1, d = 10.5, ✓ = 2⇡ 9 . The surface used is a generated Gaussian rough surface with correlation length l = and standard deviation of height is = 0.2 . The eigenvalues of the scattering matrix are listed in the first column and the second column of the table 5.1. The eigenvalues in the first column correspond to the outgoing wave and the eigenvalues in the second column correspond to the incoming wave. The values of the r n = ± p ⌫ 2 ↵2 n , where M < n < M are listed in the third column of the table 5.1: 0.0006 + 0.8504i ±0.8454i -0.0002 -0.6893i -0.0002 + 0.6893i ±0.6887i -0.0000 -0.5021i -0.0000 + 0.5021i ±0.5021i 0.0000 -0.2192i 0.0000 + 0.2192i ±0.2192i 0.3713 -0.0000i -0.3713 + 0.0000i ±0.3713 0.5529 + 0.0000i -0.5529 -0.0000i ±0.5529 0.6748 -0.0000i -0.6748 -0.0000i ±0.6748 0.7660 + 0.0000i -0.7660 -0.0000i ±0.7660 0.8368 -0.0000i -0.8368 -0.0000i ±0.8368 0.8919 -0.0000i -0.8591 -0.0000i ±0.8919 0.9341 -0.0000i -0.9342 -0.0000i ±0.9341 0.9652 -0.0000i -0.9652 + 0.0000i ±0.9651 0.9861 + 0.0000i -0.9862 -0.0000i ±0.9860 0.9975 + 0.0000i -0.9976 -0.0000i ±0.9975 0.9996 -0.0000i -0.9996 + 0.0000i ±0.9997 0.9924 + 0.0000i -0.9922 -0.0000i ±0.9929 0.9758 -0.0000i -0.9749 + 0.0000i ±0.9768 0.9478 + 0.0000i -0.9474 + 0.0000i ±0.9509 0.9044 -0.0000i -0.9094 -0.0000i ±0.9144 0.8532 + 0.0000i -0.8919 -0.0000i ±0.8660 0.7613 + 0.0000i -0.7572 -0.0000i ±0.8035 0.6781 -0.0000i -0.6757 + 0.0000i ±0.7233 0.5713 + 0.0000i -0.5711 + 0.0000i ±0.6185

Based on this observation, the expression of r (m) n or r pq (m) can be used as shifts to approximate eigenvalues of the scattering matrix. In fact, if we increase the truncation order M , we can see that the approximations are better for relatively small n and large M [START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF]. The convergence is very quick due to the good approximation. We call this the "early shift". The "early shift" is used in pairs corresponding to the double shift QR algorithm. One pair of the "early shift" will create a bulge to be chased. Wilkinson's shift can be used after the "early shift".

Parallel QR with tightly coupled bulge chasing

The parallel bulge chasing algorithm was proposed by Bai, Demmel [START_REF] Bai | On a block implementation of hessenberg multishift qr iteration[END_REF] and Braman et al. [START_REF] Braman | The multishift QR algorithm. part I: Maintaining well-focused shifts and level 3 performance[END_REF]. In order to benefit from the level 3 BLAS, they parallelize the bulge chasing procedure by performing the chasing of multiple chains of tightly coupled bulges. With the delay and accumulate technique, the main computation work become the matrixmatrix multiplications. The procedure of intrablock chasing and interblock chasing are described below.

To describe the parallel algorithm, we first introduce the data layout mapping in a distributed memory environment as follows:

• The p = p r p c processors are arranged into a p r ⇥ p c grid. Usually the values of p r and p c are set to be as close as possible.

• The N ⇥ N matrix A is partitioned in 2D block cyclic scheme [82] and is mapped on p r ⇥ p c grid as shown in table 5.2. The table shows a 4 grid with p r = p c = 2.

The four processors are denoted as (0, 0), (0, 1), (1, 0), [START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF][START_REF] Martín | Electromagnetic field theory for physicists and engineers: Fundamentals and Applications[END_REF]. The block size is M b ⇥ N b and we require the block to be square M b = N b . Generally, a processor will store a collection of non-contiguous blocks. In table 5.2, if the size of matrix N = 16, then the block size is 4 ⇥ 4 and processor (0, 0) will store the elements A(1 : 4, 1 : 4), A(1 : 4, 9 : 12), A(9 : 12, 1 : 4), A(9 : 12, 9 : 12). An array descriptor stores the details of data layout. The mapping between entries of the global matrix and their corresponding locations in the memory can be established from the array descriptor.

(0,0) (0,1) (0,0) (0,1) (1,0) (1,1) (1,0) (1,1) (0,0) (0,1) (0,0) (0,1) (1,0) (1,1) (1,0) (1,1)
Locally, each processor in the mesh may also utilize multithreading. This can be seen as adding another level of explicit parallelization by organizing all the p = p r ⇥p c processors into a three-dimension mesh.

We use the shifts that are mentioned earlier to introduce the chain of bulges into diagonal blocks. Each of the chains reside on a di↵erent diagonal block. We choose the number of shifts such that each chain covers at most half of the data layout block. The "early shift" is distributed from left-upper diagonal blocks to right-lower diagonal blocks. Each "early shift" is used once. When there are no "early shift" to distribute, we use Wilkinson's shift.

For the intrablock chasing where the chain is chased from the top left corner to the lower right corner within a contiguous diagonal block. We may use a sequence of 3 ⇥ 3

Householder transformations to chase the chain of bulges down some rows to the down right-hand corner of the contiguous diagonal block. We start from the lowest bulge of the block and chase one bulge at a time. The intrablock chasing can be performed locally on the process that own this chain and simultaneously between di↵erent diagonal contiguous blocks which saves computation time. For the interblock chasing where a chain of bulges from one contiguous diagonal block is chased to a di↵erent one on another processor, for each contiguous diagonal blocks in which the bulge chains reside, we create copies of its neighbors and it becomes similar to the case of intrablock chasing. with its neighbor as illustrated in figure 5.2 in the two large gray blocks. The o↵diagonal blocks are then updated by multiplication with the accumulated orthogonal matrix. We perform interblock chasing first for the odd-numbered blocks and then for the even-numbered blocks. This odd-even manner avoids conflicts between di↵erent tightly coupled chains [START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF].

For each diagonal block, the corresponding orthogonal transformations are accumulated into an orthogonal factor. Each orthogonal has the following shape:

U = 0 @ U 11 U 12 U 21 U 22 1 A (5.14)
where U 12 is a lower triangular matrix, and U 21 is a upper triangular matrix. So matrix multiplication by U will be broken into two dense by dense matrix multiplications and two triangular by dense matrix multiplications. Computation time is saved because of the triangular structure.

Parallel AED

The parallel aggressive early deflation (AED) algorithm was proposed in [START_REF] Braman | The multishift QR algorithm. part II: Aggressive early deflation[END_REF]. We divide the Hessenberg matrix H as follows:

H = 0 B B B @ H 11 H 12 H 13 H 21 H 22 H 23 0 H 32 H 33 1 C C C A (5.15)
where H 11 is of size (n k 1) ⇥ (n k 1) and H 33 is of size k ⇥ k. We use the pipeline parallel QR algorithm to find the Schur decomposition of H 33 :

H 33 = V T V ⇤
and perform the following similarity transformation:

0 B B B @ I 0 0 0 1 0 0 0 V 1 C C C A ⇤ 0 B B B @ H 11 H 12 H 13 H 21 H 22 H 23 0 H 32 H 33 1 C C C A 0 B B B @ I 0 0 0 1 0 0 0 V 1 C C C A = 0 B B B @ H 11 H 12 H 13 V H 21 H 22 H 23 V 0 s T 1 C C C A (5.16)
Now the matrix looks like as in figure 5.3. The spike s is denoted as the gray part as in the figure 5.3. It has been proved that it is often the case that the some of the last components of s are very small [START_REF] Braman | The multishift QR algorithm. part II: Aggressive early deflation[END_REF]. If it is the case that the trailing several components of s are negligible, they can set to be zero. The matrix is deflated. This technique often detects convergence much earlier. If it is not the case, we move the eigenvalues to the top left corner of the block H 33 .

Conclusion

In this chapter, we propose the specifically designed parallel QR algorithm for the Cmethod. We present why we propose the "early shift" and how it can be used to accelerate the convergence. We also present the techniques of parallel QR with tightly coupled bulge chasing and parallel AED. These techniques are used to reduce the computational time of the C-method.

In the next chapter, we will use these techniques including "early shift", parallel QR with tightly coupled bulge chasing and parallel AED to analyze gratings, one-dimensional and two-dimensional surfaces. It is a real novelty in the context of the C-method.

Chapter 6

Numerical experiments with parallel QR algorithm -2 Sandy Bridge E5-2670 processors -64 GB memory every node -2 Tesla K20 GPU (Cuda Capability 3.5, 4.8 GB memory every GPU)

• 4 interactive login nodes equipped with:

-2 Sandy Bridge E5-2670 processors -32 GB memory every node

We make use of the 92 nodes and login nodes. We make use of the following libraries:

F CLOADER = $(F C), CCLOADER = $(CC), F CLOADF LAGS = $(F CF LAGS), CCLOADF LAGS = $(CCF LAGS)
Our implementation is based on an simple imitation of the ScaLAPACK routine PDHSEQR we try to use this to a complex implementation with our own shift strategy. We adopt the recommended values such as the size of the deflation window and the tuning parameter NIBBLE which determines when to skip a QR sweep and perform AED in [START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF].

For all the following experiments, we use block factor M b = N b = 50.

Numerical results for one-dimensional case

For the one-dimensional case, we compare the parallel algorithm with the pipeline parallel algorithm which is implemented in ScaLAPACK routine PZLAHQR. The routine PZLAHQR is also compiled and built on the machine P oincare with the same compiler and flags as our program. We also compare the parallel algorithm with or without the "early shift". In the following experiment, we have the physical model parameters (see The accuracy of the early shift can be seen in Table 5.1, page 78, as an example. A plot will be di cult for recognition. The actual eigenvalues that are real are very well approximated, but in the picture, all one can see is blue points on the real axis. We present the Figure 6.2 and 6.3. We also performed an experiment where M and thus the size of the matrix is changed.

The number of cores used here is 16 (i.e. 4 ⇥ 4). Figure 6.6 shows the comparison between the pipeline parallelization and the parallelization with parallel multishift and AED techniques. Figure 6.6 shows that the new version of parallel QR algorithm with parallel multishift and AED is much faster than the existing parallel QR algorithm.

Numerical results for two-dimensional case

For the two-dimensional case, we performed the following experiments. To compare the computation time of sequential algorithm and the specifically parallel algorithm, we present figure 6.8. This figure is based on experimental result of one realisation.

The surface we consider is perfectly conducting. The area of the surface is 64 square wavelengths and = (1) and l x = l y = 1.41 (1) . The incident angle is chosen as ✓ 0 = 30 , ' 0 = 0 . The truncation order is M = 28, so the matrix size is N = 6498. In term of power balance criterion, the truncation order M = 28 gives good enough results (the error is smaller than 1%). [START_REF] Ishimaru | Electromagnetic wave propagation, radiation and scattering[END_REF][START_REF] Ishimaru | Electromagnetic wave propagation, radiation and scattering[END_REF], [START_REF] Petit | Electromagnetic Theory of Gratings[END_REF][START_REF] Petit | Electromagnetic Theory of Gratings[END_REF], [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface[END_REF][START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface[END_REF]. The left most point on the curve for the parallel performance represents that of one core. From this figure, we can see if we have 12800 cores, and we simulate the problem for N R = 200 times. With the naive parallel strategy (we choose N R cores and perform a simulation on each chosen core), it will cost us approximately 5.5 hours. With the new version of parallel strategy (we use 64 cores for one simulation), it will cost us approximately only 8 minutes. This shows when we have many cores, the new version of parallel strategy can be significantly more e cient than the naive parallel strategy and if we only have a few cores, the di↵erence may not so great.

We also compare the computation cost of sequential algorithm and parallel algorithm when the truncation order M is varying. Figure 6.9 shows this comparison based on one realisation. The parameters of surface are the same as above. For the parallel realisation, the number of cores is fixed to be 16, (p r , p c ) = (4, 4). These curves show that the reduction of computation time is important. For instance, the ratio is close to 25 when M = 28. We can see from figure 6.9, the computational time of the sequential code is approximately t = a(2(2M + 1) 2 ) ↵ 1 with ↵ 1 ⇡ 3.2 and the computational time of the parallel code is approximately t = b(2(2M + 1) 2 ) ↵ 2 with ↵ 2 ⇡ 2.2. While these two relations are only approximation from observing the data, they show that how the We compare the results of numerical experiments with experimental data from literature.

Comparison with experimental data for random rough surfaces

We consider an isotropic surface with a = 0.352 (1) and l x = 2.21 (1) . The optical index of the lower medium is ⌫ (2) = 1.62 0.001i. The other simulation parameters are: D = 8 (1) , ✓ 0 = 35 or ✓ 0 = 55 , ' 0 = 0 , M = 28 and the number of realizations is N R = 200. The following figures 6.10-6.13 show results of the implementation of parallel C-methods compared with the experimental data which come from [START_REF] Berginc | Small-slope approximation method: a further study of vector wave scattering from two-dimensional surfaces and comparison with experimental data[END_REF]. In these figures, the DRC which stands for di↵erential reflection coe cient is plotted versus the scattering angle. It is noteworthy that the figure 6. We then present the figures 6.14-6.17 which show results of the implementation of parallel C-methods compared with the experimental data which come from [START_REF] Johnson | Backscattering enhancement of electromagnetic waves from twodimensional perfectly conducting random rough surfaces: a comparison of Monte Carlo simulations with experimental data[END_REF]. In these figures, the bistatic coe cient is plotted versus the observation angle. The perfectly conducting surface under consideration is a very rough surface with a = (1) and l x = 1.41 (1) . The other simulation parameters are: D = 8 (1) , ✓ 0 = 20 , ' 0 = 0 , M = 28 We also consider an anisotropic surface with an = (1) , l x = 2 (1) , l y = 4 (1) . Figure 6.18 gives the co-polarized bistatic coe cient in the incidence plane for a perfectly conducting surface illuminated under ✓ 0 = 20 and ' 0 = 90 . The other simulation parameters are: D = 8 (1) , M = 28 and N R = 200. Figure 6.19 gives the co-polarized return when the incidence angles are ✓ 0 = 20 and ' 0 = 0 . Although the elementary cell area is reduced to 8l x l y , the comparison with experimental data which come from [START_REF] Phu | Copolarized and cross-polarized enhanced backscattering from two-dimensional very rough surfaces at millimeter wave frequencies[END_REF] is satisfactory. The comparison is also conclusive for other polarizations.

Conclusion

In this chapter, from numerical experiments, we have observed that some eigenvalues of the scattering matrix can be approximated e ciently by a certain formula. We designed the "early shift" algorithm to take advantage of this property. We plug this "early shift" method, together with Wilkinson's shift and exceptional shift [START_REF] Demmel | Applied Numerical Linear Algebra[END_REF], into a new parallel QR algorithm. This new QR algorithm uses multiple chains of tightly coupled bulges chasing technique to parallelize the conventional bulge chasing and the aggressive early deflation technique to detect deflation quickly. We apply this specifically designed parallel QR algorithm to the scattering matrix. We also compare the computation time with that of the sequential code. The results show a significant speed up to approximately 40 for 64 cores with our new QR algorithm. This combination of "early shift" and other shifts can also be used in the problems such as linear-quadratic optimal control problem where a large number of eigenvalues and eigenvectors are needed and background of the original problem can provides very good initial approximations.

This parallel QR algorithm can be used for analyzing crossed gratings or two-dimensional random surfaces. Comparisons with experimental data for moderate roughness and isotropic or anisotropic very rough surfaces are conclusive in both co-polarized and cross-polarized components. Comparisons allow the validity of our approach.

In the next chapter, we propose an alternative to the QR algorithm for solving the eigenvalue problem. The proposed method has better scalability than the QR algorithm.

Chapter 7

A proposal: spectral projection method as a global eigensolver

It has been shown that theoretical it is impossible for the standard QR algorithm to be scalable [START_REF] Henry | Parallelizing the qr algorithm for the unsymmetric algebraic eigenvalue problem: Myths and reality[END_REF]. We want a spectral divide-and-conquer algorithms that can provide us all the eigenvalues and eigenvectors which has a very good scalability. Two related work on spectral divide-and-conquer algorithms are [START_REF] Ballard | Minimizing communication for eigenproblems and the singular value decomposition[END_REF] and [START_REF] Demmel | Fast linear algebra is stable[END_REF]. In [START_REF] Ballard | Minimizing communication for eigenproblems and the singular value decomposition[END_REF], the authors present four versions of divide-and-conquer algorithms and present eigenvalue problem that attain lower bounds, and analyze their convergence and communication costs. Paper [START_REF] Demmel | Fast linear algebra is stable[END_REF] shows that all linear algebra operations can also be done stably in O(n !+⌘ ) operations. The authors of [START_REF] Demmel | Fast linear algebra is stable[END_REF] consider known divide-and-conquer algorithms for reducing the complexity of matrix inversion to the complexity of matrix multiplication and show that these algorithm can achieve the same forward error bound (bound on the norm of the error in the output) as a conventional backward stable algorithm. We consider a spectral divide-and-conquer algorithm which essentially transform the eigenvalue to linear system problems.

In this chapter, we propose an alternative to QR algorithm for solving the eigenvalue problem. We propose a global eigensolver by a combination of the Sakuria and Sugiura method (SSM) and multiple implicitly restarted Arnolid method with nested subspaces (MIRAMns). The first method allows the computation of interior eigenvalues while the second permits to compute the eigenvalues in the extremities of spectrum. This proposed global eigensolver allows us to calculate all or a large number of the eigenvalues of a generalized matrix. The MIRAMns [START_REF] Shahzadeh Fazeli | A key to choose subspace size in implicitly restarted arnoldi method[END_REF] is a variant of the IRAM [START_REF] Sorensen | Implicitly restarted Arnoldi/Lanczos methods for large scale eigenvalue calculations[END_REF] that is based on the projection of the eigenproblem on several nested subspaces instead of a single one.

It can therefore use the eigen-information of interest obtained in all subspaces to update the restarting vector. We will take the real matrix as an example, but this proposed global eigensolver can also be applied to a complex matrix in a similar way. Given a matrix A 2 R n⇥n , MIRAMns uses the Arnoldi method to compute the Ritz elements

of A in t nested Krylov subspaces K m i ,v , 1  i  t with K m i ,v ⇢ K m i+1 ,v
. MIRAMns selects the "best subspace" by finding which subspace contains the "best" current Ritz elements. We denote the size of the "best" subspace m best . The next step of MIRAMns is to apply the shifted QR procedure to the m best ⇥ m best matrix. By choosing the undesired eigenvalues as shifts, the information related to the desired eigenvalues are concentrated in the leading submatrix. MIRAMns then completes Arnoldi projection of t nested Krylov subspaces by restarting with this submatrix whose size is the number of wanted eigenvalues. The MIRAMns can be used to provide the extremes of the spectrum of A with good convergence properties.

The rest of the spectrum of A is thus located in a finite domain D described by the extremes provided by MIRAMns. The domain D is then divided into several sub-domains.

For each sub-domain, the contour integral based projection method projects the matrix pencil (A I) onto the subspace associated with eigenvalues that are located in the sub-domain via numerical integration. A moment-based approach can be used to find the eigenvalues in each sub-domain independently. To avoid the numerically unstable problem of the computation using explicit moments, we often use a Rayleigh-Ritz procedure instead. For the computation of the contour integral, we solve a certain number of linear systems derived from the matrices A and I. When A is large, the computational costs for solving linear systems are dominant.

Algorithms

This spectral projection method as a global eigensolver is a combination of MIRAMns and SS method. We present these methods and how they can be combined to form a global solver in this section.

MIRAMns

The MIRAMns is a variant of the implicitly restarted Arnoldi method (IRAM). Recall that IRAM allows us to compute a few eigenvalues in the extremes of the spectrum of a large matrix. For that, IRAM combines the implicitly shifted QR algorithm with a k-step Arnoldi factorization to obtain a truncated form of implicitly shifted QR iteration. This approach o↵ers a more e cient and numerically stable formulation than explicitly restarted Arnoldi method (ERAM) [START_REF] Saad | Variations on Arnoldi's method for computing eigenelements of large unsymmetric matrices[END_REF]. By using IRAM instead of ERAM, the numerical di culties and storage problems normally associated with the Arnoldi process are avoided. The algorithm is capable of computing a few (k) eigenvalues with user specified features such as largest real part, largest magnitude, smallest real part or smallest magnitude.

We start from the Arnoldi factorization of length m = p + k,

AV m = V m H m + f m e T m (7.1)
where e T m means the transpose of the vector e m . We apply p shifts µ 1 , ..., µ p implicitly

AV + m = V + m H + m + f m e T m Q (7.2)
where

Q = Q 1 Q 2 • • • Q p the product of the orthogonal matrices related to µ 1 , ..., µ p and V + m = V m Q, H + m = Q T H m Q.
From the fact that Q is the product of p (unitary) Hessenberg matrices, it is easy to see that Q has p non-zero o↵-diagonals below its main diagonal. So the first k 1 elements of the vector e T m Q are zeros. So if we discard the last p columns of equation ( 7.2), we will have

AV + m (:, 1 : k) = V + m (:, 1 : k + 1)H + m (:, 1 : k) + f m e T m Q(:, 1 : k) = V + m (:, 1 : k)H + m (:, 1 : k) + h + k+1,k v + k+1 e T k + q + m,k f m e T k = V + m (:, 1 : k)H + m (:, 1 : k) + (v + k+1 ˆ k + f m k )e T k (7.3)
where we denote

ˆ k = h + k+1,k , k = q + m,k . Equation (7.
3) can be also written as

AV + k = V + k H + k + f + k e T k (7.4)
A description of this IRAM algorithm can be found in Algorithm 3.

Algorithm 3

The implicitly restarted Arnoldi process Input: (A, V, k, m) with AV m = V m H m + f m e T m , an m-step Arnoldi factorization, with m = p + k Output: k eigenvalues with user specified features and their corresponding eigenvectors. Compute the spectrum of H m : (H m ), if convergence, stop. Otherwise, select set of p shifts µ 1 , µ 2 , ..., µ p ; 3:

q T = e T m ;
4:

for j = 1, 2, ..., p do 5:

Factor [Q j , R j ] = qr(H m µ j I);

6:

H m = Q T j H m Q j , V m = V m Q j , q T = q T Q j ; 7:
end for 8:

f k = v k+1 H m (k + 1, k) + f m q T (k); V k = V m(1:n,1:k) ; H k = H m(1:k,1:k) ; 9:
Begining with the k-step Arnoldi factorization,

AV k = V k H k + f k e T k ,
AV m = V m H m + f m e T m ;
10: end for

The MIRAMns takes advantage of IRAM by choosing a set of initial Krylov subspaces that di↵er only by their sizes. MIRAMns chooses a set of t di↵erent subspace sizes M = (m 1 , ..., m t ) with a strict order m 1 < ... < m t . MIRAMns then performs t Arnoldi projections on the subspaces K m i ,v , for 1  i  t, with K m 1 ,v ⇢ K mt,v and initial vector v. MIRAMns then chooses the subspace size m best by finding the Arnoldi factorization which o↵ers the best Ritz estimation for k desired eigenpairs.

AV m best = V m best H m best + f m best e T m best (7.5)
From this k-step Arnoldi factorization, p i = m i k, 1  i  t additional steps of Arnoldi factorizations are applied to obtain t new projections onto the updated subspaces. This procedure can go on until convergence.

A description of this MIRAMns algorithm is presented in algorithm 4.

In order to select the best results in the above algorithm, we consider that (

V m i , H m i , f m i ) is "better" than (V m j , H m j , f m j ) if r m i k < r m j k where r m k = max(⇢ 1,m , ..., ⇢ k,m ) is defined by Ritz estimates ⇢ i,m = | m e T m y (m) i |.
One advantage of the MIRAMns is that it overcomes the problem of sensitivity of convergence with respect to small perturbation of the subspace size that occurs in the normal restarted Arnoldi methods. It achieves this by choosing the "best" size among Algorithm 4 Multiple IRAM with nested subspaces Input: (A,

V m i , H m i , f m i ) with AV m i = V m i H m i + f m i e T m
i , an m i -step Arnoldi factorization, with m i = p i + k, where 1  i  t. Output: k eigenvalues with user specified features and their corresponding eigenvectors.

1: for l = 1, 2, 3, ... until convergence do 2:
Compute the spectrum of H m i : (H m i ), if convergence, stop. Otherwise compute their associated eigenvectors and residuals for 1  i  t.

3:

Select the best results in these subspaces and the associated best subspace size

m best . Set m = m best , H m = H m best , V m = V m best , f m = f m best . 4:
Select a set of p = m k shifts (µ

(m) 1 , ..., µ (m)
p ) based on (H m ) or other information.

5:

q T = e T m ;

6:

for j = 1, 2, ..., p do 7:

Factor [Q j , R j ] = qr(H m µ (m) j I);
8:

H m = Q H j H m Q j , V m = V m Q j , q = q H Q j ; 9:
end for 10:

f k = v k+1 ˆ k + f m k ; V k = V m(1:n,1:k) ; H k = H m(1:k,1:k) ; 11:
Begining with the k-step Arnoldi factorization,

AV k = V k H k + f k e T k , apply p i = m i k additional steps of the Arnoldi process to obtain t new m i -step Arnoldi factorization, AV m i = V m i H m i + f m i e T

SS method

The SS method was introduced in [71], [START_REF] Sakurai | CIRR: a Rayleigh-Ritz type method with contour integral generalized eigenvalue problems[END_REF]. Given a finite domain D, we want to calculate the eigenvalues of A that lie in it. Suppose that we cover the domain D with s subdomains D i , 1  i  s, such that D = [ i Di . Now we only need to calculate the eigenvalues that lie in each sub-domain D i , 1  i  s and these tasks can be performed in parallel.

For each subdomain D

i , 1  i  s, define f (z) = u H (zI A) 1 v (7.6)
with non-zero vectors u, v 2 R n . Define

µ k = 1 2⇡i Z @D i (z z 0 ) k f (z)dz, k = 0, 1, ... (7.7)
where z 0 is located inside D i . Suppose there are m eigenvalues lie in D i , then these eigenvalues are exactly the same as the eigenvalues of pencil

H < m H m , with the m ⇥ m Hankel matrices H m = [µ i+j 2 ] and H < m = [µ i+j 1 ], 1  i, j  m.
See [START_REF] Sakuria | A projection method for generalized eigenvalue problem using numerical integration[END_REF] for a proof, here u, v are any non-zero vectors.

The number m can be calculated from the following formula [START_REF] Senzaki | An estimation method of eigenvalues distribution with substructuring[END_REF]:

m = 1 2⇡i Z @D i tr(F (z) 1 )dz (7.8)
For the eigenvectors, if we define

s k = 1 2⇡i Z @D i (z z 0 ) k (zI A) 1 vdz, k = 0, 1, ... (7.9)
and V m is the Vandermonde matrix

V m = 0 B B B B B B @ 1 1 • • • 1 1 z 0 2 z 0 • • • m z 0 . . . . . . . . . ( 1 z 0 ) m 1 ( 2 z 0 ) m 1 • • • ( m z 0 ) m 1 1 C C C C C C A (7.10)
then the associated eigenvectors are given by the formula:

[q 1 , ..., q m ] = [s 0 , ..., s m 1 ]V T m (7.11)
Suppose that D i is a circle centered at point with radius ⇢, a description of the idea of contour integral based projection (explicit moments) can be found in algorithm 5.

Algorithm 5

The contour integral based projection method (explicit moments) Input: (u, v 2 R n , N, m, , ⇢) Output: approximated eigenvalues of A that lie in the finite domain D i : ˆ 1 , ..., ˆ m and their associated eigenvectors: q1 , ..., qm

1: Set ! j = + ⇢exp(2⇡ p 1j/N ), j = 0, ..., N 1; 2: Form y j = (! j I A) 1 v, j = 0, ..., N 1; 3: Set f j = u H y j , j = 0, ..., N 1; 4: Compute μk = 1 N P N 1 j=0 (! j ) k+1 f (! j ), k = 0, ..., 2m 1; 5
: Compute ŝk = 1 N P N 1 j=0 (! j ) k+1 y j , k = 0, ..., m 1
6: Compute the eigenvalues ⇣ 1 , ..., ⇣ m of the pencil H < m H m ;

7: Compute q1 , ..., qm given by [q 1 , ..., qm ] = [ŝ 0 , ..., ŝm 1 ] V T m ;

8: Set ˆ j = + ⇣ j , j = 1, ..., m.

IF some eigenvalues are very close to each in the contour, the Hankel matrices H < m and H m are very ill-conditioned. The Rayleigh-Ritz method can be used to avoid the explicit use of moments and improve numerical accuracy. We apply a Rayleigh-Ritz procedure by projecting the matrix A to à = ⇧ T A⇧ with an unitary basis ⇧ 2 C n⇥m . The eigenvalues of A can be approximated by the Ritz values of the projected pencil ( Ã, I).

In practice, the numerical value of m from equation (7.8) is not always an integer, so it is more convenient and more e cient to choose a number M ( m) as the size of Hankel matrices . This choice can decrease the influence of the quadrature error su↵ered from eigenvalues located outside the boundary.

A description of the contour integral based projection (Rayleigh-Ritz) can be found in algorithm 6.

Algorithm 6

The contour integral based projection method (Rayleigh-Ritz) Input: (v 2 R n , N, M, , ⇢) Output: approximated eigenvalues of A that lie in the finite domain D i : ˆ 1 , ..., ˆ m and their associated eigenvectors: x1 , ..., xm 1: Set ! j = + ⇢exp(2⇡ p 1(j + 1/2)/N ), j = 0, ..., N 1;

2: Solve (! j I A)y j = v, for y j , j = 0, ..., N 1;

3: Compute ŝk = 1 N P N 1 j=0 (! j
) k+1 y j , k = 0, ..., M 1;

4: Compute construct an unitary basis ⇧ from (ŝ 0 , ..., ŝM 1 ); 5: Form à = ⇧ T A⇧; 6: Compute eigenpairs (✓ j , w j ) with j = 1, ..., M of ( Ã, I); 7: Set p j = ⇧w j , j = 1, ..., M; 8: Select the approximated eigenpairs ( ˆ 1 , x1 ), ..., ( ˆ m , xm ) from (✓ j , p j ), j = 1, ..., M;

A global eigensolver

For the given matrix A, we first apply the MIRAMns algorithm to calculate k s eigenvalues that have the smallest magnitude and k l eigenvalues that have the largest magnitude.

Thus we get the extremes of the spectrum of A:

| 1 |  ...  | ks |  | n+1 k l |  ...  | n |. Denote R = | n+1 k l |, r = | ks |,

Parallelism analysis

The proposed method has a good parallel nature. For the MIRAMns, the computation in di↵erent subspace of multiple implicitly restarted Arnoldi method (MIRAM) can be performed in parallel. In addition to this coarse grain parallelism, the communication of the eigen-information of interest among processes can also be made asynchronously. For the SS method, the solving of each sub-problem associated with the sub-domain can be done independently. Besides, for each sub-problem, one needs to use the trapezoidal rule to calculate the contour integral. The calculation related to each quadrature points can be done in parallel. Finally, for each quadrature point, we need to solve the associated linear system, this can also be done using parallel linear solvers.

Numerical experiments

We test the matrix

A n with n = 1000, A n (i, i) = 2, A n (i, i + 1) = A n (i + 1, i) = 1, the
other entries are 0.

A n = 0 B B B B B B B B B B B B B B B B @ 2 1 1 2 1 1 2 1 . . . . . . . . . 1 2 1 1 2 1 1 2 1 C C C C C C C C C C C C C C C C A (7.
12)

The eigenvalues of A n are proved to be

k (A n ) = 2 + 2 cos( k⇡ n + 1 ), k = 1, 2, ..., n (7.13)
We choose this matrix because the eigenvalues are known and all of them are real, so we can check our algorithm easily with this matrix. as [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] and we denote this as M IRAM ns [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF]. We observed that MIRAMns converges with fewer iterations than IRAM. For the testing matrix A n , the 3 eigenvalues with largest magnitude calculated from M IRAM ns [START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] The horizontal line represents the x-axis and the vertical line represents the y-axis. The red points represent the quadrature points and the blue points represent the eigenvalues that are calculated in the domain described by the red points. The total number of return the value in 0.020527 second, which is more than 100 faster than our method with 20 cores. However, the idea is that our method may still work when some existing It can be noticed that when some of the eigenvalues are very close, the precision will decrease. To increase the precision, one can decrease the radius of the circle, thus, there will be a demand of more circles to be addressed. Similar to the test matrix A n , we also run the parallel version of this global eigensolver for the test matrix B n . Figure 7.20 shows the scalability for the matrix B n wiht n = 1000 when we increase the number of cores. We can observed that the line is not as "straight"

as that for the matrix A n . The slope of the approximated line is smaller than that for the matrix A n . This decrease of performance can be explained if one observes the location of the spectra of the two matrices. For matrix B n , there are some eigenvalues More precise error analysis of this method and the optimal division strategy, especially when the whole domain D is of two dimension, is a part of our future work. We may want to extend this method to apply it to the C-method in the future.

In the next chapter, we will propose a new approach of the curvilinear coordinate method where we do not use the translation coordinate system. The proposed method allows analyzing the complex phenomenon of incident energy absorption. The new version of the C-method could be an attractive alternative to analyze multilayered grating having parallel or non-parallel interfaces. The C-method as an initial value problem

The C-method in the previous chapter is not very e cient when we are dealing with multi-layer gratings. We want to find other solutions. Especially, we want to explore the potential parallelization of multi-layer gratings. That is we want to find a way to deal with each layer independently and then combine them. In this chapter we propose a new version of the C-method.

Eigenvalue problem and initial value problem

We will propose a new approach of the curvilinear coordinate method where we don't use the translation coordinate system. We consider two horizontal plane above and below the grating. We define a coordinate system that the grating surface and both horizontal planes correspond to surface coordinate. Similar coordinate systems have been defined

for analyzing discontinuities in rectangular waveguides [START_REF] Dusséaux | Etude de transformateur plan-E dans un système de coordonnées non orthogonales[END_REF][START_REF] Dusséaux | Analysis of rectangular waveguide Hplane junctions in nonorthogonal coordinate system[END_REF][START_REF] Dusséaux | Analyse de composants plan-E symétriques en guides d'onde à section rectangulaire[END_REF][START_REF] Dusséaux | Telegraphist's equations for rectangular waveguides and analysis in nonorthogonal coordinates[END_REF] and radiation loss of optical waveguides [START_REF] Afifi | Statistical study of radiation loss from planar optical waveguides: The curvilinear coordinate method and the small pertubation method[END_REF]. Inside the area A delimited by the two horizontal planes, in the air and the low medium, the covariant formalism of Maxwell's equations lead to a di↵erential equation system with non-constant coe cients. This system represents an initial value problem. The curvilinear coordinate method expressed in the translation coordinate system leads to an eigenvalue problem. It is the fundamental di↵erence with this new approach. The scattering matrix (S matrix) relates the amplitudes of outgoing planes waves to those of incoming waves. We show how to determine the S matrix by solving the initial value problem, by satisfying the boundary conditions on the grating interface and using the continuity relations on the two horizontal planes between covariant components of fields and Cartesian ones.

8.2 From Maxwell's equations in covariant form to an ini-

tial value problem

For simplicity, we consider only the one-dimensional case. In the Cartesian referential

Oxyz, the grating is represented by a periodic cylindrical surface y = a(x) (figure 8.1).

This surface separates the air (medium 1) from the medium with a real or complex refractive index (medium 2). The grating of period D is illuminated by a monochromatic plane wave under the incidence ✓ 0 . The incident wave vector lies in the xOy plane. For E // polarization, the electric vector is parallel to the grooves. For H // polarization, it is the case of magnetic vector. The letter m denotes indi↵erently the upper medium (m = 1) or the lower medium (m = 2). Henceforth, ⌫ (m) , Z (m) and k (m) indicate the optical index. the inpendance and the wave number of medium (m).

As show by 8.1, the space is divided into four regions, Within the regions y y 1 and y  y 2 , we consider the Cartesian coordinates (x, y, z). Outside the grooves, i.e.

when y > max(a(x)) and y < min(a(x)), the di↵racted field can be represented by a combination of elementary plane waves, the Rayleigh expansion (3.7).

Within the regions A 1 and A 2 defined by a(x)  y  y 1 and y 2  y  a(x), we consider the non-orthogonal coordinate system defined as follows:

8 > > > > > < > > > > > : x 0 = x u = y m y a(x) ym a(x) z 0 = z (8.1)
The grating surface y = a ( x) coincides with the coordinate surface u = 0 and the horizontal plane y = y m with u = y m . The problem consists in determining the S matrix by solving Maxwell's equations under covariant form expressed in the coordinate system (3.9) and by using continuity relations in planes y = y m between covariant components of fields and Cartesian ones.

The covariant components (v x 0 , v u , v z 0 ) of a vector v are obtained from the Cartesian coordinate (v x , v y , v z ) as follows:

0 B B B @ v x 0 v u v z 0 1 C C C A = A 0 B B B @ v x v y v z 1 C C C A (8.2)
Here A is the transformation matrix:

A = A i i 0 = 0 B B B @ @x @x 0 @y @x 0 @z @x 0 @x @u @y @u @z @u @x @z 0 @y @z 0 @z @z 0 1 C C C A = 0 B B B @ 1 (y m u) ȧ(x) ym 0 0 ym a(x) ym 0 0 0 1 1 C C C A (8.3)
where ȧ(x) = da(x) dx .

From the above equation, we can make several observations:

• the covariant component v z 0 (x, u, z) is equal to the Cartesian component v z (x, y, z)
and it is parallel to the surface u = 0 and u = y m .

• The covariant component v u (x, y, z) is proportional to the Cartesian component v y (x, y, z).

• The covariant component v x 0 (x, u, z) is tangential to the grating interface u = 0 and can be identified with the Cartesian component v x (x, y, z) on the plane u = y m .

The passage of covariant components (v

x 0 , v u , v z 0 ) to the contravariant components (v x 0 , v u , v z 0 ) is obtained by the metric tensor G 0 B B B @ v x 0 v u v z 0 1 C C C A = G 1 0 B B B @ v x 0 v u v z 0 1 C C C A (8.4)
The tensor G is defined by (see equation 4.6):

G = g i 0 j 0 = AG c A t = 0 B B B @ g x 0 x 0 g x 0 u g x 0 z 0 g ux 0 g uu g uz 0 g z 0 x 0 g z 0 u g z 0 z 0 1 C C C A (8.5)
where G c is the Cartesian system that is equal to the identity matrix and A t denotes the transpose matrix of A.

So, we find:

G = 0 B B B @ 1 + (ym u) 2 y 2 m ȧ2 (x) (ym a(x))(ym u) y 2 m ȧ(x) 0 (ym a(x))(ym u) y 2 m ȧ(x) (ym a(x)) 2 y 2 m 0 0 0 1 1 C C C A (8.6)
and

G 1 = g i 0 j 0 = 0 B B B @ g x 0 x 0 g x 0 u g x 0 z 0 g ux 0 g uu g uz 0 g z 0 x 0 g z 0 u g z 0 z 0 1 C C C A = 0 B B B @ 1 u ym ym a(x) ȧ(x) 0 u ym ym a(x) ȧ(x) (u ym) 2 ȧ2 (x)+y 2 m (ym a(x)) 2 0 0 0 1 1 C C C A (8.7) g = det(G) = (y m a(x)) 2 y 2 m (8.8)
From the Maxwell's equations written in the coordinate system, we find:

8 > > > > > < > > > > > : @F (m) (x,u) @u = jk (m) ym a(x) ym G (m) (x, u) + jk (m) ym u ym ȧ(x)G (m) u (x, u) @G (m) (x,u) @u = @G (m) u (x,u) @x jk (m) ym a(x) ym F (m) (x, u) @F (m) (x,u) @x = jk (m) ym u ym ȧ(x)G (m) (x, u) + j ymk (m) ym a(x) (1 + (ym u) 2 y 2 m ȧ2 (x))G (m) u (x, u) (8.9) In E // polarization, F (m) = E (m) z 0 , G (m) = Z (m) H (m) x 0 and G (m) u = Z (m) H (m) u . In H // polarization, F (m) = Z (m) H (m) z 0 , G (m) = E (m) x 0 and G (m) u = E (m)
u . The covariant components F (m) (x, u) and G (m) (x, u) are tangential to the grating interface and they appear in the boundary conditions at u = 0.

The periodic functions a(x) and ȧ(x) are expanded in Fourier series. The periodicity with respect to the variable x, as well as the excitation by a plane wave, leads to an expansion of founctions

F (m) (x, u), G (m) (x, u) and G (m) u (x, u) in terms of the quasi- periodic functions exp( j↵ n x). 8 > > > > > < > > > > > : F (m) (x, u) = P +1 n= 1 f (m) n (u)exp( j↵ n x) G (m) (x, u) = P +1 n= 1 g (m) n (u)exp( j↵ n x) G (m) u (x, u) = P +1 n= 1 g (m) u,n (u)exp( j↵ n x) (8.10)
Substituting these expansions into (8.9) and projecting on basis functions exp( j↵ n x) leads to a set of partial di↵erential equations relating f n (u):

8 > < > : d f (m) (u) du = j ym u ym D(u) Ȧ↵ f (m) (u) jk (m) D(u)g (m) (u) dg (m) (u) du = j k (m) (↵D(u)↵ k (m)2 B) f (m) (u) j ym u ym ↵D(u) Ȧg (m) (u) (8.11)
where

B = I A/y m , D = B(I + (y m u) 2 y 2 m Ȧ Ȧ) 1 (8.12)
Vector f (m) and g(m) contains the coe cients f n (u) respectively. ↵ is a diagonal matrix with the propagation coe cients ↵ n along the diagonal and I is the identity matrix. A is the Toeplitz matrix generated by the Fourier coe cients a n of functions a(x), such that its (p, q) element is a p q . Ȧ is the Toeplitz matrix generated by the Fourier ȧn of the profile derivative. The numerical solution of system (8.11) requires a truncation order M . Then, the covariant components F (m) (x, u) and G (m) (x, u) are described by only 2M + 1 expansion coe cients f

) and 2M + 1 incoming waves

(amplitude c (1 ) n , c (2+) n
), see figure 3.2.

Numerical Implementations

The di↵erential system (8.11) has non-constant coe cients and represents an initial value problem. We propose a procedure for obtaining the N dimensional S matrix (N = 4M +2). First, we define N independent vectors satisfying the boundary conditions on the grating interface u = 0. In the E // polarization, the continuity relations on the electric and magnetic components are given by: 8 > < > :

F (1) (x, u = 0) = F (2) (x, u = 0) ⌫ (1) G (1) (x, u = 0) = ⌫ (2) G (2) (x, u = 0) (8.13)
In H // polarization, we have:

8 > < > : ⌫ (1) F (1) (x, u = 0) = ⌫ (2) F (2) (x, u = 0) G (1) (x, u = 0) = G (2) (x, u = 0) ( 8
n (u = 0) = f (1)
n (u = 0)

⌫ (1) g (1)
n

(u = 0) = ⌫ (2) g (2)
n (u = 0) (8.15)

Similarly, in H // polarization, we get:

8 > < > : ⌫ (1) f (1) n (u = 0) = ⌫ (2) f (2)
n (u = 0) g

n (u = 0) = g n (u = 0). In the E // polarization, they are contained in the following matrix:

0 @ F (m) (u = 0) G (m) (u = 0) 1 A = 0 @ I I I/⌫ (m) I/⌫ (m) 1 A (8.17)
I is the (2M + 1) dimensional identity matrix. In the H // polarization, we use:

0 @ F (m) (u = 0) G (m) (u = 0) 1 A = 0 @ I/⌫ (m) I/⌫ (m) I I 1 A (8.18)
For a perfectly conducting grating, the fields inside the conductor vanish, see equation and in H // polarization, 0 @ F (m) (u = 0)

G (m) (u = 0) 1 A = 0 @ I 0 1 A (8.20)
The coupled di↵erential equation system (8.11) is solved for each column vector of initial condition matrices. For each medium, this step requires numerical integrations with

f (m) n (u = y m ) = c (m+) n exp( j (m) n y m ) + c (m ) n exp(+j (m) n y m ) g (m) n (u = y m ) = c (m+) n (m) n k (m) exp( j (m) n y m ) c (m ) n (m) n k (m) exp(+j (m) n y m ) (8.21)
For each medium, we deduce from (8.21) the amplitudes c (m±) n

. For a dielectric grating, the N dimensional S matrix is obtained as follows:

S = 0 @ C (1+) C (2 ) 1 A 0 @ C (1 ) C (2+) 1 A 1 (8.22)
For a perfectly conducting grating, we have:

S = C (1+) (C (1 ) ) 1 (8.23)

On the computational time

The used iterative algorithm is a variable order Adams-Bashforth-Moulton PECE solver (Prediction/Evaluation/Correction/Evaluation). It is a multistep solver and needs the solutions at several preceding spatial points to compute the current solution. Results presented in the next section are provided using MATLAB and the solver ODE113.

The dominant computational cost of the proposed method is due to numerical integrations and depends on the relative and absolute tolerances used by the algorithm. This relative tolerance controls the number of correct digits in all solution components, except those smaller than the absolute tolerance thresholds. The absolute tolerance is a threshold below which the value of the i th solution component is unimportant. The absolute error tolerances determine the accuracy when the solution approaches zero. For a given grating with given tolerances, the computational cost is O(N 3 ).

For a multilayer grating (n + 1 layer labeled with 1, 2, ..., n, n + 1 in sequence), we will show that it is possible to form the local scattering matrix S i,i+1 and then glue them to form the global matrix S 1,n+1 . We can also explore the parallelism in this gluing operation. For example, we can in the first setp, obtain S 1,3 , S 3,5 , S 5,7 ... in parallel, and in the second step, obtain S 1,5 , S 5,9 ,... in parallel, and so on. This needs only O(log(n)) steps to get the global scattering matrix. Moreover, for the calculation of the local scattering matrix S i,i+1 , the proposed method leads to systems of first-order linear di↵erential equations, the solution of which requires the choice of an iterative algorithm. The proposed method is based on numerical integrations with N independent initial vectors. The kernel of this computing process is the numerical integration. This approach is also particularly well adapted to large-scale parallel and distributed architectures. Indeed, in the context of a distributed system comprising a network of machines, each of the problems could be solved on a machine whose architecture can be single or multiple processors. The proposed new method has a significant degree of coarse grain parallelism and requires little communication. These features o↵er the possibility of reducing dramatically the computation time. It's an advantage compared with the conventional C-method which leads to eigenvalue problems.

Numerical results

We consider a perfectly conducting sinusoidal grating defined by a(x) = h 0 cos(2⇡x/D) with D = . Under the incidence angle ✓ 0 = 30 , the grating is in first order Littrow mounting and ✓ 1 = ✓ 0 , see paragraph 3.1.2. As a result, the grating can present a perfect blazing in the minus-first-order. Figure 8.2 gives the e ciency curves obtained with the reference C-method (C1) and based on Fourier series factorization rules and the curves derived from the new version of the C-method (C2). For the reference method, the truncation order M is 18. This value provides a very good accuracy on the e ciencies. Figure 8.2 shows superimposed curves.

Comparisons are conclusive and validate the proposed method. As shown in Figure 8.3, the new version of the C-method used with M = 9 checks the power balance with an error smaller than 10 3 whatever the groove depth. We obtain similar results in polarization E // . The new approach is well-adapted to analyze this sinusoidal grating with the peak-peak amplitude smaller or equal to two wavelengths. The proposed method only uses Rayleigh expansions outside the grooves and does not use the Rayleigh hypothesis stipulating that the scattered field away from the surface can be extended down onto the grating even though it is formed by solely up-going waves. The theoretical validity of the Rayleigh hypothesis has given rise to some works for di↵raction gratings. A classical result can be mentioned: for a perfectly conducting grating defined by a(x) = h cos(2⇡x/D) with the Dirichlet condition,the assumption does not hold if h/D > 7% [START_REF] Petit | Sur la di↵raction d'une onde plane par un rseau infiniment conducteur[END_REF]. The proposed method gives the e ciencies with a good accuracy with h/D = 1. depth and it is smaller than 10 3 whatever the profile amplitude. We obtain similar results in polarization E // .The new approach of the C-method is well-adapted to analyze this lossless dielectric grating when the peak-peak amplitude is smaller or equal to two wavelengths. Figure 8.5 also shows the error on the power balance for the conventional C-method used with M = 9. The error given by C1 is lowest when h 0 < 3 /10. On the range 3 /10 < h 0 < 3 /5, C1 and C2 methods are equivalent in terms of accuracy on the power balance. When 3 /5 < h 0 < , the proposed method leads to the lowest errors. The proposed method allows analyzing the metallic grating under consideration and the phenomenon of quasi-total absorption by surface plasmons in H // polarization.

Application to multilayer gratings with homogeneous medium

In this section, we extend the C-method as an initial value problem to multilayer gratings with homogeneous medium. We consider a n + 1 layer di↵raction grating, thus there are n interfaces separating the layers. From the uppermost to downmost, these layers are composed of medium 1 to medium n + 1. Each medium has a constant optical index.

Thus we can calculate the scattering matrix S i,i+1 which associate the incoming and outcoming waves from medium i and i + 1. Then we collect all the scattering matrix of adjacent medium, S i,i+1 , i = 1, ..., n and obtain the global matrix S 1,n+1 by combination of elementary matrices S i,i+1 , i = 1, ..., n.

We consider a n + 1 layer di↵raction grating. The interface is represented by a periodic cylindrical surface y = a i (x), 1  i  n. This surface separates the medium i from the medium i + 1. In Figure 8.8, two representative adjacent interfaces are shown. We consider separable layered grating, meaning that there exists a horizontal line y = y i+1 separates the interface y = a i (x) and the interface y = a i+1 (x) for each i. The interfaces in general have di↵erent functional forms and amplitudes, but there exists a value D such that D is the period of the function or a multiple of the period. The thickness h i is measured between the middle lines of the two boundaries. The medium between the interface y = a i (x) and the interface y = a i+1 (x) is homogeneous with optical index ⌫ (i) , impedance Z (i) and wave number k (i) .

As the interfaces are separable, we can consider each interface separately and then combine to form the global matrix. For each interface y = a i (x), we consider the problem as in the previous section. Then, we reduced the problem to the previous one.

We define the local scattering matrix(S-matrix) to relate the amplitudes of outgoing plane waves (c

(i+) n , c ((i+1) ) n
) to those of incoming waves (c

(i )
n , c

((i+1)+) n
) such that:

0 @ c (i+)
c ((i+1) ) We then combine the local scattering matrices to form the global scattering matrix S 1,n+1 . In fact, if we have local scattering matrices S p,q and S q,r such that p < q < r, then we can obtain the scattering matrix S p,r . Eliminating the vectors c (q+) and c (q ) , one glue the two scattering matrices to be one S p,r : S p,r = 0 @ S (+ ) p,q + S (++) p,q (I S (+ ) q,r S ( +) p,q ) 1 S (++) q,r S ( ) p,q S (++) p,q (I S (+ ) q,r S ( +) p,q ) 1 S (++) q,r S ( ) q,r (I S ( +) p,q S (+ ) q,r ) 1 S ( ) p,q S ( +) q,r + S ( ) q,r (I S ( +) p,q S (+ ) q,r ) 1 S ( +) p,q S (++) q,r and so on until we get the global matrix S 1,n+1 .

We perform an experiment in the article [START_REF] Popov | Theoretical study of the anomalies of coated dielectric gratings[END_REF]. The grating considered is sinusoidal and defined by a 1 (x) = a 1 cos( 2⇡ D x), a 2 (x) = a 2 cos( 2⇡ D x). Figure 8.9 shows the e ciency curves under the polarization E // , with the value of sin ✓ 0 varying from 0.24 to 0.38.

With the same parameters as in that paper, figure 8.9 is exactly the same as in the paper [START_REF] Popov | Theoretical study of the anomalies of coated dielectric gratings[END_REF]. We can also analyze structure with non parallel interfaces. We perform more experiments to see how this figure changes when we change the amplitude of function We see from Figure 8.9 that the zeroth-order e ciency changes from 0 to 1. In Figure 8.10, the e ciency can not reach 0, it changes from a small positive value to 1. In figure 8.11, the e ciency can neither reach 0 nor 1. It can be observed that when a 1 varies from 0.03µm to 0.01µm, the jump becomes steeper. One can also observe that the place of the jump moves towards the left direction. When we vary a 2 from 0.03µm to 0.01µm, no similar phenomena can be observed. In fact, the curves seem almost stay the same.

Conclusion

In this chapter, we studied our approach of the C-method as an initial value problem for the e cient calculation of the N-dimensional scattering matrix of a grating. We have shown that this formulation is an interesting tool for analyzing perfectly conducting or dielectric gratings with deep grooves. The proposed method allows analyzing the complex phenomenon of incident energy absorption. We then extend this method to multilayer with homogeneous medium. We applied this method to multilayer gratings with an arbitrary number of interfaces. We have shown how to combine the local scattering matrix to obtain the global one. We have validate our method by comparison our results with that from published paper. The proposed method has very good accuracy as well as a natural of two level parallel property. This new version of C-method is an attractive alternative to analyze multilayered grating having parallel or non-parallel interfaces. We are currently working on the extension of this method to apply it to multilayer with inhomogeneous medium and non-parallel interfaces.

Chapter 9

Conclusion

In this thesis, we study the electromagnetic di↵raction by gratings and random rough surfaces. The C-method is an exact method based on Maxwells equations under covariant form written in a nonorthogonal coordinate system. The C-method leads to the eigenvalue problem of the high dimension, dense and non-symmetric scattering matrix.

All the eigenvalues and eigenvectors of the scattering matrix are needed. The scattered field is expanded as a linear combination of eigensolutions satisfying the outgoing wave condition. The boundary conditions allow the di↵raction amplitudes to be determined.

We propose the specifically designed parallel QR algorithm for the C-method. We present why we propose the "early shift" and how it can be used to accelerate the convergence. We also present the techniques of parallel QR with tightly coupled bulge chasing and parallel AED. These techniques are used to reduce the computational time of the C-method. We apply this specifically designed parallel QR algorithm to the scattering matrix. We also compare the computation time with that of the sequential code.

The results show a significant speed up to approximately 40 for 64 cores with our new QR algorithm. This combination of early shift and other shifts can also be used in the problems such as linear-quadratic optimal control problem where a large number of eigenvalues and eigenvectors are needed and background of the original problem can provides very good initial approximations. This parallel QR algorithm can be used for analyzing crossed gratings or random rough surfaces. Comparisons with experimental data for moderate roughness and isotropic or anisotropic very rough surfaces are conclusive in both co-polarized and cross-polarized components. Comparisons allow the validity of our approach.

As a prospect, we propose a spectral projection method to solve the eigenvalue problem e ciently. We propose a global eigensolver by a combination of the SS method and MIRAMns. This proposed global eigensolver allows us to calculate a large number of (or all) the eigenvalues of a generalized matrix. Compared to QR algorithm, this method has the advantage of having very good scalability. This promising method can be continued in future work. This is the first attempt to combine MIRAMns and SS method to form a global eigensolver. Numerical experiments show this combination allows us to get all the eigenvalues and their corresponding eigenvectors. MIRAMns converges with less iterations than IRAM, and the SS method is very suitable for parallelization. The scalability of the global eigensolver is very good, we get almost linear speed up. The complexity of computation can be varying with the precision that is required. The precision can be increased with smaller sub-domain. .

For gratings, we propose a new version of C-method which leads to a di↵erential system with initial conditions. We studied the new version of C-method as an initial value problem for the e cient calculation of the N-dimensional scattering matrix of a grating.

We have shown that this formulation is an interesting tool for analyzing perfectly conducting or dielectric gratings with deep grooves. The proposed method allows analyzing the complex phenomenon of incident energy absorption. We then extend this method to multilayer with homogeneous medium. We applied this method to multilayer gratings with an arbitrary number of interfaces. We have shown how to combine the local scattering matrix to obtain the global one. We have validate our method by comparison our results with that from published paper. The proposed method has very good accuracy as well as a natural of two level parallel property. This new version of C-method is an attractive alternative to analyze multilayered grating having parallel or non-parallel interfaces.

For the future work, we plan to extend the spectral projection method to more general case. We also plan to extend our new version of C-method to multilayer with inhomogeneous medium. La transformation de la matrice avec un décalage "prématuré" des éléments diagonaux permettant d'accélérer la convergence du processus itératif.

3. 1

 1 Di↵raction by a grating where incident angle ✓ 0 and period D . . . . . . . 3.2 Scattering matrix S. Matrix associates the amplitudes of outgoing plane waves and those of incoming waves . . . . . . . . . . . . . . . . . . . . . 4.1 Computational time relative to truncation order . . . . . . . . . . . . . . 5.1 Intrablock parallel bulge chasing. The grey bulges are chased to the black bulges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Interblock parallel bulge chasing. The gray bulges are chased to the black bulges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Aggressive early deflation. The gray spike contains the vector s . . . . . . 6.1 Function log 10(|error|) relative to M . . . . . . . . . . . . . . . . . . . . 6.2 The green points represent actual eigenvalues, the blue points represent the used early shifts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 The green points represent actually eigenvalues, the blue points represent the values from equation (5.11) . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Computation time of two di↵erent parallelizations relative to number of cores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Computation time with or without early shift relative to number of cores 6.6 Computation time of two di↵erent parallelizations relative to order of matrix 6.7 Computation time with or without early shift relative to order of matrix . 6.8 Comparison of computing time, sequential code relative to parallel code with respect to number of cores . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Comparison of computing time, sequential code relative to parallel code with respect to truncation order . . . . . . . . . . . . . . . . . . . . . . . 6.10 Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 35 , polarization(hh) . . . . . . . . . . . . . . . . . . . . . . . 6.11 Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 35 , polarization(vv) . . . . . . . . . . . . . . . . . . . . . . . 6.12 Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 55 , polarization(hh) . . . . . . . . . . . . . . . . . . . . . . . 6.13 Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 55 , polarization(vv) . . . . . . . . . . . . . . . . . .

8 . 9

 89 Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are:⌫ (1) = ⌫ (3) = 1, ⌫ (3) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 1 = a 2 = 0.02µm,for E // polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.10 Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are:⌫ (1) = ⌫ (3) = 1, ⌫ (2) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 1 = 0.01µm, a 2 = 0.02µm,for E // polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.11 Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are: ⌫ (1) = ⌫ (3) = 1, ⌫ (2) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 1 = 0.03µm, a 2 = 0.02µm, for E // polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.12 Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are: ⌫ (1) = ⌫ (3) = 1, ⌫ (2) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 2 = 0.03µm, a 1 = 0.02µm, for E // polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.13 Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are: ⌫ (1) = ⌫ (3) = 1, ⌫ (2) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 2 = 0.01µm, a 1 = 0.02µm, for E // polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1

  ree charge inside any closed surf ace S (coulomb, C) J = free electric current density (ampere/square meter, A • m 2 )

  .1, which shows a light ray of wavelength incident at an angle ✓ 0 and di↵racted by a grating (of groove spacing D, also called the pitch) along a set of angles ✓ n . These angles are measured from the grating normal. The sign convention for these angles depends on whether the light is di↵racted on the same side or the opposite side of the grating as the incident light.

Figure 3 . 1 :

 31 Figure 3.1: Di↵raction by a grating where incident angle ✓ 0 and period D

n

  ) are the propagation coefficients of the wave vector k(m) n of the elementary plane wave associated with the n th di↵raction order. c (m±) n are the di↵raction amplitudes of elementary plane waves. The propagation coe cient (m) n defines the nature of the plane wave: a propagation wave if (m) n is real, and an evanescent wave if (m) n

Figure 3 . 2 :

 32 Figure 3.2: Scattering matrix S. Matrix associates the amplitudes of outgoing plane waves and those of incoming waves

  48) Equations (3.47) and (3.48) can be obtained from (3.43) and (3.44) by replace Ẽ(1) with Z (1) H(1) and Z (1) H(1) with Ẽ(1) .

  (1) r r e j ⇡ 2 ũ✓ (3.56)Substituting Ẽ(1) by Z(1) H(1) and Z (1) H(1) by Ẽ(1) in equation (3.55) and (3.56), we obtain the H // polarization components of magnetic and electric field vectors. For an incident wave in (a) polarization and a scattered wave in (b) polarization ((a, b) 2

( 3 )

 3 [START_REF] Popov | Surface-enhanced second harmonic generation in nonlinear corrugated dielectrics: new theorical approaches[END_REF]. Harris et al[START_REF] Harris | Di↵erential formalism for mulitlayer di↵raction gratings made with uniaxial meterials[END_REF][START_REF] Harris | Conical di↵raction from multicoated gratings containing uniaxial materials[END_REF], Inchaussandague and Depine[START_REF] Inchaussandague | Polarization conversion from di↵raction gratings made of uniaxial crystals[END_REF][START_REF] Inchaussandague | Rigorous vector theory for di↵raction from gratings made of biaxial crystals[END_REF] have generalized the principle of resolution with anisotropic materials. G. Granet et al investigated the di↵raction gratings with inhomogeneous materials[START_REF] Granet | Extension of the C-method to nonhomogeneous media: applications to nonhomogeneous layers with parallel modulated faces and to inclined lamellar gratings[END_REF]. L. Li et al have proposed a new formalism of C-method to study the interface with edges[START_REF] Li | Improvement of the coordinate transformation method for surface-relief gratings with sharp edges[END_REF]. This formulation is based on the factorization rules of Fourier series and has a faster numerical convergence.

  y)e j↵x e j y d↵d (4.60)

68 )

 68 With a M -th order truncated approximation, the matrices [L l ] and [L r ] are 2M sdimensional ones with M s = (2M + 1) 2 .The elementary solution of equation (4.68) is defined as follows 70) represents an eigenvalue problem, the size of which is N = 2M s . Then,

Figure 4 . 1 :

 41 Figure 4.1: Computational time relative to truncation order

1 :

 1 Set A 0 = A and U 0 = I. 2: for k = 1, 2, 3, ... do 3:

5 . 9 ) 4 .

 594 Continue similar operations to chase the bulge. In general, construct a 3 ⇥ 3 Householder transformation Q k such that the (k + 1)th and (k + 2)th entries of the (k 1)th column of H k are mapped to zero. Applying the corresponding similarity transformation to H k results the updated matrix H k+1 , where k = 2, 3, ..., N 1.

  with the angle defining a propagation direction with ↵n = sin ✓ direction of the associated eigenfunction. A similar interpretation exists for the two-dimensional case.

Figure 5 .

 5 1 shows how the intrablock chasing are performed, the grey bulges are chased to the black bulges. After the bulge chasing within the diagonal block, the accumulated unitary matrices are sent to the corresponding processors in order to update the o↵-diagonal blocks. The o↵-diagonal blocks are then updated by matrix-matrix multiplications which uses level 3 BLAS. The broadcasts are sent in parallel. In order to avoid conflicts in the intersecting parts, they are performed first in the row direction and then in the column direction. See figure5.1 for the intrablock chasing, here (p r , p c ) = (2, 2) and M b = N b = 20.

Figure 5 .

 5 2 illustrates the procedure with (p r , p c ) =[START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Jackson | Classical Electrodynamics[END_REF] and M b = N b = 20, the grey bulges are chased to the black bulges. More precisely, the processor that stores the grey bulges create a copy of the block on each side of the border. Then we can perform the chasing locally, just as in the intrablock chasing and broadcast the corresponding orthogonal factors to the blocks on both sides of the cross border. The updated neighboring block are sent to its owner. To update the corresponding o↵-diagonal blocks, we broadcast orthogonal matrix accumulated in the diagonal chasing stage to the corresponding rows/columns of processors which are involved in o↵-diagonal updating. Then each involved processor exchanges data blocks

Figure 5 . 1 :

 51 Figure 5.1: Intrablock parallel bulge chasing. The grey bulges are chased to the black bulges.

Figure 5 . 2 :

 52 Figure 5.2: Interblock parallel bulge chasing. The gray bulges are chased to the black bulges.

Figure 5 . 3 :

 53 Figure 5.3: Aggressive early deflation. The gray spike contains the vector s

mkl 11 .

 11 0, intelmpi 4.0.3, lapack 3.5 gnu47. The program is written in Fortran and compiled with the following settings: F C = mpif 90, CC = mpicc, NOOP T = 00, F CF LAGS = 03, CCF LAGS = 03,

section 3 . 1 .

 31 2 for details): ⌫ = 1, D = 200 , ✓ = 40 . The surface used is a generated Gaussian rough surface with correlation length l = 3 and standard derivation of height . We set M = 1000, so the matrix size is 4002 ⇥ 4002. In fact, we have performed experiments with di↵erent M to check the error on the power balance. If we denote error = 1 P ✏ n , see (3.12), figure 6.1 shows how the function log 10 (|error|) changes with the truncation order M . It shows if we require a precision of 10 2 , it should be enough to set M = 1000.

Figure 6 . 1 :

 61 Figure 6.1: Function log 10 (|error|) relative to M .

Figure 6 .

 6 Figure 6.4 shows the comparison between the pipeline parallelization and the parallelization with parallel multishift and AED techniques. The item Pipeline parallel represents the performance of PZLAHQR. For our parallel algorithm, the update stage takes approximately 40% of the total time, the Hessenberg translation stage takes approximately 30% to 40% of the total time, the AED stage takes approximately 15% to 25% of the total time and the chasing stage takes approximately 5% to 10% of the total time.

Figure 6 .

 6 Figure 6.5 shows that the early shift does speed up the convergence and save the computation time. With early shift, the computational time of parallel multishift and AED decreases approximately 16%.

Figure 6 . 2 :

 62 Figure 6.2: The green points represent actual eigenvalues, the blue points represent the used early shifts

Figure 6 . 3 :

 63 Figure 6.3: The green points represent actually eigenvalues, the blue points represent the values from equation (5.11)

Figure 6 . 4 :

 64 Figure 6.4: Computation time of two di↵erent parallelizations relative to number of cores

Figure 6 . 5 :

 65 Figure 6.5: Computation time with or without early shift relative to number of cores

Figure 6 . 6 :

 66 Figure 6.6: Computation time of two di↵erent parallelizations relative to order of matrix

Figure 6 . 7 :

 67 Figure 6.7: Computation time with or without early shift relative to order of matrix

Figure 6 . 8 :

 68 Figure 6.8: Comparison of computing time, sequential code relative to parallel code with respect to number of cores

Figure 6 . 9 :

 69 Figure 6.9: Comparison of computing time, sequential code relative to parallel code with respect to truncation order

  [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF] presents a minimum similar to the Brewster angle for a planar surface. (By analogy with reflection from a smooth surface, a lossless dielectric with a refractive index equal to 1.62 provides a Brewster angle close to 58 .)The comparison with experimental data is excellent.

Figure 6 . 10 :

 610 Figure 6.10: Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 35 , polarization(hh)

Figure 6 . 11 :Figure 6 . 12 :Figure 6 . 13 :Figure 6 . 14 :Figure 6 . 15 :Figure 6 . 16 :Figure 6 . 17 :

 611612613614615616617 Figure 6.11: Di↵erential reflection coe cient versus observation angle in the incidence plane, ✓ 0 = 35 , polarization(vv)

Figure 6 . 18 :Figure 6 . 19 :

 618619 Figure 6.18: Average bistatic coe cient versus observation angle in the Oyz incidence plane, very rough anisotropic surface, polarization(hh)

1 :

 1 for l = 1, 2, 3, ... until convergence do 2:

  apply p additional steps of the Arnoldi process to obtain a new m-step Arnoldi factorization,

m i ; 12 :

 12 end for the di↵erent sizes. Another advantage of MIRAMns is that it has better property of convergence with almost the same time complexity compared with IRAM.

Algorithm 7

 7 then the rest of the spectrum lies in the finite domain D = B(0, R) \ B(0, r), where B(x, y) represents the open ball centered at x with radius y in the complex plane. Cover the finite domain D by s subdomains D i , 1  i  s such that D = [ i Di . Use the contour integral based projection method to calculate the eigenvalues in each sub-domain D i to get all the eigenvalues in D. Together with the extremes of the spectrum, we get all the eigenvalues of A.Assuming that the eigenvalues of A are arranged as| 1 |  | 2 |  ...| n 1 |  | n |, the algorithm can be described as in Algorithm 7: Spectral projection method as global eigensolver Input: A, k s , k l , N Output: approximated eigenvalues of A and their associated eigenvectors.

1 : 2 : 3 : 4 :

 1234 Apply MIRAMns to calculate 1 , ..., ks , n+1 k l , ..., n and their associated eigenvectors; Divide the domain D = B(0, R) \ B(0, r) by D = [ i Di ; Apply the SS method to calculate the eigenvalues that lie in D i in parallel and their corresponding eigenvectors; Collect the information from each domain D i and obtain all the eigenvalues and eigenvectors;

Figure 7 . 1 :

 71 Figure 7.1: M IRAM ns[START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] to calculate 3 eigenvalues with largest magnitude, the result is compared with IRAM (16)

Figure 7 . 2 :Figure 7 . 1

 7271 Figure 7.2: M IRAM ns[START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] to calculate 3 eigenvalues with smallest magnitude, the result is compared with IRAM (16)

Figure 7 . 3 :

 73 Figure 7.3: SS-RR method with Center = 0.5, Radius = 0.5, 333 eigenvalues are calculated

Figure 7 .

 7 Figure 7.3 to figure 7.6 are the results of SS-RR method applied to the interval [0, 4].

Figure 7 . 5 :

 75 Figure 7.5: SS-RR method with Center = 2.5, Radius = 0.5, 167 eigenvalues are calculated

Figure 7 . 6 :Figure 7 . 7 :Figure 7 . 8

 767778 Figure 7.6: SS-RR method with Center = 3.5, Radius = 0.5, 333 eigenvalues are calculated

Figure 7 . 8 :Figure 7 . 9 :Figure 7 . 10 :

 7879710 Figure 7.8: M IRAM ns[START_REF] Millar | On the Rayleigh assumption in scattering by a periodic surface -II[END_REF][START_REF] Kleev | The convergence of point-matching techniques[END_REF][START_REF] Bourlier | Theoretical study of the Kirchho↵ integral from a two-dimensional randomly rough surface with shadowing e↵ect: application to the backscattering coe cients for a perfectly-conducting surface[END_REF] to calculate 3 eigenvalues with largest magnitude, the result is compared with IRAM (16)

Figure 7 .

 7 Figure 7.10 to figure 7.19 show how the interval [0, 16] is divided and 1000 eigenvalues are found. It can be noticed that when some of the eigenvalues are very close, the

Figure 7 .

 7 Figure 7.11: SS-RR method with Center = 6, Radius = 2, 136 eigenvalues are calculated, average residual is 5.78977861683748897 ⇥ 10 6

3 Figure 7 . 17 : 4 Figure 7 . 18 :

 37174718 Figure 7.17: SS-RR method with Center = 0.0055, Radius = 0.0045, 45 eigenvalues are calculated, average residual is 2.61507800766467365 ⇥ 10 6

Figure 7 . 20 :

 720 Figure 7.20: Scalability for test matrix B

Figure 8 . 1 :

 81 Figure 8.1: Grating illuminated by a plane wave under incidence ✓ 0 . The space is divided in four regions.

c

  (x, y) within the region y y 1 or y  y 2 by the sum of 2M + 1 outgoing plane waves (amplitude c

n

  (u = 0) are the initial values of field components F (m) (x, u) and G (m) (x, u). We define N independent vectors from di↵erent initial values f

( 2 .

 2 [START_REF] Henry | Parallelizing the qr algorithm for the unsymmetric algebraic eigenvalue problem: Myths and reality[END_REF]) -(2.90). The tangential component of the electric field is zero on the grating surface. The tangential component of the magnetic field is di↵erent to zero and gives the surface current density. For a perfect conduction, we consider 2M + 1 independent vectors and in E // polarization,

  iterative algorithm from u = 0 to u = y m and gives final values of N vectors (contained in matrices F (m) (u = y m ) and G (m) (u = y m )). In the horizontal plane y = y m , the covariant components F (m) (x, u) and G (m) (x, u) become identified with the Cartesian components F (m) c (x, y) and G (m) c (x, y). By projecting on basis functions exp( j↵ n x) of the connection relationship, we find:

Figure 8 .

 8 [START_REF] Jackson | Classical Electrodynamics[END_REF] shows the e ciency curves under the polarization H // with h 0 varying from 0 to . The truncation order M is equal to 9 and y 1 = y 2 = 1.01max(a(x)). The relative and absolute tolerances are equal to 10 6 and 10 9 , respectively. For a perfectly conducting grating illuminated under firstorder Littrow mounting, the e ciency curves when the groove depth increases oscillate between 0 and 1.

Figure 8 . 2 :

 82 Figure 8.2: Reflected e ciencies versus sinusoidal grating amplitude Perfectly conducting grating in H // polarization.

Figure 8 .

 8 Figure 8.3 also shows the error on the power balance for the conventional C-method used with the same value of truncation order. The error given by C1 is lowest when h 0 < 3 /4. When 3 /4 < h 0 < , C1 and C2 methods are equivalent in terms of accuracy on the power balance.

Figure 8 .

 8 Figure 8.4 shows the transmitted e ciencies in H // polarization for a lossless dielectric grating with a sinusoidal profile. The simulation parameters are: ✓ 0 = 15 , ⌫ (m) = 3/2, D = 3 /2 and 0  h 0  . The truncation order is equal to 9 and y 1 = y 2 = 1.01max(a(x)). The relative and absolute tolerances are equal to 10 6 and 10 9 . The grating presents four di↵raction orders: ✓ 2 = 45.8 , ✓ 1 = 15.8 , ✓ 0 = 9.94 and ✓ 1 = 38.1 . For h 0 / = 0, the zeroth-order transmitted e ciency is equal to 0.98 and all other transmitted e ciencies are null. The incident energy is distributed into di↵erent di↵raction orders when the groove depth increases. For h 0 / > 0.42 , the zeroth-order transmitted e ciency is smaller than 50%. Comparison between the reference C-method (C1) used with M = 18 and the new version (C2) used with M = 9 is conclusive and curves of e ciencies are superimposed.

Figure 8 .

 8 Figure 8.5 gives the error on the power balance. The error increases with the groove

Figure 8 . 3 :

 83 Figure 8.3: Error on the power balance versus sinusoidal grating amplitude. Perfectly conducting grating in H // polarization.

Figure 8 . 4 :

 84 Figure 8.4: Transmitted e ciencies versus sinusoidal grating amplitude. Lossless dielectric grating in H // polarization.

Figure 8 . 5 :

 85 Figure 8.5: Error on the power balance versus sinusoidal grating amplitude. Lossless dielectric grating in H // polarization.

Figure 8 .

 8 Figure 8.6 gives the reflected e ciencies versus the depth parameter p for a H-polarized metallic grating whose the profile is defined by,

Figure 8 .

 8 Figure 8.7 gives the sum of e ciencies for the two polarizations. In the E // polarization case, the conduction losses are weak. In H // , for the configuration defined by p = 0.54 , surface plasmons are excited and cause a quasi-total absorption of incident energy[START_REF] Raether | Surface plasmons on smooth and rough surfaces and on gratings[END_REF].

Figure 8 . 6 :

 86 Figure 8.6: Reflected e ciencies versus grooves depth. Metallic grating in H // polarization.

Figure 8 . 7 :

 87 Figure 8.7: Sum of reflected e ciencies versus groove depth. Metallic grating in E // and H // polarizations.

  we use c(m±) to represent a vector containing the scattering amplitudes c (m±) n , m = i, i + 1.

Figure 8 . 8 :

 88 Figure 8.8: Notation for the description of a layered grating

1 A( 8 . 27 )Chapter 8 .

 18278 The C-method as an initial value problem 123 So in this way, we glue S 1,2 and S 2,3 to get S 1,3 , and then glue S 1,3 and S 3,4 to get S 1,4

y = a 1 Figure 8 . 9 :

 189 figure 8.11 shows the curve when a 1 = 0.03.

Figure 8 .

 8 Figure 8.12 and figure 8.13 shows the curve when a 2 = 0.03 and a 2 = 0.01, respectively.

Figure 8 . 10 :Figure 8 . 11 :Figure 8 . 12 :Figure 8 . 13 :

 810811812813 Figure 8.10: Di↵raction e ciency of the zeroth reflected order of the sinusoidal grating. Parameters of the system are:⌫ (1) = ⌫ (3) = 1, ⌫ (2) = 2.3, h 1 = 0.19µm, D = 0.37µm, = 0.6328µm, a 1 = 0.01µm, a 2 = 0.02µm, for E // polarization

  de très grande taille. Cependant, elles ne fournissent pas systématiquement toutes les valeurs propres et leurs vecteurs propres correspondants. Ainsi, ces méthodes itératives sont ine caces pour la méthode C car toutes les valeurs propres et leurs vecteurs propres associés sont nécessaires. La méthode QR qui est basée sur les transformations semblables, calcule tous les éléments propres d'une matrice dense sans danger de manquer des solutions propres particulières. Nous proposons un algorithme QR parallèle conçu spécifiquement pour la méthode C pour résoudre le problème de valeur propre.Cet algorithme QR parallèle est une variante de l'algorithme QR basée sur trois techniques: "early shift" 1 , "bulge chasing" 2[START_REF] Braman | The multishift QR algorithm. part I: Maintaining well-focused shifts and level 3 performance[END_REF][START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF] parallèle et "aggressive early deflation (AED)" 3 parallèle[START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF][START_REF] Braman | The multishift QR algorithm. part II: Aggressive early deflation[END_REF]. Nous proposons la technique "early shift" pour la matrice de di↵usion en fonction des propriétés que nous avons observées. En e↵et, la méthode C et l'interprétation physique derrière la méthode C fournissent une très bonne approximation de certaines valeurs propres avant les calculs. L'utilisation de ces approximations comme "early shift" o↵re la possibilité de déflation rapide. Nous avons combiné le "early shift", le "shift" de Wilkinsons ansi que "le shift exceptionnel" afin d'accélérer la convergence de la méthode QR. Plus particulièrement, nous utilisons le "early shift" afin de déflater les valeurs propres approchées de la matrice de di↵usion et accélérer ainsi la convergence de la méthode. L'algorithme double shift QR, pour des raisons d'économie et d'accélération de convergence, combine deux itérations avec shift en une seule itération avec double shift. A chaque itération, il engendre un bulbe d'éléments non-nuls à chasser par la suite (bulge chasing). Nous utilisons la version multishift de l'algorithme QR. Ainsi, pour le "bulge chasing", au lieu de chasser un seul bulbe, contenant deux shifts, une chaîne de plusieurs bulbes étroitement couplés, contenant chacun deux shifts, est poursuivi au cours d'une itération de QR multishift. Cette idée et la technique de "retard et accumulation"[START_REF] Braman | The multishift QR algorithm. part I: Maintaining well-focused shifts and level 3 performance[END_REF][START_REF] Granat | A novel parallel QR algorithm for hybrid distributed memory HPC systems[END_REF] permet d'e↵ectuer la majeur partie des calculs en termes d'opérations BLAS de niveaux 3 (essentiellement produits matrice-matrice) et augmenter ainsi l'e cacité de l'algorithme en termes de performances. L'AED est une stratégie de la déflation qui profite des perturbations de la matrice en dehors des éléments sous-diagonaux de la matrice de Hessenberg. Elle identifie et déflate les valeurs propres convergées longtemps avant la stratégie classique de déflation et peut améliorer considérablement la convergence de l'algorithme QR. Les résultats présentés dans cette 1

  . . . . . 6.14 Average bistatic coe cient versus observation angle in the incidence plane, very rough isotropic surface, polarization(hh) . . . . . . . . . . . . . . . . 6.15 Average bistatic coe cient versus observation angle in the incidence plane, Average bistatic coe cient versus observation angle in the incidence plane, very rough isotropic surface, polarization(hv) . . . . . . . . . . . . . . . . 6.17 Average bistatic coe cient versus observation angle in the incidence plane, very rough isotropic surface, polarization(vv) . . . . . . . . . . . . . . . . 6.18 Average bistatic coe cient versus observation angle in the Oyz incidence plane, very rough anisotropic surface, polarization(hh) . . . . . . . . . . . 6.19 Average bistatic coe cient versus observation angle in the Oxz incidence plane, very rough anisotropic surface, polarization(hh) . . . . . . . . . . . 7.1 M IRAM ns(9, 12, 16) to calculate 3 eigenvalues with largest magnitude, the result is compared with IRAM (16) . . . . . . . . . . . . . . . . . . . 7.2 M IRAM ns(9, 12, 16) to calculate 3 eigenvalues with smallest magnitude, the result is compared with IRAM (16) . . . . . . . . . . . . . . . . . . . 7.3 SS-RR method with Center = 0.5, Radius = 0.5, 333 eigenvalues are calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 SS-RR method with Center = 1.5, Radius = 0.5, 167 eigenvalues are calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 SS-RR method with Center = 2.5, Radius = 0.5, 167 eigenvalues are calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 SS-RR method with Center = 3.5, Radius = 0.5, 333 eigenvalues are calculated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Scalability for test matrix A . . . . . . . . . . . . . . . . . . . . . . . . . . 7.8 M IRAM ns

very rough isotropic surface, polarization(vh) . . . . . . . . . . . . . . . . vii List of Figures viii 6.

16

  eigenvalues are calculated, average residual is 4.74249918677096298 ⇥ 10 10 . . . . . . 7.19 SS-RR method with Center = 0.00005, Radius = 0.00005, 31 eigenvalues are calculated, average residual is 4.34840427877808130 ⇥ 10 6 . . . . . . 7.20 Scalability for test matrix B . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Grating illuminated by a plane wave under incidence ✓ 0 . The space is divided in four regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Reflected e ciencies versus sinusoidal grating amplitude Perfectly conducting grating in H // polarization. . . . . . . . . . . . . . . . . . . . . . 8.3 Error on the power balance versus sinusoidal grating amplitude. Perfectly conducting grating in H // polarization. . . . . . . . . . . . . . . . . . . . . 8.4 Transmitted e ciencies versus sinusoidal grating amplitude. Lossless dielectric grating in H // polarization. . . . . . . . . . . . . . . . . . . . . . . 8.5 Error on the power balance versus sinusoidal grating amplitude. Lossless dielectric grating in H // polarization. . . . . . . . . . . . . . . . . . . . . . 8.6 Reflected e ciencies versus grooves depth. Metallic grating in H // polarization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Sum of reflected e ciencies versus groove depth. Metallic grating in E // and H // polarizations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.8 Notation for the description of a layered grating . . . . . . . . . . . . . . .

  Media in which equation (2.30) is valid are called ohmic media. A typical example of ohmic media are metals where (2.30) holds in a wide range of circumstances. However, in other materials, such as semiconductors, the

[START_REF] Chandezon | Application d'une nouvelle méthode de résolution des équations de Maxwell à l'étude de la propagation des ondes électromagnétiques dans les guides périodiques[END_REF]

Very often the relation between an electric field and the conduction current density Jc that it generates is given, at any point of the conducting material, by the phenomenological relation, called Ohm's law Jc = Ẽ (2.30) so that J is linearly related to Ẽ through the proportionality factor called the conductivity of the medium. Conductivity is measured in Siemens per meter (S • m 1 = ⌦ 1 • m 1 ) or mhos per meter (mho • m 1 ).

Table 5 . 1 :

 51 Comparison of eigenvalues and early shift

	outgoing waves	incoming waves	r n , M < n < M
	-0.1792 -1.4278i -0.1792 + 1.4278i ±1.8140i 0.1933 -1.4471i 0.1933 + 1.4471i ±1.7044i 0.1126 -1.4278i 0.1126 + 1.4278i ±1.5930i -0.1296 -1.3333i -0.1296 + 1.3333i ±1.4794i 0.0582 -1.2890i 0.0582 + 1.2890i ±1.3629i -0.0944 -1.2226i -0.0944 + 1.2226i ±1.2428i 0.0360 -1.1072i 0.0360 + 1.1072i ±1.1179i 0.0113 -0.9843i 0.0113 + 0.9843i ±0.9865i 0.0006 -0.8504i

Table 5 . 2 :

 52 2D block cyclic scheme

  11: SS-RR method with Center = 6, Radius = 2, 136 eigenvalues are calculated, average residual is 5.78977861683748897 ⇥ 10 6Figure 7.12: SS-RR method with Center = 3, Radius = 1, 95 eigenvalues are calculated, average residual is 2.18255949228685363 ⇥ 10 4 that are very close to the point zero. We can say that there is a clustered point. In fact, the smallest theoretical eigenvalue is 16 cos 4 (1000⇡/2002) = 9.7020 ⇥ 10 11 and |16 cos 4 (1000⇡/2002) 16 cos 4 (999⇡/2002)| = 1.4553 ⇥ 10 9 . These eigenvalues are so close that it is di cult to tell them apart from each other. Thus a very small sub-domain is required to keep a good precision. The initial covering strategy around this clustered point will fail which leads to a poorer performance than that for the matrix A n . Even in this case, figure 7.20 shows a good scalability of the proposed global eigensolver.Figure 7.13: SS-RR method with Center = 1.5, Radius = 0.5, 72 eigenvalues are calculated, average residual is 2.59612185434167264 ⇥ 10 6Figure 7.15: SS-RR method with Center = 0.3, Radius = 0.2, 95 eigenvalues are calculated, average residual is 2.59080105213133636 ⇥ 10 4

	0.5 0.2
	0.4
	0.15
	0.3
	0.1
	0.2
	0.05
	0.1
	0 0
	-0.1
	-0.05
	-0.2
	-0.1
	-0.3
	-0.15
	-0.4
	-0.5 -0.2
	1 0.1	1.1 0.15	1.2	0.2	1.3	1.4 0.25	1.5 0.3	1.6 0.35	1.7	0.4	1.8	1.9 0.45	2 0.5
	1 0.25 0.05
	0.8 0.2 0.04
	0.6 0.15 0.03
	0.4 0.1 0.02
	0.2 0.05 0.01
	0 0 0
	-0.2 -0.05 -0.01
	-0.4 -0.1 -0.02
	-0.6 -0.15 -0.03
	-0.8 -0.2 -0.04
	-1 -0.25 -0.05
	2 0.5 0.01	2.2 0.55 0.02	2.4 0.6 0.03	2.6 0.65 0.04	2.8 0.7 0.05 0.75 3 0.06 3.2 0.8	3.4 0.85 0.07	3.6 0.9 0.08	3.8 0.95 0.09	4 1 0.1
	Figure 7.14: SS-RR method with Center = 0.75, Radius = 0.25, 57 eigenvalues are Figure 7.16: SS-RR method with Center = 0.055, Radius = 0.045, 80 eigenvalues
	calculated, average residual is 1.45660089462261519 ⇥ 10 6 are calculated, average residual is 1.40601714058967572 ⇥ 10 5
	7.3 Conclusion This is the first attempt to combine MIRAMns and SS method to form a global eigen-
	solver. Numerical experiments show this combination allows us to get all the eigen-
	In this chapter, we propose a global eigensolver by combination of the contour integral values and their corresponding eigenvectors. MIRAMns converges with less iterations
	based projection method (SS method) and the multiple implicitly restarted Arnoldi than IRAM, and the SS method is very suitable for parallelization. The scalability of
	method with nested subspaces (MIRAMns). This proposed global eigensolver allows the global eigensolver is very good, we get almost linear speed up. The complexity of
	us to calculate a large number of (or all) eigenvalues and eigenvectors of a generalized computation can be varying with the precision that is required. The precision can be
	matrix. increased with smaller sub-domain.

  Figure 7.19: SS-RR method with Center = 0.00005, Radius = 0.00005, 31 eigenvalues are calculated, average residual is 4.34840427877808130 ⇥ 10 6

		x 10 -5
		5
		4
		3
		2
		1
		0
		-1
		-2
		-3
		-4
		-5
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
											x 10 -4
		0.22
		0.2
		0.18
		0.16
		0.14
	1/time (1/s)	0.1 0.12
		0.08
		0.06
		0.04
		0.02
		0
		0	2	4	6	8	10	12	14	16	18	20
							#cores

La chasse aux élément non-nuls introduits au cours du calcul

La dimiuntion de la taille de la matrice en fonction des valeurs propres connues et/ou déjà calculées

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisors Prof. Richard Dusséaux and Prof. Nahid Emad for the continuous support of my Ph.D study and related research, for their patience, motivation, and immense knowledge. Their guidances helped me in all the time of research and writing of this thesis. I could not have imagined having better advisors and mentors for my Ph.D study. Besides my advisors, I would like to thank the reviewers of my thesis: Porf. James Demmel and Dr. Denis Prémel, for their insightful comments and encouragement, but also for the hard question which incented me to widen my research from various perspectives. I would like also to thank the rest of my thesis committee: Dr. Christophe Bourlier, Dr. Philippe Ricoux and Prof. Valérie Ciarletti.

I thank my fellow labmates in for the stimulating discussions, for all the fun we have had in the last four years.

Chapter 4

The curvilinear coordinate method [START_REF] Ishimaru | Electromagnetic wave propagation, radiation and scattering[END_REF]

.1 Introduction

The resolution of Maxwell's equations requires to consider the continuity of certain components of the field on the interface. The continuity relation is simplified if the interface is one coordinate surface. If the interface is geometrically simple (plane, cylinder, sphere), we can use the corresponding coordinates. If it is not, in order to see clearly the continuity relations, we need to use nonorthogonal coordinates.

Several authors have adapted the Maxwell-Minkowsky equations to the three-dimensional space. By combining the electromagnetic field tensor, Minkowsky has generalized the Maxwell's equations to the space-time. E.J.Post has written the equations in the rationalized MKS system [START_REF] Post | Formal structure of electromagnetic[END_REF]. J.Chandezon et al have proposed to adapt this formulation to the three-dimensional space and valid for all the curvilinear coordinates [START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF]. Under this form, the Maxwell's equations are not a↵ected by the coordinate system which is di↵erent from classical ways.

In order to use this formalism, J.Chandezon et al introduced two systems of nonorthogonal coordinates, the translation system to study the di↵raction of a plane wave by a grating [START_REF] Chandezon | A new theoretical method for di↵raction gratings and its numerical application[END_REF] and the revolution system to study the propagation of wave in the periodical cylindrical guides [START_REF] Chandezon | Application d'une nouvelle méthode de résolution des équations de Maxwell à l'étude de la propagation des ondes électromagnétiques dans les guides périodiques[END_REF]. The principle of the proposed method is called the C-method. Nous nous concentrons sur l'aspect numérique de la méthode C, en développant une mise en oeuvre e cace de cette méthode exacte. Des méthodes itératives de recherche de valeurs propres telles que les méthodes de sous-espace de Krylov ou les méthodes de Jacobi-Davidson [START_REF] Bai | Templates for the solution of algebraic eigenvalue problems[END_REF] ont été développées pour traiter de problèmes de valeur propre thèse mettent en évidence cette amélioration de performances pour le problème considéré.

En perspective, nous proposons une méthode de projection spectrale pour résoudre le problème de valeurs propres e cacement. Cette méthode proposeé afin de palier au problème de "scalability" de la méthode QR. Elle est baseé sur une combinaison de la méthode de Sakuria et Sugiura (SSM) [START_REF] Sakuria | A projection method for generalized eigenvalue problem using numerical integration[END_REF] et "multiple implicitly restarted Arnolid method" avec des sous-espace imbriqué (MIRAMns) proposé par S. A. Shahzadeh Fazeli et al [START_REF] Shahzadeh Fazeli | A key to choose subspace size in implicitly restarted arnoldi method[END_REF]. La méthode proposée nous permet de calculer un grand nombre de (ou