Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy
email: guillaume.demesy@fresnel.fr

John Desanto

Gérard Granet

Editor Evgeny Popov

Boris Gralak

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout

Fréderic Zolla

Benjamin Vial

Guillaume Demésy

Frédéric Zolla

GratinGs: theory and numeric applications

Introduction

Finite element methods (FEM) represent a very general set of techniques to approximate solutions of partial derivative equations. Their main advantage lies in their ability to handle arbitrary geometries via unstructured meshes of the domain of interest: The discretization of oblic geometry edges is natively built in. Finite Element Methods have been widely developed in many areas of physics and engineering: mechanics, thermodynamics. . . But until the early 80's, two major drawbacks prevented them from being used in electromagnetic problems. On the one hand, existing nodal element basis did not satisfy the physical (dis)continuity of the vector fields components and lead to spurious solutions [START_REF] Bossavit | Solving Maxwell equations in a closed cavity, and the question ofspurious modes[END_REF]. On the other hand, there was no proper way to truncate unbounded regions in open wave problems.

These two major limitations were both overcome in the early 80's: Vector elements have been developed by Nédélec [START_REF] Nedelec | Mixed finite elements in R 3[END_REF][START_REF] Bossavit | Solving maxwell equations in a closed cavity, and the question ofspurious modes[END_REF], and Perfectly Matched Layers (PMLs) were discovered by Bérenger [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Since then, it has been shown that PMLs could be described in the general framework of transformation optics [START_REF] Chew | A 3d perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF][START_REF] Teixeira | General closed-form pml constitutive tensors to match arbitrary bianisotropic and dispersive linear media[END_REF][START_REF] Nicolet | Geometrical transformations and equivalent materials in computational electromagnetism[END_REF][START_REF] Agha | On the use of pml for the computation of leaky modes: An application to microstructured optical fibres[END_REF].

All the mathematical and computational ingredients now exist and the goal of this chapter is to show how to combine them to implement a general 3D numerical scheme adapted to gratings using Finite Elements. In fact, we are now facing the physical difficulties inherent to the infinite spatial characteristics of the grating problem, whereas the computation domain has to be bounded in practice: (i) Both the superstrate and the substrate are infinite regions, (ii) there is an infinite number of periods and, last but not least, (iii) the sources of the incident field (a plane wave) are located in the superstrate at an infinite distance from the grating.

In this chapter, the infinite extension of the superstrate and substrate is addressed using cartesian PMLs. In the framework of transformation optics, we demonstrate that Bérenger's original PMLs can be extended to the challenging numerical cases of grazing incidence in order to deal with extreme oblic incidences or configurations near Wood's anomalies. The second issue of infinite number of period can be addressed via Bloch conditions. Finally, we are dealing with the distant plane wave sources through an equivalence of the diffraction problem with a radiation one whose sources are localized inside the diffractive element itself. The unknown field to be approximated using Finite Elements is a radiated field with sources inside the computation box and allows to retrieve easily the total field with the plane wave source.

In a first section, we derive and implement this approach in the so-called 2D non-conical, or scalar, case. We are dealing with the infinite issues rigorously in both TE and TM polarization cases. It results in a radiation problem with sources localized in the diffractive element itself. We mathematically split the whole problem into two parts. The first one consists in the classical calculation of the total field solution of a simple interface. The second one amounts to looking for a radiated field with sources confined within the diffractive obstacles and deduced from the first elementary problem. From this viewpoint, the later radiated field can be interpreted as an exact perturbation of the total field. We show that our approach allows to tackle some kind of anisotropy without increasing the computational time or resource. Through a battery of examples, we illustrate its independence towards the geometry of the diffractive pattern. Finally, we present an Adaptative PML able to tackle grazing incidences or configurations near Wood's anomaly.

In a second section, we extend this approach to the most general configuration of vector diffraction by crossed gratings embedded in arbitrary multilayered stack. The main advantage of this method is, again, its complete independence towards the shape of the diffractive element, whereas other methods often require heavy adjustments depending on whether the geometry of the groove region presents oblique edges. This approach combined with the use of second order edge elements allows us to retrieve the few numerical academic examples found in the literature with an excellent accuracy. Furthermore, we provide a new reference case combining major difficulties: A non trivial toroidal geometry together with strong losses and a high permittivity contrast. Finally, we discuss computation time and convergence as a function of the mesh refinement as well as the choice of the direct solver.

5.2 Scalar diffraction by arbitrary mono-dimensional gratings : a Finite Element formulation

Set up of the problem and notations

We denote by x, y and z, the unit vectors of the axes of an orthogonal coordinate system Oxyz.

We deal only with time-harmonic fields; consequently, the electric and magnetic fields are represented by the complex vector fields E and H, with a time dependance in exp(-i ω t).

Besides, in this chapter, we assume that the tensor fields of relative permittivity ε and relative permeability µ can be written as follows:

ε =   ε xx εa 0 ε a ε yy 0 0 0 ε zz   and µ =   µ xx μa 0 µ a µ yy 0 0 0 µ zz   , (5.1)
where ε xx , ε a , . . . µ zz are possibly complex valued functions of the two variables x and y and where εa (resp. μa ) represents the conjugate complex of ε a (resp. µ a ). These kinds of materials are said to be z-anisotropic. It is of importance to note that with such tensor fields, lossy materials can be studied (the lossless materials correspond to tensors with real diagonal terms represented by Hermitian matrices) and that the problem is invariant along the z-axis but the tensor fields can vary continuously (gradient index gratings) or discontinuously (step index gratings). Moreover we define k 0 := ω/c. The gratings that we are dealing with are made of three regions (See Fig. 5.1 ).

• The superstratum (y > h g ) which is supposed to be homogeneous, isotropic and lossless and characterized solely by its relative permittivity ε + and its relative permeability µ + and we denote k + := k 0 ε + µ +

• The substratum (y < 0) which is supposed to be homogeneous and isotropic and therefore characterized by its relative permittivity ε -and its relative permeability µ -and we denote k -:= k 0 ε -µ -

• The groove region (0 < y < h g ) which can be heterogeneous and z-anisotropic and thus characterized by the two tensor fields ε g (x, y) and µ g (x, y). It is worth noting that the method does work irrespective of whether the tensor fields are piecewise constant. The groove periodicity along x-axis will be denoted d. This grating is illuminated by an incident plane wave of wave vector k + p = α xβ + y = k + (sin θ 0 xcos θ 0 y), whose electric field (TM case) ( resp. magnetic field (TE case)) is linearly polarized along the z-axis:

E 0 e = A 0 e exp(ik + p • r r r) z (resp. H 0 m = A 0 m exp(ik + p • r r r) z) , (5.2)
where A 0 e (resp. A 0 m ) is an arbitrary complex number and r r r = (x, y) T . In this section, a plane wave is characterized by its wave-vector denoted k {+,-} {p,c} . The subscript p (resp. c) stands for "propagative" (resp. "counter-propagative"). The superscript + (resp. -) refers to the associated wavenumber k + (resp. k -), and indicates that we are dealing with a plane wave propagating in the superstrate (resp. substrate).The magnetic (resp. electric) field derived from E 0 e (resp. H 0 m ) is denoted H 0 e (resp. E 0 m ) and the electromagnetic field associated with the incident field is therefore denoted (E 0 , H 0 ) which is equal to (E 0 e , H 0 e ) (resp. (E 0 m , H 0 m )). The diffraction problem that we address consists in finding Maxwell equation solutions in harmonic regime i.e. the unique solution (E, H) of:

curl E = i ω µ 0 µ H (5.3a) curl H = -i ω ε 0 ε E (5.3b) such that the diffracted field (E d , H d ) := (E -E 0 e , H -H 0 m
) satisfies an Outgoing Waves Condition (O.W.C. [START_REF] Petit | Electromagnetic theory of gratings[END_REF]) and where E and H are quasi-periodic functions with respect to the x coordinate.

Theoretical developments of the method

Decoupling of fields and z-anisotropy

We assume that δ (x, y) is a z-anisotropic tensor field (δ xz = δ yz = δ zx = δ zy = 0). Moreover, the left upper matrix extracted from δ is denoted δ , namely:

δ = δ xx δa δ a δ yy . (5.4)
For z-anisotropic materials, in a non-conical case, the problem of diffraction can be split into two fundamental cases (TE case and TM case). This property results from the following equality which can be easily derived:

-curl δ -1 curl (u z) = div δ T / det( δ )∇u z , (5.5)
where u is a function which does not depend on the z variable. Relying on the previous equality, it appears that the problem of diffraction in a non conical mounting amounts to looking for an electric (resp. magnetic) field which is polarized along the z-axis ; E = e(x, y) z (resp.

H = h(x, y) z).
The functions e and h are therefore solutions of similar differential equations:

L ξ ,χ (u) := div ξ ∇u + k 2 0 χ u = 0 (5.6) with u = e, ξ = μT / det( μ), χ = ε zz , (5.7)
in the TM case and

u = h, ξ = εT / det(ε), χ = µ zz , (5.8)
in the TE case.

Boiling down the diffraction problem to a radiation one

In its initial form, the diffraction problem summed up by Eq. (5.6) is not well suited to the Finite Element Method. In order to overcome this difficulty, we propose to split the unknown function u into a sum of two functions u 1 and u d 2 , the first term being known as a closed form and the latter being a solution of a problem of radiation whose sources are localized within the obstacles.

We have assumed that outside the groove region (cf. Fig. 5.1), the tensor field ξ and the function χ are constant and equal respectively to ξ -and χ -in the substratum (y < 0) and equal respectively to ξ + and χ + in the superstratum (y > h g ). Besides, for the sake of clarity, the superstratum is supposed to be made of an isotropic and lossless material and is therefore solely defined by its relative permittivity ε + and its relative permeability µ + , which leads to:

ξ + = 1
µ + Id 2 and χ + = ε + in TE case (5.9) or ξ + = 1 ε + Id 2 and χ + = µ + in TM case, (5.10) where Id 2 is the 2 × 2 identity matrix. With such notations, ξ and χ are therefore defined as follows:

ξ (x, y) :=      ξ + for y > h g ξ g (x, y) for h g > y > 0 ξ - for y < 0 , χ(x, y) :=    χ + for y > h g χ g (x, y) for h g > y > 0 χ - for y < 0 .
(5.11)

It is now apropos to introduce an auxiliary tensor field ξ 1 and an auxiliary function χ 1 :

ξ 1 (x, y) := ξ + for y > 0 ξ -for y < 0 , χ 1 (x, y) := χ + for y > 0 χ -for y < 0 , (5.12)
these quantities corresponding, of course, to a simple plane interface. Besides, we introduce the constant tensor field ξ 0 which is equal to ξ + everywhere and a constant scalar field χ 0 which is equal to χ + everywhere. Finally, we denote u 0 the function which equals the incident field u inc in the superstratum and vanishes elsewhere (see Fig. 5.1):

u 0 (x, y) := u inc for y > h g 0 for y < h g (5.13)
We are now in a position to define more precisely the diffraction problem that we are dealing with. The function u is the unique solution of:

L ξ ,χ (u) = 0 , such that u d := u -u 0 satisfies an O.W.C. ( 5.14)
In order to reduce this diffraction problem to a radiation problem, an intermediate function is necessary. This function, called u 1 , is defined as the unique solution of the equation:

L ξ 1 ,χ 1 (u 1 ) = 0 , such that u d 1 := u 1 -u 0 satisfies an O.W.C. (5.15)
The function u 1 corresponds thus to an annex problem associated to a simple interface and can be solved in closed form and from now on is considered as a known function. As written above, we need the function u d 2 which is simply defined as the difference between u and u 1 :

u d 2 := u -u 1 = u d -u d 1 .
(5.16)

The presence of the superscript d is, of course, not irrelevant: As the difference of two diffracted fields, the O.W.C. of u d 2 is guaranteed (which is of prime importance when dealing with PML cf. 5.2.2.4). As a result, the Eq. (5.14) becomes:

L ξ ,χ (u d 2 ) = -L ξ ,χ (u 1 ) , (5.17)
where the right hand member is a scalar function which may be interpreted as a known source term -S 1 (x, y) and the support of this source is localized only within the groove region. To prove it, all we have to do is to use Eq. (5.15):

S 1 := L ξ ,χ (u 1 ) = L ξ ,χ (u 1 ) -L ξ 1 ,χ 1 (u 1 ) =0 = L ξ -ξ 1 ,χ-χ 1 (u 1 ) . (5.18)
Now, let us point out that the tensor fields ξ and ξ 1 are identical outside the groove region and the same holds for χ and χ 1 . The support of S 1 is thus localized within the groove region as expected. It remains to compute more explicitly the source term S 1 . Making use of the linearity of the operator L and the equality u 1 = u d 1 + u 0 , the source term can be split into two terms

S 1 = S 0 1 + S d 1 , (5.19)
where

S 0 1 = L ξ -ξ 1 ,χ-χ 1 (u 0 ) (5.20) and S d 1 = L ξ -ξ 1 ,χ-χ 1 (u d 1 ) . (5.21)
Now, bearing in mind that u 0 is nothing but a plane wave u 0 = exp(ik + p • r) (with k + p = αxβ + y), it is sufficient to give ∇u 0 = ik + p u 0 for the weak formulation associated with Eq. (5.17):

S 0 1 = i div ξ + -ξ k + p exp(ik + p • r) + k 2 0 χ + -χ exp(ik + p • r) . (5.22)
The same holds for the term associated with the diffracted field. Since, in the superstrate, we have of course

u d 1 = ρ exp(ik + c • r) with k + c = αx + β + y, S d 1 = ρ i div ξ + -ξ k + c exp(ik + c • r) + k 2 0 χ + -χ exp(ik + c • r) , (5.23)
where ρ is simply the complex reflection coefficient associated with the simple interface:

ρ = p + -p - p + -p -with p ± =    β ± in the TM case β ±
ε ± in the TE case (5.24)

Quasi-periodicity and weak formulation

The weak formulation follows the classical lines and is based on the construction of a weighted residual of Eq. (5.6), which is multiplied by the complex conjugate of a weight function u and integrated by part to obtain:

R ξ ,χ (u, u ) = - Ω ξ ∇u • ∇u + k 2 0 χ u u dΩ + ∂ Ω u ξ ∇u • n dS (5.25)
The solution u of the weak formulation can therefore be defined as the element of the space L 2 (curl, d, α) of quasiperiodic functions (i.e. such that u(x, y) = u # (x, y) e iαx with u # (x, y) = u # (x + d, y), a d-periodic function) of L 2 (curl) on Ω such that: therefore an unknown to be determined on the boundary). A third possibility are the so-called quasi-periodicity conditions of particular importance in the modeling of gratings. Denote by Γ l and Γ r the lines parallel to the y-axis delimiting a cell of the grating (see Fig. 5.2) respectively from its left and right neighbor cell. Considering that both u and u are in

R ξ ,χ (u, u ) = 0 ∀u ∈ L 2 (curl, d, α). ( 5
L 2 (curl, d, α), the boundary term for Γ l ∪ Γ r is Γ l ∪Γ r u ξ ∇u • n dS = Γ l ∪Γ r u # e -iαx ξ ∇(u # e +iαx ) • n dS = Γ l ∪Γ r u # ξ (∇u # + iαu # x) • n dS = 0 , because the integrand u # ξ (∇u # + iαu # x)
•n is periodic along x and the normal n has opposite directions on Γ l and Γ r so that the contributions of these two boundaries have the same absolute value with opposite signs. The contribution of the boundary terms vanishes therefore naturally in the case of quasi-periodicity.

The finite element method is based on this weak formulation and both the solution and the weight functions are classically chosen in a discrete space made of linear or quadratic Lagrange elements, i.e. piecewise first or second order two variable polynomial interpolation built on a triangular mesh of the domain Ω (cf. Fig. 5.3a). Dirichlet and Neumann conditions may be used to truncate the PML domain in a region where the field (transformed by the PML) is negligible. The quasi-periodic boundary conditions are imposed by considering the u as unknown on Γ l (in a way similar to the homogeneous Neumann condition case) while, on Γ r , u is forced equal to the value of the corresponding point on Γ l (i.e. shifted by a quantity -d along x) up to the factor e iαd . The practical implementation in the finite element method is described in details in [START_REF] Zolla | Foundations of Photonic Crystal Fibres: 2nd Edition[END_REF][START_REF] Nicolet | Modelling of electromagnetic waves in periodic media with finite elements[END_REF]

Perfectly Matched Layer for z-anisotropic materials

The main drawback encountered in electromagnetism when tackling theory of gratings through the finite element method is the non-decreasing behaviour of the propagating modes in superstratum and substratum (if they are made of lossless materials): The PML has been introduced by [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] in order to get round this obstacle. The computation of PML designed for z-anisotropic gratings is the topic of what follows.

In the framework of transformation optics, a PML may be seen as a change of coordinate corresponding to a complex stretch of the coordinate corresponding to the direction along which the field must decay [START_REF] Nicolet | Leaky modes in twisted microstructured optical fibres[END_REF][START_REF] Lassas | Complex riemannian metric and absorbing boundary condition[END_REF][START_REF] Lassas | Analysis of the PML equations in general convex geometry[END_REF]. Transformation optics have recently unified various techniques in computational electromagnetics such as the treatment of open problems, helicoidal geometries or the design of invisibility cloaks ( [START_REF] Nicolet | Geometrical transformations and equivalent materials in computational electromagnetism[END_REF]). These apparently different problems share the same concept of geometrical transformation, leading to equivalent material properties. A very simple and practical rule can be set up ( [START_REF] Zolla | Foundations of Photonic Crystal Fibres: 2nd Edition[END_REF]): when changing the coordinate system, all you have to do is to replace the initial materials properties ε and µ by equivalent material properties ε s and µ s given by the following rule:

ε s = J -1 ε J -T det(J) and µ s = J -1 µ J -T det(J), (5.27)
where J is the Jacobian matrix of the coordinate transformation consisting of the partial derivatives of the new coordinates with respect to the original ones (J -T is the transposed of its inverse).

In this framework, the most natural way to define PMLs is to consider them as maps on a complex space C 3 , which coordinate change leads to equivalent permittivity and permeability tensors. We detail here the different coordinates used in this section.

• (x, y, z) are the cartesian original coordinates.

• (x s , y s , z s ) are the complex stretched coordinates. A suitable subspace Γ ⊂ C 3 is chosen (with three real dimensions) such that (x s , y s , z s ) are the complex valued coordinates of a point on Γ (e.g. x = Re(x s ), y = Re(y s ), z = Re(z s )).

• (x c , y c , z c ) are three real coordinates corresponding to a real valued parametrization of Γ ⊂ C 3 .

We use rectangular PMLs ( [START_REF] Nicolet | Leaky modes in twisted microstructured optical fibres[END_REF]) absorbing in the y-direction and we choose a diagonal matrix J = diag(1, s y (y), 1), where s y (y) is a complex-valued function of the real variable y, defined by: y s (y) = y 0 s y (y )dy .

(5.28)

The expression of the equivalent permittivity and permeability tensors are thus:

ε s =   s y ε xx ε a 0 ε a s -1 y ε yy 0 0 0 s y ε zz   and µ s =   s y µ xx µ a 0 µ a s -1 y µ yy 0 0 0 s y µ zz   .
(5.29)

Note that the equivalent medium has the same impedance than the original one as ε an µ are transformed in the same way, which guarantees that the PML is perfectly reflectionless. Now, let us define the so-called substituted field F F F s = (E s , H s ), solution of Eqs. (5.3) with ξ = ξ s and χ = χ s . It turns out that F F F s equals the field F F F in the region y b < y < y t (with y b = -h - and y t = h g + h + , see Fig. 5.3a), provided that s y (y) = 1 in this region. The main feature of this latest field F F F s is the remarkable correspondence with the first field F F F ; whatever the function s y provided that it equals 1 for y t < y < y b , the two fields F F F and F F F s are identical in the region y t < y < y b [START_REF] Agha | On the use of pml for the computation of leaky modes: An application to microstructured optical fibres[END_REF]. In other words, the PML is completely reflection-less. In addition, for complex valued functions s y (ℑm{s y } strictly positive in PML), the field F F F s converges exponentially towards zero (as y tends to ±∞, cf. Fig. 5.3c and 5.3d) although its physical counterpart F F F does not. Note that in Fig. 5.3d, the value of the computed radiated field u d 2 on each extreme boundary of the PMLs is at least 10 -8 weaker than in the region of interest. As a consequence, F F F s is of finite energy and for this substituted field a weak formulation can be easily derived which is essential when dealing with Finite Element Method.

Still remains to give a suitable function s y . Let us consider the complex coordinate mapping y(y c ), which is simply defined as the derivative of the stretching coefficient s y (y) with respect to y c . With simple stretching functions, we can obtain a reliable criterion upon proper fields decay. A classical choice is:

s y (y) =      ζ -if y < y b 1 if y b < y < y t ζ + if y > y t
(5.30)

where ζ ± = ζ ,± + iζ ,± are complex constants with ζ ,± > 0.

In that case, the complex valued function y(y c ) defined by Eq. (5.28) is explicitly given by:

y(y c ) =      y b + ζ -(y c -y b ) if y c < y b y c if y b < y c < y t y t + ζ + (y c -y t ) if y c > y t , (5.31)
Finally, let us consider a propagating plane wave in the substratum u n (x, y) := exp(i(αxβ - n y)). Its expression can be rewritten as a function of the stretched coordinates in the PML as follows:

u (5.32)

The behavior of this latest function along the y c direction is governed by the function

U sc (y c ) := e -iβ - n ζ -y c . Letting β ,- n := ℜe{β - n }, β ,- n := ℑm{β - n }, ζ ,-:= ℜe{ζ -} and ζ ,-:= ℑm{ζ -}, the non-oscillating part of the function U sc (y c ) is given by exp (β ,- n ζ ,-+ β ,- n ζ ,-)y c . Keeping in mind that β ,-
n and/or β ,- n are positive numbers, the function U sc decreases exponentially towards zero as y c tends to -∞ (Fig. 5.3d) provided that ζ -belongs to C + := {z ∈ C, ℜe{z} > 0, and ℑm{z} > 0}. In the same way, it can be shown that ζ + belongs to C + .

Let us conclude this section with two important remarks:

1. Practical choice of PML parameters. As for the complex stretch parameters, setting ζ ± = 1 + i is usually a safe choice. For computational needs, the PML has to be truncated and the other constitutive parameter of the PML is its thickness ĥ (see Fig. 5.3a). Setting ĥ± = λ 0 / √ ε ± leads to a PML thick enough to "absorb" all incident radiation. These specific values will be used in the sequel, unless otherwise specified.

2. Special cases. The reader will notice that a configuration where β ,- n is a very weak positive number compared to k 0 with β ,- n (this is precisely the case of a plane wave at grazing incidence on the bottom PML) leads to a very slow exponential decay of U sc . In such a case, close to so-called Wood's anomalies or at extreme grazing incidences, classical PML fail. We will address this tricky situation extensively in Section 5.2.4.

Synthesis of the method

In order to give a general view of the method, all information is collected here that is necessary to set up the practical Finite Element Model. First of all, the computation domain Ω (cf. Fig. 5.3a) corresponds to a truncated cell of the grating which is a finite rectangle divided into five horizontal layers. These layers are respectively from top to bottom upper PML, the superstratum, the groove region, the substratum, and the lower PML. The unknown field is the scalar function u d 2 defined in Eq. (5.16). Its finite element approximation is based on the second Lagrange elements built on a triangle meshing of the cell (cf. Fig. 5.3b). A complex algebraic system of linear equations is constructed via the Galerkin weighted residual method, i.e. the set of weight functions u is chosen as the set of shape functions of interpolation on the mesh [START_REF] Zolla | Foundations of Photonic Crystal Fibres: 2nd Edition[END_REF].

• In region 1 (upper PML, see Fig. 5.3a),

R ξ + s ,χ + s (u d 2 , u ) = 0 , (5.33)
with ξ + s and χ + s depending on the equivalent anisotropic properties of the PML given by Eq. (5.7), Eq. (5.8) and Eqs. (5.29).

• In region 2 (superstratum),

R ξ + ,χ + (u d 2 , u ) = 0 , (5.34)
with ξ + and χ + depending on the homogeneous isotropic properties of the superstratum given by Eq. (5.7), Eq. (5.8), Eq. (5.9) and Eq. (5.10).

• In region 3 (groove region),

R ξ g ,χ g (u d 2 , u ) = -R ξ g ,χ g (S 1 , u ) , (5.35)
with ξ g and χ g depending on the heterogeneous possibly anisotropic properties given by Eq. (5.7), Eq. (5.8), Eq. (5.11) and S 1 given by Eq. (5.19) , Eq. (5.22), Eq. (5.23) and Eq. (5.24).

• In region 4 (substratum),

R ξ -,χ -(u d 2 , u ) = 0 , (5.36)
with ξ -and χ -depending on the homogeneous isotropic properties of the substratum given by Eq. (5.7), Eq. (5.8), Eq. (5.9) and Eq. (5.10).

• In region 5 (lower PML),

R ξ - s ,χ - s (u d 2 , u ) = 0 , (5.37)
with ξ - s and χ - s depending on the equivalent anisotropic properties of the PML given by Eq. (5.7), Eq. (5.8) and Eqs. (5.29).

Energy balance: Diffraction efficiencies and absorption

The rough result of the FEM calculation is the complex radiated field u d 2 . Using Eq. ( 5.16), it is straightforward to obtain the complex diffracted field u d solution of Eq. ( 5.14) at each point of the bounded domain. We deduce from u d the diffraction efficiencies with the following method. The superscripts + (resp. -) correspond to quantities defined in the superstratum (resp. substratum) as previously.

On the one hand, since u d is quasi-periodic along the x-axis , it can be expanded as a Rayleigh expansion (see for instance [START_REF] Petit | Electromagnetic theory of gratings[END_REF]):

for y < 0 and y > h g , u d (x, y) = ∑ n∈Z u d n (y) e iα n x (5.38)
where

u d n (y) = 1 d d/2 -d/2 u d (x, y)e -iα n x dx with α n = α + 2π d n (5.39)
On the other hand, introducing Eq. ( 5.38) into Eq. (5.6) leads to the Rayleigh coefficients:

u d n (y) =    u + n (y) = r n e iβ + n y + a n e -iβ + n y for y > h g u - n (y) = t n e -iβ - n y + b n e iβ - n y for y < 0 with β ± 2 n = k ± 2 -α 2 n (5.40)
For a temporal dependance in e -iωt , the O.W.C. imposes a n = b n = 0. Combining Eq. (5.39) and (5.40) at a fixed y 0 altitude leads to:

           r n = 1 d d/2 -d/2 u d (x, y 0 ) e -i(α n x+β + n y 0 ) dx for y 0 > h g t n = 1 d d/2 -d/2 u d (x, y 0 ) e -i(α n x-β - n y 0 ) dx for y 0 < 0 (5.41)
We extract these two coefficients by trapezoidal numerical integration along x from a cutting of the previously calculated field map at y 0 . It is well known that the mere trapezoidal integration method is very efficient for smooth and periodic functions (integration on one period) [START_REF] Helluy | Intégrations numériques d'ordre élevé de fonctions régulières ou singulières sur un intervalle[END_REF]. Now the restriction on a horizontal straight line crossing the whole cell in homogeneous media (substratum and superstratum) is of C ∞ class. From a numerical point of view, it appears that the interpolated approximation of the unknown function, namely u d 2 preserves the good behaviour of the numerical computation of these integrals. From this we immediately deduce the reflected and transmitted diffracted efficiencies of propagative orders (T n and R n ) defined by:

     R n := r n r n β + n β + for y 0 > h g T n := t n t n β - n β - γ + γ -for y 0 < 0 with γ ± =    1 in the TM case ε ± in the TE case (5.42)
This calculation is performed at several different y 0 altitudes in the superstratum and the substratum, and the mean value found for each propagative transmitted or reflected diffraction order is presented in the numerical experiments of the following section.

Normalized losses Q can be obtained according to Poynting's theorem through the straightforward computation of the following ratio:

Q := S ω ε 0 ℑm(ε g ) E • E ds L ℜe{E 0 × H 0 } • n dl , (5.43)
The numerator in Eq. ( 5.43) clarifies losses in Watts by period of the considered grating and are computed by integrating the Joule effect losses density over the surface S of the lossy element. The denominator normalizes these losses to the incident power, i.e. the time-averaged incident Poynting vector flux across one period (a straight line L of length d in the superstrate parallel to Ox, whose normal oriented along decreasing values of y is denoted n).

Finally, combining Eqs. (5.42) and Eq. (5.43), a self consistency check of the whole numerical scheme consists in comparing the quantity B:

B := ∑ n T n + ∑ m R m + Q (5.44)
to unity. In Eq. (5.44), the summation indexed by n (resp. m) corresponds to the sum over the efficiencies of all transmitted (resp. reflected) propagative diffraction orders in the substrate (resp. superstrate). We give interpretations and concrete examples of such numerical energy balances over non trivial grating profiles in sections 5.2.3.2 and 5.2.3.3.

Numerical experiments

Numerical validation of the method

We can refer to [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF] in order to test the accuracy of our method. The studied grating is isotropic, since we lack numerical values in the literature in anisotropic cases. We compute the following problem (cf. Fig. 5.4), as described in [START_REF] Bao | Adaptive finite-element method for diffraction gratings[END_REF] and [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF]. The wavelength of the plane wave is set to 1 µm and is incoming with an angle of π/6 with respect to the normal to the grating. We present the R 0 efficiency (cf. Table 5.1) in both cases of polarization versus the mesh refinement. So we have a good agreement to the reference values, and the accuracy reached is independent from the polarization case. This grating has been studied by [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF] and is one of our points of reference [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF] 0.7342789 0.8484781

Maximum element size R TE 0 R TM 0 λ 0 /(4 √ ε) 0.
Tab. 5.1: Reflected efficiencies versus mesh refinement. Note that the efficiencies are properly computed (two significant digits) even for a rather coarse mesh.

Experiment set up based on existing materials

The method proposed in this section is adapted to z-anisotropic materials, such as transparent CaCO 3 [START_REF] Tayeb | Contribution à l'étude de la diffraction des ondes électromagnétiques par des réseaux. Réflexions sur les méthodes existantes et sur leur extension aux milieux anisotropes[END_REF], LiNbO 3 [START_REF] Ohkawa | Analysis of anisotropic dielectric grating diffraction using the finite-element method[END_REF] or Ni:YIG [START_REF] Kono | A novel finite-element method for nonreciprocal magneto-photonic crystal waveguides[END_REF] and lossy CoPt or CoPd [START_REF] Zhou | Dielectric tensor characterization for magneto-optical recording media[END_REF]. Let us now consider a trapezoidal (cf. Fig. 5.5) anisotropic grating made of aragonite (CaCO 3 ) deposited on an isotropic substratum (SiO 2 , ε SiO 2 = 2.25). Along the anisotropic crystal axis, its dielectric tensor can be written as follows [START_REF] Tayeb | Contribution à l'étude de la diffraction des ondes électromagnétiques par des réseaux. Réflexions sur les méthodes existantes et sur leur extension aux milieux anisotropes[END_REF]:

ε CaCO 3 =   2.843 0 0 0 2.341 0 0 0 2.829   and µ CaCO 3 =   µ 0 0 0 0 µ 0 0 0 0 µ 0   (5.45)
300 nm 500 nm 600 nm 600 nm Fig. 5.5: Diffractive element pattern. This element is made of aragonite for which the dielectric tensor is given by Eq. (5.46) and is deposited on a silica substrate with a period d = 600 nm. Now let's assume that the natural axis of our aragonite grating are rotated by 45 We shall here remind that our method remains strictly the same whatever the diffractive element geometry is. The 2D computational domain is bounded along the y-axis by the PMLs and along the x since we consider only one pseudo period. We propose to calculate the diffractive efficiencies at λ 0 = 633 nm in both polarization cases TE and TM, and for different incoming incidences (0 • , 20 • and 40 • ). Since both µ and ε are Hermitian, the whole incident energy is diffracted and the sum of theses efficiencies ought to be equal to the incident energy, which will stand for validation of our numerical calculation. Finally, the resulting bounded domain is meshed with a maximum mesh element side size of λ 0 /10 √ ε. Efficiencies are still post-processed in accordance with the calculation presented section 5.2.2.6. At normal incidence, the h field in the TE case (cf. Fig. 5.6b) is non symmetric whereas the e field in the TM case is (cf. Fig. 5.6a). This is illustrated by the obvious non-symmetry of T TE -1 and T TE 1 (cf. Table 5.2: 0.322510 versus 0.124722!), whereas

TM T -2 T -1 T 0 T 1 R -1 R 0 R 1 total 0 • - 0.
T -1 T 0 T 1 R -1 R 0 R 1 total 0 • - 0.
T TM -1 = T TM 1 = 0.20313. TM TE (a) ℜe{e} in V /m at θ 0 = 0 • (b) ℜe{h} in A/m at θ 0 = 0 • (c) ℜe{e} in V /m at θ 0 = 20 • (d) ℜe{h} in A/m at θ 0 = 20 • (e) ℜe{e} in V /m at θ 0 = 40 • (f) ℜe{h} in A/m at θ 0 = 40 •

A non trivial geometry

Since the beginning of this chapter, we have laid great stress upon the independence of the method towards the geometry of the pattern. But we have considered so far diffractive objects of simple trapezoidal section. Let us tackle a way more challenging case and see what this approach is made of. We can obtain an quite winding shape by extracting the contrast contour of an arbitrary image (see Fig. 5.7a-5.7b). The contour is approximated by a set of splines, and the resulting domain is finely meshed (Fig. 5.7c). Finally, as shown in Fig. 5.7b, the formed pattern (h g /λ 0 = 1.68), breathing in free space (ε substrat = 1), is supposed to be periodically repeated d/λ 0 = 1.26 on a plane ground of glass (ε SiO 2 = 2.25). The element is considered to be "made of" a lossy material of high optical index (ε marsu = 40 + 0.1 i). The response of this system to a incident s-polarized plane wave at oblic incidence (θ 0 = -30 • ) is finally calculated. The real part of the quasi-periodic total field is represented in Fig. 5.7d for three periods. Indeed, we do not have any tabulated data available to check our results. But what we do have is a pretty reliable consistency check through the computation of the energy balance described by Eqs. (5.42) and (5.43). As shown in Fig. 5.7e, we obtain at least 7 significative digits on the energetic values. The total balance of 1.00000019 is computed taking into account (i) values of the total field inside the diffractive elements, (ii) values of the diffracted field at altitudes spanning the whole (modeled) superstrate, (iii) values of the total field at altitudes spanning the entire (modeled) substrate. Finally, (iv) the calculated field u d 2 also nicely decays exponentially inside both PML. These four points allow us to check a posteriori the validity of the field everywhere in the computation cell.

Dealing with Wood anomalies using Adaptative PML

As we have noticed at the end of Section 5.2.2.4, PMLs based on "traditional coordinate stretching" are inefficient for periodic problems when dealing with grazing angles of diffracted orders, i.e. when the frequency is near a Wood's anomaly ( [START_REF] Wood | On a Remarkable Case of Uneven Distribution of Light in a Diffraction Grating Spectrum[END_REF][START_REF] Rayleigh | Note on the remarkable case of diffraction spectra described by Prof. Wood[END_REF]), leading to spurious reflexions and thus numerical pollution of the results. An important question in designing absorbing layers is thus the choice of their parameters: The PML thickness and the absorption coefficient. To this aim, adaptative formulations have already been set up, most of them employing a posteriori error estimate [START_REF] Chen | An adaptive perfectly matched layer technique for time-harmonic scattering problems[END_REF][START_REF] Bao | Adaptive finite-element method for diffraction gratings[END_REF][START_REF] Schädle | Domain decomposition method for maxwell's equations: Scattering off periodic structures[END_REF]. In this section, we propose Adaptative PMLs (APMLs) with a suitable coordinate stretching, depending both on incidence and grating parameters, capable of efficiently absorbing propagating waves with nearly grazing angles. This section is dedicated to the mathematical formulation used to determine PML parameters adapted to any diffraction orders. We provide at the end a numerical example of a dielectric slit grating showing the relevance of our approach in comparison with classical PMLs. As explained in Section 5.2.2.6, the diffracted field u d can be expanded as a Rayleigh expansion, i.e. into an infinite sum of propagating and evanescent plane waves called diffraction orders. As detailed at the end of Section 5.2.2.4, we are now in position to rewrite easily the expression of, say, a transmitted diffraction order into the substrate. Similar considerations also apply to the reflected orders in the top PML. Combining Eq. (5.32) and (5.40) lead to the expression u - n,s (y c ) of a transmitted propagative order inside the PML:

Skin depth of the PML

h g ε g (x, y), µ g (x, y) h + ε + , µ + ĥ+ ε+ (x, y), μ+ (x, y) h - ε -, µ - ĥ- ε-(x,
u - n,s (y c ) = u - n (y(y c )) = t n e -iβ - n [y t +ζ -(y c -y t )] .
The non oscillating part of this function is given by:

U - n (y) = t n exp (β ,- n ζ ,-+ β ,- n ζ ,-)y c , where β - n = β ,- n + iβ ,- n .
For a propagating order we have β ,- n > 0 and β ,- n = 0, while for an evanescent order β ,- n = 0 and β ,- n > 0. It is thus sufficient to take ζ ,-> 0 and ζ ,-> 0 to ensure the exponential decay to zero of the field inside the PML if it was of infinite extent. But, of course, for practical purposes, the thickness of the PML is finite and has to be suitably chosen. Two pitfalls must be avoided:

1. The PML thickness is chosen too small compared to the skin depth. As a consequence, the electromagnetic wave cannot be considered as vanishing: An incident electromagnetic "sees the bottom of the PML". In other words, this PML of finite thickness is no longer reflection-less.

2. The PML thickness is chosen much larger than the skin depth. In that case, a significant part of the PML is not useful, which gives rise to the resolution of linear systems of unnecessarily large dimensions.

Then remains to derive the skin depth, l - n , associated with the propagating order n. This characteristic length is defined as the depth below the PML at which the field falls to 1/e of its value near the surface:

U - n (y -l - n ) = U - n (y) e .
Finally, we find l - n = (β ,- n ζ ,-+ β ,- n ζ ,-) -1 and we define l -as the largest value among the l - n :

l -= max n∈Z l - n .
The height of the bottom PML region is set to ĥ-= 10l -.

Weakness of the classical PML for grazing diffracted angles

Let us consider the (bottom) PML adapted to the substrate. Similar conclusions will hold for the top PML. The efficiency of the classical PML fails for grazing diffracted angles, in other words when a given order appears/vanishes: this is the so-called Wood's anomaly, well known in the grating theory. In mathematical terms, there exists n 0 such that β - n 0 0. The skin depth of the PML then becomes very large. To compensate this, it is tempting to increase the value of ζ ,-, but it would lead to spurious numerical reflections due to an overdamping. For a fixed value of ĥ-, if ζ ,-is too weak, the absorption in the PMLs is insufficient and the wave is reflected on the outward boundary of the PML. To illustrate these typical behaviors (cf. Fig. 5.9), we compute the field diffracted by a grating with a rectangular cross section of height h g = 1.5 µm and width L g = 3 µm with ε g = 11.7, deposited on a substrate with permittivity ε -= 2.25. The structure is illuminated by a p-polarized plane wave of wavelength λ 0 = 10 µm and of angle of incidence θ 0 = 10 • in the air (ε + = 1). All materials are non magnetic (µ r = 1) and the periodicity of the grating is d = 4 µm. We set ĥ-= 10l - 0 and ζ ,-= 1.

Construction of an adaptative PML

To overcome the problems pointed out in the previous section, we propose a coordinate stretching that rigorously treats the problem of Wood's anomalies. The wavelengths "seen" by the system are very different depending on the order at stake:

• if the diffracted angle θ n is zero, the apparent wavelength λ 0 / cos θ n is simply the incident wavelength,

• if the diffracted angle is near ±π/2 (grazing angle), the apparent wavelength λ 0 / cos θ n is very large.

Thus if a classical PML is adapted to one diffracted order, it will not be for another, and vice versa. The idea behind the APML is to deal with each and every order when progressing in the absorbing medium.

Once again the development will be conducted only for the PML adapted to the substrate. We consider a real-valued coordinate mapping y d (y), the final complex-valued mapping is then y c (y) = ζ -y d (y), with the complex constant ζ -, with ζ ,-> 0 and ζ ,-> 0, accounting for the damping of the PML medium.

We begin with transforming the equation

β ± n 2 = k ± 2 -α ± n 2
, so that the function with integer argument n → β - n becomes a function with real argument continuously interpolated between the imposed integer values. Indeed, the geometric transformations associated to the PML has to be continuous and differentiable in order to compute its Jacobian. To that extent, we choose the parametrization:

α(y d ) = α 0 + 2π d y d λ 0 , (5.47)
so that the application β -defined by β -(y d ) 2 = k 2 0 ε -α(y d ) 2 is continuous. Thus, the propagation constant of the n th transmitted order is given by β - n = β -(nλ 0 ). The key idea is to combine the complex stretching with a real non uniform contraction (given by the continuous function y(y d ), Eq. (5.49)). This contraction is chosen in such a way that for each order n there is a depth y n d such that, around this depth, the apparent wavelength corresponding to the order in play is contracted to a value close to λ 0 . At that point of the PML, this order is perfectly absorbed thanks to the complex stretch. We thus eliminate first the orders with quasi normal diffracted angles at lowest depths up to grazing orders (near Wood's anomalies) which are absorbed at greater depths. In mathematical words, the translation of previous considerations on the real contraction can be expressed as:

exp [-iβ -(y d )y(y d )] = exp (-ik 0 y d ) (5.48)
The contraction y(y d ) is thus given by:

y(y d ) = k 0 y d β -(y d ) = y d ε --(sin θ 0 + y d /d) 2
(5.49)

The function y(y d ) has two poles, denoted

y d,± = d(± √ ε --sin θ 0 ). When y d,± = ±nλ 0 with n ∈ N , β -(y d,± ) = β -(±nλ 0 ) = β - ± = 0, i.e.
we are on a Wood's anomaly associated with the appearance/disappearance of the ±n th transmitted order. We now search for the nearest point to y * d,± associated with a Wood's anomaly, denoting: (5.51) Figure 5.10 shows an example of this coordinate mapping. Eventually, the complex stretch s y used in Eq. (5.29) is given by:

         n + / D + = min n + ∈N |y d,+ -n + λ 0 | n -/ D -= min
s y (y d ) = ζ -∂ ỹ ∂ y d (y d ).
(5.52)

Equipped with this mathematical formulation, we can tailor a layer that is doubly perfectly matched:

• to a given medium, which is the aim of the PML technique, through Eq. (5.27),

• to all diffraction orders, through the stretching coefficient s y , which depends on the characteristics of the incident wave and on opto-geometric parameters of the grating.

Numerical example

We now apply the method described in the preceding parts to design an adapted bottom PML for the same example as in part 5.2.4.2. The parameters are the same, and we choose the wavelength of the incident plane wave close to the Wood's anomaly related to the +1 transmitted order (λ 0 = 0.999y d,+ ).Moreover, we set the length of the PML ĥ-= 1.1y d,+ and choose absorption coefficients ζ + = ζ -= 1 + i. For both cases (PML and APML), parameters are alike, the only difference being the complex stretch s y .

The field maps of the norm of H z , E x and E y are plotted in logarithmic scale on Fig. 5.11, for the case of a classical PML and our APML. We can observe that the field H z that is effectively computed is clearly damped in the bottom APML (leftmost on Fig. 5.11(b)) whereas it is not in the standard case (leftmost on Fig. 5.11(a)), causing spurious reflections on the outer boundary. The fields E x and E y are deduced from H z thanks to Maxwell's equations. The high values of E y at the tip of the APML (rightmost on Fig. 5.11(b)) are due to very high values of the optical equivalent properties of the APML medium (due to high values of s y ), which does not affect the accuracy of the computed field within the domain of interest. Another feature of our approach is that it efficiently absorbs the grazing diffraction order, as illustrated on Fig. 5.12: the +1 transmitted order does not decrease in the standard PML (blue solid line), and reaches a high value at y = -ĥ-, whereas the same order tends to zero as y → -ĥin the case of the adapted PML (blue dashed line). To further validate the accuracy of the method, we compare the diffraction efficiencies computed by our FEM formulation with PML and APML to those obtained by another method. We choose the Rigorous Coupled Wave Analysis (RCWA), also known as the Fourier Modal Method (FMM, [START_REF] Li | New formulation of the fourier modal method for crossed surface-relief gratings[END_REF]). For the chosen parameters, only the 0 th order is propagative in reflexion and the orders -1, 0 and +1 are non evanescent in transmission. We can also check the energy balance B = R 0 + T -1 + T 0 + T +1 since there is no lossy medium in our example. Results are reported in Table 5.3, and show a good agreement of the FEM with APML with the results from RCWA. On the contrary, if classical PML are used, the diffraction efficiencies are less accurate 5.12: Modulus of the u n for the three propagating orders with adapted (dashed lines) and classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except for the grazing one (n = 1) as expected. This drawback is bypassed when using the adaptative PML.

compared to those computed with RCWA. Checking the energy balance leads the same conclusions: the numerical result is perturbed by the reflection of the waves at the end of the PML if it is not adapted to the situation of nearly grazing diffracted orders. Eventually, to illustrate the behavior of the adaptative PML when the incident wavelength gets closer to a given Wood's anomaly, we computed the mean value of the norm of H z along the outer boundary of the bottom PML γ = |H z (-ĥ-)| x , when λ 0 = (1 + 10 -n )y d,+ and λ 0 = (1 -10 -n )y d,+ , for n = 1, 2, ...10. The results are shown in Fig. 5.13. As the wavelength gets closer to y d,+ , γ first increases but for n > 3, it decreases exponentially. However, in all cases, the value of γ remains small enough to ensure the efficiency of the PMLs.

Concluding remarks

A novel FEM formulation was adapted to the analysis of z-anisotropic gratings relying on a rigorous treatment of the plane wave sources problem through an equivalent radiation problem with localized sources. The developed approach presents the advantage of being very general in the sense that it is applicable to every conceivable grating geometry.

Numerical experiments based on existing materials at normal and oblique incidences in both TE and TM cases showed the efficiency and the accuracy of our method. We demonstrated we could generate strongly imbalanced symmetric propagative orders in the TE polarization case and at normal incidence with an aragonite grating on a silica substratum.

We also introduced the adaptative PML for grazing incidences configurations. It based on a complex-valued coordinate stretching that deals with grazing diffracted orders, yielding an efficient absorption of the field inside the PML. We provided an example in the TM polarization case (but similar results hold for the TE case), illustrating the efficiency of our method. The value of the magnetic field on the outward boundary of the PML remains small enough to consider there is no spurious reflection. The formulation is used with the FEM but can be applied to others numerical methods. Moreover, the generalization to the vectorial three-dimensional case is straightforward: the recipes given in this last section do work irrespective of the dimension and whether the problem is vectorial.

In the next section, the scalar formulation adapted to mono-dimensional gratings is extended to the the most general case of bi-dimensional grating embedded in an arbitrary multilayered dielectric stack with arbitrary incidence.

Diffraction by arbitrary crossed-gratings : a vector Finite Element formulation

Introduction

In this section, we extend the method detailed in Sec. 5.2 to the most general case of vector diffraction by an arbitrary crossed gratings. The main advantage of the Finite Element Method lies in its native ability to handle unstructured meshes, resulting in a build-in accurate discretization of oblique edges. Consequently, our approach remains independent of the shape of the diffractive element, whereas other methods require heavy adjustments depending on whether the geometry of the groove region presents oblique edges (e.g. RCWA [START_REF] Popov | Staircase approximation validity for arbitrary-shaped gratings[END_REF], FDTD. . . ). In this section, for the sake of clarity, we recall again the rigorous procedure allowing to deal with the issue of the plane wave sources through an equivalence of the diffraction problem with a radiation one whose sources are localized inside the diffractive element itself, as already proposed in Sec. 5.2 [START_REF] Demésy | The finite element method as applied to the diffraction by an anisotropic grating[END_REF][START_REF] Demésy | Finite element method as applied to the study of gratings embedded in complementary metal-oxide semiconductor image sensors[END_REF].

This approach combined with the use of second order edge elements allowed us to retrieve with a good accuracy the few numerical academic examples found in the literature. Furthermore, we provide a new reference case combining major difficulties such as a non trivial toroidal geometry together with strong losses and a high permittivity contrast. Finally, we discuss computation time and convergence as a function of the mesh refinement as well as the choice of the direct solver.

Theoretical developments

Set up of the problem and notations

We denote by x, y and z the unit vectors of the axes of an orthogonal coordinate system Oxyz. We only deal with time-harmonic fields; consequently, electric and magnetic fields are represented by the complex vector fields E and H, with a time dependance in exp(-i ω t). Note that incident light is now propagating along the z-axis, whereas y-axis was used in the 2D case.

Besides, in this section, for the sake of simplicity, the materials are assumed to be isotropic and therefore are optically characterized by their relative permittivity ε and relative permeability µ (note that the inverse of relative permeabilities are denoted here ν). It is of importance to note that lossy materials can be studied, the relative permittivity and relative permeability being represented by complex valued functions. The crossed-gratings we are dealing with can be split into the following regions as suggested in Fig. 5.14:

• The superstrate (z > z 0 ) is supposed to be homogeneous, isotropic and lossless, and therefore characterized by its relative permittivity ε + and its relative permeability µ + (= 1/ν + ) and we denote k

+ := k 0 ε + µ + , where k 0 := ω/c,
• The multilayered stack (z N < z < z 0 ) is made of N layers which are supposed to be homogeneous and isotropic, and therefore characterized by their relative permittivity ε n , their relative permeability µ n (= 1/ν n ) and their thickness e n . We denote k n := k 0 √ ε n µ n for n integer between 1 and N.

• The groove region (z g < z < z g-1 ), which is embedded in the layer indexed g (ε g , µ g ) of the previously described domain, is heterogeneous. Moreover the method does work irrespective of whether the diffractive elements are homogeneous: The permittivity and permeability can vary continuously (gradient index gratings) or discontinuously (step index gratings). This region is thus characterized by the scalar fields ε g (x, y, z) and µ g (x, y, z)(= 1/ν g (x, y, z)). The groove periodicity along the x-axis, respectively (resp.) y-axis, is denoted d x , resp. d y , in the sequel.

• The substrate (z < z N ) is supposed to be homogeneous and isotropic and therefore characterized by its relative permittivity ε -and its relative permeability µ -(= 1/ν -) and we denote k

-:= k 0 ε -µ -,
Let us emphasize the fact that the method principles remain unchanged in the case of several diffractive patterns made of distinct geometry and/or material. The incident field on this structure is denoted:

E inc = A e 0 exp(i k + p • r) (5.53) with k + =   α 0 β 0 γ 0   = k +   -sin θ 0 cos ϕ 0 -sin θ 0 sin ϕ 0 -cos θ 0   (5.54)
and

A e 0 =   E 0 x E 0 y E 0 z   = A e  
cos ψ 0 cos θ 0 cos ϕ 0sin ψ 0 sin ϕ 0 cos ψ 0 cos θ 0 sin ϕ 0 + sin ψ 0 cos ϕ 0 cos ψ 0 sin θ 0   , (5.55)

where ϕ 0 ∈ [0, 2π], θ 0 ∈ [0, π/2] and ψ 0 ∈ [0, π] (polarization angle).

We recall here the diffraction problem: finding the solution of Maxwell equations in harmonic regime i.e. the unique solution (E, H) of:

curl E = i ω µ 0 µ H (5.56a) curl H = -i ω ε 0 ε E (5.56b)
such that the diffracted field satisfies the so-called Outgoing Waves Condition (OWC [START_REF] Zolla | Method of fictitious sources as applied to the electromagnetic diffraction of a plane wave by a grating in conical diffraction mounts[END_REF] ) and where E and H are quasi-bi-periodic functions with respect to x and y coordinates.

One can choose to calculate arbitrarily E, since H can be deduced from Eq. (5.56a). The diffraction problem amounts to looking for the unique solution E of the so-called vectorial Helmholtz propagation equation, deduced from Eqs. (5.56a,5.56b):

M ε,ν := -curl (ν curl E) + k 2 0 ε E = 0 (5.57)
such that the diffracted field satisfies an OWC and where E is a quasi-bi-periodic function with respect to x and y coordinates.

From a diffraction problem to a radiative one with localized sources

According to Fig. 5.14, the scalar relative permittivity ε and inverse permeability ν fields associated to the studied diffractive structure can be written using complex-valued functions defined by part and taking into account the notations adopted in Sec. 5.3.2.1:

υ(x, y, z) :=            υ + for z > z 0 υ n for z n-1 > z > z n with 1 ≤ n < g υ g (x, y, z) for z g-1 > z > z g υ n for z n-1 > z > z n with g < n ≤ N υ - for z < z N (5.58)
with υ = {ε, ν} , z 0 = 0 and z n = -∑ n l=1 e l for 1 ≤ n ≤ N. It is now convenient to introduce two functions defined by part ε 1 and ν 1 corresponding to the associated multilayered case (i.e. the same stack without any diffractive element) constant over Ox and Oy:

υ 1 (x, y, z) :=    υ + for z > 0 υ n for z n-1 > z > z n with 1 ≤ n ≤ N υ -for z < z N (5.59)
with υ = {ε, ν}.

We denote by E 0 the restriction of E inc to the superstrate region:

E 0 := E inc for z > z 0 0 for z ≤ z 0 (5.60)
We are now in a position to define more explicitly the vector diffraction problem that we are dealing with in this section. It amounts to looking for the unique vector field E solution of:

M ε,ν (E) = 0 such that E d := E -E 0 satisfies an OWC. ( 5.61)
In order to reduce this diffraction problem to a radiation one, an intermediary vector field denoted E 1 is necessary and is defined as the unique solution of:

M ε 1 ,ν 1 (E 1 ) = 0 such that E d 1 := E 1 -E 0 satisfies an OWC.
(5.62)

The vector field E 1 corresponds to an ancillary problem associated to the general vectorial case of a multilayered stack which can be calculated independently. This general calculation is seldom treated in the literature, we present a development in Appendix. Thus E 1 is from now on considered as a known vector field. It is now apropos to introduce the unknown vector field E d 2 , simply defined as the difference between E and E 1 , which can finally be calculated thanks to the FEM and:

E d 2 := E -E 1 = E d -E d 1 .
(5.63)

It is of importance to note that the presence of the superscript d is not fortuitous: As a difference between two diffracted fields (Eq. (5.63), E d 2 satisfies an OWC which is of prime importance in our formulation. By taking into account these new definitions, Eq. (5.61) can be written:

M ε,ν (E d 2 ) = -M ε,ν (E 1 ) , (5.64)
where the right-hand member is a vector field which can be interpreted as a known vectorial source term -S 1 (x, y, z) whose support is localized inside the diffractive element itself. To prove it, let us introduce the null term defined in Eq. ( 5.62) and make the use of the linearity of M , which leads to:

S 1 := M ε,ν (E 1 ) = M ε,ν (E 1 ) -M ε 1 ,ν 1 (E 1 ) =0 = M ε-ε 1 ,ν-ν 1 (E 1 ) .
(5.65)

Quasi-periodicity and weak formulation

The weak form is obtained by multiplying scalarly Eq. ( 5.61) by weighted vectors E chosen among the ensemble of quasi-bi-periodic vector fields of

L 2 (curl) (denoted L 2 (curl, (d x , d y ), k)) in Ω: R ε,ν (E, E ) = Ω -curl (ν curl E) • E + k 2 0 ε E • E dΩ (5.66)
Integrating by part Eq. (5.66) and making the use of the Green-Ostrogradsky theorem lead to:

R ε,ν (E, E ) = Ω -ν curl E • curl E + k 2 0 ε E • E dΩ - ∂ Ω (n × (ν curl E)) • E dS (5.67)
where n refers to the exterior unit vector normal to the surface ∂ Ω enclosing Ω.

The first term of this sum concerns the volume behavior of the unknown vector field whereas the right-hand term can be used to set boundary conditions (Dirichlet, Neumann or so called quasi-periodic Bloch-Floquet conditions).

The solution E d 2 of the weak form associated to the diffraction problem, expressed in its previously defined equivalent radiative form at Eq. (5.64), is the element of L 2 (curl, (d x , d y ), k) such that:

∀E ∈ L 2 (curl, d x , d y , k), R ε,ν (E d 2 , E ) = -R ε-ε 1 ,ν-ν 1 (E 1 , E ) . (5.68)
In order to rigorously truncate the computation a set of Bloch boundary conditions are imposed on the pair of planes defined by (y = -d y /2, y = d y /2) and (x = -d x /2, x = d x /2). One can refer to [START_REF] Nicolet | Modelling of electromagnetic waves in periodic media with finite elements[END_REF] for a detailed implementation of Bloch conditions adapted to the FEM. A set of Perfectly Matched Layers are used in order to truncate the substrate and the superstrate along z axis (see [START_REF] Ould Agha | On the use of PML for the computation of leaky modes : an application to gradient index MOF[END_REF] for practical implementation of PML adapted to the FEM). Since the proposed unknown E d 2 is quasi-bi-periodic and satisfies an OWC, this set of boundary conditions is perfectly reasonable: E d 2 is radiated from the diffractive element towards the infinite regions of the problem and decays exponentially inside the PMLs along z axis. The total field associated to the diffraction problem E is deduced at once from Eq. (5.63).

Edge or Whitney 1-form second order elements

In the vectorial case, edge elements (or Whitney forms) make a much more relevant choice [START_REF] Dular | A discrete sequence associated with mixed finite elements and itsgauge condition for vector potentials[END_REF] than nodal elements. Note that a lot of work (see for instance [START_REF] Ingelstrom | A new set of H (curl)-conforming hierarchical basis functions for tetrahedral meshes[END_REF]) has been done on higher order edge elements since their introduction by Bossavit [START_REF] Bossavit | Edge-elements for scattering problems[END_REF]. These elements are suitable to the representation of vector fields such as E d 2 , by letting their normal component be discontinuous and imposing the continuity of their tangential components. Instead of linking the Degrees Of Freedom (DOF) of the final algebraic system to the nodes of the mesh, the DOF associated to edges (resp. faces) elements are the circulations (resp. flux) of the unknown vector field along (resp. across) its edges (resp. faces).

Let us consider the computation cell Ω together with its exterior boundary ∂ Ω. This volume is sampled in a finite number of tetrahedron according to the following rules: Two distinct tetrahedrons have to either share a node, an edge or a face or have no contact. Let us denote by T the set of tetrahedrons, F the set of faces, E the set of edges and N the set of nodes. In the sequel, one will refers to the node n = {i}, the edge e = {i, j}, the face f = {i, j, k} and the tetrahedron t = {i, j, k, l}. Twelve DOF (two for each of the six edges of a tetrahedron) are classically derived from line integral of weighted projection of the field E d 2 on each oriented edge e = {i, j}:

       ϑ i j = j i E d 2 • t i j λ i dl ϑ ji = i j E d 2 • t ji λ j dl , (5.69)
where t i j is the unit vector and λ i , the barycentric coordinate of node i, is the chosen weight function.

According to Yioultsis et al. [START_REF] Yioultsis | The Mystery and Magic of Whitney Elements -An Insight in their Properties and Construction[END_REF], a judicious choice for the remaining DOF is to make the use of a tangential projection of the 1-form E d 2 on the face f = {i, j, k}.

     ϑ i jk = f E d 2 × n + i jk • grad λ j ds ϑ ik j = f E d 2 × n - i jk • grad λ k ds . ( 5.70)
The expressions for the shape functions, or basis vectors, of the second order 1-form Whitney element are given by:

w i j = (8 λ 2 i -4 λ i ) grad λ j + (-8 λ i λ j + 2 λ j ) grad λ i w i jk = 16 λ i λ j grad λ k -8 λ j λ k grad λ i -8 λ k λ i grad λ j . (5.71)
This choice of shape function ensures [START_REF] Yioultsis | Multiparametric vector finite elements: a systematic approach to the construction of three-dimensional, higher order, tangential vector shape functions[END_REF] the following fundamental property: every degree of freedom associated with a shape function should be zero for any other shape function. Finally, an approximation of the unknown E d 2 projected on the shape functions of the mesh m (E d,m 2 ) can be derived:

E d,m 2 = ∑ e∈E ϑ e w e + ∑ f ∈F ϑ f w f . (5.72)
Weight functions E (c.f. Eq. ( 5.68) are chosen in the same space than the unknown E d 2 , L 2 (curl, (d x , d y ), k). According to the Galerkin formulation, this choice is made so that their restriction to one bi-period belongs to the set of shape functions mentioned above. Inserting the decomposition of E d 2 of Eq. (5.72) in Eq. ( 5.68) leads to the final algebraic system which is solved, in the following numerical examples, thanks to direct solvers.

Energetic considerations: Diffraction efficiencies and losses

Contrarily to modal methods based on the determination of Rayleigh coefficients, the rough results of the FEM are three complex components of the vector field E d interpolated over the mesh of the computation cell. Diffraction efficiencies are deduced from this field maps as follows.

As a difference between two quasi-periodic vector fields (see Eq. ( 5.61)), E d is quasi-biperiodic and its components can be expanded as a double Rayleigh sum: Eq. ( 5.74) allows to evaluate numerically e x,c n,m (resp. e x,p n,m ) by double trapezoidal integration of a slice of the complex component E d

x at an altitude z c fixed in the superstrate (resp. substrate).

It is well known that the mere trapezoidal integration method is very efficient for smooth and periodic functions (integration on one period). The same holds for E d y and E d z components as well as their coefficients e y,{c,p} n,m and e z,{c,p} n,m . The dimensionless expression of the efficiency of each reflected and transmitted (n, m) order [START_REF] Noponen | Eigenmode method for electromagnetic synthesis of diffractive elements with three-dimensional profiles[END_REF] is deduced from Eqs. (5.75,5.76): z. Furthermore, normalized losses Q can be obtained through the computation of the following ratio:

       R n,m = 1 |A e | 2
Q = V 1 2 ω ε 0 ℑm(ε g ) E • E dV S 1 2 ℜe{E 0 × H 0 } • n dS . (5.78)
The numerator in Eq. ( 5.78) clarifies losses in watts by bi-period of the considered crossedgrating and are computed by integrating the Joule effect losses density over the volume V of the lossy element. The denominator normalizes these losses to the incident power, i.e. the timeaveraged incident Poynting vector flux across one bi-period (a rectangular surface S of area d x d y in the superstrate parallel to Oxy, whose normal oriented along decreasing values of z is denoted n). Since E 0 is nothing but the plane wave defined at Eqs. (5.54,5.55), this last term is equal to (A 2 e ε 0 /µ 0 d x d y )/(2 cos(θ 0 )). Volumes and normal to surfaces being explicitly defined, normalized losses losses Q are quickly computed once E determined and interpolated between mesh nodes.

Finally, the accuracy and self-consistency of the whole calculation can be evaluated by summing the real part of transmitted and reflected efficiencies (n, m) to normalized losses:

Q + ∑ (n,m)∈Z 2 ℜe{R n,m } + ∑ (n,m)∈Z 2
ℜe{T n,m } , quantity to be compared to 1. The sole diffraction orders taken into account in this conservation criterium correspond to propagative orders whose efficiencies have a non-null real part. Indeed, diffraction efficiencies of evanescent orders, corresponding to pure imaginary values of γ ± n,m for higher values of (n, m) (see Eq. (5.75)) are also pure imaginary values as it appears clearly in Eq. (5.77). Numerical illustrations of such global energy balances are presented in the next section.

Accuracy and convergence

Classical crossed gratings

There are only a few references in the literature containing numerical examples. For each of them, the problem only consists of three regions (superstrate, grooves and substrate) as summed up on Figure 5.16. For the four selected cases, among six found in the literature, published results are compared to ones given by our formulation of the FEM. Moreover, in each case, a satisfying global energy balance is detailed. Finally a new validation case combining all the difficulties encountered when modeling crossed-gratings is proposed: A non-trivial geometry for the diffractive pattern (a torus), made of an arbitrary lossy material leading to a large step of index and illuminated by a plane wave with an oblique incidence. Convergence of the FEM calculation as well as computation time will be discussed in Sec. 5.3.4.2.

Checkerboard grating

In this example worked out by L. Li [START_REF] Li | New formulation of the fourier modal method for crossed surface-relief gratings[END_REF], the diffractive element is a rectangular parallelepiped as shown Fig. 5.17a and the grating parameter highlighted in Fig. 5.16 are the following:

ϕ 0 = θ 0 = 0 • , ψ 0 = 45 • , d x = d y = 5 λ 0 √ 2/4, h = λ 0 , ε + = ε g = 2.
25 and ε -= ε g = 1. Our formulation of the FEM shows good agreement with the Fourier Modal Method developed by L. Li ([27], 1997) since the maximal relative difference between the array of values presented in Table 5.4 remains lower than 10 -3 . Moreover, the sum of the efficiencies of propagative orders given by the FEM is very close to 1 in spite of the addition of all errors of determination upon the efficiencies. Given in [START_REF] Derrick | Crossed gratings: A theory and its applications[END_REF] [40]

[41]

[42] FEM R -1,0 0.00254 0.00207 0.00246 0.00249 0.00251 R 0,0 0.01984 0.01928 0.01951 0.01963 0.01938 T -1,-1 0.00092 0.00081 0.00086 0.00086 0.00087 T 0,-1 0.00704 0.00767 0.00679 0.00677 0.00692 T -1,0 0.00303 0.00370 0.00294 0.00294 0.00299 T 0,0 0.96219 0.96316 0.96472 0.96448 0.96447 T 1,0 0.00299 0.00332 0.00280 0.00282 0.00290 TOTAL 0.99855 1.00001 1.00008 0.99999 1.00004 Tab. 5.5: Comparison with the results given in [START_REF] Derrick | Crossed gratings: A theory and its applications[END_REF][START_REF] Greffet | Diffraction of electromagnetic waves by crossed gratings: a series solution[END_REF][START_REF] Bräuer | Electromagnetic diffraction analysis of two-dimensional gratings[END_REF][START_REF] Granet | Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: Application to multilayer crossed gratings[END_REF].

good agreement with ones of the C method [START_REF] Derrick | Crossed gratings: A theory and its applications[END_REF][START_REF] Granet | Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: Application to multilayer crossed gratings[END_REF], the Rayleigh method [START_REF] Greffet | Diffraction of electromagnetic waves by crossed gratings: a series solution[END_REF] and the RCWA [START_REF] Bräuer | Electromagnetic diffraction analysis of two-dimensional gratings[END_REF]. Note that, in this case, some edges of the diffractive element are oblique.

Bi-sinusoidal grating In this example worked out by Bruno et al. [START_REF] Bruno | Numerical solution of diffraction problems: a method of variation of boundaries. III. doubly periodic gratings[END_REF], the surface of the grating is bi-sinusoidal (see Fig. surface, the bi-sinusoid was first sampled (15 × 15 points), then converted to a 3D file format. This sampling can account for the slight differences with the results obtained using the method of variation of boundaries developed by [START_REF] Bruno | Numerical solution of diffraction problems: a method of variation of boundaries. III. doubly periodic gratings[END_REF].

Circular apertures in a lossy layer In this example worked out by Schuster et al. [START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF], the diffractive element is a circular aperture in a lossy layer as shown Fig. 5.20a and the grating parameter highlighted in Fig. 5.16 are the following: λ 0 = 500 nm, ϕ 0 = θ 0 = 0 • , ε + = ε g = 1, ε g = 0.8125 + 5.2500 i and ε -= 2.25. In this lossy case, results obtained with the FEM show good agreement with the ones obtained with the FMM [START_REF] Li | New formulation of the fourier modal method for crossed surface-relief gratings[END_REF], the differential method [START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF][START_REF] Arnaud | Diffraction et diffusion de la lumière : modélisation tridimensionnelle et application à la métrologie de la microélectronique et aux techniques d'imagerie sélective en milieu diffusant[END_REF] and the RCWA [START_REF] Moharam | Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings[END_REF]. Joule losses inside the diffractive element can be easily calculated, which allows to provide a global energy balance for this configuration. Finally, the convergence of the value R 0,0 as a function of the mesh refinement will be examined.

Lossy tori grating

We finally propose a new test case for crossed-grating numerical methods. The major difficulty of this case lies both in the non trivial geometry (see Fig. 5.21a) of the diffractive object and in the fact that it is made of a material chosen so that losses are optimal inside it. The grating parameters highlighted in Fig. 5. [START_REF] Helluy | Intégrations numériques d'ordre élevé de fonctions régulières ou singulières sur un intervalle[END_REF] Tab. 5.8 illustrates the independence of our method towards the geometry of the diffractive element. ε g is chosen so that the skin depth has the same order of magnitude as b, which maximizes losses. Note that energy balances remain very accurate at normal and oblique incidence, in spite of both the non-triviality of the geometry and the strong losses.

Convergence and computation time

Convergence as a function of mesh refinement When using modal methods such as the RCWA or the differential method, based on the calculation of Rayleigh coefficients, a number proportional to N R have to be be determined a priori. Then, the unknown diffracted field is expanded as a Fourier serie, injected under this form in Maxwell equations, which annihilates xand y-dependencies. This leads to a system of coupled partial differential equations whose coefficients can structured in a matrix formalism. The resulting matrix is sometimes directly invertible (RCWA) depending on whether the geometry allows to suppress the z-dependance, which makes this method adapted to diffractive elements with vertically (or decomposed in staircase functions) shaped edge. In some other cases, one has to make the use of integral methods in order to solve the system, as in the pyramidal case for instance, which leads to the socalled differential method. The diffracted field map can be deduced from these coefficients. If the grating configuration only calls for a few propagative orders and if the field inside the groove region is not the main information sought for, these two close methods allow to determine the repartition of the incident energy quickly. However, if the field inside the groove region is the main piece of information, it is advisable to calculate many Rayleigh coefficients corresponding to evanescent waves which increases the computation time as (N R ) 3 or even (N R ) 4 .

FEM relies on the direct calculation of the vectorial components of the complex field. Rayleigh coefficients are determined a posteriori. The parameter limiting the computation time is the number of tetrahedral elements along which the computational domain is split up. We suppose that it is necessary to calculate at least two or three points (or mesh nodes) per period of the field (λ 0 / ℜe{ε}). Figure 5.22 shows the convergence of the efficiency R 0,0 (circular apertures case, see Fig. 5.20a) as a function of the mesh refinement characterized by the parameter N M : The maximum size of each element is set to λ 0 /(N M ℜe{ε}). It is of interest to note that even if N M < 3 the FEM still gives pertinent diffraction efficiencies: R 0,0 = 0.2334 for N M = 1 and R 0,0 = 0.2331 for N M = 2. The Galerkin method (see Eq. (5.67)) corresponds to a minimization of the error (between the exact solution and the approximation) with respect to a norm that can be physically interpreted in terms of energy-related quantities. Therefore, the finite element methods usually provide energy-related quantities that are more accurate than the local values of the fields themselves.

Computation time All the calculations were performed on a server equipped with 8 dual core Itanium1 processors and 256Go of RAM. Tetrahedral quadratic edge elements were used together with the direct solver PARDISO. Among different direct solvers adapted to sparse matrix algebra (UMFPACK, SPOOLES and PARDISO), PARDISO turned out to be the less time-consuming one as shown in Tab 5.9.

Solver

Computation time for 41720 DOF Computation time for 205198 DOF SPOOLES 15 mn 32 s 14 h44 mn UMFPACK 2 mn 07 s 1 h12 mn PARDISO 57 s 16 mn

Tab. 5.9: Computation time variations from solver to solver. It is of importance to note that for values of N M lower than 3, the problem can be solved in less than a minute on a standard laptop (4Go RAM, 2 × 2 GHz) with 3 significant digits on the diffraction efficiencies. This accuracy is more than sufficient in numerous experimental cases. Furthermore, as far as integrated values are at stake, relatively coarse meshes (N M ≈ 1) can be used trustfully, authorizing fast geometric, spectral or polarization studies.

Nowadays, the efficiency of the numerical algorithms for sparse matrix algebra together with the available power of computers and the fact that the problem reduces to a basic cell with a size of a small number of wavelengths make the finite element problem very tractable as proved here.

Concluding remarks

In this chapter, we demonstrate a general formulation of the FEM allowing to calculate the diffraction efficiencies from the electromagnetic field diffracted by arbitrarily shaped gratings embedded in a multilayered stack lightened by a plane wave of arbitrary incidence and polarization angle. It relies on a rigorous treatment of the plane wave sources problem through an equivalent radiation problem with localized sources. Bloch conditions and a new dedicated PML have been implemented in order to rigorously truncate the computational domain.

The principles of the method were discussed in detail for mono-dimensional gratings in TE/TM polarization cases (2D or scalar case) in a first part, and for the most general bidimensional or crossed gratings (3D or vector case) in a second part. Note that the very same concepts could be applied to the intermediate case of mono-dimensional gratings enlighten by an arbitrary incident plane wave (so-called conical case). The reader will find detail about the element basis relevant to this case in [START_REF] Nicolet | Modelling of electromagnetic waves in periodic media with finite elements[END_REF].

The main advantage of this formulation is its complete generality with respect to the studied geometries and the material properties, as illustrated with the lossy tori grating nontrivial case. Its principle remains independent of both the number of diffractive elements by period and number of stack layers. Its flexibility allowed us to retrieve with accuracy the few numerical academic examples found in the literature and established with independent methods. The remarkable accuracy observed in the case of coarse meshes, makes it a fast tool for the design and optimization of diffractive optical components (e.g. reflection and transmission filters, polarizers, beam shapers, pulse compression gratings. . . ). The complete independence of the presented approach towards both the geometry and the isotropic constituent materials of the diffractive elements makes it a handy and powerful tool for the study of metamaterials, finitesize photonic crystals, periodic plasmonic structures. . . The method described in this chapter has already been successfully applied to various problems, from homogenization theory [START_REF] Cabuz | Homogenization of nonlocal wire metamaterial via a renormalization approach[END_REF] or transformation optics [START_REF] Dupont | Revolution analysis of three-dimensional arbitrary cloaks[END_REF] to more applied concerns as the modeling of complex CMOS nanophotonic devices [START_REF] Demésy | Finite element method as applied to the study of gratings embedded in complementary metal-oxide semiconductor image sensors[END_REF] or ultra-thin new generation solar cells [START_REF] Demésy | Solar energy trapping with modulated silicon nanowire photonic crystals[END_REF].

5.A APPENDIX

This appendix is dedicated to the determination of the vector electric field in a dielectric stack enlightened by a plane wave of arbitrary polarization and incidence angle. This calculation, abundantly treated in the 2D scalar case, is generally not presented in the literature since, as far as isotropic cases are concerned, it is possible to project the general vectorial case on the two reference TE and TM cases. However, the presented formulation can be extended to a fully anisotropic case for which this TE/TM decoupling is no longer valid and the three components of the field have to be calculated as follows.

Let us consider the ancillary problem mentioned in Sec. 5.3.2.2, i.e. a dielectric stack made of N homogeneous, isotropic, lossy layers characterized by there relative permittivity denoted ε j and their thickness e j . This stack is deposited on a homogeneous, isotropic, possibly lossy substrate characterized by its relative permittivity denoted ε N+1 = ε -. The superstrate is air and its relative permittivity is denoted ε + = 1. Finally, we denote by z j the altitude of the interface between the j th and j + 1 th layers. The restriction of the incident field E inc to the superstrate region is denoted E 0 . The problem amounts to looking for (E 1 ,H 1 ) satisfying Maxwell equations in harmonic regime (see Eqs. (5.56a,5.56b)).

Across the interface z = z j

By projection on the main axis of the vectorial Helmholtz propagation equation (Eq. (5.57)), the total electric field inside the j th layer can be written as the sum of a propagative and a counter-propagative plane waves: What follows consists in writing the continuity of the tangential components of (E 1 , H 1 ) across the interface z = z j , i.e. the continuity of the vector field Ψ defined by:

Ψ =     E x 1 E y 1 i H x 1 i H y 1     .
(5.82)

The continuity of Ψ along Oz together with its analytical expression inside the j th and j + 1 th layers allows to establish a recurrence relation for the interface z = z j . Then, by projection of Eqs. (5.56a,5.56b) on Ox,Oy and Oz: Φ j+1 (z j )

   i β 0 H z 1 - ∂ H y 1 dz ∂ H x 1 dz -i α 0 H z 1 i α 0 H y 1 -i β 0 H x 1    = -i ω ε   E x 1 E y 1 E z 1   (5.83) and    i β 0 E z 1 - ∂ E y 1 ∂ z ∂ E x 1 ∂ z -i α 0 E z 1 i α 0 E y 1 -i β 0 E x 1    = i ω µ   H x
(5.90) Thanks to Eq. (5.90) and Eq. (5.89), a recurrence relation can be formulated for the analytical expression of E 1 in each layer: Φ j+1 (z j+1 ) = T j+1 Π -1 j+1 Π j Φ j (z j ) (5.91)

Fig. 5 . 1 :

 51 Fig. 5.1: Sketch and notations of the grating studied in this section.

Fig. 5 . 2 :

 52 Fig. 5.2: Quasi-periodicity of the field and sample of a d-periodic mesh.

  sc n (x c , y c ) := u n (x(x c ), y(y c )) = e iαx c e -iβ - n (y b +ζ -(y c -y b ))

Fig. 5 . 3 :

 53 Fig. 5.3: Example of computation of the radiated field u d 2 (TM case).

Fig. 5 . 4 :

 54 Fig.5.4: Rectangular groove grating: This pattern is repeatedly set up with a period d = 1 µm. This grating has been studied by[START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF] and is one of our points of reference

Fig. 5 . 6 :

 56 Fig. 5.6: Real part of the total calculated field depending on θ 0 and the polarization case

Fig. 5 . 7 :

 57 Fig. 5.7: (a) Initial contrasted image. (b) Proposed set up. (c) Sample mesh. (d) ℜe{E z } in V /m. (e) Energy balance of the problem.

Fig. 5 . 8 :

 58 Fig. 5.8: The basic cell used for the FEM computation of the diffracted field u d 2 .

Fig. 5 . 9 :

 59 Fig. 5.9: Zero th transmitted order by a grating with a rectangular cross section (see parameters in text, part 5.2.4.2) for different values of ζ ,-: blue line, ζ ,-= 1, correct damping; green line, ζ ,-= 0.1, underdamping; red line, ζ ,-= 20, overdamping.

  n -∈N |y d,-+ n -λ 0 | . In a second step, we look for the point y 0 d = n λ 0 such that: n / D = min n ∈{n + ,n -} (D + , D -) . (5.50) To avoid the singular behaviour at y d = y d,± , we continue the graph of the function y d (y) by a straight line tangent at y 0 d , which equation is t 0 (y d ) = s(y 0 d )(y dy 0 d ) + y(y 0 d ), where s(y d ) = ∂ y ∂ y d (y d ) is the so-called stretching coefficient. The final change of coordinate is then given by : d ) for y d ≤ y 0 d t 0 (y d ) for y d > y 0 d .

Fig. 5 . 10 :

 510 Fig. 5.10: Example of a coordinate mapping ỹ(y d ) used for the APML (black solid line). The graph of y d (y) (blue solid line) is continued by a straight line t 0 (y d ) tangent at y 0 d (red dashed line) to avoid the singular behaviour at y d = y d .

Fig. 5 . 11 :

 511 Fig. 5.11: Field maps of the logarithm of the norm of H z , E x and E y for the dielectric slit grating at λ 0 = 0.999y d,+ (same parameters as in part 5.2.4.2). (a): classical PML with inefficient damping of H z in the bottom PML. (b): APML where the H z field is correctly damped in the bottom PML. For both cases the thickness of the PML is ĥ-= 1.1y d,+ .

Fig.

  Fig.5.12: Modulus of the u n for the three propagating orders with adapted (dashed lines) and classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except for the grazing one (n = 1) as expected. This drawback is bypassed when using the adaptative PML.

Tab. 5 . 3 :

 53 0.1561 0.3959 0.1776 0.2703 0.9999 FEM + PML 0.1904 0.4118 0.1927 0.2481 1.0430 Diffraction efficiencies R 0 , T -1 , T 0 and T +1 of the four propagating orders, and energy balance B = R 0 + T -1 + T 0 + T +1 , computed by three methods: RCWA (line 1), FEM formulation with APML (line 2), FEM formulation with classical PML (line 3).

Fig. 5 . 13 :

 513 Fig. 5.13: Mean value of the norm of H z along the outer boundary of the bottom PML γ = |H z (-ĥ-)| x , for λ 0 approaching the Wood's anomaly y d,+ by inferior values (λ 0 = (1 -10 -n )y d,+ , red squares) and by superior value (λ 0 = (1 + 10 -n )y d,+ , blue circles) as a function of n.

pFig. 5 . 14 :

 514 Fig. 5.14: Scheme and notations of the studied bi-gratings.

Fig. 5 . 15 :

 515 Fig. 5.15: Degrees of freedom of a second order tetrahedral element.

2 -d x / 2 d y / 2 -

 222 , y, z) = ∑ (n,m)∈Z 2 u d,x n,m (z) e i (α n x+β m y) , (5.73)withα n = α 0 + 2 π d x n, β m = β 0 + 2 π d y m and u d,x n,m (z) = 1 d x d y d x /d y /2 E dx (x, y, z) e -i (α n x+β m y) dx dy .(5.74)By inserting the decomposition of Eq. (5.73), which is satisfied by E d x everywhere but in the groove region, into the Helmholtz propagation equation, one can express Rayleigh coefficients in the substrate and the superstrate as follows:u d,x n,m (z) = e x,p n,m e -i γ + n,m z + e x,c n,m e i γ + n,m z (5.75) with γ ± 2 n,m = k ± 2α 2 nβ 2 m, where γ n,m (or -i γ n,m ) is positive. The quantity u d,x n,m is the sum of a propagative plane wave (which propagates towards decreasing values of z, superscript p) and of a counterpropagative one (superscript c). The OWC verified by E d imposes: ∀(n, m) ∈ Z 2 e x,p n,m = 0 for z > z 0 e x,c n,m = 0 for z < z N(5.76)

T,

  (z c ) • e c n,m (z c ) for z c > z 0 (z c ) • e p n,m (z c ) for z c < z N

Fig. 5 . 16 :

 516 Fig. 5.16: Configuration of the studied cases.

Fig. 5 . 17 :

 517 Fig. 5.17: Diffractive element with vertical edges (a). ℜe{E x } in V/m (b).

FMMTab. 5 . 4 :

 54 Energy balance[START_REF] Li | New formulation of the fourier modal method for crossed surface-relief gratings[END_REF].

Fig. 5 . 18 :

 518 Fig. 5.18: Diffractive element with oblique edges (a). ℜe{E y } in V/m (b).

Fig. 5 . 19 :

 519 Fig. 5.19: Diffractive element with oblique edges (a). ℜe{E z } in V/m (b).

Tab. 5 . 6 :

 56 Energy balance[START_REF] Bruno | Numerical solution of diffraction problems: a method of variation of boundaries. III. doubly periodic gratings[END_REF].

Fig. 5 . 20 :

 520 Fig. 5.20: Lossy diffractive element with vertical edges (a). ℜe{E y } in V/m (b).

Tab. 5 . 7 :

 57 Comparison with[START_REF] Moharam | Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings[END_REF][START_REF] Li | New formulation of the fourier modal method for crossed surface-relief gratings[END_REF][START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF] and energy balance.

  and Fig. 5.21a are the following: λ 0 = 1, ϕ 0 = ψ 0 = 0 • , d x = d y = 0.3 , a = 0.1, b = 0.05, R = 0.15, h = 500 nm, ε + = ε g = 1, ε g = -21 + 20 i and ε -= 2.25.

Fig. 5 . 21 :Tab. 5 . 8 :

 52158 Fig. 5.21: Torus parameters (a). Coarse mesh of the computational domain (b).

Fig. 5 . 22 :

 522 Fig. 5.22: Convergence of R 0,0 in function of N m (circular apertures crossed-grating).

Fig. 5 . 23 :

 523 Fig. 5.23: Computation time and number of DOF as a function of N M .

Figure 5 .

 5 Figure 5.23 shows the computation time required to perform the whole FEM computational process for a system made of a number of DOF indicated on the right-hand ordinate. It is of importance to note that for values of N M lower than 3, the problem can be solved in less than a minute on a standard laptop (4Go RAM, 2 × 2 GHz) with 3 significant digits on the diffraction efficiencies. This accuracy is more than sufficient in numerous experimental cases. Furthermore, as far as integrated values are at stake, relatively coarse meshes (N M ≈ 1) can be used trustfully, authorizing fast geometric, spectral or polarization studies.Nowadays, the efficiency of the numerical algorithms for sparse matrix algebra together with the available power of computers and the fact that the problem reduces to a basic cell with a

E 1 (

 1 j (α 0 x + β 0 y + γ j z) j (α 0 x + β 0 yγ j z)

  components of H 1 can be expressed in function of tangential components of E 1 : invariance and linearity of the problem along Ox and Oy, the following notations are adopted. (5.80) and Eq. (5.84) and letting M = B -1 , it comes for the j th layer: Ψ(x, y, z) = exp(i(α 0 x + β 0 y)) continuity of Ψ at the interface z = z j leads to:Φ j+1 (z j ) = Π -1 j+1 Π j Φ j (z j ).(5.89)Normal components can be deduced using Eqs. (5.83,5.84).Traveling inside the j + 1 th layer Using Eq. (5.80), a simple phase shift allows to travel from z = z j to z = z j+1 = z je j+1 :Φ j+1 (z j+1 ) =

Reflection and transmission coefficients

The last step consists in the determination of the first term Φ 0 , which is not entirely known, since the problem definition only specifies U 0,+ x and U 0,+ y , imposed by the incident field E 0 . Let us make the use of the OWC hypothesis verified by E d 1 (see Eq. (5.62)). This hypothesis directly translates the fact that none of the components of E d 1 can either be traveling down in the superstrate or up in the substrate: U N+1,- y = U N+1,- x = 0. Therefore, the four unknowns U 0,- x , U 0,- y ,U N+1,+ y and U N+1,+ x , i.e. transverse components of the vector fields reflected and transmitted by the stack, verify the following equation system:

This allows to extend the definition of transmission and reflection widely used in the scalar case. Finally, Φ N+1 is entirely defined. Making the use of the recurrence relation of Eq. ( 5.91) and of Eq. (5.80) leads to an analytical expression for E d 1 in each layer.